A framework to evaluate IMEX schemes for atmospheric models

We present a new evaluation framework for implicit and explicit (IMEX) Runge-Kutta timestepping schemes. The new framework uses a linearized nonhydrostatic system of normal modes. We utilize the framework to investigate stability of IMEX methods and their dispersion and dissipation for 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 2nd order low-storage IMEX Runge-Kutta methods with high CFL. We show that the new framework is more selective than the 2D acoustic system previously used in literature. 5 Schemes that are stable for the 2D acoustic system are not stable for the system of normal modes.

properties. In Section 4.4 we investigate convergence of IMEX methods with respect to vertical resolution. Section 5 concludes.

Linearized system for normal modes
In this section we define the linearized system of equations that corresponds to the HOMME nonhydrostatic dynamics model.
We recover analytical and numerical frequencies for spacial discretization of the system and confirm that the discretization (which we broadly define here to include choice of prognostic variables, equation of state, boundary conditions, and staggering) 40 is nearly optimal. Therefore, later the system can be used to investigate properties of IMEX spacetime operators.

Description of the system
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ṡ = 0, following Lin (2004). In practice, the vertical levels must be remapped to a reference configuration at regular time intervals; we do not consider the effects of these vertical remap 45 operations here. After linearization and substituting single mode solutions in which each field is proportional to exp(ikx + ily − iωt), this formulation is equivalent to the system of equations (20)-(24) in isentropic coordinate from TW2005. With inclusion of the β-effect as in equations (55)-(56) of TW2005, this system is as follows: with linearized equation of state (EOS) 2 https://doi.org/10.5194/gmd-2020-178 Preprint. Discussion started: 23 July 2020 c Author(s) 2020. CC BY 4.0 License.
We also retain a version with time derivatives: Here u, v, and w are velocity components, p is pressure, ρ is density, φ is geopotential, k x and l x are horizontal wavenumbers with K 2 = k 2 x + l 2 x , g is the gravity constant, κ = R/c p is a thermodynamic constant, σ 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 60 variables defined by reference profiles of a linearized hydrostatic steady state with constant temperature T 0 . The subscript θ denotes partial differentiation with respect to potential temperature. Other variables are first-order perturbed quantities, about the reference state, as follows from linear analysis.
We note that sinceθ = 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 systems (1)-65 (5) and (6)-(7) as in TW2005, but a different top BC. We replace the rigid lid boundary condition with a constant pressure boundary condition, which for perturbed pressure variable becomes p top = 0. This BC is more typical for a mass-based vertical coordinate.
To study dispersion properties of system (1)-(5) we choose horizontal wavenumber k x = 2π/10 6 m −1 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 spacial discretization.
To form a spacetime operator using Eqs. (6)-(7) and study the stability of IMEX schemes, we set n lev = 72, to emulate the 75 default configuration of E3SM, and vary k x throughout a representative parameter space resolvable by anticipated high resolution models. In regimes where stability is controlled by Courant-Friedrichs-Lewy (CFL) condition associated with acoustics modes, we desire an IMEX method where the stability will not depend on the number of vertical levels. In Sect. 4.4 we study the stability of IMEX schemes for varying number of vertical levels.

Analytic frequencies and dispersion relations 80
The problem of finding frequencies ω in system (1)-(5) is equivalent to investigating a spectrum of an ODE. Since we replace boundary conditions at the top of the model, we obtain slightly different dispersion relation for internal modes compared to previous work. And, in contrast with TW2005, there are no external modes in our system.
To derive the dispersion relation, 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. The hydrostatic 85 equation, elimination, and use of the EOS yield the ODẼ is a change of variable, 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 ω. A, B, and C are defined as in TW2005, equation (58). In our setting, the ODE has bottom boundary condition at z = 0 and top boundary conditionp = 0 at z = D.
We first assume a > 0. Cases a < 0, a = 0 are discussed below. With m = √ a and solution of formp(z) = c 1 sin(mz) + c 2 cos(mz), we recover internal modes. Instead of internal modes with vertical wavenumber m = nπ/D, where n > 0 is the mode number, as in TW2005, we obtain solutions with wavenumber m obeying 95 tan(mD) = 2m For large m, these wavenumbers are close to nπ D + π 2D , n is a positive integer. Wavenumbers are found numerically in Matlab by solving Eq. (10) for m i ∈ {1, . . . , n lev }, n lev = 20. We recover three wave branches, acoustic, gravity, and Rossby, by solving the quintic equation for D = 40000 m and D = 50000 m and used these values in Eq. (11). We obtained five real roots with magnitudes of order 110 10 −6 , 10 −3 , and 10 −2 , as expected. We were unable to locate external modes in discretized systems with large domain sizes.
To search, we examined eigenvalues and eigenvectors (external modes have zero vertical structure in vertical velocity).

Numerical frequencies for HOMME discretization
To discretize system (6)-(7) 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 its n lev + 1 level interfaces. This staggering is denoted [wz, uvσ]. Due to the choice 115 4 https://doi.org/10.5194/gmd-2020-178 Preprint. Discussion started: 23 July 2020 c Author(s) 2020. CC BY 4.0 License. of boundary conditions, φ and w are zero at the bottom of the domain, therefore, we solve only for their n lev interface values, excluding 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 120 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 system (6)-(7) as a matrix equation, Matrix M is of size 5n lev × 5n lev . The eigenvalues of M are −iω and its eigenvectors correspond to three branches of waves, Rossby, gravity or acoustic.

125
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 = 1/3 manually. A few solutions for highest wavenumbers for Rossby and gravity waves become oscillatory and counting zeros for them is inaccurate. We diagnose them using monotonicity of numerical eigenvalues.
The numerical dispersion relation for the discretization of system (1)-(5) is plotted on Fig. 1 with blue diamonds for west-130 ward 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. This category has a near optimal dispersion relation with overestimated Rossby 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 system (1)-(5) . We now combine the spatial 135 discretization with a temporal discretization and then evaluate the resulting spacetime operator.

Spacetime operator
Similarly to Weller et al. (2013) and Lock et al. (2014), we form a spacetime operator from system (6)-(7) 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 where S and N are the stiff and nonstiff 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, Here A denotes the explicit matrix for an IMEX scheme and has no relation with constant A used in Sect. 2.2. We keep both notations as is due to their use in literature.
Later we refer to order of accuracy conditions for IMEX schemes, as defined, for example, in Rokhzadi et al. (2018). First- 150 Second-order conditions include the 1st-order conditions and contain conditions for each table and coupling conditions for explicit and implicit tables, 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 3x3 matrix and is called an amplification factor. Here, the spacetime operator is a (5n lev ) × (5n lev ) matrix.
To form a spacetime operator from system (6)- (7) we should define its stiff and nonstiff parts. For the stiff part, represented by matrix S, we consider the right-hand side terms of the equations for vertical velocity and geopotential, It can be shown analytically (Steyer et al., 2019) or numerically, using Matlab scripts for this project, that eigenvalues of S coincide with frequencies for acoustic waves. All other terms in the right-hand side of Eqs.
(1)-(5) contribute to the nonstiff matrix N. Unlike in the 2D acoustic system (Lock et al., 2014), the spectrum of N does not coincide with the slow modes of system (6)- (7) and N is not linear in k x , in the sense that N = k x N 0 for some constant operator N 0 .

Stability diagrams 160
To investigate the numerical stability of timestepping schemes it is common to refer to the 2D acoustics system (Weller et al., 2013;Lock et al., 2014;Steyer et al., 2019) with horizontal wave number k x , vertical wave number k z , k x = 2π/T x and k z = 2π/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 , λ = λ(C x , C z ), 165 for Courant numbers C x = c s k x ∆t and C z = c s k 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 λ = λ(C x , C z ) does not hold for system (6)- (7) and its corresponding spacetime operator. In this system, the eigenvalues are functions of three parameters, λ = λ(k x , ∆t, n 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 170 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 175 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 (by magnitude) eigenvalue is saved to an array which is then plotted on a stability diagram.
We declare a spacetime operator stable if its maximum eigenvalue is less than 1+ tol , with tol = 10 −12 . In our diagrams, stable regions are colored white.

Diagrams for dispersion and dissipation
Knowing the eigenpairs (−iω k , m k ) of the space operator M as in Eq. (13) and the eigenpairs (λ j , q j ) of the IMEX spacetime operator Q we recover additional properties of each IMEX scheme.
For small timesteps ∆t we expect the relationship between the space operator M and the spacetime operator Q constructed for ∆t step to be 185 q j m k and λ j = l j e −iωj ∆t (18) for some real l j > 0 and complexω j . We also expect each pair (q j , m k ) to be uniquely matched.
Ideally, l j = 1 andω j = ω k , that is, there are no dissipation or dispersion errors from timestepping. In practice, we observe at least some numerical dissipation from applying IMEX, especially for acoustic waves.
To make dissipation/dispersion diagrams, we use the Munkres algorithm (Munkres, 1957) and its Matlab implementation 190 (Cao, 2020 (accessed March 22, 2020) to uniquely match each q j with m k using the cost function − <m k ,qj > |m k ||qj | , where < ·, · > denotes an inner product. Then we examine corresponding eigenvalue λ j from the spacetime operator and compute its absolute value l j and itsω i from Eq. (18). More on dissipation/dispersion diagrams is in Sect. 4.3 where we apply them for a family of IMEX methods M2.

195
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 analysed with system (17) but is unstable for these timesteps if analysed with system of normal modes (6)-(7) . We also apply the new framework for two schemes presented in Giraldo et al. (2013) and Rokhzadi et al. (2018).

Plotting details
We plot three stability diagrams for the M1 scheme: Fig. 2(a) with axes (k x ∆t, k z ∆t) and Fig. 2(b) with (T x , ∆t) axes for 2D acoustics system (17) and Fig. 2(c) with (T x , ∆t) axes for system of normal modes (6)-(7) . In Fig. 2(a)  logarithmic spacing for 100 sampling points. For each pair of (k x , ∆t), we compute a set of spacetime operators based on k z ∈ K 0 via the same procedure as for Fig. 2(a). If for each k z operator is stable, then point (T x , ∆t) is stable in Fig. 2(b). We chose to plot T x wavelengths on x-axes of stability diagrams instead of wavenumbers to make it easier to identify horizontal resolutions.
Figure 2(c) is generated identically to Fig. 2(b) except that its results come from the system (6)-(7) . Since its spatially 220 discretized version is discretized in vertical direction, there is no need to define k z .

Stability of M1
When using stability diagrams in Figs. 2(a,b) based on the 2D acoustic system, as in Lock et al. (2014), the scheme appears stable for reasonable timesteps and resolutions as indicated by the large white (stable) regions. Both figures show that stability of the IMEX scheme is the same as the stability of its explicit table, which is defined by Courant number S M 1 ≈ 4, as follows.

Definitions
We start with the 2nd order explicit table, (20), left, from one of methods in Kinnmark, I. and Gray, W. (1984a). For the implicit topography will require storing geopotential and vertical velocity terms for each internal stage that corresponds to d j1 = 0, j 1 < ν. Therefore, we focus on methods that limit such storage space. where We consider the 1st-and 2nd-order variants of M2 methods listed in Table 1.
Variants M2a, M2b, M2c (M2c is introduced in Steyer et al. (2019)) are 2nd order methods with good stability properties; 255 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, so the scheme is expected to be the most stable but also the most dissipative method as it is 1st order accurate. Method M2cn, with Crank-Nicolson for the last stage, presumably has no dissipation for hyperbolic problems like ours. We also analyze method M2cno, where the last stage is Crank-Nicolson with off-centering, since off-centering is a common practice to stabilize timestepping schemes (Durran and 260 Blossey, 2012;Staniforth et al., 2006).

Stability diagrams and dispersion/dissipation diagrams
Stability diagrams of the M2 schemes are shown in Figs. 4(a-c) and 5(a-c) which used plotting procedures described in Sect.
4.1.1. There, we plot numerical frequencies of space operator Q and spacetime operator M to evaluate how numerical timestepping methods, IMEX, preserve frequencies ω from space discretization Q. In other words, due to hyperbolicity of our system (6)-(7) exact time integration would conserve the frequencies. Inexact IMEX time integration will introduce errors, which we evaluate below. We also evaluate numerical damping introduced by IMEX since exact time integration does not introduce damping. Note that we compare properties of the spacetime operator to properties of space operator integrated exactly in time, we do not compare solutions of the spacetime operator to analytical solutions of system (6)- (7) . This is because we want to investigate the numerical errors due solely to the timestepping methods. There, red stars, blue plus signs, and black stars are for acoustic, gravity, and Rossby waves correspondingly. Each plot shows amplification factors near 1 for gravity and Rossby waves, with additional damping of the acoustic modes. We discuss these 280 differences further in the next section.

Analysing 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 regions of smallest spatial resolutions and highest wavenumbers, since those are regions where nonhydrostatic effects are most prominent. Thus stability evaluation is easy: bigger stable (white) regions translate to 285 larger stable ∆t for those methods.
As expected, M2be scheme is the most stable. Its largest stable ∆t at T x = 2 km is approximately 4 sec, which is at least 2× larger than the largest stable ∆t for the other schemes. Its stability region in Fig. 5(a) coincides with stability region of its explicit table (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.

290
It is harder to rank schemes using dispersion/dissipation diagrams. All schemes preserve dispersion and dissipation relations for gravity and Rossby waves to a high degree. They perform very differently for acoustic waves. Method M2be has the biggest dissipation rates for acoustics waves and is the only scheme that does not have regions of negative group velocity for acoustic waves. Method M2cn is its opposite: it has no damping of acoustic waves while errors in acoustic frequencies are much larger.
M2cno, a 1st order variation of M2cn, has dispersion errors very similar to M2cn while introducing low-degree dissipation into 295 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 atmosphere is crucial for restoration of hydrostatic balance.  Among other 2nd order schemes, M2a, M2b and M2c, it is hard to declare a clear winner. Due to its smaller largest stable 300 ∆t at T x = 2 km and big dispersion errors, M2a may be less competitive. Comparing M2b and M2c, M2b has slightly larger maximum stable ∆t at T x = 2 km, its errors in dispersion for acoustics waves are smaller, but its dissipation rates are larger.
Indeed, its dissipative rates are probably what cause its better stability compared to M2c. However, depending on evaluation criteria, M2c can be viewed as a better scheme than M2b. For example, it has smaller dissipation rates and its dispersion is very similar to Crank-Nicolson's method, which is widely used for hyperbolic problems.  (20), right. In this case it is standard to form a function s(z), where z is complex and often chosen to be purely imaginary due to strong hyperbolicity of 310 systems of atmospheric dynamics. Here, since 2nd order accuracy conditions for methods M2 depend only on the last implicit stage, we make all other implicit stages backward Euler to presumably maximize stability.
As a next step, from observing the stability and dispersion/dissipation diagrams in Figs. 4, 5 we speculate that stability of any M2 method is directly related to amount of dissipation provided by the last stage coefficients d. That is, choices where |s(z)| is smaller would lead to bigger stability regions in stability diagrams.

315
In the M2 methods, dispersion and dissipation of gravity and Rossby waves do not seem to be affected by the implicit table, in particular, they seem insensitive to the most dissipative backward Euler stages 2 to 5. On the contrary, acoustic waves are affected by the implicit table and its last, 6th, implicit stage. We make a suggestion that dispersion and dissipation of acoustics waves can be tuned only by working with the implicit table of any method. For an explicit timestepping method, the most restrictive CFL condition is usually that associated with the vertically propagating acoustic waves and the stable timestep would decrease linearly with ∆z = D/n 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 325 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 sampling points. Stability diagrams are made very similarly to the ones in Fig. 4, with only difference of the horizontal axis, which now defines n lev . We vary ∆t from 1 to 10 sec 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. In contrast, for method M2b, stability is always controlled by the horizontal resolution: in Fig. 6(c) the stable region is below 335 horizontal line ∆t 0 7.2 sec. To further support this conclusion, we also computed eigenvalues of the spacetime operator for method M1, ∆t = 7 sec, 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.   6(c) up to the value of ∆t 0 . That is, stability of methods IMEX-SSP2(2,3,2) and M2 methods is controlled by the horizontal wavelengths.
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 2D acoustics system used in literature. For example, the M1 method from Sect. 4.1 appears to be stable for a large set of timesteps and resolutions when using the 2D acoustics system.
If the method is evaluated with the system of normal modes, it is unstable for the same set of timesteps and resolutions.

345
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 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 (6)-(7) but on a system of compressible Boussinesq equations as in Durran and Blossey (2012). March 22, 2020).
For this submission, we created script paper_figures.m, also archived on Zenodo Steyer, 2020 (accessed July 06, 2020), which sets parameters and launches Matlab scripts to produce all paper figures in order they appear. The script contains comments to easily identify the figures.
Sta-imex version 1.0: Copyright 2020 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of