An energy consistent discretization of the nonhydrostatic equations in primitive variables

TL;DR

This paper presents an energy-conserving discretization of nonhydrostatic equations in spherical geometry, ensuring accurate energy balance and convergence in atmospheric modeling.

Contribution

It introduces a novel energy consistent discretization in primitive variables with a discrete derivative product rule, applicable to terrain-following coordinates and various discretizations.

Findings

Conservation of discrete energy in exact time integration.

Demonstrated convergence with spectral element code.

Applicable to terrain-following and Eulerian/vertically Lagrangian discretizations.

Abstract

We derive a formulation of the nonhydrostatic equations in spherical geometry with a Lorenz staggered vertical discretization. The combination conserves a discrete energy in exact time integration when coupled with a mimetic horizontal discretization. The formulation is a version of Dubos and Tort (2014) rewritten in terms of primitive variables. It is valid for terrain following mass or height coordinates and for both Eulerian or vertically Lagrangian discretizations. The discretization relies on an extension to Simmons and Burridge (1981) vertical differencing which we show obeys a discrete derivative product rule. This product rule allows us to simplify the treatment of the vertical transport terms. Energy conservation is obtained via a term-by-term balance in the kinetic, internal and potential energy budgets, ensuring an energy-consistent discretization with no spurious sources of energy. We demonstrate convergence with respect to time truncation error in a spectral element code with a HEVI IMEX timestepping algorithm

Equations

$$\frac{\partial \pi }{\partial s}\frac{\partial s}{\partial z}=\frac{\partial \pi }{\partial z}=-\rho g$$

$$\rho=-\frac{\partial\pi}{\partial s}\Big{/}\frac{\partial\phi}{\partial s}.$$

$${\int }_{z_{\text{surf}}}^{z_{\text{top}}}\rho Xdz=\frac{1}{g}{\int }_{s_{\text{top}}}^{s_{\text{surf}}}\frac{\partial \pi }{\partial s}Xds.$$

$${\nabla }_{s}p,{\nabla }_s\times w\hat{k},{\nabla }_s\cdot \mathbf{u},{\nabla }_s\times \mathbf{u}$$

$$\frac{1}{\rho }\nabla p\cdot \hat{k}$$

$$\frac{1}{\rho }\{\nabla p{\}}_{\text{h}}$$

$$\mu=\frac{\partial p}{\partial s}\Big{/}\frac{\partial\pi}{\partial s}.$$

$$\frac{1}{\rho }{\nabla }_{s}p=c_p{\theta }_v{\nabla }_s\Pi .$$

$$\frac{\partial \mathbf{u}}{\partial t}+({\nabla }_s\times \mathbf{u}+f\hat{k})\times \mathbf{u}+\frac{1}{2}{\nabla }_s{\mathbf{u}}^2+\dot{s}\frac{\partial \mathbf{u}}{\partial s}+c_p{\theta }_v{\nabla }_s\Pi +\mu {\nabla }_s\varphi$$

$$\frac{\partial w}{\partial t}+\mathbf{u}\cdot {\nabla }_{s}w+\dot{s}\frac{\partial w}{\partial s}+g-g\mu$$

$$\frac{\partial \varphi }{\partial t}+\mathbf{u}\cdot {\nabla }_s\varphi +\dot{s}\frac{\partial \varphi }{\partial s}-gw$$

$$\frac{\partial \theta }{\partial t}+\mathbf{u}\cdot {\nabla }_s\theta +\dot{s}\frac{\partial \theta }{\partial s}$$

$$\frac{\partial }{\partial t}\left(\frac{\partial \pi }{\partial s}\right)+{\nabla }_s\cdot \left(\frac{\partial \pi }{\partial s}\mathbf{u}\right)+\frac{\partial }{\partial s}\left(\frac{\partial \pi }{\partial s}\dot{s}\right)$$

$$G_T(p_0,\theta ,q_i)$$

$$\frac{1}{\rho }$$

$$\frac{D{\theta }_v}{Dt}=0.$$

$$T=\theta {\left(\frac{p}{p_0}\right)}^{R^{\ast }/c_p^{\ast }}$$

$$\frac{\partial \mathbf{u}}{\partial t}+({\nabla }_s\times \mathbf{u}+f\hat{k})\times \mathbf{u}+\frac{1}{2}{\nabla }_s\left({\mathbf{u}}^2+w^2\right)-w{\nabla }_{s}w+\dot{s}\frac{\partial \mathbf{u}}{\partial s}+c_p{\theta }_v{\nabla }_s\Pi +\mu {\nabla }_s\varphi$$

$$\frac{\partial w}{\partial t}+\mathbf{u}\cdot {\nabla }_{s}w+\dot{s}\frac{\partial w}{\partial s}+g-g\mu$$

$$\frac{\partial \varphi }{\partial t}+\mathbf{u}\cdot {\nabla }_s\varphi +\dot{s}\frac{\partial \varphi }{\partial s}-gw$$

$$\frac{\partial \Theta }{\partial t}+{\nabla }_s\cdot (\mathbf{u}\Theta )+\frac{\partial }{\partial s}\left(\dot{s}\Theta \right)$$

$$\frac{\partial }{\partial t}\left(\frac{\partial \pi }{\partial s}\right)+{\nabla }_s\cdot \left(\frac{\partial \pi }{\partial s}\mathbf{u}\right)+\frac{\partial }{\partial s}\left(\frac{\partial \pi }{\partial s}\dot{s}\right)$$

$$\Theta =\frac{\partial \pi }{\partial s}{\theta }_v\frac{\partial \varphi }{\partial s}=-R\Theta \frac{\Pi }{p}$$

$$\dot{s}\frac{\partial \varphi }{\partial s}=gw-\mathbf{u}\cdot {\nabla }_s\varphi .$$

$$w=\mathbf{u}\cdot \frac{1}{g}{\nabla }_s\varphi$$

$$\pi =A(s)p_0+B(s)p_s$$

$$\dot{s}\frac{\partial \pi }{\partial s}=B(s){\int }_{s_{\text{top}}}^1{\nabla }_s\cdot \frac{\partial \pi }{\partial s}\mathbf{u}ds-{\int }_{s_{\text{top}}}^s{\nabla }_s\cdot \frac{\partial \pi }{\partial s}\mathbf{u}ds$$

$$w=\mathbf{u}\cdot \frac{1}{g}{\nabla }_s\varphi$$

$$K=\frac{1}{2}\frac{\partial \pi }{\partial s}{\mathbf{v}}^{2}I=c_p\Theta \Pi +\frac{\partial \varphi }{\partial s}p+\hat{p}{\varphi }_{\text{top}}P=\frac{\partial \pi }{\partial s}\varphi$$

$$\frac{\partial }{\partial t}P=-\varphi {\nabla }_s\cdot \left(\frac{\partial \pi }{\partial s}\mathbf{u}\right)-\frac{\partial \pi }{\partial s}\mathbf{u}\cdot {\nabla }_s\varphi -\varphi \frac{\partial }{\partial s}\left(\frac{\partial \pi }{\partial s}\dot{s}\right)-\frac{\partial \pi }{\partial s}\dot{s}\frac{\partial \varphi }{\partial s}+\frac{\partial \pi }{\partial s}gw$$

$$\frac{\partial }{\partial t}I=c_p\Pi \frac{\partial \Theta }{\partial t}+c_p\Theta \frac{\partial \Pi }{\partial t}+\frac{\partial }{\partial t}\left(\frac{\partial \varphi }{\partial s}p\right)+\hat{p}\frac{\partial }{\partial t}({\varphi }_{\text{top}}).$$

$$\frac{\partial }{\partial t}I=c_p\Pi \frac{\partial \Theta }{\partial t}+p\frac{\partial }{\partial s}\frac{\partial \varphi }{\partial t}+\hat{p}\frac{\partial }{\partial t}({\varphi }_{\text{top}})$$

$$\frac{\partial }{\partial t}I=c_p\Pi \frac{\partial \Theta }{\partial t}-\frac{\partial p}{\partial s}\frac{\partial \varphi }{\partial t}$$

$$\frac{\partial }{\partial t}I=-c_p\Pi {\nabla }_s\cdot \left(\Theta \mathbf{u}\right)-c_p\Pi \frac{\partial }{\partial s}\left(\Theta \dot{s}\right)+\frac{\partial p}{\partial s}\mathbf{u}\cdot {\nabla }_s\varphi +\frac{\partial p}{\partial s}\left(\dot{s}\frac{\partial \varphi }{\partial s}\right)-\frac{\partial p}{\partial s}gw$$

$$\int c_p\Pi \frac{\partial }{\partial s}\left(\Theta \dot{s}\right)ds=\int \frac{\partial p}{\partial s}\left(\dot{s}\frac{\partial \varphi }{\partial s}\right)ds.$$

Full text

\draftfalse\journalname

Journal of Advances in Modeling Earth Systems (JAMES)

Computational Science, Sandia National Laboratories, Albuquerque, New Mexico, USA. Department of Land, Air and Water Resources, University of California, Davis, California, USA. NVIDIA, Santa Clara, California, USA Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000 Grenoble, France

\correspondingauthor

Mark Taylormataylo@sandia.gov

{keypoints}

We give a discrete Hamiltonian formulation of the nonhydrostatic equations in primitive variables.

The formulation supports mass or height terrain following coordinates.

The Lorenz staggered vertical discretization obeys a derivative product rule.

Plain Language Summary

Energy consistent discretizations have proven useful in guiding the development of numerical methods for simulating fluid dynamics. They ensure that the discrete method does not have any spurious sources of energy which can lead to unstable and unrealistic simulations. Here we provide an energy consistent discretization of the equations used by global models of the Earth’s atmosphere. The discretization is written in terms of standard variables in spherical coordinates and supports a wide variety of terrain following vertical coordinates. It can be used with any horizontal discretization that has a discrete version of the integration-by-parts identity.

1 Introduction

We present a discrete formulation of the nonhydrostatic equations designed for use in a global model of the Earth’s atmosphere. It has been implemented in the version of the High Order Method Modeling Environment (HOMME), [Dennis \BOthers. (\APACyear2005), Dennis \BOthers. (\APACyear2012)] used by the Energy Exascale Earth System Model (E3SM), [Golaz \BOthers. (\APACyear2019)]. We follow a Hamiltonian approach so that the formulation will be energetically conservative and consistent in the sense of \citeAGassmannHerzog2008 when coupled with a mimetic horizontal and vertical discretizations. Energy consistency is obtained via a term-by-term balance in the discrete kinetic, internal and potential energy budgets, ensuring an energy conserving discretization with no spurious energy sources or sinks.

Recent examples of energy consistent global atmosphere models running on unstructured grids include [Taylor (\APACyear2011), Gassmann (\APACyear2013), Dubos \BOthers. (\APACyear2015)]. Here we extend the discretization from \citeATaylor11 to the nonhydrostatic equations using a formulation modeled after \citeADubosTort2014. The main difference is that we formulate the Hamiltonian in terms of the commonly used primitive variables in spherical coordinates. This may simplify implementation in existing modeling frameworks. The discretization supports terrain following mass or height coordinates [Kasahara (\APACyear1974)], which differ only in their boundary conditions and one diagnostic equation. It retains a 2D vector invariant form for horizontal transport and advective form for vertical transport, facilitating both Eulerian and vertically Lagrangian [Lin (\APACyear2004)] coordinates in a single code base. In the vertically Lagrangian (mass or height) coordinate, the vertical coordinate surfaces float in such a way as to ensure no mass flux through the surfaces. In the Eulerian (mass or height) coordinate, the transformed equations are solved in an Eulerian reference frame.

For our vertical discretization, we use \citeASB81 vertical differencing with a \citeALorenz1960 staggering. We introduce an extension to SB81 which obeys a discrete derivative product rule in addition to the well known SB81 discrete integral properties. This product rule ensures the discrete vertical transport terms remain energetically neutral despite being in non-Hamiltonian form. The discretization is energetically consistent when used with the shallow atmosphere approximation, total air mass density as a prognostic variable and general moist equations of state, but with a common approximation in the equation for virtual potential temperature.

The formulation is appropriate for any mimetic spatial discretization in spherical geometry such as \citeAthuburn2009numerical,ringler2010unified,TayFour10,cotter2012mixed,thuburn2015primal, Gassmann2018. Mimetic discretizations are those that mimic key integration and annihilator properties of the divergence, gradient and curl operators [Samarskiĭ \BOthers. (\APACyear1981), Nicolaides (\APACyear1992), Shashkov \BBA Steinberg (\APACyear1995), Bochev \BBA Hyman (\APACyear2006)]. They are commonly used in order to preserve the Hamiltonian structure of the continuum equations [Salmon (\APACyear2004), Salmon (\APACyear2007), Gassmann \BBA Herzog (\APACyear2008)]. Here we use a $p=3$ collocated spectral element discretization due to its retention of 4th order accuracy on unstructured spherical grids and its explicit mass matrix which allows for the use of efficient horizontal explicit vertically implicit (HEVI) timestepping. This efficiency comes at the cost of erratic behavior of grid scale waves [Melvin \BOthers. (\APACyear2012), Ainsworth (\APACyear2014)] which are damped as in \citeAullrich2018. We demonstrate convergence with respect to time truncation error with a 5 stage Runge-Kutta IMEX algorithm.

2 Continuum equations in terrain following coordinates

Following \citeAKasahara1974 and \citeALaprise92 we write the equations of motion in a terrain following coordinate and make use of a pseudo-density based on the hydrostatic pressure. Our physical domain will be Earth’s 3D atmosphere in spherical coordinates with latitude $\varphi$, longitude $\lambda$, and radial coordinate $r$. The altitude from an arbitrary reference level $r_{0}$ is denoted $z=r-r_{0}$. A terrain following coordinate $s$ is defined so that the bottom boundary corresponds to a surface of constant $s$. We denote the velocity $\mathbf{v}=(u,v,w)$ with $\mathbf{u}=(u,v)$ denoting the the horizontal velocity tangent to constant $z$ surfaces. For vertical velocity, we use both $\dot{z}=w=Dz/Dt$ and $\dot{s}=Ds/Dt$, with $D/Dt$ the standard material derivative. We take $\rho$ as the total air density, $\hat{k}$ the radial unit vector, and $\phi=gz$ as the geopotential with constant gravitational acceleration $g$.

We assume $z$ is a monotone function of $s$, so the pseudo-density $\partial\pi/\partial s$ can be defined in terms of the hydrostatic pressure $\pi$ and density,

[TABLE]

with

[TABLE]

Vertical integrals of mass and other density weighted quantities are transformed to terrain following coordiantes via

[TABLE]

We denote the three-dimensional physical coordinate-independent differential operators in the usual notation, $\nabla\cdot\mathbf{v},\nabla p,\nabla\times\mathbf{v}$ and use

[TABLE]

to denote the two-dimensional operators on $s$ surfaces [Ullrich \BOthers. (\APACyear2017)].

2.1 Pressure gradient formulation

We can separate the pressure gradient term into vertical and horizontal components, written in terrain following coordinates as

[TABLE]

where for compactness we introduce

[TABLE]

We use the usual definition of virtual Exner pressure $\Pi=(p/p_{0})^{R/c_{p}}$ and virtual potential temperature $\theta_{v}=p/(\rho R\Pi)$. These are defined in terms of $p$ and $\rho$ from the moist atmosphere but using the ideal gas constant for dry air, $R$, and the specific heat at constant pressure for dry air, $c_{p}$. These definitions do not imply any assumptions regarding the equation of state and are used to write the pressure gradient as

[TABLE]

2.2 Shallow atmosphere equations

We start with a standard form of the nonhydrostatic equations with the shallow atmosphere approximation in terrain following coordinates. Our goal is to derive a energetically consistent discretization for the adiabatic terms, so we initially neglect dissipation, sources and sinks of moisture and other forcing terms. In the full model, dissipation will be added through choice of timestepping and vertical remap algorithms and horizontal hyperviscosity. Following \citeALaprise92, we write the equations as

[TABLE]

with Coriolis term $f\hat{k}$.

2.3 Equation of state

To close the equations we need an equation of state in order to compute $p$, which is needed to derive $\Pi$ and $\theta_{v}$. As written above, we can support a very general equation of state via the Gibbs potential approach recommended in \citeAThuburn2017. Given a Gibbs potential $G(p,T,q_{i})$, with temperature $T$ and $q_{i}$ representing various forms of water, the algebraic system

[TABLE]

then defines $p$ and $T$ implicitly as functions of $\rho,\theta$ and $q_{i}$. This system can be solved efficiently with Newton iteration.

Here we simplify the equations substantially by assuming that

[TABLE]

This is a common assumption for moist air. It neglects small terms, the magnitude of which are discussed in \citeAstaniforth2006. Under this approximation, we can prognose $\theta_{v}$ instead of $\theta$ and then the relation $\theta_{v}=p/(\rho R\Pi)$ is sufficient to close the system. The water species $q_{i}$ then become passive tracers during the adiabatic part of the dynamics represented by (3)-(7). The equation of state is only needed when coupling with other aspects of the model, such as to compute $T$ needed by the physical parameterizations. The $D\theta_{v}/Dt=0$ approximation can be made independently of the equation of state, but it is more consistent to introduce an approximate Gibbs potential for which $D\theta_{v}/Dt=0$ and use this potential in all equation of state calculations [Thuburn (\APACyear2017)].

For the simulations presented here, we use an equation of state for dry air, $\rho_{d}$, and water vapor represented by specific humidity $q$. We assume $q$ is prognosed by the model’s transport scheme. The Gibbs potential for this equation of state, e.g. \citeA[Appendix A]Vallis2017, leads to $p=\rho R^{*}T$ with $\rho_{d}=\rho(1-q)$, $R^{*}=R+\left(R_{v}-R\right)q$, $c^{*}_{p}=c_{p}+\left(c_{pv}-c_{p}\right)q$ where $R$ and $c_{p}$ are the dry air constants defined above, $R_{v}$ is the ideal gas constants for water vapor and $c_{pv}$ is the specific heat at constant pressure for water vapor. For this equation of state, we have

[TABLE]

and $D\theta/Dt=0$ in the absence of sources or sinks of moisture ($Dq/Dt=0$).

We can introduce an approximate Gibbs potential that is consistent with $D\theta_{v}/Dt=0$ by taking $c_{pv}=(R_{v}/R)c_{p}$. This representing a 14% error in $c_{pv}$. With this $c_{pv}$ approximation, $R^{*}/c_{p}^{*}=R/c_{p}$, so that $R\theta_{v}=R^{*}\theta$ and then $D\theta_{v}/Dt=0$ since $DR^{*}/Dt=(R_{v}-R)Dq/Dt=0$.

2.4 Quasi-Hamiltonian Form

We now give the final form of the continuum equations that will be discretized in § 4. We write the thermodynamic equation in conservation form instead of advective form. We also change the formulation slightly so that it is in quasi-Hamiltonian form (See A) which requires the gradient of the full kinetic energy in the momentum equation. The result is

[TABLE]

with

[TABLE]

where in (15) we have rewritten $\theta_{v}=p/(\rho R\Pi)$ in terms of our prognostic variables.

2.5 Height coordinate

In the height coordinate, our formulation is valid for any general terrain following coordinate as long as $z$ is a monotone function of $s$ and independent of time, such as \citeAGalChen1975 and generalizations that allow for a quicker transition from terrain following to pure height such as \citeASchar2002, Leuenberger2010, klemp2011.

For boundary conditions at the model top and surface, we use a free slip condition ($\dot{s}=0$). Since $\partial\phi/\partial t=0$ throughout the domain, (12) remains a valid equation if we take diagnostic equation for $\dot{s}$ as

[TABLE]

Evaluating this equation at the model top and surface we derive the $w$ boundary condition

[TABLE]

which is combined with (11) to derive boundary conditions for $p$ as well.

2.6 Mass coordinate

For the mass coordinate we follow \citeAArakawaLamb1977,SB81 and define $s$ via

[TABLE]

with $A(s)$ zero near the surface and $B(s)$ zero near the model top. In this coordinate system, the diagnostic equation for $\dot{s}$ is obtained by integrating (14) in the vertical and making use of $\partial\pi/\partial t=B(s)\partial\pi_{s}/\partial t$ to obtain

[TABLE]

For boundary conditions at the model top and surface, we use a free slip condition ($\dot{s}=0$). Since $\partial\phi/\partial t=0$ at the model surface, from (12) we have a boundary condition on $w$,

[TABLE]

which is combined with (11) to derive boundary conditions for $p$ at the surface. At the top of the model we have a free surface with a constant pressure $p=A(s_{\text{top}})p_{0}$.

Following \citeALaprise92, the standard mass coordinate hydrostatic approximation can be obtained by taking $p=\pi$ and replacing (12) with the diagnostic equation for $\phi$ obtained by integrating (15).

2.7 Vertically Lagrangian coordinate

We also consider the vertically Lagrangian coordinate from \citeAsjlin04. We assume the vertical levels move so that there is no flow through the levels ($\dot{s}=0$) and the vertical advection terms are dropped from the equations. In both mass and height options, the diagnostic equation for $\dot{s}$ is no longer valid nor is it needed. Instead, (12) must be used as a prognostic equation for the level positions.

In practice, the floating vertical levels will eventually become too unevenly spaced and every few timesteps we must remap the solution onto a set of more uniform levels. Here we consider energy conservation only during the floating Lagrangian phase. To remap, we typically use a conservative monotone algorithm which will introduce vertical dissipation. To retain energy conservation in the vertically Lagrangian case an energy conserving remap algorithm would have to be adopted.

3 Energy conservation

We first show energy conservation and derive the energy budget in the continuum. We look at the budget term by term to check that we have formulated the equations so that most terms cancel via integration by parts, so that they will also cancel when discretized with a numerical method that has a discrete integration by parts property. As such we must be careful to only perform algebraic steps which will also hold in the discrete system. In this derivation we note a few key terms that require additional numerical properties in order to obtain energy conservation. The analysis in the discrete case will be identical but with integrals replaced by sums and additional care must be taken to account for the vertical staggering.

The kinetic, internal and potential energy of our equations in terrain following coordinates with integration with respect to $ds$ is given by

[TABLE]

with $\hat{p}=p_{\text{top}}$ in the mass coodrinate and $\hat{p}=0$ in the height coordinate. In the mass coordiante, $p_{\text{top}}$ appears in the energy due to the free surface constant pressure boundary condition.

3.1 Potential energy

For the evolution of potential energy in the mass coordinate, we multiply the continuity equation by $\phi$ and the $\phi$ equation by $\frac{\partial\pi}{\partial s}$ and sum,

[TABLE]

The first four terms will integrate to zero, requiring only integration by parts, and the remaining term is the transfer of potential energy to kinetic energy. In the Lorenz staggered case, some care will be needed in the averaging to ensure the first four terms cancel. This equation can be simplified in the Eulerian height coordinate since $\partial\phi/\partial t=0$, but since this term is not zero in the mass coordinate or the vertically Lagrangian height coordinate, we retain this general form.

3.2 Internal energy

For internal energy we have

[TABLE]

Assuming exact time integration, this can be simplified to

[TABLE]

We can eliminate the mass coordinate boundary term if we integrate the second term by parts in the vertical, leading to

[TABLE]

We then take $\Pi$ times the $\Theta$ equation, $\partial p/\partial s$ times the $\phi$ equations to derive

[TABLE]

The first, third and fifth terms represent transfer to kinetic energy and will cancel after integration by parts with similar terms in the kinetic energy equation. In order for the terms involving $\dot{s}$ to cancel, we require

[TABLE]

this relation is needed for the advection terms to be energetically neutral. It can be achieved if care is taken in how the vertical advection of $\Theta$ is defined, or if we are using a vertically Lagrangian coordinate where these terms do not appear.

3.3 Kinetic energy

For kinetic energy, take the sum of $\frac{\partial\pi}{\partial s}\mathbf{u}$ dotted with the $\mathbf{u}$ equation, $\frac{\partial\pi}{\partial s}w$ times the $w$ equation and $(\mathbf{u}^{2}+w^{2})/2$ times the continuity equation. First we consider only the terms related to advection, which should be energetically neutral. After removing two terms that cancel pointwise, we are left with

[TABLE]

The first, second, fifth and sixth terms in the above will cancel when integrating over the sphere via integration by parts. The third, forth, seventh and eight terms cancel in the continuum but rely on additional properties that are not obeyed by general mimetic discretizations. They are responsible for a small lack of conservation in \citeA[§ 4.4]Dubos2015. These terms will cancel via special integral properties satisfied by the SB81 discretization [Simmons \BBA Burridge (\APACyear1981), Taylor (\APACyear2011)]. In § 4.1 we show that this result can be obtained more directly by noting that SB81 also has a discrete derivative product rule.

Now consider the remaining terms involving the pressure gradient

[TABLE]

The first term balances the identical term in the potential energy equation and the remaining terms balance with similar terms in the internal energy equation.

3.4 Energy consistency

For energy consistency, we require that the discrete equations reproduce three key aspects of the continuum equations: (1) the transport terms are energetically neutral, meaning that they vanish to machine precision when integrated over the domain, (2) the diagnostic equations for the change in total potential, internal and kinetic energies hold in exact time integration and (3) the transfer terms in these diagnostic equations cancel to machine precision. These three requirements ensure that the discretization is not introducing any spurious sources of energy and that total energy will be conserved in exact time integration.

We denote integration over the entire domain with a double integral, representing integration in the vertical with respect to $ds$ and integration in the horizontal with respect to an area metric $dA$. The energy transfer terms in the budget equations are denoted

[TABLE]

[TABLE]

After integrating the energy budgets derived in the previous subsection, we obtain

[TABLE]

Energy conservation requires $T_{1}=-S_{1}$, which holds via integration-by-parts.

4 Vertical discretization

We use a \citeALorenz1960 staggering where the prognostic variables are defined at either level interfaces or midpoints. These levels, their indexing and the staggering of the prognostic variables is shown in Fig. 1. The placement of the prognostic variables then leads to the following natural placement for the remaining variables,

[TABLE]

where for convenience we have introduced $\dot{S}=\dot{s}\frac{\partial\pi}{\partial s}$. The formulas for $\dot{S}$ and $\mu$ mix midpoint and interface quantities, and the averaging necessary for energy consistency will be given in the next section. For the level positions, we require the midpoints to be centered within the interfaces

[TABLE]

and denote the level thickness by ${\Delta s}$ with

[TABLE]

and at the boundaries, {\Delta s}_{\text{n+{\nicefrac{{1}}{{2}}}}}={\Delta s}_{n} and ${\Delta s}_{{\nicefrac{{1}}{{2}}}}={\Delta s}_{1}$. We note that the coordinate labeling by $s$ is arbitrary. The solution will depend on the physical position of each level, but not its assigned $s$-value. For simplicity we could relabel $s$ so that $\Delta s=1$ for all levels, as is often done in numerical implementations.

Below we will define the averaging and differential operators using a generic $\phi_{i+{\nicefrac{{1}}{{2}}}}$ variable given at interface points and a generic $p_{i}$ variable given at midpoints with boundary conditions at the model top $p_{{\nicefrac{{1}}{{2}}}}$ and surface $p_{n+{\nicefrac{{1}}{{2}}}}$. We use standard ${\Delta s}$ weighted averaging to map between interface and midpoint quantities. Because of the centered midpoint levels, to average an interface quantities to midpoints we use

[TABLE]

and to average (and extrapolate) midpoint quantities to interfaces, we use

[TABLE]

This simple extrapolation formula is necessary to ensure a discrete averaging-by-parts identity holds. Vertical derivatives are defined by

[TABLE]

with the boundary conditions and one sided differencing used for computing derivatives at the surface and model top

[TABLE]

Note that for the $\partial/\partial s$ operator on midpoint quantities, these boundary conditions must be specified as part of the discretization.

Vertical integrals are approximated via quadrature using a midpoint rule. For midpoint quantities,

[TABLE]

and at interfaces

[TABLE]

where the prime on the sum $\sum^{\prime}$ is used to denote that we multiply the integrand by $\frac{1}{2}$ at the end points $i=0$ and $i=n$.

4.1 Discrete operator identities

With these definitions, our discrete operators obey integration and averaging by parts identities as well as a discrete derivative product rule. For interface quantities, averaging commutes with differentiation

[TABLE]

with the boundary values of $\phi$ used for the $\partial\bar{\phi}/\partial s$ boundary conditions. We have a discrete averaging-by-parts identity

[TABLE]

and a discrete integration-by-parts identity

[TABLE]

Most importantly for our formulation, and difficult to obtain for general discretizations [Salmon (\APACyear2007)], is the pointwise product rule. For interface quantities $a$ and $b$, we have

[TABLE]

For midpoint quantities $c$ and $d$ we have a similar formula, but we have to use unweighted averaging,

[TABLE]

where we have extra factors of ${\Delta s}$ so that the unweighted average of midpoint quantities is written in terms of the weighted average defined in (31).

4.2 SB81 vertical advection

The SB81 vertical advection for operator $\dot{s}\partial/\partial s$ is carefully constructed for midpoint quantities in order to conserve energy. We denote this operator with square brackets $[\cdot]$ and write it in terms of interface quantitiy $\dot{S}$. For midpoint quantities,

[TABLE]

For interface quantities, we introduce a generalization of this operator defined by

[TABLE]

where we have averaged the midpoint quantity $\frac{\partial\pi}{\partial s}$ to interfaces. This generalization appears to contain double averaging but it can also be written in an equivalent form without averages of averages.

5 Lorenz staggered formulation

We now write the equations in the vertically staggered discretization using the bar notation for averaging between midpoints and interfaces using (31) or (30), $\partial/\partial s$ representing (32) and (33) and $[\cdot]$ representing the SB81 transport operator. The discrete equations in their final energy consistent form are

[TABLE]

with boundary condition

[TABLE]

With this staggering, $p$ and $\Pi$ at their midpoint locations are determined from (15) without averaging,

[TABLE]

and averaging is needed for vertical momentum and $\mu$,

[TABLE]

In order for the transport terms to be energetically neutral, we need to use special averaging denoted by a tilde for two variables:

[TABLE]

For $\tilde{\mathbf{u}}$ this is the density weighted average. For $\bar{\theta}_{v}$ this represents a version of (15) at interfaces. This term will not be present in the vertically Lagrangian coordinate. In the Eulerian coordinate, (15) at interfaces is thus treated quite differently then (15) at midpoints. As an alternative one can eliminate (15) altogether by introducing an additional prognostic variable for $p$ or $\Pi$ [Gassmann (\APACyear2013)].

For the mass coordinate, we use a diagnostic integral for $\dot{S}$, $i=1\dots n-1$

[TABLE]

and then compute $\dot{s}$. For the Eulerian height coordinate, we replace the prognostic equation for $\phi$ with the diagnostic equation

[TABLE]

and then compute $\dot{S}$. In the vertically Lagrangian case, $\dot{S}=\dot{s}=0$ and all terms involving these variables are dropped. In the vertically Lagrangian height coordinate, we retain the $\phi$ equation with $\phi$ at the surface and model top held fixed.

We require that the equations for $w$ and $\phi$ hold at all interfaces including the model top and surface. This requirement is used to deduce the consistent boundary conditions for $w$ and $p$. At the surface, $\phi$ is given for both mass and height coordinates and the $\phi$ equation is satisfied by taking $gw=\tilde{\mathbf{u}}\cdot\nabla_{s}\phi_{s}$. This allows us to determine $\partial w/\partial t$ at the surface in terms of $\partial\mathbf{u}/\partial t$, which in turn determines the pressure (and equivalently $\mu$) at the surface using the prognostic equation for $w$. We omit this algebra here as it depends on details of the spatial discretization. For our purposes, we need only know that $p$ and $\mu$ are defined at the surface so that the prognostic $w$ equation holds at the surface.

At the model top, for the height coordinate the procedure is identical and $w$ and $\mu$ are defined so that the $w$ and $\phi$ equations hold at the top most interface. For the mass coordinate, we have a pressure boundary condition and a free surface at the model top. In this case, $\mu$ is well defined and $w$ and $\phi$ are evolved at the model top with their prognostic equations.

We note that we used the SB81 operator for vertical advection of $w$ and the more traditional discretization for vertical advection of $\phi$. The SB81 operator is necessary for energy consistency in the momentum equations, but it does not vanish at the surface. The traditional discretization does vanish at the surface since $\dot{s}=0$, and its use in the $\phi$ equation results in the more natural $gw=\tilde{\mathbf{u}}\cdot\nabla_{s}\phi_{s}$ boundary condition for $w$.

5.1 Discrete energy conservation

We now repeat the energy budget calculations done for the continuum formulation in § 3. For our horizontal discretization we will be using the mimetic collocated spectral element method from \citeATayFour10. In this case, the horizontal terms which cancel after integration by parts will cancel in the discrete case as well and those calculations are not repeated here. For a horizontally staggered case, such as \citeAthuburn2009numerical, more care must be taken to ensure the horizontal integration by parts properties still hold when mixed with vertical staggering [Gassmann (\APACyear2013)].

5.2 Discrete potential energy

We define discrete potential energy at midpoints

[TABLE]

For the budget, we derive the discrete analog of equations (25). We take the continuity equation, average to interfaces and multiply by $\phi$, and then add in the $\phi$ equation multiplied by $\overline{\frac{\partial\pi}{\partial s}}$. Assuming exact time integration, we have

[TABLE]

Similar to the continuum equations, four of these terms will cancel when integrated over the domain after integration by parts and making use of the definition of $\tilde{\mathbf{u}}$. The last term, after averaging by parts, becomes,

[TABLE]

5.3 Discrete internal energy

We use a discrete internal energy

[TABLE]

with $\hat{p}=p_{\nicefrac{{1}}{{2}}}$ in the mass coordinate and $\hat{p}=0$ in the height coordinate. Differentiating with respect to time, and taking $\partial\Pi/\partial t=(R/c_{p})(\Pi/p)\partial p/\partial t$ which will hold in exact time integration as well as in the the limit of small $\Delta t$, and using that $(\partial\hat{p}/\partial t)=0$ in both coordinate systems,

[TABLE]

where we have removed the terms

[TABLE]

since our discrete EOS is $\frac{\partial\phi}{\partial s}=-R(\Theta\Pi/p)_{i}.$ We then apply integration by parts to obtain

[TABLE]

Inserting the discrete prognostic equations,

[TABLE]

The first, third and fifth terms make up the transfer terms which will balance similar terms in the kinetic energy equation. For the vertical advection terms, we integrate one term by parts and use the definition of $\dot{S}$ to obtain

[TABLE]

Ensuring that these two terms cancel is the reason behind our complex formula for $\tilde{\theta}_{v}$.

5.4 Discrete kinetic energy

Starting with $K$ at midpoints,

[TABLE]

summing and applying averaging by parts we have

[TABLE]

We then derive the discrete budget following § 3.3. We dot the discrete equation for $\mathbf{u}$ by $\frac{\partial\pi}{\partial s}\mathbf{u}$ and sum with $\frac{\partial\pi}{\partial s}w$ times the $w$ equation and $(\mathbf{u}^{2}+w^{2})/2$ times the continuity equation. Starting with just the advection terms from (21), we drop four of the terms which cancel after integration by parts in the horizontal and two terms cancel after integration by parts in the vertical, using the SB81 operator for $\dot{s}\partial\mathbf{u}/\partial s$ and the product rule identity (39). The remaining advection terms are

[TABLE]

Using the definition of $\tilde{\mathbf{u}}$ from (48) the first two terms cancel. Expanding the SB81 operator and applying integration-by-parts we are left with

[TABLE]

To show these two terms cancel, we apply the SB81 product rule, average by parts, apply (36) and then integrate by parts making use of $\dot{S}=0$ at the boundaries,

[TABLE]

Finally, we are left with the discrete version of the transfer terms from (22),

[TABLE]

5.5 Discrete energy consistency

We now write the discrete integral energy budgets. Let the sum over $i$ represent the vertical integral as above. In the horizontal, we let the sum over $j$ represent the discrete quadrature for integration over spherical $s$ surfaces associated with the mimetic horizontal discretization.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

After integrating the energy budgets derived in the previous subsection, we obtain

[TABLE]

To obtain energy conservation, we note that $S_{1}=-T_{1}$ after integration by parts and $S_{2}=T_{2}$ after averaging-by-parts. To show $S_{3}=T_{3}$ we average by parts and then expand $\mu$ and $\tilde{\mathbf{u}}$. Ensuring $S_{3}=T_{3}$ is the motivation behind the averaging of $\mu\nabla_{s}\phi$ in the momentum equation. This is the only term in our formulation that contains double averaging (since $\mu$ contains an average of $\frac{\partial\pi}{\partial s}$). This is unavoidable given our definition of $\tilde{\mathbf{u}}$ needed to obtain energy neutral transport.

6 Numerical implementation

We have implemented the discretization described here in an open source version of HOMME. HOMME previously contained a shallow water and hydrostatic primitive equation model and we use HOMME-NH to refer to the new nonhydrostatic model. For the horizontal directions, HOMME-NH relies on the existing HOMME infrastructure including the mimetic spectral element discretization, monotone and conservative tracer transport [Guba \BOthers. (\APACyear2014)], and support for fully unstructured arbitrary quadrilateral grids [Guba \BOthers. (\APACyear2014b)]. Here we will be running the model on the quasi-uniform equi-angle cubed sphere grid [Ronchi \BOthers. (\APACyear1996)]. Within each spectral element we use 4th order accurate polynomial basis functions of degree $p=3$.

HOMME-NH implements (42)–(47) using method of lines. The horizontal divergence, gradient and curl operators on $s$-surfaces are evaluated with existing HOMME subroutines, and the vertical terms are computed as written. All the tendency terms are computed locally within each element and then we apply a standard finite element assembly step. Because of the spectral element method’s diagonal mass matrix, this assembly is explicit and can be computed in each element with only nearest neighbor information. For our initial HOMME-NH implementation we use the mass coordinate for with both Eulerian and vertically Lagrangian options. For the vertically Lagrangian option, we typically remap back to the initial reference levels every three timesteps using parabolic splines [Zerroukat \BOthers. (\APACyear2006)].

For timestepping, we use a HEVI approach [Satoh (\APACyear2002)] discretized with a Runge-Kutta IMEX method. The HEVI/IMEX combination has proven attractive for atmospheric models, e.g. \citeAUllrichJablonowski2012,GKC2013,Weller2013,Lock2014. HOMME-NH supports a variety of IMEX methods from the ARKODE package [Gardner \BOthers. (\APACyear2017), Gardner \BOthers. (\APACyear2018)] of the SUNDIALS library [Hindmarsh \BOthers. (\APACyear2005)]. Here we use a 5 stage 2nd order accurate high-CFL method derived in \citeASteyer2019. In the HEVI approach, the vertically implicit component can be solved independently in each vertical column. We use Newton iterations with the Jacobian computed analytically and inverted by direct LU factorization. We use a time-split approach to couple (42)–(47) with hyperviscosity and physical parameterizations following \citeADennis12.

HOMME-NH is further described in [Ullrich \BOthers. (\APACyear2017)] and results from a splitting supercell thunderstorm test case are given in [Zarzycki \BOthers. (\APACyear2019)].

7 Numerical results

We use HOMME-NH to verify the three energy consistent properties from § 3.4. Properties (1) and (3) depend on the mimetic properties of the discretization. These are verified through unit tests for all the terms in the energy budgets which should cancel when integrated over the domain. These tests are useful for ensuring all differential operators and quadrature rules are properly coded and are not presented here. To establish property (2) we verify that equations (60)-(62) are satisfied to within time truncation error when running HOMME-NH without forcing or explicit dissipation. Since these equations would be exact in exact time integration, their residual is a measure of the dissipation introduced by the time stepping algorithm. In an energy consistent discretization, this dissipation should converge to zero with $\Delta t$.

We make use of a moist baroclinic instability test case from the Dynamical Core Model Intercomparison Project (DCMIP) [Ullrich \BOthers. (\APACyear2017), Ullrich \BOthers. (\APACyear2018b)]. This test case is based on \citeAullrich2014 with the addition a prognostic variable for specific humidity and the \citeAReed2012 simple physics package. The physics includes large-scale condensation, surface fluxes and boundary layer turbulence. We use a low resolution grid with $4^{\circ}$ average grid spacing at the Equator and 30 vertical levels. This low resolution ensures energy conservation holds even in the presense of large spatial truncation errors.

To analyze the HOMME-NH discretization with a realistic flow containing energy in all resolved scales we first run the full model for 15 days. The solution at $t_{0}=$15 days is shown in Figures 2 and 3. We then disable the forcing and all explicit model dissipation and run to time $t_{1}=t_{0}+$ 2 hours. We use the two hour period $[t_{0},t_{1}]$ to study the HOMME-NH solution to the adiabatic component of the model, equations (42)–(47).

The energy conservation of the discretization as a function of $\Delta t$ is shown in Fig.4. We plot the relative change in total energy over $[t_{0},t_{1}]$

[TABLE]

for simulations with $\Delta t$ ranging from close to the CFL stability limit (1200s) down to 18s. With a second order accurate time stepping algorithm we expect second order convergence in the energy, as shown.

Next we examine the residual error for (60)-(62) instantaneously at time $t_{1}$ given by

[TABLE]

The residuals $R_{p},R_{I}$ and $R_{k}$ represent unphysical changes in potential, internal and kinetic energy in HOMME-NH as a function of $\Delta t$. For simplicity we estimate $\partial/\partial t$ with a first order accurate approximation so that $R_{p},R_{I}$ and $R_{k}$ can be computed in the code at each timestep with only data from that timestep. In Fig.5 we plot $|R_{P}|,|R_{I}|$ and $|R_{K}|$, with $\Delta t$ again ranging from the CFL limit (1200s) to 18s. From the figure we see that $R_{p}$ and $R_{I}$ decrease to zero as $O(\Delta t)$, as expected due to our first order approximation. The value of $R_{k}$ is 3 orders of magnitude smaller than $R_{p}$ and $R_{I}$, and decreases to zero at an even faster rate until reaching machine precision. For this problem, $R_{K}$ and $R_{I}$ are always negative, $R_{P}$ is always positive, $K=8.8\times 10^{5}\,J/m^{2}$, $I=1.8\times 10^{9}\,J/m^{2}$ and $P=7.2\times 10^{8}\,J/m^{2}$. At $\Delta t=300s$, $R_{k}=-1.4\times 10^{-8}\,W/m^{2}$ corresponds to change in the 12’th digit of $K$. These results verify the HOMME-NH implementation of (42)–(47), confirming that the spatial discretization does not have any spurious sources of energy and that the only changes in energy are due to the timestepping algorithm.

Most simulations with HOMME-NH are run with timesteps close to the CFL limit. In this regime, we see that the dissipation from time truncation is $\sim 0.02\,\mbox{W/m}^{2}$ at day 15. To compare this number with dissipation from the non-adiabatic terms, we note that at day 15 before the explicit dissipation was turned off, the kinetic energy dissipation from hyperviscsoity and the vertical remap is $0.24\,\mbox{W/m}^{2}$ and $0.06\,\mbox{W/m}^{2}$ respectively.

8 Conclusions

We have presented a discretization of the nonhydrostatic equations in terrain following coordinates with a Lorenz staggered vertical discretization. The discretization is energy consistent in exact time integration when coupled with a mimetic horizontal discretization, meaning the discrete diagnostic equations for potential, internal and kinetic energy are consistent and total energy is conserved to time truncation error. It supports both mass and height coordinates, which differ only in the top of model boundary condition and the diagnostic equation for $\dot{s}$. Vertical transport is kept in advective form, leading to unified support for both Eulerian or vertically Lagrangian coordinates. We use virtual potential temperature and prognose total mass to support general moist equations of state. The mass coordinate formulation has been implemented in HOMME-NH, with energy consistency verified using a moist baroclinic instability test case and a Runge-Kutta IMEX timestepping algorithm. In the design of HOMME-NH, we have found energy consistency to be a useful guiding principal that enhances the stability of the discretization. The ultimate goal is not exact energy conservation, but to preserve the Hamiltonian structure of the equations and ensure all changes in energy are through explicitly added dissipative mechanisms. In HOMME, these mechanisms include hyperviscosity in the horizontal, monotone remap in the vertical and dissipation inherent in the Runge-Kutta method. HOMME-NH is currently being evaluated with a suite of test cases from the DCMIP. Future work will be to add support for the height coordinate and coupling with the Energy Exascale Earth System Model’s suite of atmospheric parameterizations.

Appendix A Quasi-Hamiltonian form

Following \citeADubosTort2014, we show that (10)-(14) is in quasi-Hamiltonian form plus the addition of an energetically neutral vertical transport term $\mathcal{V}\cdot\dot{s}$. Let

[TABLE]

and take

[TABLE]

with $P,I,K,dA$ and $ds$ defined as in § 3. To find the functional derivatives of $\mathcal{H}$ we differentiate with respect to time as in § 3 to obtain

[TABLE]

and thus

[TABLE]

Using these functional derivatives, the equations of motion can be re-written as

[TABLE]

This can be seen to fall into the quasi-Hamiltonian form

[TABLE]

with $\mathcal{J}$ a skew symmetric operator defined via

[TABLE]

where $\eta_{s}=\nabla_{s}\times\mathbf{u}+f\hat{k}$ and

[TABLE]

We call this quasi-Hamiltonian form since no attempt has been made to prove the Jacobi identity for $\mathcal{J}$. Energy is conserved due the skew-symmetry of $\mathcal{J}$ and the condition

[TABLE]

where

[TABLE]

is the adjoint of $\mathcal{V}$. The (74) condition is a consequence of the vertical relabeling symmetry of the flow and implies the vertical transport is energetically neutral for any velocity $\dot{s}$ with boundary conditions $\dot{s}=0$.

A general mimetic discretization will preserve the skew-symmetry of $\mathcal{J}$, but additional properties are needed to obtain a discrete version of (74). To see this, we insert the functional derivatives into

[TABLE]

to obtain

[TABLE]

The product rule for $\partial/\partial s$ immediately gives cancellation for the $\mathbf{u}$, $w$ and $\phi$ terms, leaving only

[TABLE]

Now recalling that $\frac{\partial\pi}{\partial s}\mu=\frac{\partial p}{\partial s}$ and

[TABLE]

which holds in the continuum due to the definition of $\Pi$ and $\theta_{v}$. We can thus obtain a discrete version of (74) with the the SB81 product rule (38) and the $\theta_{v}$ special averaging in (48).

Acknowledgements.

HOMME-NH is part of the open source E3SM project with source code publically available from https://github.com/E3SM-Project/E3SM. Instructions for configuring and compiling standalone HOMME-NH are in the README Files in components/homme. Once the model has been configured, the checkE.sh script in dcmip_tests/dcmip2016_test1_baroclinic_wave/theta-l will generate all necessary input data and run the simulations described in this work. This research was supported as part of the Energy Exascale Earth System Model (E3SM) project, funded by the U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government. Christopher Eldred was supported by the French National Research Agency through contract ANR-14-CE23-0010 (HEAT).