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 (
-\gamma R_d \pi^t \frac{\partial \Theta^{\prime \prime, \tau}}{\partial x} + R^t_U)\\V^{\prime \prime, \tau + \delta \tau} - V^{\prime \prime, \tau} &=& \delta \tau (
-\gamma R_d \pi^t \frac{\partial \Theta^{\prime \prime, \tau}}{\partial y} + R^t_V)\end{aligned}\end{align} \]
\[ \begin{align}\begin{aligned}\begin{split}W^{\prime \prime, \tau + \delta \tau} - W^{\prime \prime, \tau} &=& \delta \tau (
-\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}}\end{split}\\ \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 )\end{aligned}\end{align} \]
\[\Theta^{\prime \prime, \tau + \delta \tau} - \Theta^{\prime \prime, \tau} = \delta \tau (
-\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 (( \beta_1 W^{\prime \prime, \tau} + \beta_2 W^{\prime \prime, \tau + \delta \tau} ) \theta^t)}{\partial z} + R^t_{\Theta} )\]
\[\rho^{\prime \prime, \tau + \delta \tau} - \rho^{\prime \prime, \tau} = \delta \tau (
- \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} )\]
where \(\beta_1 = 0.5 (1 - \beta_s)\) and \(\beta_2 = 0.5 (1 + \beta_s)\) with \(\beta_s = 0.1\).
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}.\)