Time Advance

Compressible Advance

To advance the fully compressible solution in time, ERF uses a 3rd order Runge-Kutta method. By default, acoustic substepping which solves the perturbational equations implicitly in the vertical direction is used in each Runge-Kutta stage, following the approach of Klemp, Skamarock and Dudhia (2007) However, there is a run-time option to turn off the substepping completely.

Specifically, the 3rd order Runge-Kutta method solves

\[\frac{d \mathbf{S}}{dt} = f(\mathbf{S})\]

where \(\mathbf{S}\) is the solution vector, in the following three steps:

\[ \begin{align}\begin{aligned}\mathbf{S}^{*} &=& \mathbf{S}^n + \frac{1}{3} \Delta t f(\mathbf{S}^n)\\\mathbf{S}^{**} &=& \mathbf{S}^n + \frac{1}{2} \Delta t f(\mathbf{S}^{*})\\\mathbf{S}^{n+1} &=& \mathbf{S}^n + \Delta t f(\mathbf{S}^{**})\end{aligned}\end{align} \]

Anelastic Advance

When solving the anelastic rather than fully compressible equations, ERF supports a 2nd order Runge-Kutta method and the 2nd order midpoint method (both with no substepping):

The 2nd order Runge-Kutta method solves

\[\frac{d \mathbf{S}}{dt} = f(\mathbf{S})\]

where \(\mathbf{S}\) is the solution vector, in the following two steps:

\[ \begin{align}\begin{aligned}\mathbf{S}^{*} &=& \mathbf{S}^n + \Delta t f(\mathbf{S}^n)\\\mathbf{S}^{n+1} &=& \mathbf{S}^n + \frac{\Delta t}{2} ( f(\mathbf{S}^{n}) + f(\mathbf{S}^{*}) )\end{aligned}\end{align} \]

The midpoint stages are given by:

\[ \begin{align}\begin{aligned}\mathbf{S}^{*} &=& \mathbf{S}^n + \frac{\Delta t}{2} f(\mathbf{S}^n)\\\mathbf{S}^{n+1} &=& \mathbf{S}^n + \Delta t f(\mathbf{S}^{*})\end{aligned}\end{align} \]

Acoustic Sub-stepping

When solving the fully compressible equation set, by default we substep the acoustic modes within each Runge-Kutta stage.

Recall the equations in the following form, here defining \(\mathbf{R}\) for each equation to include all additional terms that contribute to the time evolution.

\[ \begin{align}\begin{aligned}\frac{\partial \rho}{\partial t} &=& - \nabla \cdot (\rho \mathbf{u}) = R_\rho,\\\frac{\partial (\rho \mathbf{u})}{\partial t} &=& - \nabla \cdot (\rho \mathbf{u} \mathbf{u} + p^\prime I) + {\mathbf F}_\mathbf{u} = \mathbf{R}_\mathbf{u}\\\frac{\partial (\Theta)}{\partial t} &=& - \nabla \cdot (\mathbf{u} \Theta) + \rho \alpha_{T}\ \nabla^2 \theta = R_{\Theta},\\\frac{\partial (\rho C)}{\partial t} &=& - \nabla \cdot (\rho \mathbf{u} C) + \rho \alpha_{C}\ \nabla^2 C = R_{\rho C},\end{aligned}\end{align} \]

where we have defined \(\mathbf{U} = (U,V,W) = \rho \mathbf{u} = (\rho u, \rho v, \rho w)\) , \(\Theta = \rho \theta\) and \(\mathbf{F}_\mathbf{U} = (F_U, F_V, F_W) = \rho^\prime \mathbf{g} + \nabla \cdot \tau + \mathbf{F}\)

Using the relation \(\nabla p = \gamma R_d \pi \nabla \Theta,\) where the Exner function \(\pi = (p/p_0)^\kappa\) with \(\kappa = R_d / c_p,\) we can re-write the unapproximated momentum equations

\[ \begin{align}\begin{aligned}\frac{\partial U}{\partial t} &=& - \nabla \cdot (\mathbf{u} U) - \gamma R_d \pi \frac{\partial \Theta^\prime}{\partial x} + F_u\\\frac{\partial V}{\partial t} &=& - \nabla \cdot (\mathbf{u} V) - \gamma R_d \pi \frac{\partial \Theta^\prime}{\partial y} + F_v\\\frac{\partial W}{\partial t} &=& - \nabla \cdot (\mathbf{u} W) - \gamma R_d \pi \frac{\partial \Theta^\prime}{\partial z} - g (\overline{\rho} \frac{\pi^\prime}{\overline{\pi}} - \rho^\prime) + F_w\end{aligned}\end{align} \]

We then define new perturbational quantities, e.g., \(\mathbf{U}^{\prime \prime} = \mathbf{U} - \mathbf{U}^t\) where \(\mathbf{U}^t\) is the momentum at the most recent time reached by the “large” time step, which could be \(t^{n}\) or one of the intermediate Runge-Kutta steps. Then the acoustic substepping evolves the equations in the form

\[ \begin{align}\begin{aligned}U^{\prime \prime, \tau + \delta \tau} - U^{\prime \prime, \tau} &= \delta \tau \left( -\gamma R_d \pi^t \frac{\partial \Theta^{\prime \prime, \tau}}{\partial x} + R^t_U \right)\\V^{\prime \prime, \tau + \delta \tau} - V^{\prime \prime, \tau} &= \delta \tau \left( -\gamma R_d \pi^t \frac{\partial \Theta^{\prime \prime, \tau}}{\partial y} + R^t_V \right)\end{aligned}\end{align} \]

\[\begin{split}W^{\prime \prime, \tau + \delta \tau} - W^{\prime \prime, \tau} = \delta \tau \biggl( &-\gamma R_d \pi^t \frac{\partial (\beta_1 \Theta^{\prime \prime, \tau} + \beta_2 \Theta^{\prime \prime, \tau + \delta \tau} ) }{\partial z} \\ & - \frac{\gamma R_d^2}{c_v} \frac{\pi^t}{\Theta^t_m} \frac{\partial \Theta^t_m}{\partial z} (\beta_1 \Theta^{\prime \prime, \tau} + \beta_2 \Theta^{\prime \prime, \tau + \delta \tau} ) \\ & - g (\beta_1 \rho^{\prime \prime, \tau} + \beta_2 \rho^{\prime \prime, \tau + \delta \tau } ) \\ & + R^t_W \biggr)\end{split}\]

\[\Theta^{\prime \prime, \tau + \delta \tau} - \Theta^{\prime \prime, \tau} = \delta \tau \left( -\frac{\partial (U^{\prime \prime, \tau + \delta \tau} \theta^t)}{\partial x} -\frac{\partial (V^{\prime \prime, \tau + \delta \tau} \theta^t)}{\partial y} -\frac{\partial \left(( \beta_1 W^{\prime \prime, \tau} + \beta_2 W^{\prime \prime, \tau + \delta \tau} ) \theta^t\right)}{\partial z} + R^t_{\Theta} \right)\]

\[\rho^{\prime \prime, \tau + \delta \tau} - \rho^{\prime \prime, \tau} = \delta \tau \left( - \frac{\partial U^{\prime \prime, \tau + \delta \tau }}{\partial x} - \frac{\partial V^{\prime \prime, \tau + \delta \tau }}{\partial y} - \frac{\partial (\beta_1 W^{\prime \prime, \tau} + \beta_2 W^{\prime \prime, \tau + \delta \tau})}{\partial z} + R^t_{\rho} \right)\]

where \(\beta_1 = 0.5 (1 - \beta_s)\) and \(\beta_2 = 0.5 (1 + \beta_s)\), and \(\Theta_m = \Theta (1 + \frac{R_v}{R_d} q_v)\) is the moist analogue of \(\Theta\) which accounts for the contribution of water vapor to the pressure.

The second term in the \(W^{\prime \prime}\) equation above arises from differentiating \(p = R_d \pi \Theta_m\) with respect to \(z\), retaining the variation of \(\pi\) as well as that of \(\Theta_m\). Note that it is expressed entirely in terms of quantities at the most recent large time step or Runge-Kutta stage; unlike the formulation in equation (14) of Klemp, Skamarock and Dudhia (2007), it makes no reference to the hydrostatically stratified base state \((\overline{\rho}, \overline{\pi})\). This is the same choice made in the MPAS-A technical note, and yields coefficients for the double-prime terms that are exactly those obtained from the acoustic equations written in terms of the full variables rather than as perturbations from a specified reference sounding.

\(\beta_s\) is the acoustic step off-centering coefficient. When we do implicit substepping, we use the typical WRF value of 0.1. This off-centering is intended to provide damping of both horizontally and vertically propagating sound waves by biasing the time average toward the future time step. This may be controlled with the erf.beta_s parameter.

Over steep terrain more off-centering may be needed. With erf.use_variable_beta_s = true, \(\beta_s\) is instead defined on each \(w\)-face from the slopes of the coordinate surface there,

\[\beta_s = \max\left(\beta_{s,min}, \min\left(\beta_{s,max}, h_\xi^2 + h_\eta^2\right)\right),\]

with \(\beta_{s,min}\) and \(\beta_{s,max}\) set by erf.beta_s_min (default 0.1) and erf.beta_s_max (default 0.8). The vertical flux on a face is then weighted by \(\beta_1\) and \(\beta_2\) of that face, and the \(W^{\prime \prime}\) equation on a face by \(\beta_2\) of that face, so that the entries of the tridiagonal system for the face \(k\) carry \(\beta_2\) of face \(k\) times \(\beta_2\) of face \(k-1\), \(k\) or \(k+1\). Without terrain, \(\beta_s = \beta_{s,min}\) everywhere.

To solve the coupled system, we first evolve the equations for \(U^{\prime \prime, \tau + \delta \tau}\) and \(V^{\prime \prime, \tau + \delta \tau}\) explicitly using \(\Theta^{\prime \prime, \tau}\) which is already known. We then solve a tridiagonal system for \(W^{\prime \prime, \tau + \delta \tau}\), and once \(W^{\prime \prime, \tau + \delta \tau}\) is known, we update \(\rho^{\prime \prime, \tau + \delta \tau}\) and \(\Theta^{\prime \prime, \tau + \delta \tau}.\)

In addition to the acoustic off-centering, divergence damping is also applied to control horizontally propagating sound waves.

\[p^{\prime\prime,\tau*} = p^{\prime\prime,\tau} + \beta_d \left( p^{\prime\prime,\tau} - p^{\prime\prime,\tau-\delta\tau} \right)\]

where \(\tau*\) is the forward projected value used in RHS of the acoustic substepping equations for horizontal momentum. According to Skamarock et al, this is equivalent to including a horizontal diffusion term in the continuity equation. A typical damping coefficient of \(\beta_d = 0.1\) is used by default, as in WRF; it can be set with erf.beta_d, which must lie in \([0, 1]\).