# Exact spatial and temporal balance of energy exchanges within a horizontally
explicit/vertically implicit non-hydrostatic atmosphere

David Lee<EMAIL_ADDRESS>Artur Palha  Bureau of Meteorology,
Melbourne, Australia  Faculty of Aerospace Engineering, Delft University of
Technology, Kluyverweg 2, 2629 HS Delft, The Netherlands

###### Abstract

A new horizontally explicit/vertically implicit (HEVI) time splitting scheme
for atmospheric modelling is introduced, for which the horizontal divergence
terms are applied within the implicit vertical substep. The new HEVI scheme is
implemented in conjunction with a mixed mimetic spectral element spatial
discretisation and semi-implicit vertical time stepping scheme that both
preserve the skew-symmetric structure of the non-canonical Hamiltonian form of
the equations of motion. Within this context the new HEVI scheme allows for
the exact balance of all energetic exchanges in space and time. However since
the choice of horizontal fluxes for which this balance is satisfied is not
consistent with the horizontal velocity at the end of the time level the
scheme still admits a temporal energy conservation error. Linearised
eigenvalue analysis shows that similar to a fully implicit method, the new
HEVI scheme is neutrally stable for all buoyancy modes, and unlike a second
order trapezoidal HEVI scheme is stable for all acoustic modes below a certain
horizontal CFL number. The scheme is validated against standard test cases for
both planetary and non-hydrostatic regimes. For the planetary scale baroclinic
instability test case, the new formulation exhibits a secondary oscillation in
the potential to kinetic energy power exchanges, with a temporal frequency
approximately four times that exhibited by a horizontally third order,
vertically second order trapezoidal scheme. For the non-hydrostatic test case,
the vertical upwinding of the potential temperature diagnostic equation is
shown to reduce spurious oscillations without altering the energetics of the
solution, since this upwinding is performed in an energetically consistent
manner. For this test case, which is configured on an affine geometry, the
exact balance of energy exchanges allows the model to run stably without any
form of dissipation.

Horizontally explicit/vertically implicit (HEVI) schemes are a popular
approach to time stepping in atmospheric models. This is on account of the
importance of removing the CFL limit imposed by the explicit representation of
vertical dynamics, and the desire to avoid the computational expense of a
three dimensional implicit solve at each nonlinear iteration. Numerous HEVI
schemes have been implemented [1, 2, 3, 4, 5] in previous models. These
schemes are typically motived by a desire to improve upon the stability and
dispersion properties of the second order trapezoidal HEVI scheme through the
use of additional sub steps within each time step [6, 7]. The application of
additional substeps qualitatively improves the stepping scheme by ensuring the
linear stability of the horizontal buoyancy and acoustic models below some CFL
limit.

In the present article a new HEVI scheme is introduced, motivated by a desire
to preserve the exact balance of energetic exchanges. The scheme is
implemented within the context of a mixed mimetic spectral element spatial
discretisation [8] and energy conserving implicit vertical time stepping
method [9] that together allow for the exact balance of all energy exchanges
in space and time. This is achieved by evaluating the horizontal mass and
temperature fluxes within the vertically implicit sub-step, after first
computing a provisional horizontal velocity from which to derive a second
order temporal representation of the horizontal fluxes. Unlike other HEVI
schemes, the new method is composed of only two horizontal and one vertical
substep. In contrast to the standard second order trapezoidal splitting
scheme, linear analysis shows the new method to be unconditionally stable for
all acoustic modes below some horizontal CFL number, and all buoyancy modes
independent of CFL number. While the energetic exchanges are all exactly
balanced, since the provisional horizontal velocity does not exactly match the
final horizontal velocity at the end of the time step, the scheme admits a
horizontal energy conservation error, as is the case for explicit schemes used
to solve compressible flow problems.

In section 1 the continuous 3D compressible Euler equations for atmospheric
flows are introduced in variational form. The discrete analogue of these
equations and a discussion of the energy conserving time stepping are then
presented in section 2, and the new HEVI splitting scheme is described in
section 3. Section 4 describes the energetically consistent stabilisation of
spurious oscillations in potential temperature via the upwinding of the
corresponding test space. The reproduction of standard test cases is described
in section 5. Finally the conclusions are discussed in section 6.

##  1 Skew-symmetric variational structure of the compressible Euler equations

The compressible Euler equations for atmospheric flows may be expressed in
skew-symmetric form for the evolution of the velocity, $\boldsymbol{u}$,
density, $\rho$, and density weighted potential temperature, $\Theta$, as:

$\displaystyle\frac{\partial\boldsymbol{u}}{\partial t}$
$\displaystyle=-\boldsymbol{q}\times\boldsymbol{U}-\nabla\Phi-\theta\nabla\Pi$
(1a) $\displaystyle\frac{\partial\rho}{\partial t}$
$\displaystyle=-\nabla\cdot\boldsymbol{U}$ (1b)
$\displaystyle\frac{\partial\Theta}{\partial t}$
$\displaystyle=-\nabla\cdot(\theta\boldsymbol{U})$ (1c)

where $\boldsymbol{q}:=(\nabla\times\boldsymbol{u}+\boldsymbol{f})/\rho$ is
the potential vorticity (and $\boldsymbol{f}=f\boldsymbol{e}_{z}$ is the
Coriolis term), $\boldsymbol{U}:=\rho\boldsymbol{u}$ is the mass flux,
$\Phi:=\frac{1}{2}\boldsymbol{u}\cdot\boldsymbol{u}+\rho gz$ is the Bernoulli
function (with $g$ the gravitational acceleration and $z$ the height),
$\theta:=\Theta/\rho$ is the potential temperature, and
$\Pi:=c_{p}\left(R\Theta/p_{0}\right)^{\frac{R}{c_{v}}}$ is the Exner pressure
(where $c_{v}$ and $c_{p}$ are the specific heats at constant volume and
temperature respectively, $R$ is the ideal gas constant and $p_{0}$ is the
reference pressure). The corresponding Hamiltonian (total energy) is given as:

${H}=\int_{\Omega}\overbrace{\frac{1}{2}\rho\boldsymbol{u}\cdot\boldsymbol{u}}^{{K}}+\overbrace{\vphantom{\frac{1}{2}}\rho
gz}^{{P}}+\overbrace{\frac{c_{v}}{c_{p}}\Theta\Pi}^{{I}}\mathrm{d}\Omega,$ (2)

where ${K}$ is the kinetic energy, ${P}$ is the potential energy, and ${I}$ is
the internal energy. In this work the domain, $\Omega$ is configured as either
the volume on the surface of a sphere in physical units, or a high resolution
Cartesian geometry which is periodic in the two horizontal dimensions. The
variational derivatives of the Hamiltonian with respect to the dependent
variables ($\boldsymbol{u}$, $\rho$, $\Theta$) are

$\displaystyle\frac{\delta{H}}{\delta\boldsymbol{u}}$
$\displaystyle=\rho\boldsymbol{u}=:\boldsymbol{U},$ (3a)
$\displaystyle\frac{\delta{H}}{\delta\rho}$
$\displaystyle=\frac{1}{2}\boldsymbol{u}\cdot\boldsymbol{u}+gz=:\Phi,$ (3b)
$\displaystyle\frac{\delta{H}}{\delta\Theta}$
$\displaystyle=c_{p}\Bigg{(}\frac{R\Theta}{p_{0}}\Bigg{)}^{\frac{R}{c_{v}}}=:\Pi,$
(3c)

which give the mass flux, $\boldsymbol{U}$, the Bernoulli function, $\Phi$ and
the Exner pressure, $\Pi$ respectively.

The weak form of (1) is given for $\boldsymbol{u},\boldsymbol{U}\in
H(\mathrm{div},\Omega)$ and $\rho,\Theta,\Phi,\Pi\in L^{2}(\Omega)$ as:

$\displaystyle\Bigg{\langle}\boldsymbol{\beta},\frac{\partial\boldsymbol{u}}{\partial
t}\Bigg{\rangle}$
$\displaystyle=-\langle\boldsymbol{\beta},\boldsymbol{q}\times\boldsymbol{U}\rangle+\langle\nabla\cdot\boldsymbol{\beta},\Phi\rangle+\langle\nabla\cdot\mathbb{P}^{\mathrm{div}}\left[\theta\boldsymbol{\beta}\right],\Pi\rangle,\qquad\forall\boldsymbol{\beta}\in
H(\mathrm{div},\Omega)$ (4a)
$\displaystyle\Bigg{\langle}\gamma,\frac{\partial\rho}{\partial
t}\Bigg{\rangle}$
$\displaystyle=-\langle\gamma,\nabla\cdot\boldsymbol{U}\rangle,\qquad\forall\gamma\in
L^{2}(\Omega)$ (4b)
$\displaystyle\Bigg{\langle}\sigma,\frac{\partial\Theta}{\partial
t}\Bigg{\rangle}$
$\displaystyle=-\langle\sigma,\nabla\cdot\mathbb{P}^{\mathrm{div}}\left[\theta\boldsymbol{U}\right]\rangle,\qquad\forall\sigma\in
L^{2}(\Omega)$ (4c)

where we have invoked inner product notation as $\langle
a,b\rangle=\int_{\Omega}ab\mathrm{d}\Omega$, the integration-by-parts relation
$\langle\boldsymbol{\beta},\nabla\alpha\rangle=-\langle\nabla\cdot\boldsymbol{\beta},\alpha\rangle$
(assuming periodic or homogeneous boundary conditions) in (4), and introduced
the $L^{2}$ projection to $H(\mathrm{div},\Omega)$

$\mathbb{P}^{\mathrm{div}}:\left[L^{2}(\Omega)\right]^{3}\rightarrow
H(\mathrm{div},\Omega)\,.$ (5)

Specifically for $\boldsymbol{v}\in\left[L^{2}(\Omega)\right]^{3}$, the
projection $\mathbb{P}^{\mathrm{div}}\left[\boldsymbol{v}\right]\in
H(\mathrm{div},\Omega)$ is given by

$\langle\boldsymbol{\beta},\mathbb{P}^{\mathrm{div}}\left[\boldsymbol{v}\right]\rangle=\langle\boldsymbol{\beta},\boldsymbol{v}\rangle\,,\qquad\forall\boldsymbol{\beta}\in
H(\mathrm{div},\Omega)\,.$ (6)

Energy conservation is assured for the choice of
$\boldsymbol{\beta}=\boldsymbol{U}$, $\gamma=\Phi$, $\sigma=\Pi$ by summing
all equations in (4), such that

$\Bigg{\langle}\frac{\delta{H}}{\delta\boldsymbol{u}},\frac{\partial\boldsymbol{u}}{\partial
t}\Bigg{\rangle}+\Bigg{\langle}\frac{\delta{H}}{\delta\rho},\frac{\partial\rho}{\partial
t}\Bigg{\rangle}+\Bigg{\langle}\frac{\delta{H}}{\delta\Theta},\frac{\partial\Theta}{\partial
t}\Bigg{\rangle}=\frac{\mathrm{d}H}{\mathrm{d}t}=0.$ (7)

where the key ingredient for energy conservation is the skew-symmetry of (4).

##  2 Discrete formulation

In this section the mixed mimetic spectral element spatial discretisation used
in the model is introduced, as well as an analysis of the geometric properties
of the energy conserving time integration scheme.

###  2.1 Spatial discretisation

For the purposes of this article the salient feature of the mixed mimetic
spectral element method is that it allows for the preservation of the skew-
symmetric structure of the compressible Euler equations, and thus the
conservation of energy and energetic exchanges in the discrete form. For a
more detailed discussion the reader is referred to previous work on the use of
this method for geophysical flow modelling [10, 11, 8], as well as more
foundational works on the subject [12, 13, 14, 15]. In order to begin this
discussion we introduce the finite dimensional subspaces
$\mathcal{P}_{h}\subset H^{1}(\Omega)$, $\mathcal{W}_{h}\subset
H(\mathrm{curl},\Omega)$, $\mathcal{U}_{h}\subset H(\mathrm{div},\Omega)$ and
$\mathcal{Q}\subset L^{2}(\Omega)$. These subspaces are spanned by basis
functions

$\mathrm{span}\\{\psi_{1},\dots,\psi_{\mathrm{d}_{\mathcal{P}}}\\}=\mathcal{P}_{h}\,,\quad\mathrm{span}\\{\boldsymbol{\alpha}_{1},\dots,\boldsymbol{\alpha}_{\mathrm{d}_{\mathcal{W}}}\\}=\mathcal{W}_{h}\,,\quad\mathrm{span}\\{\boldsymbol{\beta}_{1},\dots,\boldsymbol{\beta}_{\mathrm{d}_{\mathcal{U}}}\\}=\mathcal{U}_{h}\,,\quad\text{and}\quad\mathrm{span}\\{\gamma_{1},\dots,\gamma_{\mathrm{d}_{\mathcal{Q}}}\\}=\mathcal{Q}_{h}\,,\quad$
(8)

with $\mathrm{d}_{\mathcal{P}}$, $\mathrm{d}_{\mathcal{W}}$,
$\mathrm{d}_{\mathcal{U}}$, and $\mathrm{d}_{\mathcal{Q}}$, the number of
degrees of freedom of the spaces $\mathcal{P}_{h}$, $\mathcal{W}_{h}$,
$\mathcal{U}_{h}$, and $\mathcal{Q}_{h}$, respectively. The basis functions
themselves are constructed by tensor product combinations of nodal and edge
polynomials [12, 16]. As usual, elements of these spaces can be represented as
a linear combination of the basis functions, for example, for
$p_{h}\in\mathcal{P}_{h}$, $\boldsymbol{w}_{h}\in\mathcal{W}_{h}$,
$\boldsymbol{u}_{h}\in\mathcal{U}_{h}$, and $c_{h}\in\mathcal{Q}_{h}$ we have
that

$p_{h}=\sum_{i=1}^{\mathrm{d}_{\mathcal{P}}}\hat{p}_{i}\psi_{i}\,,\quad\boldsymbol{w}_{h}=\sum_{i=1}^{\mathrm{d}_{\mathcal{W}}}\hat{w}_{i}\boldsymbol{\alpha}_{i}\,,\quad\boldsymbol{u}_{h}=\sum_{i=1}^{\mathrm{d}_{\mathcal{U}}}\hat{u}_{i}\boldsymbol{\beta}_{i}\,,\quad\text{and}\quad
c_{h}=\sum_{i=1}^{\mathrm{d}_{\mathcal{Q}}}\hat{c}_{i}\gamma_{i}\,.$ (9)

We now introduce an equivalent, but more compact, notation for the expansion
of a function as a linear combination of the basis functions

$p_{h}=\psi_{h}\hat{p}_{h}\,,\quad\boldsymbol{w}_{h}=\boldsymbol{\alpha}_{h}\hat{\boldsymbol{w}}_{h}\,,\quad\boldsymbol{u}_{h}=\boldsymbol{\beta}_{h}\hat{\boldsymbol{u}}_{h}\,,\quad\text{and}\quad
c_{h}=\gamma_{h}\hat{c}_{h}\,,$ (10)

where
$\hat{p}_{h}:=[\hat{p}_{1},\dots,\hat{p}_{\mathrm{d}_{\mathcal{P}}}]^{\top}$,
$\hat{\boldsymbol{w}}_{h}:=[\hat{w}_{1},\dots,\hat{w}_{\mathrm{d}_{\mathcal{W}}}]^{\top}$,
$\hat{\boldsymbol{u}}_{h}:=[\hat{u}_{1},\dots,\hat{u}_{\mathrm{d}_{\mathcal{U}}}]^{\top}$,
and
$\hat{c}_{h}:=[\hat{c}_{1},\dots,\hat{c}_{\mathrm{d}_{\mathcal{Q}}}]^{\top}$
are the column vectors containing the coefficients of the discrete
representation of the associated functions, and
$\psi_{h}:=[\psi_{1},\dots,\psi_{\mathrm{d}_{\mathcal{P}}}]$,
$\boldsymbol{\alpha}_{h}:=[\boldsymbol{\alpha}_{1},\dots,\boldsymbol{\alpha}_{\mathrm{d}_{\mathcal{W}}}]$,
$\boldsymbol{\beta}_{h}:=[\boldsymbol{\beta}_{1},\dots,\boldsymbol{\beta}_{\mathrm{d}_{\mathcal{U}}}]$,
and $\gamma_{h}:=[\gamma_{1},\dots,\gamma_{\mathrm{d}_{\mathcal{Q}}}]$ are the
row vectors containing the basis functions.

###### Remark 1

Throughout this work we will extensively use the basis functions $\psi_{h}$,
$\boldsymbol{\alpha}_{h}$, $\boldsymbol{\beta}_{h}$, and $\gamma_{h}$. For
clarity in the notation, these symbols will be solely used to refer to the
basis functions:

$\displaystyle\psi_{h}:=[\psi_{1},\dots,\psi_{\mathrm{d}_{\mathcal{P}}}]\qquad\mathrm{such\
that}$
$\displaystyle\qquad\mathrm{span}\\{\psi_{1},\dots,\psi_{\mathrm{d}_{\mathcal{P}}}\\}=\mathcal{P}_{h}$
(11)
$\displaystyle\boldsymbol{\alpha}_{h}:=[\boldsymbol{\alpha}_{1},\dots,\boldsymbol{\alpha}_{\mathrm{d}_{\mathcal{W}}}]\qquad\mathrm{such\
that}$
$\displaystyle\qquad\mathrm{span}\\{\boldsymbol{\alpha}_{1},\dots,\boldsymbol{\alpha}_{\mathrm{d}_{\mathcal{W}}}\\}=\mathcal{W}_{h}$
$\displaystyle\boldsymbol{\beta}_{h}:=[\boldsymbol{\beta}_{1},\dots,\boldsymbol{\beta}_{\mathrm{d}_{\mathcal{U}}}]\qquad\mathrm{such\
that}$
$\displaystyle\qquad\mathrm{span}\\{\boldsymbol{\beta}_{1},\dots,\boldsymbol{\beta}_{\mathrm{d}_{\mathcal{U}}}\\}=\mathcal{U}_{h}$
$\displaystyle\gamma_{h}:=[\gamma_{1},\dots,\gamma_{\mathrm{d}_{\mathcal{Q}}}]\qquad\mathrm{such\
that}$
$\displaystyle\qquad\mathrm{span}\\{\gamma_{1},\dots,\gamma_{\mathrm{d}_{\mathcal{Q}}}\\}=\mathcal{Q}_{h}\,.$

These discrete function spaces satisfy a compatibility relation that is
expressed as a discrete de Rham complex of the form

$\mathbb{R}\longrightarrow\mathcal{P}_{h}\stackrel{{\scriptstyle\nabla}}{{\longrightarrow}}\mathcal{W}_{h}\stackrel{{\scriptstyle\nabla\times}}{{\longrightarrow}}\mathcal{U}_{h}\stackrel{{\scriptstyle\nabla\cdot}}{{\longrightarrow}}\mathcal{Q}_{h}\longrightarrow
0.$ (12)

The differential operators $\nabla$, $\nabla\times$, and $\nabla\cdot$ are
concisely represented at the discrete level by incidence matrices,
$\boldsymbol{\mathsf{E}}^{1,0}$, $\boldsymbol{\mathsf{E}}^{2,1}$,
$\boldsymbol{\mathsf{E}}^{3,2}$, such that for $p_{h}\in\mathcal{P}_{h}$,
$\boldsymbol{w}_{h}\in\mathcal{W}_{h}$, and
$\boldsymbol{u}_{h}\in\mathcal{U}_{h}$ we have that

$\nabla
p_{h}=\boldsymbol{\alpha}_{h}\boldsymbol{\mathsf{E}}^{1,0}\hat{p}_{h},\quad\nabla\times\boldsymbol{w}_{h}=\boldsymbol{\beta}_{h}\boldsymbol{\mathsf{E}}^{2,1}\hat{\boldsymbol{w}}_{h},\quad\text{and}\quad\nabla\cdot\boldsymbol{u}_{h}=\gamma_{h}\boldsymbol{\mathsf{E}}^{3,2}\hat{\boldsymbol{u}}_{h}\,.$
(13)

In a similar way, we may apply the differential operators to the basis
functions $\psi_{h}$, $\boldsymbol{\alpha}_{h}$, $\boldsymbol{\beta}_{h}$, and
$\gamma_{h}$ yielding

$\nabla\psi_{h}=\boldsymbol{\alpha}_{h}\boldsymbol{\mathsf{E}}^{1,0},\quad\nabla\times\boldsymbol{\alpha}_{h}=\boldsymbol{\beta}_{h}\boldsymbol{\mathsf{E}}^{2,1},\quad\text{and}\quad\nabla\cdot\boldsymbol{\beta}_{h}=\gamma_{h}\boldsymbol{\mathsf{E}}^{3,2}\,.$
(14)

These incidence matrices satisfy the identities

$\boldsymbol{\mathsf{E}}^{2,1}\boldsymbol{\mathsf{E}}^{1,0}=\boldsymbol{0}^{2,0},\quad\text{and}\quad\boldsymbol{\mathsf{E}}^{3,2}\boldsymbol{\mathsf{E}}^{2,1}=\boldsymbol{0}^{3,1}\,,$
(15)

where $\boldsymbol{0}^{2,0}$ is the
$\mathrm{d}_{\mathcal{U}}\times\mathrm{d}_{\mathcal{P}}$ zero matrix and
$\boldsymbol{0}^{3,1}$ is the
$\mathrm{d}_{\mathcal{Q}}\times\mathrm{d}_{\mathcal{W}}$ zero matrix. Note
that these two identities are directly related to the well known vector
calculus identities $\nabla\times\nabla=0$ and $\nabla\cdot\nabla\times=0$.
For a more detailed discussion of the incidence matrices and their properties
the reader is directed to [14].

Additionally, these compatible function spaces satisfy an integration-by-parts
property, such that the reverse mappings ($\tilde{\nabla}\cdot$,
$\tilde{\nabla}\times$, $\tilde{\nabla}$), satisfy

$0\longleftarrow\mathcal{P}_{h}\stackrel{{\scriptstyle\tilde{\nabla}\cdot}}{{\longleftarrow}}\mathcal{W}_{h}\stackrel{{\scriptstyle\tilde{\nabla}\times}}{{\longleftarrow}}\mathcal{U}_{h}\stackrel{{\scriptstyle\tilde{\nabla}}}{{\longleftarrow}}\mathcal{Q}_{h}\longleftarrow\mathbb{R}.$
(16)

and are represented by weak form adjoint relations. Assuming periodic or
homogeneous boundary conditions we have that for
$\boldsymbol{w}_{h}\in\mathcal{W}_{h}$,
$\boldsymbol{u}_{h}\in\mathcal{U}_{h}$, and $c_{h}\in\mathcal{Q}_{h}$

$\displaystyle\langle\boldsymbol{\beta}_{h},\tilde{\nabla}c_{h}\rangle$
$\displaystyle=-\langle\nabla\cdot\boldsymbol{\beta}_{h},c_{h}\rangle,$ (17)
$\displaystyle\langle\boldsymbol{\alpha}_{h},\tilde{\nabla}\times\boldsymbol{u}_{h}\rangle$
$\displaystyle=\langle\nabla\times\boldsymbol{\alpha}_{h},\boldsymbol{u}_{h}\rangle,$
$\displaystyle\langle\psi_{h},\tilde{\nabla}\cdot\boldsymbol{w}_{h}\rangle$
$\displaystyle=-\langle\nabla\psi_{h},\boldsymbol{w}_{h}\rangle\,.$

###### Remark 2

In (17) we have used an implicit matrix notation. As an example, consider the
expression $\langle\boldsymbol{\beta}_{h},\boldsymbol{\beta}_{h}\rangle$. This
expression corresponds to

$\langle\boldsymbol{\beta}_{h},\boldsymbol{\beta}_{h}\rangle:=\int_{\Omega}\boldsymbol{\beta}_{h}^{\top}\boldsymbol{\beta}_{h}\,\mathrm{d}\Omega\,,$
(18)

which is a matrix of dimensions
$\mathrm{d}_{\mathcal{U}}\times\mathrm{d}_{\mathcal{U}}$ with element $(i,j)$,
$\langle\boldsymbol{\beta}_{h},\boldsymbol{\beta}_{h}\rangle_{i,j}$, given by

$\langle\boldsymbol{\beta}_{h},\boldsymbol{\beta}_{h}\rangle_{i,j}=\langle\boldsymbol{\beta}_{i},\boldsymbol{\beta}_{j}\rangle\,.$

As a more complex example, consider the expression
$\langle\nabla\psi_{h},\boldsymbol{w}_{h}\rangle$, with
$\boldsymbol{w}_{h}\in\mathcal{W}_{h}$. In a similar way, this expression
corresponds to

$\langle\nabla\psi_{h},\boldsymbol{w}_{h}\rangle=\int_{\Omega}\nabla\psi_{h}^{\top}\boldsymbol{w}_{h}\,\mathrm{d}\Omega\,,$

which is a column vector of $\mathrm{d}_{\mathcal{W}}$ elements with element
$j$, $\langle\nabla\psi_{h},\boldsymbol{w}_{h}\rangle_{j}$, given by [17]

$\langle\nabla\psi_{h},\boldsymbol{w}_{h}\rangle_{j}\stackrel{{\scriptstyle\eqref{eq:matrix_inner_product}}}{{=}}\left(\int_{\Omega}\nabla\psi_{h}^{\top}\boldsymbol{w}_{h}\right)_{j}\stackrel{{\scriptstyle\eqref{eq_incidence_matrices_basis}}}{{=}}\left(\int_{\Omega}\left(\boldsymbol{\mathsf{E}}^{1,0}\right)^{\top}\boldsymbol{\alpha}^{\top}_{h}\boldsymbol{w}_{h}\right)_{j}=\left(\left(\boldsymbol{\mathsf{E}}^{1,0}\right)^{\top}\langle\boldsymbol{\alpha}_{h},\boldsymbol{w}_{h}\rangle\right)_{j}=\sum_{k=1}^{\mathrm{d}_{\mathcal{W}}}\boldsymbol{\mathsf{E}}^{1,0}_{k,j}\langle\boldsymbol{\alpha}_{k},\boldsymbol{w}_{h}\rangle\,.$

This leads to the following discrete representations of the dual operators
applied to $\boldsymbol{w}_{h}\in\mathcal{P}_{h}$,
$\boldsymbol{u}_{h}\in\mathcal{U}_{h}$, and $c_{h}\in\mathcal{Q}_{h}$

$\displaystyle\tilde{\nabla}c_{h}$
$\displaystyle=-\boldsymbol{\beta}_{h}\,\langle\boldsymbol{\beta}_{h},\boldsymbol{\beta}_{h}\rangle^{-1}\left(\boldsymbol{\mathsf{E}}^{3,2}\right)^{\top}\langle\gamma_{h},\gamma_{h}\rangle\,\hat{c}_{h}\,,$
(19) $\displaystyle\tilde{\nabla}\times\boldsymbol{u}_{h}$
$\displaystyle=\boldsymbol{\alpha}_{h}\,\langle\boldsymbol{\alpha}_{h},\boldsymbol{\alpha}_{h}\rangle^{-1}\left(\boldsymbol{\mathsf{E}}^{2,1}\right)^{\top}\langle\boldsymbol{\beta}_{h},\boldsymbol{\beta}_{h}\rangle\,\hat{\boldsymbol{u}}_{h}\,,$
$\displaystyle\tilde{\nabla}\cdot\boldsymbol{w}_{h}$
$\displaystyle=-\psi_{h}\,\langle\psi_{h},\psi_{h}\rangle^{-1}\left(\boldsymbol{\mathsf{E}}^{1,0}\right)^{\top}\langle\boldsymbol{\alpha}_{h},\boldsymbol{\alpha}_{h}\rangle\,\hat{\boldsymbol{w}}_{h}\,.$

The discrete form of the variational 3D compressible Euler equations (4) may
then be expressed as: Given
$\left.\boldsymbol{u}_{h}\right|_{t=0}\in\mathcal{U}$,
$\left.\rho\right|_{t=0},\left.\Theta\right|_{t=0}\in\mathcal{Q}$, find
$\boldsymbol{u}_{h}\in\mathcal{U}$, $\rho,\Theta\in\mathcal{Q}$ such that

$\displaystyle\Bigg{\langle}\boldsymbol{\beta}_{h},\frac{\partial\boldsymbol{u}_{h}}{\partial
t}\Bigg{\rangle}$
$\displaystyle=-\langle\boldsymbol{\beta}_{h},\boldsymbol{q}_{h}\times\boldsymbol{U}_{h}\rangle+\langle\nabla\cdot\boldsymbol{\beta}_{h},\Phi_{h}\rangle+\langle\nabla\cdot\mathbb{P}^{\mathrm{div}}\left[\theta_{h}\boldsymbol{\beta}_{h}\right],\Pi_{h}\rangle\,,$
(20a) $\displaystyle\Bigg{\langle}\gamma_{h},\frac{\partial\rho_{h}}{\partial
t}\Bigg{\rangle}$
$\displaystyle=-\langle\gamma_{h},\nabla\cdot\boldsymbol{U}_{h}\rangle\,,$
(20b)
$\displaystyle\Bigg{\langle}\gamma_{h},\frac{\partial\Theta_{h}}{\partial
t}\Bigg{\rangle}$
$\displaystyle=-\langle\gamma_{h},\nabla\cdot\mathbb{P}^{\mathrm{div}}\left[\theta_{h}\boldsymbol{U}_{h}\right]\rangle\,,$
(20c)

and $\boldsymbol{U}_{h}\in\mathcal{U}$, $\Phi_{h},\Pi_{h}\in\mathcal{Q}$ such
that

$\displaystyle\langle\boldsymbol{\beta}_{h},\frac{\delta
H_{h}}{\delta\boldsymbol{u}_{h}}\rangle$
$\displaystyle=\langle\boldsymbol{\beta}_{h},\boldsymbol{U}_{h}\rangle=\langle\boldsymbol{\beta}_{h},\rho_{h}\boldsymbol{u}_{h}\rangle\,,$
(21a) $\displaystyle\langle\gamma_{h},\frac{\delta
H_{h}}{\delta\rho_{h}}\rangle$
$\displaystyle=\langle\gamma_{h},\Phi_{h}\rangle=\langle\gamma_{h},\frac{1}{2}\boldsymbol{u}_{h}\cdot\boldsymbol{u}_{h}+gz\rangle\,,$
(21b) $\displaystyle\langle\gamma_{h},\frac{\delta
H_{h}}{\delta\Theta_{h}}\rangle$
$\displaystyle=\langle\gamma_{h},\Pi_{h}\rangle=\langle\gamma_{h},c_{p}\Bigg{(}\frac{R\Theta_{h}}{p_{0}}\Bigg{)}^{\frac{R}{c_{v}}}\rangle.$
(21c)

The system of equations (20), obtained after spatial discretisation,
constitutes a system of ordinary differential equations. We now wish to devise
a discrete time integration formulation that preserves the energy conservation
properties of the spatial discretisation in (20). To introduce the energy
conserving time integration scheme employed in this work, we recall the total
energy $H$, (2),

${H}_{h}=\int_{\Omega}\overbrace{\frac{1}{2}\rho_{h}\boldsymbol{u}_{h}\cdot\boldsymbol{u}_{h}}^{{K_{h}}}+\overbrace{\vphantom{\frac{1}{2}}\rho_{h}gz}^{{P_{h}}}+\overbrace{\frac{c_{v}}{c_{p}}\Theta_{h}\Pi_{h}}^{{I_{h}}}\mathrm{d}\Omega,$
(22)

and its time derivative, (7),

$\Bigg{\langle}\frac{\delta
H_{h}}{\delta\boldsymbol{u}_{h}},\frac{\partial\boldsymbol{u}_{h}}{\partial
t}\Bigg{\rangle}+\Bigg{\langle}\frac{\delta
H_{h}}{\delta\rho_{h}},\frac{\partial\rho_{h}}{\partial
t}\Bigg{\rangle}+\Bigg{\langle}\frac{\delta
H_{h}}{\delta\Theta_{h}},\frac{\partial\Theta_{h}}{\partial
t}\Bigg{\rangle}=\Bigg{\langle}\boldsymbol{U}_{h},\frac{\partial\boldsymbol{u}_{h}}{\partial
t}\Bigg{\rangle}+\Bigg{\langle}\Phi_{h},\frac{\partial\rho_{h}}{\partial
t}\Bigg{\rangle}+\Bigg{\langle}\Pi_{h},\frac{\partial\Theta_{h}}{\partial
t}\Bigg{\rangle}=\frac{\mathrm{d}H_{h}}{\mathrm{d}t}=0.$ (23)

We now note that (20) can be written in matrix form as

$\left[\begin{array}[]{>{\displaystyle}c}\langle\boldsymbol{\beta}_{h},\boldsymbol{\beta}_{h}\rangle\frac{\mathrm{d}\hat{\boldsymbol{u}}_{h}}{\mathrm{d}t}\\\
\langle\gamma_{h},\gamma_{h}\rangle\frac{\mathrm{d}\hat{\rho}_{h}}{\mathrm{d}t}\\\
\langle\gamma_{h},\gamma_{h}\rangle\frac{\mathrm{d}\hat{\Theta}_{h}}{\mathrm{d}t}\end{array}\right]=\left[\begin{array}[]{ccc}-\langle\boldsymbol{\beta}_{h},\boldsymbol{q}_{h}\times\boldsymbol{\beta}_{h}\rangle&\langle\nabla\cdot\boldsymbol{\beta}_{h},\gamma_{h}\rangle&\langle\nabla\cdot\mathbb{P}^{\mathrm{div}}\left[\theta_{h}\boldsymbol{\beta}_{h}\right],\gamma_{h}\rangle\\\
-\langle\gamma_{h},\nabla\cdot\boldsymbol{\beta}_{h}\rangle&0&0\\\
-\langle\gamma_{h},\nabla\cdot\mathbb{P}^{\mathrm{div}}\left[\theta_{h}\boldsymbol{\beta}_{h}\right]\rangle&0&0\end{array}\right]\left[\begin{array}[]{c}\hat{\boldsymbol{U}}_{h}\\\
\hat{\Phi}_{h}\\\ \hat{\Pi}_{h}\end{array}\right]\,.$ (24)

###  2.2 Temporal discretisation

Equation (24) is essentially a system of ordinary differential equations
(ODEs) with the general form

$\boldsymbol{\mathsf{M}}\frac{\mathrm{d}\hat{\boldsymbol{y}}(t)}{\mathrm{d}t}=\boldsymbol{\mathsf{S}}(\hat{\boldsymbol{y}}(t))\,\hat{\boldsymbol{x}}(t),$
(25)

where

$\hat{\boldsymbol{y}}(t):=\left[\begin{array}[]{c}\hat{\boldsymbol{u}}_{h}\\\
\hat{\rho}_{h}\\\
\hat{\Theta}_{h}\end{array}\right]\,,\qquad\hat{\boldsymbol{x}}(t):=\left[\begin{array}[]{c}\hat{\boldsymbol{U}}_{h}\\\
\hat{\Phi}_{h}\\\ \hat{\Pi}_{h}\end{array}\right]\,,$ (26)

$\boldsymbol{\mathsf{S}}(\hat{\boldsymbol{y}}(t))$ is a time dependent skew-
symmetric matrix and $\boldsymbol{\mathsf{M}}$ is a constant in time symmetric
positive definite matrix. The variable $\hat{\boldsymbol{y}}$ (26) is a column
vector containing the coefficients of the discrete representation of the
fields $\boldsymbol{u}_{h}$, $\rho_{h}$, and $\Theta_{h}$, as introduced in
(9) and (10). In this way we have that

$\hat{y}_{i}=\hat{u}_{i},\quad\mathrm{for}\quad
i=1,\dots,\mathrm{d}_{\mathcal{U}},\quad\hat{y}_{i+\mathrm{d}_{\mathcal{U}}}=\hat{\rho}_{i},\quad\mathrm{for}\quad
i=1,\dots,\mathrm{d}_{\mathcal{Q}},\quad\mathrm{and}\quad\hat{y}_{i+\mathrm{d}_{\mathcal{U}}+\mathrm{d}_{\mathcal{Q}}}=\hat{\Theta}_{i},\quad\mathrm{for}\quad
i=1,\dots,\mathrm{d}_{\mathcal{Q}}\,.$ (27)

These coefficients $\hat{\boldsymbol{y}}$ are time dependent and their
evolution is governed by the system of ODEs (25).

To solve this system of ODEs in time, we follow a similar approach to that for
the spatial discretisation. We approximate the time dependent coefficients
$\hat{\boldsymbol{y}}(t)$ on the time step interval $[t^{n},t^{n+1}]$ with a
temporal discrete polynomial function space of degree $s$,
$\mathcal{T}^{0}_{h}([t^{n},t^{n+1}]):=\mathrm{span}\\{l_{0},\dots,l_{s}\\}$
where $l_{i}(t)$, $i=0,\dots s$, are the ($s+1$) nodal polynomials of degree
$s$, as mentioned before for the spatial discretisation. In this way the
polynomial temporal expansion of the spatial coefficients $\hat{y}_{i}^{h}$ is

$\hat{y}_{i}^{h}=l_{h}\,\hat{\hat{y}}_{i}^{h}:=\sum_{j=0}^{s}l_{j}(t)\hat{\hat{y}}^{j}_{i}\,,\qquad
t\in([t^{n},t^{n+1}])\,,$ (28)

where, as in (9) and (10),
$\hat{\hat{y}}_{i}^{h}:=[\hat{\hat{y}}_{i}^{1},\dots\hat{\hat{y}}_{i}^{s+1}]$.
The double hat notation $\hat{\hat{y}}_{i}^{h}$ in (28) highlights the double
discretization taking place (spatial and temporal). The terms
$\hat{y}_{i}^{h}$ with
$i=1,\dots,\mathrm{d}_{\mathcal{U}}+2\mathrm{d}_{\mathcal{Q}}$ correspond to
the time dependent polynomial expansion of the $i$-th coefficient of the
spatial polynomial expansion
$\hat{\boldsymbol{y}}:=[\hat{\boldsymbol{u}}_{h},\hat{\rho}_{h},\hat{\Theta}_{h}]^{\top}$.
The terms $\hat{\hat{y}}^{j}_{i}$ correspond to the $j$-th coefficient of the
temporal polynomial expansion associated to the $i$-th coefficient of the
spatial polynomial expansion.

The time dependent coefficients $\hat{\boldsymbol{x}}(t)$ are approximated
with a temporal discrete polynomial function space of degree ($s-1$),
$\mathcal{T}^{1}_{h}([t^{n},t^{n+1}]):=\mathrm{span}\\{e_{1},\dots,e_{s}\\}$
where $e_{i}(t)$, $i=1,\dots s$, are the $s$ edge polynomials of degree
($s-1$), as mentioned before for the spatial discretisation. In this way

$\hat{x}_{i}^{h}=e_{h}\,\hat{\hat{x}}_{i}^{h}:=\sum_{j=1}^{s}e_{j}(t)\hat{\hat{x}}^{j}_{i}\,,\qquad
t\in([t^{n},t^{n+1}])\,.$ (29)

As before, these two function spaces constitute a discrete de Rham complex

$\mathbb{R}\longrightarrow\mathcal{T}^{0}_{h}([t^{n},t^{n+1}])\stackrel{{\scriptstyle\frac{\mathrm{d}}{\mathrm{d}t}}}{{\longrightarrow}}\mathcal{T}^{1}_{h}([t^{n},t^{n+1}])\longrightarrow
0\,,$ (30)

and we have that

$\frac{\mathrm{d}\hat{y}_{i}^{h}}{\mathrm{d}t}=e_{h}\boldsymbol{\mathsf{E}}^{1,0}_{t}\hat{\hat{y}}_{i}^{h}\,,\qquad\mathrm{and}\qquad\frac{\mathrm{d}l_{h}}{\mathrm{d}t}=e_{h}\boldsymbol{\mathsf{E}}^{1,0}_{t}\,.$
(31)

###### Remark 3

Note that here we added the subscript $t$ to the incidence matrix,
$\boldsymbol{\mathsf{E}}^{1,0}_{t}$, to distinguish it from the spatial ones.
Also, as mentioned before, the double hat highlights the fact that we have
performed both discretisation in space and time and these coefficients are
associated to this double discretisation. When an explicit reference to the
temporal indices is required they are presented as superscript, as in (28).

###### Remark 4

The discrete de Rham complex (30) is associated to a single time step (for
simplicity of exposition and computational efficiency). Naturally, as for the
spatial discretisation, considering multiple time steps is a straightforward
extension. In this work we perform the discretisation in a per step fashion.

The temporal discretisation then follows a similar mixed Galerkin approach as
presented for the spatial discretisation: Given
$\hat{x}_{i}^{h}\in\mathcal{T}^{1}_{h}([t^{n},t^{n+1}])$, with
$i=1,\dots,(\mathrm{d}_{\mathcal{U}}+2\mathrm{d}_{\mathcal{Q}})$, find
$\hat{y}_{i}^{h}\in\mathcal{T}^{0}_{h}([t^{n},t^{n+1}])$, with
$i=1,\dots,(\mathrm{d}_{\mathcal{U}}+2\mathrm{d}_{\mathcal{Q}})$, such that

$\boldsymbol{\mathsf{M}}\,\langle
e_{h},\frac{\mathrm{d}\hat{\boldsymbol{y}}_{h}(t)}{\mathrm{d}t}\rangle_{[t^{n},t^{n+1}]}=\langle
e_{h},\boldsymbol{\mathsf{S}}(\hat{\boldsymbol{y}}_{h}(t))\,\hat{\boldsymbol{x}}_{h}(t)\rangle_{[t^{n},t^{n+1}]}\,.$
(32)

###### Remark 5

The time-related inner products follow the same rules outlined before, namely

$\langle
e_{h},\frac{\mathrm{d}\hat{\boldsymbol{y}}_{h}(t)}{\mathrm{d}t}\rangle_{[t^{n},t^{n+1}]}:=\int_{t^{n}}^{t^{n+1}}e_{h}^{\top}\frac{\mathrm{d}\hat{\boldsymbol{y}}_{h}(t)}{\mathrm{d}t}\,\mathrm{d}t\,.$
(33)

Therefore, (32) is

$\boldsymbol{\mathsf{M}}\,\int_{t^{n}}^{t^{n+1}}e_{h}^{\top}\frac{\mathrm{d}\hat{\boldsymbol{y}}_{h}(t)}{\mathrm{d}t}\,\mathrm{d}t=\int_{t^{n}}^{t^{n+1}}e_{h}^{\top}\boldsymbol{\mathsf{S}}(\hat{\boldsymbol{y}}_{h}(t))\,\hat{\boldsymbol{x}}_{h}(t)\,\mathrm{d}t\,.$
(34)

Note that we implicitly separate the spatial from the temporal coefficients,
for simplicity of the notation, i.e., the constant in time matrix
$\boldsymbol{\mathsf{M}}$ and the time dependent matrix
$\boldsymbol{\mathsf{S}}(\hat{\boldsymbol{y}}_{h}(t))$ apply at the spatial
level also.

This choice of time integration is notable because it is energy conserving. To
prove this, we recall the expression for the time variation of the energy of
this system of equations (the Hamiltonian, $H_{h}$) (23)

$\int_{t^{n}}^{t^{n+1}}\frac{\mathrm{d}H_{h}}{\mathrm{d}t}\,\mathrm{d}t=H_{h}(\boldsymbol{y}_{h}(t^{n+1}))-H_{h}(\boldsymbol{y}_{h}(t^{n}))=\int_{t^{n}}^{t^{n+1}}\Bigg{\langle}\frac{\delta{H_{h}}}{\delta\boldsymbol{y}_{h}},\frac{\partial\boldsymbol{y}_{h}}{\partial
t}\Bigg{\rangle}\,\mathrm{d}t.$ (35)

We may now use (21) to rewrite (35) as

$\int_{t^{n}}^{t^{n+1}}\frac{\mathrm{d}H_{h}}{\mathrm{d}t}\,\mathrm{d}t=H_{h}(\boldsymbol{y}_{h}(t^{n+1}))-H_{h}(\boldsymbol{y}_{h}(t^{n}))=\int_{t^{n}}^{t^{n+1}}\Bigg{\langle}\boldsymbol{x}_{h},\frac{\partial\boldsymbol{y}_{h}}{\partial
t}\Bigg{\rangle}\,\mathrm{d}t=\int_{t^{n}}^{t^{n+1}}\hat{\boldsymbol{x}}_{h}^{\top}\boldsymbol{\mathsf{M}}\frac{\mathrm{d}\hat{\boldsymbol{y}}_{h}}{\mathrm{d}t}\,\mathrm{d}t=\hat{\hat{\boldsymbol{x}}}_{h}^{\top}\boldsymbol{\mathsf{M}}\int_{t^{n}}^{t^{n+1}}e_{h}^{\top}\frac{\mathrm{d}\hat{\boldsymbol{y}}_{h}}{\mathrm{d}t}\,\mathrm{d}t.$
(36)

Finally, using the proposed time integration scheme, (34), we obtain energy
conservation

$H_{h}(\boldsymbol{y}_{h}(t^{n+1}))-H_{h}(\boldsymbol{y}_{h}(t^{n}))=\hat{\hat{\boldsymbol{x}}}_{h}^{\top}\boldsymbol{\mathsf{M}}\int_{t^{n}}^{t^{n+1}}e_{h}^{\top}e_{h}\,\mathrm{d}t(\hat{\hat{\boldsymbol{y}}}_{h}^{n+1}-\hat{\hat{\boldsymbol{y}}}_{h}^{n})=\hat{\hat{\boldsymbol{x}}}_{h}^{\top}\int_{t^{n}}^{t^{n+1}}e_{h}^{\top}\boldsymbol{\mathsf{S}}(\hat{\boldsymbol{y}}_{h}(t))\,\hat{\boldsymbol{x}}_{h}(t)\,\mathrm{d}t=\hat{\hat{\boldsymbol{x}}}_{h}^{\top}\left(\int_{t^{n}}^{t^{n+1}}e_{h}^{\top}\boldsymbol{\mathsf{S}}(\hat{\boldsymbol{y}}_{h}(t))\,e_{h}\,\mathrm{d}t\right)\,\hat{\hat{\boldsymbol{x}}}_{h}.$
(37)

Since the term inside the parenthesis on the right hand side is skew-
symmetric, the quadratic term cancels and we obtain energy conservation at the
discrete level

$H_{h}(\boldsymbol{y}_{h}(t^{n+1}))-H_{h}(\boldsymbol{y}_{h}(t^{n}))=0.$ (38)

###### Remark 6

In equation (37), the term
$\int_{t^{n}}^{t^{n+1}}e_{h}^{\top}e_{h}\,\mathrm{d}t$ is a temporal mass
matrix.

###### Remark 7

Two key ingredients need to be highlighted. The first is the construction of a
system of equations of the form (24) and (25) where $\boldsymbol{\mathsf{S}}$
is skew-symmetric. This was only possible to obtain by a judicious choice of
discrete function spaces and an equally careful selection of physical field
quantities to employ in the construction of the system of equations. The
second key ingredient is the construction of the time integration scheme also
based on a sequence of polynomial spaces together with a Galerkin projection
employing exact integration. Exact temporal integration is fundamental to
guarantee the equality needed to perform the last step, equation (37).

This temporal discretisation is very closely related to the works by Hairer et
al. [18, 19]. For the lowest order case (discussed below) this approach
results in time stepping scheme identical to [19].

In this formulation any polynomial degree $s$ may be employed to approximate
the solution in time. In this work we choose $s=1$ since this choice will
greatly simplify the temporal discretisation. If this choice is made we will
have

$l_{0}(t)=\frac{t^{n+1}-t}{\Delta t}\,,\quad l_{1}(t)=\frac{t-t^{n}}{\Delta
t}\,,\quad\mathrm{and}\quad e_{1}(t)=\frac{1}{\Delta t}\,.$ (39)

where $\Delta t:=t^{n+1}-t^{n}$. In this case, since
$e_{h}=e_{1}(t)=\frac{1}{\Delta t}$, the coefficients of the temporal
polynomial expansion of $\hat{\boldsymbol{x}}_{h}(t)$ become

$\hat{\hat{\boldsymbol{x}}}_{h}(t)=[\hat{\hat{U}}^{1}_{1},\dots,\hat{\hat{U}}^{1}_{\mathrm{d}_{\mathcal{U}}},\hat{\hat{\Phi}}^{1}_{1},\dots,\hat{\hat{\Phi}}^{1}_{\mathrm{d}_{\mathcal{Q}}},\hat{\hat{\Pi}}^{1}_{1},\dots,\hat{\hat{\Pi}}^{1}_{\mathrm{d}_{\mathcal{Q}}}]^{\top}\,e_{1}(t)\,.$
(40)

We recall now the properties of the edge polynomials, [12, 16], and its
associated coefficients. Just as the nodal polynomial expansion interpolates a
function, the edge polynomial expansion histopolates (interpolates the
integral of) a function. While the coefficients of the nodal polynomial
expansion correspond to pointwise sampling of a function, the coefficients of
the edge polynomial expansion correspond to sampling of the integral of the
function over intervals [12, 16]. Note that in the temporal discretization
considered here there is only one polynomial basis in the edge polynomial
expansion and therefore the temporal expansion of each spatial coefficient has
only one coefficient. Given the integral interpolation properties of the edge
basis functions, the single coefficient of this expansion is the temporal
integral

$\hat{\hat{\boldsymbol{x}}}_{i}^{1}=\hat{\hat{U}}^{1}_{i}=\int_{t^{n}}^{t^{n+1}}\hat{U}_{i}^{h}(t)\,\mathrm{d}t,\quad\hat{\hat{\boldsymbol{x}}}_{j+\mathrm{d}_{\mathcal{U}}}^{1}=\hat{\hat{\Phi}}^{1}_{j}=\int_{t^{n}}^{t^{n+1}}\hat{\Phi}_{j}^{h}(t)\,\mathrm{d}t,\quad\hat{\hat{\boldsymbol{x}}}_{k+\mathrm{d}_{\mathcal{U}}+\mathrm{d}_{\mathcal{Q}}}^{1}=\hat{\hat{\Pi}}^{1}_{k}=\int_{t^{n}}^{t^{n+1}}\hat{\Pi}_{k}^{h}(t)\,\mathrm{d}t\,,$
(41)

with $i=1,\dots,\mathrm{d}_{\mathcal{U}}$, and
$j,k=1,\dots,\mathrm{d}_{\mathcal{Q}}$. We now introduce

$\overline{\hat{\hat{\boldsymbol{x}}}_{i}^{h}}:=\frac{\hat{\hat{\boldsymbol{x}}}_{i}^{1}}{\Delta
t}=\frac{\int_{t^{n}}^{t^{n+1}}\hat{\hat{\boldsymbol{x}}}_{i}^{h}\,\mathrm{d}t}{\Delta
t}\,,$ (42)

the exact time integral of the temporal approximation between time levels $n$
and $n+1$. Following the notation introduced before, (9), we will use

$\overline{\hat{\hat{\boldsymbol{x}}}}_{h}:=[\overline{\hat{\hat{\boldsymbol{x}}}_{1}^{h}},\dots,\overline{\hat{\hat{\boldsymbol{x}}}^{h}}_{\mathrm{d}_{\mathcal{U}}+2\mathrm{d}_{\mathcal{Q}}}]^{\top}\stackrel{{\scriptstyle\eqref{eq::x_coefficients}+\eqref{eq:mean_01}+\eqref{eq:mean_02}}}{{=}}\frac{1}{\Delta
t}[\hat{\hat{U}}^{1}_{1},\dots,\hat{\hat{U}}^{1}_{\mathrm{d}_{\mathcal{U}}},\hat{\hat{\Phi}}^{1}_{1},\dots,\hat{\hat{\Phi}}^{1}_{\mathrm{d}_{\mathcal{Q}}},\hat{\hat{\Phi}}^{1}_{1},\dots,\hat{\hat{\Pi}}^{1}_{\mathrm{d}_{\mathcal{Q}}}]^{\top}\,.$
(43)

If we substitute this into (34) we obtain

$\boldsymbol{\mathsf{M}}\,\left(\hat{\hat{\boldsymbol{y}}}^{n+1}_{h}-\hat{\hat{\boldsymbol{y}}}^{n}_{h}\right)=\left(\int_{t^{n}}^{t^{n+1}}\boldsymbol{\mathsf{S}}(\hat{\boldsymbol{y}}_{h}(t))\,\mathrm{d}t\right)\,\frac{\hat{\hat{\boldsymbol{x}}}_{h}}{\Delta
t}\stackrel{{\scriptstyle\eqref{eq:mean_02}}}{{=}}\left(\int_{t^{n}}^{t^{n+1}}\boldsymbol{\mathsf{S}}(\hat{\boldsymbol{y}}_{h}(t))\,\mathrm{d}t\right)\,\overline{\hat{\hat{\boldsymbol{x}}}}_{h}\,.$
(44)

Since $\boldsymbol{\mathsf{S}}(\hat{\boldsymbol{y}}_{h}(t))$ depends only
linearly on time, we may use lowest order Gauss integration (midpoint rule) to
exactly integrate the term involving
$\boldsymbol{\mathsf{S}}(\hat{\boldsymbol{y}}_{h}(t))$, yielding

$\boldsymbol{\mathsf{M}}\,\left(\hat{\hat{\boldsymbol{y}}}^{n+1}_{h}-\hat{\hat{\boldsymbol{y}}}^{n}_{h}\right)=\Delta
t\,\boldsymbol{\mathsf{S}}\left(\frac{\hat{\hat{\boldsymbol{y}}}^{n+1}_{h}+\hat{\hat{\boldsymbol{y}}}^{n}_{h}}{2}\right)\,\overline{\hat{\hat{\boldsymbol{x}}}}_{h}\,.$
(45)

Therefore this system of equations may be cast in a fully discrete form as:

$\left[\begin{array}[]{c}\langle\boldsymbol{\beta}_{h},\boldsymbol{\beta}_{h}\rangle\,\hat{\boldsymbol{u}}_{h}^{n+1}\\\
\langle\gamma_{h},\gamma_{h}\rangle\,\hat{\rho}_{h}^{n+1}\\\
\langle\gamma_{h},\gamma_{h}\rangle\,\hat{\Theta}_{h}^{n+1}\end{array}\right]=\left[\begin{array}[]{c}\langle\boldsymbol{\beta}_{h},\boldsymbol{\beta}_{h}\rangle\,\hat{\boldsymbol{u}}_{h}^{n}\\\
\langle\gamma_{h},\gamma_{h}\rangle\,\hat{\rho}_{h}^{n}\\\
\langle\gamma_{h},\gamma_{h}\rangle\,\hat{\Theta}_{h}^{n}\end{array}\right]-\\\
\Delta
t\left[\begin{array}[]{ccc}\langle\boldsymbol{\beta}_{h},\boldsymbol{q}^{n+\frac{1}{2}}_{h}\times\boldsymbol{\beta}_{h}\rangle&-\left(\boldsymbol{\mathsf{E}}^{3,2}\right)^{\top}\langle\gamma_{h},\gamma_{h}\rangle&-\langle\boldsymbol{\beta}_{h},\theta^{n+\frac{1}{2}}_{h}\boldsymbol{\beta}_{h}\rangle\langle\boldsymbol{\beta}_{h},\boldsymbol{\beta}_{h}\rangle^{-1}\left(\boldsymbol{\mathsf{E}}^{3,2}\right)\langle\gamma_{h},\gamma_{h}\rangle\\\
\langle\gamma_{h},\gamma_{h}\rangle\boldsymbol{\mathsf{E}}^{3,2}&0&0\\\
\langle\gamma_{h},\gamma_{h}\rangle\boldsymbol{\mathsf{E}}^{3,2}\langle\boldsymbol{\beta}_{h},\boldsymbol{\beta}_{h}\rangle^{-1}\langle\boldsymbol{\beta}_{h},\theta^{n+\frac{1}{2}}_{h}\boldsymbol{\beta}_{h}\rangle&0&0\end{array}\right]\left[\begin{array}[]{c}\overline{\hat{\boldsymbol{U}}}_{h}\\\
\overline{\hat{\Phi}}_{h}\\\ \overline{\hat{\Pi}}_{h}\end{array}\right]\,.$
(46)

As mentioned before, $\overline{{\hat{\boldsymbol{U}}}}_{h}$,
$\overline{\hat{{\Phi}}}_{h}$, and $\overline{\hat{{\Pi}}}_{h}$ are the exact
time integrals of the temporal polynomial reconstruction of
$\hat{\boldsymbol{U}}_{h}$, $\hat{\Phi}_{h}$, and $\hat{\Pi}_{h}$. An
important point to note is the following. Since $\hat{U}_{h}$,
$\hat{\Phi}_{h}$, and $\hat{\Pi}_{h}$ are functions of $\hat{\boldsymbol{u}}$,
$\hat{\rho}$, and $\hat{\Theta}$, (3), we employ discrete piecewise linear
approximations between time levels $n$ and $n+1$ to be described in the
following section, and
${\theta}^{n+\frac{1}{2}}_{h}=(\theta_{h}^{n}+\theta_{h}^{n+1})/2$,
${\boldsymbol{q}}^{n+\frac{1}{2}}_{h}=(\boldsymbol{q}_{h}^{n}+\boldsymbol{q}_{h}^{n+1})/2$
are the time centered potential temperature and potential vorticity.

##  3 HEVI splitting

The above system (46) may be dimensionally split into an implicit solve for
the vertical dynamics (incorporating the horizontal divergence terms) and
explicit momentum advection by defining the horizontal and vertical velocity
components respectively as $\boldsymbol{v}_{h}\in\mathcal{U}_{h}^{\parallel}$,
$\boldsymbol{w}_{h}\in\mathcal{U}_{h}^{\perp}$, such that
$\boldsymbol{u}_{h}=\boldsymbol{v}_{h}+\boldsymbol{w}_{h}\in\mathcal{U}_{h}$
[8]. Similarly the discrete vorticity vector,
$\boldsymbol{q}_{h}\in\mathcal{W}_{h}$ may be partitioned into its respective
vertical and horizontal components (in global coordinates),
$q_{h}^{\perp}=\boldsymbol{q}_{h}\cdot\boldsymbol{e}_{z}$,
$\boldsymbol{q}_{h}^{\parallel}=(-\boldsymbol{q}_{h}\cdot\boldsymbol{e}_{\phi},\boldsymbol{q}\cdot\boldsymbol{e}_{\lambda})$,
where $\boldsymbol{e}_{\lambda}$, $\boldsymbol{e}_{\phi}$,
$\boldsymbol{e}_{z}$ are unit vectors in the zonal, meridional and vertical
global coordinates. The discrete strong form divergence operator, as given by
the incidence matrix $\boldsymbol{\mathsf{E}}^{3,2}$, may itself be decomposed
into its horizontal, $\boldsymbol{\mathsf{E}}^{3,2}_{\parallel}$ and vertical
$\boldsymbol{\mathsf{E}}^{3,2}_{\perp}$ components, such that

$\nabla\cdot\boldsymbol{u}_{h}\stackrel{{\scriptstyle\eqref{eq_incidence_matrices}}}{{=}}\gamma_{h}\boldsymbol{\mathsf{E}}^{3,2}\hat{\boldsymbol{u}}_{h}=:\gamma_{h}\left(\boldsymbol{\mathsf{E}}^{3,2}_{\parallel}\hat{\boldsymbol{v}}_{h}+\boldsymbol{\mathsf{E}}^{3,2}_{\perp}\hat{\boldsymbol{w}}_{h}\right)\,.$
(47)

The second order HEVI splitting is then given as:

_Step 1: Explicit horizontal momentum solve_

$\langle\boldsymbol{\beta}_{h}^{\parallel},\boldsymbol{\beta}_{h}^{\parallel}\rangle{\hat{v}_{h}}^{\prime}=\langle\boldsymbol{\beta}_{h}^{\parallel},\boldsymbol{\beta}_{h}^{\parallel}\rangle{\hat{v}_{h}}^{n-1}-2\Delta
t\langle\boldsymbol{\beta}_{h}^{\parallel},q^{\perp,n}_{h}\times\boldsymbol{\beta}_{h}^{\parallel}\rangle{\hat{V}_{h}}^{n}+2\Delta
t\langle\boldsymbol{\beta}_{h}^{\parallel}\cdot\boldsymbol{q}_{h}^{\parallel,n},\beta_{h}^{\perp}\rangle\hat{W}_{h}^{n}+\\\
2\Delta
t\Big{(}\boldsymbol{\mathsf{E}}^{3,2}_{\parallel}\Big{)}^{\top}\langle\gamma_{h},\gamma_{h}\rangle{\hat{\Phi}_{h}}^{n}+2\Delta
t\langle\boldsymbol{\beta}_{h}^{\parallel},\theta_{h}^{n}\boldsymbol{\beta}_{h}^{\parallel}\rangle\langle\boldsymbol{\beta}_{h}^{\parallel},\boldsymbol{\beta}_{h}^{\parallel}\rangle^{-1}\Big{(}\boldsymbol{\mathsf{E}}^{3,2}_{\parallel}\Big{)}^{\top}\langle\gamma_{h},\gamma_{h}\rangle{\hat{\Pi}_{h}}^{n}$
(48)

_Step 2: Implicit vertical solve (including horizontal divergence)_

$\begin{bmatrix}\langle\beta_{h}^{\perp},\beta_{h}^{\perp}\rangle{\hat{w}_{h}}^{n+1}\\\
\langle\gamma_{h},\gamma_{h}\rangle{\hat{\rho}_{h}}^{n+1}\\\
\langle\gamma_{h},\gamma_{h}\rangle{\hat{\Theta}_{h}}^{n+1}\\\
\end{bmatrix}=\begin{bmatrix}\langle\beta_{h}^{\perp},\beta_{h}^{\perp}\rangle{\hat{w}_{h}}^{n}\\\
\langle\gamma_{h},\gamma_{h}\rangle{\hat{\rho}_{h}}^{n}\\\
\langle\gamma_{h},\gamma_{h}\rangle{\hat{\Theta}_{h}}^{n}\\\ \end{bmatrix}-\\\
\Delta
t\begin{bmatrix}\boldsymbol{\mathsf{0}}&-\Big{(}\boldsymbol{\mathsf{E}}^{3,2}_{\perp}\Big{)}^{\top}\langle\gamma_{h},\gamma_{h}\rangle&-\langle\beta_{h}^{\perp},{\theta}_{h}^{n+1/2}\beta_{h}^{\perp}\rangle\langle\beta_{h}^{\perp},\beta_{h}^{\perp}\rangle^{-1}\Big{(}\boldsymbol{\mathsf{E}}^{3,2}_{\perp}\Big{)}^{\top}\langle\gamma_{h},\gamma_{h}\rangle\\\
\langle\gamma_{h},\gamma_{h}\rangle\boldsymbol{\mathsf{E}}^{3,2}_{\perp}&\boldsymbol{\mathsf{0}}&\boldsymbol{\mathsf{0}}\\\
\langle\gamma_{h},\gamma_{h}\rangle\boldsymbol{\mathsf{E}}^{3,2}_{\perp}\langle\beta_{h}^{\perp},\beta_{h}^{\perp}\rangle^{-1}\langle\beta_{h}^{\perp},{\theta}_{h}^{n+1/2}\beta_{h}^{\perp}\rangle&\boldsymbol{\mathsf{0}}&\boldsymbol{\mathsf{0}}\\\
\end{bmatrix}\begin{bmatrix}\overline{{\hat{W}}}_{h}\\\
\overline{{\hat{\Phi}}}_{h}\\\ \overline{{\hat{\Pi}}}_{h}\\\ \end{bmatrix}-\\\
\Delta
t\begin{bmatrix}\langle\beta_{h}^{\perp},\boldsymbol{q}_{h}^{\parallel,n+1/2}\cdot\boldsymbol{\beta}_{h}^{\parallel}\rangle\overline{\hat{V}}_{h}\\\
\langle\gamma_{h},\gamma_{h}\rangle\boldsymbol{\mathsf{E}}^{3,2}_{\parallel}\overline{{\hat{V}}}_{h}\\\
\langle\gamma_{h},\gamma_{h}\rangle\boldsymbol{\mathsf{E}}^{3,2}_{\parallel}\langle\boldsymbol{\beta}_{h}^{\parallel},\boldsymbol{\beta}_{h}^{\parallel}\rangle^{-1}\langle\boldsymbol{\beta}_{h}^{\parallel},{\theta}_{h}^{n+1/2}\boldsymbol{\beta}_{h}^{\parallel}\rangle\overline{{\hat{V}}}_{h}\end{bmatrix}$
(49)

_Step 3: Explicit horizontal momentum solve_

$\langle\boldsymbol{\beta}_{h}^{\parallel},\boldsymbol{\beta}_{h}^{\parallel}\rangle{\hat{v}_{h}}^{n+1}=\langle\boldsymbol{\beta}_{h}^{\parallel},\boldsymbol{\beta}_{h}^{\parallel}\rangle{\hat{v}_{h}}^{n}-\Delta
t\langle\boldsymbol{\beta}_{h}^{\parallel},q^{\perp,n+1/2}_{h}\times\boldsymbol{\beta}_{h}^{\parallel}\rangle\overline{\hat{V}}_{h}+\Delta
t\langle\boldsymbol{\beta}_{h}^{\parallel}\cdot\boldsymbol{q}_{h}^{\parallel,n+1/2},\beta_{h}^{\perp}\rangle\overline{\hat{W}}_{h}+\\\
\Delta
t\Big{(}\boldsymbol{\mathsf{E}}^{3,2}_{\parallel}\Big{)}^{\top}\langle\gamma_{h},\gamma_{h}\rangle\overline{{\hat{\Phi}}}_{h}+\Delta
t\langle\boldsymbol{\beta}_{h}^{\parallel},\theta_{h}^{n+1/2}\boldsymbol{\beta}_{h}^{\parallel}\rangle\langle\boldsymbol{\beta}_{h}^{\parallel},\boldsymbol{\beta}_{h}^{\parallel}\rangle^{-1}\Big{(}\boldsymbol{\mathsf{E}}^{3,2}_{\parallel}\Big{)}^{\top}\langle\gamma_{h},\gamma_{h}\rangle\overline{{\hat{\Pi}}}_{h}$
(50)

Left multiplication of (49) by
$[\overline{{\hat{W}}}_{h}^{\top}\quad\overline{{\hat{\Phi}}}_{h}^{\top}\quad\overline{{\hat{\Pi}}}_{h}^{\top}]$
and (50) by $\overline{{\hat{V}}}_{h}^{\top}$ leads to the cancellation of all
forcing terms in both equations, resulting in the expression

$\overline{{\hat{V}}}_{h}\langle\boldsymbol{\beta}_{h}^{\parallel},\boldsymbol{\beta}_{h}^{\parallel}\rangle({\hat{v}}_{h}^{n+1}-{\hat{v}}_{h}^{n})+\overline{{\hat{W}}}_{h}\langle\boldsymbol{\beta}_{h}^{\parallel},\boldsymbol{\beta}_{h}^{\parallel}\rangle({\hat{w}}_{h}^{n+1}-{\hat{w}}_{h}^{n})+\overline{{\hat{\Phi}}}_{h}\langle\gamma_{h},\gamma_{h}\rangle({\hat{\rho}}_{h}^{n+1}-{\hat{\rho}}_{h}^{n})+\overline{{\hat{\Pi}}}_{h}\langle\gamma_{h},\gamma_{h}\rangle({\hat{\Theta}}_{h}^{n+1}-{\hat{\Theta}}_{h}^{n})\neq
0$ (51)

Note that since ${\boldsymbol{v}_{h}}^{\prime}$ and
${\boldsymbol{v}_{h}}^{n+1}$ differ, energy is not strictly conserved in the
temporal discretisation of the horizontal momentum equation. There is a
conservation error of magnitude
$\overline{{\hat{V}}}_{h}\langle\boldsymbol{\beta}_{h}^{\parallel},\boldsymbol{\beta}_{h}^{\parallel}\rangle({\hat{v}}_{h}^{\prime}-{\hat{v}}_{h}^{n+1})$,
such that (51) is not strictly a discrete analogue of the continuous
expression

$\Bigg{(}\frac{\delta{H}}{\delta\boldsymbol{y}}\Bigg{)}^{\top}\cdot\frac{\partial\boldsymbol{y}}{\partial
t}=0.$ (52)

The above formulation is similar to the “u-forward, pressure-backward”
formulation [6], by which the horizontal divergence term is treated implicitly
(using the most recent horizontal velocity update), and the horizontal
pressure gradient term is treated explicitly. However in the present case the
horizontal pressure gradient term at the final step is evaluated at the same
time level as the (implicit) horizontal divergence term.

The exact, second order time integrals of the variational derivatives are then
given as:

$\displaystyle\langle\boldsymbol{\beta}_{h}^{\parallel},\boldsymbol{\beta}_{h}^{\parallel}\rangle\overline{{\hat{V}}}_{h}=$
$\displaystyle\frac{1}{3}\langle\boldsymbol{\beta}_{h}^{\parallel},\rho_{h}^{n}\boldsymbol{\beta}_{h}^{\parallel}\rangle{\hat{v}_{h}}^{n}+\frac{1}{6}\langle\boldsymbol{\beta}_{h}^{\parallel},\rho_{h}^{n+1}\boldsymbol{\beta}_{h}^{\parallel}\rangle{\hat{v}_{h}}^{n}+\frac{1}{6}\langle\boldsymbol{\beta}_{h}^{\parallel},\rho_{h}^{n}\boldsymbol{\beta}_{h}^{\parallel}\rangle{\hat{v}_{h}}^{\prime}+\frac{1}{3}\langle\boldsymbol{\beta}_{h}^{\parallel},\rho_{h}^{n+1}\boldsymbol{\beta}_{h}^{\parallel}\rangle{\hat{v}_{h}}^{\prime}$
(53a)
$\displaystyle\langle\beta_{h}^{\perp},\beta_{h}^{\perp}\rangle\overline{{\hat{W}}}_{h}=$
$\displaystyle\frac{1}{3}\langle\beta_{h}^{\perp},\rho_{h}^{n}\beta_{h}^{\perp}\rangle{\hat{w}_{h}}^{n}+\frac{1}{6}\langle\beta_{h}^{\perp},\rho_{h}^{n+1}\beta_{h}^{\perp}\rangle{\hat{w}_{h}}^{n}+\frac{1}{6}\langle\beta_{h}^{\perp},\rho_{h}^{n}\beta_{h}^{\perp}\rangle{\hat{w}_{h}}^{n+1}+\frac{1}{3}\langle\beta_{h}^{\perp},\rho_{h}^{n+1}\beta_{h}^{\perp}\rangle{\hat{w}_{h}}^{n+1}$
(53b)
$\displaystyle\langle\gamma_{h},\gamma_{h}\rangle\overline{{\hat{\Phi}}}_{h}=$
$\displaystyle\frac{1}{6}\langle\gamma_{h},\boldsymbol{v}_{h}^{n}\cdot\boldsymbol{\beta}_{h}^{\parallel}\rangle{\hat{v}_{h}}^{n}+\frac{1}{6}\langle\gamma_{h},\boldsymbol{v}^{\prime}_{h}\cdot\boldsymbol{\beta}_{h}^{\parallel}\rangle{\hat{v}_{h}}^{n}+\frac{1}{6}\langle\gamma_{h},\boldsymbol{v}^{\prime}_{h}\cdot\boldsymbol{\beta}_{h}^{\parallel}\rangle{\hat{v}_{h}}^{\prime}+$
$\displaystyle\frac{1}{6}\langle\gamma_{h},w_{h}^{n}\beta_{h}^{\perp}\rangle{\hat{w}_{h}}^{n}+\frac{1}{6}\langle\gamma_{h},w^{n+1}_{h}\beta_{h}^{\perp}\rangle{\hat{w}_{h}}^{n}+\frac{1}{6}\langle\gamma_{h},w^{n+1}_{h}\beta_{h}^{\perp}\rangle{\hat{w}_{h}}^{n+1}+g\langle\gamma_{h},\gamma_{h}\rangle{\hat{z}_{h}}$
(53c)
$\displaystyle\langle\gamma_{h},\gamma_{h}\rangle\overline{{\hat{\Pi}}}_{h}=$
$\displaystyle\frac{1}{2}\langle\gamma_{h},\gamma_{h}\rangle{\hat{\Pi}_{h}}^{n}+\frac{1}{2}\langle\gamma_{h},\gamma_{h}\rangle{\hat{\Pi}_{h}}^{n+1}$
(53d)

Note that while a second order “leap-frog” method has been used to determine
the provisional velocity, ${\hat{v}_{h}}^{\prime}$ in (48), a first order,
forward Euler step is probably also acceptable, since in either case the
resulting mass flux (53a) will be second order in time.

We have not as yet made mention of the spaces used to represent the quantities
found within the operator in (49): $\boldsymbol{q}_{h}$, $\theta_{h}$.
Conservation of energy is satisfied for any representation these quantities,
provided that the skew-symmetry of this operator is preserved. We choose to
represent $\boldsymbol{q}_{h}\in\mathcal{W}_{h}$, such that this may be
derived from $\boldsymbol{u}_{h}$ and $\rho_{h}$ via a weak curl operator [8],
while $\theta_{h}$ is represented in the vertical component of the
$H(\text{div},\Omega)$ space, $\mathcal{U}_{h}^{\perp}$ (such that it is
$C^{0}$ continuous in the vertical direction only). Previous work has shown
that this representation of $\theta_{h}$ yields an improved representation of
the dispersion relation for linearised buoyancy modes [20].

The implicit vertical system in (49) is solved in serial using a direct LU
method, while the projections onto the horizontal
$\mathcal{U}_{h}^{\parallel}$ space in (48), (50), as well as to derive the
horizontal pressure gradients and temperature fluxes in (49) and the
horizontal mass flux in (53a) are solved in parallel using the conjugate
gradient method. These are both done using the PETSc library [22, 21, 23].

###  3.1 Stability analysis

In order to study the stability of the HEVI splitting scheme described above
this is applied to the two dimensional linearised compressible Boussinesq
equations [24], which are given as:

$\displaystyle\frac{\partial u}{\partial t}+\frac{\partial p}{\partial x}$
$\displaystyle=0$ (54a) $\displaystyle\frac{\partial w}{\partial
t}+\frac{\partial p}{\partial z}-b$ $\displaystyle=0$ (54b)
$\displaystyle\frac{\partial p}{\partial t}+c^{2}\Bigg{(}\frac{\partial
u}{\partial x}+\frac{\partial w}{\partial z}\Bigg{)}$ $\displaystyle=0$ (54c)
$\displaystyle\frac{\partial b}{\partial t}+N^{2}w$ $\displaystyle=0$ (54d)

where $u$ and $w$ are the horizontal and vertical velocities, $p$ is the
pressure, $b$ is the buoyancy, $c$ is the speed of sound and $N$ is the Brunt-
Väisälä frequency. Using the splitting scheme detailed above (but with a first
order forward Euler step for the initial provisional velocity, $u^{\prime}$)
and assuming periodic solutions of the form $(u,w,p,b)(t)e^{ikx+ilz}$ these
may be time-stepped as:

$\displaystyle u^{\prime}$ $\displaystyle=u^{n}-\Delta t\ ikp^{n}$ (55a)
$\displaystyle w^{n+1}+\frac{\Delta t\ il}{2}p^{n+1}-\frac{\Delta
t}{2}b^{n+1}$ $\displaystyle=w^{n}-\frac{\Delta t\ il}{2}p^{n}+\frac{\Delta
t}{2}b^{n}$ (55b) $\displaystyle p^{n+1}+\frac{\Delta t\ ilc^{2}}{2}w^{n+1}$
$\displaystyle=p^{n}-\frac{\Delta t\ ilc^{2}}{2}w^{n}-\frac{\Delta t\
ikc^{2}}{2}u^{n}-\frac{\Delta t\ ikc^{2}}{2}u^{\prime}$ (55c) $\displaystyle
b^{n+1}+\frac{\Delta t\ N^{2}}{2}w^{n+1}$ $\displaystyle=b^{n}-\frac{\Delta t\
N^{2}}{2}w^{n}$ (55d) $\displaystyle u^{n+1}+\frac{\Delta t\ ik}{2}p^{n+1}$
$\displaystyle=u^{n}-\frac{\Delta t\ ik}{2}p^{n}$ (55e)

Note that a temporal description of the HEVI scheme that omits the spatial
discretisation is provided in Appendix A. Substituting the expression for
$u^{\prime}$ into that for $p^{n+1}$ and expressing in matrix form gives

$\begin{bmatrix}1&0&\frac{\Delta t\ ik}{2}&0\\\ 0&1&\frac{\Delta t\
il}{2}&-\frac{\Delta t}{2}\\\ 0&\frac{\Delta t\ ilc^{2}}{2}&1&0\\\
0&\frac{\Delta t\ N^{2}}{2}&0&1\\\ \end{bmatrix}\begin{bmatrix}u^{n+1}\\\
w^{n+1}\\\ p^{n+1}\\\ b^{n+1}\\\
\end{bmatrix}=\begin{bmatrix}1&0&-\frac{\Delta t\ ik}{2}&0\\\
0&1&-\frac{\Delta t\ il}{2}&\frac{\Delta t}{2}\\\ -\Delta t\
ikc^{2}&-\frac{\Delta t\ ilc^{2}}{2}&1-\frac{\Delta t^{2}k^{2}c^{2}}{2}&0\\\
0&-\frac{\Delta t\ N^{2}}{2}&0&1\\\ \end{bmatrix}\begin{bmatrix}u^{n}\\\
w^{n}\\\ p^{n}\\\ b^{n}\\\ \end{bmatrix}.$ (56)

In order to analyse the stability of the splitting scheme we inspect the
eigenvalues of the operator

$\boldsymbol{\mathsf{A}}=\begin{bmatrix}1&0&\frac{\Delta t\ ik}{2}&0\\\
0&1&\frac{\Delta t\ il}{2}&-\frac{\Delta t}{2}\\\ 0&\frac{\Delta t\
ilc^{2}}{2}&1&0\\\ 0&\frac{\Delta t\ N^{2}}{2}&0&1\\\
\end{bmatrix}^{-1}\begin{bmatrix}1&0&-\frac{\Delta t\ ik}{2}&0\\\
0&1&-\frac{\Delta t\ il}{2}&\frac{\Delta t}{2}\\\ -\Delta t\
ikc^{2}&-\frac{\Delta t\ ilc^{2}}{2}&1-\frac{\Delta t^{2}k^{2}c^{2}}{2}&0\\\
0&-\frac{\Delta t\ N^{2}}{2}&0&1\\\ \end{bmatrix}$ (57)

For a given eigenvalue, $\lambda_{i}$, this will be unstable if
$\lambda_{i}\lambda_{i}^{*}>1$. The four eigenvalues correspond to two
acoustic modes, which differ in the sign of their imaginary component only,
and two gravity (buoyancy) modes which similarly differ only in the sign of
the imaginary term.

The results for the new HEVI scheme are compared to the standard second order
trapezoidal HEVI scheme, for which the vertical terms are all treated
implicitly, and the horizontal terms explicitly. This scheme may be written
as:

$\begin{bmatrix}1&0&0&0\\\ 0&1&\frac{\Delta t\ il}{2}&-\frac{\Delta t}{2}\\\
0&\frac{\Delta t\ ilc^{2}}{2}&1&0\\\ 0&\frac{\Delta t\ N^{2}}{2}&0&1\\\
\end{bmatrix}\begin{bmatrix}u^{\prime}\\\ w^{\prime}\\\ p^{\prime}\\\
b^{\prime}\\\ \end{bmatrix}=\begin{bmatrix}1&0&-\Delta t\ ik&0\\\
0&1&-\frac{\Delta t\ il}{2}&\frac{\Delta t}{2}\\\ -\Delta t\
ikc^{2}&-\frac{\Delta t\ ilc^{2}}{2}&1&0\\\ 0&\frac{\Delta t\ N^{2}}{2}&0&1\\\
\end{bmatrix}\begin{bmatrix}u^{n}\\\ w^{n}\\\ p^{n}\\\ b^{n}\\\ \end{bmatrix}$
(58)

$\begin{bmatrix}u^{n+1}\\\ w^{n+1}\\\ p^{n+1}\\\ b^{n+1}\\\
\end{bmatrix}=\begin{bmatrix}1&0&-\frac{\Delta t\ ik}{2}&0\\\
0&1&-\frac{\Delta t\ il}{2}&\frac{\Delta t}{2}\\\ -\frac{\Delta t\
ikc^{2}}{2}&-\frac{\Delta t\ ilc^{2}}{2}&1&0\\\ 0&-\frac{\Delta t\
N^{2}}{2}&0&1\\\ \end{bmatrix}\begin{bmatrix}u^{n}\\\ w^{n}\\\ p^{n}\\\
b^{n}\\\ \end{bmatrix}+\begin{bmatrix}0&0&-\frac{\Delta t\ ik}{2}&0\\\
0&0&-\frac{\Delta t\ il}{2}&\frac{\Delta t}{2}\\\ -\frac{\Delta t\
ikc^{2}}{2}&-\frac{\Delta t\ ilc^{2}}{2}&0&0\\\ 0&-\frac{\Delta t\
N^{2}}{2}&0&0\\\ \end{bmatrix}\begin{bmatrix}u^{\prime}\\\ w^{\prime}\\\
p^{\prime}\\\ b^{\prime}\\\ \end{bmatrix}$ (59)

Combining the above expressions gives a linear operator by which the solution
at time level $n+1$ may be determined from the solution at time level $n$,
from which the corresponding eigenvalues may be derived.

For completeness we also compare against the Crank-Nicolson scheme, which is
semi-implicit in both horizontal and vertical dimensions, and neutrally stable
for all vertical and horizontal wave numbers for both buoyancy and acoustic
modes. The amplification factors for the Crank-Nicolson scheme are given as:

$\boldsymbol{\mathsf{A}}=\begin{bmatrix}1&0&\frac{\Delta t\ ik}{2}&0\\\
0&1&\frac{\Delta t\ il}{2}&-\frac{\Delta t}{2}\\\ \frac{\Delta t\
ikc^{2}}{2}&\frac{\Delta t\ ilc^{2}}{2}&1&0\\\ 0&\frac{\Delta t\
N^{2}}{2}&0&1\\\ \end{bmatrix}^{-1}\begin{bmatrix}1&0&-\frac{\Delta t\
ik}{2}&0\\\ 0&1&-\frac{\Delta t\ il}{2}&\frac{\Delta t}{2}\\\ -\frac{\Delta t\
ikc^{2}}{2}&-\frac{\Delta t\ ilc^{2}}{2}&1&0\\\ 0&-\frac{\Delta t\
N^{2}}{2}&0&1\\\ \end{bmatrix}$ (60)

Figure 1: Amplification factor for the acoustic waves using the Crank-Nicolson
(left), new (center) and trapezoidal HEVI (right) splitting. The black line on
the center plot indicates an amplification factor of $1.0+10^{-12}$,
indicating the onset of unstable modes. The Crank-Nicolson scheme is neutrally
stable for all modes, the new scheme is neutrally stable for all modes below a
finite horizontal wave number and the trapezoidal scheme is unstable for all
modes.

The amplification factors for the acoustic and gravity modes for the three
different time splitting schemes are given in figs. 1 and 2 respectively for a
periodic domain of size 1000.0m in both dimensions, with
$c=340.0\mathrm{ms}^{-1}$, $N=0.01\mathrm{s}^{-1}$. While the Crank-Nicolson
scheme is neutrally stable, $|\lambda|=1$ for both modes, and the trapezoidal
HEVI scheme is unstable, $|\lambda|>1$ for both modes (albeit with a very
small amplification factor for the gravity modes), the energetically balanced
splitting scheme is neutrally stable for all acoustic modes below some
horizontal wave number, and unconditionally stable for all gravity modes. This
represents a qualitative improvement on the standard second order HEVI
splitting.

Figure 2: Amplification factor for the gravity waves using the Crank-Nicolson
(left), new (center) and trapezoidal HEVI (right) splitting. The Crank-
Nicolson scheme and the new scheme are both neutrally stable for all modes,
while the trapezoidal scheme is unstable for all modes.

##  4 Potential temperature upwinding

This section describes the suppression of spurious oscillations in potential
temperature via upwinding of the potential temperature diagnostic equation in
an energetically consistent manner. Since the potential vorticity and
potential temperature are represented within the skew-symmetric operator for
the discrete variational form of the compressible Euler equations (46), these
quantities may be upwinded so as to suppress oscillations without breaking
energy conservation (so long as skew-symmetry is preserved). This strategy has
been used previously to suppress oscillations for the shallow water equations
by upwinding potential vorticity so as to dissipate potential enstrophy [25,
26, 27], and for the compressible Euler equations by upwinding potential
temperature so as to generate entropy [28]. In the present context the
upwinding of potential temperature is applied not to the material form
advection equation, but rather to the diagnostic equation by which the
potential temperature is derived from the prognostic variables of density (1b)
and density-weighted potential temperature (1c).

While this is typically achieved for variational methods by augmenting the
test space with a linear correction which allows for dissipation in the
direction of the flow [29], in the present case this is achieved by
integrating the test functions backwards along velocity characteristics [27].
Here the upwinding is applied in the vertical direction only. Since the
discrete potential temperature is represented as
$\theta_{h}\in\mathcal{U}_{h}^{\perp}$ [8], for which the vertical component
of the tensor product basis is represented by piecewise linear polynomials,
for this case the upwinded bases are effectively equivalent to the SUPG method
[29, 28]. However in the present case this upwinding is applied to the
diagnostic equation (the discrete form of the relation $\theta:=\Theta/\rho$,
(62)), and not the material form of the advection equation,
$\partial\theta/\partial t+\boldsymbol{u}\cdot\nabla\theta=0$, as is customary
in the SUPG method. Since $C^{0}$ continuity for the $\mathcal{U}_{h}^{\perp}$
space is enforced only in the vertical direction, upwinding of $\theta_{h}$ in
the horizontal directions will have no effect.

The upwinded trial functions, $\beta_{h}^{\perp,u}\in\mathcal{U}_{h}^{\perp}$
are computed by evaluating the quadrature points at vertically _downwind_
locations, $\zeta^{d}$, which are computed to first order as [27]:

$\zeta^{d}=\zeta+0.5\Delta
t\sum_{i}\hat{w}_{i}\beta_{i}^{\perp}(\xi,\eta,\zeta),$ (61)

where $\xi$, $\eta$ and $\zeta$ are the local element coordinates. Note that
in the above expression the velocity is interpolated to its local element
value, and not its global physical value. Note also that the choice to
integrate the vertical quadrature points downwind for a time of $0.5\Delta t$
is arbitrary, and may be tuned to optimise the balance between reduced
oscillation and less diffusive solutions. The upwinded trial functions are
then evaluated as $\beta_{h}^{\perp,u}=\beta_{h}^{\perp}(\xi,\eta,\zeta^{d})$.
The upwinded potential temperature is then diagnosed as:

$\langle\beta_{h}^{\perp,u},\rho_{h}\beta_{h}^{\perp}\rangle\hat{\theta}_{h}=\langle\beta_{h}^{\perp,u},\gamma_{h}\rangle\hat{\Theta}_{h}.$
(62)

Note that in the above expression only the test functions for the potential
temperature diagnostic equation are evaluated at upstream locations, while the
trial functions remain static. Moreover since continuity is enforced only in
the vertical direction for the trial space of the potential temperature,
$\mathcal{U}_{h}^{\perp}$, it is only in this direction that upwinding has any
material effect.

While the variational form of the compressible Euler equations (4) conserve
entropy, $s=c_{p}\log(\theta)+s_{0}$ (see Appendix B), this is not true in the
discrete form, since transcendental functions such as logarithms cannot be
integrated exactly in the discrete form. As such, while the upwinding of the
test functions as described above will generate additional entropy, there are
already additional entropy sources at the discrete level.

##  5 Results

###  5.1 Baroclinic instability

In order to verify the energetically balanced HEVI scheme, this is applied to
the solution of a standard test case for baroclinic instability on the sphere
[30]. This test involves a small perturbation to an otherwise hydrostatically
and geostrophically balanced atmosphere, which over several days generates a
baroclinic wave. The model is stabilised using a biharmonic viscosity term on
both the horizontal potential temperature equation and momentum equation
(which does not conserve energy due to the explicit time stepping), and a
Rayleigh friction term in the top layers of the atmosphere in the vertical
(which is not strictly necessary however helps to accelerate the initial
hydrostatic adjustment process) [8]. Figures 3, 4 and 5 show the surface level
Exner pressure, as well as the vertical vorticity component and potential
temperature at a height of $z=1.57\mathrm{km}$ at days 8 and 10 respectively.
As an initial validation exercise these compare well against previously
published results [30, 9], however they do not allow us to draw any particular
conclusion as to the superiority of any particular scheme. The signature of
the baroclinic instability is also observed in the perturbation to the steady
hydrostatic pressure profile, for which a vertical cross section at
$50^{\circ}\mathrm{N}$ is shown at days 8 and 10 in fig 6.

Figure 3: Baroclinic instability test case: Bottom level Exner pressure, days
8 and 10.

Figure 4: Baroclinic instability test case: Vertical vorticity component at
$z=1.57km$, days 8 and 10.

Figure 5: Baroclinic instability test case: Potential temperature at
$z=1.57km$, days 8 and 10.

Figure 6: Baroclinic instability test case: Zonal-vertical pressure slice at
$50^{\circ}N$, days 8 and 10.

####  5.1.1 Comparison to TRAP(2,3,2)

The energetically balanced HEVI scheme is compared to the horizontally third
order, vertically second order trapezoidal scheme, TRAP(2,3,2) [6, 7]. In most
respects the energetics of the two schemes are similar, despite the fact that
the TRAP(2,3,2) scheme is not energetically balanced, but is higher order in
the horizontal (also a time step of $\Delta t=60\mathrm{s}$ is used for the
energetically balanced scheme, while the TRAP(2,3,2) simulation was run with a
time step of $\Delta t=120\mathrm{s}$). While the TRAP(2,3,2) scheme is stable
for longer time steps than the energetically balanced scheme, it also involves
an additional substep in both the vertical and horizontal directions. Both
schemes are stabilised using biharmonic viscosity in the horizontal and
Rayleigh friction in the top layers in the vertical [8, 9].

The one respect in which the two schemes differ markedly is the potential to
kinetic power exchanges, the transition of potential to kinetic energy over
time. This is computed (assuming homogeneous Dirichlet boundary conditions for
the vertical velocity) as the discrete equivalent of the term $-\partial
P/\partial t=g\int z\partial(\rho w)/\partial z\mathrm{d}\Omega$ as $-\partial
P_{h}/\partial
t=g\hat{z}_{h}\langle\gamma_{h},\gamma_{h}\rangle\boldsymbol{\mathsf{E}}^{3,2}_{\perp}\overline{\hat{W}}_{h}$
[8]. For the potential to kinetic power exchange the energetically balanced
scheme exhibits a higher frequency oscillation not present in the TRAP(2,3,2)
scheme, as shown in fig. 7. Fourier analysis of these power exchanges, taken
between days 6 and 10, shows that while the amplitudes of the low-frequency
modes ($k<10$) match each other closely, both solutions exhibit a secondary
peak, which for the TRAP(2,3,2) scheme occurs at $k=33$ (2.91 hours), while
for the energetically balanced scheme this peak occurs at $k=131$ (0.73
hours). Curiously, the secondary peak for the energetically balanced scheme is
almost exactly four times the frequency observed for the TRAP(2,3,2) scheme.
Since the secondary peak for the TRAP(2,3,2) scheme is closer in frequency to
the low frequency modes, which exhibit very similar behaviour for both
schemes, there is a possibility that this may alias onto the low frequency
time scales of the baroclinic motions.

Figure 7: Baroclinic instability test case: Comparison of potential to kinetic
energy power exchanges for the energetically balanced and TRAP(2,3,2) HEVI
schemes; left: real space representation from day 9 to 10, right: amplitude of
the corresponding Fourier modes from days 6 to 10.

###  5.2 3D Rising bubble

In order to verify the scheme for non-hydrostatic dynamics this is applied to
a standard test case for a 3D rising bubble, which is initialised as a small
perturbation in potential temperature against a constant background value,
with the Exner pressure specified so as to satisfy hydrostatic balance [2, 31,
32]. The domain is configured as a horizontally periodic box of size
$1000\mathrm{m}\times 1000\mathrm{m}\times 1500\mathrm{m}$ with a flat bottom
and top with homogeneous boundary conditions, and discretised using $24\times
24$ elements of degree $p=3$ in the horizontal and 150 vertical levels, with a
time step of $\Delta t=0.01\mathrm{s}$. No dissipation of any kind is applied
in the vertical, or to the potential temperature advection equation. However
the energy conserving upwinding is applied to the potential temperature
diagnostic equation (62) in order to suppress the development of oscillations
associated with non-hydrostatic motions. This is not applied to the baroclinic
instability test case, which evolves in a predominantly hydrostatic regime.

Since the implicit vertical time stepping scheme conserves energy [9], and the
energetic exchanges are exactly balanced in the horizontal, no stabilisation
is required of the energetically balanced HEVI scheme, provided that the
horizontal CFL condition is satisfied. This is in contrast to the TRAP(2,3,2)
scheme, which requires a horizontal biharminic viscosity in order suppress the
onset of numerical instability. In order to directly compare to the
TRAP(2,3,2) scheme, simulations of the energetically balanced HEVI scheme have
been run both with and without horizontal biharmonic viscosity on the momentum
equation [11], with a coefficient of $624.78\mathrm{m}^{4}\mathrm{s}^{-1}$.

Figure 8: 3D rising bubble: Potential temperature cross section at $y=0$ at
times 200s (left) and 400s (right).

While the results qualitatively match those previously published in terms of
both the position and deformation of the bubble [31, 32], in the absence of
potential temperature upwinding a secondary oscillation develops in the wake
of the bubble as it steepens, as shown in fig. 8. This oscillation is also
observed for the TRAP (2,3,2) scheme [9]. The evolution of the energy and
power exchanges (computed as in [8]) are shown in fig 9 , and clearly show the
signature of the high frequency internal gravity wave, as well as the increase
in vertical kinetic energy and decrease in potential energy which represent
the ascension of the bubble.

Figure 9: 3D rising bubble: Difference in energy from initial value (left),
and power exchanges associated with the changes in potential to kinetic and
internal to kinetic energy with time (right). The power exchanges are computed
as the discrete analogues of the relations $g\int z\partial(\rho w)/\partial
z\mathrm{d}\Omega$ and
$\int\nabla\cdot(\rho\boldsymbol{u}\theta)\Pi\mathrm{d}\Omega$ respectively.

####  5.2.1 Upwinding of potential temperature

In order to suppress the spurious potential temperature oscillations observed
in fig. 8 in an energetically consistent manner, the upwinding formulation
presented in section 4 is applied to the rising bubble test case. As observed
in fig. 10, this upwinding formulation effectively clears up the spurious
overshoots and oscillations. Since the upwinding is performed within the skew-
symmetric operator, as detailed in (46), the energy conservation errors are
almost identical to those for the original formulation, as shown in fig. 11.

Fig. 11 also shows the energy conservation errors for the TRAP(2,3,2) scheme
(with biharmonic viscosity) and the energetically balanced scheme in the
absence of viscosity. While horizontal biharmonic viscosity is necessary to
stabilise the TRAP(2,3,2) scheme, the energetically balanced scheme may be run
entirely without dissipation. As observed, the difference in the energy
conservation error of the energetically balanced scheme with and without
biharmonic viscosity is small compared to the overall energy conservation
error, indicating that most of the energy conservation error arises from the
explicit time stepping of the horizontal velocity, and not the dissipation
term. In all cases the total energy of the system decreases with time.

This figure also shows the power exchanges for the TRAP(2,3,2) scheme. These
are almost identical to those for the energetically balanced scheme shown in
fig. 9, such that the difference in the time scale of the internal gravity
mode observed for the baroclinic test case is not observed for the high
resolution rising bubble test case, for which the ratio of the horizontal and
vertical resolutions is closer to unity, and the vertical time scales are
effectively resolved by the time step.

Figure 10: 3D rising bubble: Upwinded potential temperature cross section at
$y=0$ at times 200s (left) and 400s (right).

Figure 11: 3D rising bubble: Comparison of energy conservation errors for the
original and upwinded potential temperature formulations, as well as the
upwinded inviscid solution and the TRAP(2,3,2) scheme (left). Power exchanges
associated with the changes in potential to kinetic and internal to kinetic
energy with time for the TRAP (2,3,2) scheme (right).

##  6 Conclusions

A new horizontally explicit/vertically implicit time splitting scheme for non-
hydrostatic atmospheric dynamics is introduced. The scheme allows for the
exact balance of energetic exchanges when coupled with a spatial and vertical
time integration scheme that preserve the skew-symmetric property of the
Hamiltonian formulation of the equations of motion. However since the choice
of horizontal velocity for which exchanges balance differs from the horizontal
velocity at the end of the time step, the scheme still permits an energy
conservation error. Linear eigenvalue analysis shows that the second order
scheme improves upon the existing second order trapezoidal HEVI splitting by
ensuring the neutral stability of all buoyancy modes, and all acoustic modes
below a certain horizontal CFL number. Comparison to a horizontally third
order, vertically second order trapezoidal splitting for the full 3D
compressible Euler equations on the sphere shows that both schemes exhibit a
secondary oscillation in the potential to kinetic power exchanges, associated
with an internal gravity mode, however this oscillation has a temporal
frequency approximately four times faster for the new energetically balances
scheme, such that the time scales of the oscillation for the trapezoidal
scheme are closer to those of the baroclinic motions.

An energetically consistent formulation for potential temperature upwinding is
also presented, by which the test functions are evaluated at vertically
downstream locations. This ensures that spurious oscillations in potential
temperature can be removed for vertical motions at non-hydrostatic regimes.
This upwinding is applied within the skew-symmetric formulation, so as to
suppress oscillations associated with vertical motions without altering the
energetic properties of the model, for which the exchanges remain balanced.

Notably, for the rising bubble test case, which is run on an affine geometry,
the energetically balanced HEVI scheme may be run stably without dissipation
of any kind. For the baroclinic instability test case on the sphere however,
horizontal biharmonic viscosity is still required in order to stabilise the
simulation. Since the model uses inexact integration, and the sphere is
represented as a non-affine geometry, it is possible that aliasing errors may
provide an additional source of instability. As part of future work exact
integration will be employed in order to determine if the scheme can be run
entirely without dissipation for the baroclinic instability test case on the
sphere also.

The code used in this article is publicly available at:
`https://github.com/davelee2804/MiMSEM`.

## Appendix A: Temporal discreisation of the energetically balanced HEVI
scheme

For the benefit of readers who are not familiar with compatible finite element
methods, in this appendix we present a formulation of the energetically
balanced HEVI scheme that is discretised in the temporal dimension only. Note
that in order to preserve energetic balance a spatial discretisation that
respects the skew-symmetric structure of the equations of motion must also be
employed. The temporal scheme is expressed as follows:

_Step 1: Explicit horizontal momentum solve_

$\boldsymbol{v}^{\prime}=\boldsymbol{v}^{n}-\Delta t\
{q}_{\perp}^{n}\times\boldsymbol{V}^{n}+\Delta t\
\boldsymbol{q}_{\parallel}^{n}W^{n}-\Delta t\ \nabla_{h}\Phi^{n}-\Delta t\
\theta^{n}\nabla_{h}\Pi^{n}$ (63)

_Step 2: Implicit vertical solve (including horizontal divergence)_

$\displaystyle w^{n+1}$ $\displaystyle=w^{n}-\Delta t\
\partial_{z}\overline{\Phi}-\Delta t\
\theta^{n+1/2}\partial_{z}\overline{\Pi}-\Delta t\
\boldsymbol{q}_{\parallel}^{n+1/2}\cdot\overline{\boldsymbol{V}}$ (64a)
$\displaystyle\rho^{n+1}$ $\displaystyle=\rho^{n}-\Delta t\
\partial_{z}\overline{W}-\Delta t\ \nabla_{h}\cdot\overline{\boldsymbol{V}}$
(64b) $\displaystyle\Theta^{n+1}$ $\displaystyle=\rho^{n}-\Delta t\
\partial_{z}(\theta^{n+1/2}\overline{W})-\Delta t\
\nabla_{h}\cdot(\theta^{n+1/2}\overline{\boldsymbol{V}})$ (64c)

with

$\displaystyle\overline{V}$
$\displaystyle=\frac{1}{6}\Big{(}2\rho^{n}\boldsymbol{v}^{n}+\rho^{n+1}\boldsymbol{v}^{n}+\rho^{n}\boldsymbol{v}^{\prime}+2\rho^{n+1}\boldsymbol{v}^{\prime}\Big{)}$
(65a) $\displaystyle\overline{W}$
$\displaystyle=\frac{1}{6}\Big{(}2\rho^{n}w^{n}+\rho^{n+1}w^{n}+\rho^{n}w^{n+1}+2\rho^{n+1}w^{n+1}\Big{)}$
(65b) $\displaystyle\overline{\Phi}$
$\displaystyle=\frac{1}{6}\Big{(}\boldsymbol{v}^{n}\boldsymbol{v}^{n}+\boldsymbol{v}^{n}\boldsymbol{v}^{\prime}+\boldsymbol{v}^{\prime}\boldsymbol{v}^{\prime}+w^{n}w^{n}+w^{n}w^{n+1}+w^{n+1}w^{n+1}\Big{)}+gz$
(65c) $\displaystyle\overline{\Pi}$
$\displaystyle=\frac{1}{2}\Big{(}\Pi^{n}+\Pi^{n+1}\Big{)}$ (65d)

_Step 3: Explicit horizontal momentum solve_

$\boldsymbol{v}^{n+1}=\boldsymbol{v}^{n}-\Delta t\
{q}_{\perp}^{n+1/2}\times\overline{\boldsymbol{V}}+\Delta t\
\boldsymbol{q}_{\parallel}^{n}\overline{W}-\Delta t\
\nabla_{h}\overline{\Phi}-\Delta t\ \theta^{n+1/2}\nabla_{h}\overline{\Pi}$
(66)

where $\nabla_{h}$ is the horizontal gradient operator.

Note that if energetic balance is not of concern, a standard time centered
construction of the mass flux terms, $\boldsymbol{V}$, $W$ and Bernoulli
function, $\Phi$ at $n+1/2$ should also preserve second order accuracy. Also
note that a simple first order Euler integration is used here in step 1, which
is sufficient to construct a second order mass flux and Bernoulli function.

## Appendix B: Entropy conservation for the variational form of the
compressible Euler equations

The time derivative of the internal energy is derived by multiplying the
density-weighted potential temperature advection equation,

$\frac{\partial\rho\theta}{\partial
t}+\nabla\cdot(\rho\boldsymbol{u}\theta)=0,$ (67)

by the Exner pressure,

$\frac{\delta
I}{\delta\Theta}=\Pi=c_{p}\Bigg{(}\frac{R\rho\theta}{p_{0}}\Bigg{)}^{R/c_{v}},$
(68)

such that the internal energy, ${I}$ evolves as

$\frac{\partial{I}}{\partial t}=\Bigg{\langle}\frac{\delta
I}{\delta\Theta},\frac{\partial\Theta}{\partial
t}\Bigg{\rangle}=\Bigg{\langle}\Pi,\frac{\partial\rho\theta}{\partial
t}\Bigg{\rangle}=\frac{c_{v}}{c_{p}}\frac{\partial}{\partial
t}\int\rho\theta\Pi\mathrm{d}\Omega=-\langle\Pi,\nabla\cdot(\boldsymbol{u}\rho\theta)\rangle,$
(69)

since $R=c_{p}-c_{v}$ and $\Theta=\rho\theta$. This expression may then be
expanded as:

$\displaystyle\frac{c_{v}}{c_{p}}\frac{\partial}{\partial
t}\int\rho\theta\Pi\mathrm{d}\Omega$
$\displaystyle=-\langle\Pi,\theta\nabla\cdot(\rho\boldsymbol{u})+\rho\boldsymbol{u}\cdot\nabla\theta\rangle$
(70a)
$\displaystyle=-\Bigg{\langle}\rho\theta\Pi,\frac{1}{\rho}\nabla\cdot(\rho\boldsymbol{u})+\frac{\boldsymbol{u}}{\theta}\cdot\nabla\theta\Bigg{\rangle}$
(70b)
$\displaystyle=\Bigg{\langle}\rho\theta\Pi,\frac{\partial\log(\rho)}{\partial
t}-\boldsymbol{u}\cdot\nabla\log(\theta)\Bigg{\rangle}$ (70c)

The left hand side of this expression may be re-arranged as:

$\frac{c_{v}}{c_{p}}\int\frac{\partial\rho\theta\Pi}{\partial
t}\mathrm{d}\Omega=\frac{c_{v}}{c_{p}}\Bigg{\langle}\frac{\rho\theta\Pi}{\rho\theta\Pi},\frac{\partial\rho\theta\Pi}{\partial
t}\Bigg{\rangle}=\frac{c_{v}}{c_{p}}\Bigg{\langle}\rho\theta\Pi,\frac{\partial\log(\rho\theta\Pi)}{\partial
t}\Bigg{\rangle}=\Bigg{\langle}\rho\theta\Pi,\frac{\partial\log(\rho)}{\partial
t}-\boldsymbol{u}\cdot\nabla\log(\theta)\Bigg{\rangle}$ (71)

For $\rho\theta\Pi\in L^{2}(\Omega)$, cancellation of this factor from both
sides yields the variational expression

$\frac{c_{v}}{c_{p}}\Bigg{\langle}\gamma,\frac{\partial\log(\rho\theta\Pi)}{\partial
t}\Bigg{\rangle}=\Bigg{\langle}\gamma,\frac{\partial\log(\rho)}{\partial
t}-\boldsymbol{u}\cdot\nabla\log(\theta)\Bigg{\rangle}\qquad\forall\gamma\in
L^{2}(\Omega)$ (72)

Recalling the equation of state gives

$\rho\theta\Pi=c_{p}\Bigg{(}\frac{R}{p_{0}}\Bigg{)}^{R/c_{v}}(\rho\theta)^{c_{p}/c_{v}},$
(73) $\Bigg{\langle}\gamma,\frac{\partial\log(\rho\theta)}{\partial
t}\Bigg{\rangle}=\Bigg{\langle}\gamma,\frac{\partial\log(\rho)}{\partial
t}-\boldsymbol{u}\cdot\nabla\log(\theta)\Bigg{\rangle}\qquad\forall\gamma\in
L^{2}(\Omega)$ (74)

which is further simplified as:

$\Bigg{\langle}\gamma,\frac{\partial\log(\theta)}{\partial
t}+\boldsymbol{u}\cdot\nabla\log(\theta)\Bigg{\rangle}=0\qquad\forall\gamma\in
L^{2}(\Omega)$ (75)

Recalling the original definition of the entropy then gives

$\Bigg{\langle}\gamma,\frac{\partial s}{\partial t}+\boldsymbol{u}\cdot\nabla
s\Bigg{\rangle}=0\qquad\forall\gamma\in L^{2}(\Omega)$ (76)

such that entropy is materially conserved. Setting the test function as the
density $\gamma=\rho$ in the above expression, and multiplying the continuity
equation by the entropy, we have

$\displaystyle\Bigg{\langle}\rho,\frac{\partial s}{\partial
t}+\boldsymbol{u}\cdot\nabla s\Bigg{\rangle}$ $\displaystyle=0$ (77a)
$\displaystyle\Bigg{\langle}s,\frac{\partial\rho}{\partial
t}+\nabla\cdot\boldsymbol{u}\rho\Bigg{\rangle}$ $\displaystyle=0.$ (77b)

Adding this to the continuity equation, and recalling that
$\boldsymbol{U}:=\rho\boldsymbol{u}$ gives

$\int\frac{\partial\rho s}{\partial t}+\nabla\cdot(\boldsymbol{u}\rho
s)\mathrm{d}\Omega=0,$ (78)

which yields an additional conservation law for the entropy function-entropy
flux pair $\rho s,\boldsymbol{u}\rho s$ [33].

##  7 Acknowledgments

David Lee would like to thank Dr Justin Freeman for his continued
encouragement, and Drs. Marcus Thatcher and John McGregor for their support
and access to computing resources. This project was supported by resources and
expertise provided by CSIRO IMT Scientific Computing. We are also grateful to
the two anonymous reviewers, whose comments and insights helped to improve the
quality of this article.

## References

  * [1] P. Ullrich, C. Jablonowski, Operator-split Runge–Kutta–Rosenbrock methods for nonhydrostatic atmospheric models, Mon. Wea. Rev. 140 (2012) 1257–1284.
  * [2] F. X. Giraldo, J. F. Kelly, E. M. Constantinescu, Implicit-explicit formulations of a three-dimensional nonhydrostatic unified model of the atmosphere (NUMA), SIAM J. Sci. Comput. 35 (2013) B1162–B1194.
  * [3] L. Bao, R. Klöfkorn, R. D. Nair, Horizontally explicit and vertically implicit (HEVI) time discretization scheme for a discontinuous Galerkin nonhydrostatic model, Mon. Wea. Rev. 143 (2015) 972–990.
  * [4] D. J. Gardner, J. E. Guerra, F. P. Hamon, D. R. Reynolds, P. A. Ullrich, C. S. Woodward, Implicit-explicit (IMEX) Runge-Kutta methods for non-hydrostatic atmospheric models, Geosci. Model Dev. 11 (2018) 1497–1515.
  * [5] A. Steyer, C. J. Vogl, M. Taylor, O. Guba, Efficient IMEX Runge-Kutta Methods for Nonhydrostatic dynamics, arXiv:1906.07219 (2019)
  * [6] H. Weller, S. J. Lock, N. Wood, Runge–Kutta IMEX schemes for the Horizontally Explicit/Vertically Implicit (HEVI) solution of the wave equations, J. Comp. Phys. 252 (2013) 365–381.
  * [7] S.-J. Lock, N. Wood, H. Weller, Numerical analyses of Runge–Kutta implicit–explicit schemes for horizontally explicit, vertically implicit solutions of atmospheric models, Q. J. R. Meteorol. Soc. 140 (2014) 1654–1669.
  * [8] D. Lee, A. Palha, A mixed mimetic spectral element model of the 3D compressible Euler equations on the cubed sphere, J. Comp. Phys. 401 (2020) 108993.
  * [9] D. Lee, An energetically balanced, quasi-Newton integrator for non-hydrostatic vertical atmospheric dynamics, J. Comp. Phys. (2021) 109988.
  * [10] D. Lee, A. Palha, M. Gerritsma, Discrete conservation properties for shallow water flows using mixed mimetic spectral elements, J. Comp. Phys. 357 (2018) 282–304.
  * [11] D. Lee, A. Palha, A mixed mimetic spectral element model of the rotating shallow water equations on the cubed sphere, J. Comp. Phys. 375 (2018) 240–262.
  * [12] M. Gerritsma, Edge functions for spectral element methods, in: Spectral and high order methods for partial differential equations, Lecture Notes in Computational Science and Engineering, Springer 76 (2011) 199–207.
  * [13] J. Kreeft, M. Gerritsma, Mixed mimetic spectral element method for Stokes flow: A pointwise divergence-free solution, J. Comp. Phys. 240 (2013) 284–309.
  * [14] A. Palha, P. P. Rebelo, R. Hiemstra, J. Kreeft, M. Gerritsma, Physics-compatible discretization techniques on single and dual grids, with application to the Poisson equation of volume forms, J. Comp. Phys. 257 (2014) 1394–1422.
  * [15] R. Hiemstra, D. Toshniwal, R. H. M Huijsmans, M. I. Gerritsma, High order geometric methods with exact conservation properties, J. Comp. Phys. 257 (2014) 1444–1471.
  * [16] V. Jain, Y. Zhang, A. Palha, M. Gerritsma, Construction and application of algebraic dual polynomial representations for finite element methods on quadrilateral and hexahedral meshes, Comput. Math. Appl. in press (2020).
  * [17] M. Gerritsma, P. Bochev, A spectral mimetic least-squares method for the Stokes equations with no-slip boundary conditions, Comput. Math. Appl. 71 (2016) 2285–2300.
  * [18] E. Hairer, Energy-Preserving Variant of Collocation Methods, JNAIAM J. Numer. Anal. Ind. Appl. Math. 5 (2010) 7–84.
  * [19] D. Cohen, E. Hairer, Linear energy-preserving integrators for Poisson systems, BIT Numer. Math. 51 (1) (2011) 91–101.
  * [20] T. Melvin, T. Benacchio, J. Thuburn, C. Cotter, Choice of function spaces for thermodynamic variables in mixed finite‐element methods, Q. J. Royal. Meteorol. Soc. 144 (2018) 900–916.
  * [21] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, D. A. May, L. C. McInnes, K. Rupp, B. F. Smith, S. Zampini, H. Zhang, PETSc Web page, http://www.mcs.anl.gov/petsc (2017).
  * [22] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, D. A. May, L. C. McInnes, K. Rupp, P. Sanan, B. F. Smith, S. Zampini, H. Zhang, PETSc users manual, Tech. Rep. ANL-95/11 - Revision 3.8, Argonne National Laboratory (2017).
  * [23] S. Balay, W. D. Gropp, L. C. McInnes, B. F. Smith, Efficient management of parallelism in object oriented numerical software libraries, in: E. Arge, A. M. Bruaset, H. P. Langtangen (Eds.), Modern Software Tools in Scientific Computing, Birkhäuser Press, 1997, pp. 163–202.
  * [24] D. Durran, Numerical Methods for Fluid Dynamics: With Applications in Geophysics, Springer (2010) 516pp.
  * [25] R. Sadourny, C. Basdevant, Parameterization of Subgrid Scale Barotropic and Baroclinic Eddies in Quasi-geostrophic Models: Anticipated Potential Vorticity Method, J. Atmos. Sci. 42 (1985) 1353–1363
  * [26] A. T. T. McRae, C. J. Cotter, Energy- and enstrophy conserving schemes for the shallow water equations, based on mimetic finite elements, Q. J. R. Meteorol. Soc. 140 (2014) 2223–2234.
  * [27] D. Lee, Petrov-Galerkin flux upwinding for mixed mimetic spectral elements, and its application to geophysical flow problems, Comput. Math. Appl. 89 (2021) 68–77
  * [28] G. A. Wimmer, C.J. Cotter, W. Bauer, Energy conserving SUPG methods for compatible finite element schemes in numerical weather prediction, arXiv:2001.09590 (2020)
  * [29] A. N. Brooks, T. J. R. Hughes, Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations, Comput. Meth. Appl. Mech. Engrg. 32 (1982) 199–259
  * [30] P. A. Ullrich, T. Melvin, C. Jablonowski, A. Staniforth, A proposed baroclinic wave test case for deep‐ and shallow‐atmosphere dynamical cores, Q. J. R. Meteorol. Soc. 140 (2014) 1590–1602.
  * [31] D. S. Abdi, F. X. Giraldo, Efficient construction of unified continuous and discontinuous Galerkin formulations for the 3D Euler equations, J. Comp. Phys. 320 (2016) 46–68.
  * [32] T. Melvin, T. Benacchio, B. Shipway, N. Wood, J. Thuburn, C. Cotter, A mixed finite-element, finite-volume, semi-implicit discretisation for atmospheric dynamics: Cartesian geometry, Q. J. R. Meteorol. Soc. (2019) 1–19.
  * [33] E. Tadmor, Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems, Acta Numerica 12 (2003) 451–512.