A Mixed Mimetic Spectral Element Model of the 3D Compressible Euler Equations on the Cubed Sphere

TL;DR

This paper introduces a novel 3D compressible Euler equations model on the cubed sphere using mixed mimetic spectral elements, enabling energy conservation and efficient vertical-horizontal decoupling for atmospheric simulations.

Contribution

It presents a new discretization approach combining mimetic spectral elements with dimensional splitting, improving energy conservation and computational efficiency in atmospheric modeling.

Findings

Successfully validated with baroclinic instability tests

Accurately modeled non-hydrostatic gravity waves

Achieved energy exchanges consistent with physical laws

Abstract

A model of the three-dimensional rotating compressible Euler equations on the cubed sphere is presented. The model uses a mixed mimetic spectral element discretization which allows for the exact exchanges of kinetic, internal and potential energy via the compatibility properties of the chosen function spaces. A Strang carryover dimensional splitting procedure is used, with the horizontal dynamics solved explicitly and the vertical dynamics solved implicitly so as to avoid the CFL restriction of the vertical sound waves. The function spaces used to represent the horizontal dynamics are discontinuous across vertical element boundaries, such that each horizontal layer is solved independently so as to avoid the need to invert a global 3D mass matrix, while the function spaces used to represent the vertical dynamics are similarly discontinuous across horizontal element boundaries, allowing for the serial solution of the vertical dynamics independently for each horizontal element. The model is validated against standard test cases for baroclinic instability within an otherwise hydrostatically and geostrophically balanced atmosphere, and a non-hydrostatic gravity wave as driven by a temperature perturbation.

Equations

$$\displaystyle\frac{\partial\bm{u}}{\partial t}+(\bm{\omega}+\bm{f})\times\bm{u}+\nabla\Bigg{(}\frac{1}{2}\|u\|^{2}+gz\Bigg{)}+\theta\nabla\Pi$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{u})$$

$$\frac{\partial \rho \theta }{\partial t}+\nabla \cdot (\rho \theta \mathbf{u})$$

$$\Pi$$

$$\theta$$

$$p$$

$$\Pi$$

$$\theta$$

$$\frac{1}{\rho}\nabla p\stackrel{{\scriptstyle\eqref{eq:pressure_definition}}}{{=}}\frac{RT}{p}\nabla p\stackrel{{\scriptstyle\eqref{eq:exener_definition}+\eqref{eq:potential_temperature_definition}}}{{=}}\frac{R\theta}{p}\Bigg{(}\frac{p}{p_{0}}\Bigg{)}^{\frac{R}{c_{p}}}\nabla p=\frac{R\theta}{p_{0}^{R/c_{p}}}p^{R/c_{p}-1}\nabla p=\frac{R\theta}{p_{0}^{R/c_{p}}}\frac{c_{p}}{R}\nabla p^{R/c_{p}}=c_{p}\theta\nabla\Bigg{(}\frac{p}{p_{0}}\Bigg{)}^{\frac{R}{c_{p}}}\stackrel{{\scriptstyle\eqref{eq:exener_definition}+\eqref{eq:potential_temperature_definition}}}{{=}}\theta\nabla\Pi.$$

$${w|}_{z=0}={w|}_{z=z^{\mathrm{top}}}=0,$$

$$\frac{\partial\Pi}{\partial z}\Bigg{|}_{z=0}=\frac{\partial\Pi}{\partial z}\Bigg{|}_{z=z^{top}}=0.$$

$$K:=\frac{1}{2}⟨\mathbf{u},\rho \mathbf{u}⟩=\frac{1}{2}{\int }_{\Omega }\rho \Vert u{\Vert }^2,$$

$$⟨f,g⟩:={\int }_{\Omega }fg\mathrm{d}\Omega ,$$

$$⟨\mathbf{u},\mathbf{v}⟩:={\int }_{\Omega }\mathbf{u}\cdot \mathbf{v}\mathrm{d}\Omega ,$$

$$\frac{\partial K}{\partial t}=-⟨g,\rho w⟩-⟨\rho \mathbf{u},\theta \nabla \Pi ⟩,$$

$$P:=⟨\rho ,gz⟩={\int }_{\Omega }\rho gz\mathrm{d}\Omega ,$$

$$\frac{\partial P}{\partial t}=⟨gz,\frac{\partial \rho }{\partial t}⟩\stackrel{\text{eqref}eq::dens}{=}-⟨gz,\nabla \cdot (\rho \mathbf{u})⟩=⟨g,\rho w⟩,$$

$$I:={\int }_{\Omega }{c}_v\rho T\mathrm{d}\Omega \stackrel{\text{eqref}eq:potentia{l}_{t}emperatur{e}_{d}efinition}{=}{\int }_{\Omega }\frac{c_v}{c_p}\rho \theta \Pi \mathrm{d}\Omega \stackrel{\text{eqref}eq:p{i}_{f}ro{m}_{t}het{a}_{r}ho}{=}{\int }_{\Omega }{c}_v\rho \theta {\left(\frac{R\rho \theta }{p_0}\right)}^{\frac{R}{c_v}}\mathrm{d}\Omega ={\int }_{\Omega }{c}_v{\left(\frac{R}{p_0}\right)}^{\frac{R}{c_v}}{\left(\rho \theta \right)}^{\frac{c_p}{c_v}}\mathrm{d}\Omega .$$

$$\frac{\partial I}{\partial t}=-⟨\nabla \cdot (\rho \theta \mathbf{u}),\Pi ⟩=⟨\rho \mathbf{u},\theta \nabla \Pi ⟩.$$

$$\mathcal{H}:=K+P+I={\int }_{\Omega }\frac{1}{2}\rho u^2\mathrm{d}\Omega +{\int }_{\Omega }\rho gz\mathrm{d}\Omega +{\int }_{\Omega }\frac{c_v}{c_p}\Theta \Pi \mathrm{d}\Omega .$$

$$\bm{\mathsf{a}}:=\left[\begin{array}[]{ccc}\bm{u}&\rho&\Theta\end{array}\right]^{\top}\,,$$

$$\bm{\mathsf{h}}:=\left[\begin{array}[]{ccc}\bm{U}&\Phi&\Pi\end{array}\right]^{\top}\,,$$

$$\bm{\mathsf{B}}:=\left[\begin{array}[]{ccc}-\bm{q}\times(\cdot)&-\nabla(\cdot)&-\theta\nabla(\cdot)\\ -\nabla\cdot(\cdot)&0&0\\ -\nabla\cdot(\theta\cdot)&0&0\\ \end{array}\right]\,,$$

$$\frac{\partial \mathsf{a}}{\partial t}=\mathsf{B}\mathsf{h}.$$

$$\frac{\delta \mathcal{H}}{\delta \mathbf{u}}=\rho \mathbf{u}=\mathbf{U},\frac{\delta \mathcal{H}}{\delta \rho }=\frac{1}{2}{u}^2+gz=\Phi ,\frac{\delta \mathcal{H}}{\delta \Theta }=\Pi ,$$

$$\mathsf{h}=\frac{\delta \mathcal{H}}{\delta \mathsf{a}}.$$

$$\frac{\partial \mathsf{a}}{\partial t}=\mathsf{B}\frac{\delta \mathcal{H}}{\delta \mathsf{a}}.$$

$$\frac{\partial\mathcal{H}}{\partial t}=\frac{\delta\mathcal{H}}{\delta\bm{\mathsf{a}}}\cdot\frac{\partial\bm{\mathsf{a}}}{\partial t}\stackrel{{\scriptstyle\eqref{eq:da_dt_B_H}}}{{=}}\frac{\delta\mathcal{H}}{\delta\bm{\mathsf{a}}}\cdot\left(\bm{\mathsf{B}}\frac{\delta\mathcal{H}}{\delta\bm{\mathsf{a}}}\right)=0,$$

$$\frac{\delta \mathcal{H}}{\delta \mathbf{a}}\cdot \frac{\partial \mathbf{a}}{\partial t}={\int }_{\Omega }\sum _{i=1}^3\frac{\delta \mathcal{H}}{\delta {\mathsf{a}}_i}\frac{\partial {\mathsf{a}}_i}{\partial t}\mathrm{d}\Omega .$$

$$l_i^p({\xi }_j)={\delta }_{i,j},i,j=0,\dots ,p,$$

$$l_i^p(\xi )=\prod _{\begin{array}{c}k=0\\ k\ne i\end{array}}^p\frac{\xi -{\xi }_k}{{\xi }_i-{\xi }_k}.$$

$$q_h(\xi ):=\sum _{i=0}^p{q}_i{l}_i^p(\xi ).$$

$$e_i^p(\xi ):=-\sum _{k=0}^{i-1}\frac{\mathrm{d}{l}_k^p(\xi )}{\mathrm{d}\xi },i=1,\dots ,p.$$

Full text

††journal: Journal of Computational Physics

1 Introduction

Mimetic finite element families are an appealing choice for the discretization of geophysical flow problems. This is on account of their capacity to preserve both conservation laws and leading order balance relations in the discrete form [1, 2, 3, 4], due to the compatibility properties of the chosen function spaces, as well as their ability to represent complex geometries such as the surface of the sphere [5, 6, 7].

The use of mimetic discretizations to represent the solution variables, and the adjoint properties of the differential operators implied by those spaces, allows for the conservation of energy via the exact balance of energetic exchanges, as well as the orthogonality of vorticity evolution to those exchanges [2, 3, 4]. By satisfying exactly the balance between energetic exchanges, it is hoped that these methods may improve the statistical behaviour of climate simulations over long time scales by mitigating against internal biases in the representation of dynamical processes.

Several mimetic finite volume models of the compressible Euler equations on the sphere have previously been presented [8, 9, 10], and the Raviart-Thomas family of compatible finite elements has been chosen to form the basis of the LFric atmospheric model [11, 12]. Collocated spectral element methods, which in contrast to the above methods represent all solution variables on the same function space, are a popular choice for the simulation of non-hydrostatic atmospheric flows using both continuous [13] and discontinuous Galerkin [14, 15] formulations. These models typically satisfy some mimetic properties, such as the divergence theorem and the adjoint property between the gradient and divergence operators [16], however the HOMME spectral element model [17, 18] is capable of preserving a more complete set of mimetic properties so as to satisfy both the divergence and circulation theorems through the careful consideration and construction of the metric terms for the various differential operators.

In the present formulation we make use of the mixed mimetic spectral element method [19, 20], a compatible family of function spaces with spectral error convergence. We use this method to build on previous work on the rotating shallow water equations [3, 5] in order to develop a solver for the three-dimensional rotating compressible Euler equations on the cubed sphere. In contrast to collocated methods, weak formulations using mixed function spaces provide a clearer mathematical formalism for the analysis and construction of schemes which preserve internal properties of the governing equations. Like other weak form mimetic discretizations such as mixed finite elements [1, 2, 7] and mimetic Galerkin differences [4], the mixed mimetic spectral element method satisfies by construction an exact mapping between function spaces via the various differential operators, as well as adjoint relationships between strong form mappings to higher function spaces, and weak form mappings to lower spaces, as defined within an appropriate Hilbert space.

The remainder of this article proceeds as follows: In Section 2 the rotating compressible Euler equations, and their energetic properties will be introduced in the continuous form. Section 3 will provide a brief introduction to the mixed mimetic spectral element method. Readers are referred to references therein for more detailed discussions. Section 4 will discuss the construction of discrete function spaces built off the mixed mimetic spectral element method required to solve the compressible Euler equations, as well as the use of those spaces to ensure consistent energetic exchanges in the discrete form and the associated metric transformations for these spaces. The details of the time stepping scheme, including the implicit vertical solver will be discussed in Section 5, and the results for a standard baroclinic instability test case and a high resolution gravity wave will be presented in Section 6. Section 7 will discuss the conclusions of this work and the future directions we intend to pursue.

2 The rotating compressible Euler equations

The compressible Euler equations for a shallow atmosphere may be expressed as [21, 9]

[TABLE]

where $\bm{u}=u\bm{e}_{\lambda}+v\bm{e}_{\varphi}+w\bm{e}_{z}$ are the zonal, meridional and vertical velocity components respectively, $\rho$ is the density, $\bm{f}=f\bm{e}_{z}$ is the Coriolis term, $g$ is the acceleration due to gravity, $\theta$ is the potential temperature, and $\Pi$ is the Exner pressure (including the specific heat at constant pressure). The last two are defined with respect to the standard thermodynamic variables of temperature, $T$, and pressure, $p$, as

[TABLE]

For these identities we used $c_{p}$ for the specific heat at constant pressure, $p_{0}$ for the reference pressure, and $R$ for the ideal gas constant. We may remove the direct dependence on pressure from the system by simply substituting expression (2c) for pressure into (2a) and (2b), obtaining

[TABLE]

where $c_{v}=c_{p}-R$ is the specific heat at constant volume.

The potential temperature/Exner pressure form of the pressure gradient term, $\theta\nabla\Pi$, in (1a) is equivalent to the standard density/pressure form since

[TABLE]

One advantage of the Exner pressure/potential temperature representation of the thermodynamics is the formulation of the temperature equation in flux form (1c), which allows us to exploit the adjoint relationship between gradient and divergence in the mimetic framework in order to preserve energetic exchanges.

To obtain a closed system for the solution of the compressible Euler equations, (1) and (3) must be supplemented by Dirichlet and Neumann boundary conditions. The following identities are imposed as Dirichlet boundary conditions on the $z$-component of velocity, $w$

[TABLE]

where $z^{\mathrm{top}}$ corresponds to the $z$-coordinates of the top boundary of the domain. For Neumann boundary conditions the following identities are imposed

[TABLE]

Note that in this formulation we have invoked the shallow atmosphere approximation, for which gravity is constant throughout the fluid column, the height of the fluid column is negligible with respect to the earth’s radius, and the horizontal components of the Coriolis term are omitted [22].

2.1 Energetics

Before introducing the discrete form of the Euler equations, we analyse the energetics of the continuous system. This will help to guide our choice of function spaces for the various solution variables for the discrete form.

2.1.1 Kinetic, potential, and internal energy

The kinetic energy, $K$, is defined as

[TABLE]

where $\|u\|:=\langle\bm{u},\bm{u}\rangle$, and $\langle\cdot,\cdot\rangle$ is the $L^{2}$ inner product given as usual as

[TABLE]

for scalar fields, and as

[TABLE]

for vector fields.

The time variation of kinetic energy is obtained by summing the $L^{2}$ inner product, between the momentum equation, (1a), and $\rho\bm{u}$, and between the continuity equation, (1b), and $\frac{1}{2}\|u\|^{2}$

[TABLE]

where again $w$ is the $z$-component of the velocity field, $\bm{u}$.

The potential energy, $P$, is given by

[TABLE]

and its time derivative follows directly

[TABLE]

where we have used integration by parts on the last identity and assumed periodic boundary conditions in the horizontal directions, together with homogeneous boundary conditions for the vertical component of the velocity field, (5).

The internal energy, $I$, is defined as

[TABLE]

After some manipulation, the time variation of internal energy is given by

[TABLE]

where integration by parts was used on the last identity, together with homogeneous boundary conditions for $\bm{u}$ and periodic boundary conditions on the horizontal directions.

2.1.2 Conservation of total energy

Following [9], the total energy of the system, $\mathcal{H}$, is given as the sum of kinetic, $K$, potential, $P$, and internal, $I$, energy

[TABLE]

where we used $\Theta:=\rho\theta$.

For the proof of conservation of total energy $\mathcal{H}$, (15), first consider the column vectors

[TABLE]

and

[TABLE]

where $\bm{U}:=\rho\bm{u}$, and $\Phi:=\frac{1}{2}u^{2}+gz$. Introducing the skew-symmetric operator

[TABLE]

where $\bm{q}=(\bm{\omega}+\bm{f})/\rho$ is the potential vorticity, the original prognostic equations, (1a)-(1c), may be rewritten as

[TABLE]

Note now that the variational derivatives of $\mathcal{H}$ with respect to the prognostic variables $\bm{u}$, $\rho$, and $\Theta$, are

[TABLE]

and therefore

[TABLE]

Substituting (21) into (19) yields

[TABLE]

Conservation of total energy follows directly since [23]

[TABLE]

where the last identity follows from the skew-symmetry of $\bm{\mathsf{B}}$. Note that here the dot product involves not only a summation over the elements of the vectors but also an integration over $\Omega$, e.g.

[TABLE]

3 Mimetic polynomial basis functions

3.1 1D mimetic polynomial function spaces

The mixed mimetic spectral element method is built off two types of one-dimensional polynomials: one associated with nodal interpolation, and the other with integral interpolation (histopolation) [24, 19]. Subsequently, these two types of polynomials will be combined to generate the family of three-dimensional basis functions used to discretize the system.

3.1.1 Nodal polynomial basis functions

Consider the canonical interval $I=[-1,1]\subset\mathbb{R}$ and the Legendre polynomials, $L_{p}(\xi)$ of degree $p$ with $\xi\in I$. The $p+1$ roots, $\xi_{i}$, of the polynomial $(1-\xi^{2})\frac{\mathrm{d}L_{p}}{\mathrm{d}\xi}$ are called Gauss-Lobatto-Legendre (GLL) nodes and satisfy $-1=\xi_{0}<\xi_{1}<\dots<\xi_{p-1}<\xi_{p}=1$. Let $l^{p}_{i}(\xi)$ be the Lagrange polynomial of degree $p$ through the GLL $\xi_{i}$, such that

[TABLE]

where $\delta_{i,j}$ is the Kronecker delta. The explicit form of these Lagrange polynomials is given by

[TABLE]

Let $q_{h}(\xi)$ be a polynomial of degree $p$ defined on $I=[-1,1]$ and $q_{i}=q_{h}(\xi_{i})$, then the expansion of $q_{h}(\xi)$ in terms of Lagrange polynomials is given by

[TABLE]

Because the expansion coefficients in (26) are given by the value of $q_{h}$ in the nodes $\xi_{i}$, we refer to this interpolation as a nodal interpolation and we will denote the Lagrange polynomials in (25) by nodal polynomials.

3.1.2 Histopolant polynomial basis functions

In addition to the nodal basis functions defined above, a second set of basis functions is required in order to discretize integral based quantities across edges. Together these nodal and histopolant basis functions will then be used to construct a compatible family of finite element function spaces in multiple dimensions. Using the nodal polynomials we can define this set of basis polynomials, $e^{p}_{i}(\xi)$, as

[TABLE]

These polynomials $e^{p}_{i}(\xi)$ have polynomial degree $p-1$ and satisfy,

[TABLE]

Using (27) the integral of $e^{p}_{i}(\xi)$ becomes [24, 19]

[TABLE]

Let $g_{h}(\xi)$ be a polynomial of degree $(p-1)$ defined on $I=[-1,1]$ and $g_{i}=\int_{\xi_{i-1}}^{\xi_{i}}g_{h}(\xi)\,\mathrm{d}\xi$, then its expansion in terms of the polynomials $e_{i}^{p}(\xi)$ is given by

[TABLE]

Because the expansion coefficients in (29) are the integral values of $g_{h}(\xi)$, we denote the polynomials in (27) by histopolant polynomials and refer to (29) as histopolation. It can be shown, [24, 19], that if $q_{h}(\xi)$ is expanded in terms of nodal polynomials, as in (26), then the expansion of its derivative $\frac{\mathrm{d}q_{h}(\xi)}{\mathrm{d}\xi}$ in terms of histopolant, or edge polynomials is

[TABLE]

where $\mathsf{E}^{1,0}_{i,j}$ are the coefficients of the $p\times(p+1)$ matrix $\bm{\mathsf{E}}^{1,0}$ for the one dimensional case, hereafter referred to as an incidence matrix (see A for details and [3, 25, 42] for an extensive discussion). The following identity holds (Commuting property)

[TABLE]

For an example of the one-dimensional basis polynomials corresponding to $p=4$, see Figure 1.

3.2 3D mimetic polynomial function spaces

A fundamental element in the proposed discretization for the compressible Euler equations, (1), is the de Rham sequence of function spaces in the domain $\Omega\subset\mathbb{R}^{3}$

[TABLE]

where, as usual, the space $H^{1}(\Omega)$ represents square integrable functions over $\Omega$ whose gradient is also square integrable, the function spaces $H(\mathrm{curl},\Omega)$ and $H(\mathrm{div},\Omega)$ contain square integrable vector fields over $\Omega$ with square integrable curl and divergence, respectively, and the function space $L^{2}(\Omega)$ contains square integrable functions.

More specifically, this work relies on approximating polynomial spaces $\mathcal{P}_{h}(\Omega)\subset H^{1}(\Omega)$, $\mathcal{W}_{h}(\Omega)\subset H(\mathrm{curl},\Omega)$, $\mathcal{U}_{h}(\Omega)\subset H(\mathrm{div},\Omega)$, and $\mathcal{Q}_{h}(\Omega)\subset L^{2}(\Omega)$ such that

[TABLE]

These 3-dimensional polynomial function spaces may be constructed from tensor products of the 1-dimensional function spaces presented in Section 3.1. Moreover, each of these polynomial function spaces has an associated finite set of basis functions $\epsilon^{\mathcal{P}}_{i}$, $\bm{\epsilon}^{\mathcal{W}}_{i}$, $\bm{\epsilon}^{\mathcal{U}}_{i}$, and $\epsilon^{\mathcal{Q}}_{i}$, such that

[TABLE]

As previously discussed, e.g. [25], these basis functions are given by

[TABLE]

where for compactness, the subscripts $+$ and $-$ mean addition and subtraction of 1, e.g. $p_{+}:=p+1$ and $i_{-}:=i-1$.

Moreover, the basis functions satisfy the following identities, see for example [25]:

[TABLE]

where $\bm{\mathsf{E}}^{1,0}$, $\bm{\mathsf{E}}^{2,1}$, and $\bm{\mathsf{E}}^{3,2}$ are the incidence matrices corresponding to the discrete versions of the differential operators grad, curl, and div, respectively (see Appendix for details).

4 Numerical discretization

4.1 Splitting into horizontal and vertical contributions

Consider the following splitting into the horizontal, $\bm{u}_{\parallel}$, and vertical, $\bm{u}_{\perp}$, components of the velocity field $\bm{u}=u\bm{e}_{\lambda}+v\bm{e}_{\varphi}+w\bm{e}_{z}$

[TABLE]

Moreover, let $\nabla_{\parallel}$ and $\nabla_{\perp}$ represent the horizontal and vertical components of the gradient operator of a scalar field $\rho$

[TABLE]

In a similar way, $\nabla_{\parallel}\times$ and $\nabla_{\perp}\times$ represent, respectively, the horizontal and vertical components of the curl of a vector field $\bm{u}=u\bm{e}_{\lambda}+v\bm{e}_{\varphi}+w\bm{e}_{z}$

[TABLE]

Note that we have assumed the shallow atmosphere approximation of constant radius, $r$, in the above expressions. From (37) and (38) follows directly that

[TABLE]

With (36) and (38) it is possible to rewrite the definition of vorticity, $\bm{\omega}:=\nabla\times\bm{u}$, as

[TABLE]

Using (36), (37), (38), and (39), it is possible to split the compressible Euler equations, (1), into horizontal and vertical components

[TABLE]

where, as in (39), $\bm{\omega}_{\parallel,\parallel}:=\nabla_{\parallel}\times\bm{u}_{\parallel}$ and $\bm{\omega}_{\parallel,\perp}:=\nabla_{\parallel}\times\bm{u}_{\perp}$.

The same splitting into horizontal and vertical components may be done for the basis functions, Section 3.2,

[TABLE]

and

[TABLE]

Recalling the definition of $\bm{\epsilon}^{\mathcal{W}}_{j}$ and $\bm{\epsilon}^{\mathcal{U}}_{j}$, Section 3.2, we can explicitly write their horizontal and vertical components as

[TABLE]

Using this splitting of the basis functions we can also split the discrete function spaces such that (34) becomes

[TABLE]

where

[TABLE]

Remark 1

An important point relevant in the derivations that follow is that

[TABLE]

4.2 Unsplit discretization

To discretize the unsplit form of the compressible Euler equations, (1), we first introduce the weak form: Given a domain $\Omega\subset\mathbb{R}^{3}$ and a Coriolis term $\bm{f}=f\bm{e}_{z}\in H(\mathrm{curl},\Omega)$, find $\bm{u},\bm{U},\bm{F},\bm{P}\in H(\mathrm{div},\Omega)$, $\bm{\omega}\in H(\mathrm{curl},\Omega)$, and $\Pi,\rho,\theta,\Theta\in L^{2}(\Omega)$ for the prognostic equations

[TABLE]

as well as the associated diagnostic equations

[TABLE]

Consider now the domain $\Omega\subset\mathbb{R}^{3}$ and its tessellation $\mathcal{T}(\Omega)$ consisting of $M$ arbitrary quadrilaterals (curved), $\Omega_{m}$, with $m=1,\dots,M$. We assume that all quadrilateral elements $\Omega_{m}$ can be obtained from a map $\Phi_{m}:(\xi,\eta,\zeta)\in I^{3}\mapsto(\lambda,\phi,z)\in\Omega_{m}$. Then the pushforward $\Phi_{m,*}$ maps functions in the reference element $I^{3}$ to functions in the physical element $\Omega_{m}$, see for example [26, 27]. For this reason it suffices to explore the analysis on the reference domain $I^{3}$. Additionally, the multi-element case follows the standard approach in finite elements.

Remark 2

If a differential geometry formulation was used, the physical quantities would be represented by differential k-forms and the map $\Phi_{m}:(\xi,\eta,\zeta)\in I^{3}\mapsto(\lambda,\phi,z)\in\Omega_{m}$ would generate a pullback, $\Phi_{m}^{*}$, mapping $k$-forms in physical space, $\Omega_{m}$, to $k$-forms in the reference element, $I^{3}$, [25].

The discrete weak formulation can be stated as: Given $\Omega=I^{3}$, the polynomial degree $N$ and a Coriolis term $\bm{f}_{h}\in\mathcal{W}_{h,\perp}(\Omega)$, for any time $t\in(0,t_{F}]$ find $\bm{u}_{h},\bm{U}_{h},\bm{F}_{h},\bm{P}_{h}\in\mathcal{U}_{h}(\Omega)$, $\theta_{h}\in\mathcal{U}_{h,\perp}(\Omega)$, $\bm{\omega}_{h}\in\mathcal{W}_{h}(\Omega)$, and $\Pi_{h},\rho_{h},\Theta_{h}\in\mathcal{Q}_{h}(\Omega)$ such that

[TABLE]

and

[TABLE]

Using the expansions for all unknowns, (48), (49) may be written as: Find $\bm{\mathsf{u}},\bm{\mathsf{U}},\bm{\mathsf{F}},\bm{\mathsf{P}}\in\mathbb{R}^{d_{\mathcal{U}}}$, $\mathsf{\theta}\in\mathbb{R}^{d_{\mathcal{U}_{\perp}}}$, $\bm{\mathsf{\omega}}\in\mathbb{R}^{d_{\mathcal{W}}}$, and $\mathsf{\Pi},\mathsf{\rho},\mathsf{\Theta}\in\mathbb{R}^{d_{\mathcal{Q}}}$ such that

[TABLE]

and

[TABLE]

4.3 Split discretization

It is important to note that functions in $\mathcal{W}_{\perp,h}$ and $\mathcal{U}_{\parallel,h}$ are discontinuous across vertical element boundaries, while those in $\mathcal{U}_{\perp,h}$ are discontinuous across horizontal boundaries. Similarly functions in $\mathcal{Q}_{h}$ are discontinuous across both vertical and horizontal boundaries. These properties allow us to split the three dimensional problem presented in (50), (51) into separate horizontal and vertical problems, and in doing so avoid solving any global three dimensional implicit systems.

The split discretization is obtained by using (36), (37), (38), and (39), together with (41) and (42) in (50), (51) in order to obtain a discrete version of the split Euler equations (40). The horizontal discrete equations are given by

[TABLE]

where we have introduced $\omega_{h,\parallel,\parallel}:=\sum_{i=0}^{d_{\mathcal{W}_{\parallel}}-1}\mathsf{\omega}_{i,\parallel,\parallel}\bm{\epsilon}^{\mathcal{W}_{\parallel}}_{i}$, $\omega_{h,\parallel,\perp}:=\sum_{i=0}^{d_{\mathcal{W}_{\parallel}}-1}\mathsf{\omega}_{i,\parallel,\perp}\bm{\epsilon}^{\mathcal{W}_{\parallel}}_{i}$. In the same way, the vertical discrete equations are

[TABLE]

where we have introduced $\omega_{h,\perp}:=\sum_{i=0}^{d_{\mathcal{W}_{\perp}}-1}\mathsf{\omega}_{i,\perp}\bm{\epsilon}^{\mathcal{W}_{\perp}}_{i}$. Note that for efficiency we do not assemble the second term in (53a) for simulations at hydrostatic resolutions (and as such we do not apply the diagnostic equation (53e) for $\omega_{h,\perp}$). This term represents the horizontal advection of vertical velocity, and as such is not significant for hydrostatic motions where the vertical scales are small with respect to the horizontal scales, and so the horizontal gradients of vertical terms have minimal impact on the dynamics. While this term is also omitted for hydrostatic formulations [22], our motivation here is based on simple scale analysis. Additionally the rotational terms have no projection onto the energy of the system within the mimetic discretization [3] (though they do re-arrange kinetic energy). At non-hydrostatic resolutions however we do include this term, since the horizontal and vertical scales are closer to unity, and failure to incorporate this term means that vertical motions are not being properly transported with the horizontal flow, such that upstream vertical oscillations may be exaggerated.

The only equations in the above system that cannot be effectively split between the horizontal and vertical systems are (52e) and (53e), since these involve vertical gradients of a horizontally continuous field and vice versa. These terms are required for the vertical advection of horizontal velocity and the horizontal advection of vertical velocity respectively. However by employing a horizontal velocity space which is piecewise constant in the vertical, and a vertical velocity that is piecewise linear in the vertical, we can avoid the need to diagnose these terms through the use of global matrices by employing a direct differencing at the vertical layer interfaces.

Additionally, we also have the flux form equations for density and density weighted potential temperature transport that contain both vertical and horizontal components. While we have not included these in the split systems described in (52) and (53), since doing so incurs a temporal splitting error, in practice these equations are also split between their horizontal and vertical components. For the Strang carryover scheme detailed in Section 5.3 this results in a second order temporal error for the full system. These equations are given as

[TABLE]

Note that the diagnostic equations for potential temperature (51e) and the equation of state (51f) are also included in both the horizontal and vertical systems. We also note that the mass matrices all cancel in (54), such that these flux form transport equations may effectively be solved in the strong form [3, 5] for the point-wise conservation of mass and mass weighted potential temperature.

Remark 3

In (52), (53), and (54), we have used the fact that the incidence matrices can be written as

[TABLE]

where $\bm{\mathsf{E}}^{1,0}_{\parallel}$ is a $d_{\mathcal{W}_{\parallel}}\times d_{\mathcal{P}}$ matrix, $\bm{\mathsf{E}}^{1,0}_{\perp}$ is a $d_{\mathcal{W}_{\perp}}\times d_{\mathcal{P}}$ matrix, $\bm{\mathsf{E}}^{2,1}_{\parallel,\parallel}$ is a $d_{\mathcal{U}_{\parallel}}\times d_{\mathcal{W}_{\parallel}}$ matrix, $\bm{\mathsf{E}}^{2,1}_{\parallel,\perp}$ is a $d_{\mathcal{U}_{\parallel}}\times d_{\mathcal{W}_{\perp}}$ matrix, $\bm{\mathsf{E}}^{2,1}_{\perp}$ is a $d_{\mathcal{U}_{\perp}}\times d_{\mathcal{W}_{\parallel}}$ matrix, $\bm{\mathsf{E}}^{3,2}_{\parallel}$ is a $d_{\mathcal{Q}}\times d_{\mathcal{U}_{\parallel}}$ matrix, and $\bm{\mathsf{E}}^{3,2}_{\perp}$ is a $d_{\mathcal{Q}}\times d_{\mathcal{U}_{\perp}}$ matrix.

Remark 4

In (52), (53), and (54), two important points to note are

[TABLE]

and

[TABLE]

In compact matrix notation we can write (52) as

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

In a similar manner, the vertical equations, (53), can be written in compact matrix notation as

[TABLE]

with

[TABLE]

[TABLE]

In (57a) the vector $\bm{\mathsf{z}}$ represents the discrete vertical coordinate projected onto $\mathcal{Q}_{h}$. While it may seem counter-intuitive to take a discrete vertical gradient of the vertical coordinate itself, rather than just representing this as unity, this weak representation is required in order to preserve the exact balance of kinetic and potential energy exchanges, as shown below in Section 4.4.

Finally, (54) may be written in compact matrix notation as

[TABLE]

with

[TABLE]

4.4 Discrete energetics

The conservation of energy for the rotating shallow water equations via balanced kinetic-potential exchanges has previously been analysed [3] and experimentally verified [5] for the mixed mimetic spectral element method. In terms of energetics, the qualitative difference between the rotating shallow water equations and the compressible Euler equations is the presence of kinetic and internal energy exchanges. As such we here extend the previous analysis to demonstrate that these exchanges may be balanced in the discrete from.

As seen before, we have that the discrete velocity field, $\bm{u}_{h}$, density, $\rho_{h}$, and density weighted potential temperature, $\Theta_{h}$, are

[TABLE]

The discrete Hamiltonian $\mathcal{H}_{h}:=\mathcal{H}[\bm{u}_{h},\rho_{h},\Theta_{h}]$ is then given by

[TABLE]

Using the definition of the variational derivative, see for example [28], we can compute the variational derivative of the Hamiltonian with respect to the velocity, $\frac{\delta\mathcal{H}}{\delta\bm{u}_{h}}$

[TABLE]

The left hand side of this expression may be evaluated to yield

[TABLE]

Since this expression is valid for all $\bm{v}_{h}\in\mathcal{U}_{h}$, then

[TABLE]

Following the same procedure we may obtain the weak equations for the variational derivative with respect to $\rho_{h}$

[TABLE]

and the variational derivative with respect to $\Theta_{h}$

[TABLE]

The discrete compressible Euler equations (48), (49) may then be formulated as a skew-symmetric system for the discrete analogue of (16)-(19) as

[TABLE]

Multiplying both sides by $\begin{bmatrix}\bm{\mathsf{U}}^{\top}&\mathsf{\Phi}^{\top}&\mathsf{\Pi}^{\top}\end{bmatrix}$, gives

[TABLE]

where

[TABLE]

Note that for $\mathsf{R}_{q,ij}:=\left\langle\bm{\epsilon}^{\mathcal{U}}_{i}\,,\,\bm{q}_{h}\times\bm{\epsilon}^{\mathcal{U}}_{j}\right\rangle_{\Omega}$, where $\bm{q}_{h}$ is the potential vorticity [3], this is itself a skew-symmetric operator such that $\bm{\mathsf{U}}^{\top}\bm{\mathsf{R}}_{q}\bm{\mathsf{U}}=\bm{\mathsf{U}}^{\top}\bm{\mathsf{R}}\bm{\mathsf{u}}=\bm{\mathsf{0}}$. As such neither $\bm{\mathsf{R}}_{q}$ nor $\bm{\mathsf{R}}$ projects onto the energy in the discrete form.

Note also that the pressure gradient diagnostic equation (51a) and the temperature flux diagnostic equation (51c) appear within the skew-symmetric operator in (66) within the top right and bottom left blocks respectively. The discrete energy exchanges are therefore given as

[TABLE]

The right hand side terms of (69) exactly balance those of (70) and (71), thus allowing for the exact balances of kinetic to potential and kinetic to internal energy respectively. This holds for both the horizontal and vertical discretisations presented above, assuming periodic boundary conditions in the horizontal, homogeneous Dirichlet boundary conditions for the vertical velocity (5) and Neumann conditions for the Exner pressure (6) in the vertical. As an aside we note that energetic consistency is satisfied independent of the choice of function space for $\theta_{h}$, since this only appears within $\bm{\mathsf{S}}^{\mathcal{U}}$, thus justifying our choice to represent $\theta_{h}\in\mathcal{U}_{\perp}$.

4.5 Metric terms

The Jacobian matrix between local element coordinates $\bm{\xi}:=(\xi,\eta,\zeta)$ and global coordinates $\bm{x}:=(\lambda,\phi,z)$ is given as

[TABLE]

where the subscripts represent derivatives with respect to local element coordinates and we have assumed that all horizontal layers are perfectly flat, such that the projection of vertical local coordinates onto horizontal global coordinates and horizontal local coordinates onto vertical global coordinates are zero. The $H(\mathrm{curl},\Omega)$, $H(\mathrm{div},\Omega)$ and $L^{2}(\Omega)$ forms of the Piola transformation are given respectively as [11, 12]

[TABLE]

where $J$ is the determinant of the Jacobian matrix. The metric transformations for the respective mass matrices are therefore $(\mathsf{J}^{-\top})^{\top}\mathsf{J}^{-\top}$, $\frac{1}{J^{2}}\mathsf{J}^{\top}\mathsf{J}$ and $\frac{1}{J^{2}}$. Since the horizontal and vertical components of both the $H(\mathrm{curl},\Omega)$ and $H(\mathrm{div},\Omega)$ transformations are orthogonal, these metric transformations are further simplified for these spaces as

[TABLE]

As a special case, we also note that an inner product between basis functions in $\mathcal{U}_{h,\perp}$ and functions in $\mathcal{Q}_{h}$ has the metric term

[TABLE]

As stated these metric terms do not account for any projection of horizontal vector components onto vertical components and vice versa. As such we limit ourselves here to the case where the horizontal and vertical degrees of freedom are strictly orthogonal, and we are unable to account for the tilting of layers or bottom topography in the current implementation. We note however that the tilting of horizontal layers in order to represent topography is naturally incorporated into the finite element formulation via these cross terms, and these may be implemented at a later date with minimal disruption to the current formulation.

Note that the vorticity term $\langle\bm{\sigma}_{h},\bm{\omega}_{h,\parallel,\parallel}\times u_{h,\perp}\rangle$ is alternatively formulated as $\langle\bm{\sigma}_{h},u_{h,\perp}\nabla_{\perp}\bm{u}_{h,\parallel}\rangle$ (with the vertical derivative derived in the weak form). In this form the vertical velocity derivatives may be interpreted as being oriented normal to the edges, rather than tangent to them, and as such the $H(\mathrm{div},\Omega)$ form of the Piola transformation is used to construct the metric term for this operator as well as for (52e).

5 Time stepping

In this section we introduce the horizontally-explicit/vertically implicit (HEVI) time stepping scheme employed in the model. Such schemes are popular in non-hydrostatic atmospheric modelling since they negate the explicit time step restrictions associated with the fast time scales of the vertical gravity and acoustic waves [29, 13, 30, 15, 31]. These schemes are often constructed by perturbing the vertical dynamics around a leading order state of hydrostatic balance [13, 15] and/or further splitting the vertical dynamics so as to only solve for the terms responsible for the fast dynamics implicitly [31]. Here we avoid these formulations due to our concern that they may break the skew-symmetric structure of the discrete system (66) that is central to the conservation of energy in the spatial discretisation. Nevertheless we note that while our current implicit scheme detailed below preserves the exact balance kinetic to potential energy exchanges, balance is not satisfied for the vertical kinetic to internal energy exchanges due to the formulation of the pressure gradient term. Restoring this balance within the implicit vertical scheme is a subject of ongoing research. We further emphasise that both our time splitting scheme as well as our individual implicit vertical and explicit horizontal time integration schemes will conserve energy only to truncation order in time.

5.1 Directional splitting

We use a Strang carryover splitting scheme to partition the horizontal and vertical dynamics [29, 30]. Consider the previously introduced splitting of the spatial operator (here denoted by $\bm{L}$) into a vertical ($\bm{L}^{v}$) and horizontal ($\bm{L}^{h}$) components: equations (56) (horizontal), and equations (57) together with equations (58) (vertical). With this splitting we can write the full system of equations as

[TABLE]

where $\bm{b}:=[\bm{u},\bm{U},\bm{F},\bm{P},\bm{\omega},\rho,\theta,\Theta,\Pi]$. Then, the time stepping procedure to evolve from time step $t^{n}$ to $t^{n+1}$ follows the sequence

[TABLE]

This splitted time evolution procedure is then numerically integrated. The first vertical half time step (77) is explicitly integrated using a forward Euler scheme (carried over) as

[TABLE]

The horizontal full time step (78) is numerically integrated using an explicit third order stiffly stable Runge-Kutta scheme, see Section 5.2 for more details. Finally, the second vertical half step (79) is integrated using an implicit Euler scheme. A nonlinear Picard iteration, $m$, is applied to solve this vertical system, which is assumed to converge once $|\bm{b}^{m+1}(t^{n+1})-\bm{b}^{m}(t^{n+1})|_{2}/|\bm{b}^{m+1}(t^{n+1})|_{2}<10.0^{-8}$, where $|\cdot|_{2}$ is the $L^{2}$ norm. For a detailed discussion see Section 5.3. For the 30 models levels of the baroclinic test case described below in Section 6.1, 13 Picard iterations are required to reach convergence.

This scheme is formally second order accurate in time due to the Strang splitting, with the horiztonal scheme being third order accurate, and the vertical scheme being equivalent to a second order trapezoidal scheme [30]. Our anecdotal experience is that if a first order horizontal-vertical splitting is used, as opposed to the second order splitting described here, then the horizontal-vertical coupling is too weak to allow for the correct transfer of potential and internal energy to vertical kinetic energy required to properly simulate the baroclinic instability described below in Section 6.1.

Since all the solution variables in $\bm{b}$ required to solve for the vertical dynamics are discontinuous across horizontal element boundaries, the inner linear system above is solved in serial for each horizontal element, via a direct LU solve using the PETSc library [33, 34, 35]. Similarly, since all the solution variables involved in the vertical dynamics are discretized on function spaces that are discontinuous across vertical element boundaries, this system may be solved independently, in parallel for each horizontal layer. The mass matrices are solved using the PETSc GMRES solver, with a block Jacobi preconditioner for each element, as was done for the shallow water equations in our previous work [5].

5.2 Explicit horizontal solve

The horizontal dynamics are solved using an explicit stiffly stable Runge-Kutta scheme of the form [32, 29]

[TABLE]

where $b^{\prime}(t)$ and $b^{\prime\prime}(t)$ are as defined in (77) and (78), respectively.

5.3 Implicit vertical solve

The implicit vertical solve involves a half step, $\textstyle\frac{1}{2}\Delta t$, from the state at the end of the horizontal solve, $b^{\prime\prime}(t^{n}+\Delta t)$, to the end of the time level $b^{n+1}$. Following [9] we begin by taking the logarithm of the equation of state (3) as

[TABLE]

Differentiating both sides with respect to time over a half step via the chain rule, while recalling that the second and third terms on the right hand side are simply constants then gives

[TABLE]

Substituting in the potential temperature equation (1c), we then express the evolution of the Exner pressure as

[TABLE]

Note the similarity between (87) and the internal energy evolution equation (14). The vertical dynamics may then be discretized in time (with the vorticity components omitted, as discussed in Section 4.3) as

[TABLE]

We define the additional matrix operators $\bm{\mathsf{M}}^{\mathcal{Q}}_{\Theta}$, $\bm{\mathsf{M}}^{\mathcal{Q}}_{\Pi}$ and $\bm{\mathsf{M}}^{\mathcal{U}_{\perp}}_{\Theta}$, for which

[TABLE]

Dropping the $n+1$ superscripts for variables at the end of the vertical half step (but keeping those for variables at time level ′′), the discrete form of (87) is then given as:

[TABLE]

and the discrete form of (88) as:

[TABLE]

Substituting (92) into (93) gives

[TABLE]

Finally, rearranging gives an expression for the vertical velocity at each fixed point Picard iteration as:

[TABLE]

In order to incorporate the pressure gradient term into the left hand side in (95) so as to ensure the stable implicit solution of vertical motions, we have sacrificed the energetically consistent formulation of the vertical pressure gradient term as presented in (53a), (53b) and (66). As such we do not expect the vertical kinetic and internal exchanges to exactly balance, as we do for the horizontal explicit discretization. However we note that unlike the perturbation and operator splitting approaches to HEVI discretisations discussed above, our current approach leaves open the possibility of recovering the exact balance of kinetic and internal energy exchanges. We are actively exploring preconditioning strategies to address this problem. We also note that our implicit vertical solve may potentially dampen the fast vertical motions for which the time scales are not resolved, in comparison to sub-cycled vertically explicit formulations.

The Exner pressure at the current Picard iteration is then derived from (92), and the corresponding Picard iteration solves for the other variables are then given (omitting the horizontal terms) as

[TABLE]

5.4 Dissipative terms

In order to stabilise the model we also include various dissipative terms. These include a biharmonic viscosity on both the horizontal momentum and temperature equations [36, 5] with a value of $0.072\Delta x^{3.2}$, where $\Delta x$ is the average spacing between GLL nodes. A Rayleigh friction term with a coefficient of $0.2$ is also applied to the top layer of the vertical momentum equation, as is often used in atmospheric models to suppress orographically forced gravity waves [37]. This term is added to (53a) as $\sum_{i=0}^{d_{\mathcal{U}_{\perp}}-1}0.2\left\langle\bm{\epsilon}^{\mathcal{U}_{\perp}}_{j}\delta_{j,l}\,,\,\delta_{i,k}\bm{\epsilon}^{\mathcal{U}_{\perp}}_{i}\right\rangle_{\Omega}\mathsf{u}_{i,\perp}$, where $k,l$ are degrees of freedom for the trial and test functions in the top layer only, and $\delta_{i,k}$ is the standard delta function. While this Rayleigh friction term is not strictly necessary for the stability of the simulation presented here, it greatly reduces the noise in the energetic profiles as the model adjusts to a state of hydrostatic balance from its initial conditions. No viscosity is required to stabilize the vertical solution, perhaps because this is only second order accurate, so the internal dissipation of this low order discretization is sufficient to prevent nonlinear instabilities.

We emphasise that these dissipative terms will necessarily remove energy from the system, and in the case of the horizontal viscous terms are necessary to stabilise the model by arresting the nonlinear cascades at the grid scale. While upwinding terms may be used instead of viscosity as a means of suppressing grid scale oscillations without dissipating energy, by systematically adding these to the skew-symmetric formulation [38], here we limit ourselves to a consideration only of the internal energetic processes, and not the external forcings.

6 Results

6.1 Baroclinic instability

We validate the model using a dry baroclinic instability test case [39] with the shallow atmosphere approximation. The appeal of this test case is that the initial condition is specified for a $z-$level vertical coordinate, whereas other such test cases that are defined on pressure level vertical coordinates require the solution of a nonlinear problem in order to compute the corresponding $z-$level configuration. The initial state is one of geostrophic horizontal and hydrostatic vertical balance, overlaid with a small, $\mathcal{O}(1\mathrm{m/s})$, perturbation to the zonal and meridional velocity components.

The model was run with $24\times 24$ elements of degree $p=3$ on each face of the cubed sphere (and linear elements in the vertical), for an averaged resolution of $\Delta x\approx 128\mathrm{km}$ and 30 vertical levels on 96 processors ($4^{2}=16$ per face of the cubed sphere) with a time step of $\Delta t=120\mathrm{s}$. While the vertical dynamics are all solved implicitly and so do not limit the time step size, the explicit horizontal dynamics present both diffusive and advective CFL restrictions due to the biharmonic viscosity and sound waves respectively.

Figures 2 and 3 show the zonal averages of density $\rho$, Exner pressure $\Pi$, potential temperature $\theta$ and zonal velocity $u$ at day 10 (solid lines), as well as the differences between the final and initial states. These profiles show little difference between the initial and final states, with the exception of the zonal velocity, which exhibits a small kink near the bottom boundary where the baroclinic instability occurs, demonstrating that the leading order geostrophic and hydrostatic balances in the horizontal and vertical are well satisfied. The potential temperature (Fig. 3), is approximately $15^{\circ}\mathrm{K}$ cooler in the top layer at day 10, which is outside the range of the color bar. This is due to the hydrostatic adjustment of the mean state, as discussed below.

Figures 4 and 5 show the evolution of the kinetic (horizontal and vertical), potential and internal energy with time, and the associated exchanges. These are shown on both logarithmic scales for their normalised absolute difference from initial value (Fig. 4), and as a direct difference between their current and initial value (Fig. 5).

The growth in the baroclinic instability is evident in the increase in kinetic energy, and the reduction in both potential and internal energy as isopycnals flatten in the region of the instability. Note that the total amounts of potential and internal energy are approximately $3.6\times 10^{23}$ and $9.2\times 10^{23}\mathrm{kg\cdot m^{2}s^{-2}}$ respectively, and so are several orders of magnitude greater than the amounts of horizontal and vertical kinetic energy (approximately $4.0\times 10^{20}$ and $2.5\times 10^{13}\mathrm{kg\cdot m^{2}s^{-2}}$ respectively). As such the flattening of the density contours from which the baroclinic instability draws energy are barely evident in Fig. 2.

Figure 5 also shows the normalised sum of the globally integrated kinetic to potential (69) and potential to kinetic (70) energy exchanges, as well as the (horizontal) globally integrated kinetic to internal (69) and internal to kinetic (71) exchanges. The kinetic to potential and potential to kinetic exchanges balance to machine precision, owing to the exact skew-symmetry of these operators, as shown in (66). The kinetic to internal exchanges balance only to approximately $\mathcal{O}(10^{-6})$, and are sometimes out by a factor of $\mathcal{O}(10^{-3})$. This is due to the fact that both terms involve the inverse of a $H(\mathrm{div},\Omega)$ mass matrix, the action of which is approximated via an iterative Krylov method. We therefore anticipate pointwise errors in these exchanges, which may be compounded by the fact that these inverses are approximated twice, once to determine the Exner pressure gradient (51a), and a second time to determine the density weighted potential temperature flux (51c). Nevertheless, these errors exhibit no systematic drift throughout the simulation. We do not show the vertical kinetic to internal balance errors, as our vertical implicit solve does not satisfy this balance.

The observation that potential energy is greater than its initial value for most of the simulation, while internal energy is smaller, as shown in Figs. 4 and 5, is explained by the fact that the initial condition is not precisely in hydrostatic balance. As such the initial adjustment process leads to a slight rise in the fluid, resulting in an increase in potential energy, and a corresponding reduction in pressure, leading to a reduction of internal energy via the equation of state. These changes are of $\mathcal{O}(10^{-4})$ compared to the total amounts of potential and internal energy in the system. The high frequency oscillation in the internal and potential energies observed over the first 24 hours of the simulation during this adjustment process is reduced by the application of the Rayleigh friction to the top layer of the vertical momentum equation.

We present the bottom level pressure $p$ and temperature, $T$, and the vertical component of the relative vorticity, $\omega$, at $z\approx 1.5\mathrm{km}$ at days $8$ and $10$ in Figs. 6, 7 and 8, as well as a meridional cross section of the pressure perturbation at $50^{\circ}\mathrm{N}$ in Fig. 9. In the cases of the pressure and temperature, these are reconstructed from the model variables as $p=p_{0}(\Pi/c_{p})^{c_{p}/R}$ and $T=\theta\Pi/c_{p}$. The pressure perturbation in Fig. 9 is derived by removing the average pressure at the corresponding vertical level at $50^{\circ}\mathrm{S}$. These results compare well with the previously published test case results [39] in both shape and magnitude, and clearly show the signal of the baroclinic instability.

For completeness we also show the results at days 8 and 10 for the original model variables of bottom level Exner pressure and potential temperature and relative vertical vorticity at $z\approx 1.5\mathrm{km}$ in Figs. 10, 11 and 12 respectively. These are presented for the northern hemisphere only, looking down from the north pole. These results are perhaps slightly sharper than the reconstructed temperature and pressure fields since they are interpolated directly from the degrees of freedom.

6.2 Non-hydrostatic gravity wave

The baroclinic instability test case detailed above provides an excellent means of validating both the leading order balance relations and the horizontal and vertical coupling required to correctly simulate a secondary nonlinear three-dimensional instability. However the limitation of this test case is that it is configured for spatial resolutions at which the non-hydrostatic dynamics are negligible. As such we also test the behaviour of the model for the propagation of a non-hydrostatic gravity wave driven by a potential temperature perturbation on a planet with a reduced radius $125$ times smaller than that of the earth. This test was originally proposed as part of the 2012 DCMIP workshop, for which numerous hydrostatic and non-hydrostatic dynamical cores presented results [40]. Specific details of the initial configuration can be found within the DCMIP test case document on the web site.

As for the baroclinic instability test, we run the simulation with a resolution of $24\times 24$ elements of degree $p=3$ in each cubed sphere panel. We use 16 evenly spaced vertical levels over a total height of $10,000.0$m, and a time step of $\Delta t=0.5$s for a total simulation time of $3600$s. No Rayleigh damping or viscosity is applied in the vertical, however we rescale the horizontal biharmonic viscosity by a factor of $2.0$ for both the momentum and temperature equations for a value of $0.144\Delta x^{3.2}$. Since the horizontal and vertical scales are of equal order in this configuration, we have also included the additional vorticity term in (53a) and the associated diagnostic term (53e), both of which were omitted from the baroclinic test case, where the vertical scales were $\mathcal{O}(10^{-3})$ with respect to the horizontal scales.

This test case is especially challenging for our higher order spectral element formulation, firstly since we have not applied any sort of upwinding or monotonicity preservation method to either the continuity or temperature equation, and secondly because unlike other models we formulate our energy equation as a flux form equation for the density weighted potential temperature, $\Theta=\rho\theta$, from which the potential temperature, $\theta$, is then diagnosed. This is in contrast to a material form of the potential temperature advection, in which upwinding is directly applied in the flux reconstruction [12]. For this reason we configure our model with slightly higher resolution than that specified [40].

Figure 13 shows the longitude-height equatorial ($\phi=0^{\circ}$) cross section of the potential temperature perturbation, $\theta^{\prime}(\lambda,0,z)=\theta(\lambda,0,z)-\bar{\theta}(z)$, where $\bar{\theta}(z)$ is the mean potential temperature at a given height, after 30 minutes and 1 hour. While the structure and evolution of the perturbation qualitatively match the results presented for other non-hydrostatic models, the absence of a monotone upwinded advection scheme means that our results are slightly more oscillatory than those of other models. Moreover our results also show the ejection of a small hot bubble which rises to the top of the domain, which is also most likely an artefact of our non-monotone scheme. In future work we will investigate energetically consistent methods for stabilising the advective terms, as has previously been addressed for the shallow water equations [38], and hence recover less oscillatory solutions at non-hydrostatic scales. Fortunately the energetically consistent flux form of the potential temperature equation presented here (48c, 49c), matches with the formulation necessary to ensure that tracer concentration fluxes are consistent with respect to the mass flux in monotone schemes [41].

We also use this test to study the conservation properties of the model, as shown in Fig. 14. The third order Runge-Kutta scheme used for horizontal advection, as given in (82-84) involves a linear extrapolation from previous to current states. It would appear that this convex operation leads to a small loss of mass that grows linearly with time. If we replace the third order scheme with a second order scheme of the form $\bm{b}^{(1)}=\bm{b}^{\prime}+\Delta tL^{h}(\bm{b}^{n});\bm{b}^{\prime\prime}=\bm{b}^{\prime}+\Delta t(L^{h}(\bm{b}^{n})+L^{h}(\bm{b}^{(1)}))/2$, which does not involve a linear extrapolation, then exact energy conservation is recovered.

Figure 14 also shows the power (in units of $\mathrm{kg\cdot m^{2}s^{-3}}$) associated with the horizontal and vertical energetic exchanges, and the sum of vertical exchanges. Since the vertical gravity wave oscillates at a high frequency compared to the horizontal dynamics, we only show these results for the first 12 minutes of the simulation. The gravity wave is expressed through the periodic exchange of potential to internal energy, via vertical motions. The power associated with these exchanges is of $\mathcal{O}(10^{13}\mathrm{kg\cdot m^{2}s^{-3}})$, which is approximately $\mathcal{O}(10^{4})$ times as large as both the sum of these exchanges, as well as the horizontal kinetic to internal exchanges, which are observed to occur on a longer time scale than the gravity wave.

7 Conclusions

A model of the rotating compressible Euler equations on the cubed sphere using the mixed mimetic spectral element method is presented. The model uses a Strang carryover directional splitting scheme with an implicit Picard solver for the vertical dynamics in order to negate the CFL condition of the vertical acoustic and gravity waves. The discontinuities in the discrete function spaces are exploited so as to solve for each horizontal layer and each vertical element independently in order to avoid the need to invert a global 3D mass matrix. The forcing and flux terms are constructed so as to take advantage of the adjoint relations between the discrete gradient and divergence operators and balance the exchanges of kinetic, potential and internal energy. The exception to this is the construction of the vertical pressure gradient operator, which violates this principal in order to allow for a stable implicit solve of the vertical dynamics, with the consequence that the vertical kinetic to internal energy exchanges do not exactly balance.

In future work we will explore fully implicit formulations so as to avoid the horizontal vertical splitting and restore the balance of vertical kinetic internal energy exchanges to the discrete system. We will also explore energetically consistent flux reconstructions so as to recover less oscillatory solutions at non-hydrostatic scales.

We also intend to optimise the parallel formulation of the code as part of additional future work. For the baroclinic instability simulation presented here at a resolution of $128\mathrm{km}$ on 30 vertical levels over 96 processors on a set of Dell PowerEdge M630 servers, the model only runs approximately 4.5 times faster than real time. By re-formulating the parallel decomposition of the different function spaces on the cubed sphere we hope to better align these decompositions with the native PETSc decompositions of vectors and matrices so as to reduce parallel communication.

8 Acknowledgments

David Lee would like to thank Drs. Marcus Thatcher and John McGregor for their continued support and access to computing resources, and Dr. Thomas Melvin for several helpful discussions. We are also grateful to the three anonymous reviewers, each of whom contributed thoughtful and insightful comments that helped improve the quality of this article.

Appendix A Incident matrix

A.1 Discrete gradient (2D)

The continuous equation

[TABLE]

can be represented on the mesh depicted on Figure 15 by the following algebraic system of equations

[TABLE]

where the incidence matrix $\bm{\mathsf{E}}^{1,0}$ is given by

[TABLE]

We note that we can easily split into a vertical component ($z$-direction) and a parallel component ($x$-direction). Therefore we can rewrite (100) as

[TABLE]

In the same way, we can rewrite the discrete counterpart of this equation, (101), as

[TABLE]

where

[TABLE]

and

[TABLE]

A.2 Discrete curl (2D)

The continuous equation

[TABLE]

can be represented on the mesh depicted on Figure 16 by the following system of equations

[TABLE]

where the incidence matrix $\bm{\mathsf{E}}^{2,1}$ is given by

[TABLE]

We note that, as done before, it is possible to split the curl operator into a vertical component ($z$-direction) and a parallel component ($x$-direction). Therefore we can rewrite (107) as

[TABLE]

Similarly, we can rewrite the discrete counterpart of this equation, (108), as

[TABLE]

Note that since the mesh under consideration, Figure 16, is a two dimensional mesh, $\bm{\omega}$ will only have a component in the $y$-direction. For this reason, $\bm{\mathsf{E}}^{2,1}_{\perp}=\bm{\mathsf{0}}$, and for compactness we suppress this term (for three dimensional meshes this term must be included)

[TABLE]

where

[TABLE]

and

[TABLE]

A.3 Discrete divergence (2D)

For the divergence operator, the continuous equation

[TABLE]

can be expressed on the mesh depicted on Figure 17 by the following system of equations

[TABLE]

where the incidence matrix is given by

[TABLE]

As discussed previously, it is possible to split the divergence operator into a vertical component ($z$-direction) and a parallel component ($x$-direction). Therefore we can rewrite (115) as

[TABLE]

In the same fashion, we can rewrite the discrete version of this equation, (116), as

[TABLE]

where

[TABLE]

and

[TABLE]