Hydrostatic balance

During the setup-phase of a run, CGCAM will compute the thermodynamic state such that the hydrostatic balance is satisfied numerically on a particular point distribution in \(z\). This point distribution corresponds to the cell center locations for a Cartesian grid type (i_grid=1), the \(z\) point distribution contained in the grid file for a terrain grid (i_grid=3), and a uniform point distribution for a general mesh (i_grid=0). With the thermodynamic value sampling points being equal to the cell centers in the Cartesian grid case, the hydrostatic balance will be satisfied at all points within the domain. For the other cases, however, interpolation errors, both when the thermodynamic variables are interpolated to the cell centers, and when the pressure is interpolated from the cell centers to the cell faces, will upset the hydrostatic balance. While the severity of these errors are problem dependent, they will result in non-physical flows at early times during the simulation since pressure and buoyancy forces are out of balance. Currently this problem can be avoided by specifying the input i_hydrostat=1, in which case the conservation laws are modified by analytically subtracting the hydrostatic force balance. While this modification is straightforward for the momentum equation is is somewhat more involved for the energy equation. For the momentum equation we simply decompose the pressure gradient and buoyancy force terms into background and fluctuating components and then subtract the background (hydrostatic) part:

\[\begin{split}\begin{eqnarray*} \frac{\partial \rho u_i}{\partial t} & = & -\frac{\partial}{\partial x_j}\Big( u_i \; \rho u_j \Big) + \frac{\partial \sigma_{ij}}{\partial x_j} - \frac{\partial p}{\partial x_i} - \rho g\delta_{i3}\\[0.1in] & = & -\frac{\partial}{\partial x_j}\Big( u_i \; \rho u_j \Big) + \frac{\partial \sigma_{ij}}{\partial x_j} - \;\underbrace{\frac{\partial \bar{p}}{\partial x}\delta_{i1} - \frac{\partial \bar{p}}{\partial y}\delta_{i2}} _{{\rm usually}~0} - \;\underbrace{ \left( \frac{\partial \bar{p}}{\partial z} + \bar{\rho} g\right)}_{=0} - \frac{\partial p^\prime}{\partial x_i} - \rho^\prime g\delta_{i3}\\[0.1in] & = & -\frac{\partial}{\partial x_j}\Big( u_i \; \rho u_j \Big) + \frac{\partial \sigma_{ij}}{\partial x_j} - \frac{\partial p^\prime}{\partial x_i} - \rho^\prime g\delta_{i3}. \end{eqnarray*}\end{split}\]

Here \(\sigma_{ij}\) is the viscous stress as defined in section Governing Equations. Note that the horizontal background pressure gradient forces have been dropped along with the vertical hydrostatic balance. While these terms should be zero analytically in most cases (but see Data Assimilation), interpolation errors along with the approximation that the pressure be constant over a cell face typically renders these terms non-zero for non-Cartesian grid cells. By removing the horizontal background pressure gradient terms we are guaranteed that a simple background pressure field (i.e. \(\bar{p}(z)\)) exerts no force imbalance in any direction.

Since one component of the energy equation involves \(u_i({\rm momentum})_i\), the hydrostatic balance reappears in the energy equation as a work term. The situation is more subtle, however, since no work exists when the velocity is zero. Thus the energy equation is already in hydrostatic balance, resulting in \(\partial\rho E/\partial t = 0\) when \(u_i=0\), provided there is no net heat transfer. The situation is more complicated when the velocity is non-zero, however, as erroneous work will be associated with any imbalance in the background pressure forces. We can demonstrate this fact if we decompose the pressure work and buoyancy terms into background and fluctuating components, and if we expand the pressure work term, as follows:

\[\begin{split}\begin{eqnarray*} \frac{\partial \rho E}{\partial t} & = & -\frac{\partial}{\partial x_j}\Big( E \; \rho u_j \Big) + \frac{\partial u_i \sigma_{ij}}{\partial x_j} - \frac{\partial u_j p}{\partial x_j} - w \rho g - \frac{\partial q_j}{\partial x_j}\\[0.1in] & = & -\frac{\partial}{\partial x_j}\Big( E \; \rho u_j \Big) + \frac{\partial u_i \sigma_{ij}}{\partial x_j} - \frac{\partial u_j p^\prime}{\partial x_j} - \frac{\partial u_j \bar{p}}{\partial x_j} - w \bar{\rho} g - w \rho^\prime g - \frac{\partial q_j}{\partial x_j}\\[0.1in] & = & -\frac{\partial}{\partial x_j}\Big( E \; \rho u_j \Big) + \frac{\partial u_i \sigma_{ij}}{\partial x_j} - \frac{\partial u_j p^\prime}{\partial x_j} - \bar{p}\frac{\partial u_j }{\partial x_j} - u\!\!\underbrace{\frac{\partial \bar{p}}{\partial x}}_{{\rm usually}~0} - v\!\!\underbrace{\frac{\partial \bar{p}}{\partial y}}_{{\rm usually}~0} - w\underbrace{ \left(\frac{\partial \bar{p}}{\partial z} - \bar{\rho} g\right)}_{=0} - \\& & w \rho^\prime g - \frac{\partial q_j}{\partial x_j}\\[0.1in] & = & -\frac{\partial}{\partial x_j}\Big( E \; \rho u_j \Big) + \frac{\partial u_i \sigma_{ij}}{\partial x_j} - \frac{\partial u_j p^\prime}{\partial x_j} - \bar{p}\frac{\partial u_j }{\partial x_j} - w \rho^\prime g - \frac{\partial q_j}{\partial x_j}\\[0.1in] \end{eqnarray*}\end{split}\]

Note that an extra term, the background pressure multiplied by the velocity divergence, appears in the modified energy equation. This term is somewhat problematic since it is not in pure divergence form, and thus could potentially lead to global energy conservation errors. One way to guard against this possibility is to write the original pressure-work term as

\[\frac{\partial pu_j}{\partial x_j} \simeq \frac{\delta \widetilde{pu}_j}{\delta x_j},\]

where

\[\widetilde{pu}\vert_{i+1/2} = 1/2\left(p_{i}u_{i+1} + p_{i+1}u_{i}\right),\]

where \(i\) is the grid index. Then, according to Morinishi et al. (JCP 143, 1998)

\[\frac{\delta_1 \widetilde{pu_j}}{\delta_1 x_j} = p \frac{\delta_2 u_j}{\delta_2 x_j} + u_j\frac{\delta_2 p }{\delta_2 x_j},\]

where \(\delta_1/\delta_1 x\) and \(\delta_2/\delta_2 x\) are one and two-point finite differences, respectively. Thus if we use the special interpolation \(\widetilde{pu}_j\), the decomposition of the pressure work term done above should still result in a globally conservative form.

Hydrostataic balance removal in conjunction with data assimilation

The data assimilation forms are derived by subtracting the background balances from the flow equations. If the hydrostatic balance are also to be removed we can start with the results above to derive the corresponding data assimilation forms. Since the background equivalent of any prime term is null, only the advective and viscous stress terms in the momentum equation above have a background equivalent. Subtracting these, adding the measured momentum time rate of change, and using the notation \(\sigma^\prime_{ij} = \sigma_{ij} - \bar{\sigma}_{ij}\), we get

\[ \frac{\partial \rho u_i}{\partial t} = -\frac{\partial}{\partial x_j}\Big( u_i \; \rho u_j - \bar{u}_i \; \bar{\rho}\bar{u}_j \Big) + \frac{\partial \sigma_{ij}^\prime}{\partial x_j} - \frac{\partial p^\prime}{\partial x_j} - \rho^\prime g\delta_{i3} + \left(\frac{\partial \bar{\rho} \bar{u}_i}{\partial t}\right)_D\]

A similar procedure is used to derive the data assimilation form of the energy equation

\[\begin{split}\frac{\partial \rho E}{\partial t} = -\frac{\partial}{\partial x_j}\Big( E \; \rho u_j - \bar{E} \; \bar{\rho}\bar{u}_j \Big) + \frac{\partial}{\partial x_j}\left( u_i\sigma_{ij} - \bar{u}_i\bar{\sigma}_{ij}\right) - \\[0.1in] \frac{\partial u_j p^\prime}{\partial x_j} - \bar{p}\frac{\partial u^\prime_j }{\partial x_j} - w \rho^\prime g - \frac{\partial q^\prime_j}{\partial x_j} + \left(\frac{\partial \bar{\rho}\bar{E}}{\partial t}\right)_D\end{split}\]

Note that the viscous work term can be written as

\[u_i\sigma_{ij} - \bar{u}_i\bar{\sigma}_{ij} \; = \; u_i(\sigma^\prime_{ij} + \bar{\sigma}_{ij}) - \bar{u}_i\bar{\sigma}_{ij} \; = \; u_i\sigma^\prime_{ij} + (u_i - \bar{u}_i)\bar{\sigma}_{ij} \; = \; u_i\sigma^\prime_{ij} + u^\prime_i\bar{\sigma}_{ij}\]

This result may be used to rewrite the energy equation as

\[\frac{\partial \rho E}{\partial t} = -\frac{\partial}{\partial x_j}\Big( E \; \rho u_j - \bar{E} \; \bar{\rho}\bar{u}_j \Big) + \frac{\partial}{\partial x_j}\Big( u_i\Sigma^\prime_{ij} + u^\prime_i\bar{\sigma}_{ij} \Big) - \bar{p}\frac{\partial u^\prime_j }{\partial x_j} - w \rho^\prime g - \frac{\partial q^\prime_j}{\partial x_j} + \left(\frac{\partial \bar{\rho}\bar{E}}{\partial t}\right)_D\]

where \(\Sigma^\prime_{ij}=\sigma^\prime_{ij}-p^\prime\delta_{ij}\).