# Entropy and energy conservation for thermal atmospheric dynamics using mixed
compatible finite elements

Kieran Ricardo  David Lee<EMAIL_ADDRESS>Kenneth Duru  Mathematical
Sciences Institute, Australian National University, Canberra, Australia
Bureau of Meteorology, Melbourne, Australia

###### Abstract

Atmospheric systems incorporating thermal dynamics must be stable with respect
to both energy and entropy. While energy conservation can be enforced via the
preservation of the skew-symmetric structure of the Hamiltonian form of the
equations of motion, entropy conservation is typically derived as an
additional invariant of the Hamiltonian system, and satisfied via the exact
preservation of the chain rule. This is particularly challenging since the
function spaces used to represent the thermodynamic variables in compatible
finite element discretisations are typically discontinuous at element
boundaries. In the present work we negate this problem by constructing our
equations of motion via weighted averages of skew-symmetric formulations using
both flux form and material form advection of thermodynamic variables, which
allow for the necessary cancellations required to conserve entropy without the
chain rule. We show that such formulations allow for the stable simulation of
both the thermal shallow water and 3D compressible Euler equations on the
sphere using mixed compatible finite elements without entropy damping.

Entropy stability is an essential feature of numerical schemes for the
solution hyperbolic systems involving thermally forced dynamics. In recent
years great advances have been made in the formulation of entropy stable
schemes using finite difference [1], spectral element [2] and discontinuous
Galerkin [3, 4] methods for the solution of hyperbolic systems, including the
3D compressible Euler equations used for the simulation of atmospheric
dynamics [5]. These formulations are based on the preservation of the
summation by parts property between gradient and divergence operators for flux
form representations of these systems. In a separate line of research, a
finite volume scheme for the rotating thermal shallow water equations has been
developed using global integrals of the souce terms in order to derive well
balanced fluxes that preserve equilibrium solutions [6].

One draw back of these flux form energy and entropy conserving schemes is that
they do not naturally conserve vorticity, as do vector invariant formulations
using compatible finite element schemes [7, 8, 9]. More recently, an entropy
stable discontinuous Galerkin formulation has been developed for the vector
invariant thermal shallow water equations [10] which builds off previous work
on the shallow water equations [11] to incorporate the properties of
geostrophic balance and vorticity conservation, making this an appealing
scheme for the simulation of atmospheric dynamics. In addition to the
summation by parts property, this vector invariant formulation negates the
need to preserve the chain rule in order to conserve entropy for the case of
thermodynamic variables represented by discontinuous function spaces, for
which the chain rule is not defined at element boundaries.

In the present work we extend this approach to the mixed finite element
method, using compatible function spaces. Moreover we show that the negation
of the chain rule can be derived via an averaging of two different skew-
symmetric Hamiltonian formulations, one involving the material form transport
of the thermodynamic variable, and the other using flux form transport. In
doing so we show that thermal atmospheric systems including the thermal
shallow water equations and the 3D compressible Euler equations can run stably
without any entropy dissipation, as is typically required for the stable
simulation of these systems. Note that in the present article we are referring
specifically to a _mathematical_ entropy, by which we mean a conserved moment
that remains convex for all possible solutions due the presence of strictly
positive eigenvalues for its Hessian matrix [1], and not a _physical_ entropy,
as may be derived from the functional derivatives of a thermodynamic potential
[12].

The remainder of this article proceeds as follows: In Section 1 we will
introduce the thermal shallow water and 3D compressible Euler equations in
continuous form, and well as their entropy functions. In Section 2 we describe
discrete, mixed finite element formulations of these systems for which energy
and entropy are conserved. We present these formulations first as an average
of discrete skew-symmetric Hamiltonian systems, and secondly using a coupled
approach using flux form transport only, which is similar to the approach
previously derived for the discontinuous Galerkin method [10]. Results will be
presented in Section 3, and finally we will discuss the conclusions from this
work in Section 4.

##  1 Atmospheric systems with thermal dynamics

###  1.1 Thermal shallow water equations

The thermal shallow water equations may be considered as an extension of the
shallow water equations for which the buoyancy, $b$ is allowed to vary with
respect to changes in density, $\rho$ as $b=g\rho/\overline{\rho}$, where $g$
is gravity and $\overline{\rho}$ is the mean density. If there is no variation
in density, such that $\rho=\overline{\rho}$, then the original shallow water
equations are recovered for a constant buoyancy. The Hamiltonian structure of
the thermal shallow water equations has been studied in continuous form in
[13, 14], and in a discrete mixed Galerkin variational form in [9]. These may
be given for the velocity, $\boldsymbol{u}$, fluid depth, $h$ and buoyancy $b$
as

$\displaystyle\frac{\partial\boldsymbol{u}}{\partial
t}+q\times\boldsymbol{F}+\nabla\Bigg{(}\frac{1}{2}\boldsymbol{u}\cdot\boldsymbol{u}+hb\Bigg{)}-\frac{h}{2}\nabla
b$ $\displaystyle=0$ (1a) $\displaystyle\frac{\partial h}{\partial
t}+\nabla\cdot\boldsymbol{F}$ $\displaystyle=0$ (1b)
$\displaystyle\frac{\partial b}{\partial t}+\frac{\nabla
b}{h}\cdot\boldsymbol{F}$ $\displaystyle=0,$ (1c)

where $\boldsymbol{F}=h\boldsymbol{u}$ is the mass flux and
$q=(\nabla\times\boldsymbol{u}+f)/h$ is the potential vorticity (with $f$
being the Coriolis term due to the earth’s rotation). Combining (1b) and (1c)
we may alternatively express the buoyancy transport as a flux form equation
for the density weighted buoyancy $B=hb$ as

$\frac{\partial B}{\partial t}+\nabla\cdot(b\boldsymbol{F})=0.$ (2)

Both the material transport of $b$ (1c) and flux form transport of $B$ (2)
support skew-symmetric Hamiltonian formulations of the thermal shallow water
equations. For the material form we express the global energy over the domain
$\Omega$ as

$\mathcal{H}_{m}^{tsw}=\int\frac{1}{2}h\boldsymbol{u}^{2}+\frac{1}{2}h^{2}b\mathrm{d}\Omega,$
(3)

for which the variational derivatives are expressed as

$\nabla\mathcal{H}_{m}^{tsw}=\begin{bmatrix}h\boldsymbol{u}\\\
\frac{1}{2}\boldsymbol{u}\cdot\boldsymbol{u}+hb\\\
\frac{1}{2}h^{2}\end{bmatrix}=\begin{bmatrix}\boldsymbol{F}\\\ \Phi^{m}\\\
T^{m}\end{bmatrix}.$ (4)

This leads to a skew-symmetric representation of the thermal shallow water
equations as

$\frac{\partial}{\partial t}\begin{bmatrix}\boldsymbol{u}\\\ h\\\
b\end{bmatrix}=-\begin{bmatrix}q\times&\nabla&-\frac{\nabla b}{h}\\\
\nabla\cdot&0&0\\\ \frac{\nabla
b}{h}&0&0\end{bmatrix}\begin{bmatrix}\boldsymbol{F}\\\ \Phi^{m}\\\
T^{m}\end{bmatrix}.$ (5)

Conversely, for the flux form, we have the global energy as

$\mathcal{H}_{f}^{tsw}=\int\frac{1}{2}h\boldsymbol{u}\cdot\boldsymbol{u}+\frac{1}{2}hB\mathrm{d}\Omega$
(6)

with the corresponding variational derivatives as

$\nabla\mathcal{H}_{f}^{tsw}=\begin{bmatrix}h\boldsymbol{u}\\\
\frac{1}{2}\boldsymbol{u}\cdot\boldsymbol{u}+\frac{1}{2}B\\\
\frac{1}{2}h\end{bmatrix}=\begin{bmatrix}\boldsymbol{F}\\\ \Phi^{f}\\\
T^{f}\end{bmatrix}.$ (7)

For which the skew-symmetric form of the thermal shallow water equations is
expressed as

$\frac{\partial}{\partial t}\begin{bmatrix}\boldsymbol{u}\\\ h\\\
B\end{bmatrix}=-\begin{bmatrix}q\times&\nabla&\frac{B}{h}\nabla\\\
\nabla\cdot&0&0\\\
\nabla\cdot(\frac{B}{h}\cdot)&0&0\end{bmatrix}\begin{bmatrix}\boldsymbol{F}\\\
\Phi^{f}\\\ T^{f}\end{bmatrix}.$ (8)

For both (5) and (8) energy conservation arises from the skew-symmetric nature
of the right hand side operator. In either case multiplying both sides by
$(\nabla\mathcal{H}^{tsw})^{\top}$ leads to the expression

$\begin{bmatrix}\frac{\delta\mathcal{H}^{tsw}}{\delta\boldsymbol{u}}&\frac{\delta\mathcal{H}^{tsw}}{\delta
h}&\frac{\delta\mathcal{H}^{tsw}}{\delta
a}\end{bmatrix}\begin{bmatrix}\frac{\partial\boldsymbol{u}}{\partial t}\\\
\frac{\partial h}{\partial t}\\\ \frac{\partial a}{\partial
t}\end{bmatrix}=\frac{\mathrm{d}\mathcal{H}^{tsw}}{\mathrm{d}t}=0$ (9)

for $a\in b,B$.

Additionally, both skew-symmetric operators have a null-space, which contains
the variational derivatives of the additional invariants of the system. These
include the mass, $\mathcal{M}^{tsw}=\int h\mathrm{d}\Omega$, the vorticity,
$\mathcal{W}^{tsw}=\int hq-f\mathrm{d}\Omega$, the density weighted buoyancy,
$\mathcal{B}^{tsw}=\int hb\mathrm{d}\Omega$ [9], as well as the higher order
moments of the buoyancy [13] such as the entropy
$\mathcal{S}^{tsw}=\frac{1}{2}\int
hb^{2}\mathrm{d}\Omega=\int\frac{B^{2}}{2h}\mathrm{d}\Omega$. While the
preservation of all of these quantities is dynamically important, we are
principally concerned in the present study with the conservation of entropy,
since this is know to have eigenvalues which are all $\geq 0$, which ensures
the stability of thermal processes.

The variational derivatives of the entropy for the material form are given as
$\nabla\mathcal{S}_{m}^{tsw}=[\boldsymbol{0},\frac{1}{2}b^{2},hb]^{\top}$,
such that this is conserved as

$\frac{\mathrm{d}\mathcal{S}_{m}^{tsw}}{\mathrm{d}t}=\begin{bmatrix}\boldsymbol{0}&\frac{1}{2}b^{2}&hb\end{bmatrix}\begin{bmatrix}\frac{\partial\boldsymbol{u}}{\partial
t}\\\ \frac{\partial h}{\partial t}\\\ \frac{\partial a}{\partial
t}\end{bmatrix}=-\begin{bmatrix}\boldsymbol{0}&\frac{1}{2}b^{2}&hb\end{bmatrix}\begin{bmatrix}q\times&\nabla&-\frac{\nabla
b}{h}\\\ \nabla\cdot&0&0\\\ \frac{\nabla
b}{h}&0&0\end{bmatrix}\begin{bmatrix}\boldsymbol{F}\\\ \Phi^{m}\\\
T^{m}\end{bmatrix}=-\begin{bmatrix}\frac{1}{2}\nabla b^{2}-b\nabla
b&0&0\end{bmatrix}\begin{bmatrix}\boldsymbol{F}\\\ \Phi^{m}\\\
T^{m}\end{bmatrix}=0.$ (10)

Similarly for the flux form, the variational derivatives for the corresponding
entropy are given as
$\nabla\mathcal{S}_{f}^{tsw}=[\boldsymbol{0},\frac{-B^{2}}{2h^{2}},\frac{B}{h}]^{\top}$,
for which the entropy is conserved as

$\frac{\mathrm{d}\mathcal{S}_{f}^{tsw}}{\mathrm{d}t}=\begin{bmatrix}\boldsymbol{0}&-\frac{B^{2}}{2h^{2}}&\frac{B}{h}\end{bmatrix}\begin{bmatrix}\frac{\partial\boldsymbol{u}}{\partial
t}\\\ \frac{\partial h}{\partial t}\\\ \frac{\partial a}{\partial
t}\end{bmatrix}=-\begin{bmatrix}\boldsymbol{0}&-\frac{B^{2}}{2h^{2}}&\frac{B}{h}\end{bmatrix}\begin{bmatrix}q\times&\nabla&\frac{B}{h}\nabla\\\
\nabla\cdot&0&0\\\
\nabla\cdot(\frac{B}{h}\cdot)&0&0\end{bmatrix}\begin{bmatrix}\boldsymbol{F}\\\
\Phi^{f}\\\
T^{f}\end{bmatrix}=-\begin{bmatrix}-\frac{1}{2}\nabla\Big{(}\frac{B^{2}}{h^{2}}\Big{)}+\frac{B}{h}\nabla\Big{(}\frac{B}{h}\Big{)}&0&0\end{bmatrix}\begin{bmatrix}\boldsymbol{F}\\\
\Phi^{f}\\\ T^{f}\end{bmatrix}=0.$ (11)

In both cases, entropy conservation is derived from the chain rule. This is
problematic in the discrete form if $h$ and $b$ or $B$ are represented by
discontinuous function spaces, for which the chain rule is not defined at
element boundaries.

In order to negate this issue and ensure that entropy is conserved independent
of the chain rule, we may alternatively construct our solution as an average
of (5) and (8) as

$\displaystyle\frac{\partial\boldsymbol{u}}{\partial
t}+q\times\boldsymbol{F}+\frac{1}{2}\nabla\Bigg{(}\boldsymbol{u}\cdot\boldsymbol{u}+hb+\frac{1}{2}B\Bigg{)}+\frac{1}{4}\Bigg{(}-\frac{1}{h}\nabla\Bigg{(}\frac{B}{h}\Bigg{)}h^{2}+b\nabla
h\Bigg{)}$ $\displaystyle=0$ (12a) $\displaystyle\frac{\partial h}{\partial
t}+\nabla\cdot\boldsymbol{F}$ $\displaystyle=0$ (12b)
$\displaystyle\frac{\partial b}{\partial
t}+\frac{1}{h}\nabla\Bigg{(}\frac{B}{h}\Bigg{)}\cdot\boldsymbol{F}$
$\displaystyle=0$ (12c) $\displaystyle\frac{\partial B}{\partial
t}+\nabla\cdot\Big{(}{b\boldsymbol{F}}\Big{)}$ $\displaystyle=0.$ (12d)

Note that the above formulation differs from a weighted average of (5) and (8)
in one respect: the contents of the skew-symmetric matrix for the material
form, (5) are instead of $b$ using $\frac{B}{h}$, while the skew-symmetric
operator in (8) are instead of $B$ using $hb$. This will become important in
the discrete form, since it will allow for a tighter coupling of the material
and flux form transport equations. Now pre-multiplying the material form terms
by $(\nabla\mathcal{S}_{m}^{tsw})^{\top}$, and the flux form terms by
$(\nabla\mathcal{S}_{f}^{tsw})^{\top}$ and recalling that $b=\frac{B}{h}$
leads to cancellations between (10) and (11) such that
$\mathcal{S}_{m}^{tsw}=\mathcal{S}_{f}^{tsw}$ is conserved independently of
the chain rule.

More succinctly, we may combine (1b) and (1c) to eliminate the material
advection equation for $b$ from the above system (12c) to give

$\displaystyle\frac{\partial\boldsymbol{u}}{\partial
t}+q\times\boldsymbol{F}+\frac{1}{2}\nabla\Bigg{(}\boldsymbol{u}\cdot\boldsymbol{u}+\frac{3}{2}B\Bigg{)}+\frac{1}{4}\Bigg{(}\frac{B}{h}\nabla
h-h\nabla\Big{(}\frac{B}{h}\Big{)}\Bigg{)}$ $\displaystyle=0$ (13a)
$\displaystyle\frac{\partial h}{\partial t}+\nabla\cdot\boldsymbol{F}$
$\displaystyle=0$ (13b) $\displaystyle\frac{\partial B}{\partial
t}+\frac{1}{2}\nabla\cdot\Bigg{(}{\frac{B}{h}\boldsymbol{F}}\Bigg{)}+\frac{1}{2}\frac{B}{h}\nabla\cdot\boldsymbol{F}+\frac{1}{2}\boldsymbol{F}\cdot\Big{(}\frac{B}{h}\Big{)}$
$\displaystyle=0,$ (13c)

as has been previously done for an entropy conserving discretisation of the
thermal shallow water equations using a discontinuous Galerkin method [10].
Once again, pre-multiplying (13) by $(\nabla\mathcal{H}_{f}^{tsw})^{\top}$
allows for the conservation of energy, while pre-multiplying by
$(\nabla\mathcal{S}_{f}^{tsw})^{\top}$ allows for the conservation of entropy,
without necessitating the invocation of the chain rule.

###  1.2 3D compressible Euler equations

Extending the above formulation to the 3D compressible Euler equations is
remarkably simple. This is because the 3D compressible Euler equations have
the same skew-symmetric structure as the thermal shallow water equations as
given in (5) and (8). The only thing that differs is the energy and is
variational derivatives, which for the flux form is given as

$\mathcal{H}_{f}^{ce}=\int\frac{1}{2}\rho\boldsymbol{u}\cdot\boldsymbol{u}+\rho
gz+\frac{c_{v}}{c_{p}}\Theta\Pi\mathrm{d}\Omega,$ (14)

where $g$ is gravity, $z$ is height, $\rho$ is the density, $\Theta$ the
density weighted weighted potential temperature (and
$\theta:=\frac{\Theta}{\rho}$ is the potential temperature) and $c_{v}$,
$c_{p}$ are the specific heats at constant volume and pressure respectively.
The Exner pressure is denoted by $\Pi$, which is derived from the ideal gas
law as $\Pi=c_{p}(R\Theta/p_{0})^{R/c_{v}}$, with $p_{0}$ being the reference
pressure at the surface and $R=c_{p}-c_{v}$ is the ideal gas constant.

The variational derivatives of the energy for the 3D compressible Euler
equations with the flux form transport of density weighted potential
temperature are then given as

$\nabla\mathcal{H}_{f}^{ce}=\begin{bmatrix}\rho\boldsymbol{u}\\\
\frac{1}{2}\boldsymbol{u}\cdot\boldsymbol{u}+gz\\\
c_{p}\Big{(}\frac{R\Theta}{p_{0}}\Big{)}^{\frac{R}{c_{v}}}\end{bmatrix}=\begin{bmatrix}\boldsymbol{U}\\\
\Phi\\\ \Pi\end{bmatrix}.$ (15)

Note that we have omitted superscripts for the equation set and subscripts for
flux form from the above definitions for the sake of brevity. The 3D
compressible Euler equations can then be expressed in a direct analogue of (8)
as

$\frac{\partial}{\partial t}\begin{bmatrix}\boldsymbol{u}\\\ \rho\\\
\Theta\end{bmatrix}=-\begin{bmatrix}q\times&\nabla&\frac{\Theta}{\rho}\nabla\\\
\nabla\cdot&0&0\\\
\nabla\cdot(\frac{\Theta}{\rho}\cdot)&0&0\end{bmatrix}\begin{bmatrix}\boldsymbol{U}\\\
\Phi\\\ \Pi\end{bmatrix}.$ (16)

We may then construct an analogue of the coupled thermal shallow water system
(13) for the 3D compressible Euler equations as

$\displaystyle\frac{\partial\boldsymbol{u}}{\partial
t}+q\times\boldsymbol{U}+\nabla\Bigg{(}\frac{1}{2}\boldsymbol{u}\cdot\boldsymbol{u}+gz+\frac{1}{2}\frac{\Theta}{\rho}\Pi\Bigg{)}+\frac{1}{2}\Bigg{(}\frac{\Theta}{\rho}\nabla\Pi-\Pi\nabla\Big{(}\frac{\Theta}{\rho}\Big{)}\Bigg{)}$
$\displaystyle=0$ (17a) $\displaystyle\frac{\partial\rho}{\partial
t}+\nabla\cdot\boldsymbol{U}$ $\displaystyle=0$ (17b)
$\displaystyle\frac{\partial\Theta}{\partial
t}+\frac{1}{2}\nabla\cdot\Bigg{(}{\frac{\Theta}{\rho}\boldsymbol{U}}\Bigg{)}+\frac{1}{2}\frac{\Theta}{\rho}\nabla\cdot\boldsymbol{U}+\frac{1}{2}\boldsymbol{U}\cdot\nabla\Big{(}\frac{\Theta}{\rho}\Big{)}$
$\displaystyle=0.$ (17c)

The corresponding entropy is given as
$\mathcal{S}^{ce}=\int\frac{\Theta^{2}}{2\rho}\mathrm{d}\Omega$, and its
variational derivatives as
$\nabla\mathcal{S}^{ce}=[\boldsymbol{0},\frac{-\Theta^{2}}{2\rho^{2}},\frac{\Theta}{\rho}]^{\top}$.
Once again, pre-multiplication of (17) by the transpose of this vector allows
for the conservation of $\mathcal{S}^{ce}$ without resorting to the chain
rule.

##  2 Mixed finite element discretisations

In this study we apply a mixed compatible finite element discretisation in
order to preserve the energy and entropy conserving properties of the
continuous formulations discussed in the previous section to the discrete
form.

###  2.1 Thermal shallow water equations

For the thermal shallow water equations we discretise the two dimensional
domain $\Omega^{2}$ using the finite dimensional subspaces
$\mathbb{V}_{0}\subset H^{1}(\Omega^{2})$, $\mathbb{V}_{1}\subset
H(\mathrm{div},\Omega^{2})$, $\mathbb{V}_{2}\subset L^{2}(\Omega^{2})$. These
subspaces are spanned by the basis functions $\psi_{h}\in\mathbb{V}_{0}$,
$\boldsymbol{v}_{h}\in\mathbb{V}_{1}$ and $\phi_{h}\in\mathbb{V}_{2}$
respectively. The compatible relationships between these subspaces may then be
expressed via a discrete de Rham complex as [15]

$\mathbb{R}\longrightarrow\mathbb{V}_{0}\stackrel{{\scriptstyle\nabla^{\perp}}}{{\longrightarrow}}\mathbb{V}_{1}\stackrel{{\scriptstyle\nabla\cdot}}{{\longrightarrow}}\mathbb{V}_{2}\longrightarrow
0,$ (18)

such that the vector identity $\nabla\cdot\nabla^{\perp}:=0$ is preserved
discretely. The discrete analogues of the _curl_ and _gradient_ operators are
then applied as adjoints of the strong form $\nabla^{\perp}$ and $\nabla\cdot$
operators [7, 8, 9].

In order to discretely solve the entropy conserving thermal shallow water
equations we first introduce the operators:

$\displaystyle\boldsymbol{\mathsf{M}}_{0}$
$\displaystyle:\mathbb{V}_{0}\rightarrow\mathbb{V}_{0}^{\prime}$
$\displaystyle\quad\mathrm{for\
which}\quad\boldsymbol{\mathsf{M}}_{0}\psi_{h}$
$\displaystyle:=\langle\psi_{h},\cdot\rangle,\quad$
$\displaystyle\forall\psi_{h}\in\mathbb{V}_{0},$ (19a)
$\displaystyle\boldsymbol{\mathsf{M}}_{1}$
$\displaystyle:\mathbb{V}_{1}\rightarrow\mathbb{V}_{1}^{\prime}$
$\displaystyle\quad\mathrm{for\
which}\quad\boldsymbol{\mathsf{M}}_{1}\boldsymbol{v}_{h}$
$\displaystyle:=\langle\boldsymbol{v}_{h},\cdot\rangle,\quad$
$\displaystyle\forall\boldsymbol{v}_{h}\in\mathbb{V}_{1},$ (19b)
$\displaystyle\boldsymbol{\mathsf{M}}_{2}$
$\displaystyle:\mathbb{V}_{2}\rightarrow\mathbb{V}_{2}^{\prime}$
$\displaystyle\quad\mathrm{for\
which}\quad\boldsymbol{\mathsf{M}}_{2}\phi_{h}$
$\displaystyle:=\langle\phi_{h},\cdot\rangle,\quad$
$\displaystyle\forall\phi_{h}\in\mathbb{V}_{2},$ (19c)
$\displaystyle\boldsymbol{\mathsf{D}}_{2}$
$\displaystyle:\mathbb{V}_{1}\rightarrow\mathbb{V}_{2}^{\prime}$
$\displaystyle\quad\mathrm{for\
which}\quad\boldsymbol{\mathsf{D}}_{2}\boldsymbol{v}_{h}$
$\displaystyle:=\langle\nabla\cdot\boldsymbol{v}_{h},\cdot\rangle,\quad$
$\displaystyle\forall\boldsymbol{v}_{h}\in\mathbb{V}_{1},$ (19d)
$\displaystyle\boldsymbol{\mathsf{R}}_{2}$
$\displaystyle:\mathbb{V}_{0}\rightarrow\mathbb{V}_{1}^{\prime}$
$\displaystyle\quad\mathrm{for\
which}\quad\boldsymbol{\mathsf{R}}_{2}\psi_{h}$
$\displaystyle:=\langle\nabla^{\perp}\psi_{h},\cdot\rangle,\quad$
$\displaystyle\forall\psi_{h}\in\mathbb{V}_{0},$ (19e)
$\displaystyle\boldsymbol{\mathsf{C}}_{2}$
$\displaystyle:\mathbb{V}_{0}\otimes\mathbb{V}_{1}\rightarrow\mathbb{V}_{1}^{\prime}$
$\displaystyle\quad\mathrm{for\
which}\quad\boldsymbol{\mathsf{C}}_{2}(\psi_{h},\boldsymbol{v}_{h})$
$\displaystyle:=\langle\psi_{h}\times\boldsymbol{v}_{h},\cdot\rangle,\quad$
$\displaystyle\forall\psi_{h}\in\mathbb{V}_{0},\boldsymbol{v}_{h}\in\mathbb{V}_{1},$
(19f) $\displaystyle\boldsymbol{\mathsf{M}}_{0*}$
$\displaystyle:\mathbb{V}_{0}\otimes\mathbb{V}_{2}\rightarrow\mathbb{V}_{0}^{\prime}$
$\displaystyle\quad\mathrm{for\
which}\quad\boldsymbol{\mathsf{M}}_{0*}(\psi_{h},\phi_{h})$
$\displaystyle:=\langle\phi_{h}\psi_{h},\cdot\rangle,\quad$
$\displaystyle\forall\psi_{h}\in\mathbb{V}_{0},\phi_{h}\in\mathbb{V}_{2},$
(19g) $\displaystyle\boldsymbol{\mathsf{M}}_{1*}$
$\displaystyle:\mathbb{V}_{1}\otimes\mathbb{V}_{2}\rightarrow\mathbb{V}_{1}^{\prime}$
$\displaystyle\quad\mathrm{for\
which}\quad\boldsymbol{\mathsf{M}}_{1*}(\boldsymbol{v}_{h},\phi_{h})$
$\displaystyle:=\langle\phi_{h}\boldsymbol{v}_{h},\cdot\rangle,\quad$
$\displaystyle\forall\boldsymbol{v}_{h}\in\mathbb{V}_{1},\phi_{h}\in\mathbb{V}_{2},$
(19h) $\displaystyle\boldsymbol{\mathsf{M}}_{2*}$
$\displaystyle:\mathbb{V}_{2}\otimes\mathbb{V}_{2}\rightarrow\mathbb{V}_{2}^{\prime}$
$\displaystyle\quad\mathrm{for\
which}\quad\boldsymbol{\mathsf{M}}_{2*}(\phi_{h},\phi_{k})$
$\displaystyle:=\langle\phi_{h}\phi_{k},\cdot\rangle,\quad$
$\displaystyle\forall\phi_{h},\phi_{k}\in\mathbb{V}_{2},$ (19i)
$\displaystyle\boldsymbol{\mathsf{K}}_{2}$
$\displaystyle:\mathbb{V}_{1}\otimes\mathbb{V}_{1}\rightarrow\mathbb{V}_{2}^{\prime}$
$\displaystyle\quad\mathrm{for\
which}\quad\boldsymbol{\mathsf{K}}_{2}(\boldsymbol{v}_{h},\boldsymbol{v}_{k})$
$\displaystyle:=\langle\boldsymbol{v}_{h}\cdot\boldsymbol{v}_{k},\cdot\rangle,\quad$
$\displaystyle\forall\boldsymbol{v}_{h},\boldsymbol{v}_{k}\in\mathbb{V}_{1},$
(19j)

where primes denote corresponding dual spaces and angled brackets denote inner
products as $\langle a,b\rangle=\int ab\mathrm{d}\Omega^{2}$, and likewise for
tri-linear operators.

####  2.1.1 Thermal shallow water equations: mixed formulation

Having defined the discrete function spaces and operators, we now proceed to
describe the semi-discrete analogue of the mixed material-flux form of the
thermal shallow water equations (12) for the discrete solution variables
$\boldsymbol{u}_{h}\in\mathbb{V}_{1}$, $h_{h},B_{h},b_{h}\in\mathbb{V}_{2}$.
Firstly the diagnostic equations for the discrete potential vorticity,
$q_{h}\in\mathbb{V}_{0}$ and variational derivatives
$\boldsymbol{F}_{h}\in\mathbb{V}_{1}$,
$\Phi^{m}_{h},\Phi^{f}_{h},T^{m}_{h},T^{f}_{h}\in\mathbb{V}_{2}$ are given as

$\displaystyle\boldsymbol{\mathsf{M}}_{0*}(q_{h},h_{h})$
$\displaystyle=-\boldsymbol{\mathsf{R}}_{2}^{\top}\boldsymbol{u}_{h}+\boldsymbol{\mathsf{M}}_{0}f_{h}$
(20a) $\displaystyle\boldsymbol{\mathsf{M}}_{1}\boldsymbol{F}_{h}$
$\displaystyle=\boldsymbol{\mathsf{M}}_{1*}(\boldsymbol{u}_{h},h_{h})$ (20b)
$\displaystyle\boldsymbol{\mathsf{M}}_{2}\Phi^{m}_{h}$
$\displaystyle=\frac{1}{2}\boldsymbol{\mathsf{K}}_{2}(\boldsymbol{u}_{h},\boldsymbol{u}_{h})+\boldsymbol{\mathsf{M}}_{2*}(h_{h},b_{h})$
(20c) $\displaystyle\boldsymbol{\mathsf{M}}_{2}\Phi^{f}_{h}$
$\displaystyle=\frac{1}{2}\boldsymbol{\mathsf{K}}_{2}(\boldsymbol{u}_{h},\boldsymbol{u}_{h})+\frac{1}{2}\boldsymbol{\mathsf{M}}_{2}B_{h}$
(20d) $\displaystyle\boldsymbol{\mathsf{M}}_{2}T^{m}_{h}$
$\displaystyle=\frac{1}{2}\boldsymbol{\mathsf{M}}_{2*}(h_{h},h_{h})$ (20e)
$\displaystyle\boldsymbol{\mathsf{M}}_{2}T^{f}_{h}$
$\displaystyle=\frac{1}{2}\boldsymbol{\mathsf{M}}_{2}h_{h},$ (20f)

where $f_{h}\in\mathbb{V}_{0}$ is the discrete representation of the Coriolis
term. We also compute a discrete analogue of the buoyancy gradient per unit
depth, $\boldsymbol{G}:=\frac{1}{h}\nabla b^{\prime}$, where
$b^{\prime}=\frac{B}{h}$, as appears in the skew symmetric operator for the
material form of the thermal shallow water equations (5) for
$\boldsymbol{G}_{h}\in\mathbb{V}_{1}$, $b_{h}^{\prime}\in\mathbb{V}_{2}$ as

$\displaystyle\boldsymbol{\mathsf{M}}_{2*}(h_{h},b_{h}^{\prime})$
$\displaystyle=\boldsymbol{\mathsf{M}}_{2}B_{h}$ (21a)
$\displaystyle\boldsymbol{\mathsf{M}}_{1*}(\boldsymbol{G}_{h},h_{h})$
$\displaystyle=-\boldsymbol{\mathsf{D}}_{2}^{\top}b_{h}^{\prime}.$ (21b)

Once we have computed the necessary diagnostic quantities, we may solve for
the prognostic variables in the semi-discrete form of (12) as

$\displaystyle\boldsymbol{\mathsf{M}}_{1}\frac{\partial\boldsymbol{u}_{h}}{\partial
t}$
$\displaystyle=-\boldsymbol{\mathsf{C}}_{2}(q_{h},\boldsymbol{F}_{h})+\frac{1}{2}\boldsymbol{\mathsf{D}}_{2}^{\top}(\Phi^{m}_{h}+\Phi^{f}_{h})+\frac{1}{2}\boldsymbol{\mathsf{K}}_{2}^{\top}(\boldsymbol{G}_{h},T_{h}^{m})+\frac{1}{2}\boldsymbol{\mathsf{M}}_{1*}(\boldsymbol{\mathsf{M}}_{1}^{-1}\boldsymbol{\mathsf{D}}_{2}^{\top}T_{h}^{f},b_{h})$
(22a) $\displaystyle\boldsymbol{\mathsf{M}}_{2}\frac{\partial h_{h}}{\partial
t}$ $\displaystyle=-\boldsymbol{\mathsf{D}}_{2}\boldsymbol{F}_{h}$ (22b)
$\displaystyle\boldsymbol{\mathsf{M}}_{2}\frac{\partial b_{h}}{\partial t}$
$\displaystyle=-\boldsymbol{\mathsf{K}}_{2}(\boldsymbol{G}_{h},\boldsymbol{F}_{h})$
(22c) $\displaystyle\boldsymbol{\mathsf{M}}_{2}\frac{\partial B_{h}}{\partial
t}$
$\displaystyle=-\boldsymbol{\mathsf{D}}_{2}\boldsymbol{\mathsf{M}}_{1}^{-1}\boldsymbol{\mathsf{M}}_{1*}(\boldsymbol{F}_{h},b_{h})$
(22d)

Left-multiplication of (22) by the vector
$[\boldsymbol{F}_{h}^{\top},\frac{1}{2}(\Phi_{h}^{m\top}+\Phi_{h}^{f\top}),\frac{1}{2}T_{h}^{m\top},\frac{1}{2}T_{h}^{f\top}]$
leads to the cancellation of all right hand side terms as

$\displaystyle-\boldsymbol{F}_{h}^{\top}\boldsymbol{\mathsf{C}}_{2}(q_{h},\boldsymbol{F}_{h})+\frac{1}{2}\boldsymbol{F}_{h}^{\top}\boldsymbol{\mathsf{D}}_{2}^{\top}(\Phi^{m}_{h}+\Phi^{f}_{h})+\frac{1}{2}\boldsymbol{F}_{h}^{\top}\boldsymbol{\mathsf{K}}_{2}^{\top}(\boldsymbol{G}_{h},T_{h}^{m})+\frac{1}{2}\boldsymbol{F}_{h}^{\top}\boldsymbol{\mathsf{M}}_{1*}(\boldsymbol{\mathsf{M}}_{1}^{-1}\boldsymbol{\mathsf{D}}_{2}^{\top}T_{h}^{f},b_{h})$
$\displaystyle-\frac{1}{2}(\Phi_{h}^{m\top}+\Phi_{h}^{f\top})\boldsymbol{\mathsf{D}}_{2}\boldsymbol{F}_{h}-\frac{1}{2}T_{h}^{m\top}\boldsymbol{\mathsf{K}}_{2}(\boldsymbol{G}_{h},\boldsymbol{F}_{h})-\frac{1}{2}T_{h}^{f\top}\boldsymbol{\mathsf{D}}_{2}\boldsymbol{\mathsf{M}}_{1}^{-1}\boldsymbol{\mathsf{M}}_{1*}(\boldsymbol{F}_{h},b_{h})$
$\displaystyle=$ (23a)
$\displaystyle-\langle\boldsymbol{F}_{h},\boldsymbol{q}_{h}\times\boldsymbol{F}_{h}\rangle+\frac{1}{2}\langle\nabla\cdot\boldsymbol{F}_{h},\Phi_{h}^{m}+\Phi_{h}^{f}\rangle+\frac{1}{2}\langle\boldsymbol{F}_{h}\cdot\boldsymbol{G}_{h},T_{h}^{m}\rangle+\frac{1}{2}\langle\nabla\cdot(b_{h}\boldsymbol{F}_{h}),T_{h}^{f}\rangle$
$\displaystyle-\frac{1}{2}\langle\Phi_{h}^{m}+\Phi_{h}^{f},\nabla\cdot\boldsymbol{F}_{h}\rangle-\frac{1}{2}\langle
T_{h}^{m},\boldsymbol{G}_{h}\cdot\boldsymbol{F}_{h}\rangle-\frac{1}{2}\langle
T_{h}^{f},\nabla\cdot(b_{h}\boldsymbol{F}_{h})\rangle$ $\displaystyle=0,$
(23b)

and results in the expression

$\frac{1}{2}\Bigg{(}\boldsymbol{F}_{h}^{\top}\boldsymbol{\mathsf{M}}_{1}\frac{\partial\boldsymbol{u}_{h}}{\partial
t}+\Phi_{h}^{m\top}\boldsymbol{\mathsf{M}}_{2}\frac{\partial h_{h}}{\partial
t}+T_{h}^{m\top}\boldsymbol{\mathsf{M}}_{2}\frac{\partial b_{h}}{\partial
t}+\boldsymbol{F}_{h}^{\top}\boldsymbol{\mathsf{M}}_{1}\frac{\partial\boldsymbol{u}_{h}}{\partial
t}+\Phi_{h}^{f\top}\boldsymbol{\mathsf{M}}_{2}\frac{\partial h_{h}}{\partial
t}+T_{h}^{f\top}\boldsymbol{\mathsf{M}}_{2}\frac{\partial B_{h}}{\partial
t}\Bigg{)}=\frac{1}{2}\Bigg{(}\frac{\partial\mathcal{H}_{m,h}^{tsw}}{\partial
t}+\frac{\partial\mathcal{H}_{f,h}^{tsw}}{\partial t}\Bigg{)}=0$ (24)

In order to further conserve energy in time as well as space, the functional
derivatives on the left hand side in the above expression should be exactly
integrated in time up to the order of the temporal discretisation [16, 9, 17,
18].

The discrete functional derivatives of the entropy are expressed as

$\displaystyle\nabla\mathcal{S}_{m,h}^{tsw}$
$\displaystyle=\begin{bmatrix}\boldsymbol{0}\\\
\frac{1}{2}\Big{(}\boldsymbol{\mathsf{M}}_{2*}(b_{h},b_{h})\Big{)}^{\top}\boldsymbol{\mathsf{M}}_{2}^{-1}\\\
\Big{(}\boldsymbol{\mathsf{M}}_{2*}(h_{h},b_{h})\Big{)}^{\top}\boldsymbol{\mathsf{M}}_{2}^{-1}\end{bmatrix}$
(25a) $\displaystyle\nabla\mathcal{S}_{f,h}^{tsw}$
$\displaystyle=\begin{bmatrix}\boldsymbol{0}\\\
-\frac{1}{2}\Big{(}\boldsymbol{\mathsf{M}}_{2*}(B_{h},B_{h})\Big{)}^{\top}\boldsymbol{\mathsf{M}}_{2*}^{-1}(h_{h}^{2},\cdot)\\\
\Big{(}\boldsymbol{\mathsf{M}}_{2}B_{h}\Big{)}^{\top}\boldsymbol{\mathsf{M}}_{2*}^{-1}(h_{h},\cdot)\end{bmatrix}.$
(25b)

For the case of exact integration where
$\int\phi_{h}h_{h}^{2}b_{h}^{2}\mathrm{d}\Omega^{2}=\int\phi_{h}B_{h}^{2}\mathrm{d}\Omega^{2}$
$\forall\phi_{h}\in\mathbb{V}_{2}$, such that
$b_{h}^{2}=\boldsymbol{\mathsf{M}}_{2*}^{-1}(h_{h}^{2},\boldsymbol{\mathsf{M}}_{2}B_{h}^{2})$,
we have

$\frac{1}{2}\Big{(}\boldsymbol{\mathsf{M}}_{2*}(b_{h},b_{h})\Big{)}^{\top}\boldsymbol{\mathsf{M}}_{2}^{-1}\boldsymbol{\mathsf{D}}_{2}\boldsymbol{F}_{h}-\frac{1}{2}\Big{(}\boldsymbol{\mathsf{M}}_{2*}(B_{h},B_{h})\Big{)}^{\top}\boldsymbol{\mathsf{M}}_{2*}^{-1}(h_{h}^{2},\boldsymbol{\mathsf{D}}_{2}\boldsymbol{F}_{h})=\frac{1}{2}\Big{\langle}b_{h}^{2},\nabla\cdot\boldsymbol{F}_{h}\Big{\rangle}-\frac{1}{2}\Big{\langle}\frac{B_{h}^{2}}{h_{h}^{2}},\nabla\cdot\boldsymbol{F}_{h}\Big{\rangle}=0,$
(26)

such that the left multiplication by (22b) by an average of the second rows in
(25a) and (25b) results in no change in entropy. Left multiplication of (22c)
by the bottom row in (25a) and (22d) by the bottom row in (25b) and adding
then gives

$\displaystyle-\Big{(}\boldsymbol{\mathsf{M}}_{2*}(h_{h},b_{h})\Big{)}^{\top}\boldsymbol{\mathsf{M}}_{2}^{-1}\boldsymbol{\mathsf{K}}_{2}(\boldsymbol{F}_{h},\boldsymbol{\mathsf{M}}_{1*}^{-1}\boldsymbol{\mathsf{D}}_{2}^{\top}\boldsymbol{\mathsf{M}}_{2*}^{-1}(h_{h},\boldsymbol{\mathsf{M}}_{2}B_{h}),h_{h}))+\Big{(}\boldsymbol{\mathsf{M}}_{2}B_{h}\Big{)}^{\top}\boldsymbol{\mathsf{M}}_{2*}^{-1}(h_{h},\boldsymbol{\mathsf{D}}_{2}\boldsymbol{\mathsf{M}}_{1}^{-1}\boldsymbol{\mathsf{M}}_{1*}(\boldsymbol{F}_{h},b_{h}))$
$\displaystyle=$ (27a)
$\displaystyle-\Big{\langle}\nabla\cdot\Bigg{(}\frac{h_{h}b_{h}\boldsymbol{F}_{h}}{h_{h}}\Bigg{)},\frac{B_{h}}{h_{h}}\Big{\rangle}+\Big{\langle}\frac{B_{h}}{h_{h}},\nabla\cdot(b_{h}\boldsymbol{F}_{h})\Big{\rangle}$
$\displaystyle=0,$ (27b)

such that the material and flux form transport terms taken together also
result in no change of entropy, under the assumption of exact spatial
integration, but without resorting to the chain rule.

Note that since $b_{h}$ and $B_{h}$ are time stepped independently, while
there is no spatial error in the entropy conservation under the assumption of
exact spatial integration, we do anticipate a temporal error in the entropy
conservation. This error is controlled by ensuring that the temporal evolution
of $b_{h}$ is a function of $B_{h}$ in (22c), while the temporal evolution of
$B_{h}$ is a function of $b_{h}$ in (22d).

As for other compatible discretisations of both the shallow water [7, 8] and
thermal shallow water [9] equations, discrete conservation of mass
$\mathcal{M}_{h}^{tsw}=\int h_{h}\mathrm{d}\Omega^{2}$, and vorticity,
$\mathcal{W}_{h}^{tsw}=\int h_{h}q_{h}-f_{h}\mathrm{d}\Omega^{2}$ is derived
by setting $\phi_{h}=1$ in (22b) and
$\boldsymbol{v}_{h}=\nabla^{\perp}\psi_{h}$ in (22a) respectively. This also
holds for the discrete coupled formulation presented below.

####  2.1.2 Thermal shallow water equations: coupled formulation

In order to avoid solving separately for both $b_{h}$ and $B_{h}$, we may
construct an alternative energy and entropy conserving formulation as a
discrete analogue of (13) as

$\displaystyle\boldsymbol{\mathsf{M}}_{1}\frac{\partial\boldsymbol{u}_{h}}{\partial
t}$
$\displaystyle=-\boldsymbol{\mathsf{C}}_{2}(q_{h},\boldsymbol{F}_{h})+\boldsymbol{\mathsf{D}}_{2}^{\top}\Phi_{h}^{f}+\frac{1}{4}\boldsymbol{\mathsf{D}}_{2}^{\top}B_{h}+\frac{1}{4}\boldsymbol{\mathsf{M}}_{1*}(\boldsymbol{\mathsf{M}}_{1}^{-1}\boldsymbol{\mathsf{D}}_{2}^{\top}h_{h},b_{h}^{\prime})-\frac{1}{4}\boldsymbol{\mathsf{M}}_{1*}(\boldsymbol{\mathsf{M}}_{1}^{-1}\boldsymbol{\mathsf{D}}_{2}^{\top}b_{h}^{\prime},h_{h})$
(28a) $\displaystyle\boldsymbol{\mathsf{M}}_{2}\frac{\partial h_{h}}{\partial
t}$ $\displaystyle=-\boldsymbol{\mathsf{D}}_{2}\boldsymbol{F}_{h}$ (28b)
$\displaystyle\boldsymbol{\mathsf{M}}_{2}\frac{\partial B_{h}}{\partial t}$
$\displaystyle=-\frac{1}{2}\boldsymbol{\mathsf{D}}_{2}\boldsymbol{\mathsf{M}}_{1}^{-1}\boldsymbol{\mathsf{M}}_{1*}(\boldsymbol{F}_{h},b_{h}^{\prime})-\frac{1}{2}\boldsymbol{\mathsf{M}}_{2*}(b_{h}^{\prime},\boldsymbol{\mathsf{M}}_{2}^{-1}\boldsymbol{\mathsf{D}}_{2}\boldsymbol{F}_{h})+\frac{1}{2}\boldsymbol{\mathsf{K}}_{2}(\boldsymbol{F}_{h},\boldsymbol{\mathsf{M}}_{1}^{-1}\boldsymbol{\mathsf{D}}_{2}^{\top}b_{h}^{\prime}).$
(28c)

Recalling the discrete definition of $b_{h}^{\prime}$ (21a) and again assuming
exact spatial integration, left multiplication of (28) by
$[\boldsymbol{F}_{h}^{\top},\Phi_{h}^{f\top},T_{h}^{f\top}]$ leads to the
cancellation of all right hand side terms, ensuring that there is no spurious
generation of energy from the spatial discretisation as

$\displaystyle-\boldsymbol{F}_{h}^{\top}\boldsymbol{\mathsf{C}}_{2}(q_{h},\boldsymbol{F}_{h})+\boldsymbol{F}_{h}^{\top}\boldsymbol{\mathsf{D}}_{2}^{\top}\Phi_{h}^{f}+\frac{1}{4}\boldsymbol{F}_{h}^{\top}\boldsymbol{\mathsf{D}}_{2}^{\top}B_{h}+\frac{1}{4}\boldsymbol{F}_{h}^{\top}\boldsymbol{\mathsf{M}}_{1*}(\boldsymbol{\mathsf{M}}_{1}^{-1}\boldsymbol{\mathsf{D}}_{2}^{\top}h_{h},b_{h}^{\prime})-\frac{1}{4}\boldsymbol{F}_{h}^{\top}\boldsymbol{\mathsf{M}}_{1*}(\boldsymbol{\mathsf{M}}_{1}^{-1}\boldsymbol{\mathsf{D}}_{2}^{\top}b_{h}^{\prime},h_{h})-\Phi_{h}^{f\top}\boldsymbol{\mathsf{D}}_{2}\boldsymbol{F}_{h}$
$\displaystyle-\frac{1}{2}T_{h}^{f\top}\boldsymbol{\mathsf{D}}_{2}\boldsymbol{\mathsf{M}}_{1}^{-1}\boldsymbol{\mathsf{M}}_{1*}(\boldsymbol{F}_{h},b_{h}^{\prime})-\frac{1}{2}T_{h}^{f\top}\boldsymbol{\mathsf{M}}_{2*}(b_{h}^{\prime},\boldsymbol{\mathsf{M}}_{2}^{-1}\boldsymbol{\mathsf{D}}_{2}\boldsymbol{F}_{h})+\frac{1}{2}T_{h}^{f\top}\boldsymbol{\mathsf{K}}_{2}(\boldsymbol{F}_{h},\boldsymbol{\mathsf{M}}_{1}^{-1}\boldsymbol{\mathsf{D}}_{2}^{\top}b_{h}^{\prime})$
$\displaystyle=$ (29a)
$\displaystyle+\frac{1}{4}\boldsymbol{F}_{h}^{\top}\boldsymbol{\mathsf{D}}_{2}^{\top}B_{h}+\frac{1}{4}\boldsymbol{F}_{h}^{\top}\boldsymbol{\mathsf{M}}_{1*}(\boldsymbol{\mathsf{M}}_{1}^{-1}\boldsymbol{\mathsf{D}}_{2}^{\top}h_{h},b_{h}^{\prime})-\frac{1}{4}\boldsymbol{F}_{h}^{\top}\boldsymbol{\mathsf{M}}_{1*}(\boldsymbol{\mathsf{M}}_{1}^{-1}\boldsymbol{\mathsf{D}}_{2}^{\top}b_{h}^{\prime},h_{h})$
$\displaystyle-\frac{1}{4}h_{h}^{\top}\boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{M}}_{1}^{-1}\boldsymbol{\mathsf{M}}_{1*}(\boldsymbol{F}_{h},b_{h}^{\prime})-\frac{1}{4}h_{h}^{\top}\boldsymbol{\mathsf{M}}_{2*}(b_{h}^{\prime},\boldsymbol{\mathsf{M}}_{2}^{-1}\boldsymbol{\mathsf{D}}_{2}\boldsymbol{F}_{h})+\frac{1}{4}h_{h}^{\top}\boldsymbol{\mathsf{K}}_{2}(\boldsymbol{F}_{h},\boldsymbol{\mathsf{M}}_{1}^{-1}\boldsymbol{\mathsf{D}}_{2}^{\top}b_{h}^{\prime})$
$\displaystyle=$ (29b)
$\displaystyle\frac{1}{4}\langle\nabla\cdot\boldsymbol{F}_{h},B_{h}\rangle+\frac{1}{4}\langle\nabla\cdot(b_{h}^{\prime}\boldsymbol{F}_{h}),h_{h}\rangle-\frac{1}{4}\langle\nabla\cdot(h_{h}^{\prime}\boldsymbol{F}_{h}),b_{h}^{\prime}\rangle-$
$\displaystyle\frac{1}{4}\langle
h_{h},\nabla\cdot(b_{h}^{\prime}\boldsymbol{F}_{h})\rangle-\frac{1}{4}\langle
h_{h}b_{h}^{\prime},\nabla\cdot\boldsymbol{F}_{h}\rangle+\frac{1}{4}\langle\nabla\cdot(h_{h}\boldsymbol{F}_{h}),b_{h}^{\prime}\rangle$
$\displaystyle=0.$ (29c)

Defining the entropy for the coupled system as
$\mathcal{S}_{c,h}^{tsw}=\frac{1}{2}\int h_{h}b_{h}^{\prime
2}\mathrm{d}\Omega^{2}$, the variational derivatives may be expressed as
$(\nabla\mathcal{S}_{c,h}^{tsw})^{\top}=[\boldsymbol{0}^{\top},-\frac{1}{2}(b_{h}^{\prime
2})^{\top},b_{h}^{\prime\top}]$. Left multiplication of (28) by this vector
leads to the cancellation of all right hand side terms (again, without
invocation of the chain rule), and hence the conservation of entropy as

$\displaystyle\frac{1}{2}(b_{h}^{\prime
2})^{\top}\boldsymbol{\mathsf{D}}_{2}\boldsymbol{F}_{h}-\frac{1}{2}b_{h}^{\prime\top}\boldsymbol{\mathsf{D}}_{2}\boldsymbol{\mathsf{M}}_{1}^{-1}\boldsymbol{\mathsf{M}}_{1*}(\boldsymbol{F}_{h},b_{h}^{\prime})-\frac{1}{2}b_{h}^{\prime\top}\boldsymbol{\mathsf{M}}_{2*}(b_{h}^{\prime},\boldsymbol{\mathsf{M}}_{2}^{-1}\boldsymbol{\mathsf{D}}_{2}\boldsymbol{F}_{h})+\frac{1}{2}b_{h}^{\prime\top}\boldsymbol{\mathsf{K}}(\boldsymbol{F}_{h},\boldsymbol{\mathsf{M}}_{1}^{-1}\boldsymbol{\mathsf{D}}_{2}^{\top}b_{h}^{\prime})$
$\displaystyle=$ (30a) $\displaystyle\frac{1}{2}\langle b_{h}^{\prime
2},\nabla\cdot\boldsymbol{F}_{h}\rangle-\frac{1}{2}\langle
b_{h}^{\prime},\nabla\cdot(b_{h}^{\prime}\boldsymbol{F}_{h})\rangle-\frac{1}{2}\langle
b_{h}^{\prime
2},\nabla\cdot\boldsymbol{F}_{h}\rangle+\frac{1}{2}\langle\nabla\cdot(b_{h}^{\prime}\boldsymbol{F}_{h}),b_{h}^{\prime}\rangle$
$\displaystyle=0.$ (30b)

If we once again assume exact integration such that
$\int\phi_{h}B_{h}\mathrm{d}\Omega^{2}=\int\phi_{h}h_{h}b_{h}^{\prime}\mathrm{d}\Omega^{2}$
yields a unique solution $\forall\phi_{h}\in\mathbb{V}_{2}$, then this
conserved entropy can be re-expressed in terms of the prognostic variables as
$\frac{1}{2}\int B_{h}^{2}/h_{h}\mathrm{d}\Omega^{2}$. Note that unlike the
mixed formulation, there is no temporal error in the relation between
$b_{h}^{\prime}$ and $B_{h}$, such that this entropy is conserved independent
of the time stepping scheme.

###  2.2 3D compressible Euler equations

Our three dimensional compatible discretisation over the domain $\Omega^{3}$
is a natural extension of the two dimensional discretisation presented in the
previous section, for which we introduce the finite dimensional subspaces
$\mathbb{W}_{0}\subset H^{1}(\Omega^{3})$, $\mathbb{W}_{1}\subset
H(\mathrm{curl},\Omega^{3})$, $\mathbb{W}_{2}\subset
H(\mathrm{div},\Omega^{3})$, $\mathbb{W}_{3}\subset L^{2}(\Omega^{3})$, for
which we have a corresponding de Rham complex [15]

$\mathbb{R}\longrightarrow\mathbb{W}_{0}\stackrel{{\scriptstyle\nabla}}{{\longrightarrow}}\mathbb{W}_{1}\stackrel{{\scriptstyle\nabla\times}}{{\longrightarrow}}\mathbb{W}_{2}\stackrel{{\scriptstyle\nabla\cdot}}{{\longrightarrow}}\mathbb{W}_{3}\longrightarrow
0.$ (31)

This complex discretely preserves the identities $\nabla\times\nabla:=0$ and
$\nabla\cdot\nabla\times:=0$. As for the two dimensional complex in the
previous section, we also recover weak form analogues of the _div_ , _curl_
and _grad_ operators as adjoints of the strong form $\nabla$, $\nabla\times$
and $\nabla\cdot$ operators [19, 18].

In a slight abuse of notation, we will label the basis functions that span
these discrete subspaces as $\alpha_{h}\in\mathbb{W}_{0}$,
$\boldsymbol{\psi}_{h}\in\mathbb{W}_{1}$,
$\boldsymbol{v}_{h}\in\mathbb{W}_{2}$, $\phi_{h}\in\mathbb{W}_{3}$. In a
further abuse of notation we will introduce the operators integrated over the
domain $\Omega^{3}$ as

$\displaystyle\boldsymbol{\mathsf{M}}_{1}$
$\displaystyle:\mathbb{W}_{1}\rightarrow\mathbb{W}_{1}^{\prime}$
$\displaystyle\quad\mathrm{for\
which}\quad\boldsymbol{\mathsf{M}}_{1}\boldsymbol{\psi}_{h}$
$\displaystyle:=\langle\boldsymbol{\psi}_{h},\cdot\rangle,\quad$
$\displaystyle\forall\boldsymbol{\psi}_{h}\in\mathbb{W}_{1},$ (32a)
$\displaystyle\boldsymbol{\mathsf{M}}_{2}$
$\displaystyle:\mathbb{W}_{2}\rightarrow\mathbb{W}_{2}^{\prime}$
$\displaystyle\quad\mathrm{for\
which}\quad\boldsymbol{\mathsf{M}}_{2}\boldsymbol{v}_{h}$
$\displaystyle:=\langle\boldsymbol{v}_{h},\cdot\rangle,\quad$
$\displaystyle\forall\boldsymbol{v}_{h}\in\mathbb{W}_{2},$ (32b)
$\displaystyle\boldsymbol{\mathsf{M}}_{3}$
$\displaystyle:\mathbb{W}_{3}\rightarrow\mathbb{W}_{3}^{\prime}$
$\displaystyle\quad\mathrm{for\
which}\quad\boldsymbol{\mathsf{M}}_{3}\phi_{h}$
$\displaystyle:=\langle\phi_{h},\cdot\rangle,\quad$
$\displaystyle\forall\phi_{h}\in\mathbb{W}_{3},$ (32c)
$\displaystyle\boldsymbol{\mathsf{D}}_{3}$
$\displaystyle:\mathbb{W}_{2}\rightarrow\mathbb{W}_{3}^{\prime}$
$\displaystyle\quad\mathrm{for\
which}\quad\boldsymbol{\mathsf{D}}_{3}\boldsymbol{v}_{h}$
$\displaystyle:=\langle\nabla\cdot\boldsymbol{v}_{h},\cdot\rangle,\quad$
$\displaystyle\forall\boldsymbol{v}_{h}\in\mathbb{W}_{2},$ (32d)
$\displaystyle\boldsymbol{\mathsf{R}}_{3}$
$\displaystyle:\mathbb{W}_{1}\rightarrow\mathbb{W}_{2}^{\prime}$
$\displaystyle\quad\mathrm{for\
which}\quad\boldsymbol{\mathsf{R}}_{3}\boldsymbol{\psi}_{h}$
$\displaystyle:=\langle\nabla\times\boldsymbol{\psi}_{h},\cdot\rangle,\quad$
$\displaystyle\forall\boldsymbol{\psi}_{h}\in\mathbb{W}_{1},$ (32e)
$\displaystyle\boldsymbol{\mathsf{C}}_{3}$
$\displaystyle:\mathbb{W}_{1}\otimes\mathbb{W}_{2}\rightarrow\mathbb{W}_{2}^{\prime}$
$\displaystyle\quad\mathrm{for\
which}\quad\boldsymbol{\mathsf{C}}_{3}(\boldsymbol{\psi}_{h},\boldsymbol{v}_{h})$
$\displaystyle:=\langle\boldsymbol{\psi}_{h}\times\boldsymbol{v}_{h},\cdot\rangle,\quad$
$\displaystyle\forall\boldsymbol{\psi}_{h}\in\mathbb{W}_{1},\boldsymbol{v}_{h}\in\mathbb{W}_{2},$
(32f) $\displaystyle\boldsymbol{\mathsf{M}}_{1*}$
$\displaystyle:\mathbb{W}_{1}\otimes\mathbb{W}_{3}\rightarrow\mathbb{W}_{1}^{\prime}$
$\displaystyle\quad\mathrm{for\
which}\quad\boldsymbol{\mathsf{M}}_{1*}(\boldsymbol{\psi}_{h},\phi_{h})$
$\displaystyle:=\langle\phi_{h}\boldsymbol{\psi}_{h},\cdot\rangle,\quad$
$\displaystyle\forall\boldsymbol{\psi}_{h}\in\mathbb{W}_{1},\phi_{h}\in\mathbb{W}_{3},$
(32g) $\displaystyle\boldsymbol{\mathsf{M}}_{2*}$
$\displaystyle:\mathbb{W}_{2}\otimes\mathbb{W}_{3}\rightarrow\mathbb{W}_{3}^{\prime}$
$\displaystyle\quad\mathrm{for\
which}\quad\boldsymbol{\mathsf{M}}_{2*}(\boldsymbol{v}_{h},\phi_{h})$
$\displaystyle:=\langle\phi_{h}\boldsymbol{v}_{h},\cdot\rangle,\quad$
$\displaystyle\forall\boldsymbol{v}_{h}\in\mathbb{W}_{2},\phi_{h}\in\mathbb{W}_{3},$
(32h) $\displaystyle\boldsymbol{\mathsf{M}}_{3*}$
$\displaystyle:\mathbb{W}_{3}\otimes\mathbb{W}_{3}\rightarrow\mathbb{W}_{3}^{\prime}$
$\displaystyle\quad\mathrm{for\
which}\quad\boldsymbol{\mathsf{M}}_{3*}(\phi_{h},\phi_{k})$
$\displaystyle:=\langle\phi_{h}\phi_{k},\cdot\rangle,\quad$
$\displaystyle\forall\phi_{h},\phi_{k}\in\mathbb{W}_{3},$ (32i)
$\displaystyle\boldsymbol{\mathsf{K}}_{3}$
$\displaystyle:\mathbb{W}_{2}\otimes\mathbb{W}_{2}\rightarrow\mathbb{W}_{3}^{\prime}$
$\displaystyle\quad\mathrm{for\
which}\quad\boldsymbol{\mathsf{K}}_{3}(\boldsymbol{v}_{h},\boldsymbol{v}_{k})$
$\displaystyle:=\langle\boldsymbol{v}_{h}\cdot\boldsymbol{v}_{k},\cdot\rangle,\quad$
$\displaystyle\forall\boldsymbol{v}_{h},\boldsymbol{v}_{k}\in\mathbb{W}_{2},$
(32j) $\displaystyle\boldsymbol{\mathsf{D}}_{3*}$
$\displaystyle:\mathbb{W}_{2}\otimes\mathbb{W}_{3}\rightarrow\mathbb{W}_{3}^{\prime}$
$\displaystyle\quad\mathrm{for\
which}\quad\boldsymbol{\mathsf{D}}_{3*}(\boldsymbol{v}_{h},\phi_{h})$
$\displaystyle:=\langle\nabla\cdot\boldsymbol{v}_{h},\phi_{h}\cdot\rangle,\quad$
$\displaystyle\forall\boldsymbol{v}_{h}\in\mathbb{W}_{2},\phi_{h}\in\mathbb{W}_{3}.$
(32k)

Having defined the various operators, the discrete potential vorticity,
$\boldsymbol{q}_{h}\in\mathbb{W}_{1}$ and variational derivatives
$\boldsymbol{U}_{h}\in\mathbb{W}_{2}$, $\Phi_{h},\Pi_{h}\in\mathbb{W}_{3}$ are
given as

$\displaystyle\boldsymbol{\mathsf{M}}_{1*}(\boldsymbol{q}_{h},h_{h})$
$\displaystyle=-\boldsymbol{\mathsf{R}}_{3}^{\top}\boldsymbol{u}_{h}+\langle\boldsymbol{\psi}_{h},f_{h}\boldsymbol{e}_{z}\rangle$
(33a) $\displaystyle\boldsymbol{\mathsf{M}}_{2}\boldsymbol{U}_{h}$
$\displaystyle=\boldsymbol{\mathsf{M}}_{2*}(\boldsymbol{u}_{h},\rho_{h})$
(33b) $\displaystyle\boldsymbol{\mathsf{M}}_{3}\Phi_{h}$
$\displaystyle=\frac{1}{2}\boldsymbol{\mathsf{K}}_{3}(\boldsymbol{u}_{h},\boldsymbol{u}_{h})+\langle\phi_{h},gz\boldsymbol{e}_{z}\rangle$
(33c) $\displaystyle\boldsymbol{\mathsf{M}}_{3}\Pi_{h}$
$\displaystyle=c_{p}\Big{(}\frac{R}{p_{0}}\Big{)}^{\frac{R}{c_{v}}}\Big{\langle}\phi_{h},\Theta_{h}^{\frac{R}{c_{v}}}\Big{\rangle},$
(33d)

where $\boldsymbol{e}_{z}$ is the unit vector in the vertical direction.
Similarly to (21a), we define the discrete potential temperature,
$\theta_{h}\in\mathbb{W}_{3}$ via

$\boldsymbol{\mathsf{M}}_{3*}(\rho_{h},\theta_{h})=\boldsymbol{\mathsf{M}}_{3}\Theta_{h}.$
(34)

The semi-discrete equivalent of the prognostic equations (17) for
$\boldsymbol{u}_{h}\in\mathbb{W}_{2}$, $\rho_{h},\Theta_{h}\in\mathbb{W}_{3}$
are then expressed as

$\displaystyle\boldsymbol{\mathsf{M}}_{2}\frac{\partial\boldsymbol{u}_{h}}{\partial
t}$
$\displaystyle=-\boldsymbol{\mathsf{C}}_{3}(\boldsymbol{q}_{h},\boldsymbol{U}_{h})+\boldsymbol{\mathsf{D}}_{3}^{\top}\Phi_{h}+\frac{1}{2}\boldsymbol{\mathsf{D}}_{3*}^{\top}(\theta_{h},\Pi_{h})+\frac{1}{2}\boldsymbol{\mathsf{M}}_{2*}(\boldsymbol{\mathsf{M}}_{2}^{-1}\boldsymbol{\mathsf{D}}_{3}^{\top}\Pi_{h},\theta_{h})-\frac{1}{2}\boldsymbol{\mathsf{M}}_{2*}(\boldsymbol{\mathsf{M}}_{2}^{-1}\boldsymbol{\mathsf{D}}_{3}^{\top}\theta_{h},\Pi_{h})$
(35a)
$\displaystyle\boldsymbol{\mathsf{M}}_{3}\frac{\partial\rho_{h}}{\partial t}$
$\displaystyle=-\boldsymbol{\mathsf{D}}_{3}\boldsymbol{U}_{h}$ (35b)
$\displaystyle\boldsymbol{\mathsf{M}}_{3}\frac{\partial\Theta_{h}}{\partial
t}$
$\displaystyle=-\frac{1}{2}\boldsymbol{\mathsf{D}}_{3}\boldsymbol{\mathsf{M}}_{2}^{-1}\boldsymbol{\mathsf{M}}_{2*}(\boldsymbol{U}_{h},\theta_{h})-\frac{1}{2}\boldsymbol{\mathsf{M}}_{3*}(\theta_{h},\boldsymbol{\mathsf{M}}_{3}^{-1}\boldsymbol{\mathsf{D}}_{3}\boldsymbol{U}_{h})+\frac{1}{2}\boldsymbol{\mathsf{K}}_{3}(\boldsymbol{U}_{h},\boldsymbol{\mathsf{M}}_{2}^{-1}\boldsymbol{\mathsf{D}}_{3}^{\top}\theta_{h}).$
(35c)

As for the coupled formulation of the thermal shallow water equations, energy
conservation is assured via the left-multiplication of (35) by
$[\boldsymbol{U}_{h}^{\top},\Phi_{h}^{\top},\Pi_{h}^{\top}]$, such that the
right hand side cancels as

$\displaystyle-\boldsymbol{U}_{h}^{\top}\boldsymbol{\mathsf{C}}_{3}(\boldsymbol{q}_{h},\boldsymbol{U}_{h})+\boldsymbol{U}_{h}^{\top}\boldsymbol{\mathsf{D}}_{3}^{\top}\Phi_{h}+\frac{1}{2}\boldsymbol{U}_{h}^{\top}\boldsymbol{\mathsf{D}}_{3*}^{\top}(\theta_{h},\Pi_{h})+\frac{1}{2}\boldsymbol{U}_{h}^{\top}\boldsymbol{\mathsf{M}}_{2*}(\boldsymbol{\mathsf{M}}_{2}^{-1}\boldsymbol{\mathsf{D}}_{3}^{\top}\Pi_{h},\theta_{h})-\frac{1}{2}\boldsymbol{U}_{h}^{\top}\boldsymbol{\mathsf{M}}_{2*}(\boldsymbol{\mathsf{M}}_{2}^{-1}\boldsymbol{\mathsf{D}}_{3}^{\top}\theta_{h},\Pi_{h})$
$\displaystyle-\Phi_{h}^{\top}\boldsymbol{\mathsf{D}}_{3}\boldsymbol{U}_{h}-\frac{1}{2}\Pi_{h}^{\top}\boldsymbol{\mathsf{D}}_{3}\boldsymbol{\mathsf{M}}_{2}^{-1}\boldsymbol{\mathsf{M}}_{2*}(\boldsymbol{U}_{h},\theta_{h})-\frac{1}{2}\Pi_{h}^{\top}\boldsymbol{\mathsf{M}}_{3*}(\theta_{h},\boldsymbol{\mathsf{M}}_{3}^{-1}\boldsymbol{\mathsf{D}}_{3}\boldsymbol{U}_{h})+\frac{1}{2}\Pi_{h}^{\top}\boldsymbol{\mathsf{K}}_{3}(\boldsymbol{U}_{h},\boldsymbol{\mathsf{M}}_{2}^{-1}\boldsymbol{\mathsf{D}}_{3}^{\top}\theta_{h})$
$\displaystyle=$ (36a)
$\displaystyle\frac{1}{2}\langle\nabla\cdot\boldsymbol{U}_{h},\theta_{h}\Pi_{h}\rangle+\frac{1}{2}\langle\nabla\cdot(\theta_{h}\boldsymbol{U}_{h}),\Pi_{h}\rangle-\frac{1}{2}\langle\nabla\cdot(\Pi_{h}\boldsymbol{U}_{h}),\theta_{h}\rangle$
$\displaystyle-\frac{1}{2}\langle\Pi_{h},\nabla\cdot(\theta_{h}\boldsymbol{U}_{h})\rangle-\frac{1}{2}\langle\theta_{h}\Pi_{h},\nabla\cdot\boldsymbol{U}_{h}\rangle+\frac{1}{2}\langle\nabla\cdot(\Pi_{h}\boldsymbol{U}_{h}),\theta_{h}\rangle$
$\displaystyle=0.$ (36b)

Conservation of entropy is established in an almost identical fashion as for
the discrete coupled thermal shallow water equations in (30a). For the entropy
$\mathcal{S}^{ce}_{h}=\frac{1}{2}\int\rho_{h}\theta_{h}^{2}\mathrm{d}\Omega^{3}$,
left multiplication of (35) by
$(\nabla\mathcal{S}^{ce}_{h})^{\top}=[\boldsymbol{0}^{\top},-\frac{1}{2}\theta_{h}^{2\top},\theta_{h}^{\top}]$
gives

$\displaystyle\frac{1}{2}\theta_{h}^{2\top}\boldsymbol{\mathsf{D}}_{3}\boldsymbol{U}_{h}-\frac{1}{2}\theta_{h}^{\top}\boldsymbol{\mathsf{D}}_{3}\boldsymbol{\mathsf{M}}_{2}^{-1}\boldsymbol{\mathsf{M}}_{2*}(\boldsymbol{U}_{h},\theta_{h})-\frac{1}{2}\theta_{h}^{\top}\boldsymbol{\mathsf{M}}_{3*}(\theta_{h},\boldsymbol{\mathsf{M}}_{3}^{-1}\boldsymbol{\mathsf{D}}_{3}\boldsymbol{U}_{h})+\frac{1}{2}\theta_{h}^{\top}\boldsymbol{\mathsf{K}}_{3}(\boldsymbol{U}_{h},\boldsymbol{\mathsf{M}}_{2}^{-1}\boldsymbol{\mathsf{D}}_{3}^{\top}\theta_{h})$
$\displaystyle=$ (37a)
$\displaystyle\frac{1}{2}\langle\theta_{h}^{2},\nabla\cdot\boldsymbol{U}_{h}\rangle-\frac{1}{2}\langle\theta_{h},\nabla\cdot(\theta_{h}\boldsymbol{U}_{h})\rangle-\frac{1}{2}\langle\theta_{h}^{2},\nabla\cdot\boldsymbol{U}_{h}\rangle+\frac{1}{2}\langle\nabla\cdot(\theta_{h}\boldsymbol{U}_{h}),\theta_{h}\rangle$
$\displaystyle=0.$ (37b)

As for the coupled formulation of the thermal shallow water equations, if we
assume exact spatial integration, then
$\int\phi_{h}\rho_{h}\theta_{h}\mathrm{d}\Omega^{3}=\int\phi_{h}\Theta_{h}\mathrm{d}\Omega^{3}$
yields a unique solution $\forall\phi_{h}\in\mathbb{W}_{3}$, such that this
entropy can be expressed in terms of the prognostic variables.

##  3 Results

###  3.1 Thermal shallow water equations: thermogeostrophic balance

In order to first demonstrate the correctness of the entropy conserving
construction, the thermal shallow water equations are solved in both mixed
(22) and coupled (28) form using a cubed sphere discretisation [20] for an
analytic test case for thermogeostropic balance. This represents an extension
of a similar planar test case [9] extended to the sphere [10] for a steady
solution of the form

$\displaystyle u(\lambda,\chi)$ $\displaystyle=u_{0}\cos(\chi)$ (38a)
$\displaystyle v(\lambda,\chi)$ $\displaystyle=0$ (38b) $\displaystyle
h(\lambda,\chi)$
$\displaystyle=h_{0}-\frac{u_{0}}{g}\Big{(}r_{e}\omega_{e}+\frac{1}{2}u_{0}\Big{)}\sin^{2}(\chi)$
(38c) $\displaystyle b(\lambda,\chi)$
$\displaystyle=g\Bigg{(}1+0.05\Bigg{(}\frac{h_{0}}{h}\Bigg{)}^{2}\Bigg{)},$
(38d)

where $u_{0}=2\pi r_{e}/T$, with $r_{e}=6371220.0m$ being the earth’s radius
and $T=12\ days$, $h_{0}=2.94\times 10^{4}/g$, $\omega_{e}=7.292\times
10^{-5}s^{-1}$ being the earth’s angular frequency and $(\lambda,\chi)$ being
the longitudinal and latitudinal coordinates respectively in the physical
domain $\Omega^{2}$.

As for all test cases presented here, we use third order elements and six
point Gauss-Lobatto-Legendre quadrature, which is exact for polynomials of
degree nine and lower. However since our Jacobian involves trigonometic
functions [21, 20] which cannot be integrated exactly, we do not preserve the
properties of exact spatial integration in our results.

Figure 1: $L^{2}$ errors (left), energy conservation error (center) and
entropy conservation error (right) convergence at day 5 of the
thermogeostrophic test case using the mixed and coupled formulations of the
thermal shallow water equations.

Figure 1 shows the convergence of errors in the $L^{2}$ norm for $h_{h}$,
$\boldsymbol{u}_{h}$ and $B_{h}$ for both the mixed (22) and coupled (28)
formulations (in both cases using a stiffly-stable third order Runge-Kutta
time integrator [22]) after 5 days of integration, using $6\times 8^{2}$,
$6\times 16^{2}$ and $6\times 32^{2}$ third order elements on the cubed sphere
(and corresponding time steps of $\Delta t=120s$, $\Delta t=60s$ and $\Delta
t=30s$). In all cases the errors converge with the anticipated third power of
the average grid spacing, $\Delta x$, with the exception of the velocities,
$\boldsymbol{u}_{h}$, which exhibit a steeper rate of convergence due to the
higher polynomial degree in the normal direction for the basis functions
spanning the $\mathbb{V}_{1}\subset H(\mathrm{div},\Omega^{2})$ subspace [20].

This figure also shows the energy and entropy conservation errors for both
formulations. In terms of energy, this converges at a higher rate for the
mixed formulation, with the temporal errors being associated with the inexact
representation of the temporal chain rule (24). As mentioned previously,
energy can be conserved in time also if the temporal chain rule is preserved
exactly using a semi-implicit time integrator with the variational derivatives
integrated exactly in time [16, 9, 17, 18]. For the entropy, this is conserved
at close to machine precision for the coupled scheme, whereas the mixed scheme
involves a (convergent) temporal error due to the independent time stepping of
both $b_{h}$ and $B_{h}$.

###  3.2 Thermal shallow water equations: shear flow instability

In order to study the conservation properties of the two different energy and
entropy conserving thermal shallow water formulations in a more dynamically
rich setting we present results for a modified shear flow test instability
case [23]. This configuration is extended from its original formulation for
the rotating shallow water equations with the addition of a perturbed buoyancy
field of the form

$b(\lambda,\chi)=g\Big{(}1-0.1\cos(\chi)e^{-(3\lambda)^{2}-(15(\pi/4-\chi))^{2}}\Big{)},$
(39)

which is similar to the Gaussian perturbation to the otherwise geostrophically
balanced depth field used to trigger the shear flow instability in the
original shallow water configuration. For both the mixed and coupled
formulations the model was configured with $6\times 32^{2}$ third order
elements and a time step of $\Delta t=30s$, and run for 20 days so as to
ensure model stability for a mature turbulent state.

No energy or entropy damping is applied in either case, and the only
stabilisation present is applied to the potential enstrophy, via the
anticipated potential vorticity method [24] (using an upwinding coefficient of
$\tau=\frac{\Delta t}{2}$). Since this preserves the skew symmetric structure
of the discrete system, and leaves the buoyancy gradient and transport terms
un-altered, it has no effect on the energy or entropy conservation, and acts
solely to smooth the vorticity field.

Figure 2 shows the conservation errors for the shear flow test case for the
mass, vorticity and depth weighted buoyancy. Both the mixed and coupled
formulations exhibit a machine precision loss of mass at each time step, as
has been observed previously for third order Runge-Kutta time integration
[25]. The un-normalised vorticity conservation errors are consistent with
those observed for previous compatible finite element discretisations of the
rotating shallow water equations [20, 26] for the shear flow instability test
case [23]. While the depth weighted buoyancy conservation errors are of
comparable magnitude to the mass conservation errors for the mixed
formulation, these are $\mathcal{O}(100)$ times larger for the coupled
formulation. The energy and entropy conservation errors for the two
formulations are presented in Fig. 3. In contrast to the depth weighted
buoyancy conservation errors, the energy and entropy conservation errors are
in both cases $\mathcal{O}(100)$ times smaller for the coupled formation.

Figure 2: Conservation errors for the mixed and coupled formulations of the
thermal shallow water equations for the mass (left), un-normalised vorticity
(center) and buoyancy (right).

Figure 3: Conservation errors for the mixed and coupled formulations of the
thermal shallow water equations for the energy (left) and entropy (right).

The northern-hemisphere vorticity and buoyancy fields are shown at day 6 for
both schemes in Figs. 4 and 5. The results are similar for both schemes. Since
no energy or entropy damping is applied, both develop a large amount of grid
scale noise due to nonlinear aliasing. However since both energy and entropy
are conserved by the spatial discretisation, this noise does not trigger any
instability in the model, and both are able to run stably for the duration of
the simulation.

Figure 4: Vorticity field (northern hemisphere) at day 6 for the thermal
shallow water equations using the mixed (left) and coupled (right)
formulations.

Figure 5: Buoyancy field (northern hemisphere) at day 6 for the thermal
shallow water equations using the mixed (left) and coupled (right)
formulations.

###  3.3 3D compressible Euler equations: baroclinic instability

We validate the entropy and energy conserving formulation for the 3D
compressible Euler equations (35) for a standard baroclinic instability test
case on the sphere [27], using $6\times 24^{2}$ elements of third order in the
horizontal dimensions and first order in the vertical dimension, with 30
vertical levels and a time step of $\Delta t=75s$ for 12 days. The model uses
a second order horizontally explicit/vertically implicit time splitting that
exactly balances all energy exchanges [18], and an exact energy conserving
implicit solver for the vertical [17], with the energy and entropy conserving
formulation (35) applied in the horizontal only. In the vertical, the model
uses a flux form transport representation for the density weighted potential
temperature, and a corresponding pressure gradient formulation which conserves
energy but not entropy [17]. As part of this formulation the potential
temperature is projected onto the vertical component of the $\mathbb{W}_{2}$
space, so as to reduce errors associated with spurious buoyancy modes [28].
The model uses a biharmonic viscosity in the horizontal for the momentum
equation [21] and a Rayleigh damping term on the top three levels in the
vertical momentum equation in order to suppress vertical oscillations
associated with the hydrostatic adjustment at the beginning of the simulation.

The model is compared to a near identical formulation [17, 18], for which the
horizontal density weighted potential temperature is transported in flux form
only, which necessitates the inclusion of an additional biharmonic viscosity
term for the flux form density weighted potential temperature equation in
order to ensure long term stability. In the present case neither the flux form
nor the coupled form is run with biharmonic viscosity on the density weighted
potential temperature equation so as to demonstrate the loss of entropy
stability for the original flux formulation. While the new coupled entropy
conserving formulation is able to run stably for the duration of the
simulation with no sign of instability, the original flux formulation goes
unstable due to unbounded entropy conservation errors after 11 days of
simulation time.

Figure 6 shows the evolution of the kinetic (vertical and horizontal),
potential and internal energies (difference from initial values) for the flux
form and coupled formulations. Note the difference in scale for the vertical
kinetic energy evolution (right hand side axis on the first figure), for which
the energy changes are $\mathcal{O}(10^{4})$ times smaller than for the other
energy components. While the horizontal kinetic, potential and internal energy
changes are negligible between the two formulations, there is a sharp uptick
in the vertical kinetic energy for the flux form shortly before the simulation
goes unstable. This is not the case for the coupled entropy conserving
formulation, which remains stable, despite accurately capturing the evolution
of the vertical kinetic energy due to the baroclinic process.

The instability of the flux formulation is preceeded by a steady growth in
entropy, which is not observed for the entropy conserving coupled formulation.
This is despite the fact that the coupled formulation also exhibits a slower
decay in energy, presumably since there is less grid scale noise to be
smoothed out by the biharmonic viscosity on the momentum equation.

Figure 6: Temporal evolution of the energy components for the 3D compressible
Euler equations for the flux and coupled form potential temperature
formulations (left), and the total energy (center) and entropy (right)
conservation errors.

The lowest level Exner pressure is presented at day 11 for both the flux form
and coupled form simulations, as well as for the coupled form at day 12 (after
the flux form has gone unstable) in Fig. 7. While the onset of instability for
the flux form is visible from the somewhat noisy structure of the Exner
pressure contours, this is not the case for the coupled form, which remain
smooth and coherent after 12 days of simulation.

Figure 7: Surface level Exner pressure for the original (flux) formulation and
entropy conserving (coupled) formulation at day 11 (left and center), and
entropy conserving formulation at day 12.

The onset of instability is also observed for the flux form in the vertical
component of the vorticity field, presented at day 11 at a vertical height of
$z\approx 1.5km$ in Fig. 8. Again, no such noise is observed for the coupled
form vorticity field, either at day 11 or 12.

Figure 8: Vertical vorticity component at $z\approx 1.5km$ for the original
(flux) formulation and entropy conserving (coupled) formulation at day 11
(left and center), and entropy conserving formulation at day 12.

##  4 Conclusion

This article presents two alternative energy and entropy conserving spatial
formulations for thermal atmospheric systems in vector invariant form using
compatible finite elements. Based on the insight that the chain rule is not
satisfied for thermodynamic variables on discontinuous function spaces, this
requirement is negated for the conservation of entropy via the averaging of
flux form and material form transport of thermodynamic variables (buoyancy for
the thermal shallow water equations and potential temperature for the 3D
compressible Euler equations). Experiments for the thermal shallow water
equations show that the coupled entropy conserving formulation (an extension
of a recent entropy conserving discontinuous Galerkin formulation [10])
exhibits superior conservation properties to an alternative mixed
material/flux form transport formulation. Corresponding representations of
pressure gradient terms are presented which preserve the skew-symmetric
property, and hence energy conservation, with respect to the thermodynamic
transport terms.

Experiments for the 3D compressible Euler equations show that in contrast to a
non-entropy-conserving potential temperature formulation, the new entropy
conserving coupled formulation is able to run stably without the need for a
dissipative term on the thermodynamic transport.

## References

  * [1] T. C. Fisher, M. H. Carpenter, High-order entropy stable finite difference schemes for nonlinear conservation laws: Finite domains, J. Comp. Phys. 252 (2013) 518–557.
  * [2] M. H. Carpenter, T. C. Fisher, E. J. Nielsen, S. H. Frankel, Entropy Stable Spectral Collocation Schemes for the Navier–Stokes Equations: Discontinuous Interfaces, SIAM J. of Sci. Comput. 36 (2014)
  * [3] G. J. Gassner, A. R. Winters, D. A. Kopriva, A well balanced and entropy conservative discontinuous Galerkin spectral element method for the shallow water equations, Applied Mathematics and Computation 272 (2016) 291–308.
  * [4] J. Chan, On discretely entropy conservative and entropy stable discontinuous Galerkin methods, J. Comp. Phys. 362 (2018) 346–374.
  * [5] M. Waruszewski, J. E. Kozdon, L. C. Wilcox, T. H. Gibson, F. X. Giraldo, Entropy stable discontinuous Galerkin methods for balance laws in non-conservative form: Applications to the Euler equations with gravity, J. Comp. Phys. 468 (2022) 111507.
  * [6] A. Kurganov, Y. Liu, V. Zeitlin. A well-balanced central-upwind scheme for the thermal rotating shallow water equations, J. Comp. Phys. 411 (2020) 109414.
  * [7] 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.
  * [8] 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.
  * [9] C. Eldred, T. Dubos, E. Kritsikis, A quasi-Hamiltonain discretization of the thermal shallow water equations, J. Comp. Phys. 379 (2019) 1–31.
  * [10] K. Ricardo, K. Duru, D. Lee, An entropy stable discontinuous Galerkin method for the vector invariant thermal shallow water equations, (2023) _in prep_.
  * [11] K. Ricardo, D. Lee, K. Duru, Conservation and stability in a discontinuous Galerkin method for the vector invariant spherical shallow water equations, (2023) arXiv:2303.17120.
  * [12] J. Thuburn, Use of the Gibbs thermodynamic potential to express the equation of state in atmospheric models, Q. J. R. Meteorol. Soc. 143 (2017) 1185–1196.
  * [13] P. J. Dellar, Common Hamiltonian structure of the shallow water equations with horizontal temperature gradients and magnetic fields, Phys. Fluids 15 (2003) 292–297.
  * [14] E. S. Warneford, P. J. Dellar, The quasi-geostrophic theory of the thermal shallow water equations, J. Fluid Mech. 723 (2013) 374–403
  * [15] D. N. Arnold, R. S. Falk, R. Winther. Finite element exterior calculus, homological techniques, and applications. Acta Numerica. 15 (2006) 1–155.
  * [16] W. Bauer, C. J. Cotter, Energy-enstrophy conserving compatible finite element schemes for the rotating shallow water equations with slip boundary conditions, J. Comp. Phys. 373 (2018) 171–187.
  * [17] D. Lee, An energetically balanced, quasi-Newton integrator for non-hydrostatic vertical atmospheric dynamics, J. Comp. Phys. 429 (2021) 109988.
  * [18] D. Lee, A Palha, Exact spatial and temporal balance of energy exchanges within a horizontally explicit/vertically implicit non-hydrostatic atmosphere, J. Comp. Phys. 440 (2021) 110432.
  * [19] A. Natale, J. Shipton, C. J. Cotter, Compatible finite element spaces for geophysical fluid dynamics, Dyn. Stat. Climate Sys. 1 (2016) 1–31.
  * [20] 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.
  * [21] O. Guba, M. A. Taylor, P. A. Ullrich, J. R. Overfelt, M. N. Levy, The spectral element method (SEM) on variable-resolution grids: evaluating grid sensitivity and resolution-aware numerical viscosity, Geosci. Model. Dev. 7 (2014) 2803–2816.
  * [22] C. W. Shu, S. Osher, Efficient implementation of essentially nonoscillatory shock-capturing schemes, J. Comp. Phys. 77 (1988) 439–471.
  * [23] J. Galewsky, R. K. Scott, L. M. Polvani, An initial-value problem for testing numerical models of the global shallow-water equations, Tellus A., 56 (2004) 429–440.
  * [24] 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.
  * [25] 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.
  * [26] D. Lee, A. Martin, C. Bladwell, S. Badia, A comparison of variational upwinding schemes for geophysical fluids, and their application to potential enstrophy conserving discretisations in space and time, arXiv:2203.04629 (2022).
  * [27] 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.
  * [28] 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.