To advance the solution in time, ERF uses a 3rd order Runge-Kutta method with acoustic sub-stepping
in each Runge-Kutta stage, following the approach of Klemp, Skamarock and Dudhia (2007)
\[ \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} \]
Acoustic Sub-stepping
We sub-step the acoustic modes within each Runge-Kutta stage.
We first 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} \\
& - g \overline{\rho} \frac{R_d}{c_v} \frac{\pi^t}{\overline{\pi}}
\frac{ (\beta_1 \Theta^{\prime \prime, \tau} +
\beta_2 \Theta^{\prime \prime, \tau + \delta \tau} )}{\Theta^t} \\
& + 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)\) with \(\beta_s = 0.1\).
\(\beta_s\) is the acoustic step off-centering coefficient and 0.1 is the typical WRF value. 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.
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, as in
WRF.