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)

Specifically, for

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

where $$\mathbf{S}$$ is the solution vector, we solve

\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}.$$