characterised in Holm et al. (1998) as Euler-Poicaré reduction by symmetry. In
that work, this framework was extended to quantities that are advected by the
flow, including densities solving the continuity equation, tracers solving the
advection equation, etc., leading to the derivation of the full family of
geophysical fluid dynamics models and beyond: quasigeostrophic approximations
(Holm and Zeitlin 1998), complex fluids (Holm 2002, Gay-Balmaz and Ratiu
2009), vertical slice models (Cotter and Holm 2013a), pseudocompressible and
anelastic approximations (Cotter and Holm 2014), the Boussinesq-alpha model
(Badin et al. 2018), models with moisture and irreversible processes (Gay-
Balmaz 2019). The derivation of conservation laws associated with particle
relabelling symmetries through Noether’s theorem applied to this framework was
presented in Cotter and Holm (2013b), Cotter and Cullen (2019).

###  6.2 Koopman representation: continuous theory

The difficulty with adapting Hamilton’s Principle as described above to
discretisations is that there does not exist a finite dimensional subspace of
flow maps that closes appropriately under composition. A solution to this,
proposed by (Pavlov et al. 2011), is to use the Koopman representation of flow
maps. We shall only briefly describe their approach here, but will provide
more detail about the extension to compatible finite elements shortly.

In the Koopman framework, flow maps $\chi$ in the group
$\operatorname{Diff}(\Omega)$ of diffeomorphisms on $\Omega$ are represented
by elements of $GL(L^{2}(\Omega))$, the invertible linear maps from
$L^{2}(\Omega)$ to $L^{2}(\Omega)$. In particular, flow maps $\chi$
represented by linear maps $\rho_{\chi}\in GL(L^{2}(\Omega))$, defined by
$\rho_{\chi}\cdot a\equiv\rho_{\chi}(a)=a\circ\chi^{-1}$. Of course, there are
many maps in $GL(L^{2}(\Omega))$ that can not be written this way, and so
there is not an isomorphism. In fact, the map defines a subgroup of
$GL(L^{2}(\Omega))$, which we call $G(L^{2}(\Omega))$, which is isomorphic to
the group $\operatorname{Diff}(\Omega)$.

To specialise to incompressible flows, we use $GL_{0}(L^{2}(\Omega))$, defined
as

$\displaystyle GL_{0}(L^{2}(\Omega))=\Bigg{\\{}\rho\in
GL(L^{2}(\Omega)):(\rho\cdot a,\rho\cdot b)_{\Omega}=(a,b)_{\Omega},$
$\displaystyle\quad\forall a,b,\in L^{2}(\Omega)\mbox{ and }\rho\cdot
c=c,\forall c\in\mathbb{R}\Bigg{\\}}.$ (187)

The $\rho_{\chi}$ representation defines an isomorphism of a subgroup of
$GL_{0}(L^{2}(\Omega))$, which we call $G_{0}(L^{2}(\Omega))$, to the group of
volume preserving diffeomorphisms
$\operatorname{Diff}_{\operatorname{vol}}(\Omega)$. The group
$G_{0}(L^{2}(\Omega))$ was approximated in numerical discretisations of the
incompressible Euler equations in Pavlov et al. (2011); the extension to the
full diffeomorphism group was used in the extension to compressible models of
geophysical fluid dynamics in Desbrun et al. (2014), Bauer and Gay-Balmaz
(2017), Brecht et al. (2019), Bauer and Gay-Balmaz (2019).

The Koopman representation also provides an isomorphism between the Lie
algebra $\mathfrak{X}(\Omega)$, corresponding to vector fields on $\Omega$,
and a subspace of the Lie algebra $\mathfrak{gl}(L^{2}(\Omega))$, which we
call $\mathfrak{g}(L^{2}(\Omega))$. For a one parameter family
$\chi_{s}\in\operatorname{Diff}(L^{2}(\Omega))$ of maps with $\chi_{0}=Id$ and
$\frac{\operatorname{d}}{\operatorname{d}s}|_{s=0}\chi_{s}=v$,

$\frac{\operatorname{d}}{\operatorname{d}s}|_{s=0}\rho_{\chi_{s}}a=\frac{\operatorname{d}}{\operatorname{d}s}|_{s=0}(a\circ\chi_{s}^{-1}):=-L_{v}a=v\cdot\nabla
a,$ (188)

where $L_{v}$ is called the _Lie derivative_. Corresponding definitions follow
for the restrictions
$\operatorname{Diff}(\Omega)\to\operatorname{Diff}_{\operatorname{vol}}(\Omega)$
and $GL(L^{2}(\Omega))$ to $GL_{0}(L^{2}(\Omega))$.

###  6.3 Previous work on variational discretisations using Koopman
representation

Pavlov et al. (2011) used the Koopman representation to derive a discrete
variational principle that considered finite subspaces $GL(V)$, where $V$ is
the space of cellwise constant functions defined on a mesh. The subset
$G_{h}\subset GL(V)$ approximating the flow maps was identified as being
generated by maps that only allow instantaneous fluxes between neighbouring
cells (thus approximating the diffeomorphism property). These fluxes can be
described within the framework of Discrete Exterior Calculus (DEC) (Hirani
2003). The difficulty with trying to discretise this structure is that it is
not possible to find a subset $G_{h}$ that can be generated by a finite
subspace of $\mathfrak{gl}(V)$. This was addressed by introducing a
_nonholonomic constraint_ on the time derivative of the flow map (and
corresponding Koopman representative of $GL(V)$). Remarkably, Hamilton’s
principle under these constraints is still reducible (i.e., the $GL(V)$
elements can be eliminated in favour of their generating vector fields),
leading to a discretisation that can be solved entirely in terms of Eulerian
velocities, yielding a spatial discretisation that corresponds to a known
marker-and-cell scheme when a structured grid of square cells is used. This
spatial discretisation is then combined with quadrature approximation in time
to produce a fully discrete variational integrator for fluid dynamics.

The principle goal of the variational discretisation is to obtain schemes with
discrete conservation laws. Before time discretisation, the discrete action is
invariant under time translations, leading to conservation of energy. Time
discretisation breaks this symmetry, but there is the potential to apply
backward error analysis to obtain exponentially accurate conservation of a
modified energy. Whilst numerical results from these variational schemes do
exhibit long time approximate energy conservation as would be expected from
this backward error analysis, it has not yet been adapted to the type of
nonholonomic constraints occurring in this framework. The variational scheme
also contains an echo of the circulation theorem. Instead of considering
loops, we consider currents. These are objects dual to velocities, defined by
the duality pairing

$c[v]=\int_{\Gamma}v\cdot\operatorname{d}x,\quad\forall
v\in\mathfrak{X}(\Omega),$ (189)

for some curve $\Gamma$. Currents are transported by a flow map $\chi$ via
$c\mapsto((\nabla\chi)c)\circ\chi^{-1}$, allowing Kelvin’s circulation theorem
to be reformulated (in the case of incompressible Euler equations) as

$\frac{\operatorname{d}}{\operatorname{d}t}\left(((\nabla\chi)c)\circ\chi^{-1}\right)[u]=0.$
(190)

There is a discrete analogue of this formula when currents are approximated by
objects dual to velocity fields on the discrete grid. These discrete currents
can be considered to be loops that have been smoothed out over a finite area.
This can be used to define a discrete circulation that is conserved along the
solution, including the discrete time solution. However, since the
approximation of a discrete current to a continuous current corresponding to
an advected loop gets noisier as time progresses, it is not clear how or
whether these discrete conservation laws constrain the discrete fluid dynamics
in the same way that circulation conservation does in the continuous case
(which occurs through the link to Casimirs on the Hamiltonian side). This is
another important open question about this framework.

###  6.4 Compatible finite element variational discretisation of
incompressible flow

Inspired by the observation of Dmitry Pavlov that one could extend this
framework to other discretisation methods by simply selecting a discrete space
for velocity fields and a discrete approximation of their Lie group action on
scalar fields, and noting the links between DEC and FEEC, Natale and Cotter
(2018) developed such an extension to compatible finite element methods for
the incompressible Euler case. In this case, we use the space of velocity
fields $\mathring{W}_{h}^{r}$ defined by

$\mathring{W}_{h}^{r}=\left\\{u\in W_{h}^{r}:\nabla\cdot u=0,\,u\cdot
n|_{\partial\Omega}=0\right\\},$ (191)

where $W_{h}^{r}$ is a degree $r$ $BDM$ or $RT$ space (they both have the same
divergence free subspace so the distinction is not important here). For the
discrete Koopman representation of the flows generated by these velocity
fields, we select $V^{s}_{h}$, the space of discontinuous piecewise
polynomials of degree $s$. Following the discrete Lie derivative framework
(the Eulerian version in particular) of Heumann and Hiptmair (2011), given
$u\in\mathring{W}_{h}^{r}$, we define the discrete advection operator
$X_{u}\in\mathfrak{gl}(V_{h}^{s})$ defined by

$\left\langle X_{u}a,b\right\rangle=-\left\langle
a,\nabla\cdot(ub)\right\rangle+\left\llangle\left[\\!\\!\left[au\right]\\!\\!\right],\\{b\\}\right\rrangle_{\Gamma},\quad\forall
a,b\in V_{h}^{s},$ (192)

where $\\{b\\}=(b^{+}+b^{-})/2$. This is a centred approximation of the
advection operator on $V_{h}^{s}$, meaning that $X_{u}$ is linear in $u$ (an
upwinded approximation would break this).

To develop the nonholonomic constraint that enforces a dynamics that converges
to fluid motion by diffeomorphism, we seek a Koopman representation of a time
dependent flow map $\hat{\chi}(\cdot,t)\in G(V_{h}^{s})$ that transports
advected tracers according to $a=\hat{\chi}\cdot a_{0}$. If we require that
all such advected tracers satisfy the equation,

$\frac{\operatorname{d}}{\operatorname{d}t}a+X_{u}a_{0}=0,$ (193)

for some time dependent velocity field $u\in\mathring{W}^{r}_{h}$, then we
obtain the _nonholonomic constraint_ on $\hat{\chi}$ that

$\frac{\operatorname{d}}{\operatorname{d}t}{\hat{\chi}}_{h}\cdot
a+X_{u}\hat{\chi}\cdot a=0,\quad\forall a\in V_{h},$ (194)

for some $X_{u}\in\mathfrak{gl}(V_{h}^{s})$, i.e.,

$\frac{\operatorname{d}}{\operatorname{d}t}\hat{\chi}+X_{u}\hat{\chi}=0.$
(195)

This constraint describes the approximation $\hat{\chi}\in G_{h}$ to the
Koopman representation $\hat{\chi}\in G(L^{2}(\Omega))$ of the flow map
$\chi$. It is nonholonomic because it constrains the time derivative of
$\hat{\chi}$, not $\hat{\chi}$ itself, and this constraint cannot be
integrated to obtain such a constraint $F(\hat{\chi})=0$. This is because the
subspace $S_{h}^{r}(V_{h}^{s})\subset\mathfrak{gl}(V_{h}^{s})$ defined by the
image of the map $u\in\mathring{W}_{h}^{r}\mapsto X_{u}$ is not closed under
Lie brackets. In other words $[X_{u},X_{v}]$ is not guaranteed to be in
$S_{h}^{r}(V_{h}^{s})$ for all $X_{u},X_{v}\in S_{h}^{r}(V_{h}^{s})$.

Natale and Cotter (2018) proved that if $r\geq s$, then the map between $u\in
W^{r}_{h}$ and $X_{u}\in S_{h}^{r}(V_{h}^{s})\subset\mathfrak{gl}(V_{h}^{s})$
is an isomorphism. Gawlik and Gay-Balmaz (2021b) took this further, by
considering the extension of $X_{u}$ to the whole of
$H(\operatorname{div})\cap L^{p}(\Omega)^{n}$ (with some $p>2$; this
technicality ensures the existence of traces on individual facets), using the
same formula (192). They consider the space
$\hat{S}_{h}(V_{h}^{s})\subset\mathfrak{gl}(V_{h}^{s})$, defined by

$\hat{S}_{h}(V_{h}^{s})=\\{X_{u}:u\in H(\operatorname{div})\cap
L^{p}(\Omega)^{n}\\},$ (196)

and proved that $\hat{S}_{h}(V_{h}^{s})$ is isomorphic to $RT^{2r}$ (via the
isomorphism $u\mapsto X_{u}$). This shows that subspaces of $RT^{2r}$ are a
necessary choice for $W^{r}_{h}$. When discontinuous finite element spaces are
chosen for $V$, the compatible finite element framework is thus a necessity
rather than a choice.

To continue the derivation of the discrete incompressible Euler equations, we
need to form the Lagrangian, which is just the kinetic energy. Defined in
terms of the original flow map, the undiscretised Lagrangian is

$\int_{\Omega}\frac{1}{2}|u|^{2}\operatorname{d}x.$ (197)

In the finite element framework, we need to write this Lagrangian as a
functional of $X_{u}\in{S}_{h}(V_{h}^{s})$, but we need to be able to recover
$u$ from $X_{u}$ to substitute into the kinetic energy. As proposed by Dmitry
Pavlov, Natale and Cotter (2018) obtained an approximation to this Lagrangian
by applying $X_{u}$ to each of the Cartesian coordinates, i.e.

$u=\overline{X_{u}}=\sum_{i=1}^{N}X_{u}(x_{i})e_{i},$ (198)

where $e_{i}$ is the unit vector in the direction of increasing coordinate
$x_{i}$.444In fact, this was used to prove the isomorphism between $X_{u}$ and
$u$. This presentation assumes a Cartesian metric, and suitable changes need
to be made for solution of the equations on manifolds such as the surface of a
sphere.

To properly define Hamilton’s principle, we need to define the Lagrangian on
the whole tangent bundle $TG_{h}(V)$ of $G_{h}(V)$ (consisting of pairs
$(\frac{\partial}{\partial t}{\hat{\chi}},\hat{\chi})$ with $\hat{\chi}\in
G_{h}(V)$). Thus we extend the definition of the “overbar map” (198) to the
whole of $\mathfrak{gl}(V_{h})$,

$\overline{\xi}=\sum_{i=1}^{N}e_{i}\xi\cdot
x_{i},\quad\forall\xi\in\mathfrak{gl}(V_{h}),$ (199)

and write the Lagrangian $L_{h}:TG_{h}(V)\to\mathbb{R}$, as

$L_{h}=\int_{\Omega}\frac{1}{2}\left\|\overline{\left(\frac{\partial}{\partial
t}{\hat{\chi}}\circ\hat{\chi}^{-1}\right)}\right\|^{2}\operatorname{d}x.$
(200)

We can write this as a reduced Lagrangian on $\mathfrak{gl}(V_{h}^{s})$,

$\ell_{h}[X]=\int_{\Omega}\frac{1}{2}\left\|\overline{X}\right\|^{2}\operatorname{d}x.$
(201)

The Lagrange-D’Alembert principle (Hamilton’s principle with nonholonomic
constraints) seeks $\hat{\chi}$ with $\frac{\partial}{\partial t}{\hat{\chi}}$
satisfying (195) such that

$\delta\int_{0}^{T}L_{h}(\frac{\partial}{\partial
t}{\hat{\chi}},\hat{\chi})\operatorname{d}t=0,$ (202)

for all variations $\delta\hat{\chi}$ satisfying

$\delta\hat{\chi}+X_{w}\hat{\chi}=0,$ (203)

for some time dependent $w\in\mathring{W}_{h}^{r}$. If we have a Lagrangian
(such as (200)) that is reducible by right action, i.e.

$L_{h}(\frac{\partial}{\partial
t}{\hat{\chi}},\hat{\chi}^{-1})=\ell_{h}(\frac{\partial}{\partial
t}\hat{\chi}\circ\hat{\chi}^{-1}),$ (204)

then we can perform Euler-Poincarè reduction, taking care with the
nonholonomic constraint. This proceeds much as in the unapproximated case
described above, as follows. If we have $g$ such that

$\frac{\partial}{\partial
t}\hat{\chi}+X\hat{\chi}=0,\,\delta\hat{\chi}+Y\hat{\chi}=0,$ (205)

for $X,Y\in\mathfrak{gl}(V_{h}^{s})$, then calculations identical to the above
lead to

$(\delta X)\hat{\chi}+XY\hat{\chi}-\frac{\partial}{\partial
t}Y\hat{\chi}-YX\hat{\chi}=0,$ (206)

i.e.

$\delta X=\frac{\partial}{\partial t}Y+[X,Y],$ (207)

where $[X,Y]=XY-YX$ is the usual commutator for linear operators. Thus the
Lagrange-D’Alembert principle can be reduced to the corresponding reduced
D’Alembert principle, as follows. Find $X\in\mathfrak{gl}(V_{h}^{s})$ subject
to the constraint

$X=X_{u},$ (208)

for some $u\in\mathring{W}_{h}^{r}$, such that

$\delta\int_{0}^{T}\ell_{h}(X)\operatorname{d}t=0,$ (209)

for variations of the form (207) for all time dependent
$Y\in\mathfrak{gl}(V_{h}^{s})$, subject to the constraint $Y=X_{w}$ for some
time dependent $w\in\mathring{W}_{h}^{r}$. The nonclosure of
$S^{r}_{h}(V_{h}^{s})$ under Lie brackets together with the appearance of the
Lie bracket in (207) is the reason why we defined a Lagrangian $TG_{h}(V)$ and
not just for $\hat{\chi}$ satisfying the constraint.555In fact, it is only
necessary to define the reduced Lagrangian $\ell_{h}$ for $X=X_{u}$ and
$X=[X_{u},X_{v}]$ for $u,v\in\mathring{W}_{h}^{r}$. This is how the problem
was approached for the discrete exterior calculus formulation Pavlov et al.
(2011). However, the finite element framework and the bar map makes it easy
enough to extend to the whole of $\mathfrak{gl}(V_{h}^{s})$. From this reduced
principle, we can derive the Euler-Poincarè-D’Alembert equation,

$\left\langle\frac{\operatorname{d}}{\operatorname{d}t}\frac{\delta
l_{h}}{\delta X},Y\right\rangle+\left\langle\frac{\delta l_{h}}{\delta
X},[X,Y]\right\rangle=0,$ (210)

for all $Y$ satisfying (207), where

$\left\langle\frac{\delta l_{h}}{\delta X},Y\right\rangle:=\delta
l_{h}[X;\delta X]=\lim_{\epsilon\to 0}\frac{L[X+\epsilon\delta
X]-L[X]}{\epsilon}.$ (211)

For our reduced Lagrangian (201), we have

$\left\langle\frac{\delta\ell_{h}}{\delta X},\delta
X\right\rangle=\delta\ell_{h}[X;\delta
X]=\left\langle\delta\bar{X},\bar{X}\right\rangle=\left\langle\bar{\delta
X},\bar{X}\right\rangle,\quad\forall\delta X\in\mathfrak{gl}(V_{h}^{s}).$
(212)

This means that solving (210) is equivalent to finding $A_{h}\in
S^{r}_{h}(V_{h})$ such that

$\left\langle\frac{\operatorname{d}}{\operatorname{d}t}\bar{A}_{h},\bar{B}_{h}\right\rangle_{\Omega}+\left\langle\bar{A}_{h},\overline{[A_{h},B_{h}]}\right\rangle_{\Omega}=0,$
(213)

for all $B_{h}\in S^{r}_{h}(V_{h})$, where $[\cdot,\cdot]$ is the commutator
bracket for linear operators. Natale and Cotter (2018) then showed that (210)
is equivalent to finding $u\in\mathring{W}_{s}$ such that

$\left\langle\frac{\partial}{\partial
t}{u},v\right\rangle_{\Omega}+\left\langle
X_{u}u,v\right\rangle=0,\quad\forall v\in\mathring{W}_{h}^{r},$ (214)

where $X_{u}:\mathring{W}_{h}^{r}\to\mathring{W}_{h}^{r}$ is defined as

$\left\langle X_{u}a,b\right\rangle=\left\langle a,\nabla_{h}\times(u\times
b)-u\nabla\cdot
b\right\rangle+\left\llangle\\{u\\},\left[\\!\\!\left[n\times(u\times
b)\right]\\!\\!\right]\right\rrangle_{\Gamma}.$ (215)

Surprisingly, given all of the complexity in the formulation, this takes the
form of a conventional finite element approximation without ever needing to
calculate $X_{u}$. Some further manipulation shows that this discretisation is
in fact identical to the centred flux discretisation described in Guzmán et
al. (2017), which emerged around the same time (but without the variational
derivation).

###  6.5 Compatible finite element discretisation for compressible fluids

The framework was extended to compressible fluid equations in Gawlik and Gay-
Balmaz (2021b), which applied the programme of Desbrun et al. (2014), Bauer
and Gay-Balmaz (2019) to the compatible finite element case. This involves the
introduction of advected quantities such as temperature (which satisfies a
scalar advection equation) and density (which satisfies a continuity
equation). In the unapproximated equations, this enables us to treat
Lagrangians of the form

$L(\frac{\partial\chi}{\partial t},\chi)=\ell(u,a_{1},a_{2},\ldots,a_{n}),$
(216)

where $a_{i}$ are advected quantities satisfying

$\frac{\partial}{\partial t}a_{i}+\mathcal{L}_{u}a_{i}=0,$ (217)

where $\mathcal{L}_{u}$ is a Lie derivative of an appropriate type e.g. for
advected scalars $a$,

$\mathcal{L}_{u}a=u\cdot\nabla a,$ (218)

as before, whilst for advected densities $D$,

$\mathcal{L}_{u}D=\nabla\cdot(uD).$ (219)

In that case, Hamilton’s principle leads to the Euler-Poincaré equation with
advected quantities (Holm et al. 1998),

$\frac{\partial\frac{\delta l}{\delta u}}{\partial t}+u\cdot\nabla\frac{\delta
l}{\delta u}+(\nabla u)^{T}\frac{\delta l}{\delta
u}=\sum_{i}a_{i}\diamond\frac{\delta l}{\delta a_{i}},$ (220)

where the diamond operator $\diamond$ is defined by

$\left\langle a_{i}\diamond\frac{\delta l}{\delta
a_{i}},w\right\rangle=-\left\langle\mathcal{L}_{w}a_{i},\frac{\delta l}{\delta
a_{i}}\right\rangle,$ (221)

for all vector fields $w$. This allows for the relaxation to arbitrary
diffeomorphisms instead of volume preserving ones, and enables the variational
derivation of the full range of compressible fluid models arising in
geophysical fluid dynamics and beyond. In particular, the incompressible Euler
equation can be recovered by introducing a Lagrange multiplier (the pressure)
to enforce constant density $D$.

In the Koopman operator framework, discrete advected densities $D\in
V_{h}^{s}$ are treated by defining their transport equation as being dual to
that of scalar functions $f\in V_{h}^{s}$ i.e. if $f=\hat{\chi}f_{0}\in
V_{h}^{s}$.

$\left\langle D,f\right\rangle=\left\langle D_{0},f_{0}\right\rangle.$ (222)

Hence,

$\displaystyle 0$ $\displaystyle=\frac{\partial}{\partial t}\left\langle
D,f\right\rangle=\left\langle\frac{\partial}{\partial
t}D,f\right\rangle+\left\langle D,\frac{\partial}{\partial t}f\right\rangle,$
(223) $\displaystyle=\left\langle\frac{\partial}{\partial
t}D,f\right\rangle-\left\langle D,Xf\right\rangle.$ (224)

Therefore, we conclude that

$\left\langle\frac{\partial}{\partial t}D,\phi\right\rangle-\left\langle
X^{*}D,\phi\right\rangle=0,\quad\forall\phi\in V_{h}^{s},$ (225)

where $X=(\frac{\partial}{\partial t}g)\circ
g^{-1}\in\mathfrak{gl}(V_{h}^{s})$. When $X$ satisfies the nonholonomic
constraint $X=X_{u}$ for some $u\in{W}_{h}^{r}$, this becomes

$\displaystyle\left\langle\frac{\partial}{\partial t}D,\phi\right\rangle$
$\displaystyle=\left\langle X_{u}^{*}D,\phi\right\rangle=\left\langle
D,X_{u}\phi\right\rangle,$ (226)
$\displaystyle=-\left\langle\nabla_{h}\cdot(uD)\phi\right\rangle+\left\llangle\left[\\!\\!\left[Du\right]\\!\\!\right],\\{\phi\\}\right\rrangle_{\Gamma},$
(227) $\displaystyle=\left\langle
u\cdot\nabla_{h}\phi,D\right\rangle+\int_{\Gamma}u\cdot\Bigg{(}\frac{1}{2}(D^{+}n^{+}+D^{-}n^{-})(\phi^{+}+\phi^{-})$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad-n^{+}D^{+}\phi^{+}-n^{-}D^{-}\phi^{-}\Bigg{)}\operatorname{d}S,$
(228) $\displaystyle=\left\langle
u\cdot\nabla_{h}\phi,D\right\rangle+\int_{\Gamma}u\cdot\Bigg{(}\frac{1}{2}(-D^{+}n^{+}\phi^{+}+D^{+}n^{+}\phi^{-}+$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad
D^{-}n^{-}\phi^{+}-D^{-}n^{-}\phi^{-}\Bigg{)}\operatorname{d}S,$ (229)
$\displaystyle=\left\langle
u\cdot\nabla_{h}\phi,D\right\rangle+\int_{\Gamma}u\cdot\Bigg{(}\frac{1}{2}(-D^{+}n^{+}\phi^{+}-D^{+}n^{-}\phi^{-}-$
(230) $\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad
D^{-}n^{+}\phi^{+}-D^{-}n^{-}\phi^{-}\Bigg{)}\operatorname{d}S,$ (231)
$\displaystyle=\left\langle
u\cdot\nabla_{h}\phi,D\right\rangle-\left\llangle\left[\\!\\!\left[u\phi\right]\\!\\!\right],\\{D\\}\right\rrangle_{\Gamma},$
(232)

which is the standard Discontinuous Galerkin centred flux scheme for the
continuity equation for advected densities.

Similarly, by considering $\delta\left\langle D,f\right\rangle=0$, we obtain

$\delta D-Y^{*}D=0,$ (233)

where $Y$ is such that $\delta\hat{\chi}+Y\hat{\chi}=0$.

We can again use the overbar map to build Lagrangians with advected quantities
defined on the whole of $GL(V_{h}^{s},V_{h}^{s})$. For example, the shallow
water equations (with flat topography) has the (reduced) Lagrangian

$\ell=\int_{\Omega}D\frac{|u|^{2}}{2}-\frac{gD^{2}}{2}\operatorname{d}x.$
(234)

The discrete Lagrangian can then be written as

$\displaystyle L(\frac{\partial}{\partial
t}\hat{\chi},\hat{\chi})=\int_{\Omega}(\hat{\chi}^{*}D_{0})\frac{|\overline{\frac{\partial}{\partial
t}\hat{\chi}\circ\hat{\chi}^{-1}}|^{2}}{2}-\frac{g(\hat{\chi}^{*}D_{0})^{2}}{2}\operatorname{d}x.$
(235)

This has discrete reduced Lagrangian

$\ell_{h}[D,X]=\int_{\Omega}D\frac{|\overline{X}|^{2}}{2}-\frac{gD^{2}}{2}\operatorname{d}x.$
(236)

The reduced Hamilton’s principle with nonholonomic constraints becomes

$\delta\int_{0}^{T}\ell_{h}(X,D)\operatorname{d}t=0,$ (237)

with variations

$\delta X=\frac{\partial}{\partial t}Y+[X,Y],\quad\delta D=Y^{*},$ (238)

and constraints $X=X_{u}$ for some $u\in{W}_{h}^{r}$, and $Y=Y_{w}$ for all
$w\in{W}_{h}^{r}$.

This gives the equation

$\left\langle\frac{\operatorname{d}}{\operatorname{d}t}\frac{\delta l}{\delta
X_{u}},X_{w}\right\rangle-\left\langle\frac{\delta l}{\delta
X_{u}},[X,X_{w}]\right\rangle+\left\langle\frac{\delta l}{\delta
X)_{u}},X_{w}^{*}D\right\rangle,\quad\forall w\in{W}_{h}^{r},$ (239)

which can be written as

$\frac{\operatorname{d}}{\operatorname{d}t}\frac{\delta l}{\delta
X_{u}}+\operatorname{ad}^{*}_{X_{u}}\frac{\delta l}{\delta X_{u}}-\frac{\delta
l}{\delta D}\diamond D\in(S^{r}_{h})^{*},$ (240)

where

$\left\langle\operatorname{ad}^{*}_{X}Y,Z\right\rangle=-\left\langle
Y,[X,Z]\right\rangle,\,-\left\langle X\diamond Y,Z\right\rangle=\left\langle
X,Y^{*}Z\right\rangle,$ (241)

and $(S^{r}_{h})^{*}$ is the dual space to $S^{r}_{h}$ in
$\mathfrak{gl}(V_{h})$. This is the discrete analogue of (220).

###  6.6 Compatible finite element variational discretisations: discrete
conservation laws

The goal of variational discretisations is to derive numerical methods that
provide discrete analogues of conservation laws of the unapproximated
equations. Noether’s theorem can still be applied to nonholonomic variational
principles without change, provided that the constraints are invariant under
the relevant symmetry as well as the Lagrangian. The discrete varational
principles discussed here are invariant under time translation, leading to the
conservation of energy as usual. This is straightforward to check directly in
the incompressible case of Natale and Cotter (2018) since the energy is

$E=\left\langle\bar{A}_{h},\bar{A}_{h}\right\rangle,$ (242)

so the energy equation is obtained in (210) by taking
$\bar{B}_{h}=\bar{A}_{h}$ and using antisymmetry of the bracket.

Regarding Kelvin’s circulation theorem, Natale and Cotter (2018) demonstrated
the same “echo” in (190) for the compatible finite element case, which also
emerges from Gawlik and Gay-Balmaz (2021b) in the case where the Lagrangian
depends only on velocity and density, as expected. In the incompressible case,
the conservation law takes the form

$\frac{\operatorname{d}}{\operatorname{d}t}\left\langle
u,\overline{\hat{\chi}X_{c}\hat{\chi}^{-1}}\right\rangle=0,$ (243)

for all time independent $c\in\mathring{W}_{h}^{r}$.

Gawlik and Gay-Balmaz (2021b) further extended the framework by introducing
advected tracers (which can represent potential temperature, and salinity in
the ocean), which are discretised as elements $\theta$ of $V_{h}$ that can be
acted on by velocity via $X_{u}\theta$ as above. This unlocks the variational
discretisation of all of the main variational models of geophysical fluid
dynamics. The framework was further extended in Gawlik and Gay-Balmaz (2021a)
to accommodate advected transported fluxes represented in ${W}^{r}_{h}$,
leading to variational discretisations of magnetohydrodynamics.

###  6.7 Variational time integrators

So far in the this section we have only discussed the variational
discretisation in space, leading to a system that is still continuous in time.
As set out in the original vision of Pavlov et al. (2011), the idea is to also
discretise Hamilton’s principle in time. Variational integrators arising from
time discretisation in Hamilton’s principle have a long history, with the
programme being formally set out in Marsden and West (2001). Variational time
integration using a finite difference discretisation in time was investigated
in Gawlik and Gay-Balmaz (2021b). Since the kinetic energy is a nonlinear
function of both the coordinate in $V_{h}$ and its rate of change, any
variational integrator will result in a system that requires the solution of
an implicit nonlinear system (using Newton’s method, for example). This
contrasts with the case of classical mechanics, where the Lagrangian
$L(z,\frac{\partial}{\partial t}{z})$ splits into a kinetic energy depending
only on $\frac{\partial}{\partial t}{z}$ and a potential energy depending only
on $z$. This split makes explicit variational integrators possible in that
case. Following Pavlov et al. (2011), Gawlik and Gay-Balmaz (2021b) made the
choice of replacing $A(t)=\frac{\partial}{\partial
t}{\hat{\chi}}\hat{\chi}^{-1}$ with
$A_{k}=\tau^{-1}(\hat{\chi}_{k+1}\hat{\chi}^{-1}_{k})/\Delta t$, where
$\tau:\mathfrak{gl}(V_{h})\to G_{(}V_{h})$ is the Cayley transform

$\tau(A)=\left(I-\frac{A}{2}\right)^{-1}\left(I+\frac{A}{2}\right),$ (244)

with other possibilities for $\tau$ being discussed in Bou-Rabee and Marsden
(2009). After this replacement, the variational integrator is derived by
finding the stationary point of the resulting discrete action principle
depending on $\hat{\chi}_{k+1}$ and $\hat{\chi}_{k}$. It was found through
numerical experiments that the resulting scheme is only conditional stable,
requiring condition $\delta t<Ch$. This is disappointing given that intensive
computation is required to advance the solution by a small step. The nature of
this stepsize requirement is an open problem in the area, as is the question
of whether a variational integrator can be found that allows larger timesteps.
Natale and Cotter (2018) used the implicit midpoint rule to discretise the
semidiscrete variational discretisation. The implicit midpoint rule is not
known to be a variational integrator for such systems, but does have the
property that it preserves any quadratic invariants of the time continuous
system, which includes the energy in the incompressible case.

##  7 Almost-Poisson brackets

An alternative, but related, route to structure preserving discretisations is
found via Poisson bracket formulations.

Poisson brackets are bilinear, skew symmetric maps that take pairs of
functionals on some space where the solutions of the PDE reside, which also
satisfy the Jacobi identity,

$\\{A,\\{B,C\\}\\}+\\{B,\\{C,A\\}\\}+\\{C,\\{A,B\\}\\}=0,$ (245)

for all functionals $A,B,C$. In the process of discretisation of Poisson
brackets for fluid dynamics, the Jacobi identity is lost, for similar reasons
that the nonholonomic constraints are required for the variational
discretisations of Section 6. Poisson brackets without the Jacobi identity are
called “almost Poisson brackets” but we shall just use the term Poisson
bracket in the rest of this article for brevity. Poisson brackets can have
special functionals called Casimirs, which make the Poisson bracket vanish,
i.e. $C$ is a Casimir if $\\{F,C\\}=0$ for all functionals $F$. Then, Casimirs
are conserved by the dynamics since

$\frac{\partial}{\partial t}{C}=-\\{C,H\\}=0.$ (246)

The goal of building discretisations using Poisson brackets is that they
automatically conserve the Hamiltonian, and if any Casimirs survive the
discretisation process then they will be conserved as well. We shall discuss
specific examples later.

The aim of building structure preserving discretisations for fluid PDEs using
Poisson brackets was introduced in Morrison (1982), Salmon (1983), although in
fact the energy and enstrophy preserving discretisations in Arakawa (1966) are
the first instance of a Poisson bracket discretisation (but not presented that
way), which was extended to the rotating shallow water equations and beyond in
subsequent work by Arakawa, Sadourny and others (Sadourny et al. 1968,
Sadourny 1972, 1975, Arakawa and Lamb 1977, 1981, Arakawa and Hsu 1990)).

In the 2000s, when attention was focussed much more on triangular and
polygonal grids to provide a more uniform coverage of the sphere, the idea of
using Poisson brackets was revived to produce energy conserving schemes (or to
at least guide the design of practical schemes that are as energy consistent
as possible). This took place in a number of groups (Ringler et al. 2010,
Skamarock et al. 2012, Eldred and Randall 2017, Gassmann and Herzog 2008,
Gassmann 2013, Tort et al. 2015, Dubos et al. 2015). There was also work on
extending Poisson brackets to Nambu brackets in pursuit of additional
conserved quantities in the method (Sommer and Névir 2009).

There have also been some interesting studies about the relevance of
conservation for geophysical models. Dubinkina and Frank (2007) demonstrated
the benefits of using energy-enstrophy conserving schemes to obtain correct
statistical equilibria, and Thuburn et al. (2014b) demonstrated that energy
conservation is important to obtain realistic backscatter in underresolved
simulations of two dimensional turbulence (provided that enstrophy is
dissipated at the small scale). Dubinkina (2018) demonstrated that conserving
both energy and enstrophy is important in the context of data assimilation.
Even when structure preserving discretisations lead to systems of equations
that are challenging to solve efficiently, it is useful to consider how they
relate to more standard discretisations to see where conservation errors are
being committed, and to see when they are likely to be large or small.

In this section, we describe how Poisson brackets can be used to construct
energy (and enstrophy) conserving schemes. This work has been heavily informed
by previous works using finite difference methods, such as those cited above.
We shall start by briefly discussing the 2D incompressible Euler equations,
the the rotating shallow water equations and then vertical slice and three
dimensional models. We initially assume that we are solving the equations on a
closed manifold (the surface of the sphere, or periodic boundary conditions,
for example), and will return to the treatment of boundary conditions later.

###  7.1 Incompressible Euler equations: continuous theory

Returning to Arnold’s variational formulation for incompressible flow in
(182), we can reformulate as

$\left\langle\frac{\partial m}{\partial
t},v\right\rangle+\left\langle\left[\frac{\delta H}{\delta
m},v\right],m\right\rangle=0,\,\forall
v\in\mathfrak{X}_{\operatorname{vol}}(\Omega),$ (247)

where

$H(m)=\left\langle m,u\right\rangle-\ell(u),$ (248)

and we take $m=\frac{\delta l}{\delta
u}\in\mathfrak{X}_{\operatorname{vol}}(\Omega)$, now inverting the
relationship so that $u$ is considered as an operator applied to $m$. Here we
use the variational derivative $\frac{\delta F}{\delta
u}\in\mathfrak{X}_{\operatorname{vol}}(\Omega)$ defined by

$\left\langle\frac{\delta F}{\delta u},v\right\rangle=\lim_{\epsilon\to
0}\frac{1}{\epsilon}\left(F[u+\epsilon v]-F[u]\right).$ (249)

Equation (247) is equivalent to the Poisson formulation

$\frac{\partial}{\partial t}{F}[m]+\left\\{F,H\right\\}=0,$ (250)

for arbitrary functionals
$F:\mathfrak{X}_{\operatorname{vol}}(\Omega)\to\mathbb{R}$, with Poisson
bracket

$\left\\{F,G\right\\}=\int_{\Omega}\left[\frac{\delta G}{\delta
m},\frac{\delta F}{\delta m}\right]\cdot m\operatorname{d}x.$ (251)

Since the Poisson bracket is antisymmetric, this leads immediately to
conservation of the Hamiltonian,

$\frac{\operatorname{d}}{\operatorname{d}t}{H}=-\\{H,H\\}+\frac{\partial
H}{\partial t}=0,$ (252)

provided that $H$ has no explicit dependence on time $t$.

For the 2D incompressible Euler equations, a clear path towards deriving
enstrophy conservation requires modification of the bracket (251) by changing
variables $m\to u$. This produces the equivalent Poisson formulation

$\frac{\partial}{\partial t}{F}[u]+\left\\{F,H\right\\}[u]=0,$ (253)

for all functionals
$F:\mathfrak{X}_{\operatorname{vol}}(\Omega)\to\mathbb{R}$, where

$\\{F,G\\}=\int_{\Omega}\omega\left(\frac{\delta F}{\delta
u}\right)\cdot\frac{\delta G}{\delta u}^{\perp}\operatorname{d}x,$ (254)

$\omega=\nabla^{\perp}\cdot u$, and

$H=\frac{1}{2}\int_{\Omega}|u|^{2}\operatorname{d}x.$ (255)

The derivation of this bracket formulation from (251) is discussed in Morrison
(1982), Marsden and Weinstein (1983). Here, we just directly demonstrate that
it leads to the incompressible Euler equations by computation. First, we
compute $\frac{\delta H}{\delta u}=u$. For a linear functional
$F[u]=\left\langle u,w\right\rangle$ with
$w\in\mathfrak{X}_{\operatorname{vol}}(\Omega)$, (253) becomes

$\left\langle u_{t},w\right\rangle+\left\langle\omega,w\cdot
u^{\perp}\right\rangle=0,\quad\forall
w\in\mathfrak{X}_{\operatorname{vol}}(\Omega),$ (256)

which is a weak formulation of the equation

$u_{t}+\omega u^{\perp}+\nabla P=0,$ (257)

where $P$ is some potential chosen so that
$u_{t}\in\mathfrak{X}_{\operatorname{vol}}(\Omega)$ (because we only test
against functions from $\mathfrak{X}_{\operatorname{vol}}(\Omega)$ in (257),
so the equation is projected into
$\mathfrak{X}_{\operatorname{vol}}(\Omega)$). Writing
$P=\frac{1}{2}|u|^{2}+p$, we obtain

$u_{t}+u^{\perp}\nabla^{\perp}\cdot u+\nabla\frac{1}{2}|u|^{2}+\nabla p=0,$
(258)

which becomes recognisable as the incompressible Euler equations after
recalling the identity

$(u\cdot\nabla)u=\omega u^{\perp}+\frac{1}{2}\nabla|u|^{2}.$ (259)

Returning to the Poisson bracket (254), we find that it has an infinite number
of Casimirs of the form

$C_{n}[u]=\int_{\Omega}\omega^{n}\operatorname{d}x,\quad n=1,2,\ldots$ (260)

To verify that $C_{n}$ is a Casimir, we compute

$\displaystyle\left\langle\frac{\delta C_{n}}{\delta u},v\right\rangle$
$\displaystyle=n\int_{\Omega}\omega^{n-1}\nabla^{\perp}v\operatorname{d}x,\quad\forall
v\in\mathfrak{X}_{\operatorname{vol}}(\Omega),$ (261)
$\displaystyle=\int_{\Omega}\nabla^{\perp}(-n\omega^{n-1})\cdot
v\operatorname{d}x,$ (262)

having integrated by parts (we assume for now that there are now boundaries
and so may ignore the surface term). Hence, we conclude formally that

$\frac{\delta C_{n}}{\delta u}=-n\nabla^{\perp}\omega^{n-1}.$ (263)

Inserting into the Poisson bracket then gives

$\displaystyle\left\\{F,C_{n}\right\\}$
$\displaystyle=-n\int_{\Omega}\omega\frac{\delta F}{\delta
u}\cdot\left(\nabla^{\perp}\omega^{n-1}\right)^{\perp}\operatorname{d}x,$
(264) $\displaystyle=n\int_{\Omega}\omega\frac{\delta F}{\delta
u}\cdot\nabla\omega^{n-1}\operatorname{d}x,$ (265)
$\displaystyle=(n-1)\int_{\Omega}\frac{\delta F}{\delta
u}\cdot\nabla\omega^{n}\operatorname{d}x,$ (266)
$\displaystyle=-(n-1)\int_{\Omega}\omega^{n}\underbrace{\nabla\cdot\frac{\delta
F}{\delta u}}_{=0}\operatorname{d}x=0,$ (267)

having integrated by parts again, for any functional $F$ on
$\mathfrak{X}_{\operatorname{vol}}(\Omega)$. We have $\nabla\cdot\frac{\delta
F}{\delta u}=0$ since $\frac{\delta F}{\delta
u}\in\mathfrak{X}_{\operatorname{vol}}(\Omega)$. Hence, $C_{n}$ is a conserved
quantity for the Poisson dynamics from any Hamiltonian. In particular, $C_{1}$
is the _total vorticity_ , whilst $C_{2}$ is the _enstrophy_ , both of which
provide strong constraints on two dimensional incompressible turbulence.

###  7.2 Incompressible Euler equations: compatible finite element
discretisation

To make our compatible finite element discretisation, we restrict $u$ and $w$
to the divergence-free subspace $\zeta_{h}$ of some chosen H(div) finite
element space $\mathbb{V}_{h}^{1}$ (such as Raviart-Thomas or Brezzi-Douglas-
Marini on triangles). We have to make a further approximation since
$\nabla^{\perp}\cdot u$ is not defined for H(div) spaces, and so we define
$\omega_{h}\in\mathbb{V}_{h}^{0}$, such that

$\left\langle\gamma,\omega_{h}\right\rangle=-\left\langle\nabla^{\perp}\gamma,u\right\rangle,\quad\forall\gamma\in\mathbb{V}_{h}^{0},$
(268)

i.e. $\omega_{h}$ is defined from the approximated weak curl of $u$, the dual
of the $\nabla^{\perp}$ operator restricted to $\mathbb{V}_{h}^{0}$ and
$\mathbb{V}_{h}^{1}$. Note that this is where we have used the absence of
boundary, otherwise there would be a boundary term causing complications that
we shall discuss later. Having defined $\omega_{h}$, we write the discrete
Poisson bracket as

$\left\\{F,G\right\\}=\int_{\Omega}\omega_{h}\frac{\delta F}{\delta
u}\cdot\frac{\delta G}{\delta u}^{\perp}\operatorname{d}x,$ (269)

where $F$ and $G$ are now functionals on $\mathbb{V}_{h}^{1}$. With the same
Hamiltonian (255) now restricted to $\mathbb{V}_{h}^{1}$, the Poisson
formulation implies the following dynamics for $u\in\zeta_{h}$,

$\left\langle w,u_{t}\right\rangle+\left\langle\omega_{h}w\cdot
u^{\perp}\right\rangle=0,\quad\forall w\in\zeta_{h}.$ (270)

The discrete Hamiltonian is conserved as usual through the antisymmetry of the
bracket and the time-independence of $H$. Concerning Casimirs of the bracket,
we can repeat the earlier calculation for functionals $C_{n}$ computed with
$\omega_{h}$ substituted for $\omega$, i.e.

$C_{n,h}[u]=\int_{\Omega}\omega_{h}^{n}\operatorname{d}x.$ (271)

Following (262), we obtain

$\left\langle\frac{\delta C_{n,h}}{\delta
v}\right\rangle=\int_{\Omega}\nabla^{\perp}(-n\omega^{n-1})\cdot
v\operatorname{d}x,\quad\forall v\in\zeta_{h},$ (272)

but now can only conclude that $\frac{\delta C_{n,h}}{\delta
v}=P_{1}\left(-n\nabla^{\perp}\omega^{n-1}\right)$, where $P_{1}$ is the
$L^{2}$ projection into $\zeta_{h}$, which prevents us from showing that
$C_{n,h}$ is a Casimir for $n>2$. However, when $n=1$, we obtain that
$\frac{\delta C_{1,h}}{\delta v}=0$ (so $C_{1,h}$ is trivially conserved, just
as $C_{1}$ is for the unapproximated case). When $n=2$, we have
$\nabla^{\perp}\omega\in\mathbb{V}_{h}^{1}$ by the embedding property of the
discrete de Rham complex, and then a calculation identical to (264-267) shows
that $\\{C_{2,h},G\\}=0$ for any functional $G$ on $\zeta_{h}$, and hence the
numerical enstrophy $C_{2,h}$ is a Casimir and is conserved for dynamics
generated from any Hamiltonian. Hence, this scheme conserves energy, total
vorticity, and enstrophy.

To make a practical implementation of the scheme, one can follow two
approaches. The main hurdle is that the scheme is defined on $\zeta_{h}$, and
not the whole of $\mathbb{V}_{h}^{1}$. Since $\nabla\cdot$ maps from
$\mathbb{V}_{h}^{1}$ onto $\mathbb{V}_{h}^{2}$, the divergence free subspace
$\zeta_{h}$ is equivalently represented as

$\zeta_{h}=\left\\{u\in\mathbb{V}_{h}^{1}:\int_{\Omega}\phi\nabla\cdot
u\operatorname{d}x=0,\quad\forall\phi\in\mathbb{V}_{h}^{2}\right\\}.$ (273)

Hence, we can equivalently write the following system for
$(\omega,u,P)\in\mathbb{V}_{h}^{0}\times\mathbb{V}_{h}^{1}\times\mathbb{V}_{h}^{2}$
such that

$\displaystyle\left\langle\gamma,\omega\right\rangle-\left\langle\nabla^{\perp}\gamma,u\right\rangle$
$\displaystyle=0,\quad\forall\gamma\in\mathbb{V}_{h}^{0},$ (274)
$\displaystyle\left\langle w,u_{t}\right\rangle+\left\langle w,\omega
u^{\perp}\right\rangle-\left\langle\nabla\cdot w,P\right\rangle$
$\displaystyle=0,\quad\forall w\in\mathbb{V}_{h}^{1},$ (275)
$\displaystyle\left\langle\phi,\nabla\cdot u\right\rangle$
$\displaystyle=0,\quad\forall\phi\in\mathbb{V}_{h}^{2}.$ (276)

It can easily be checked that the solution satisfies $\nabla\cdot u=0$ in
$L^{2}$. Selecting $w\in\zeta_{h}\subset\mathbb{V}_{h}^{1}$ makes the $P$ term
disappear and we recover (270). This formulation builds a bridge to the
shallow water and compressible systems that we shall look at later.

On the other hand, for $u\in\zeta_{h}$ we can directly parameterise
$u=\nabla^{\perp}\psi$ for $\psi\in\mathbb{V}_{h}^{0}$ and choose
$w=\nabla^{\perp}\beta$ in (270) for all $\beta\in\mathbb{V}_{h}^{0}$ and we
obtain

$\displaystyle\left\langle\gamma,\omega\right\rangle+\left\langle\nabla\gamma,\nabla\psi\right\rangle$
$\displaystyle=0,\quad\forall\gamma\in\mathbb{V}_{h}^{0},$ (277)
$\displaystyle\left\langle\nabla\beta,\nabla\psi_{t}\right\rangle+\left\langle\nabla\beta,\omega\nabla^{\perp}\psi\right\rangle$
$\displaystyle=0,\quad\forall\beta\in\mathbb{V}_{h}^{0},$ (278)

which is a discretisation of the incompressible Euler equation in vorticity
streamfunction form,

$\omega_{t}+\nabla\cdot(\omega\nabla^{\perp}\psi)=0,\quad-\nabla^{2}\psi=\omega.$
(279)

This formulation was presented in Liu et al. (2001), where it was analysed in
the viscous case, but the energy conservation of the inviscid equations did
form an important part of the proof. A related scheme was presented in Liu and
Shu (2000), but with the vorticity in a discontinuous space, with appropriate
jump terms defining the fluxes between cells. With an average flux, the energy
and enstrophy are both conserved. If a Lax Friedrichs flux is used, then
energy is still conserved but enstrophy is dissipated (as is appropriate for
long time simulations of cascading 2D vortex dynamics). It is possible to
modify the Poisson bracket so that enstrophy is dissipated in this way, with
the antisymmetric formulation still conserving energy, as we shall see later.
Liu and Shu (2000) proved convergence for both types of fluxes. This
formulation and analysis was extended by Bernsen et al. (2006) to the
quasigeostrophic model of large scale rotating geophysical fluid dynamics,
including the case with islands in the flow (considering the fluid as an
ocean).

To extend the conservation properties to a full discrete method, one can use
the implicit midpoint rule, which conserves all quadratic invariants of the
continuous time system. This includes energy and enstrophy (or just energy
where enstrophy is dissipated) in the incompressible Euler case.

###  7.3 Rotating shallow water equations: continuous theory

Moving on to the rotating shallow water equations, it is tempting to continue
working with the vorticity-streamfunction formulation above. However, this
places a limitation on the possibilities of extension to three dimensional
models, complicates the boundary conditions, and is not preferred by
practitioners since the prognostic variables are not quantities that are
directly measurable. Hence, we must address the challenge of finding a
compatible finite element discretisation of the rotating shallow water
equations using velocity $u$ and layer depth $D$.

We start from the Lie-Poisson formulation of the rotating shallow water
equations given by

$\displaystyle H$
$\displaystyle=\int_{\Omega}\frac{1}{2D}|m|^{2}+gD\left(\frac{D}{2}+b\right)\operatorname{d}x,$
(280) $\displaystyle\\{F,G\\}$ $\displaystyle=\left\langle[\frac{\delta
G}{\delta m},\frac{\delta F}{\delta
m}],m\right\rangle+\left\langle\frac{\delta F}{\delta
D},\nabla\cdot(D\frac{\delta G}{\delta
m})\right\rangle-\left\langle\frac{\delta G}{\delta
D},\nabla\cdot(D\frac{\delta F}{\delta m})\right\rangle,$ (281)

where $m$ is defined by

$\int_{\Omega}m\cdot v\operatorname{d}x=\int_{\Omega}D\left({u}+R\right)\cdot
v\operatorname{d}x,\quad\forall v\in\mathfrak{X}(\Omega),$ (282)

where $\nabla^{\perp}\cdot R=f$, the Coriolis parameter. These equations
emerge from the reduced Hamilton’s principle with advected density $D$ after
applying the Legendre transform as described in Holm et al. (1998). Following
the incompressible case, if we want enstrophy conservation to emerge then we
need to change variables to $(u,D)$, this leads to the equivalent Poisson
bracket formulation

$\displaystyle H$
$\displaystyle=\frac{1}{2}\int_{\Omega}D\|u\|^{2}+gD\left(\frac{D}{2}+b\right)\operatorname{d}x,$
(283) $\displaystyle\\{F,G\\}$ $\displaystyle=\left\langle q,\frac{\delta
F}{\delta u}\cdot\frac{\delta G}{\delta
u}^{\perp}\right\rangle-\left\langle\nabla\cdot\frac{\delta F}{\delta
u},\frac{\delta G}{\delta D}\right\rangle+\left\langle\nabla\cdot\frac{\delta
G}{\delta u},\frac{\delta F}{\delta D}\right\rangle,$ (284)

where $q=\frac{\nabla^{\perp}\cdot u+f}{D}$ is the potential vorticity. After
computing the variational derivatives

$\displaystyle\frac{\delta H}{\delta u}$ $\displaystyle=Du,$ (285)
$\displaystyle\frac{\delta H}{\delta D}$
$\displaystyle=\frac{1}{2}\|u\|^{2}+g(D+b),$ (286)

and substituting into the Poisson bracket equation $F_{t}+\\{F,H\\}$, we
formally obtain

$\displaystyle u_{t}+qDu^{\perp}+\nabla\left(\frac{1}{2}|u|^{2}+g(D+b)\right)$
$\displaystyle=0,$ (287) $\displaystyle D_{t}+\nabla\cdot(uD)$
$\displaystyle=0,$ (288)

which we recognise as the rotating shallow water equations in vector invariant
form. We define functionals $C_{n}$ (which will turn out to be Casimirs) by

$C_{n}[u,D]=\int_{\Omega}Dq^{n}\operatorname{d}x.$ (289)

To compute the variational derivatives, given $\epsilon>0$ and $(v,\phi)\in
H(\operatorname{div})\times L^{2}$, writing $D_{\epsilon}=D+\epsilon\phi$,
$u_{\epsilon}=u+\epsilon v$, we define $q_{\epsilon}\in H^{1}$ such that

$\left\langle\gamma,D_{\epsilon}q_{\epsilon}\right\rangle=-\left\langle\nabla^{\perp}\gamma,u_{\epsilon}\right\rangle+\left\langle\gamma,f\right\rangle,\quad\forall\gamma\in
H^{1},$ (290)

noting that $q_{\epsilon}|_{\epsilon=0}=q$ (after integrating by parts and
using the lack of surface term in a domain without boundary). Then,

$\displaystyle\left\langle\frac{\delta C_{n}}{\delta
v},v\right\rangle+\left\langle\frac{\delta C_{n}}{\delta D},\phi\right\rangle$
$\displaystyle=\lim_{\epsilon\to 0}\frac{1}{\epsilon}\left(C_{n}[u+\epsilon
v,D+\epsilon\phi]-C_{n}[u,D]\right),$ (291)
$\displaystyle=\int_{\Omega}\lim_{\epsilon\to
0}\frac{D_{\epsilon}q_{\epsilon}^{n}-Dq^{n}}{\epsilon}\operatorname{d}x,$
(292)
$\displaystyle=\int_{\Omega}\frac{\operatorname{d}}{\operatorname{d}\epsilon}|_{\epsilon=0}D_{\epsilon}q_{\epsilon}^{n}\operatorname{d}x,$
(293) $\displaystyle=\left\langle
q_{\epsilon}^{n-1}|_{\epsilon=0},\frac{\operatorname{d}}{\operatorname{d}\epsilon}|_{\epsilon=0}D_{\epsilon}q_{\epsilon}\right\rangle$
$\displaystyle\qquad+\left\langle(n-1)q_{\epsilon}^{n-2}|_{\epsilon=0}\frac{\operatorname{d}}{\operatorname{d}\epsilon}|_{\epsilon=0}q_{\epsilon},D_{\epsilon}q_{\epsilon}\right\rangle,$
(294) $\displaystyle=n\left\langle
q_{\epsilon}^{n-1}|_{\epsilon=0},\frac{\operatorname{d}}{\operatorname{d}\epsilon}|_{\epsilon=0}D_{\epsilon}q_{\epsilon}\right\rangle$
$\displaystyle\qquad-\left\langle(n-1)q_{\epsilon}^{n}|_{\epsilon=0},\frac{\operatorname{d}}{\operatorname{d}\epsilon}|_{\epsilon=0}D_{\epsilon})\right\rangle,$
(295)

From (290) we have

$\left\langle\gamma,\frac{\operatorname{d}}{\operatorname{d}\epsilon}|_{\epsilon=0}D_{\epsilon}q_{\epsilon}\right\rangle=-\left\langle\nabla^{\perp}\gamma,v\right\rangle,\quad\forall\gamma\in
H^{1},$ (296)

(correcting a typographic error in Equation 58 of McRae and Cotter (2014)),
and hence we have

$\displaystyle\left\langle\frac{\delta C_{n}}{\delta
v},v\right\rangle+\left\langle\frac{\delta C_{n}}{\delta D},\phi\right\rangle$
$\displaystyle=n\left\langle\nabla^{\perp}q^{n-1},v\right\rangle-\left\langle(n-1)q^{n},\phi\right\rangle,$
(297)

i.e.

$\frac{\delta C_{n}}{\delta v}=n\nabla^{\perp}q^{n-1},\quad\frac{\delta
C_{n}}{\delta D}=-(n-1)q^{n}.$ (298)

Then, inserting into the Poisson bracket (284), we obtain

$\displaystyle\\{C_{n},G\\}$
$\displaystyle=\left\langle\underbrace{nq\nabla^{\perp}q^{n-1}}_{=(n-1)\nabla^{\perp}q^{n}},\frac{\delta
G}{\delta u}^{\perp}\right\rangle$
$\displaystyle\qquad-\left\langle\underbrace{\nabla\cdot
n\nabla^{\perp}q^{n-1}}_{=0},\frac{\delta G}{\delta
D}\right\rangle+\left\langle(n-1)q^{n},\nabla\cdot\frac{\delta G}{\delta
u}\right\rangle,$ (299) $\displaystyle=\left\langle(n-1)\nabla
q^{n},\frac{\delta G}{\delta
u}\right\rangle-\left\langle(n-1)q^{n},\nabla\cdot\frac{\delta G}{\delta
u}\right\rangle=0,$ (300)

where the last line is obtained by integrating by parts, and hence $C_{n}$ is
a Casimir of the bracket (284), and so it is a conserved quantity of the
dynamics for any Hamiltonian. Another, more simple, Casimir is the mass

$M[u,D]=\int_{\Omega}D\operatorname{d}x,$ (301)

with variational derivatives

$\frac{\delta M}{\delta u}=0,\quad\frac{\delta M}{\delta D}=1,$ (302)

and hence

$\\{M,G\\}=\left\langle\nabla\cdot\frac{\delta G}{\delta
u},1\right\rangle=\int_{\Omega}\nabla\cdot\frac{\delta G}{\delta
D}\operatorname{d}x=0,$ (303)

using the divergence theorem (assuming no boundary currently), so mass is also
a Casimir of the bracket and is conserved for any Hamiltonian.

###  7.4 Rotating shallow water equations: Poisson bracket discretisation

To produce a discretisation of this structure, McRae and Cotter (2014) simply
took the structure (283-284), restricted $(u,D)$ to
$\mathbb{V}_{h}^{1}\times\mathbb{V}_{h}^{2}$, and replaced $q$ with the
discrete approximation $q\in\mathbb{V}_{h}^{0}$, with

$\left\langle\gamma,Dq\right\rangle=-\left\langle\nabla^{\perp}\gamma,u\right\rangle+\left\langle\gamma,f\right\rangle,\quad\forall\gamma\in\mathbb{V}_{h}^{0},$
(304)

following the exposition for the discretisation for incompressible Euler
equations above. Then, we obtain

$\displaystyle\left\langle w,\frac{\delta H}{\delta u}-Du\right\rangle$
$\displaystyle=0,\,\forall w\in\mathbb{V}_{h}^{1},$ (305)
$\displaystyle\left\langle\phi,\frac{\delta H}{\delta
D}-\frac{1}{2}|u|^{2}-g(D+b)\right\rangle$
$\displaystyle=0,\quad\forall\phi\in\mathbb{V}_{h}^{2},$ (306)

i.e.,

$\frac{\delta H}{\delta u}=P_{1}(Du),\quad\frac{\delta H}{\delta
D}=P_{2}\left(\frac{1}{2}|u|^{2}+g(D+b)\right),$ (307)

where $P_{1}$ and $P_{2}$ are the $L^{2}$ projections into
$\mathbb{V}_{h}^{1}$ and $\mathbb{V}_{h}^{2}$, respectively. To derive the
equations of motion, we take $F[u,D]=\left\langle
w,u\right\rangle+\left\langle\phi,D\right\rangle$ for
$w,\phi\in\mathbb{V}_{h}^{1}\times\mathbb{V}_{h}^{2}$, so that $\frac{\delta
F}{\delta u}=w$, $\frac{\delta F}{\delta D}=\phi$, and substitute into the
Poisson dynamics to obtain

$\displaystyle\left\langle
w,u_{t}\right\rangle+\left\langle\phi,D_{t}\right\rangle$
$\displaystyle=\frac{\partial}{\partial t}{F}$ (308)
$\displaystyle=-\\{F,H\\},$ (309) $\displaystyle=-\left\langle q,w\cdot
P_{1}(Du)^{\perp}\right\rangle+\left\langle\nabla\cdot
w,P_{2}\left(\frac{1}{2}|u|^{2}+g(D+b)\right)\right\rangle$
$\displaystyle\qquad-\left\langle\nabla\cdot P_{1}(Du),\phi\right\rangle,$
(310) $\displaystyle=-\left\langle q,w\cdot
P_{1}(Du)^{\perp}\right\rangle+\left\langle\nabla\cdot
w,\frac{1}{2}|u|^{2}+gD\right\rangle-\left\langle\nabla\cdot
P_{1}(Du),\phi\right\rangle,$ $\displaystyle\qquad\qquad\quad\forall
w,\phi\in\mathbb{V}_{h}^{1}\times\mathbb{V}_{h}^{2},$ (311)

where we were able to drop the $P_{2}$ since the result of the projection was
in an $L^{2}$ inner product with $\nabla\cdot w\in P_{2}$ (so the discrete de
Rham complex is crucial here). Writing $m=P_{1}(Du)$ we put everything
together as
$(u,D,q,m)\in\mathbb{V}_{h}^{1}\times\mathbb{V}_{h}^{2}\times\mathbb{V}_{h}^{0}\times\mathbb{V}_{h}^{1}$,
such that

$\displaystyle\left\langle w,u_{t}\right\rangle+\left\langle
qw,m^{\perp}\right\rangle-\left\langle\nabla\cdot
w,\frac{1}{2}|u|^{2}+g(D+b)\right\rangle$ $\displaystyle=0,\,\quad\forall
w\in\mathbb{V}_{h}^{1},$ (312)
$\displaystyle\left\langle\phi,D_{t}+\nabla\cdot m\right\rangle$
$\displaystyle=0,\,\forall\phi\in\mathbb{V}_{h}^{2},$ (313)
$\displaystyle\left\langle\gamma,qD\right\rangle+\left\langle\nabla^{\perp}\gamma,u\right\rangle-\left\langle\gamma,f\right\rangle$
$\displaystyle=0,\,\forall\gamma\in\mathbb{V}_{h}^{0},$ (314)
$\displaystyle\left\langle v,m-uD\right\rangle$ $\displaystyle=0,\,\forall
v\in\mathbb{V}_{h}^{1}.$ (315)

This is a set of coupled equations which must be solved together, but since
Equations (314-315) do not contain time derivatives, $q$ and $m$ may be
reconstructed at any time from $u$ and $D$. Since the equations are derived
from a Poisson bracket formulation, we can immediately deduce that they
conserve the Hamiltonian. Concerning the Casimirs, we make the same
calculations for mass $M$ as for the undiscretised case, obtaining

$\displaystyle\frac{\delta M}{\delta u}=0,\quad\frac{\delta M}{\delta
D}=P_{1}(D)=1.$ (316)

For $C_{n}$, we can only make use of (290) when $n=1$ (so that $\gamma=1$) or
$n=2$ (so that $\gamma=q_{\epsilon}$), leading to

$\frac{\delta C_{n}}{\delta
u}=n\left(\nabla^{\perp}q^{n-1}\right),\quad\frac{\delta C_{n}}{\delta
D}=-(n-1)P_{2}\left(q^{n}\right).$ (317)

(300) then follows but only for $n=1,2$,

$\displaystyle\\{C_{n},G\\}$ $\displaystyle=\left\langle
nq\nabla^{\perp}q^{n-1},\frac{\delta G}{\delta u}^{\perp}\right\rangle$
$\displaystyle\qquad-\left\langle
n\underbrace{\nabla\cdot\nabla^{\perp}q^{n-1}}_{=0},\frac{\delta G}{\delta
u}\right\rangle+\left\langle P_{2}((n-1)q^{n}),\nabla\cdot\frac{\delta
G}{\delta u}\right\rangle,$ (318) $\displaystyle=\left\langle
nq\nabla^{\perp}q^{n-1},\frac{\delta G}{\delta
u}^{\perp}\right\rangle+\left\langle(n-1)q^{n},\nabla\cdot\frac{\delta
G}{\delta u}\right\rangle=0,$ (319)

where we may integrate by parts since $\frac{\delta G}{\delta
u}\in\mathbb{V}_{h}^{1}\subset H(\operatorname{div})$ and
$q^{n}\in\mathbb{V}_{h}^{0}\subset H^{1}$. Hence, $C_{n}$ is a Casimir for
$n=1$ (total vorticity) and $n=2$ (enstrophy).

McRae and Cotter (2014) verified these conservation properties in numerical
experiments, and showed second order convergence with $h$ for the scheme with
the $BDFM_{1}$-$DG_{1}$ finite element spaces on triangles ($BDFM_{1}$ is a
slightly more exotic variant which has an intermediate number of degrees of
freedom between $BDM_{1}$ and $RT_{1}$, which results in a 2:1 ratio of
velocity to pressure degrees of freedom).

The Poisson bracket approach has been extended to other compatible spaces with
various motivations. Eldred et al. (2019) used finite element spaces built
around splines to form higher order discrete de Rham complexes. These spaces
have the same degrees of freedom as the lowest order $Q_{1}$-$RT_{0}$-$DG_{0}$
complex on quadrilaterals, and achieve higher order by making use of degrees
of freedom from a patch of neighbouring cells. The advantage is that this
removes the jump in the dispersion relation for gravity waves in higher order
spaces on quadrilaterals, as discussed in Section 3.5. The price to be paid is
that there is increased interelement coupling, and that there are some
technicalities at the boundaries between patches of structured quadrilaterals,
e.g. at the edges and vertices of the cube upon which a cubed sphere mesh is
constructed. Since Hamiltonian and the Poisson brackets are the same, the only
thing that has changed is the finite element spaces, which still satisfy the
discrete de Rham complex so energy and enstrophy conservation follows
directly.

Lee et al. (2018) extended the method to mixed mimetic spectral elements,
which are a variant of mixed elements using spectral element histopolation
functions to construct high order spaces. The usual spectral element technique
of using incomplete quadrature then leads to diagonal mass matrices for the
continuous space $\mathbb{V}_{h}^{0}$ without losing the discrete de Rham
complex property. Again, since the Hamiltonian and Poisson brackets are the
same (excepting some details on quadrature rules, where care must be taken),
and the new finite element spaces still satisfy the discrete de Rham complex,
the energy and enstrophy conservation follows directly. Lee and Palha (2018)
extended these spaces to the surface of the sphere.

###  7.5 Poisson integrators

To extend these conservation properties to a fully discrete method after time
discretisation, we need to look beyond the implicit midpoint rule into the
more general case of Poisson integrators. To make this generalisation, we
write the Poisson bracket as

$\\{F,G\\}=A\left(\frac{\delta F}{\delta z},\frac{\delta G}{\delta
z};z\right),$ (320)

where $z\in W$ comprises the dynamic fields (i.e., $z=(u,D)$ and
$W=\mathbb{V}_{h}^{1}\times\mathbb{V}_{h}^{2}$ for the case of the shallow
water equations). We use this notation to express that Poisson brackets are
bilinear in $(\frac{\delta F}{\delta z},\frac{\delta G}{\delta z})$ but with
possibly arbitrary additional dependence on $z$, which acts as a coefficient.
When there is no explicit dependence on $z$, we obtain linear dynamics, and
the $z$ dependency encodes nonlinear dynamics. From the properties of the
Poisson bracket, $A$ is bilinear and antisymmetric in $\frac{\delta F}{\delta
z}$ and $\frac{\delta G}{\delta z}$. To derive one particular Poisson
integrator, we write $z(s)=z^{n}+s(z^{n}-z^{n+1})$, and seek $z^{n+1}$ such
that

$\int_{0}^{\Delta t}\left(\frac{\partial}{\partial t}{F}[z(s)]+A(\frac{\delta
F}{\delta z}[z(s)],\frac{\delta H}{\delta
z}[z(s)];z^{n+1/2})\right)\beta(s)\operatorname{d}s=0,$ (321)

for all linear functions $\beta(s)$. In other words, we replace the bracket
$\\{F,G\\}$ by

$\\{F,G\\}=A\left(\frac{\delta F}{\delta z},\frac{\delta G}{\delta
z};z^{n+1/2}\right),$ (322)

where $z^{n+1/2}=(z^{n}+z^{n+1})/2$, and project the equation onto linear
dynamics in time. We observe energy conservation since taking $F=H$ and
$\beta=1$ leads to

$\displaystyle\frac{H[z^{n+1}]-H[z^{n}]}{\Delta t}$
$\displaystyle=\int_{0}^{\Delta t}\frac{\partial}{\partial
s}{H}[z(s)]\operatorname{d}s,$ (323) $\displaystyle=-\int_{0}^{\Delta
t}\underbrace{A\left(\frac{\delta H}{\delta z}[z(s)],\frac{\delta H}{\delta
z}[z(s)];z^{n+1/2}\right)}_{=0}\operatorname{d}s=0,$ (324)

by antisymmetry.

This formulation leads to a practical method since taking $F[z]=\left\langle
w,z\right\rangle$ for $w\in W$ and $\beta=1$ gives

$\left\langle w,z^{n+1}-z^{n}\right\rangle+A\left(w,\int_{0}^{\Delta
t}\frac{\delta H}{\delta z}[z(s)]\operatorname{d}s;z^{n+1/2}\right)=0,$ (325)

by linearity in the second argument. This scheme was introduced along with
higher order variants as a larger set of Poisson integrators in Hairer (2010),
Cohen and Hairer (2011). The easiest way to obtain implementable formulae for
the scheme is to choose a quadrature rule for the time integral in (325) such
that the integral is exact. This is possible whenever the Hamiltonian is
polynomial. For example, when the Hamiltonian is quadratic, $\frac{\delta
H}{\delta z}$ is linear, and the midpoint rule

$\int_{0}^{\Delta t}\frac{\delta H}{\delta
z}[z(s)]\operatorname{d}s=\frac{1}{2}\left(\frac{\delta H}{\delta
z}\left[z^{n}+\frac{1}{2}\left(z^{n+1}-z^{n}\right)\right]\right),$ (326)

(evaluated at $s=1/2$ with weight 1), is exact; the scheme is then equivalent
to the implicit midpoint rule.

For the rotating shallow water equation scheme described above, the
Hamiltonian is cubic, so a two point quadrature must be used to compute the
time integral involving the quadratic derivatives of the Hamiltonian exactly.
We obtain the scheme

$\displaystyle\left\langle
w,u^{n+1}\right\rangle+\left\langle\phi,D^{n+1}\right\rangle=$
$\displaystyle\left\langle
w,u^{n}\right\rangle+\left\langle\phi,D^{n}\right\rangle$
$\displaystyle\qquad-\Delta tA\left((w,\phi),\left(\overline{\frac{\delta
H}{\delta u}},\overline{\frac{\delta H}{\delta
D}}\right);(u^{n+1/2},D^{n+1/2})\right),$
$\displaystyle\qquad\forall(w,\phi)\in\mathbb{V}_{h}^{1}\times\mathbb{V}_{h}^{2},$
(327)

where

$\displaystyle\overline{\frac{\delta H}{\delta u}}$
$\displaystyle=P_{1}[m^{n+1/2}],$ (328)
$\displaystyle=P_{1}\frac{1}{3}\left(D^{n}u^{n}+\frac{1}{2}D^{n}u^{n+1}+\frac{1}{2}D^{n+1}u^{n}+D^{n+1}u^{n+1}\right),$
(329) $\displaystyle\overline{\frac{\delta H}{\delta D}}$
$\displaystyle=P_{2}(\pi^{n+1/2}),$ (330)
$\displaystyle=P_{2}\left(\frac{1}{6}\left(|u^{n}|^{2}+u^{n}\cdot
u^{n+1}+|u^{n+1}|^{2}\right)+\frac{g}{2}(D^{n+1}+D^{n})+b\right).$ (331)

This can be implemented as
$(u^{n+1},D^{n+1},q^{n+1/2},\overline{m}^{n+1/2})\in\mathbb{V}_{h}^{1}\times\mathbb{V}_{h}^{2}\times\mathbb{V}_{h}^{0}\times\mathbb{V}_{h}^{1}$
such that

$\displaystyle\left\langle w,u^{n+1}-u^{n}\right\rangle+\Delta t\left\langle
w,q^{n+1/2}(\overline{m}^{n+1/2})^{\perp}\right\rangle$
$\displaystyle\qquad-\Delta t\left\langle\nabla\cdot
w,\pi^{n+1/2}\right\rangle$ $\displaystyle=0,\quad\forall
w\in\mathbb{V}_{h}^{1},$ (332)
$\displaystyle\left\langle\phi,D^{n+1}-D^{n}\right\rangle+\Delta
t\left\langle\phi,\nabla\cdot\overline{m}^{n+1/2}\right\rangle$
$\displaystyle=0,\quad\forall\phi\in\mathbb{V}_{h}^{2},$ (333)
$\displaystyle\left\langle\gamma,D^{n+1/2}q^{n+1/2}\right\rangle+\left\langle\nabla^{\perp}\gamma,u^{n+1/2}\right\rangle-\left\langle\gamma,f\right\rangle$
$\displaystyle=0,\quad\forall\gamma\in\mathbb{V}_{h}^{0},$ (334)
$\displaystyle\left\langle v,\overline{m}^{n+1/2}-m^{n+1/2}\right\rangle$
$\displaystyle=0,\quad\forall v\in\mathbb{V}_{h}^{1},$ (335)

with $m^{n+1/2}$, $\pi^{n+1/2}$ defined as above. This scheme is similar, but
not identical to, the implicit midpoint rule, but results in energy
conservation. Cohen and Hairer (2011) showed that this class of Poisson
integrators also preserves Casimirs that are at most quadratic functions of
the state space variables. This covers mass and total vorticity, but not
enstrophy, which is a nonpolynomial function of $u$ and $D$. Finding Poisson
integrators that preserve enstrophy for this discrete formulation of the
rotating shallow water equations is an open problem.

###  7.6 Enstrophy conservation on domains with boundaries

Something that we have neglected from our discussion so far is the case when
the domain $\Omega$ has an exterior boundary. This is important when extending
these tools to ocean applications (where there are coastlines). When
boundaries are present, we consider the subcomplex

$\begin{CD}\mathring{H}^{1}@>{\nabla^{\perp}}>{}>\mathring{H}(\textrm{div})@>{\nabla\cdot}>{}>L^{2}\\\
@V{}V{\pi_{0}}V@V{}V{\pi_{1}}V@V{}V{\pi_{2}}V\\\
\mathring{\mathbb{V}}_{h}^{0}@>{\nabla^{\perp}}>{}>\mathring{\mathbb{V}}_{h}^{1}@>{\nabla\cdot}>{}>\mathbb{V}_{h}^{2}\\\
\end{CD}$ (336)

where

$\displaystyle\mathring{H}_{1}$ $\displaystyle=\\{\phi\in
H^{1}:\operatorname{tr}_{\partial\Omega}{\phi}=0\\},$ (337)
$\displaystyle\mathring{\mathbb{V}}_{h}^{0}$
$\displaystyle=\\{\phi\in\mathbb{V}_{h}^{0}:\operatorname{tr}_{\partial\Omega}{\phi}=0\\},$
(338) $\displaystyle\mathring{H}(\operatorname{div})$ $\displaystyle=\\{u\in
H(\operatorname{div}):\operatorname{tr}_{\partial\Omega}{u\cdot n}=0\\},$
(339) $\displaystyle\mathring{\mathbb{V}}_{h}^{1}$
$\displaystyle=\\{u\in\mathbb{V}_{h}^{1}:\operatorname{tr}_{\partial\Omega}{u\cdot
n}=0\\},$ (340)

and where $\operatorname{tr}$ is the boundary trace operator. When a boundary
is present, we must modify (304) to incorporate a boundary integral,

$\left\langle\gamma,Dq\right\rangle=-\left\langle\nabla^{\perp}\gamma,u\right\rangle+\left\llangle\gamma,n^{\perp}\cdot
u\right\rrangle+\left\langle\gamma,f\right\rangle,\quad\forall\gamma\in\mathbb{V}_{h}^{0},$
(341)

where $\left\llangle\cdot,\cdot\right\rrangle$ is the usual $L^{2}$ inner
product on $\partial\Omega$. This is necessary because $q$ does not vanish on
the boundary in general. One can then try to proceed using this modified
definition of $q$ in the Poisson bracket, now defined over
$\mathring{\mathbb{V}}_{1}\times\mathbb{V}_{h}^{2}$. The proof of enstrophy
conservation then fails because although
$\nabla^{\perp}q\in\mathbb{V}_{h}^{1}$,
$\nabla^{\perp}q\notin\mathring{\mathbb{V}}_{1}$ in general, because $q$ is
not constant on the boundary. One possible solution is to restrict
$q\in\mathring{\mathbb{V}}_{h}^{0}$, and use Equation (304) with test
functions in $\mathring{\mathbb{V}}_{h}^{0}$, so that
$\nabla^{\perp}q\in\mathring{\mathbb{V}}_{h}^{1}$. This recovers enstrophy
conservation, but at the expense of consistency as it commits a first-order
error in forcing $q$ to be zero on the boundary.

An alternative solution presented in Bauer and Cotter (2018) is to split
$\mathbb{V}_{h}^{0}=\mathring{\mathbb{V}}_{h}^{0}\oplus(\mathring{\mathbb{V}}_{h}^{0})^{\perp}$
where $(\mathring{\mathbb{V}}_{h}^{0})^{\perp}$ is the $L^{2}$-orthogonal
complement of $\mathring{\mathbb{V}}_{h}^{0}$ in $\mathbb{V}_{h}^{0}$. We then
extend the solution space to
$(u,D,Z^{\prime})\in\mathring{\mathbb{V}}_{h}^{1}\times\mathbb{V}_{h}^{2}\times(\mathring{\mathbb{V}}_{h}^{0})^{\perp}$,
and define the bracket

$\displaystyle\\{F,G\\}$ $\displaystyle=\left\langle q,\frac{\delta F}{\delta
u}\cdot\frac{\delta G}{\delta
u}^{\perp}\right\rangle-\left\langle\nabla\cdot\frac{\delta F}{\delta
u},\frac{\delta G}{\delta D}\right\rangle+\left\langle\nabla\cdot\frac{\delta
G}{\delta u},\frac{\delta F}{\delta D}\right\rangle$
$\displaystyle\qquad+\left\langle\nabla\frac{\delta F}{\delta
Z^{\prime}},q\frac{\delta G}{\delta
u}\right\rangle-\left\langle\nabla\frac{\delta G}{\delta
Z^{\prime}},q\frac{\delta F}{\delta u}\right\rangle,$ (342)

with the same Hamiltonian, where $q\in\mathbb{V}_{h}^{0}$ such that

$\displaystyle\left\langle\gamma,qD\right\rangle-\left\langle\gamma,f+\mathring{Z}+Z^{\prime}\right\rangle$
$\displaystyle=0,\forall\gamma\in\mathbb{V}_{h}^{0},$ (343)
$\displaystyle\left\langle\gamma,\mathring{Z}\right\rangle-\left\langle\nabla^{\perp}\gamma,u\right\rangle$
$\displaystyle=0,\quad\forall\gamma\in\mathring{\mathbb{V}}_{h}^{0}.$ (344)

Here, $Z^{\prime}$ represents the contribution to the vorticity from the
boundary, which would normally be given by

$\left\langle\gamma,Z^{\prime}\right\rangle=-\left\langle\nabla^{\perp}\gamma,u\right\rangle+\left\llangle\gamma,n^{\perp}\cdot
u\right\rrangle+\left\langle\gamma,f\right\rangle,\quad\forall\gamma\in(\mathring{\mathbb{V}}_{h}^{0})^{\perp}.$
(345)

However, here we just ensure that $Z^{\prime}$ satisfies this condition
initially, and that at future times $Z^{\prime}$ has its own dynamics
consistent with the potential vorticity equation as we shall see later (and
see 8 for further discussion of this). Since the Hamiltonian does not depend
on $Z^{\prime}$, we have $\frac{\delta H}{\delta u}$ and $\frac{\delta
H}{\delta D}$ as before, and $\frac{\delta H}{\delta Z^{\prime}}=0$. Thus the
$u_{t}$ and $D_{t}$ equations are unchanged, and we have

$\left\langle\gamma,Z^{\prime}_{t}\right\rangle=\left\langle\nabla\gamma,q\frac{\delta
H}{\delta
u}\right\rangle=\left\langle\nabla\gamma,mq\right\rangle,\forall\gamma\in(\mathring{\mathbb{V}}_{h}^{0})^{\perp}.$
(346)

The bracket is antisymmetric, so the Hamiltonian is conserved.

To check whether this has repaired the Casimirs, we recompute the derivatives
of $C_{n}$ for this extended phase space, writing
$D_{\epsilon}=D+\epsilon\phi$ for $\phi\in\mathbb{V}_{h}^{2}$,
$u_{\epsilon}=u+\epsilon v$ for $v\in\mathring{\mathbb{V}}_{h}^{1}$,
$Z^{\prime}_{\epsilon}=Z^{\prime}+\epsilon\gamma^{\prime}$ for
$\gamma^{\prime}\in(\mathring{\mathbb{V}}_{h}^{0})^{\perp}$, and defining
$q_{\epsilon}\in\mathbb{V}_{h}^{0}$ and
$\mathring{Z}\in\mathring{\mathbb{V}}_{h}^{0}$ such that

$\displaystyle\left\langle\gamma,q_{\epsilon}D_{\epsilon}\right\rangle-\left\langle\gamma,f+Z^{\prime}_{\epsilon}+\mathring{Z}_{\epsilon}\right\rangle$
$\displaystyle=0,\forall\gamma\in\mathbb{V}_{h}^{0},$ (347)
$\displaystyle\left\langle\gamma,\mathring{Z}_{\epsilon}\right\rangle-\left\langle\nabla^{\perp}\gamma,u_{\epsilon}\right\rangle$
$\displaystyle=0,\quad\forall\gamma\in\mathring{\mathbb{V}}_{h}^{0}.$ (348)

Then, for $n\in(1,2)$,

$\displaystyle\delta C_{n}$ $\displaystyle=\left\langle\frac{\delta
C_{n}}{\delta u},v\right\rangle+\left\langle\frac{\delta C_{n}}{\delta
D},\phi\right\rangle+\left\langle\frac{\delta C_{n}}{\delta
Z^{\prime}},\gamma^{\prime}\right\rangle$ (349)
$\displaystyle=\lim_{\epsilon\to 0}\frac{1}{\epsilon}\left(C_{n}[u+\epsilon
v,D+\epsilon\phi,Z^{\prime}+\epsilon\gamma^{\prime}]-C_{n}[u,D]\right),$ (350)
$\displaystyle=\int_{\Omega}\frac{\operatorname{d}}{\operatorname{d}\epsilon}|_{\epsilon=0}D_{\epsilon}q_{\epsilon}^{n}\operatorname{d}x,$
(351) $\displaystyle=\left\langle
q_{\epsilon}^{n-1}|_{\epsilon=0},\frac{\operatorname{d}}{\operatorname{d}\epsilon}|_{\epsilon=0}D_{\epsilon}q_{\epsilon}\right\rangle$
$\displaystyle\qquad+\left\langle(n-1)q_{\epsilon}^{n-2}|_{\epsilon=0}\frac{\operatorname{d}}{\operatorname{d}\epsilon}|_{\epsilon=0}q_{\epsilon},D_{\epsilon}q_{\epsilon}\right\rangle,$
(352) $\displaystyle=n\left\langle
q_{\epsilon}^{n-1}|_{\epsilon=0},\frac{\operatorname{d}}{\operatorname{d}\epsilon}|_{\epsilon=0}D_{\epsilon}q_{\epsilon}\right\rangle$
$\displaystyle\qquad-\left\langle(n-1)q_{\epsilon}^{n}|_{\epsilon=0},\frac{\operatorname{d}}{\operatorname{d}\epsilon}|_{\epsilon=0}D_{\epsilon}\right\rangle,$
(353) $\displaystyle=n\left\langle
q^{n-1},\gamma^{\prime}+\frac{\operatorname{d}}{\operatorname{d}\epsilon}|_{\epsilon=0}\mathring{Z}_{\epsilon}\right\rangle-\left\langle(n-1)q^{n},\phi\right\rangle,$
(354)
$\displaystyle=n\left\langle\mathring{P}_{0}^{\perp}q^{n-1},\gamma^{\prime}\right\rangle+n\left\langle\mathring{P}_{0}q^{n-1},\frac{\operatorname{d}}{\operatorname{d}\epsilon}|_{\epsilon=0}\mathring{Z}_{\epsilon}\right\rangle-\left\langle(n-1)q^{n},\phi\right\rangle,$
(355)
$\displaystyle=n\left\langle\mathring{P}_{0}^{\perp}q^{n-1},\gamma^{\prime}\right\rangle+n\left\langle\nabla^{\perp}\mathring{P}_{0}q^{n-1},v\right\rangle-\left\langle(n-1)q^{n},\phi\right\rangle,$
(356)

where $\mathring{P}_{0}^{\perp}$ is the $L^{2}$ projection onto the orthogonal
subspace to $\mathring{\mathbb{V}}_{h}^{0}$. Hence,

$\frac{\delta C_{n}}{\delta
u}=n\nabla^{\perp}\mathring{P}_{0}q^{n-1},\quad\frac{\delta C_{n}}{\delta
Z^{\prime}}=n\mathring{P}_{0}^{\perp}q^{n-1},\quad\frac{\delta C_{n}}{\delta
D}=-(n-1)P_{2}(q^{n}),$ (357)

for $n\in(1,2)$. Then,

$\displaystyle\\{C_{n},G\\}$ $\displaystyle=\left\langle
q,n\nabla^{\perp}\left(\mathring{P}_{0}q^{n-1}\right)\cdot\frac{\delta
G}{\delta u}^{\perp}\right\rangle$
$\displaystyle\qquad-\left\langle\underbrace{\nabla\cdot
n\nabla^{\perp}\left(\mathring{P}_{0}q^{n-1}\right)}_{=0},\frac{\delta
G}{\delta D}\right\rangle+\left\langle\nabla\cdot\frac{\delta G}{\delta
u},-(n-1)P_{2}(q^{n})\right\rangle$ $\displaystyle\qquad+\left\langle\nabla
n\mathring{P}_{0}^{\perp}q^{n-1},q\frac{\delta G}{\delta
u}\right\rangle-\left\langle\nabla\frac{\delta G}{\delta
Z^{\prime}},n\nabla^{\perp}\left(\mathring{P}_{0}q^{n-1}\right)\right\rangle,$
(358) $\displaystyle=\left\langle
q,n\nabla\underbrace{\left(\mathring{P}_{0}q^{n-1}+\mathring{P}_{0}^{\perp}q^{n-1}\right)}_{=q^{n-1}}\cdot\frac{\delta
G}{\delta u}\right\rangle$
$\displaystyle\qquad-\left\langle\underbrace{\nabla\cdot
n\nabla^{\perp}q^{n-1}}_{=0},\frac{\delta G}{\delta
D}\right\rangle+\left\langle\nabla\cdot\frac{\delta G}{\delta
u},-(n-1)P_{2}(q^{n})\right\rangle$
$\displaystyle\qquad+\left\langle\underbrace{\nabla^{\perp}\cdot\nabla\frac{\delta
G}{\delta
Z^{\prime}}}_{=0},n\left(\mathring{P}_{0}q^{n-1}\right)\right\rangle,$ (359)
$\displaystyle=-\left\langle(n-1)\nabla q^{n},\frac{\delta G}{\delta
u}\right\rangle+\left\langle\nabla\cdot\frac{\delta G}{\delta
u},-(n-1)q^{n}\right\rangle$
$\displaystyle=\left\langle(n-1)q^{n},\nabla\cdot\frac{\delta G}{\delta
u}\right\rangle-\left\langle\nabla\cdot\frac{\delta G}{\delta
u},(n-1)q^{n}\right\rangle=0,$ (360)

where we have repeatedly used $q^{n}\in\mathbb{V}_{h}^{0}$ for $n\in(1,2)$.
The surface integral vanished in the penultimate line since
$\mathring{P}_{0}q^{n-1}\in\mathring{\mathbb{V}}_{h}^{0}$ vanishes on the
boundary, and the surface integral vanished in the final line since the normal
component of $\frac{\delta G}{\delta u}\in\mathring{\mathbb{V}}_{h}^{1}$
vanishes on the boundary. Hence the total vorticity $C_{1}$ and the enstrophy
$C_{2}$ are both Casimirs and hence are conserved by the dynamics.

Bauer and Cotter (2018) demonstrated through numerical experiments that this
scheme produces convergent solutions. This idea is closely related to the
approach of Ketefian and Jacobson (2009), who introduced vorticity variables
at the boundary to make an energy enstrophy conserving staggered finite
difference method.

###  7.7 Thermal shallow water equations

The thermal shallow water equations provide a useful stepping stone between
the rotating shallow water equations and three dimensional models. This is
because they incorporate an additional advected tracer, the temperature,
whilst remaining in the two dimensional setting. Additionally, they provide an
interesting reduced model for describing some atmospheric processes,
especially when further augmented with a moisture variable, as discussed in
the excellent book by Zeitlin (Zeitlin 2018). The variational derivation of
these equations originates from Ripa (1993), and can be placed in the
framework of Euler-Poincaré equations by treating the buoyancy $s$ as an
additional advected quantity satisfying

$s_{t}+u\cdot\nabla s=0.$ (361)

Then, the equations take the form

$\displaystyle\frac{\partial u}{\partial t}+u\cdot\nabla
u+fu^{\perp}+s\nabla(D+b)+\frac{D}{2}\nabla s$ $\displaystyle=0,$ (362)
$\displaystyle\frac{\partial D}{\partial t}+\nabla\cdot(Du)$
$\displaystyle=0,$ (363) $\displaystyle\frac{\partial s}{\partial
t}+u\cdot\nabla s$ $\displaystyle=0,$ (364)

where $u$ is the horizontal velocity and $D$ is the layer thickness as before.
Instead of $s$, one can instead work with the bouyancy density $S=sD$, which
satisfies

$S_{t}+\nabla\cdot(uS)=0.$ (365)

As a consequence of the variational derivation, these equations have Poisson
bracket formulations. When $S$ is used as the prognostic variable, the
equations can be obtained from the following Poisson bracket,

$\displaystyle\\{F,G\\}_{S}$
$\displaystyle=\\{F,G\\}_{0}+\left\langle\frac{\delta F}{\delta
S},\nabla\cdot\left(s\frac{\delta G}{\delta
u}\right)\right\rangle-\left\langle\frac{\delta G}{\delta
S},\nabla\cdot\left(s\frac{\delta F}{\delta u}\right)\right\rangle,$ (366)

where $\\{\cdot,\cdot\\}_{0}$ is the Poisson bracket presented in (284),
combined with the Hamiltonian

$H=\int_{\Omega}\frac{D|u|^{2}}{2}+S\left(\frac{D}{2}+b\right)\operatorname{d}x.$
(367)

This Poisson bracket formulation has Casimirs of the form

$C[D,u,S]=\int_{\Omega}DqA\left(\frac{S}{D}\right)+DB\left(\frac{S}{D}\right)\operatorname{d}x,$
(368)

where $A$ and $B$ are arbitrary functions.

Alternatively, when $s$ is used as the prognostic variable, we have the
following Poisson bracket formulation,

$\displaystyle\\{F,G\\}_{s}$
$\displaystyle=\\{F,G\\}_{0}-\left\langle\nabla\frac{\delta F}{\delta
s},s\frac{\delta G}{\delta u}\right\rangle+\left\langle\nabla\frac{\delta
G}{\delta s},s\frac{\delta F}{\delta u}\right\rangle,$ (369) $\displaystyle H$
$\displaystyle=\int_{\Omega}\frac{D|u|^{2}}{2}+Ds\left(\frac{D}{2}+b\right)\operatorname{d}x.$
(370)

Similarly, this Poisson bracket formulation has Casimirs of the form

$C[D,u,s]=\int_{\Omega}DqA\left(s\right)+DB\left(s\right)\operatorname{d}x.$
(371)

The advantage of the $S$ formulation is that conservation of total buoyancy

$B=\int_{\Omega}S\operatorname{d}x,$ (372)

is naturally incorporated, and local conservation is possible when choosing
$S\in\mathbb{V}_{h}^{2}$. In both cases, for a conforming discretisation,
$s\in H^{1}$ is required, i.e. we should take $s\in\mathbb{V}_{h}^{0}$. In the
case of the $S$ formulation, following the approach to $q$ in McRae and Cotter
(2014), Eldred et al. (2019) proposed to introduce $s\in\mathbb{V}_{h}^{0}$ as
a diagnostic quantity defined by

$\int_{\Omega}\gamma Ds\operatorname{d}x=\int_{\Omega}\gamma
S\operatorname{d}x,\quad\forall s\in\mathbb{V}_{h}^{0}.$ (373)

Alternatively, a nonconforming discretisation can be obtained by introducing
additional facet terms into the Poisson bracket; we shall discuss this further
in section 7.13.

In the case of the conforming $S$ formulation with $S\in\mathbb{V}_{h}^{2}$
and prognostic $s\in\mathbb{V}_{h}^{0}$, the variational derivatives of $H$
with respect to $u\in\mathbb{V}_{h}^{1}$, $D\in\mathbb{V}_{h}^{2}$ and
$S\in\mathbb{V}_{h}^{2}$ become

$\frac{\delta H}{\delta u}=m:=P_{1}(uD),\,\frac{\delta H}{\delta
D}=P_{2}(|u|^{2}/2)+S/2,\,\frac{\delta H}{\delta S}=\frac{D}{2}+b.$ (374)

Picking $F=\left\langle
u,w\right\rangle+\left\langle\phi,D\right\rangle+\left\langle
S,\alpha\right\rangle$ and using the Poisson bracket (366) then gives the
system of equations

$\displaystyle\left\langle w,u_{t}\right\rangle+\left\langle
w,qm^{\perp}\right\rangle$ $\displaystyle\quad-\left\langle\nabla\cdot
w,\frac{|u|^{2}}{2}+\frac{S}{2}\right\rangle-\left\langle\frac{D}{2}+b,\nabla\cdot\left(sw\right)\right\rangle$
$\displaystyle=0,\quad\forall w\in\mathbb{V}_{h}^{1},$ (375)
$\displaystyle\left\langle\phi,D_{t}\right\rangle+\left\langle\phi,\nabla\cdot
m\right\rangle=0,\quad\forall\phi\in\mathbb{V}_{h}^{2},$ (376)
$\displaystyle\left\langle\alpha,S_{t}\right\rangle+\left\langle\alpha,\nabla\cdot(sD)\right\rangle$
$\displaystyle=0,\quad\forall\alpha\in\mathbb{V}_{h}^{2},$ (377)
$\displaystyle\left\langle\gamma,qD\right\rangle-\left\langle\nabla^{\perp}\gamma,u\right\rangle-\left\langle\gamma,f\right\rangle$
$\displaystyle=0,\quad\forall\gamma\in\mathbb{V}_{h}^{0},$ (378)
$\displaystyle\left\langle\beta,sD\right\rangle-\left\langle\beta,S\right\rangle$
$\displaystyle=0,\quad\forall\beta\in\mathbb{V}_{h}^{0}.$ (379)

The introduction of the thermal variable breaks the symmetry that makes
$C_{n}$ a Casimir for $n>1$ for the shallow water equations. However,
straightforward calculations show that mass $M$, total buoyancy $B$ and total
vorticity $Z$ are all Casimirs for this discrete bracket.

A similar conforming discretisation obtained from the Poisson bracket (369)
also preserves all three of these quantities, with the variation that now $B$
is a quadratic functional

$B=\int_{\Omega}sD\operatorname{d}x,$ (380)

so requires a time integrator that preserves quadratic Casimirs.

Eldred et al. (2019) introduced these formulations plus a number of
nonconforming versions with facet integrals. They introduced time integration
methods that conserve the relevant Casimirs and demonstrated all of these
schemes in convergence tests and other benchmarks.

###  7.8 Rotating compressible Euler equations

Using the approaches described in this section, Poisson bracket
discretisations are possible for any variational fluid model with a
Hamiltonian being a function of velocity $u$, density $D$, and a thermal field
$\theta$ (as well as extensions to e.g. magnetic flux $B$ etc.). In this
section we briefly discuss such discretisations for the compressible Euler
equations that are the basis for atmospheric dynamical cores in weather and
climate models. These equations are given (in the “$\theta$-$\Pi$”
formulation) by

$\displaystyle u_{t}+(u\cdot\nabla)u+2\Omega\times u+c_{p}\theta\nabla\Pi$
$\displaystyle=-g\hat{z},$ (381) $\displaystyle\theta_{t}+u\cdot\nabla\theta$
$\displaystyle=0,$ (382) $\displaystyle D_{t}+\nabla\cdot(uD)$
$\displaystyle=0,$ (383) $\displaystyle\Pi^{\frac{1-\kappa}{\kappa}}$
$\displaystyle=\frac{R}{p_{0}}D\theta,$ (384)

where $u$ is the velocity, $\Omega$ is the rotation vector for the Earth,
$\theta$ is the potential temperature (a scaling of temperature that absorbs
the changes in temperature due to changes in pressure), $\Pi$ is the Exner
function, $g$ is the acceleration due to gravity, $\hat{z}$ is the unit vector
pointing away from the centre of the Earth, $D$ is the density,
$\kappa=R/c_{p}$, $R$ is the ideal gas constant, $c_{p}=R+c_{v}$ is the
specific heat at constant pressure, $c_{v}$ is the specific heat at constant
volume, and $p_{0}$ is a reference pressure used to define $\theta$.

One Poisson bracket formulation for these equations based around the three
dimensional vorticity vector $\omega=\nabla\times u+2\Omega$ is

$\displaystyle\\{F,G\\}$ $\displaystyle=\left\langle\frac{\delta F}{\delta
u},\omega\times\frac{\delta G}{\delta u}\right\rangle+\left\langle\frac{\delta
F}{\delta D},\nabla\cdot\frac{\delta G}{\delta
u}\right\rangle+\left\langle\frac{1}{D}\frac{\delta
F}{\delta\theta}\nabla\theta,\frac{\delta G}{\delta u}\right\rangle$
$\displaystyle\qquad-\left\langle\frac{\delta G}{\delta
D},\nabla\cdot\frac{\delta F}{\delta
u}\right\rangle-\left\langle\frac{1}{D}\frac{\delta
G}{\delta\theta}\nabla\theta,\frac{\delta F}{\delta u}\right\rangle,$ (385)

with Hamiltonian

$H=\int_{\Omega}\frac{D|u|^{2}}{2}+Dgz+\underbrace{c_{v}D\theta\Pi}_{\mbox{thermal
energy}}\operatorname{d}x,$ (386)

where $z$ is the height above some reference altitude. Other Poisson bracket
formulations are also possible, notably with $\Theta=D\theta$ instead of
$\theta$, but we do not intend to be encyclopaedic here.

Similarly to the treatment of the thermal shallow water equations, a
conforming discretisation requires that $\theta\in H^{1}$, i.e. we should take
$\theta\in\mathbb{W}_{h}^{0}$. We take $u\in\mathbb{W}_{h}^{2}$,
$D\in\mathbb{W}_{h}^{3}$. The discrete variational derivatives of $H$ are then

$\displaystyle\frac{\delta H}{\delta u}$ $\displaystyle=m:=P_{2}(Du),$ (387)
$\displaystyle\frac{\delta H}{\delta D}$
$\displaystyle=P_{3}\left(\frac{1}{2}|u|^{2}+gz+c_{p}\theta\Pi\right),$ (388)
$\displaystyle\frac{\delta H}{\delta\theta}$
$\displaystyle=c_{p}P_{0}\left(D\Pi\right),$ (389)

where the derivatives of the internal energy, the third term in $H$, require a
little algebra. Using the Poisson bracket (7.8) as the basis for a
discretisation requires an approximation of $\nabla\times u$ since
$\mathbb{W}_{h}^{2}$ is not a curl-conforming space. Analogously to (304), we
approximate $\omega\in\mathbb{W}_{h}^{1}$ such that

$\left\langle v,\omega\right\rangle-\left\langle\nabla\times
v,u\right\rangle=0,\quad\forall v\in\mathbb{W}_{h}^{1},$ (390)

for the case of a domain without boundaries. When boundaries are present, as
indeed they must be since the gravitational potential energy term $Dgz$ does
not work if the domain is periodic in the vertical, then we must again define
$\mathring{\mathbb{W}}_{h}^{1}$ as

$\mathring{\mathbb{W}}_{h}^{1}=\left\\{\omega\in\mathbb{W}_{h}^{1}:\omega\times
n=0\mbox{ on }\partial\Omega\right\\},$ (391)

where $n$ is the outward pointing normal to $\Omega$. Then the state space
must be extended to include
$Z^{\prime}\in(\mathring{\mathbb{W}}_{h}^{1})^{\perp}$, the $L^{2}$-orthogonal
complement to $\mathring{\mathbb{W}}_{h}^{1}$ in $\mathbb{W}_{h}^{1}$, which
represents vorticity components on the boundary which have their own dynamics
consistent with the conservation of total vorticity. We ignore this aspect for
now, but return to it in Section 8.6.

This construction leads to the discretisation

$\displaystyle\left\langle w,u_{t}\right\rangle+\left\langle w,\omega\times
m\right\rangle$ $\displaystyle\qquad-\left\langle\nabla\cdot
w,\frac{1}{2}|u|^{2}+gz+c_{p}\theta\Pi\right\rangle-\left\langle
w,\frac{1}{D}s\nabla\theta\right\rangle$ $\displaystyle=0,\quad\forall
w\in\mathbb{W}_{h}^{2},$ (392)
$\displaystyle\left\langle\phi,D_{t}+\nabla\cdot m\right\rangle$
$\displaystyle=0,\quad\forall\phi\in\mathbb{W}_{h}^{3},$ (393)
$\displaystyle\left\langle\gamma,\theta_{t}\right\rangle+\left\langle\gamma\frac{1}{D}\alpha\nabla\theta,m\right\rangle$
$\displaystyle=0,\quad\forall\gamma\in\mathbb{W}_{h}^{0},$ (394)
$\displaystyle\left\langle r,m-uD\right\rangle$ $\displaystyle=0,\quad\forall
r\in\mathbb{W}_{h}^{2},$ (395) $\displaystyle\left\langle
v,\omega\right\rangle-\left\langle\nabla\times v,u\right\rangle$
$\displaystyle=0,\quad\forall v\in\mathbb{W}_{h}^{1},$ (396)
$\displaystyle\left\langle\alpha,s-c_{p}D\Pi\right\rangle$
$\displaystyle=0,\quad\forall\alpha\in\mathbb{W}_{h}^{0},$ (397)
$\displaystyle\Pi^{\frac{1-\kappa}{\kappa}}$
$\displaystyle=\frac{R}{p_{0}}D\theta.$ (398)

In this formulation, there is no approximation in the definition of $\Pi$. The
equations have nonpolynomial terms due to the fractional powers in the
definition of $\Pi$, which cannot be integrated exactly. This can be dealt
with by replacing the exact integral in the definition of the variational
derivatives with a quadrature rule. Care must be taken to use this quadrature
rule consistently in all of the terms to obtain an energy conserving
formulation. In contrast to $\frac{\delta H}{\delta D}$, it is not possible to
remove the projection $P_{\theta}$ in the definition of $\frac{\delta
H}{\delta\theta}$ from the equations, because $\frac{\delta H}{\delta\theta}$
does not appear in an inner product with a function from $\mathbb{V}_{h}^{0}$.
Hence, we have to introduce a third auxiliary variable $s$. Lee and Palha
(2021) used a related formulation to build a discretisation using mimetic
spectral elements.

Focussing on the approximation of the pressure gradient term
$c_{p}\theta\nabla\Pi$, the relevant terms are

$\displaystyle-\left\langle\nabla\cdot
w,c_{p}\theta\Pi\right\rangle-\left\langle
w,\frac{1}{D}s\nabla\theta\right\rangle$
$\displaystyle=-\left\langle\nabla\cdot(\theta w),c_{p}\Pi\right\rangle$
$\displaystyle\qquad+\left\langle
w,\left(c_{p}\Pi-\frac{1}{D}s\right)\nabla\theta\right\rangle,$ (399)

which is a consistent approximation to $c_{p}\theta\nabla\Pi$ since $s/D$ only
differs from $c_{p}\Pi$ by multiplication by $D$, projection to
$\mathbb{W}_{h}^{0}$ and division by $D$ again.

As discussed earlier in this article, it can be preferable to use the
temperature space $\mathbb{W}_{\theta}$ for $\theta$, which is more compatible
with hydrostatic balance. Since $\mathbb{W}_{\theta}$ allows discontinuities
in the horizontal direction, we need to modify the Poisson bracket formulation
to incorporate the nonconforming discretisation. This is done by focussing on
the following term in the Poisson bracket,

$\left\langle\frac{1}{D}\frac{\delta F}{\delta\theta}\nabla\theta,\frac{\delta
G}{\delta u}\right\rangle=\left\langle\frac{\delta
F}{\delta\theta},\frac{1}{D}\frac{\delta G}{\delta
u}\cdot\nabla\theta\right\rangle,$ (400)

which leads to the term in (394) approximating $u\cdot\nabla\theta$. To adapt
this to the partially discontinuous space $\mathbb{W}_{\theta}$, we replace
with the discontinuous Galerkin discretisation,

$L\left[\frac{\delta G}{\delta u},\frac{\delta
F}{\delta\theta};\theta\right]=-\left\langle\nabla_{h}\cdot\left(\frac{1}{D}\frac{\delta
G}{\delta u}\frac{\delta
F}{\delta\theta}\right),\theta\right\rangle+\left\llangle\left[\\!\\!\left[\frac{1}{D}\frac{\delta
G}{\delta u}\frac{\delta
F}{\delta\theta}\right]\\!\\!\right],\\{\theta\\}\right\rrangle_{\Gamma},$
(401)

using the discontinuous Galerkin notation as introduced in Section 4. Here we
have chosen a centred flux $\\{\theta\\}$ but will discuss upwind fluxes in
7.13. As usual, this is a consistent approximation with the facet integrals
vanishing when $u$, $D$, $\theta$, etc are all smooth functions. Then, we use
$L$ in a modified Poisson bracket,

$\displaystyle\\{F,G\\}$ $\displaystyle=\left\langle\frac{\delta F}{\delta
u},\omega\times\frac{\delta G}{\delta u}\right\rangle+\left\langle\frac{\delta
F}{\delta D},\nabla\cdot\frac{\delta G}{\delta
u}\right\rangle+L\left[\frac{\delta G}{\delta u},\frac{\delta
F}{\delta\theta};\theta\right]$ $\displaystyle\qquad-\left\langle\frac{\delta
G}{\delta D},\nabla\cdot\frac{\delta F}{\delta
u}\right\rangle-L\left[\frac{\delta F}{\delta u},\frac{\delta
G}{\delta\theta};\theta\right].$ (402)

This leads to the modified velocity equation,

$\displaystyle\left\langle w,u_{t}\right\rangle+\left\langle w,\omega\times
m\right\rangle-\left\langle\nabla\cdot
w,\frac{1}{2}|u|^{2}+gz+c_{p}\theta\Pi\right\rangle$
$\displaystyle\quad+\left\langle\nabla_{h}\cdot\left(\frac{1}{D}ws\right),\theta\right\rangle-\left\llangle\left[\\!\\!\left[\frac{1}{D}ws\right]\\!\\!\right],\\{\theta\\}\right\rrangle_{\Gamma}$
$\displaystyle=0,\quad\forall w\in\mathbb{W}_{h}^{2},$ (403)

in which we recognise another consistent approximation of the pressure
gradient term. This concept for introducing any chosen $\theta$ advection
scheme into a Poisson bracket formulation was introduced in Gassmann and
Herzog (2008), with application to a global finite difference model in
Gassmann (2013).

Here we should note that formulating a Poisson time integrator for this system
is challenging, because the integrands are nonpolynomial. As noted to us by
Chris Eldred, and implemented in Wimmer et al. (2021), a Poisson integrator
must be approximated by using an incomplete quadrature rule in the time
averaged variational derivatives. This can be done to high order at the
expense of a more complicated assembly.

###  7.9 Upwinding for incompressible Euler: SUPG

One of the interesting and useful features of the Poisson bracket formulation
is that it can be modified to incorporate stabilisation of the transport
schemes whilst remaining antisymmetric and hence energy conserving. In the
context of two dimensional incompressible turbulence (and geostrophic
turbulence), this is useful because energy cascades to large scales (where
functions can be well approximated by finite element functions), whilst
enstrophy cascades to small scales (where they can not). In energy enstrophy
conserving schemes, vorticity features pile up at the grid scale, leading to
unphysical noise, when really they should be cascading to scales below the
grid scale. Here a scheme that conserves energy whilst dissipating enstrophy
at the small scale through upwind stabilisation is appropriate. In three
dimensional isotropic turbulence, energy is also cascading towards small
scales, and so additional dissipative gridscale closures or parameterisations
are necessary. If these are added to a Poisson bracket formulation with upwind
stabilisation, we know that there are no spurious energy transfers between
scales and between potential, kinetic and internal energy, and the only energy
changes are due to the additional dissipative closures and parameterisations.
If desired, the energy dissipated from those terms can be collected and
recycled into subgrid closures, as is done in Gassmann (2013).

First we discuss energy conserving upwinding techniques for the advection term
in the velocity equation. For incompressible quasigeostrophic models, Sadourny
and Basdevant (1985) proposed a subgrid closure within the Arakawa Jacobian
finite difference formulation by replacing $q\to q-\tau u\cdot\nabla u$ in the
Poisson bracket, where $\tau$ is a chosen timescale. This provides upwinding
by approximating the value of $q$ taken upstream along the streamline passing
through the gridpoint, hence the name Anticipated Potential Vorticity Method
(APVM). This was included in a rotating shallow water formulation by Arakawa
and Hsu (1990), and has been included in more recent unstructured grid
formulations in Ringler et al. (2010), Chen et al. (2012).

When regarded as a numerical scheme instead of a turbulence closure, this
modification appears as a $\mathcal{O}(\tau)$ consistency error. Instead, in
the context of incompressible Euler, one can replace
$\omega\to\omega-\tau(\omega_{t}+u\cdot\nabla\omega)$. Here, the idea is that
this term vanishes when the approximation of the solution is accurate and
smooth, since then the vorticity equation $\omega_{t}+u\cdot\nabla\omega=0$ is
well approximated. Hence, the approximation is consistent. Following this
idea, we can modify the incompressible Euler bracket (269) to become

$\\{F,G\\}=\int_{\Omega}\left(\omega-\tau\left(\omega_{t}+u\cdot\nabla\omega\right)\right)\frac{\delta
F}{\delta u}\cdot\frac{\delta G}{\delta u}^{\perp}\operatorname{d}x,$ (404)

where $\omega$ is obtained from (268) as usual. The discretisation becomes

$\left\langle
w,u_{t}\right\rangle+\left\langle\left(\omega-\tau\left(\omega_{t}+u\cdot\nabla\omega\right)\right)w,u^{\perp}\right\rangle=0,\quad\forall
w\in\zeta_{h}.$ (405)

Writing $u=\nabla^{\perp}\psi$, $w=\nabla^{\perp}\phi$ for
$\psi,\phi\in\mathbb{V}_{h}^{0}$, we get

$\underbrace{\left\langle\nabla\phi,\nabla\psi_{t}\right\rangle}_{=-\left\langle\phi,\omega_{t}\right\rangle}+\left\langle\left(\omega-\tau\left(\omega_{t}+u\cdot\nabla\omega\right)\right)\nabla\phi,u\right\rangle=0,\quad\forall\phi\in\mathbb{V}_{h}^{0},$
(406)

which we rewrite as

$\left\langle\phi+\tau
u\cdot\nabla\phi,\omega_{t}\right\rangle+\left\langle\phi+\tau
u\cdot\nabla\phi,u\cdot\nabla\omega\right\rangle=0,\quad\forall\phi\in\mathbb{V}_{h}^{0},$
(407)

after integrating by parts in the term
$-\left\langle\omega\nabla\phi,u\right\rangle$, which is permissible since
$\omega,\phi\in H^{1}(\Omega)$ and $u\in H(\operatorname{div})$. This is the
SUPG discretisation of the incompressible Euler equation, which is obtained by
replacing the test function $\phi$ with $\phi+\tau u\cdot\nabla\phi$. The
additional term leads to streamwise diffusion of the vorticity $\omega$
without harming the consistency of the scheme.

Since the scheme as written here is derived from a Poisson bracket
formulation, it conserves energy by construction. Regarding Casimirs $C_{n}$,
we now get

$\displaystyle\left\\{F,C_{n}\right\\}$
$\displaystyle=-n\int_{\Omega}(\omega+\tau u\cdot\nabla\omega)\frac{\delta
F}{\delta
u}\cdot\left(\nabla^{\perp}\omega^{n-1}\right)^{\perp}\operatorname{d}x,$
(408) $\displaystyle=n\int_{\Omega}(\omega+\tau
u\cdot\nabla\omega)\frac{\delta F}{\delta
u}\cdot\nabla\omega^{n-1}\operatorname{d}x,$ (409)
$\displaystyle=(n-1)\int_{\Omega}\frac{\delta F}{\delta
u}\cdot\nabla\omega^{n}\operatorname{d}x-n\tau\int_{\Omega}u\cdot\nabla\omega\frac{\delta
F}{\delta
u}\cdot\left(\nabla^{\perp}\omega^{n-1}\right)^{\perp}\operatorname{d}x,$
(410)

and this latter term only vanishes when $n=1$, hence we have conservation of
total vorticity but not enstrophy. By substituting $\phi=\omega$ into (407),
we can obtain the enstrophy dynamics

$\displaystyle\frac{\operatorname{d}}{\operatorname{d}t}\int_{\Omega}\frac{1}{2}\omega^{2}\operatorname{d}x$
$\displaystyle=\left\langle\omega,\omega_{t}\right\rangle,$ (411)
$\displaystyle=\left\langle\left(\omega-\tau\left(\omega_{t}+u\cdot\nabla\omega\right)\right)\nabla\omega,u\right\rangle=0,$
(412) $\displaystyle=\underbrace{\left\langle
u,\frac{1}{2}\nabla\omega^{2}\right\rangle}_{=0}-\tau\left\langle\omega_{t},u\cdot\nabla\omega\right\rangle-\tau\left\langle
u\cdot\nabla\omega,u\cdot\nabla\omega\right\rangle,$ (413)

where the first term in the last line vanishes after integration by parts and
noting that $u$ is divergence free. The last term is negative semidefinite,
and corresponds to diffusion of enstrophy along streamlines, as occurs in
APVM. The middle term, which is what we get if we change from APVM to SUPG to
ensure consistency of the scheme, is indefinite. However, it only contains one
derivative, so it is a lower order term compared to the streamwise diffusion,
and hence the streamwise component of $\omega$ is kept smooth.

###  7.10 Upwinding for incompressible Euler: vorticity free formulation

An alternative Poisson bracket with upwinding for the two dimensional
incompressible Euler equations stems from the variational formulation of
(Natale and Cotter 2018), discussed in the previous section. Since that
discretisation conserves energy, it should not come as a surprise that it has
a Poisson bracket formulation, given by

$\displaystyle\\{F,G\\}$ $\displaystyle=-\left\langle
u,\nabla^{\perp}\left(\frac{\delta F}{\delta u}^{\perp}\cdot\frac{\delta
G}{\delta
u}\right)\right\rangle+\left\llangle\\{u\\},\left[\\!\\!\left[n^{\perp}\cdot\left(\frac{\delta
F}{\delta u}^{\perp}\cdot\frac{\delta G}{\delta
u}\right)\right]\\!\\!\right]\right\rrangle_{\Gamma},$ (414) $\displaystyle H$
$\displaystyle=\frac{1}{2}\int_{\Omega}|u|^{2}\operatorname{d}x.$ (415)

This formulation can be thought of as an alternative way to obtain an
approximation to $\omega$ when $u\in\mathbb{V}_{h}^{1}$. Instead of using an
auxiliary equation to define $\omega$, here we integrate the curl by parts in
each cell and choose an approximation to $u$ on the facets. This approximation
is necessary because although $u$ has continuous normal components, it does
not have continuous tangential components; the tangential component of $u$ is
multivalued on the boundary. The variational derivation leads to a centred
approximation $\\{u\\}$, but we can equally take an upwind approximation
$\tilde{u}$ (where $\tilde{u}$ is the value of $u$ on the upwind side of the
facet). This leads to

$\left\langle w,u_{t}\right\rangle-\left\langle
u,\nabla^{\perp}\left(w^{\perp}\cdot
u\right)\right\rangle+\left\llangle\\{u\\},\left[\\!\\!\left[n^{\perp}\cdot\left(w^{\perp}\cdot
u\right)\right]\\!\\!\right]\right\rrangle_{\Gamma}=0,\quad\forall
w\in\zeta^{\perp}.$ (416)

It is more difficult to diagnose the enstrophy budget for this scheme than it
was for the SUPG scheme. However, numerical experiments in Natale and Cotter
(2018) showed that this scheme does indeed tend to reduce enstrophy whilst
exactly conserving energy. They also proved convergence of the upwinded
scheme, albeit at a suboptimal rate; numerical experiments showed convergence
at the rate expected given the degree of the polynomials (_i.e._ , second
order $L^{2}$ convergence for $BDM_{1}$).

###  7.11 Scale selective dissipation

Natale and Cotter (2017) investigated the ability of the schemes in Sections
7.9 and 7.10 to produce energy backscatter from small to large scales
consistently with two dimensional turbulence in the forced dissipative
setting.

For the scheme of 7.10 they provided a multiscale interpretation of the
discretisation, in the case of the $BDM_{1}$ space (the approach is general
for $BDM$ spaces, but discussion is simplified by just describing the lowest
order case). Since $\mathbb{V}_{h}^{1}=$$BDM_{1}$ contains
$\mathbb{V}_{h}^{1,l}=$P12 (the space of vector valued continuous piecewise
linear functions) as a subspace, we can write an $L^{2}$ orthogonal
decomposition,

$\mathbb{V}_{h}^{1}=\mathbb{V}_{h}^{1,l}\oplus\mathbb{V}_{h}^{1,s},$ (417)

where $\mathbb{V}_{h}^{1,s}$ is the $L^{2}$ orthogonal complement of
$\mathbb{V}_{h}^{1,l}$ in $\mathbb{V}_{h}^{1}$. Since $\mathbb{V}_{h}^{1,l}$
is a continuous finite element space and therefore contains functions that are
smoother than $\mathbb{V}_{h}^{1,s}$, we can consider $\mathbb{V}_{h}^{1,l}$
to be a subspace of larger scale fields whilst $\mathbb{V}_{h}^{1,s}$ contains
the small scales. This decomposition is not compatible with the decomposition
$\mathbb{V}_{h}^{1}=\zeta\oplus\zeta^{\perp}$, so we have to use the mixed
formulation where the divergence free condition is enforced explicitly _via_
the pressure gradient term. The formulation may be written as

$\displaystyle\left\langle v,u_{t}\right\rangle+a(u;u,v)+s(u;u,v)-\left\langle
P,\nabla\cdot v\right\rangle$ $\displaystyle=0,\quad\forall
v\in\mathbb{V}_{h}^{1},$ (418) $\displaystyle\left\langle\nabla\cdot
u,\phi\right\rangle$ $\displaystyle=0,\quad\forall\phi\in\mathbb{V}_{h}^{0},$
(419)

where

$\displaystyle a(\hat{u};u,v)$
$\displaystyle=\left\langle\hat{u}^{\perp},\nabla(u^{\perp}\cdot
v)\right\rangle-\left\llangle\\{\hat{u}^{\perp}\\}\cdot
n_{+}\left[\\!\\!\left[u^{\perp}\cdot{v}\right]\\!\\!\right]\right\rrangle_{\Gamma},$
(420) $\displaystyle s(\hat{u};u,v)$ $\displaystyle=-\left\llangle
c_{+}\left[\\!\\!\left[\hat{u}^{\perp}\right]\\!\\!\right]\cdot
n_{+}\left[\\!\\!\left[u^{\perp}\cdot{v}\right]\\!\\!\right]\right\rrangle_{\Gamma},$
(421)

noting that we use $\hat{u}=u$ in (418), and where $c_{+}$ being equal to 1 if
$u\cdot n_{+}\geq 0$ and 0 otherwise.

If we now write $u=u^{l}+u^{s}$, with $u^{l}\in\mathbb{V}_{h}^{1,l}$ and
$u^{s}\in\mathbb{V}_{h}^{1,s}$, we observe that

$\displaystyle s(u;u,u^{l})$ $\displaystyle=-\left\llangle
c_{+}\left[\\!\\!\left[u^{\perp}\right]\\!\\!\right]\cdot
n_{+},\left[\\!\\!\left[u^{\perp}\cdot{u^{l}}\right]\\!\\!\right]\right\rrangle_{\Gamma},$
(422) $\displaystyle=-\left\llangle
c_{+}\left[\\!\\!\left[u^{\perp}\right]\\!\\!\right]\cdot
n_{+},\left[\\!\\!\left[u^{\perp}\right]\\!\\!\right]\cdot
u^{l}\right\rrangle_{\Gamma},$ (423) $\displaystyle=-\left\llangle
c_{+}\left[\\!\\!\left[u^{\perp}\cdot
n_{+}\right]\\!\\!\right],\left[\\!\\!\left[u^{\perp}\cdot
n_{+}\right]\\!\\!\right]n_{+}\cdot u^{l}\right\rrangle_{\Gamma},$ (424)
$\displaystyle=-\left\llangle c_{+}\left[\\!\\!\left[u^{\perp}\cdot
n_{+}\right]\\!\\!\right]n_{+},n_{+}\left[\\!\\!\left[u^{\perp}\cdot
n_{+}\right]\\!\\!\right]n_{+}\cdot u^{l}\right\rrangle_{\Gamma},$ (425)
$\displaystyle=-\left\llangle
c_{+}\left[\\!\\!\left[u^{\perp}\right]\\!\\!\right],\left[\\!\\!\left[u^{\perp}\right]\\!\\!\right]n_{+}\cdot
u^{l}\right\rrangle_{\Gamma},$ (426) $\displaystyle=-\left\llangle
c_{+}\left[\\!\\!\left[u\right]\\!\\!\right],\left[\\!\\!\left[u\right]\\!\\!\right]n_{+}\cdot
u^{l}\right\rrangle_{\Gamma},$ (427)

where we used
$\left[\\!\\!\left[u^{\perp}\right]\\!\\!\right]=\left[\\!\\!\left[u^{\perp}\cdot
n_{+}\right]\\!\\!\right]n_{+}$ since $u^{\perp}$ has continuous tangential
components. Consequently,

$\displaystyle s(u;u,u^{s})$ $\displaystyle=-\left\llangle
c_{+}\left[\\!\\!\left[u^{\perp}\right]\\!\\!\right]\cdot
n_{+},\left[\\!\\!\left[u^{\perp}\cdot{u^{s}}\right]\\!\\!\right]\right\rrangle_{\Gamma},$
(428) $\displaystyle=-\left\llangle
c_{+}\left[\\!\\!\left[u^{\perp}\right]\\!\\!\right]\cdot
n_{+},\left[\\!\\!\left[u^{\perp}\right]\\!\\!\right]\cdot(u-u^{l})\right\rrangle_{\Gamma},$
(429) $\displaystyle=\left\llangle
c_{+}\left[\\!\\!\left[u^{\perp}\right]\\!\\!\right]\cdot
n_{+},\left[\\!\\!\left[u^{\perp}\right]\\!\\!\right]\cdot
u^{l}\right\rrangle_{\Gamma}-\left\llangle
c_{+}\left[\\!\\!\left[u^{\perp}\right]\\!\\!\right]\cdot
n_{+},\underbrace{\left[\\!\\!\left[u^{\perp}\cdot
u\right]\\!\\!\right]}_{=0}\right\rrangle_{\Gamma},$ (430)
$\displaystyle=\left\llangle
c_{+}\left[\\!\\!\left[u\right]\\!\\!\right],\left[\\!\\!\left[u\right]\\!\\!\right]n_{+}\cdot
u^{l}\right\rrangle_{\Gamma}.$ (431)

Since $\left\langle u^{l},u^{s}\right\rangle=0$, we can obtain equations for
the evolution of $E^{l}=\|u^{l}\|^{2}/2$ and $E^{s}=\|u^{s}\|/2$ by setting
$v=u^{l}$ and $v=u^{s}$, respectively,

$\displaystyle\frac{\operatorname{d}E^{l}}{\operatorname{d}t}+a(u;u,u^{l})-\left\langle
p,\nabla\cdot u^{l}\right\rangle$ $\displaystyle=\left\llangle c_{+}u^{l}\cdot
n_{+}\left[\\!\\!\left[u\right]\\!\\!\right],\left[\\!\\!\left[u\right]\\!\\!\right]\right\rrangle_{\Gamma},$
(432)
$\displaystyle\frac{\operatorname{d}E^{s}}{\operatorname{d}t}+a(u;u,u^{s})-\left\langle
p,\nabla\cdot u^{s}\right\rangle$ $\displaystyle=-\left\llangle
c_{+}u^{l}\cdot
n_{+}\left[\\!\\!\left[u\right]\\!\\!\right],\left[\\!\\!\left[u\right]\\!\\!\right]\right\rrangle_{\Gamma}.$
(433)

Since $c_{+}u^{l}\cdot n_{+}\geq 0$, the upwinding creates an energy transfer
from small to large scales.

Natale and Cotter (2017) demonstrated that this leads to energy backscatter in
practice, by using the analysis technique of Thuburn et al. (2014b). This
technique involves simulating two dimensional incompressible turbulence with
Newtonian damping and wavenumber 16 forcing. The instantaneous rate of change
of local energy is computed, and then Fourier transformed to obtain
$\frac{\partial}{\partial t}{E}(k)$, the rate of change of energy at
wavenumber $E$. Then the same solution is filtered by removing all wavenumbers
above a cutoff $k_{T}$, and the rate of change of energy is recomputed using
this filtered solution, to obtain $\frac{\partial}{\partial t}{E}_{T}(k)$.
Then, the rate of change of subgrid energy $E_{SG}(k)=\frac{\partial}{\partial
t}{E}(k)-\frac{\partial}{\partial t}{E}_{T}(k)$ is computed. This shows the
rate of change of energy at wavenumber $k$ due to wavenumbers $>k_{T}$. This
can then be compared with a high resolution reference solution which has a
much larger range of scales to support backscatter from. This computation is
then repeated using enstrophy instead of energy, obtaining $Z_{SG}(k)$. Natale
and Cotter (2017) examined the upwind scheme of Natale and Cotter (2018),
together with the SUPG scheme (405), using this technique, as well as
comparing the upwind momentum flux formulation of Guzmán et al. (2017), which
does not conserve energy. The experiments showed that all three schemes
exhibit the trough in $Z_{SG}(k)$ near $k=k_{T}$, which demonstrates that
enstrophy is being transported to smaller scales, consistent with the
enstrophy cascade. However, the upwind momentum flux formulation showed no
peak at low $k$ in $E_{SG}(k)$ that is indicative of the energy inverse
cascade in the reference solution. Both the Natale and Cotter (2018) and SUPG
schemes showed such a peak, although it is stronger and closer to the
reference solution for the SUPG scheme. When the upwinding is replaced by
centred approximation in the Natale and Cotter (2018) scheme, and when the
$\tau$ parameter is set to zero in the SUPG scheme, no statistical equilibrium
is reached and the implicit solvers for the systems eventually fail. Hence, we
conclude that energy conservation and some form of stabilisation by upwinding
or SUPG are both critical to obtaining these important features in two
dimensional forced dissipative turbulence.

###  7.12 Upwinding for rotating shallow water equations using potential
vorticity

Following Arakawa and Hsu (1990), McRae and Cotter (2014) demonstrated that an
energy conserving enstrophy dissipating scheme for the rotating shallow water
equations is possible using the APVM technique where $q$ is replaced by
$q-\tau u\cdot\nabla q$ in (284), extending the APVM idea discussed above to
the rotating shallow water equations. This was demonstrated to have a
beneficial effect on the smoothness of the solution whilst still preserving
energy. Natale and Cotter (2018) proposed to replace $q$ in (284) by
$q-\tau((Dq)_{t}+\nabla\cdot(mq))/D$ (or equivalently, by
$q-\tau(q_{t}+m\cdot\nabla q)$, since $D_{t}+\nabla\cdot m=0$ in $L^{2}$) in
order to obtain streamwise stabilisation within a consistent scheme, since
smooth solutions of the rotating shallow water equations satisfy
$(Dq)_{t}+\nabla\cdot(Duq)=0$. This leads to the system

$\displaystyle\left\langle
w,u_{t}\right\rangle+\left\langle\left(q-{\tau}(q_{t}+m\cdot\nabla
q\right)w,m^{\perp}\right\rangle$ $\displaystyle\qquad-\left\langle\nabla\cdot
w,\frac{1}{2}|u|^{2}+g(D+b)\right\rangle$ $\displaystyle=0,\,\quad\forall
w\in\mathbb{V}_{h}^{1},$ (434)
$\displaystyle\left\langle\phi,D_{t}\right\rangle+\left\langle\nabla\cdot
m\right\rangle$ $\displaystyle=0,\,\forall\phi\in\mathbb{V}_{h}^{2},$ (435)
$\displaystyle\left\langle\gamma,qD\right\rangle+\left\langle\nabla^{\perp}\gamma,u\right\rangle-\left\langle\gamma,f\right\rangle$
$\displaystyle=0,\,\forall\gamma\in\mathbb{V}_{h}^{0},$ (436)
$\displaystyle\left\langle v,m-uD\right\rangle$ $\displaystyle=0,\,\forall
v\in\mathbb{V}_{h}^{1}.$ (437)

To obtain the enstrophy dynamics, we use $w=\nabla^{\perp}q$ in (7.12) and
substitute with $\gamma=q$ in (436) to obtain

$\displaystyle\frac{\operatorname{d}}{\operatorname{d}t}\int_{\Omega}\frac{1}{2}q^{2}D\operatorname{d}x$
$\displaystyle=\left\langle\left(q-\tau(q_{t}+m\cdot\nabla
q)\right),m\cdot\nabla q\right\rangle,$ (438)
$\displaystyle=-\left\langle{\tau}(qD)_{t}m\cdot\nabla
q\right\rangle-\left\langle\nabla\cdot(qm)),m\cdot\nabla q\right\rangle,$
(439)

where we again have an indefinite consistency term and a streamwise diffusion
term that is always $\leq 0$. The energy conservation and stabilisation of
enstrophy, with decay of enstrophy when gridscale features arise through
vortex stretching, was demonstrated for this upwind stabilised rotating
shallow water scheme in Bauer and Cotter (2018) using numerical experiments.

The enrichment of this scheme with upwind stabilisation presents an
opportunity for a stability analysis which seems attainable at the time of
writing but is presently open. In particular, it would be interesting to
consider the stability of a backward Euler step for these equations.

###  7.13 Upwind discontinuous Galerkin methods for active tracers

Inspired by the observation of Gassmann and Herzog (2008) that upwinding for
advected quantities (such as layer depth, density, temperature etc.) can be
incorporated into an energy conserving scheme by simply ensuring that the
antisymmetry is maintained in the bracket, Wimmer et al. (2020) and Wimmer et
al. (2021) examined energy conserving tracer upwinding using upwind
discontinuous Galerkin schemes and SUPG schemes, respectively. Wimmer et al.
(2020) considered upwind discontinuous Galerkin schemes for the rotating
shallow water equations, applied to both the layer depth (which is in
$\mathbb{V}_{h}^{2}$ which allows arbitrary discontinuities between cells) and
the velocity (which is in $\mathbb{V}_{h}^{1}$, and hence allows discontinuity
in the tangential components, which is sufficient to allow for some
dissipation of small scale enstrophy).

The goal is to obtain an energy conserving formulation for which the layer
depth equation takes the upwind discontinuous Galerkin form,

$\left\langle\phi,D_{t}\right\rangle-\left\langle\nabla\phi,Du\right\rangle+\left\llangle\left[\\!\\!\left[\phi
u\right]\\!\\!\right],\tilde{D}\right\rrangle_{\Gamma}=0,\quad\forall\phi\in\mathbb{V}_{h}^{2}.$
(440)

To do this, we have to realise it as a modification of the component of the
Poisson bracket given by

$\left\langle\frac{\delta F}{\delta D},\nabla\cdot\frac{\delta H}{\delta
u}\right\rangle,$ (441)

where $\frac{\delta F}{\delta u}=\phi$ and $\frac{\delta H}{\delta
u}=P_{1}(Du)$. The problem is that $u$ explicitly appears in (440), whilst
$\frac{\delta H}{\delta u}$ involves a projection of $uD$. Wimmer et al.
(2020) solved this problem by introducing a recovery operator
$\mathcal{U}:(D,m)\in\mathbb{V}_{h}^{2}\times\mathbb{V}_{h}^{1}\to
u\in\mathbb{V}_{h}^{1}$ defined by

$\left\langle Dv,u\right\rangle=\left\langle v,m\right\rangle,\quad\forall
v\in\mathbb{V}_{h}^{1},$ (442)

which is well defined provided that $D>0$ (breaking this condition will cause
the scheme to fail anyway). Note in particular that if $m=P_{1}(Dw)$ for any
$w\in\mathbb{V}_{h}^{1}$, then $\mathcal{U}(D,m)=w$.

Thus, we can rewrite (440) as

$\displaystyle\left\langle\frac{\delta F}{\delta
D},D_{t}\right\rangle\underbrace{-\left\langle\nabla\frac{\delta F}{\delta
D},D\mathcal{U}\left(D,\frac{\delta H}{\delta
u}\right)\right\rangle+\left\llangle\left[\\!\\!\left[\frac{\delta F}{\delta
D}\mathcal{U}\left(D,\frac{\delta H}{\delta
u}\right)\right]\\!\\!\right],\tilde{D}\right\rrangle_{\Gamma}}=0,$
$\displaystyle\quad\forall\frac{\delta F}{\delta D}\in\mathbb{V}_{h}^{2},$
(443)

with the term marked with the underbrace replacing (441). In order to keep the
bracket antisymmetric, we also make the same substitutions in the
corresponding term with $F$ and $G$ exchanged, and the Poisson bracket becomes

$\displaystyle\\{F,G\\}=$ $\displaystyle\left\langle\frac{\delta F}{\delta
u},q\frac{\delta G}{\delta
u}^{\perp}\right\rangle{-\left\langle\nabla\frac{\delta F}{\delta
D},D\mathcal{U}\left(D,\frac{\delta G}{\delta
u}\right)\right\rangle+\left\llangle\left[\\!\\!\left[\frac{\delta F}{\delta
D}\mathcal{U}\left(D,\frac{\delta G}{\delta
u}\right)\right]\\!\\!\right],\tilde{D}\right\rrangle_{\Gamma}}$
$\displaystyle\quad{+\left\langle\nabla\frac{\delta G}{\delta
D},D\mathcal{U}\left(D,\frac{\delta G}{\delta
u}\right)\right\rangle-\left\llangle\left[\\!\\!\left[\frac{\delta G}{\delta
D}\mathcal{U}(D,\frac{\delta F}{\delta
u})\right]\\!\\!\right],\tilde{D}\right\rrangle_{\Gamma}}.$ (444)

The velocity upwinding used is a development of the Natale and Cotter (2018)
scheme, which does not require an auxiliary potential vorticity variable. To
achieve this, we apply the modification in the $q$ term of the Poisson
bracket,

$\left\langle\frac{\delta F}{\delta u},q\frac{\delta G}{\delta
u}\right\rangle\mapsto\left\langle\frac{\delta F}{\delta
u},\frac{1}{D}\nabla^{\perp}\cdot u\frac{\delta G}{\delta
u}^{\perp}\right\rangle+\left\langle\frac{\delta F}{\delta
u},\frac{f}{D}\frac{\delta G}{\delta u}^{\perp}\right\rangle.$ (445)

Following a discontinuous Galerkin methodology, we then integrate by parts
separately in each cell, and select $\tilde{u}$, the value of $u$ from the
upwind side, in the corresponding facet integral, to obtain

$-\left\langle u,\nabla^{\perp}\left(\frac{\delta F}{\delta
u}\cdot\frac{1}{D}\frac{\delta G}{\delta
u}^{\perp}\right)\right\rangle+\left\llangle\tilde{u},\left[\\!\\!\left[n^{\perp}\cdot\left(\frac{\delta
F}{\delta u}\cdot\frac{1}{D}\frac{\delta G}{\delta
u}^{\perp}\right)\right]\\!\\!\right]\right\rrangle_{\Gamma}.$ (446)

Since we are already using $\mathcal{U}$ in the layer depth terms, we might as
well avoid additional projections and use it to replace $\frac{\delta
G}{\delta u}/D$. Putting all of this together gives the

$\displaystyle\\{F,G\\}=$ $\displaystyle-\left\langle
u,\nabla^{\perp}\left(D\mathcal{U}\left(D,\frac{\delta F}{\delta
u}\right)\cdot\mathcal{U}\left(D,\frac{\delta G}{\delta
u}\right)^{\perp}\right)\right\rangle$
$\displaystyle\quad+\left\llangle\tilde{u},\left[\\!\\!\left[n^{\perp}\cdot\left(D\mathcal{U}\left(D,\frac{\delta
F}{\delta u}\right)\cdot\mathcal{U}\left(D,\frac{\delta G}{\delta
u}\right)^{\perp}\right)\right]\\!\\!\right]\right\rrangle_{\Gamma}$
$\displaystyle\quad{-\left\langle\nabla\frac{\delta F}{\delta
D},D\mathcal{U}\left(D,\frac{\delta G}{\delta
u}\right)\right\rangle+\left\llangle\left[\\!\\!\left[\frac{\delta F}{\delta
D}\mathcal{U}\left(D,\frac{\delta G}{\delta
u}\right)\right]\\!\\!\right],\tilde{D}\right\rrangle_{\Gamma}}$
$\displaystyle\quad{+\left\langle\nabla\frac{\delta G}{\delta
D},D\mathcal{U}\left(D,\frac{\delta F}{\delta
u}\right)\right\rangle-\left\llangle\left[\\!\\!\left[\frac{\delta G}{\delta
D}\mathcal{U}\left(D,\frac{\delta F}{\delta
u}\right)\right]\\!\\!\right],\tilde{D}\right\rrangle_{\Gamma}}.$ (447)

Now, if we use this Poisson bracket to generate the dynamical equations for
$D$ and $u$, we get

$\displaystyle\left\langle w,u_{t}\right\rangle-\left\langle
u,\nabla^{\perp}\left(D\mathcal{U}(D,w)\cdot u^{\perp}\right)\right\rangle$
$\displaystyle\qquad+\left\llangle\tilde{u},\left[\\!\\!\left[Dn^{\perp}\cdot\left(\mathcal{U}(D,w)\cdot
u^{\perp}\right)\right]\\!\\!\right]\right\rrangle_{\Gamma}$
$\displaystyle\quad{+\left\langle\nabla\frac{\delta H}{\delta
D},D\mathcal{U}(D,w)\right\rangle-\left\llangle\left[\\!\\!\left[\frac{\delta
H}{\delta
D}\mathcal{U}(D,w)\right]\\!\\!\right],\tilde{D}\right\rrangle_{\Gamma}}$
$\displaystyle=0,\quad\forall w\in\mathbb{V}_{h}^{1},$ (448)
$\displaystyle\left\langle\phi,D_{t}\right\rangle-\left\langle\nabla\phi,Du\right\rangle+\left\llangle\left[\\!\\!\left[\phi
u\right]\\!\\!\right],\tilde{D}\right\rrangle_{\Gamma}$
$\displaystyle=0,\quad\forall\phi\in\mathbb{V}_{h}^{2},$ (449)

which conserves energy by construction, despite the presence of the upwind
terms. As written, it is not clear how to implement this scheme, due to the
presence of the application of $\mathcal{U}$ to the test function $w$. It is
not scalable to compute and store $\mathcal{U}(D,w)$ for each test function
$w$, since $\mathcal{U}(D,w)$ is globally supported in general. Wimmer et al.
(2020) solved this problem by introducing an auxiliary variable
$r\in\mathbb{V}_{h}^{1}$ such that

$\displaystyle\left\langle r,Dv\right\rangle$ $\displaystyle=-\left\langle
u,\nabla^{\perp}\left(Dv\cdot u^{\perp}\right)\right\rangle$
$\displaystyle\quad+\left\llangle\tilde{u},\left[\\!\\!\left[Dn^{\perp}\cdot\left(v\cdot
u^{\perp}\right)\right]\\!\\!\right]\right\rrangle_{\Gamma},\quad\forall
v\in\mathbb{V}_{h}^{1}.$ (450)

Then,

$\displaystyle\left\langle r,w\right\rangle$ $\displaystyle=\left\langle
r,\mathcal{U}(D,w)D\right\rangle$ (451) $\displaystyle=-\left\langle
u,\nabla^{\perp}\left(D\mathcal{U}(D,w)\cdot u^{\perp}\right)\right\rangle$
$\displaystyle\quad+\left\llangle\tilde{u},\left[\\!\\!\left[Dn^{\perp}\cdot\left(\mathcal{U}(D,w)\cdot
u^{\perp}\right)\right]\\!\\!\right]\right\rrangle_{\Gamma},\quad\forall
w\in\mathbb{V}_{h}^{1},$ (452)

as required. Hence, we obtain the three coupled equations for
$(u,r,D,\frac{\delta H}{\delta
D})\in\mathbb{V}_{h}^{1}\times\mathbb{V}_{h}^{1}\times\mathbb{V}_{h}^{2}\times\mathbb{V}_{h}^{2}$,
where

$\displaystyle\left\langle w,u_{t}\right\rangle-\left\langle
Dw,r\right\rangle$ $\displaystyle=0,\forall w\in\mathbb{V}_{h}^{1},$ (453)
$\displaystyle\left\langle r,v\right\rangle-\left\langle
u,\nabla^{\perp}\left(Dv\cdot
u^{\perp}\right)\right\rangle+\left\llangle\tilde{u},\left[\\!\\!\left[Dn^{\perp}\cdot\left(v\cdot
u^{\perp}\right)\right]\\!\\!\right]\right\rrangle_{\Gamma}$
$\displaystyle\quad{+\left\langle\nabla\frac{\delta H}{\delta
D},Dv\right\rangle-\left\llangle\left[\\!\\!\left[\frac{\delta H}{\delta
D}v\right]\\!\\!\right],\tilde{D}\right\rrangle_{\Gamma}}$
$\displaystyle=0,\quad\forall v\in\mathbb{V}_{h}^{1},$ (454)
$\displaystyle\left\langle\phi,D_{t}\right\rangle-\left\langle\nabla\phi,Du\right\rangle+\left\llangle\left[\\!\\!\left[\phi
u\right]\\!\\!\right],\tilde{D}\right\rrangle_{\Gamma}$
$\displaystyle=0,\quad\forall\phi\in\mathbb{V}_{h}^{2},$ (455)
$\displaystyle\left\langle s,\frac{\delta H}{\delta
D}-\frac{1}{2}|u|^{2}-gD\right\rangle$ $\displaystyle=0,\quad\forall
s\in\mathbb{V}_{h}^{2}.$ (456)

Note that the use of upwinding in the $D$ equation has altered the
$\frac{\delta H}{\delta D}$ term in the $u$ equation, and it now appears in
terms other than an inner product with $\nabla\cdot w$. This means we have to
solve for the projection to $\mathbb{V}_{h}^{1}$ in an additional equation,
which complicates the solution of this system after discretisation with an
implicit Poisson integrator (which is necessary for exact energy conservation
for the fully discrete scheme). Wimmer et al. (2020) addressed the solution by
using a Picard iteration on $(u,D)$, keeping the standard linearisation about
a state of rest as an update equation and forming the nonlinear residual by
first computing $\frac{\delta H}{\delta D}$ and $r$ and substituting them into
the residuals for $u$ and $D$. The linear system can then be solved for using
hybridisation techniques as described in Section 5. Another possible approach
is to apply Newton’s method to the full four component system. Then, when
solving the Jacobian linear system for the update, $\delta r$ and
$\delta\frac{\delta H}{\delta D}$ can be eliminated as part of a Schur
complement preconditioner; the Schur complement in $u$ and $D$ can then be
approximated by the corresponding Schur complement arising from a more
standard discretisation such as those discussed in earlier sections.

Wimmer et al. (2020) demonstrated that this scheme combined with a Poisson
integrator produces relative energy conservation error of size $10^{-9}$ after
4 Picard iterations when applied to a standard test case used in the
development of numerical schemes for numerical weather prediction (the flow
over a mountain testcase 5 from Williamson et al. (1992)). In fact, similar
sized energy errors were even obtained when using the implicit midpoint rule
(also approximated using 4 Picard iterations), showing that the most important
aspect of energy conservation is in the spatial discretisation, at least for
this test problem. Additionally, if a midpoint rule is used then there is no
need for the additional auxiliary variable $r$ above. This is an example where
the reduction in energy error through space and time discretisation tricks is
lost due to the energy error from doing a finite number of nonlinear
iterations666Thanks to Golo Wimmer for pointing out this observation..

###  7.14 SUPG methods for active tracers

Wimmer et al. (2021) examined similar approaches when SUPG schemes are applied
to the advected quantities. The motivation for this is that an SUPG scheme is
required to stabilise the vertical transport of temperature when it is
approximated in the $\mathbb{W}_{h}^{\theta}$ space proposed in Section 2.
SUPG schemes have additional complications as they modify the test function
throughout the equation, including in the time derivative term. Wimmer et al.
(2021) addressed this problem as follows. First, we adopt the notation that
SUPG makes the modification $\gamma\mapsto\gamma+\tau S(u;\gamma)$ to the test
function $\gamma\in\mathbb{W}_{h}^{\theta}$. This general notation is to cover
different possibilities; for standard SUPG for continuous finite element
spaces the modification takes the form

$S(u;\gamma)=u\cdot\nabla\gamma.$ (457)

We note that $S(u;\gamma)$ will always be linear in the test function
$\gamma$, but we do not require linearity in the other variables. Then, if the
plain version of the discrete transport equation takes the form

$\left\langle\gamma,\theta_{t}\right\rangle+L\left(u,\theta;\gamma\right)=0,\quad\forall\gamma\in\mathbb{W}_{h}^{\theta},$
(458)

then the SUPG version takes the form

$\left\langle s+\tau S(u;s),\theta_{t}\right\rangle+L\left(u,\theta;s+\tau
S(u;s)\right)=0,\quad\forall s\in\mathbb{W}_{h}^{\theta}.$ (459)

For example, with a standard continuous Galerkin approximation, we would have

$L(u,\theta;\gamma)=-\left\langle\nabla\gamma,u\theta\right\rangle,\quad\forall\gamma\in\mathbb{W}_{h}^{\theta}.$
(460)

We then define $s(u;\gamma)\in\mathbb{W}_{h}^{\theta}$ according to

$\left\langle s(u;\gamma)+\tau
S(u;s(u,\gamma)),\sigma\right\rangle=\left\langle\gamma,\sigma\right\rangle,\quad\forall\sigma\in\mathbb{W}_{h}^{\theta}.$
(461)

Wimmer et al. (2021) proved the well posedness of this definition. Our SUPG
transport equation then becomes

$\left\langle\gamma,\theta_{t}\right\rangle+L\left(u,\theta;s(u;\gamma)+\tau
S(u;s(u;\gamma))\right)=0,\quad\forall\gamma\in\mathbb{W}_{h}^{\theta}.$ (462)

Following the principle of Gassmann and Herzog (2008) again, we substitute
this form into the relevant terms in the Poisson bracket, which are

$\displaystyle\\{F,G\\}=$
$\displaystyle\ldots-\left\langle\frac{1}{D}\frac{\delta
G}{\delta\theta}\nabla\theta,\frac{\delta F}{\delta
u}\right\rangle+\left\langle\frac{1}{D}\frac{\delta
F}{\delta\theta}\nabla\theta,\frac{\delta G}{\delta u}\right\rangle.$ (463)

Since we know that $\frac{\delta H}{\delta u}/D=u$ in the unapproximated case,
the substitution gives

$\displaystyle\\{F,G\\}=$ $\displaystyle\ldots-L\left(\frac{1}{D}\frac{\delta
G}{\delta u},\theta;s\left(u;\frac{\delta F}{\delta\theta}\right)+\tau
S\left(u;s\left(u;\frac{\delta F}{\delta\theta}\right)\right)\right)$
$\displaystyle\quad+L\left(\frac{1}{D}\frac{\delta F}{\delta
u},\theta;s\left(u;\frac{\delta G}{\delta\theta}\right)+\tau
S\left(u;s\left(u;\frac{\delta G}{\delta\theta}\right)\right)\right).$ (464)

This produces dynamical equations of the form

$\displaystyle\left\langle\gamma+\tau S\left(\frac{1}{D}\frac{\delta H}{\delta
u};\gamma\right),\theta_{t}\right\rangle+L\left(\frac{1}{D}\frac{\delta
H}{\delta u},\theta;\gamma+\tau S\left(\frac{1}{D}\frac{\delta H}{\delta
u};\gamma\right)\right)$
$\displaystyle=0,\quad\forall\gamma\in\mathbb{W}_{\theta},$ (465)
$\displaystyle\left\langle w,u_{t}\right\rangle+\ldots$
$\displaystyle\qquad-L\left(\frac{1}{D}w,\theta;s\left(u;\frac{\delta
H}{\delta\theta}\right)+\tau S\left(u;s\left(u;\frac{\delta
H}{\delta\theta}\right)\right)\right)$ $\displaystyle=0,\quad\forall
w\in\mathbb{W}_{h}^{1},$ (466)

where “$\ldots$” represents the terms coming from the other parts of the
Poisson bracket. The latter term in the $u_{t}$ equation is an approximation
of the term

$-c_{p}\theta\nabla\Pi,$ (467)

in the case of the compressible Euler equations. To replace $\frac{\delta
H}{\delta u}/D$ in the $\theta_{t}$ equation back with $u$ again, Wimmer et
al. (2021) also used the $\mathcal{U}$ operator as described above for the
shallow water equations, but we will not incorporate that additional
complexity here.

Wimmer et al. (2021) demonstrated robust upwind stabilisation combined with
energy conservation in various test problems using this method. For the
thermal shallow water equations, the upwinded version demonstrated much
smoother temperature fields than in the standard energy conserving version.
This is significant because the thermal shallow water equations exhibit very
fine structures in the temperature field which lead to the accumulation of
numerical noise at the gridscale, if upwinding is not used. For the
compressible Euler equations in a vertical slice formulation in a falling
bubble configuration, they showed that the energy conserving form of the
upwinded scheme leads to the appearance of secondary Kelvin Helmholtz vortices
that appear in much higher resolution simulations of the same problem (but do
not appear with upwinding schemes on the same resolution that do not conserve
energy). This latter result seems to suggest that the energy conserving
formulation is transferring potential energy dissipated from the upwind
transport scheme and injecting it into the kinetic energy in a manner that is
consistent with subscale processes in the higher resolution simulation. Wimmer
et al. (2021) also developed formulae that showed that the energy conservation
is indeed maintained by the transfer of dissipated potential energy into
kinetic energy. Since the potential energy dissipation occurs at the
gridscale, this raises the concern that the energy conservation leads to the
production of noise in the velocity field. By careful measurement of the
gridscale component of the velocity field in numerical experiments, it was
shown that the energy injection is at a larger scale (similar to what was
observed by Natale and Cotter (2017)). The naïve intuition that the energy
conserving formulation balances upwinding diffusion with antidiffusion in the
velocity term is incorrect, as the additional terms are not second order
derivatives in velocity.

##  8 Consistent vorticity and potential vorticity transport