Governing Equations

ERF can be run in two different modes: in the first, ERF solves the fully compressible fluid equations, in the second, ERF solves a modified set of equations which approximates the density field with the hydrostatic density and imposes the anelastic constraint on the velocity field.

In compressible mode, ERF solves partial differential equations expressing conservation of mass, momentum, potential temperature, and scalars (such as moisture variables) subject to an equation of state.

In anelastic mode, ERF solves partial differential equations expressing conservation of momentum, potential temperature, and scalars (such as moisture variables), as well the anelastic constraint on the velocity.

Below \(\rho_d, T, \theta_{d}\), and \(p\) are the dry-air density, temperature, dry potential temperature, and pressure, respectively; these variables are all defined at cell centers. \(\phi\) is an advected scalar, also defined at cell centers. \(\mathbf{u}\) and \((\rho_d \mathbf{u})\) are the velocity and momentum, respectively, and are defined on faces.

In the compressible moist formulation, the prognostic density is the dry-air density \(\rho_d\). When total moist density is needed, ERF forms it from the dry density and the water mixing ratios. In the EOS utility functions, rho denotes \(\rho_d\), and rhotheta denotes \(\rho_d \theta_d\).

Compressible Equations

The first three equations governing fully compressible flow are

\[ \begin{align}\begin{aligned}\frac{\partial \rho_d}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u}),\\\frac{\partial (\rho_d \mathbf{u})}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u} \mathbf{u}) - \frac{1}{1 + q_t} ( \nabla p^{\prime} - \delta_{i,3}\mathbf{B} ) - \nabla \cdot \boldsymbol{\tau} + \mathbf{F}_{u},\\\frac{\partial (\rho_d \theta_d)}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u} \theta_d) + \nabla \cdot (\rho_d \alpha_{\theta}\ \nabla \theta_d) + F_{\theta} + H_{n} + H_{p},\end{aligned}\end{align} \]

supplemented with the equation of state as given below.

Anelastic Equations

The first two equations for the anelastic formulation are

\[ \begin{align}\begin{aligned}\frac{\partial (\rho_0 \mathbf{u})}{\partial t} &= - \nabla \cdot (\rho_0 \mathbf{u} \mathbf{u}) - \frac{1}{1 + q_t} ( \nabla p^\prime + \delta_{i,3}\mathbf{B} ) - \nabla \cdot \boldsymbol{\tau} + \mathbf{F}_{u},\\\frac{\partial (\rho_0 \theta_d)}{\partial t} &= - \nabla \cdot (\rho_0 \mathbf{u} \theta_d) + \nabla \cdot (\rho_0 \alpha_{\theta}\ \nabla \theta_d) + F_{\theta} + H_{n} + H_{p},\end{aligned}\end{align} \]

supplemented with the constraint

\[\nabla \cdot (\rho_0 \mathbf{u}) = 0\]

(Dry and Moist) Scalars

We supplement the above equations with the following equations for advected scalars (\(\phi\)) and precipitating (\(\mathbf{q_{p}}\)) and non-precipitating (\(\mathbf{q_{n}}\)) moisture variables (identical for compressible and anelastic)

\[ \begin{align}\begin{aligned}\frac{\partial (\rho_d \boldsymbol{\phi})}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u} \boldsymbol{\phi}) + \nabla \cdot ( \rho_d \alpha_{\phi}\ \nabla \boldsymbol{\phi}) + \mathbf{F}_{\phi},\\\frac{\partial (\rho_d \mathbf{q_{n}})}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u} \mathbf{q_{n}}) + \nabla \cdot (\rho_d \alpha_{q} \nabla \mathbf{q_{n}}) + \mathbf{F_{n}} + \mathbf{G_{p}},\\\frac{\partial (\rho_d \mathbf{q_{p}})}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u} \mathbf{q_{p}}) + \partial_{z} \left( \rho_d \mathbf{w_{t}} \mathbf{q_{p}} \right) + \mathbf{F_{p}}.\end{aligned}\end{align} \]

The non-precipitating water mixing ratio vector \(\mathbf{q_{n}} = \left[ q_v \;\; q_c \;\; q_i \right]\) includes water vapor, \(q_v\), cloud water, \(q_c\), and cloud ice, \(q_i\), although some microphysical moisture models may not include cloud ice; similarly, the precipitating water mixing ratio vector \(\mathbf{q_{p}} = \left[ q_r \;\; q_s \;\; q_g \right]\) involves rain, \(q_r\), snow, \(q_s\), and graupel, \(q_g\), though some models may not include these terms. The source terms for moisture variables, \(\mathbf{F_{p}}\), \(\mathbf{F_{n}}\), \(\mathbf{G_{p}}\), and their corresponding impact on potential temperature, \(H_{n}\) and \(H_{p}\), and the terminal velocity, \(\mathbf{w_{t}}\) are specific to the employed model. See the Microphysics section for more details.

Height-Following Terrain Coordinates

Consider two coordinate systems that correspond to a terrain-following grid, \(\mathbf{X}\), and a flat cartesian grid, \(\mathbf{Z}\), with axes given by

\[\mathbf{X} = \left[ x \; y \; z \right]^{\intercal}, \quad \quad \mathbf{\Xi} = \left[ \xi \; \eta \; \zeta \right]^{\intercal},\]

and

\[x = \xi, \quad \quad y = \eta, \quad \quad z = h \left(\xi, \, \eta, \, \zeta \right).\]

Only the vertical coordinate in the physical domain is deformed by the terrain-fitting. To account for isotropic lateral grid stretching as represented by “map factors” \(m_x = m_y = m\) as in WRF, we augment the coordinate transform above with stretching in the lateral directions only.

These combined transformations yield the following Jacobian, \(\bar{\mathbf{J}}\), and inverse Jacobian, \(\bar{\mathbf{T}}\), matrices

\[\begin{split}\bar{\mathbf{J}} = \begin{bmatrix} \frac{1}{m} & 0 & 0 \\ 0 & \frac{1}{m} & 0\\ h_{\xi} & h_{\eta} & h_{\zeta} \\ \end{bmatrix}, \quad \quad \bar{\mathbf{T}} = \mathbf{J}^{-1} = \frac{m^2}{h_{\zeta}} \begin{bmatrix} \frac{h_{\zeta}}{m} & 0 & 0 \\ 0 & \frac{h_{\zeta}}{m} & 0\\ -\frac{h_{\xi}}{m} & -\frac{h_{\eta}}{m} & \frac{1}{m^2} \\ \end{bmatrix} = \begin{bmatrix} m & 0 & 0 \\ 0 & m & 0\\ -\frac{h_{\xi}}{h_\zeta}m & -\frac{h_{\eta}}{h_\zeta}m & \frac{1}{h_\zeta} \\ \end{bmatrix}.\end{split}\]

In the above, \(J = \left| \bar{\mathbf{J}} \right |= h_{\zeta} / m^2\) is the Jacobian determinant. To explicitly close the governing equations in terrain-following coordinates, we provide relations for the gradient of a scalar (\(f\)) and divergence of a vector (\(\mathbf{F}\)):

\[ \begin{align}\begin{aligned}\nabla_{\mathbf{X}} f &= \bar{\mathbf{T}}^{\intercal} \nabla_{\mathbf{Z}} f,\\\nabla_{\mathbf{X}} \cdot \left( \mathbf{F} \right) &= \frac{1}{J} \nabla_{\mathbf{Z}} \cdot \left( J \bar{\mathbf{T}} \mathbf{F}\right).\end{aligned}\end{align} \]

Vector rotation of the fluid velocity yields \(J \bar{\mathbf{T}} \mathbf{u} = \left[h_{\zeta}u/m, \;\; h_{\zeta}v/m, \;\; \omega/m^2 \right]^{\intercal}\), where \(\omega = w -h_{\xi} u m - h_{\eta} v m\) is the vertical velocity that is normal to the top/bottom faces of the grid cells.

Background (reference) state

Pressure and density perturbations are defined with respect to a hydrostatically stratified background state, i.e.

\[p = p_{0}(z) + p^\prime \hspace{24pt} \rho = \rho_{0}(z) + \rho^\prime\]

with

\[\frac{d p_{0}}{d z} = - \rho_{0} g\]

Equation of state (compressible only)

In the fully compressible formulation, the total pressure is computed as

\[p = P_{00} \left( \frac{R_d \rho_d \theta_m}{P_{00}} \right)^\gamma\]

where \(\gamma = c_{p} / (c_{p} - R_{d})\) and

\[\theta_m = \theta_d (1 + \frac{R_v}{R_d} q_v)\]

is the moist potential temperature. This is the only place \(\theta_m\) is used; we evolve \(\theta_d\) above. In the above, \(R_d\), \(c_p\), \(P_{00} = 1\times10^{5}\) are the gas constant, specific heat capacity for dry air, and reference pressure, respectively.

Define

\[M(q_v) = 1 + \frac{R_v}{R_d} q_v.\]

Then the EOS can be written equivalently as

\[p = P_{00} \left( \frac{R_d \rho_d \theta_d M(q_v)}{P_{00}} \right)^\gamma.\]

Here \(q_v\) is the vapor mixing ratio per unit dry-air mass. The quantity \(\rho_d \theta_d\) is the rhotheta state variable used by the EOS utility functions, and EOS pressure arguments and return values are in Pa.

The same EOS implies

\[p = \rho_d R_d T M(q_v),\]

and therefore

\[\rho_d = \frac{p}{R_d T M(q_v)}.\]

This matches getRhogivenTandPress in Source/Utils/ERF_EOS.H.

The EOS utilities rely on the constant relation

\[\kappa \equiv \frac{R_d}{c_p} = \frac{\gamma - 1}{\gamma}, \qquad \frac{1}{\gamma} = 1 - \kappa.\]

This identity makes the pressure, temperature, Exner, and density inverse relations in Source/Utils/ERF_EOS.H algebraically consistent.

The functions in Source/Utils/ERF_EOS.H implement these relations:

Function	Relation
getThgivenTandP	\(\theta_d = T(P_{00}/p)^{R_d/c_p}\)
getTgivenPandTh	\(T = \theta_d(p/P_{00})^{R_d/c_p}\)
getPgivenRTh	\(p = P_{00}(R_d\rho_d\theta_d M/P_{00})^\gamma\)
getRhoThetagivenP	inverse of getPgivenRTh
getRhogivenTandPress	\(\rho_d = p/(R_d T M)\)
getRhogivenThetaPress	density from \(\theta_d\), \(p\), and \(q_v\)
getExnergivenP	\(\Pi = (p/P_{00})^{R_d/c_p}\)
getdPdRgivenConstantTheta	\((\partial p/\partial\rho_d)_{\theta_d,q_v} = \gamma p/\rho_d\)

Additional terms

• \(\boldsymbol{\tau}\) is the viscous stress tensor,

  \[\tau_{ij} = -2\mu \sigma_{ij},\]

with \(\sigma_{ij} = S_{ij} -D_{ij}\) being the deviatoric part of the strain rate, and

\[S_{ij} = \frac{1}{2} \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right), \hspace{24pt} D_{ij} = \frac{1}{3} S_{kk} \delta_{ij} = \frac{1}{3} (\nabla \cdot \mathbf{u}) \delta_{ij},\]

• \(\mathbf{F}_{u}\) and \(F_{\theta_d}\) are the forcing terms described in Physical Forcings,

• \(\mathbf{B} = -(\rho - \rho_{0})\mathbf{g}\) is the buoyancy term described in Buoyancy,

• \(\mathbf{g} = (0,0,-g)\) is the gravity vector,

• The dry potential temperature \(\theta_d\) is defined from temperature \(T\), pressure \(p\), and reference pressure \(P_{00} = 10^{5}\) Pa as

\[\theta_d = T \left( \frac{P_{00}}{p} \right)^{R_d / c_p}.\]

(In the anelastic case, \(p\) is replaced by \(p_0\) in the relationship between \(\theta_d\) and \(T\).)

Assumptions

The assumptions involved in deriving these equations from first principles are:

• Continuum behavior

• Ideal gas behavior with constant specific heats (\(c_p,c_v\)). In dry configurations, \(p = \rho_d R_d T\). In moist configurations, ERF uses dry density and vapor mixing ratio per dry-air mass:

  \[p = \rho_d R_d T \left(1 + \frac{R_v}{R_d}q_v\right).\]

  This is equivalent to the moist potential temperature factor in the compressible EOS above. The dry shorthand \(p = \rho R_d T\) is recovered only when \(q_v = 0\).

• Viscous heating is negligible

• No chemical reactions, second order diffusive processes or radiative heat transfer

• Newtonian viscous stress with no bulk viscosity contribution (i.e., \(\kappa S_{kk} \delta_{ij}\))

• Depending on the simulation mode, the transport coefficients \(\mu\), \(\rho\alpha_{\phi}\), and \(\rho\alpha_{\theta}\) may correspond to the molecular transport coefficients, turbulent transport coefficients computed from an LES or PBL model, or a combination. See the sections on DNS vs. LES modes and PBL schemes for more details.