## Abstract

Steep orography can cause noisy solutions and instability in models of the atmosphere. A new technique for modeling flow over orography is introduced that guarantees curl-free gradients on arbitrary grids, implying that the pressure gradient term is not a spurious source of vorticity. This mimetic property leads to better hydrostatic balance and better energy conservation on test cases using terrain-following grids. Curl-free gradients are achieved by using the covariant components of velocity over orography rather than the usual horizontal and vertical components.

In addition, gravity and acoustic waves are treated implicitly without the need for mean and perturbation variables or a hydrostatic reference profile. This enables a straightforward description of the implicit treatment of gravity waves.

Results are presented of a resting atmosphere over orography and the curl-free pressure gradient formulation is advantageous. Results of gravity waves over orography are insensitive to the placement of terrain-following layers. The model with implicit gravity waves is stable in strongly stratified conditions, with *N*Δ*t* up to at least 10 (where *N* is the Brunt–Väisälä frequency). A warm bubble rising over orography is simulated and the curl-free pressure gradient formulation gives much more accurate results for this test case than a model without this mimetic property.

## 1. Introduction

As the resolution of atmospheric models increases, the orography resolved becomes steeper, which leads to pressure gradient errors (Gary 1973), which can lead to noisy solutions (e.g., Hoinka and Zangl 2004) or even instability (e.g., Webster et al. 2003). A variety of techniques for avoiding this problem have been proposed, which will be discussed. However, none of them solves the problem that existing discretizations of the pressure gradient over orography are not curl free. This means that pressure gradients can be spurious sources of vorticity, which may lead to noisy vorticity fields away from the surface such as that reported by, for example, Hoinka and Zangl (2004).

While resolution is increasing, it is still necessary to create models that can run stably with long time steps in the presence of high stratification. This means that gravity waves, as well as acoustic waves, should be treated implicitly (at least in the vertical direction, in which resolution is higher). A variety of methods for treating gravity waves implicitly have been described (e.g., Cullen 1990; Smolarkiewicz et al. 2014) that involve separating atmospheric variables into mean and perturbation quantities and linearizing. These will be discussed, which will motivate an alternative approach that does not rely on an explicit linearization.

The introduction of orography into atmosphere models is usually done using terrain-following coordinates, so that the grid does not intersect with the ground, grid boxes are arranged exactly in vertical columns, and high resolution of the planetary boundary layer is maintained (e.g., Schär et al. 2002; White 2003; Melvin et al. 2010). Whether the equations in the transformed coordinates are discretized on a uniform grid, or the equations in Cartesian coordinates are discretized on a curvilinear terrain-following grid defined by the terrain-following coordinates, existing models do not have curl-free discretizations of the gradient operator, which is likely to lead to problems over steep orography. There are a variety of approaches to alleviating this problem that will be discussed.

Smoothing of orography has been used in order to avoid noisy solutions and instabilities associated with steep orography (e.g., Kanamitsu et al. 2002; Webster et al. 2003) but smoothed orography can lead to problems such as reduced barrier heights and raised sea levels (Rutt et al. 2006) or elevated heat sources (Kanamitsu et al. 2002). A popular alternative is to use terrain-following coordinates (or layers), which rapidly become smooth with height (e.g., Schär et al. 2002; Klemp 2011) so that the pressure gradient errors are reduced away from the ground. Hoinka and Zangl (2004) found that this approach avoided the spurious potential vorticity (pv) fields near the tropopause over steep orography in the fifth-generation Pennsylvania State University–National Center for Atmospheric Research Mesoscale Model (MM5). However, this smooth-layers approach leads to very thin model layers over mountain peaks, which can lead to instability, while layers adjacent to the mountain slopes will not be smooth and so will still have large numerical errors that can be detrimental for predictions of mountain weather (Fast 2003).

A complementary approach is to improve the accuracy of the pressure gradient calculation. In atmospheric models, the prognostic velocity variables are usually the vertical velocity and two components of horizontal velocity. To solve the components of the momentum equation, the pressure gradient is needed in the same direction as the velocity components. This is straightforward for the vertical velocity because the prognostic pressure variables will also be aligned in vertical columns and so the vertical pressure gradient can be accurately calculated in a straightforward manner. However, around steep orography, horizontal pressure gradients will be more difficult to calculate because the pressure is not known along constant horizontal surfaces but along terrain-following surfaces. Consequently, much work has gone into accurate evaluations of horizontal pressure gradients using pressure data from different layers (e.g., Zängl 2012). The increased accuracy will reduce the curl of the pressure gradient but is not guaranteed to remove it. It is also possible to eliminate pressure gradient errors in the absence of stratification (Botta et al. 2004).

To eliminate errors associated with sloping coordinate surfaces, cut cells can be used adjacent to the orography (Adcroft et al. 1997; Bonaventura 2000; Steppeler et al. 2002; Good et al. 2013) so that horizontal grid layers intersect with the orography. However, it is difficult to maintain the resolution of the boundary layer at mountain peaks with cut cells and nonorthogonal distortions will still exist between-cut and non-cut cells next to the ground, meaning that pressure gradients will still not be curl free.

The common approach of using vertical and horizontal velocity components as prognostic variables with terrain-following coordinates implies that the vertical velocity is a covariant component of the velocity whereas the horizontal velocity is a contravariant component. [An exception is examined by Simarro and Hortal (2012), who use contravariant velocity components in all directions.] On horizontal, nonorthogonal grids, regardless of using Arakawa B or C grids (Rančić et al. 1996; Adcroft et al. 2004; Thuburn et al. 2014; Weller 2014), the CD grid (a blend between the Arakawa C and D grids; Putman 2007; Harris and Lin 2013), or discontinuous Galerkin (Nair et al. 2005), the covariant velocity is used as the prognostic variable. This means that, on nonorthogonal horizontal grids, pressure gradients can be curl free (Thuburn and Cotter 2012). In this paper, we will explore the use of covariant velocity components as prognostic variables in all directions in combination with terrain-following grids in Cartesian space. This will enable us to calculate pressure gradients that are curl free and consequently not a spurious vorticity source. This follows recent mimetic discretizations on nonorthogonal horizontal grids (e.g., Thuburn and Cotter 2012; Thuburn et al. 2014; Weller 2014). This work entails applying the horizontal discretization described by Weller (2014) in a vertical slice rather than in the horizontal plane in order to achieve some of the same mimetic properties.

Implicit treatment of gravity waves is necessary for using a long time step for strongly stratified flow. If gravity waves are treated explicitly, there will be a time-step restriction based on the stratification. The semi-implicit method including the implicit treatment of gravity waves, as described by Cullen (1990) and Tanguay et al. (1990), involves separating the thermodynamic variables into hydrostatically balanced and perturbation variables. The use of hydrostatically balanced reference profiles that are uniform in time and in the horizontal directions leads to the cancellation of various terms, which consequently simplifies the algorithm. But the perturbation parts can be large and, as a consequence, if linearization assumptions are made, these will not always be accurate.

To avoid large deviations from reference profiles, Davies et al. (2005) and Melvin et al. (2010) use a reference profile consisting of the profile from the previous time step and so the profile about which the model is linearized is no longer in hydrostatic balance. This means that fewer approximations are made but the semi-implicit technique is more complicated since all terms are retained. The retention of all mean and perturbation terms and the description involving the semi-Lagrangian method makes the presentation of the technique complicated. The description of the semi-implicit, semi-Lagrangian (SISL) algorithm employed by Qian et al. (1998) is also very complicated and we conjecture that the semi-implicit solution of the fully compressible equations has not been taken up so widely because these descriptions are so complex.

Gravity waves have also been treated implicitly in models of various simplified equation sets, such as soundproof or pseudo-incompressible (e.g., Smolarkiewicz et al. 2001; Smolarkiewicz and Szmelter 2011; Durran and Blossey 2012; Weller 2014). The use of simplified equation sets often implies that a global Poisson distribution must be solved rather than a global Helmholtz problem, which does not reduce the computational cost. However, an understanding of simplified equation sets can inform the design of solution algorithms for the fully compressible Euler equations, since the large, stiff terms are the same. To move away from the complication of using mean and perturbation variables, Benacchio et al. (2014) describe a method of treating sound but not gravity waves implicitly in which a blend between fully compressible and pseudo-incompressible dynamics can be made.

This article describes a new discretization of the fully compressible Euler equations suitable for strongly stratified flow over orography. The discretization has exactly curl-free pressure gradients, implying that the pressure gradient term is not a spurious source of vorticity. A new technique for treating gravity waves implicitly is presented that does not rely on a background mean state, a hydrostatic mean state, or perturbation variables, and that works on a Lorenz C grid. This simplified approach enables more clarity in ensuring the conservation of mass. The numerical method is described in section 2, some test cases and results demonstrating the properties of the method are presented in section 3, and conclusions are drawn in section 4.

## 2. Numerical method

The numerical method comprises the following elements:

solution of the nonlinear, fully compressible Euler equations in flux form;

semi-implicit treatment of acoustic and gravity waves and explicit treatment of advection;

no explicitly defined reference profile or hydrostatic profile and no reliance on perturbation variables;

exact conservation of mass;

curl-free pressure gradients over orography, following the technique of Weller (2014);

a split space–time (method of lines) multidimensional cubic upwind advection scheme;

Lorenz staggering of

*θ*and Π (with some Charney–Phillips elements within each time step); andthe C-grid finite-volume method for spatial discretization.

This numerical method has been implemented using OpenFOAM 2.3 (OpenFOAM 2014). The implementation described in this paper is available for download online [http://www.met.rdg.ac.uk/~sws02hs/AtmosFOAM/ExnerFoam.tar.gz].

### a. The fully compressible Euler equations

The fully compressible, nonrotating Euler equations in flux form (and advective form for potential temperature) are

where *ρ* is the density, **u** is the velocity, **g** is the acceleration due to gravity, *c*_{p} is the heat capacity at constant pressure, *θ* = *T*(*p*_{0}/*p*)^{κ} is the potential temperature, *T* is the temperature, *p* is the pressure, *p*_{0} is a reference pressure, Π = (*p*/*p*_{0})^{κ} is the Exner function of pressure, and is the ratio of the gas constant to the heat capacity. Both forms of the potential temperature equation will be used in this discretization.

In this *θ*–Π form, a curl-free discretization of **∇**Π does not automatically lead to a curl-free discretization of **∇***p* and consequently pressure gradients may still be spurious sources of vorticity. In the continuous equations, pressure gradients should only be a source of vorticity if pressure gradients are not parallel to density gradients; that is, the solenoidal term, **∇***p ***× ****∇***ρ*, is not zero. If we discretize *c*_{p}*ρθ***∇**Π so that **∇**Π is curl free, it does not follow that there will be no spurious vorticity source. However, there should, at least, be no spurious vorticity source due to the discretization of **∇**Π.

The thermodynamic variables of *θ* and Π are used in order to treat gravity waves implicitly following Davies et al. (2005). The *θ* and Π in the *c*_{p}*ρθ***∇**Π term are treated implicitly but *ρ* in *ρ***g** and in *c*_{p}*ρθ***∇**Π is treated explicitly. The important point is that the same *ρ* is used for both of these terms, which define the hydrostatic balance.

### b. Spatial discretization

The spatial discretization is a C-grid staggered finite volume with Lorenz staggering of thermodynamic variables using covariant velocity components as prognostic variables at the faces between cells. None of the spatial discretization described assumes a structured grid and the implementation is for an arbitrarily structured 3D grid. However, all of the test cases described in section 3 use 2D, terrain-following, structured grids.

For most interpolations, the arithmetic mean is used. The exception is for advection where an upwind multidimensional cubic fit is used. The arithmetic mean is second-order accurate only on uniform grids. For nonuniform grids, alternatives will be needed in order to maintain second-order accuracy but care will be needed to maintain balance and conservative energy transfers. For example, in some situations, volume-weighted interpolation may be preferred to linear or to higher order.

#### 1) Notation

A variable *ψ* located at a cell center is given a subscript *c*: *ψ*_{c}, where *c* is the cell number. A variable, *ψ* located on a face is given subscript *f*: *ψ*_{f}, where *f* is a face number. A variable without a subscript implies an array of all of the cell or face values over the entire grid. Interpolation of cell center values to face values is denoted with subscript *F*: *ψ*_{F}. Reconstruction of cell values from face values is denoted with subscript *C*: *ψ*_{C}; *f* ∈ *c* means the faces of cell *c* and *c* ∈ *f* means the (two) cells either side of face *f*.

#### 2) Prognostic variables

The prognostic variables are the cell center Exner function, Π_{c}; the cell center potential temperature, *θ*_{c} (hence Lorenz staggering); and the momentum at the cell faces in the cell center to cell center direction, *V*_{f} = *ρ*_{f}**u**_{f} · **d**_{f}, where the vector **d**_{f} is defined for each face and is the vector between the cell centers on either side of the face. These variables and vectors are shown in Fig. 1.

#### 3) Cell center and normal velocities from prognostic velocities (operator *H*)

To solve the continuity equation using Gauss’s divergence theorem, we will need the mass flux over every cell face as a diagnostic variable. This is denoted *U*_{f} = *ρ*_{f}**u**_{f} · **S**_{f}, where the face area vector, **S**_{f}, is normal to each face with the magnitude of the face area. To find the field of *U* values from the field of *V* values, we need operator *H* [following the notation of Thuburn and Cotter 2012)]:

Thuburn et al. (2014) define a symmetric, positive-definite *H* for two-dimensional grids with centroidal duals. To use an *H* suitable for three-dimensional, arbitrary grids, we use an *H* similar to that defined by Weller (2014), which is asymmetric and so does not guarantee energy conservation:

where (*ρ***u**)_{F} is the momentum vector interpolated from cell centers onto faces using arithmetic mean interpolation, . The second term in Eq. (6) is a correction to ensure that *H* is diagonal wherever the grid is orthogonal. The cell center momentum, (*ρ***u**)_{C}, is reconstructed from surrounding values of *V*_{f′}:

where is a 3 × 3 tensor and so the inversion of the tensor sum is a local operation that can be calculated once for each cell of the grid rather than at each time step. Equation (7) is a least squares fit that reconstructs uniform vector fields exactly and so it is first-order accurate on arbitrary grids (and second order on uniform grids). To prove the consistency of Eq. (7), we can assume that *ρ***u** = (*ρ***u**)_{f} = (*ρ***u**)_{c} is uniform and see if Eq. (7) reconstructs this uniform velocity field exactly. So we move the inverted tensor to the lhs to give . Each term in the sum on the lhs is equal to , which is identical to the terms in the sum on the rhs only if , which is in fact the definition of *V*_{f}.

The use of *V* (covariant momentum component) rather than *U* (contravariant component) as a prognostic variable was recommended by Thuburn and Cotter (2012) for nonorthogonal horizontal grids in order to achieve a combination of mimetic properties, including curl-free pressure gradients. Although the asymmetric *H* has not been proved to conserve energy, Weller (2014) showed that it gives the same unity amplification factors for the solution of the linearized shallow-water equations as the symmetric *H*.

#### 4) Gradients

For a cell-centered, scalar field, Ψ_{c}, two different types of gradients are defined. For the C-grid-staggered method with *V* = *ρ***u** · **d** as the prognostic variable, the gradient at the face in direction **d** is required:

where *n*_{f} = 1 if **S**_{f} points outward from the cell and −1 otherwise. This simple two-point gradient leads to curl-free pressure gradients. For the solution of the advective form potential temperature equation, the gradient at the cell center is also needed and is defined using Gauss’s theorem:

where the cell has volume *V*_{c}. The interpolation of Ψ from cell centers onto faces to calculate Ψ_{F} in Eq. (9) uses an arithmetic mean interpolation. For solving the potential temperature equation in advective form, the potential temperature gradient at the face in the plane normal to **d** is needed. This is interpolated from the cell-centered potential temperature gradient using arithmetic mean interpolation, , and the component parallel to **d** is not used. In general, **d**_{f} · (**∇**_{c}Ψ)_{F} ≠ |**d**_{f}|**∇**_{d}Ψ but changes to equality for linearly varying fields, for which this discretization would be perfect.

#### 5) Divergence

Divergences are calculated at cell centers using Gauss’s divergence theorem, for example, for scalar field Ψ and vector field **v**, both defined at cell centers,

or since momentum component *U* = *ρ***u** · **S** is defined at the face, then , which is simply denoted **∇** · *U*. Similarly, , which is denoted **∇** · *U*Ψ. A multidimensional cubic least squares fit over an upwind biased stencil of cells is used to calculate Ψ_{F}, which is described in section 8. For solving the momentum equation on the face, the nonlinear advection term is needed on the face. This is interpolated from the cell-centered values using arithmetic mean interpolation: .

#### 6) Perpendicular component of velocity

For the implicit treatment of gravity waves using the advective form of the potential temperature equation, the component of the velocity perpendicular to **d**_{f} will be needed:

where **u**_{F} is calculated using arithmetic mean interpolation from **u**_{C}, which is reconstructed from *V*/*ρ*_{F} as in Eq. (7).

#### 7) Interpolations for Lorenz staggering

Using Lorenz staggering, *θ*, Π, and *ρ* are all stored at cell centers and, where needed, interpolated onto faces using the arithmetic mean: and . However, as will be described in section 2 below, in the course of one time step, *θ* is also advanced on the face using the advective form of the potential temperature in Eq. (4). This is denoted *θ*_{f}. At the beginning of the time step, *θ*_{f} is set to *θ*_{F} by interpolating from *θ*_{c} but then, during a time step, *θ*_{f} is advanced independently from *θ*_{c}. Charney–Phillips staggering could be achieved by setting *θ*_{c} from *θ*_{f} at the beginning of the time step instead but this has not yet been done and care would be needed to maintain the same level of energy conservation.

#### 8) Advection of momentum and potential temperature

The interpolation operations, **u**_{F} and *θ*_{F}, in the terms and control the advection of momentum and potential temperature and so should be undertaken using an upwind-biased interpolation scheme. We have used a least squares fit to a multidimensional cubic using an upwind-biased stencil of cells (Weller et al. 2009). In two dimensions, the multidimensional cubic is

omitting terms in *y*^{3}, where *x* is the direction normal to a cell face and *y* is perpendicular to *x*. Coefficients *a* to *i* are set from a least squares fit to the cell data in the stencil. The least squares problem involves a 9 × *m* matrix singular-value decomposition for every face where *m* is the size of the stencil. However, this is purely a geometric calculation and is therefore a preprocessing activity since the grid is fixed. This generates a set of weights for calculating Ψ_{F} from the cell values in the stencil, leaving *m* multiplies for each face for each call of the advection operator. The stencils are found for three-dimensional, arbitrarily structured grids by finding the face(s) closest to upwind of the face we are interpolating onto, taking the two cells on either side of the upwind face(s) and then taking the vertex neighbors of those central cells. For a two-dimensional structured grid, this gives the stencil shown in Fig. 2a.

The advection scheme is not an important part of the algorithm described. Other good advection schemes, monotonic and/or forward in time, could be used instead.

#### 9) Sponge layer

Following Melvin et al. (2010), a damping term is added to the momentum equation to suppress wave reflections at the rigid lid. This term is −*μρ***u**, where *μ* is nonzero only for vertical velocities near the model top. The distribution of the sponge layer is

where is the distance from position **x** to the surface, *z*_{B} is the height of the bottom of the sponge layer, and *z*_{t} is the height of the model top.

### c. Semi-implicit solution technique

Terms involving acoustic and gravity waves are solved using Crank–Nicholson (trapezoidal) time stepping with no off centering. Advection is treated explicitly with no splitting between explicit and implicit terms (details below). Two outer iterations are performed for each time step so that terms treated explicitly are updated for the second iteration. We will first describe the advection of *ρ* and *θ*, then the derivation of the discretized Helmholtz equation, and finally give a summary of the whole solution procedure.

#### 1) Advection of *ρ* and *θ*

The first task, at the beginning of each outer iteration of each time step, is to solve the continuity and potential temperature equations explicitly, with identical fluxes *U* for each so that they are transported consistently:

where Δ*t* is the time step, superscript *n* represents values from the previous time level, *n* + 1 gives values at the new time level, and ℓ represents lagged values. At the beginning of the time step, values at level ℓ are set to values at level *n* and then these lagged values are updated as soon as new values are available. So at convergence, values at ℓ and *n* + 1 are the same (if enough outer iterations are taken). Crank–Nicholson time stepping uses .

Since the advection is treated explicitly, the time step is limited according to the multidimensional definition of Courant number for cell *c*:

so that Co < 1 is required.

Next in the outer iteration, the Helmholtz equation is solved for Π^{n+1}.

#### 2) Derivation of the discretized Helmholtz equation

A simultaneous solution in all of the prognostic variables together is needed in order to treat acoustic and gravity waves implicitly. To construct a Helmholtz equation in just one variable (Π^{n+1}), the momentum, continuity, and potential temperature equations are combined by hand. First, the potential temperature equation is substituted into the momentum equation to replace *θ* in the *c*_{p}*ρθ***∇**Π term with *V*^{n+1} and then the momentum equation is substituted into the continuity equation to replace *V*^{n+1} with Π^{n+1}. Finally, *ρ*^{n+1} must be replaced by Π^{n+1} on the left-hand side of the continuity equation using a linearization of the equation of state in order to create a Helmholtz equation for Π^{n+1}.

First, we take the dot product of the momentum Eq. (1) with **d** and discretize in time to get an equation for *V*^{n+1}:

where (∂*V*/∂*t*)^{n} is the term in curly brackets in Eq. (14) but from time-level *n*. We have not yet said how we define *ρ*_{f} and *θ*_{f}. This will be done below.

Following the semi-implicit solution technique of Davies et al. (2005), in Eq. (14) is calculated from the advective form of the potential temperature equation,

so that in Eq. (14) can be replace by *V*^{n+1} (all other terms being lagged or from the previous time level). Note that is used rather than . That is, *θ* on the face from the previous time step is interpolated from the prognostic, cell center *θ*_{c} rather than storing *θ*_{f} from one time step to the next, which would result in an overspecification of *θ*. Equation (15) can be rewritten

so that the part

is calculated explicitly. From Eq. (16), can now be substituted into the discretized momentum equation [Eq. (14)]. This can be rearranged so that terms involving *V*^{n+1} are on the lhs. Additionally, one instance of **∇**_{d}Π^{n+1} is replaced by **∇**_{d}Π^{ℓ} so that the equation is linear in implicit terms:

where *G* takes a form similar to that defined by Davies et al. (2005):

and

Fixed flow-rate boundary conditions are imposed on *V*′ and the remaining terms of *V*^{n+1}, the gravity and pressure gradient, are set to cancel exactly on boundaries. Setting *V*′ to zero and at rigid boundaries gives no flow across the boundaries as long as the boundary faces are orthogonal (**d × S** = 0 on the boundary). This can always be enforced by setting **d** to be parallel to **S** on the boundaries.

It may be counterintuitive that *ρ*_{f} is treated explicitly in Eq. (18) since we are treating gravity waves implicitly. However, this follows from Davies et al. (2005), who solve the advective form of the momentum equation. The important detail is to use the same density in the gravity and pressure gradient terms of Eq. (18).

To make this into a Helmholtz equation for Π^{n+1}, we need to replace *ρ*^{n+1} on the lhs with a linear function of Π^{n+1}. This can be done using the equation of state [Eq. (5)]:

where

Because of the low powers of *ρ* and *θ* in Ψ (assuming that *κ* = 0.288), Ψ varies less than *ρ* and *θ*, so the above linearization is useful and leads to convergent outer iterations (in the tests so far undertaken but analysis is needed). Substituting Eq. (22) into Eq. (21) gives the Helmholtz equation for Π^{n+1}:

Given the spatial discretization defined, Eq. (23) is a sparse matrix equation that could be solved to find Π^{n+1}. However, to simplify the construction of the matrix, the operator *H* is split into its diagonal and off-diagonal components and only the diagonal components are treated implicitly: *H* = *H*_{d} + *H*_{off}. So the final term in Eq. (23) becomes

This version is not considered better, just simpler to implement. This version would not be stable for long time steps for highly nonorthogonal grids since too much of the pressure gradient would be treated explicitly. However, it is stable for the test cases described in this paper.

This leads to a sparse matrix that is solved using the conjugate gradient solver from OpenFOAM (2014) with incomplete Cholesky preconditioning.

We now come to how is defined in Eqs. (14)–(23). The algorithm, as defined so far, has too many prognostic variables: *ρ*, Π, *θ*, and *V*, and the continuity equation is used to advance both *ρ* and Π independently. The overspecification is removed by setting using *ρ* advanced from the continuity equation and then setting

This ensures that, over the course of a long simulation, *ρ* advanced from the continuity equation and ΨΠ do not drift.

#### 3) Summary of the semi-implicit solution procedure

The entire update procedure for one time step is given in Algorithm 1. Note that, while the mathematical description talks about values at time levels *n*, *n* + 1, and ℓ, only values at levels *n* and *n* + 1 need storage. In addition, primed variables, *θ*′ and *V*′, use the same storage as *θ*^{n+1} and *V*^{n+1}.

Once Eq. (23) is solved for Π^{n+1}, *V*^{n+1} is updated from Eq. (18) (the back substitution). Unlike in ENDGAME (Davies et al. 2005; Melvin et al. 2010), there is no back substitution for *θ*_{f}. Instead, final solutions for Eqs. (12) and (13) are calculated for the beginning of the next time step.

Regardless of the level of convergence, this solution algorithm will always give the exact local mass conservation since *ρ*_{c} is advanced using fluxes over cell faces from the continuity Eq. (2). However, only at convergence will the density calculated from the continuity equation equal the density calculated from the equation of state (ΨΠ).

### d. Alternative model formulations

To demonstrate the value of the novel aspects of the discretization presented, two alternative approaches are presented and have been implemented for comparisons.

#### 1) Horizontal pressure gradient (∂*p*/∂*x*)

Most models of the global atmosphere use horizontal winds as prognostic variables and require the reconstruction of horizontal pressure gradients (e.g., Klemp 2011; Zängl 2012). A similar approach is presented in order to compare with the new version that uses the *H* operator. This version is called ∂*p*/∂*x*. In this form, *U* is the prognostic variable rather than *V* and *V* is in fact not defined. The momentum equation is formulated in direction **S** rather than **d** (i.e., in the horizontal direction and in the near-vertical direction, normal to the cell faces). The derivation in direction **S** rather than **d** is very similar apart from the gradient at the face in direction **S**, **∇**_{s}Ψ = **∇**Ψ · **S**, which is given by a least squares fit to the linear polynomial:

where the coefficients *a*, *b*, and *c* are set using values of Ψ in stencils around each edge, as shown in Fig. 2. This approach is not as sophisticated a form as that used by Zängl (2012) but it is of a similar complexity to the *H* version so as to make a like-for-like comparison.

#### 2) Explicit solution of gravity waves

To treat gravity waves explicitly and acoustic waves implicitly, and to make no other changes to the formulation, *θ*_{f} = *θ*_{F} is used instead of Eq. (16) and

is used instead of Eq. (19).

## 3. Results

A number of test cases from Skamarock and Doyle (2013) are used to demonstrate the value of the curl-free pressure gradient and the implicit treatment of gravity waves. All of the test cases use a reference pressure of *p*_{0} = 10^{5} Pa in the definition of the Exner pressure, zero viscosity, zero hyperviscosity, and *κ* = 0.287 698. The test cases use different reference temperatures for the definition of the potential temperature.

All of the test cases require hydrostatically balanced initial conditions. To find Π in discrete hydrostatic balance with an initial *θ* field, we numerically solve the Poisson equation,

subject to the boundary conditions *θ***∇**Π = **g** at the ground and lateral boundaries and a fixed Π at the flat upper boundary. The upper boundary value of Π is iteratively adjusted to get Π = 1 at *z* = 0. This is found to be more stable and reliable than setting Π = 1 at *z* = 0. For test cases with prescribed *θ*, one solution of Eq. (27) is needed per value of Π at the upper boundary. One of the test cases specifies uniform *T*, so outer iterations are needed, setting *θ* = *T*/Π between each solution of the Poisson equation.

For some of the test cases, results on different grids are compared. For example, solutions are compared with high-resolution reference solutions. This is done by mapping the reference solution onto the target grid assuming that the reference solution is piecewise constant on each cell and using volume weighting. This requires calculating overlapping volumes between the two grids. This is done using the OpenFOAM mapFields utility.

### a. Resting atmosphere

#### 1) Test case setup

We start with the simulation of a resting stratified atmosphere over a range of hills (Schär et al. 2002) using the test case setup of Klemp (2011). The lower boundary takes the form used by Schär et al. (2002):

where *a* = 5 km, *λ* = 4 km, and *h*_{m} = 1000 m. By specifying the Brunt–Väisälä frequency, we can set *θ*:

and *θ*(*z* = 0) = 288 K. Through initialization, Π is in discrete hydrostatic balance. The resolution is set to Δ*x* = 500 m and Δ*z* = 500 m away from the terrain. The depths of the terrain-following layer are set in two ways to compare with the results of Klemp (2011). First, the *z* coordinates of the gridpoint locations are set using the basic terrain-following (BTF) coordinate definition:

where *z*_{t} is the domain top and *ζ* is the layer height before orography is added. Next, the layer depths are set to follow the SLEVE vertical coordinate (Schär et al. 2002) with decay functions specified following Leuenberger et al. (2010) with *s*_{1} = 4 km and *s*_{2} = 1 km and *n* = 1.35. The BTF and SLEVE grids are shown in Fig. 3. The domain is 20 km in the *x* direction and 20 km in the *z* direction with rigid boundaries at the top, bottom, and sides. This domain is smaller than that used by Klemp (2011) in the *x* direction to reduce run times and the boundaries are rigid rather than open for simplicity. No sponge layer or diffusion is used. All test case results shown have implicit treatments of gravity waves, but for the time step used, the results are almost indistinguishable when using the explicit gravity wave formulation.

#### 2) Test case results

Potential temperature contours after 5 h from the model using the new *H* formulation and the model using the ∂*p*/∂*x* formulation on the BTF and SLEVE grids are presented in Fig. 3, all with implicit treatment of gravity waves. The advection scheme is not very diffusive and no explicit diffusion is used, so the simulation on the BTF grid using the ∂*p*/∂*x* formulation has distorted *θ* contours. The distortions can be reduced by either using the *H* formulation or by using the SLEVE grid and using both makes the contours very flat. This demonstrates the improved hydrostatic balance from using the *H* model formulation.

The maximum (spurious) vertical velocities generated for each of the model runs are shown in Fig. 4, where they are compared with the maximum vertical velocity from Fig. 4 of Klemp (2011) (note different *y* scales). This shows that the ∂*p*/∂*x* formulation on the BTF grid generates the largest spurious vertical velocities, and in the first hour, the erroneous velocities are higher than those of Klemp (2011). This could be due to different initializations or to the higher-order treatment of boundaries by Klemp (2011). However, on the BTF grid, the Klemp (2011) errors grow, unlike the ∂*p*/∂*x* and *H* errors on both grids. Use of either the SLEVE grid or the *H* version reduces the errors to a level similar to the results of Klemp (2011) on the SLEVE grid. The *H* model results are less sensitive to the choice of grid than are the ∂*p*/∂*x* results or the results of Klemp (2011), again demonstrating the value of the *H* model formulation.

The discretization described in this paper does not give exact energy conservation; therefore, it is worth examining the energy conservation and the influence of using the *H* operator on energy conservation. The normalized energy changes from the initial conditions (normalized by the initial total energy) for the simulations using the BTF grid and the SLEVE grid and for the simulations using the ∂*p*/∂*x* formulation and the *H* operator are shown in Fig. 5. The energy conservation using the *H* formulation is at least an order of magnitude better than that using the ∂*p*/∂*x* version on the BTF grid. The dashed lines on the left of Fig. 5 show negative changes, which implies that the *H* formulation mostly loses energy, which is desirable for stability. On the BTF grid, the contributing terms to the energy conservation are shown for both model versions in Fig. 5. Although the *H* version does not conserve energy to machine precision, the transfers between internal and potential energy on short time scales are represented whereas they both increase and decrease in tandem for the ∂*p*/∂*x* version, leading to large energy changes.

There are a few reasons why energy is not conserved precisely in any of the models presented. We are solving the flux form rather than a vector-invariant momentum equation and so the transfer between potential and kinetic energy is not conservative; the advection scheme is upwind biased with an amplification factor less than 1 and so destroys kinetic energy and there are inconsistencies between the *θ* that is advected by a conservative advection scheme and the *θ* that appears in the pressure gradient term, *c*_{p}*ρθ***∇**Π, of the momentum equation.

### b. Schär et al.’s (2002) mountain waves test case

#### 1) Test case setup

The test case described by Schär et al. (2002) simulates flow over mountains with small and large features that are lower and less steep in comparison to those described in section 3a. The lower boundary has the same form [Eq. (28)] and again uses *a* = 5 km and *λ* = 4 km, but for this test case *h*_{m} = 250 m. The initial conditions are set using *N* = 0.01 s^{−1} and a mean wind of *U* = 10 m s^{−1}. We follow Schär et al. (2002) and Klemp et al. (2003) and use a time step of 8 s. Following Melvin et al. (2010), we use a domain of 100 km × 30 km, Δ*x* = 0.5 km, Δ*z* = 300 m, surface temperature of 288 K, and an absorbing layer in the top 10 km of the domain with . The terrain is followed using both the BTF grid and the SLEVE grid. The top and bottom boundaries are rigid with zero flow. The inlet boundary has the prescribed *θ* profile and winds of 10 m s^{−1}, and the outlet boundary is zero gradient. The boundary condition for Π is hydrostatic all around.

This is not a good case for testing the implementation of the implicit gravity waves since *N*Δ*t* = 0.08, which is stable even if gravity waves are treated explicitly (Cullen 1990). The time step could be increased to around 40 s while treating advection explicitly, but this would still not raise *N*Δ*t* above 1.

#### 2) Test case results

The vertical velocities generated by the mountain are shown in Fig. 6 using both model versions (∂*p*/∂*x* and *H*) with implicit gravity waves on the BTF and SLEVE grids. These are compared with the simulations using the *H* model version at 4 times the resolution, a quarter the time step on the SLEVE grid, and with result from Melvin et al. (2010). The four simulations using the ∂*p*/∂*x* and *H* models on both grids are similar. The results on both grids are similar because the advection scheme used accounts properly for the distortions in the grid. The *H* and ∂*p*/∂*x* model versions give similar results for this test case because the small-scale gravity waves generated by the orography are evanescent and so their structure, whether realistic or not, is not present at a few kilometers above the ground.

Differences with the high-resolution solution are not presented because the differences are dominated by boundary reflections and the varying strength of the sponge layer with the time step. This case demonstrates that the advection scheme accounts properly for the grid distortions but is not useful for testing the curl-free pressure gradients or for the implicit treatment of gravity waves.

### c. Linear hydrostatic flow over a hill

#### 1) Test case setup

To demonstrate the value of having implicit gravity waves, it is necessary to go to coarser horizontal resolution to allow for running with a longer time step and hence achieving larger *N*Δ*t*. Decreasing the horizontal resolution brings the resolved flow closer to hydrostatic. The simulations of near-hydrostatic flow are done with the same nonhydrostatic model. The test case of hydrostatic mountain waves from Skamarock and Doyle (2013) and used by Melvin et al. (2010) has a witch of Agnesi hill profile:

with *h*_{m} = 1 m (to ensure that the solution is close to linear), *a* = 10 km, a mean wind speed of 20 m s^{−1}, and an initial isothermal temperature of *T* = 250 K in discrete hydrostatic balance. Following Melvin et al. (2010), an absorbing layer is applied in the top 20 km with . The domain is 240 km wide, centered on the hill, and 50 km deep with grid spacing Δ*x* = 2 km and Δ*z* = 250 m. This resolution is used with time steps of 20 and 50 s, giving Courant numbers of 0.2 and 0.5 and *N*Δ*t* of 0.4 and 1. Coarser resolutions are also used with larger time steps. The boundary conditions for all simulations are as described in section 1.

#### 2) Test case results

Vertical velocity contours for the near-hydrostatic flow over a hill are shown in Fig. 7 for model formulations using explicit and implicit gravity waves and the different time steps and spatial resolutions. The *H* and ∂*p*/∂*x* versions of the model are almost identical for this test case because the grids are practically orthogonal, with a hill height of only 1 m. The well-resolved numerical solutions are similar to the linear analytic, anelastic, nonhydrostatic solution (from Melvin et al. 2010). The version with explicit gravity waves is stable for a time step of 20 s (corresponding to Co = 0.2, maximum *N*Δ*t* = 0.4) but, unlike the version with implicit gravity waves, is not quite stable for a time step of 50 s (corresponding to Co = 0.5, maximum *N*Δ*t* = 1). The version with implicit gravity waves can be run stably at much longer time steps at coarser resolution so that the advection Courant number remains below 1 but with *N*Δ*t* = 2, *N*Δ*t* = 4, and *N*Δ*t* = 10. (Larger values have not been tried.) For the coarser resolutions the accuracy is reduced but, as long as gravity waves are treated implicitly, the simulations remain stable.

Results of this test case demonstrate that gravity waves are treated implicitly, as required, and the model is stable in the presence of strong stratification.

### d. Rising bubble over orography

To test the model in nonlinear flow regimes and to further test the representation of orography, we use the rising bubble test case of Bryan and Fritsch (2002), modified by Good et al. (2013) so that the bubble is rising over orography. This tests the representation of the nonhydrostatic buoyancy and pressure gradient terms on distorted grids such as those associated with terrain-following layers. The nonlinear terms are more important in this case than those with orographically produced gravity waves since there is no mean wind.

#### 1) Test case setup

The rising bubble test case of Bryan and Fritsch (2002) consists of an initially stationary atmosphere with no stratification (*θ* = 300 K) with pressure in hydrostatic balance with this temperature profile. A warm bubble is then placed, centered at (*x*_{c}, *y*_{c}) with temperature perturbation

where and *A*_{x} = *A*_{z} = 2 km. Good et al. (2013) place the bubble higher than do Bryan and Fritsch (2002) to allow for the orography under the initial bubble. The bubble is placed at (*x*_{c}, *y*_{c}) = (0, 4.5 km), directly above the orography. The domain is 20 km wide and 20 km high with Δ*x* = Δ*z* = 100 m. For this resolution, a time step of 2 s is used, which gives a final (maximum) Courant number of 0.47 and *N*Δ*t* of 0.03, and the model is run for 1000 s. Free-slip boundaries are placed all around [unlike Good et al. (2013), but this is not expected to be critical]. The model is run with both no orography and also with a witch of Agnesi hill profile:

with *h*_{m} = 1000 m and *a* = 1000 m. A BTF grid is used over the orography in order to accentuate the differences between the models with and without orography and to accentuate the differences between the *H* and ∂*p*/∂*x* model versions. The hill is sufficiently far beneath the flow generated by the rising bubble that it should not significantly affect the bubble (Good et al. 2013). The no-orography case is compared with a high-resolution simulation that uses Δ*x* = Δ*z* = 31.25 m and Δ*t* = 0.5 s.

#### 2) Test case results

The potential temperature and contours of vertical velocity for the rising bubble over flat terrain and over orography using both model versions (*H* and ∂*p*/∂*x*) are shown in the top row of Fig. 8. The potential temperature and vertical velocity are similar to those shown in Bryan and Fritsch (2002). In particular, the levels of unboundedness (values less that 300 K and greater than 302 K) are similar to those of Bryan and Fritsch (2002). However, the central vertical jet is not as sharp as that of Bryan and Fritsch (2002) and the bubble is developing a nipple (it is beginning to burst), unlike that of Bryan and Fritsch (2002). The differences are not surprising since Bryan and Fritsch (2002) use fifth-order spatial derivatives for the advection terms whereas our advection scheme is second order with cubic interpolations. The differences between the results using Δ*x* = Δ*z* = 100 m (without orography) and the results using Δ*x* = Δ*z* = 31.25 m (without orography) also shown Fig. 8, bottom left). The differences are between −0.7 and 0.5 K. The inclusion of orography below the bubble makes very little difference when using the *H* model version, but differences from the no-orography case ranging from −0.3 to 0.2 K yield larger differences when using the ∂*p*/∂*x* model (−1 to 0.3 K). The extrema are not much altered by the orography but the maximum *θ* are now on either side of the center for the ∂*p*/∂*x* model. The differences between cases with and without orography are larger than those presented by Good et al. (2013) when they used cut cells but smaller than their differences when they used a terrain-following grid (errors up to 1.67 K). Our terrain-following model results are likely better than theirs because we are using an improved advection scheme and curl-free pressure gradients.

For the rising bubble test case, we have also inspected convergence with space and time steps. Normalized ℓ_{2}(*θ*) and ℓ_{∞}(*θ*) errors are calculated for a range of spatial and temporal resolutions in comparison to reference solutions. The error norms are defined as

where *θ*_{T} is the reference solution. When looking at convergence with spatial resolution, the reference solution uses Δ*x* = Δ*z* = 31.25 m and Δ*t* = 0.5 s. When looking at convergence with the time step, all solutions use Δ*x* = 100 m and the reference solution uses Δ*t* = 0.1 s. Convergence with spatial resolution is shown in the top left of Fig. 9 at 400 s after initialization for simulations using Δ*x* = 250 m, Δ*t* = 4 s; Δ*x* = 125 m, Δ*t* = 2 s; and Δ*x* = 62.5 m, Δ*t* = 1 s. Convergence is second order at 400 s after initialization but drops to first order at 1000 s (not shown). The drop to first order is likely to be due to the insufficient resolution of the very sharp gradients, meaning that the theoretical convergence is not met. The convergence with time step (Fig. 9, top right) also mostly shows second-order convergence apart from at the longest time step, where insufficient temporal resolution reduces the accuracy more sharply. This reduced accuracy at the longest time step is the reason why the simulations presented above did not use the longest stable time step.

For the rising bubble test case, Norman et al. (2011) also show the maximum *θ* and vertical velocity for each time step as a function of resolution. Similar plots to theirs are shown in the bottom of Fig. 9, using the same spatial resolutions but with much longer time steps because Norman et al. (2011) use entirely explicit time stepping. There are similarities between our results: for the finer resolutions, the maximum *θ* increases toward the end of the simulation and, after about 800 s, there is a dramatic acceleration in the maximum vertical velocity as the bubble starts to burst.

Results for this test case demonstrate the second-order accuracy of the model and the benefits of the *H* model formulation.

## 4. Discussion and conclusions

A new semi-implicit model of the fully compressible Euler equations has been presented that offers an implicit treatment of gravity waves and the use of covariant components of velocity over orography that permits the calculation of curl-free pressure gradients. This is achieved by solving all of the flux form equations in a finite-volume model without mean and perturbation variables. These properties have enabled the following test case results:

Simulation of a resting, stratified atmosphere over steep terrain with covariant rather than contravariant prognostic velocities has led to smaller spurious velocities, better energy conservation, and a realistic transfer between potential and internal energy.

Simulations of nonhydrostatic gravity waves over orography are not dependent on the type of terrain-following grid.

Simulations with strong stratification and long time steps using a formulation applicable to arbitrary grids, which are not necessarily aligned in the vertical.

An insensitivity to grid distortions when simulating a rising warm bubble is seen

## Acknowledgments

HW acknowledges support from NERC Grant NE/H015698/1. AS acknowledges support from NERC Grant NE/K006797/1. We also acknowledge useful discussions and ideas from John Thuburn, Nigel Wood, Colin Cotter, and Tom Melvin.

## REFERENCES

*Quart. J. Roy. Meteor. Soc.,*

**116,**1253–1258, doi:.

*Quart. J. Roy. Meteor. Soc.,*

**131,**1759–1782, doi:.

*Mon. Wea. Rev.,*

**140,**1307–1325, doi:.

**15,**44–49, doi:.

*Mon. Wea. Rev.,*

**132,**1860–1867, doi:.

*Bull. Amer. Meteor. Soc.,*

**83,**1019–1037, doi:.

*Quart. J. Roy. Meteor. Soc.,*

**132B,**1795–1813, doi:.

*Acta Geophys.,*

**59,**1109–1134, doi:.

**7,**909–929, doi:.

*Quart. J. Roy. Meteor. Soc.,*

**129B,**1989–2010, doi:.

*Geosci. Model Dev.,*

**7,**779–797, doi:.