A framework to evaluate IMEX schemes for atmospheric models

We present a new evaluation framework for implicit and explicit (IMEX) Runge–Kutta time-stepping schemes. The new framework uses a linearized nonhydrostatic system of normal modes. We utilize the framework to investigate the stability of IMEX methods and their dispersion and dissipation of gravity, Rossby, and acoustic waves. We test the new framework on a variety of IMEX schemes and use it to develop and analyze a set of second-order low-storage IMEX Runge–Kutta methods with a high Courant–Friedrichs–Lewy (CFL) number. We show that the new framework is more selective than the 2-D acoustic system previously used in the literature. Schemes that are stable for the 2-D acoustic system are not stable for the system of normal modes.

1 Introduction

Differences in phase speeds between slow and fast waves in atmospheric models motivate development of time-stepping schemes with an implicit component to avoid expensive time-step restrictions imposed by fast waves on explicit methods. The nonlinearity of the equations and the expense of solving globally coupled linear equations often impose a prohibitive cost on the solvers required by fully implicit methods and hybrid implicit–explicit (IMEX) schemes that leverage the strengths of both have become common. Here, we develop a new framework for evaluating IMEX methods for atmospheric modeling.

We follow approaches of Durran and Blossey (2012), Weller et al. (2013), Lock et al. (2014), Rokhzadi et al. (2018), and others to present an evaluation framework that is simpler than a full 3-D model while still containing the challenges associated with the presence of both slow and fast modes. Our framework is based on the normal mode analysis for systems introduced in Thuburn et al. (2002a, b) and Thuburn and Woollings (2005).

We focus on IMEX Runge–Kutta (RK) methods and their use in our primary application, the High-Order Method Modeling Environment (HOMME) dynamical core (Dennis et al., 2012; Taylor et al., 2020). HOMME is the nonhydrostatic atmospheric dynamical core of the US Department of Energy Exascale Earth System Model's (E3SM) (Rasch et al., 2019) atmosphere component. HOMME is formulated in horizontally explicit, vertically implicit (HEVI) form and is well suited for IMEX RK schemes in which terms that carry vertically propagating acoustic waves are treated implicitly.

We adapt the techniques of Thuburn and Woollings (2005), hereafter TW2005, to the specific system of equations and prognostic variables in HOMME as well as other dynamical cores. In particular, we use a system of normal modes for a mass-based vertically Lagrangian coordinate system with a Lorenz-staggered vertical discretization. We construct a spacetime operator for this system and study its properties, including stability, dispersion, and dissipation. Compared with the previously used 2-D acoustic system and the compressible Boussinesq equations (Durran and Blossey, 2012; Weller et al., 2013; Lock et al., 2014; Rokhzadi et al., 2018), this system provides more complexity and more closely resembles the equations used in modern dynamical cores. It contains a full set of modes: east- and west-propagating acoustic and gravity waves and westward-propagating Rossby waves. It is linearized about a hydrostatic reference state and has the commonly used constant-pressure boundary condition at the model top.

Using the new framework, we develop a family of second-order, high-CFL (Courant–Friedrichs–Lewy), low-storage IMEX RK schemes and analyze their suitability for operational use in E3SM's high-resolution science campaigns.

The remainder of this paper is structured as follows. In Sect. 2, we present the linearized system of equations associated with our formulation of the nonhydrostatic dynamics equations and compute its spatial numerical dispersion properties. Section 3 introduces the spacetime operator and describes our analysis of its numerical stability properties. In Sect. 4, we compare the stability diagrams of several schemes and develop a new family of IMEX RK methods with desirable stability and storage properties. In Sect. 4.4, we investigate convergence of IMEX methods with respect to vertical resolution. Section 5 concludes the paper.

2 Linearized system for normal modes

In this section, we define the linearized system of equations that corresponds to the HOMME nonhydrostatic dynamics model. Then we present a discretization – which we broadly define here to include the choices of prognostic variables, equation of state, boundary conditions, and staggering – of this system corresponding to HOMME's discretization. We obtain analytical and numerical frequencies for, respectively, the analytical and discretized systems and confirm that the discretization is nearly optimal. Therefore, we can subsequently use the system to investigate properties of IMEX spacetime operators.

2.1 Description of the system

In Thuburn et al. (2002a, b) and TW2005, the Euler equations for a dry adiabatic atmosphere are simplified to study normal modes. Various approximations about the geometry and Coriolis terms are made, and the systems are linearized about a hydrostatic reference state at rest. Furthermore, TW2005 presents such systems for different choices for thermodynamic variables, vertical coordinates, and equations of state. We use a vertical coordinate based on hydrostatic pressure (Laprise, 1992), where hybrid pressure levels are located on constant s surfaces, and s is the vertical Lagrangian coordinate satisfying $\dot{s}=\mathrm{0}$, following Lin (2004). Therefore, we adopt the system of Eqs. (20)–(24) in TW2005 for the shallow atmosphere approximation and a Lagrangian vertical coordinate:

$$\begin{array}{}\text{(1)}& {\displaystyle }{u}_{\mathrm{t}}=fv-\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}^{\mathrm{r}}}}{\displaystyle \frac{\partial p}{\partial x}}+{\displaystyle \frac{\partial \mathit{\varphi }}{\partial x}}\right)\text{(2)}& {\displaystyle }{v}_{\mathrm{t}}=-fu-\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}^{\mathrm{r}}}}{\displaystyle \frac{\partial p}{\partial y}}+{\displaystyle \frac{\partial \mathit{\varphi }}{\partial y}}\right)\text{(3)}& {\displaystyle }{w}_{\mathrm{t}}=-g{\displaystyle \frac{\mathit{\sigma }}{{\mathit{\sigma }}^{\mathrm{r}}}}-{\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}^{\mathrm{r}}}}{\displaystyle \frac{\partial p}{\partial \mathit{\theta }}}\text{(4)}& {\displaystyle }{\mathit{\varphi }}_{\mathrm{t}}=gw\text{(5)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{t}}=-{\mathit{\sigma }}^{\mathrm{r}}\left({\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial v}{\partial y}}\right).\end{array}$$

Here, u, v, and w are velocity components, p is pressure, ρ is density, ϕ is geopotential, g is the gravity constant, f is the Coriolis parameter, σ is pseudo-density defined with respect to the vertical coordinate (see Taylor et al., 2020, for details), and θ is potential temperature. The superscript “r” denotes variables defined by reference profiles of a linearized hydrostatic steady state with constant temperature T₀. The subscript “t” denotes partial differentiation with respect to time. Variables u, v, w, ϕ, p, and σ are first-order perturbed quantities, about the reference state, as follows from linear analysis.

After substituting single-mode solutions in which each field is proportional to $\exp (i{k}_{x}x+i{l}_{y}y-i\mathit{\omega }t)$, this formulation is equivalent to the system of Eqs. (20)–(24) in the isentropic coordinate from TW2005. With inclusion of the β effect as in Eqs. (55)–(56) of TW2005, this system is as follows:

$$\begin{array}{}\text{(6)}& {\displaystyle }-i\mathit{\omega }u=fv+{\displaystyle \frac{i{k}_x}{K^{\mathrm{2}}}}\mathit{\beta }u-i{k}_x\left({\displaystyle \frac{p}{{\mathit{\rho }}^{\mathrm{r}}}}+\mathit{\varphi }\right)\text{(7)}& {\displaystyle }-i\mathit{\omega }v=-fu+{\displaystyle \frac{i{k}_x}{K^{\mathrm{2}}}}\mathit{\beta }v-i{l}_y\left({\displaystyle \frac{p}{{\mathit{\rho }}^{\mathrm{r}}}}+\mathit{\varphi }\right)\text{(8)}& {\displaystyle }-i\mathit{\omega }w=-g{\displaystyle \frac{\mathit{\sigma }}{{\mathit{\sigma }}^{\mathrm{r}}}}-{\displaystyle \frac{p_{\mathit{\theta }}}{{\mathit{\sigma }}^{\mathrm{r}}}}\text{(9)}& {\displaystyle }-i\mathit{\omega }\mathit{\varphi }=gw\text{(10)}& {\displaystyle }-i\mathit{\omega }\mathit{\sigma }=-{\mathit{\sigma }}^{\mathrm{r}}(i{k}_{x}u+i{l}_{y}v),\end{array}$$

with linearized equation of state (EOS)

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{p}{p^{\mathrm{r}}}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-\mathit{\kappa }}}{\displaystyle \frac{\mathit{\sigma }}{{\mathit{\sigma }}^{\mathrm{r}}}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-\mathit{\kappa }}}{\displaystyle \frac{{\mathit{\varphi }}_{\mathit{\theta }}}{{\mathit{\varphi }}_{\mathit{\theta }}^{\mathrm{r}}}}.\end{array}$$

We also retain a version of the system of Eqs. (6)–(10) with time derivatives on the left-hand side:

$$\begin{array}{}\text{(12)}& {\displaystyle }{u}_{\mathrm{t}}=fv+{\displaystyle \frac{i{k}_x}{K^{\mathrm{2}}}}\mathit{\beta }u-i{k}_x\left({\displaystyle \frac{p}{{\mathit{\rho }}^{\mathrm{r}}}}+\mathit{\varphi }\right)\text{(13)}& {\displaystyle }{v}_{\mathrm{t}}=-fu+{\displaystyle \frac{i{k}_x}{K^{\mathrm{2}}}}\mathit{\beta }v-i{l}_y\left({\displaystyle \frac{p}{{\mathit{\rho }}^{\mathrm{r}}}}+\mathit{\varphi }\right)\text{(14)}& {\displaystyle }{w}_{\mathrm{t}}=-g{\displaystyle \frac{\mathit{\sigma }}{{\mathit{\sigma }}^{\mathrm{r}}}}-{\displaystyle \frac{p_{\mathit{\theta }}}{{\mathit{\sigma }}^{\mathrm{r}}}}\text{(15)}& {\displaystyle }{\mathit{\varphi }}_{\mathrm{t}}=gw\text{(16)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{t}}=-{\mathit{\sigma }}^{\mathrm{r}}(i{k}_{x}u+i{l}_{y}v).\end{array}$$

In addition to the variables and constants defined above, $\mathit{\kappa }=R/c_p$ is a thermodynamic constant, and k_x and l_y are horizontal wavenumbers with $K^{\mathrm{2}}=k_x^{\mathrm{2}}+l_y^{\mathrm{2}}$. Here and later in the text, the subscripts x and y in k_x and l_y do not denote differentiation in x or y; we use this notation to be consistent with the symbols for horizontal and vertical wavenumbers introduced in Weller et al. (2013) and Lock et al. (2014). The subscript θ denotes partial differentiation with respect to potential temperature. In Eqs. (6)–(10), (11), and (12)–(16), the variables ρ^r, σ^r, p^r, and derivative ${\mathit{\varphi }}_{\mathit{\theta }}^{\mathrm{r}}$ are variables defined by the reference profile of a linearized hydrostatic steady state with constant temperature T₀. The variables u, v, w, ϕ, p, and σ are first-order perturbed quantities about the reference state. All variables are scalar quantities.

We note that since $\dot{\mathit{\theta }}=\mathrm{0}$, the linear system associated with the TW2005 isentropic model is equivalent to the linear system derived with a vertically Lagrangian coordinate. We use the same bottom boundary condition (BC) ϕ=0 for the systems of Eqs. (6)–(10) and (12)–(16) as in TW2005. We replace the rigid lid boundary condition with a constant-pressure boundary condition, which for the perturbed-pressure variable becomes p_top=0. This BC is more typical for a mass-based vertical coordinate.

We define meridional wavenumber l_y=0, temperature of reference state T₀=250 K, depth of domain in the vertical direction D=10⁵ m, Coriolis parameters $\mathit{\beta }=\mathrm{1.619}\times {\mathrm{10}}^{-\mathrm{11}}$ ${\mathrm{s}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{1}}$ and $f=\mathrm{1.031}\times {\mathrm{10}}^{-\mathrm{4}}$ s⁻¹, gravitational acceleration g=9.80616 m s⁻², and thermodynamic constants R=287.05 $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ and c_p=1005.0 $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$.

To study the dispersion properties of the system of Eqs. (6)–(10), we choose horizontal wavenumber $k_x=\mathrm{2}\mathit{\pi }/{\mathrm{10}}^{\mathrm{6}}$ m⁻¹ and set the number of vertical levels n_lev=20. Dispersion and dissipation diagrams of the spacetime operators are also computed with the same n_lev and k_x to match frequencies and eigenvectors of the spacetime operators with ones from the spatial discretization.

To form a spacetime operator using Eqs. (12)–(16) and study the stability of IMEX schemes, we set n_lev=72, to emulate the default configuration of E3SM, and vary k_x throughout a representative parameter space resolvable by anticipated high-resolution models. In regimes in which stability is controlled by the CFL condition associated with acoustics modes, we desire an IMEX method whose stability will not depend on the number of vertical levels. In Sect. 4.4, we study the stability of IMEX schemes for a varying number of vertical levels.

2.2 Analytical frequencies and dispersion relations

The problem of finding frequencies ω in the system of Eqs. (6)–(10) is equivalent to investigating the spectrum of a differential operator. Since we replace the boundary condition at the top of the model, we obtain a slightly different dispersion relation for internal modes compared with previous work. Additionally, in contrast to TW2005, there are no external modes in our system.

To derive the dispersion relation from Eqs. (6)–(10), we follow Sect. 3 of TW2005 and Thuburn et al. (2002b). The dispersion relation is independent of the choice of vertical coordinate and is most easily found using the height coordinate, z. In TW2005, the hydrostatic equation, elimination, and use of the EOS yield the ODE (Eq. 57):

$$\begin{array}{}\text{(17)}& ({\partial }_z+A)({\partial }_z+B)p+C=\mathrm{0},\end{array}$$

where the constants A and B are related to the static stability and sound speed, respectively, of the isothermal reference state and C(ω) is a cubic function of the frequency ω. Expressions for A, B, and C are defined as in TW2005 (Eq. 58). As in TWS2002b, we make the change of variable $\stackrel{\mathrm{̃}}{p}=p\exp \left(\frac{(A+B)z}{\mathrm{2}}\right)$ to obtain

$$\begin{array}{}\text{(18)}& {\stackrel{\mathrm{̃}}{p}}_{zz}+a\stackrel{\mathrm{̃}}{p}=\mathrm{0},a\left(\mathit{\omega }\right)=C\left(\mathit{\omega }\right)-{\displaystyle \frac{(B-A{)}^{\mathrm{2}}}{\mathrm{4}}}.\end{array}$$

In our setting, the ODE (Eq. 18) has bottom boundary condition

$$\begin{array}{}\text{(19)}& {\stackrel{\mathrm{̃}}{p}}_z+{\displaystyle \frac{B-A}{\mathrm{2}}}\stackrel{\mathrm{̃}}{p}=\mathrm{0}\end{array}$$

at z=0 and top boundary condition $\stackrel{\mathrm{̃}}{p}=\mathrm{0}$ at z=D.

We first assume a>0. The cases a<0 and a=0 are discussed below. With $m=\sqrt{a}$ and a solution of form $\stackrel{\mathrm{̃}}{p}\left(z\right)=c_{\mathrm{1}}\sin \left(mz\right)+c_{\mathrm{2}}\cos \left(mz\right)$, we obtain the internal modes. From the top boundary condition, we recover

$$\mathrm{0}=\stackrel{\mathrm{̃}}{p}\left(D\right)=c_{\mathrm{1}}\sin \left(mD\right)+c_{\mathrm{2}}\cos \left(mD\right)\Rightarrow c_{\mathrm{2}}=-c_{\mathrm{1}}\tan \left(mD\right).$$

From the bottom boundary condition, we recover

$$\begin{array}{rl}\mathrm{0}& ={\stackrel{\mathrm{̃}}{p}}_z\left(\mathrm{0}\right)+{\displaystyle \frac{B-A}{\mathrm{2}}}\stackrel{\mathrm{̃}}{p}\left(\mathrm{0}\right)=c_{\mathrm{1}}m+c_{\mathrm{2}}{\displaystyle \frac{B-A}{\mathrm{2}}}\\ & \Rightarrow c_{\mathrm{2}}=-c_{\mathrm{1}}{\displaystyle \frac{\mathrm{2}m}{B-A}}.\end{array}$$

Combining these, we obtain a condition on wavenumber m:

$$\begin{array}{}\text{(20)}& \tan \left(mD\right)={\displaystyle \frac{\mathrm{2}m}{B-A}}.\end{array}$$

In TW2005, the internal modes obey $m=n\mathit{\pi }/D$, where n>0 is the mode number. In Eq. (20), for large m, wavenumbers are close to $\frac{n\mathit{\pi }}{D}+\frac{\mathit{\pi }}{\mathrm{2}D}$, where n is a positive integer.

Due to the nonlinearity of Eq. (20) with respect to m, wavenumbers m obeying Eq. (20) are found numerically in MATLAB by solving Eq. (20) for the first n_lev values of m_i, $i\in \mathit{\{}\mathrm{1},\mathrm{\dots },n_{\text{lev}}\mathit{\}}$, n_lev=20. We recover three wave branches (acoustic, gravity, and Rossby) by solving the quintic equation,

$$\begin{array}{}\text{(21)}& a\left(\mathit{\omega }\right)=C\left(\mathit{\omega }\right)-{\displaystyle \frac{(B-A{)}^{\mathrm{2}}}{\mathrm{4}}}=m_i^{\mathrm{2}},\end{array}$$

that follows from substituting $a=m_i^{\mathrm{2}}$ in Eq. (18) and solving for ω for each m_i. Three branches of internal waves are plotted as solid blue lines in Fig. 1.

External modes are derived assuming a<0 in Eq. (18). Solutions are then represented by $\stackrel{\mathrm{̃}}{p}\left(z\right)=c_{\mathrm{1}}{e}^{mz}+c_{\mathrm{2}}{e}^{-mz}$, $m=\sqrt{-a}$. This leads to the equation

$$\begin{array}{}\text{(22)}& \tanh \left(mD\right)={\displaystyle \frac{\mathrm{2}m}{B-A}},\end{array}$$

which does not have a solution, as follows. First, rewrite it as $\tanh \left(\stackrel{\mathrm{̃}}{m}\right)=\frac{\mathrm{2}\stackrel{\mathrm{̃}}{m}}{(B-A)D}$ with $\stackrel{\mathrm{̃}}{m}=mD$. Since $\frac{\mathrm{2}}{(B-A)D}>\mathrm{3}$, the line $\frac{\mathrm{2}\stackrel{\mathrm{̃}}{m}}{(B-A)D}$ and the curve $\tanh \stackrel{\mathrm{̃}}{m}$ do not intersect except at the origin. The origin is not a solution since we assumed a<0 and thus m≠0. Similarly, the choice a=0 cannot have solutions satisfying the boundary conditions for our particular value of D.

Analytically, one can recover external modes if the depth of the domain, D, is larger. We searched for values m in Eq. (22) for D=40 000 m and D=50 000 m and used these values in Eq. (21). We obtained five real roots with magnitudes of the order of 10⁻⁶, 10⁻³, and 10⁻², as expected. We were unable to locate external modes in discretized systems with large domain sizes. To search, we examined eigenvalues and eigenvectors using the fact that external modes have zero vertical structure in the vertical velocity coordinate.

2.3 Numerical frequencies in the HOMME discretization

To discretize the right-hand side of the systems of Eqs. (12)–(16) and (6)–(10) vertically in space, we use a Lorenz staggering and place u, v, and σ at the midpoints of the model's n_lev vertical levels, and ϕ and w at its n_lev+1 level interfaces. This staggering is denoted [wz,uvσ]. Due to the choice of boundary conditions, ϕ and w are zero at the bottom of the domain; therefore, we solve only for their n_lev interface values, excluding the bottom interface. The total vector length in the discretized system is 5n_lev.

Our placement of variables requires four discrete operators: one to interpolate ϕ from interfaces to midlevels, one to approximate the derivative ϕ_θ at midlevels, one to interpolate σ from midlevels to interfaces, and one to approximate the derivative p_θ on interfaces. Derivatives are formed using second-order finite differencing with constant level spacing Δθ. Interpolation to and from midlevels is implemented via simple averaging of neighbor values. Applying these operators at each level and interface, we can now write the discretization of the system of Eqs. (12)–(16) as the matrix equation

$$\begin{array}{}\text{(23)}& {\mathit{U}}_{\mathrm{t}}=\mathbf{M}\mathit{U}.\end{array}$$

Matrix M has a size of 5n_lev×5n_lev. The eigenvalues of M are discrete representations of quantities −iω in Eqs. (6)–(10), and the eigenvectors of M correspond to the three branches of waves: Rossby, gravity, and acoustic.

We compute the numerical eigenvalues of M with MATLAB, then match a vertical mode to each numerical eigenvalue. To find a vertical mode, we wrote a routine to count zeros in an imaginary eigenvector part that corresponds to w. For the five smallest wavenumbers, we diagnose $n=\mathrm{1}/\mathrm{3}$ manually. A few solutions for the highest wavenumbers for Rossby and gravity waves become oscillatory. Counting zeros for these is inaccurate; instead, we diagnose them using the monotonicity of numerical eigenvalues.

The numerical dispersion relation for the discretization of the system of Eqs. (6)–(10) is plotted in Fig. 1, with blue diamonds for westward-propagating waves with ω<0 and red stars for eastward-propagating waves with ω>0. As in TW2005, system [wz,uvσ], which is characterized by its staggering, choice of prognostic variables, and EOS, is in category 2b. Categories for discretizations are defined in TW2005. The optimal category is category 1; methods in this category have numerical dispersion very close to analytical. The next most optimal methods belong to the set of category 2 methods. Category 2b methods have a nearly optimal dispersion relation. Rossby frequencies are overestimated, as shown in Fig. 1, where numerical frequencies for the Rossby branch for large mode numbers are larger by absolute value (Rossby frequencies are negative) than their analytical counterparts.

Figure 1 Analytical and numerical dispersion relation for the system of Eqs. (6)–(10). Solid blue curves, from top to bottom, are acoustic, gravity, and Rossby branches of analytical solutions for ω. Blue diamonds are negative numerical eigenvalues and red stars are positive numerical frequencies ω.

3 Stability of IMEX methods from eigenvalues of a spacetime operator

In the previous section, we evaluated the properties of the spatial discretization for the system of Eqs. (6)–(10). We now combine the spatial discretization with a temporal discretization and then evaluate the resulting spacetime operator.

3.1 Spacetime operator

Similarly to Weller et al. (2013) and Lock et al. (2014), we form a spacetime operator from the system of Eqs. (12)–(16) and compute its spectrum numerically. To be stable, the eigenvalues of the spacetime operator should lie on or inside of the unit circle.

The spacetime operator is defined by the underlying IMEX scheme. Given a linear ODE

$$\begin{array}{}\text{(24)}& {\mathit{y}}_{\mathrm{t}}=\mathbf{S}\mathit{y}+\mathbf{N}\mathit{y},\end{array}$$

where S and N are the stiff and non-stiff parts already discretized in space, the spacetime operator Q can be formed from the double Butcher tableau of explicit (left) and implicit (right) tables associated with a particular IMEX scheme:

$$\begin{array}{cc}\mathit{c}& \mathbf{A}\\ & {\mathit{b}}^T\end{array}\begin{array}{cc}\widehat{\mathit{c}}& \widehat{\mathbf{A}}\\ & {\widehat{\mathit{b}}}^T\end{array}.$$

Here, A denotes the explicit matrix for an IMEX scheme and has no relation with the constant A used in Sect. 2.2. We follow the literature in using this notation.

The matrices $\mathbf{A}=\mathit{\left\{}{a}_{j_{\mathrm{1}}{j}_{\mathrm{2}}}\mathit{\right\}}$ and $\widehat{\mathbf{A}}=\mathit{\left\{}{\widehat{a}}_{j_{\mathrm{1}}{j}_{\mathrm{2}}}\mathit{\right\}}$, $j_{\mathrm{1}},j_{\mathrm{2}}=\mathrm{1},\mathrm{\dots },\mathit{\nu }$, where ν is the number of stages, and vectors c={c_j1} and $\widehat{\mathit{c}}=\mathit{\left\{}{\widehat{c}}_{j\mathrm{1}}\mathit{\right\}}$, which determine the location of internal stages, obey the constraints $c_{j_{\mathrm{1}}}={\sum }_{j_{\mathrm{2}}}{a}_{j_{\mathrm{1}}{j}_{\mathrm{2}}}$ and ${\widehat{c}}_{j_{\mathrm{1}}}={\sum }_{j_{\mathrm{2}}}{\widehat{a}}_{j_{\mathrm{1}}{j}_{\mathrm{2}}}$. The weight vectors are b={b_j2} and $\widehat{\mathit{b}}=\mathit{\left\{}{\widehat{b}}_{j\mathrm{2}}\mathit{\right\}}$. Upper-diagonal and diagonal coefficients of the explicit matrix by definition are zero: $a_{j_{\mathrm{1}}{j}_{\mathrm{2}}}=\mathrm{0}$, j₁≤j₂. We are only interested in diagonally implicit Runge–Kutta (DIRK) methods; therefore, in the implicit matrix, ${\widehat{a}}_{j_{\mathrm{1}}{j}_{\mathrm{2}}}=\mathrm{0}$ for j₁<j₂.

Later we refer to order of accuracy conditions for IMEX schemes, as defined, for example, in Rokhzadi et al. (2018). First-order conditions are

$$\begin{array}{}\text{(25)}& \sum _{j_{\mathrm{1}}}{b}_{j_{\mathrm{1}}}=\sum _{j_{\mathrm{1}}}{\widehat{b}}_{j_{\mathrm{1}}}=\mathrm{1}.\end{array}$$

Second-order conditions include the first-order conditions, conditions for each table, and the following coupling conditions for the explicit and implicit tables:

$$\begin{array}{}\text{(26)}& \sum _{j_{\mathrm{1}}}{b}_{j_{\mathrm{1}}}{c}_{j_{\mathrm{1}}}=\sum _{j_{\mathrm{1}}}{\widehat{b}}_{j_{\mathrm{1}}}{\widehat{c}}_{j_{\mathrm{1}}}=\sum _{j_{\mathrm{1}}}{\widehat{b}}_{j_{\mathrm{1}}}{c}_{j_{\mathrm{1}}}=\sum _{j_{\mathrm{1}}}{b}_{j_{\mathrm{1}}}{\widehat{c}}_{j_{\mathrm{1}}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}.\end{array}$$

For details on how to construct a spacetime operator, see Lock et al. (2014), Eq. (29) or (41), where the spacetime operator is either a scalar or a 3×3 matrix and is called an amplification factor. Here, the spacetime operator is a (5n_lev)×(5n_lev) matrix.

To form a spacetime operator from the system of Eqs. (12)–(16), we should define its stiff and non-stiff parts. For the stiff part, represented by matrix S, we consider the right-hand side terms of the equations for vertical velocity and geopotential:

$$-g{\displaystyle \frac{\mathit{\sigma }}{{\mathit{\sigma }}^{\mathrm{r}}}}-{\displaystyle \frac{p_{\mathit{\theta }}}{{\mathit{\sigma }}^{\mathrm{r}}}}\text{ and }gw.$$

It can be shown analytically (Steyer et al., 2019) or numerically, using MATLAB scripts for this project, that the eigenvalues of S coincide with frequencies for acoustic waves. All other terms on the right-hand side of Eqs. (6)–(10) contribute to the non-stiff matrix N. Unlike in the 2-D acoustic system (Lock et al., 2014), the spectrum of N does not coincide with the slow modes of the system of Eqs. (12)–(16), and N is not linear in k_x, in the sense that N≠k_xN₀ for some constant operator N₀.

3.2 Stability diagrams

To investigate the numerical stability of time-stepping schemes, it is common to refer to the 2-D acoustics system (Weller et al., 2013; Lock et al., 2014; Steyer et al., 2019):

$$\begin{array}{}\text{(27)}& {\mathit{y}}_{\mathrm{t}}=-i{k}_x\left(\begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{1}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ c_s^{\mathrm{2}}& \mathrm{0}& \mathrm{0}\end{array}\right)\mathit{y}-i{k}_z\left(\begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\\ \mathrm{0}& c_s^{\mathrm{2}}& \mathrm{0}\end{array}\right)\mathit{y},\end{array}$$

with horizontal wavenumber k_x, vertical wavenumber k_z, and $k_x=\mathrm{2}\mathit{\pi }/T_x$ and $k_z=\mathrm{2}\mathit{\pi }/T_z$ for wavelengths T_x and T_z. In this system, the spacetime operator has three eigenvalues, which we denote by λ. They are functions of C_x and C_z, $\mathit{\lambda }=\mathit{\lambda }(C_x,C_z)$, for Courant numbers C_x=c_sk_xΔt and C_z=c_sk_zΔt, where c_s is the speed of sound, and k_x and k_z are horizontal and vertical wavenumbers, respectively. The full stability diagram can then be plotted as a function of C_x and C_z or related quantities.

The relation $\mathit{\lambda }=\mathit{\lambda }(C_x,C_z)$ does not hold for the system of Eqs. (12)–(16) and its corresponding spacetime operator. In this system, the eigenvalues are functions of three parameters, $\mathit{\lambda }=\mathit{\lambda }(k_x,\mathrm{\Delta }t,n_{\text{lev}})$. For each k_x and Δt, there are 5n_nlev eigenvalues. To study the stability properties in this three-dimensional parameter space, we consider two regimes. We first set n_lev=72, to emulate the default configuration of E3SM, and vary k_x throughout a representative parameter space resolvable by anticipated high-resolution models. In the second regime, we fix k_x to the highest frequency resolvable by a model with 3 km grid spacing and consider a range of vertical levels. As we are interested in IMEX methods that treat vertical acoustic waves implicitly, an ideal method should remain stable for all choices of n_lev.

In the first regime, we plot stability diagrams with horizontal wavelength T_x on the horizontal axis and Δt on the vertical axis. We vary T_x from approximately 2 to 220 km and the time-step range from 0.5 to 400 s. For each k_x and Δt, we compose a spacetime operator that corresponds to a particular IMEX method. The operator's eigenvalues are computed numerically using one of MATLAB's solvers. The largest-magnitude eigenvalue is saved to an array that is then plotted on a stability diagram. We declare a spacetime operator stable if its largest-magnitude eigenvalue has magnitude of less than 1+ϵ_tol, with ${\mathit{ϵ}}_{\text{tol}}={\mathrm{10}}^{-\mathrm{12}}$. In our diagrams, stable regions are colored white.

3.3 Diagrams for dispersion and dissipation

Knowing the eigenpairs $(-i{\mathit{\omega }}_k,{\mathit{m}}_k)$ of the space operator M as in Eq. (23) and the eigenpairs (λ_j,q_j) of the IMEX spacetime operator Q, we recover additional properties of each IMEX scheme.

For small time steps Δt, we expect the relationship between the space operator M and the spacetime operator Q constructed for Δt step to be

$$\begin{array}{}\text{(28)}& {\mathit{q}}_j\simeq {\mathit{m}}_k\text{ and }{\mathit{\lambda }}_j=l_j{e}^{-i{\stackrel{\mathrm{̃}}{\mathit{\omega }}}_j\mathrm{\Delta }t}\end{array}$$

for some real l_j>0 and real ${\stackrel{\mathrm{̃}}{\mathit{\omega }}}_j$. We also expect each pair (q_j,m_k) to be uniquely matched.

Ideally, l_j=1 and ${\stackrel{\mathrm{̃}}{\mathit{\omega }}}_j={\mathit{\omega }}_k$; that is, there are no dissipation or dispersion errors from time stepping. In practice, we observe at least some numerical dissipation from applying an IMEX scheme, especially for acoustic waves.

To make dissipation and dispersion diagrams, we use the Munkres algorithm (Munkres, 1957) and its MATLAB implementation (Cao, 2020) to uniquely match each q_j with m_k using the cost function $-\frac{<{\mathit{m}}_k,{\mathit{q}}_j>}{\left|{\mathit{m}}_k\right|\left|{\mathit{q}}_j\right|}$, where $<\cdot ,\cdot >$ denotes an inner product. Then we examine the corresponding eigenvalue λ_j from the spacetime operator and compute its absolute value l_j and its ${\stackrel{\mathrm{̃}}{\mathit{\omega }}}_i$ from Eq. (28). We discuss dissipation and dispersion diagrams further in Sect. 4.3, where we apply them for the family of IMEX M2 methods.

4 Selectiveness of new framework

In Sect. 4.1, we provide an example of a scheme that appears to be stable for many practical choices of time steps if it is analyzed with the system of Eq. (27) but is unstable for these time steps if analyzed with the system of normal modes (Eqs. 12–16). We also apply the new framework for two schemes presented in Giraldo et al. (2013) and Rokhzadi et al. (2018).

4.1 Scheme M1

In tables (29), we present a six-stage IMEX scheme based on one of the explicit Runge–Kutta methods in Kinnmark and Gray (1984b). The explicit table (Eq. 29, top) is a low-storage, third-order method with a high CFL of $\sqrt{\mathrm{15}}\approx \mathrm{4}$. We construct the implicit table (Eq. 29, bottom) using a backward-Euler method for all implicit stages except the last one. The last stage is constructed to have three positive coefficients, including a nonzero coefficient on the main diagonal of matrix $\widehat{\mathbf{A}}$, and to obey the second-order convergence conditions for IMEX methods, Eqs. (25)–(26). The method has the same time locations of explicit and implicit internal stages and is second-order accurate. It satisfies a stiffly accurate condition; that is, the last row of $\widehat{\mathbf{A}}$ matches components of $\widehat{\mathit{b}}$.

$$\begin{array}{}\text{(29)}& \begin{array}{rl}& \begin{array}{ccccccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}/\mathrm{5}& \mathrm{1}/\mathrm{5}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}/\mathrm{5}& \mathrm{0}& \mathrm{1}/\mathrm{5}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}/\mathrm{3}& \mathrm{0}& \mathrm{0}& \mathrm{1}/\mathrm{3}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}/\mathrm{2}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}/\mathrm{2}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ & \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\end{array}\\ & \begin{array}{ccccccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}/\mathrm{5}& \mathrm{0}& \mathrm{1}/\mathrm{5}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}/\mathrm{5}& \mathrm{0}& \mathrm{0}& \mathrm{1}/\mathrm{5}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}/\mathrm{3}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}/\mathrm{3}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}/\mathrm{2}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}/\mathrm{2}& \mathrm{0}\\ \mathrm{1}& \mathrm{5}/\mathrm{18}& \mathrm{5}/\mathrm{18}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{8}/\mathrm{18}\\ & \mathrm{5}/\mathrm{18}& \mathrm{5}/\mathrm{18}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{8}/\mathrm{18}\end{array}\end{array}\end{array}$$

4.1.1 Plotting details

We plot three stability diagrams for the M1 scheme: Fig. 2a, with axes (k_xΔt,k_zΔt); Fig. 2b, with axes (T_x,Δt) for the 2-D acoustics the system of Eq. (27); and Fig. 2c, with axes (T_x,Δt) for the system of normal modes (Eqs. 12–16). In Fig. 2a, the x axis k_xΔt varies from 10⁻⁴ to 10^0.1 and the y axis k_zΔt varies from 10⁻⁴ to 10². In each dimension, we use 100 logarithmically spaced points. As in Lock et al. (2014) and Steyer et al. (2019), for each pair of k_xΔt and k_zΔt, the spacetime operator from Eq. (27) is computed using IMEX M1 method given by tables (29).

Figure 2b and c use the horizontal wavenumber k_x that corresponds to the wavelength T_x. In the figures, T_x varies from approximately 2 to 220 km, with 100 logarithmically spaced points. Since the acoustic system of Eq. (27) requires k_z, for Fig. 2b we make k_z span the range $K_{\mathrm{0}}=k_x\times [{\mathrm{10}}^{-\mathrm{2}},{\mathrm{10}}^{\mathrm{4}}]$. On the y axis, the time step varies from 0.5 to 400 s, again with 100 logarithmically spaced points. For each pair of (k_x,Δt), we compute a set of spacetime operators based on k_z∈K₀ via the same procedure as for Fig. 2a. If for each k_z the operator is stable, then point (T_x,Δt) is stable in Fig. 2b. We chose to use wavelength on the x axes of the stability diagrams, rather than wavenumber, to make it easier to identify horizontal resolutions.

Figure 2c is generated identically to Fig. 2b except that its results come from the system of Eqs. (12)–(16). Since the spatially discretized system is discretized in the vertical direction, there is no need to define k_z.

These stability diagrams are not scaled by the number of stages in the IMEX methods.

4.1.2 Stability of the M1 schemes

When using the stability diagrams in Fig. 2a, b based on the 2-D acoustic system, as in Lock et al. (2014), the scheme appears stable for reasonable time steps and resolutions as indicated by the large white (stable) regions. Both figures show that the stability of the IMEX scheme is the same as the stability of its explicit table, which is defined by Courant number S_M1≈4, as follows. In Fig. 2a, a straight vertical line going through a point k_xΔt=S_M1 remains in the white region, and in Fig. 2b the stable (white) region lies below the straight line with slope 1 that goes through the point that corresponds to values $(S_{\mathrm{M}\mathrm{1}}/(c_{\mathrm{s}}{k}_x),\mathrm{\Delta }t)$. Indeed, the approximate values T_x=2000 m, Δt=4 s, k_x=0.0031 m⁻¹, and c_s=317 m s⁻¹, which is the speed of sound of the constant reference state in the system of Eqs. (12)–(16), satisfy the last condition.

However, in Fig. 2c, time steps based solely on the stability of the explicit table in Eq. (29) are not stable. That is, the 2-D acoustic system in Eq. (27) does not have enough complexity to indicate that the method can be unstable in practice. Compared to Eq. (27), the system of normal modes contains a full set of modes: east- and west-propagating acoustic and gravity waves and westward-propagating Rossby waves. It is linearized about a non-constant hydrostatic reference state and has the commonly used constant-pressure boundary condition at the model top.

Figure 2 Example of selectiveness of the new framework using the M1 method: the stability diagram (c) based on the system of normal modes which determines scheme M1 is unstable for practical applications, while time steps that are stable in the stability diagrams (a, b) based on the 2-D acoustic system look acceptable.

4.2 Schemes ARK2(2,3,2) (Giraldo et al., 2013) and IMEX-SSP2(2,3,2) (Rokhzadi et al., 2018)

In Rokhzadi et al. (2018), one of the ARK2(2,3,2) methods from Giraldo et al. (2013) is compared with a new scheme, IMEX-SSP2(2,3,2). The family of ARK2(2,3,2) schemes is characterized by parameter a₃₂ in the explicit table (Giraldo et al., 2013). In Rokhzadi et al. (2018), the authors choose the method with $a_{\mathrm{32}}=\frac{\mathrm{1}}{\mathrm{6}}(\mathrm{3}+\mathrm{2}\sqrt{\mathrm{2}})$, which we denote here as ARK2(2,3,2)(1). Rokhzadi et al. (2018) apply optimization to derive an ARK2 method with improved accuracy, stability, and strong stability preserving (SSP) properties relative to ARK2(2,3,2)(1) for a linear wave equation, the 2-D acoustics system, the compressible Boussinesq equations, and the van der Pol equation as in Durran and Blossey (2012), Weller et al. (2013), and Lock et al. (2014). We compare these two methods and method ARK2(2,3,2) with a₃₂=0.85, which we denote ARK2(2,3,2)(2), using our system of normal modes (Eqs. 12–16). We conclude that ARK2(2,3,2)(2) and IMEX-SSP2(2,3,2) have very similar stability properties, as shown in Fig. 3b and c, but the stable (white) region for ARK2(2,3,2)(1) is significantly smaller, as shown in Fig. 3a.

Figure 3 Schemes analyzed in Rokhzadi et al. (2018) and the ARK2(2,3,2)(2) scheme: (a) ARK2(2,3,2)(1), (b) ARK2(2,3,2)(2), and (c) IMEX-SSP2(2,3,2).

4.3 A set of low-storage, high-CFL IMEX schemes

We develop a set of methods we denote M2 using a second-order, explicit, low-storage, CFL-4 Runge–Kutta scheme from Kinnmark and Gray (1984a). Low-storage, high-CFL methods developed in Kinnmark and Gray (1984a) and Kinnmark and Gray (1984b) are used in HOMME, the nonhydrostatic atmospheric dynamical core of E3SM's atmosphere component. It is practical to extend existing explicit RK schemes to IMEX RK methods. We analyze the stability, dispersion, and dissipation properties of these M2 schemes using the system of normal modes (Eqs. 12–16).

4.3.1 Definitions

We start with the second-order explicit table (Eq. 30, top) from one of the methods in Kinnmark and Gray (1984a). For the implicit table (Eq. 30, bottom), we choose the same times for internal stages and make all but the last implicit stage backward Euler. Internal backward-Euler stages provide stability and do not affect the second-order accuracy conditions for IMEX given by Eqs. (25)–(26).

The M2 methods vary only in their last implicit stage. We require the last implicit stage to obey a stiffly accurate condition and have only non-negative entries in its table. The last stage is defined by the vector d, whose entries correspond to the last row of the implicit Butcher tableau, as shown in the table of Eq. (30), bottom. Moreover, here we consider only schemes with at most three nonzero entries in d to limit storage. In practice, using an IMEX method in a 3-D model with topography will require storing geopotential and vertical velocity terms for each internal stage that corresponds to $d_{j_{\mathrm{1}}}\ne \mathrm{0}$, j₁<ν. Therefore, we focus on methods that limit such storage space.

$$\begin{array}{}\text{(30)}& \begin{array}{rl}& \begin{array}{ccccccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}/\mathrm{4}& \mathrm{1}/\mathrm{4}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}/\mathrm{6}& \mathrm{0}& \mathrm{1}/\mathrm{6}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{3}/\mathrm{8}& \mathrm{0}& \mathrm{0}& \mathrm{3}/\mathrm{8}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}/\mathrm{2}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}/\mathrm{2}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ & \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\end{array}\\ & \begin{array}{ccccccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}/\mathrm{4}& \mathrm{0}& \mathrm{1}/\mathrm{4}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}/\mathrm{6}& \mathrm{0}& \mathrm{0}& \mathrm{1}/\mathrm{6}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{3}/\mathrm{8}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{3}/\mathrm{8}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}/\mathrm{2}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}/\mathrm{2}& \mathrm{0}\\ \mathrm{1}& d_{\mathrm{1}}& d_{\mathrm{2}}& d_{\mathrm{3}}& d_{\mathrm{4}}& d_{\mathrm{5}}& d_{\mathrm{6}}\\ & d_{\mathrm{1}}& d_{\mathrm{2}}& d_{\mathrm{3}}& d_{\mathrm{4}}& d_{\mathrm{5}}& d_{\mathrm{6}}\end{array}\end{array}\end{array}$$

where

$$(d_{\mathrm{1}},d_{\mathrm{2}},d_{\mathrm{3}},d_{\mathrm{4}},d_{\mathrm{5}},d_{\mathrm{6}})=\mathit{d}.$$

We consider the first- and second-order variants of the M2 methods listed in Table 1.

Table 1 The M2 schemes.

Variants M2a, M2b, and M2c, where M2c was introduced in Steyer et al. (2019), are second-order methods with good stability properties; their dispersive and dissipative characteristics are different, as shown below in Fig. 4. Variants M2be and M2cn are the two extremes of the M2 family. In M2be, the last stage is the backward-Euler method. Thus, it is first-order accurate, and we expect this scheme to be the most stable but also the most dissipative. We expect method M2cn, which has the Crank–Nicolson method for the last stage, to have no dissipation for hyperbolic problems like ours. We also analyze method M2cno, whose last stage is the Crank–Nicolson method with off-centering, since off-centering is a common approach to stabilize time-stepping schemes (Durran and Blossey, 2012; Staniforth et al., 2006).

4.3.2 Stability diagrams and dispersion and dissipation diagrams

Stability diagrams for the M2 schemes are shown in Figs. 4a–c and 5a–c, which follow the plotting procedures described in Sect. 4.1.1. We plot numerical frequencies of the space operator Q and the spacetime operator M to evaluate how the numerical time-stepping methods preserve the frequencies ω from the space discretization Q. Hyperbolicity of the system of Eqs. (12)–(16) implies exact time integration will conserve the frequencies. Inexact time integration will introduce errors, which we evaluate below. We also evaluate numerical damping since exact time integration has none. Note that we compare the properties of the spacetime operator with those of the space operator integrated exactly in time, but we do not compare solutions of the spacetime operator with analytical solutions of the system of Eqs. (12)–(16). This choice of comparison focuses on the numerical errors due solely to the time-stepping methods.

Dispersion and dissipation plots are shown below the stability diagrams, for the spacetime operator Q with eigenpairs (λ_j,q_j), in Figs. 4d–f and 5d–f. In each figure, Δt=50 s and n_lev=20. The top plots show numerical frequencies λ_j vs. vertical mode number. Red diamonds are numerical frequencies of the space operator for east- and west-propagating acoustic waves. Blue squares represent east- and west-propagating gravity waves for the space operator. Black diamonds are frequencies for west-propagating Rossby waves for the space operator. Red stars, blue plus signs, and black stars are for corresponding branches of the spacetime operator. Vertical mode number and wave characterization are obtained by uniquely matching eigenvector q_j to its counterpart, eigenvector m_k, of the space operator M, also computed using 20 vertical levels.

The bottom plots in Figs. 4d–f and 5d–f show the amplification factors of eigenvalues, $\left|{\mathit{\lambda }}_j\right|$, for the spacetime operator. The red stars, blue plus signs, and black stars are for acoustic, gravity, and Rossby waves, respectively. Each plot shows amplification factors near 1 for gravity and Rossby waves, with additional damping of the acoustic modes. We discuss these differences further in the next section.

Figure 4 Properties of M2 schemes: stability, dispersion, and dissipation diagrams.

Figure 5 Properties of M2 schemes: stability, dispersion, and dissipation diagrams.

4.3.3 Analyzing the M2 schemes

Due to their different final stages, the M2 schemes have different stability properties and dispersive and dissipative characteristics. To evaluate stability, we focus on the regions having smallest spatial resolution and highest wavenumbers, since these are the regions where nonhydrostatic effects are most prominent. Evaluation of stability is easy: larger stable (white) regions imply larger stable Δt for those methods.

As expected, due to the presence of the last backward-Euler stage in the implicit table, the M2be scheme is the most stable. Recall that analytically for hyperbolic problems, the backward-Euler method is unconditionally stable and is very dissipative. For the M2be method, the largest stable Δt at T_x=2 km is approximately 4 s, which is at least twice as large as the largest stable Δt for the other schemes.

It is desirable to have an IMEX method with stability properties similar to those of an explicit method used for the non-stiff part of the system of Eq. (24). In other words, it is desirable to be able to integrate a nonhydrostatic system using an IMEX method with a time step as large as the time step used to integrate a hydrostatic system using an explicit Runge–Kutta method. Therefore, we compare the stability of the IMEX method M2be with the stability of the Runge–Kutta method consisting of the explicit table in M2be, a method we denote MExplicit, when it is applied to the non-stiff part of Eq. (24). The stability region of M2be in Fig. 5a is almost as big as the stability region of MExplicit (not shown here) up to a minor difference at approximate wavelength T_x=220 km. That is, the stability region of M2be is the biggest region we could possibly get from an IMEX scheme whose explicit table is part of the M2 set.

It is harder to rank schemes using dispersion and dissipation diagrams. All schemes preserve the dispersion and dissipation relations for gravity and Rossby waves to a high degree. But they perform very differently for acoustic waves. Method M2be has the largest dissipation for acoustic waves and is the only scheme that does not have regions of negative group velocity for acoustic waves. Method M2cn has opposite characteristics: it does not dampen acoustic waves, and the errors in the acoustic frequencies are much larger. M2cno, a first-order variation of M2cn, has dispersion errors that are very similar to M2cn while introducing low-degree dissipation into acoustic waves.

Since acoustic waves can be considered insignificant for atmospheric applications due to their low energy, one is tempted to discard numerical errors in the dispersion and dissipation of acoustic waves. However, there is an argument (Thuburn, 2012) that correct representation of even energetically weak waves in the atmosphere is crucial for accurate restoration of hydrostatic balance.

Among other second-order schemes (M2a, M2b, and M2c), it is hard to declare a clear winner. Due to its smaller largest stable Δt at T_x=2 km and large dispersion errors, M2a may be less competitive. Compared with M2c, M2b has slightly larger maximum stable Δt at T_x=2 km, and its errors in dispersion for acoustic waves are smaller, but its dissipation is larger. Indeed, its larger dissipation is probably what causes its better stability compared with M2c. However, depending on the evaluation criteria, M2c can be viewed as a better scheme than M2b. For example, it has less dissipation and its dispersion is very similar to that of the Crank–Nicolson method, which is widely used for hyperbolic problems.

4.3.4 Role of the implicit table

We chose to limit our search for a suitable M2 method by varying only the vector d in the implicit table. Since the second-order accuracy conditions for the M2 family depend only on the last implicit stage, we make all other implicit stages (stages 2–5) use the backward-Euler method to presumably maximize stability. One could also try to use Crank–Nicolson or off-centered Crank–Nicolson methods for implicit stages 2–5.

To understand how the implicit stages influence dissipation of acoustic waves, we consider the expression for the final solution of Eq. (24) using Lock et al. (2014), Eq. (15), and the definition of the M2 family in Eq. (30):

$$\begin{array}{rl}{\mathit{y}}^{n+\mathrm{1}}& ={\mathit{y}}^n+\mathrm{\Delta }t\mathbf{N}({\mathit{y}}^{\left(\mathrm{5}\right)},t^n+\mathrm{\Delta }t)+d_{\mathrm{1}}\mathrm{\Delta }t\mathbf{S}\left({\mathit{y}}^n,t^n\right)\\ & +d_{\mathrm{2}}\mathrm{\Delta }t\mathbf{S}\left({\mathit{y}}^{\left(\mathrm{2}\right)},t^n+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\mathrm{\Delta }t\right)\\ & +d_{\mathrm{3}}\mathrm{\Delta }t\mathbf{S}\left({\mathit{y}}^{\left(\mathrm{3}\right)},t^n+{\displaystyle \frac{\mathrm{1}}{\mathrm{6}}}\mathrm{\Delta }t\right)\\ & +d_{\mathrm{4}}\mathrm{\Delta }t\mathbf{S}\left({\mathit{y}}^{\left(\mathrm{4}\right)},t^n+{\displaystyle \frac{\mathrm{3}}{\mathrm{8}}}\mathrm{\Delta }t\right)\\ & +d_{\mathrm{5}}\mathrm{\Delta }t\mathbf{S}\left({\mathit{y}}^{\left(\mathrm{5}\right)},t^n+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\Delta }t\right)\\ & +d_{\mathrm{6}}\mathrm{\Delta }t\mathbf{S}\left({\mathit{y}}^{\left(\mathrm{6}\right)},t^n+\mathrm{\Delta }t\right).\end{array}$$

Scheme M2cn, given by its final implicit stage $(d_{\mathrm{1}},\mathrm{\dots },d_{\mathrm{6}})=(\mathrm{1}/\mathrm{2},\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{1}/\mathrm{2})$, does not have dissipation. For M2cn, the solution yⁿ⁺¹ is directly influenced by intermediate implicit stages 2–5 only by the non-stiff term. Also, its final implicit stage is represented by the Crank–Nicolson method, known to be non-dissipative for hyperbolic problems. We conclude that both of these facts contribute to the lack of dissipation in M2cn. Scheme M2be has final implicit stage $(d_{\mathrm{1}},\mathrm{\dots },d_{\mathrm{6}})=(\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{1})$, which gives the backward-Euler method. Similarly to the backward-Euler method for hyperbolic problems, M2be is very dissipative.

We suggest that the dispersion and dissipation of acoustics waves can be tuned by working only with the implicit table of any method.

Figure 6 Stability diagrams with respect to varying n_lev.

4.4 Stability properties with respect to vertical resolution

For an explicit time-stepping method, the most restrictive CFL condition is usually that associated with the vertically propagating acoustic waves and the stable time step decreases linearly with $\mathrm{\Delta }z=D/n_{\text{lev}}$. Ideally, with an implicit treatment of vertical acoustic waves, an IMEX method should remain stable as n_lev is increased, and the stability should be controlled only by the CFL condition associated with the horizontal resolution.

To analyze this aspect of various IMEX methods, we fix k_x to the highest frequency resolvable by a model with 3 km grid spacing and vary the number of vertical levels from n_lev=20 to n_lev=100. We plot the method's stability as a function of n_lev using a logarithmic scale (up to rounding to the nearest integer) and 50 sample points. The stability diagrams are made very similarly to the ones in Fig. 4; the only difference is the horizontal axis, which now is n_lev. We vary Δt from 1 to 10 s with logarithmic spacing and 100 samples. Note that the horizontal axis is not defined by a vertical wavenumber, k_z, because for any fixed resolution Δz the model supports waves with many vertical wavenumbers.

Figure 6 contains stability diagrams for schemes M1, ARK2(2,3,2)(1), ARK2(2,3,2)(2), and M2b. For schemes M1 and ARK2(2,3,2)(1), the stability is independent of Δz, as desired, only for up to approximately n_lev=57 (Δz≃175 m). In Fig. 6a–b, the stable region for the approximate interval $n_{\text{lev}}\in [\mathrm{20},\mathrm{57}]$ ($\mathrm{\Delta }z\in [\mathrm{175},\mathrm{500}]$ m) is under a straight line for some Δt=Δt₀. For finer Δz, the stable regions lie below a line with a constant slope for both schemes.

In contrast, for methods ARK2(2,3,2)(2) and M2b, stability is always controlled by the horizontal resolution: in Fig. 6c, the stable region is below the horizontal line Δt₀≃7.2 s. To further support this conclusion, we also computed eigenvalues of the spacetime operator for method M2b, Δt=7 s, and a few large values of n_lev, up to 600. The spacetime operator for all large n_lev was stable.

We do not present stability diagrams for Δz studies for other methods from this paper because they are identical to Fig. 6c up to the value of Δt₀. That is, the stability of methods IMEX-SSP2(2,3,2) and M2 methods is controlled by the horizontal wavelengths.

5 Conclusions

We developed a new framework to evaluate IMEX RK methods for atmospheric modeling. The framework uses a system of normal modes and is proven to be simple but more selective than the 2-D acoustics system used in the literature. For example, the M1 method from Sect. 4.1 appears to be stable for a large set of time steps and resolutions when using the 2-D acoustics system. If the method is evaluated with the system of normal modes, it is unstable for the same set of time steps and resolutions.

The new framework gives us insight to develop a set of second-order, low-storage, high-CFL IMEX RK methods to use in atmospheric dynamical cores. Furthermore, we use the spacetime operator built with the system of normal modes to investigate dispersion and dissipation of IMEX RK schemes for three types of waves: gravity, Rossby, and acoustic.

One extension of this work would be to investigate selectiveness of the framework based not on the system of normal modes (Eqs. 12–16) but on a system of compressible Boussinesq equations as in Durran and Blossey (2012).