CMlargesymbols”00 CMlargesymbols”01 Compatible finite element methods for
geophysical fluid dynamics–References

# Compatible finite element methods for geophysical fluid dynamics

Colin J. Cotter
Department of Mathematics  Imperial College London
South Kensington Campus  London SW7 2AZ
United Kingdom of Great Britain and Northern Ireland
E-mail<EMAIL_ADDRESS>
###### Abstract

This article surveys research on the application of compatible finite element
methods to large scale atmosphere and ocean simulation. Compatible finite
element methods extend Arakawa’s C-grid finite difference scheme to the finite
element world. They are constructed from a discrete de Rham complex, which is
a sequence of finite element spaces which are linked by the operators of
differential calculus. The use of discrete de Rham complexes to solve partial
differential equations is well established, but in this article we focus on
the specifics of dynamical cores for simulating weather, oceans and climate.
The most important consequence of the discrete de Rham complex is the Hodge-
Helmholtz decomposition, which has been used to exclude the possibility of
several types of spurious oscillations from linear equations of geophysical
flow. This means that compatible finite element spaces provide a useful
framework for building dynamical cores. In this article we introduce the main
concepts of compatible finite element spaces, and discuss their wave
propagation properties. We survey some methods for discretising the transport
terms that arise in dynamical core equation systems, and provide some example
discretisations, briefly discussing their iterative solution. Then we focus on
the recent use of compatible finite element spaces in designing structure
preserving methods, surveying variational discretisations, Poisson bracket
discretisations, and consistent vorticity transport.

###### doi:

XXXXXXXX

###### Contents

  1. 1 Introduction
  2. 2 Compatible finite element spaces
  3. 3 Wave propagation properties
  4. 4 Transport and stabilisation
  5. 5 Example discretisations and iterative solution strategies
  6. 6 Variational discretisations
  7. 7 Almost-Poisson brackets
  8. 8 Consistent vorticity and potential vorticity transport
  9. 9 Structure preserving schemes on non-affine meshes
  10. 10 Summary and outlook

##  1 Introduction

Atmosphere and ocean models used in weather forecasting and climate simulation
are built around dynamical cores, which predict the fundamental quantities of
fluid motion: fluid velocity, pressure, density, temperature, and, in the case
of the ocean, salinity. Dynamical cores are computer implementations of
numerical discretisations of partial differential equation models of
geophysical fluid dynamics in the absence of viscosity. These models are then
coupled with physics parameterisations that describe additional physics as
well as fluid dynamical processes involving scales that are too small and fast
to represent explicitly in the dynamical core. In the case of atmosphere
models, this can include radiation processes, cloud models, moisture and
precipitation in various phases, models of unresolved convection, boundary
layers, and momentum transfer due to unresolved internal waves. In the case of
ocean models, this includes vertical mixing due to convection, and
parameterisation of unresolved eddy motions. Atmosphere and ocean models can
be coupled together using parameterisations of air-sea interaction processes,
and with other process models to form climate models. Other process models
include models of land albedo, land ice, sea ice, atmospheric chemistry etc.
Additionally, in operational forecasting systems, atmosphere and ocean models
are blended with observed data using data assimilation algorithms.

A very important aspect of global atmosphere and ocean models is the wide
range of timescales. In the rotating compressible Euler equations that model
atmospheric flow (although usually with additional approximations), high
frequency motions are possible in the form of pressure waves and internal
waves. However, in the solutions relevant to atmosphere models these
frequencies are observed to have very low amplitude, with the velocity
remaining very close to being horizontal and divergence-free. Since the
numerical solution is never anywhere near being resolved in an atmosphere
model, and due to the coupling with all of the other modelled processes
described above, plus the modifications made to the solution by the data
assimilation algorithms, it is critically important that unphysical large
amplitude wave motions are avoided. Similar issues arise in the ocean model
context. All of this means that it is critically important to pay attention to
what happens in the dynamical core at the gridscale, to avoid numerical error
triggering these motions. These concerns lie at the heart of the decisions
about which discretisations to use in building a dynamical core.

In the past decades, atmosphere and ocean dynamical cores have been mainly
built around finite difference, finite volume and pseudospectral methods. In
this article, we describe a more recent approach to building dynamical cores
using compatible finite element methods. Finite element methods have the
benefit that they do not depend strongly on the structure of the underlying
mesh for their consistency and rate of convergence under mesh refinement. They
can use polygonal cells of different types, with triangulations in the
horizontal being particularly flexible. This allows mesh refinement in regions
of focus, and the constuction of meshes that conform to coastlines and areas
of high topography curvature. It also allows adaptive mesh refinement although
the advantages of this are less clear in large scale weather and climate
simulation.

The relaxation of mesh structure also allows for more uniform meshes in
spherical geometry. Instead of using a latitude longitude grid, which has very
thin cells at the poles due to the convergence of lines of latitude, a cubed
sphere mesh (quadrilateral refinement of a cube mapped to a sphere) or an
icosahedral mesh (triangular refinement of an icosahedron mapped to a sphere)
can be used. These pseudouniform grids avoid the parallel scalability
bottlenecks (Zängl et al. 2015, Adams et al. 2019).

Finally, finite element methods allow for increasing the degree of polynomials
used in each cell, leading to higher order accuracy. Higher order finite
element methods have a denser, more structured data layout in each cell, which
can be exploited to try to achieve a higher computational intensity, doing
more computational work whilst fetching data from nearby cells (Dennis et al.
2012, Giraldo et al. 2013).

Compatible finite elements, the subject of this survey, address the problem of
spurious numerical waves that cause problems when coupled with numerical
errors, physics parameterisations and data assimilation schemes at the
gridscale (Staniforth and Thuburn 2012). They can be seen as an extension of
the C grid finite difference method (Arakawa and Lamb 1977), which avoids many
spurious wave issues when used with quadrilateral cells (and slightly less so
with triangular cells). The C grid method achieves this by placing different
fields (velocity components, pressure, _etc,_) on different entities of the
grid (cells, edges, vertices) so that the discretised vector calculus
operators (div, grad, curl) retain kernels of the appropriate size. This was
formalised as the Discrete Exterior Calculus (DEC) (Hirani 2003). Compatible
finite element methods do the same thing, but at the level of the finite
element spaces and the vector calculus mappings between them. These methods
have been unified in the Finite Element Exterior Calculus (Arnold et al.
2006)111The author uses the term “compatible finite elements” as a way to
discuss them with practitioners when using the language of vector calculus
instead of differential forms, which are more unifying but require more
background material to discuss., with a long history going back to the 1970s
of applications to porous media, elasticity and fluid dynamics. Their
application to geophysical fluid dynamics also has a long history, especially
in the ocean (Le Roux et al. 1998, Walters 2005), and their representation of
exact geostrophic balance has been well known from numerical dispersion
analysis for some time (Le Roux et al. 2007, Rostand and Le Roux 2008).
However, it was Cotter and Shipton (2012) that first noticed that it was the
compatible structure that was behind this property, and that this structure
can be used to understand the numerical dispersion properties of compatible
finite element schemes. Following earlier work by the C grid and DEC
communities, compatible finite element methods have also been used to build
structure preserving discretisations that embed conservation laws at the
discrete level.

This article will introduce compatible finite element methods for geophysical
fluid dynamics and their properties, discuss how to build atmosphere and ocean
dynamical cores out of them, before focussing on work on structure preserving
discretisations. The rest of the article is organised as follows. Section 2
will review compatible finite element spaces and their fundamental properties.
Section 3 will discuss their application to linear wave models in geophysical
fluid dynamics and the use of the compatible structure to understand discrete
wave propagation properties. Moving towards nonlinear dynamical cores, Section
4 will discuss the discretisation of transport schemes for compatible finite
elements, and Section 5 will use them to describe discretisations for
nonlinear geophysical fluid dynamics models and survey their use. One common
theme is that the analysis of the convergence and stability of these schemes
is very scarse and there is a lot of opportunity for future research in this
area. The next four sections concentrate on the use of compatible finite
element methods in structure preserving discretisations of various types.
Section 6 surveys methods constructed from a discrete Hamilton’s Principle.
Section 7 covers methods constructed from almost Poisson brackets that ensure
conservation of energy and other quantities. Section 8 discusses methods that
embed a diagnosed vorticity or potential vorticity field with its own
consistent dynamics. Finally, Section 9 discusses some issues common to all
structure preserving methods related to nonaffine meshes such as meshes
approximating the sphere as well as terrain following meshes. We then end the
survey with a very brief summary in Section 10.

##  2 Compatible finite element spaces

In this section we describe compatible finite element spaces as they are used
in geophysical fluid dynamics applications. Compatible finite element spaces
are sequences of spaces that form a discrete differential complex (which we
shall describe below). In applications to geophysical fluid dynamics, the
focus is on the de Rham complex. There is a growing body of work for
compatible finite elements spaces built around the Stokes complex (Tai and
Winther 2006, Falk and Neilan 2013, Neilan 2020, Hu et al. 2022, for example),
the elasticity complex, and the Regge complex for general relativity
(Christiansen 2011, Christiansen et al. 2020, for example), but we shall not
discuss those spaces here.

###  2.1 Preliminary notation

First we establish some notation. Having defined the usual space
$L^{2}(\Omega)$ of square integrable scalar functions on $\Omega$ (and writing
$L^{2}(\Omega)^{N}$ as the space of vector functions on $\Omega$ in dimension
$N$), for a domain $\Omega$ in three dimensions we have

$\displaystyle H^{1}(\Omega)$ $\displaystyle=\left\\{\phi\in
L^{2}(\Omega):\nabla\phi\in L^{2}(\Omega)\right\\},$ (1) $\displaystyle
H(\operatorname{curl};\Omega)$ $\displaystyle=\left\\{u\in
L^{2}(\Omega)^{3}:\nabla\times u\in L^{2}(\Omega)^{3}\right\\},$ (2)
$\displaystyle H(\operatorname{div};\Omega)$ $\displaystyle=\left\\{u\in
L^{2}(\Omega)^{3}:\nabla\cdot u\in L^{2}(\Omega)\right\\},$ (3)

where $\nabla$, $\nabla\cdot$ and $\nabla\times$ are the appropriately defined
weak derivative operators (see a textbook on analysis of PDEs for details, or
proceed just assuming that these are the usual derivatives for smooth
functions and the cellwise derivatives for finite element functions with the
appropriate continuity, as we shall discuss below). These spaces are
accompanied by norms defined as

$\displaystyle\|\phi\|^{2}_{L^{2}(\Omega)}$
$\displaystyle=\int_{\Omega}\phi^{2}\operatorname{d}x,$ (4)
$\displaystyle\|u\|^{2}_{L^{2}(\Omega)^{N}}$
$\displaystyle=\int_{\Omega}|u|^{2}\operatorname{d}x,$ (5)
$\displaystyle\|\phi\|^{2}_{H^{1}(\Omega)}$
$\displaystyle=\|\phi\|^{2}_{L^{2}(\Omega)}+\|\nabla\phi\|^{2}_{L^{2}(\Omega)^{N}},$
(6) $\displaystyle\|u\|^{2}_{H(\operatorname{curl};\Omega)}$
$\displaystyle=\|u\|^{2}_{L^{2}(\Omega)^{N}}+\|\nabla\times
u\|^{2}_{L^{2}(\Omega)^{N}},$ (7)
$\displaystyle\|u\|^{2}_{H(\operatorname{div};\Omega)}$
$\displaystyle=\|u\|^{2}_{L^{2}(\Omega)^{N}}+\|\nabla\cdot
u\|^{2}_{L^{2}(\Omega)}.$ (8)

In two dimensions, we similarly define

$\displaystyle H(\operatorname{curl};\Omega)$ $\displaystyle=\left\\{u\in
L^{2}(\Omega)^{2}:\nabla^{\perp}\cdot u\in L^{2}(\Omega)\right\\},$ (9)
$\displaystyle H(\operatorname{div};\Omega)$ $\displaystyle=\left\\{u\in
L^{2}(\Omega)^{2}:\nabla\cdot u\in L^{2}(\Omega)\right\\},$ (10)

having defined the operators
$\nabla^{\perp}\phi=(-\partial_{x_{2}}\phi,\partial_{x_{1}}\phi)$ and
$\nabla^{\perp}\cdot u=-\partial_{x_{2}}u_{1}+\partial_{x_{1}}u_{2}$.
Frequently we drop the $\Omega$ from this notation when the meaning is clear.
We will also use the $L^{2}$ inner product notation,

$\displaystyle\left\langle p,q\right\rangle$
$\displaystyle=\int_{\Omega}pq\operatorname{d}x,\,\forall p,q\in
L^{2}(\Omega),$ (11) $\displaystyle\left\langle u,v\right\rangle$
$\displaystyle=\int_{\Omega}u\cdot v\operatorname{d}x,\,\forall u,v\in
L^{2}(\Omega)^{N},$ (12)

noting that

$\|p\|_{L^{2}(\Omega)}^{2}=\left\langle
p,p\right\rangle,\quad\|u\|_{L^{2}(\Omega)^{N}}^{2}=\left\langle
u,u\right\rangle.$ (13)

We also define $L^{2}(\Omega)$ inner products for functions on boundaries,

$\displaystyle\left\llangle p,q\right\rrangle$
$\displaystyle=\int_{\partial\Omega}pq\operatorname{d}S,\,\forall p,q\in
L^{2}(\partial\Omega),$ (14) $\displaystyle\left\llangle u,v\right\rrangle$
$\displaystyle=\int_{\partial\Omega}u\cdot v\operatorname{d}S,\,\forall u,v\in
L^{2}(\partial\Omega)^{N}.$ (15)

###  2.2 Discrete de Rham complexes

In three dimensions, discrete de Rham complexes on a domain $\Omega$ consist
of subspaces $\mathbb{W}_{0}\subset H_{1}(\Omega)$, $\mathbb{W}_{1}\subset
H(\operatorname{curl};\Omega)$, $\mathbb{W}_{2}\subset
H(\operatorname{div};\Omega)$, $\mathbb{W}_{3}\subset L^{2}(\Omega)$,
satisfying the following commutation relations,

$\begin{CD}\mathbb{W}^{0}={H}^{1}@>{\operatorname{d}^{1}=\nabla}>{}>\mathbb{W}^{1}={H}(\operatorname{curl})@>{\operatorname{d}^{2}=\nabla\times}>{}>\mathbb{W}^{2}={H}(\textrm{div})@>{\operatorname{d}^{3}=\nabla\cdot}>{}>L^{2}\\\
@V{}V{\pi_{0}}V@V{}V{\pi_{1}}V@V{}V{\pi_{2}}V@V{}V{\pi_{3}}V\\\
{\mathbb{W}}^{0}_{h}@>{\operatorname{d}^{1}=\nabla}>{}>{\mathbb{W}}^{1}_{h}@>{\operatorname{d}^{2}=\nabla\times}>{}>\mathbb{W}^{2}_{h}@>{\operatorname{d}^{3}=\nabla\cdot}>{}>\mathbb{W}^{3}_{h},\\\
\end{CD}$ (16)

where $\pi_{i}$, $i=0,1,2,3$ are surjective projections, satisfying a bound
appropriate to their domain, e.g.

$\|\pi_{1}u\|_{H(\operatorname{curl})}\leq
C\|u\|_{H(\operatorname{curl})},\,\forall u\in H(\operatorname{curl}).$ (17)

The commutation property means that
$\operatorname{d}^{k+1}\pi_{k}u=\pi_{k+1}\operatorname{d}^{k+1}u$, $k=0,1,2$.
These projection operators (see Arnold et al. (2006) for a guide to their
construction) do not play a role in computations, but they ensure that the
finite element spaces are compatible in the sense that the differential
operators $\nabla$, $\nabla\times$ and $\nabla\cdot$ map surjectively onto the
kernel of the next operator in the sequence, just as is the case for the
infinite dimensional spaces at the top of the diagram.

In two dimensions, the de Rham complex is shorter, and there are two possible
ways to write it,

$\displaystyle\begin{CD}\mathbb{V}^{0}={H}^{1}@>{\operatorname{d}^{1}=\nabla^{\perp}}>{}>\mathbb{V}^{1}={H}(\operatorname{div})@>{\operatorname{d}^{2}=\nabla\cdot}>{}>\mathbb{V}^{2}=L^{2}\\\
@V{}V{\pi_{0}}V@V{}V{\pi_{1}}V@V{}V{\pi_{2}}V\\\
{\mathbb{V}}^{0}_{h}@>{\operatorname{d}^{1}=\nabla^{\perp}}>{}>{\mathbb{V}}^{1}_{h}@>{\operatorname{d}^{2}=\nabla\cdot}>{}>\mathbb{V}^{2}_{h},\\\
\end{CD}$ (18)
$\displaystyle\begin{CD}\tilde{\mathbb{V}}^{0}={H}^{1}@>{\operatorname{d}^{1}=\nabla}>{}>\tilde{\mathbb{V}}^{1}={H}(\operatorname{curl})@>{\operatorname{d}^{2}=\nabla^{\perp}\cdot}>{}>\tilde{\mathbb{V}}^{2}=L^{2}\\\
@V{}V{\pi_{0}}V@V{}V{\pi_{1}}V@V{}V{\pi_{2}}V\\\
\tilde{{\mathbb{V}}}^{0}_{h}<EMAIL_ADDRESS>(19)

The equivalence stems from the fact that any vector field in
$H(\operatorname{div};\Omega)$ in two dimensions can be transformed into a
vector field in $H(\operatorname{curl};\Omega)$ by rotating it by $\pi/2$. In
compatible finite element methods for geophysical fluid dynamics we tend to
use (18), since it allows for local mass conservation and exact application of
flux boundary conditions.

###  2.3 Discrete Hodge-Helmholtz decomposition

Crucially for the geophysical fluid dynamics setting, the bounded, commuting,
surjective projections ensure a discrete version of the Hodge-Helmholtz
decompositions. At the infinite dimensional level, these decompositions are

$\mathbb{W}^{k}=B^{k}\oplus\mathfrak{h}^{k}\oplus(B^{*})^{k},$ (20)

where

$\displaystyle B^{k}$
$\displaystyle=\left\\{u\in\mathbb{W}^{k}:\exists\phi\in\mathbb{W}^{k-1},\,\mbox{with}\,u=\operatorname{d}^{k}\phi\right\\},$
(21) $\displaystyle\mathfrak{h}^{k}$
$\displaystyle=\left\\{u\in\mathbb{W}^{k}:\operatorname{d}^{k+1}u=0,\,\delta^{k}u=0\right\\},$
(22) $\displaystyle(B^{*})^{k}$
$\displaystyle=\left\\{u\in\mathbb{W}^{k}:\exists\phi\in\mathbb{W}^{k+1},\,\mbox{with}\,u=\delta^{k}\phi\right\\},$
(23)

defining the dual operator $\delta^{k}:\mathbb{W}^{k}\to\mathbb{W}^{k-1}$ such
that

$\left\langle\phi,\delta^{k}u\right\rangle=-\left\langle\operatorname{d}^{k+1}\phi,u\right\rangle.$
(24)

When $u\in\mathbb{W}^{k}$ is appropriately constructed so that boundary
integrals vanish under integration by parts (and is sufficiently smooth that
integration by parts is well defined) then we can make the following
associations,

$\delta^{0}=-\nabla\cdot,\,\delta^{1}=\nabla\times,\,\delta^{2}=-\nabla.$ (25)

For example, considering $\delta^{1}$, we have

$\left\langle v,\delta^{1}u\right\rangle=\left\langle\nabla\times
v,u\right\rangle=\left\langle v,\nabla\times u\right\rangle-\left\llangle
v,n\times u\right\rrangle,$ (26)

where $n$ is the outward pointing normal to $\partial\Omega$. If we choose
$u\in\mathring{\mathbb{W}}^{2}$, where

$\mathring{\mathbb{W}}^{2}=\left\\{u\in\mathring{\mathbb{W}}^{2}:u\cdot
n=0\mbox{ on }\partial\Omega\right\\},$ (27)

then we obtain that $\delta^{1}u=\nabla\times u$ for
$u\in\mathring{\mathbb{W}}^{2}\cap\mathbb{W}^{1}$. Similarly, we define

$\displaystyle\mathring{\mathbb{W}}^{0}$
$\displaystyle=\left\\{\phi\in\mathbb{W}^{0}:\phi=0\mbox{ on
}\partial\Omega\right\\},$ (28) $\displaystyle\mathring{\mathbb{W}}^{1}$
$\displaystyle=\left\\{u\in\mathbb{W}^{1}:u\times n=0\mbox{ on
}\partial\Omega\right\\},$ (29) $\displaystyle\mathring{\mathbb{W}}^{3}$
$\displaystyle=\mathbb{W}^{0}.$ (30)

See Chapter 3 of Arnold (2018) for an accessible description of the full
functional setting for these aspects.

Returning to (20), we explain the $\oplus$ notation. It means that any
$u\in\mathbb{W}^{k}$ can be written uniquely as $b+h+c$ with $b\in B^{k}$,
$h\in\mathfrak{h}^{k}$, and $c\in(B^{*})^{k}$. In fact, further we have that
the spaces $B^{k}$, $\mathfrak{h}^{k}$ and $(B^{*})^{k}$ are mutually
orthogonal under the $L^{2}$ inner product. When $k=2$, we recognise this as
the Helmholtz decomposition for vector fields, which says that a vector field
can be written uniquely as $u=\nabla\times v+h+\nabla\phi$, with $\nabla\cdot
h=0$ and $\nabla\times h=0$; h is referred to as a “harmonic vector field”.
This decomposition is crucial to the understanding of rapidly rotating fluid
dynamics.

As described in Arnold et al. (2006), three dimensional compatible finite
element spaces satisfy a discrete Hodge-Helmholtz decomposition,

$\mathbb{W}^{k}_{h}=B^{k}_{h}\oplus\mathfrak{h}^{k}_{h}\oplus(B^{*})^{k}_{h},$
(31)

where

$\displaystyle\mathfrak{h}^{k}_{h}$
$\displaystyle=\left\\{u\in\mathbb{W}^{k}_{h}:\operatorname{d}^{k+1}u=0,\,\delta^{k}_{h}u=0\right\\},$
(32) $\displaystyle(B^{*})^{k}_{h}$
$\displaystyle=\left\\{u\in\mathbb{W}^{k}_{h}:\exists\phi\in\mathbb{W}^{k+1}_{h},\,\mbox{with}\,u=\delta^{k}_{h}\phi\right\\},$
(33)

and we have the discrete dual operator
$\delta^{k}_{h}:\mathbb{W}^{k+1}_{h}\to\mathbb{W}^{k}_{h}$, such that

$\left\langle\phi,\delta^{k}_{h}u\right\rangle=\left\langle\operatorname{d}^{k+1}\phi,u\right\rangle,\quad\forall\phi\in\mathbb{W}^{k}_{h},u\in\mathbb{W}^{k+1}_{h}.$
(34)

We note the asymmetry under discretisation: $\operatorname{d}^{k}$ has the
same definition at the discrete level, but $\delta^{k}$ is replaced by the
approximation $\delta^{k}_{h}$.

The discrete Hodge-Helmholtz decomposition inherits some important properties
from the infinite dimensional decomposition. First, we have $B^{k}_{h}\subset
B^{k}$. Second, although $\mathfrak{h}^{k}_{h}\neq\mathfrak{h}^{k}$, we do
have that $\dim(\mathfrak{h}^{k}_{h})=\dim(\mathfrak{h}^{k})$. Further,
$\mathfrak{h}^{k}_{h}$ converges to $\mathfrak{h}^{k}$ as the mesh is refined.
These two properties have important consequences for the correct
representation of inertial oscillations in geophysical fluid dynamics models.
We do not have $(B^{*})^{k}_{h}\subset(B^{*})^{k}$.

For later discussion we also define

$\zeta^{k}=\left\\{u\in\mathbb{W}^{k}:d^{k+1}u=0\right\\},$ (35)

noting that $\zeta^{k}=B^{k}\oplus\mathfrak{h}^{k}$. Similarly, we write
$\zeta^{k}_{h}=B^{k}_{h}\oplus\mathfrak{h}^{k}_{h}$.

The two dimensional discrete de Rham complexes also provide an analogous
Helmholtz-Hodge decomposition for $\mathbb{V}_{h}^{k}$, which is most
significant for $\mathbb{V}_{h}^{1}$.

###  2.4 Two dimensional compatible finite element spaces used in geophysical
fluid dynamics

For compatible finite element methods for geophysical fluid dynamics, we are
mainly focussed on the goal of producing three dimensional fluid models of the
atmosphere and ocean. However, model development usually starts by
consideration of the rotating shallow water equations, since these equations
encompass many of the challenges of designing numerical schemes for atmosphere
and ocean but without the additional computational challenges of three
dimensional models. On quadrilateral meshes (such as the “cubed sphere”
meshing of the sphere), the most commonly used spaces are
$\mathbb{V}_{h}^{0}=Q_{k}$ (tensor product Lagrange elements, e.g. bilinear,
biquadratic, _etc_), $\mathbb{V}_{h}^{1}=RT_{k-1}$ (Raviart-Thomas elements on
quadrilaterals, noting here that we use the traditional numbering convention
according to the largest complete polynomial space contained by the finite
element, not the UFL/FIAT numbering according to the highest degree polynomial
in the space), and $\mathbb{V}_{h}^{2}=DQ_{1}$ (discontinuous tensor product
Lagrange elements). However, it may be interesting to consider the trimmed
serendipity family of elements, especially at higher order Gillette and
Kloefkorn (2019). On triangular meshes (which are much more flexible in
allowing local mesh refinement), the main possibilities are
$\mathbb{V}_{h}^{0}=P_{k+1}$ (Lagrange elements), $\mathbb{V}_{h}^{1}=BDM_{k}$
(Brezzi-Douglas-Marini elements on triangles), and
$\mathbb{V}_{h}^{2}=DG_{k-1}$, or $\mathbb{V}_{h}^{0}=P_{k+1}$,
$\mathbb{V}_{h}^{1}=RT_{k}$ (Raviart-Thomas elements on triangles), and
$\mathbb{V}_{h}^{0}=DG_{k}$.Whilst in most applications, the second grouping
is preferred, because it requires less degrees of freedom for the same
accuracy, in geophysical fluid dynamics this grouping has wave propagation
issues related to the Coriolis term, and we tend to prefer the first grouping
based on $BDM$ elements (we shall discuss this later). In all of these
discrete de Rham complexes, we see the reduction in interelement continuity
properties moving across the discrete de Rham complex: $\mathbb{V}_{h}^{0}$
contains only continuous functions, $\mathbb{V}_{h}^{1}$ requires continuity
of normal components its vector valued functions across cell facets (but not
tangential components), and $\mathbb{V}_{h}^{2}$ has functions with no
interelement continuity constraints. See Boffi et al. (2013) for definitions
of the finite elements introduced in this paragraph.

###  2.5 Three dimensional compatible finite element spaces used in
geophysical fluid dynamics

When moving to three dimensional models, since the Earth’s atmosphere and
ocean are much larger in horizontal extent than the vertical, good
preservation of hydrostatic balance (balance between gravitational
acceleration and vertical pressure gradient) requires the use of prismatic
meshes that are constructed by extruding a two dimensional base mesh into
layers. When the base mesh is constructed from quadrilaterals, this produces
hexahedra, and when the base mesh is constructed from triangles, this produces
triangular prisms. Three dimensional discrete de Rham complexes are
constructed on these spaces by a tensor product of a two dimensional de Rham
complex on the base mesh with a one dimensional de Rham complex for the
vertical direction, given by

$\begin{CD}\mathbb{U}^{0}={H}^{1}@>{\operatorname{d}^{1}=\partial_{x}}>{}>\mathbb{V}^{1}=L^{2}\\\
@V{}V{\pi_{0}}V@V{}V{\pi_{1}}V\\\
<EMAIL_ADDRESS>$ (36)

Excluding splines _etc._ , the main family of options for these one
dimensional spaces is continuous Lagrange elements of degree $k+1$ for
$\mathbb{U}_{h}^{0}$ and discontinuous Lagrange elements of degree $k$ for
$\mathbb{U}_{h}^{1}$. The 3D discrete de Rham complex is then formed on the
reference cell $\hat{K}_{3}=\hat{K}_{2}\times\hat{K}_{1}$ follows,

$\displaystyle\mathbb{W}_{h}^{0}(\hat{K}_{3})$
$\displaystyle=\mathbb{V}_{h}^{0}(\hat{K}_{2})\otimes\mathbb{U}_{h}^{0}(\hat{K}_{1}),$
(37) $\displaystyle\mathbb{W}_{h}^{1}(\hat{K}_{3})$
$\displaystyle=\underbrace{\hat{k}\mathbb{V}_{h}^{0}(\hat{K}_{2})\otimes\mathbb{U}_{h}^{1}(\hat{K}_{1})}_{=\mathbb{W}_{h}^{1,V}}\oplus\underbrace{\iota(\mathbb{V}_{h}^{1}(\hat{K}_{2}))^{\perp}\otimes\mathbb{U}_{h}^{0}(\hat{K}_{1})}_{=\mathbb{W}_{h}^{1,H}},$
(38) $\displaystyle\mathbb{W}_{h}^{2}(\hat{K}_{3})$
$\displaystyle=\underbrace{\mathbb{V}_{h}^{2}(\hat{K}_{2})\otimes\hat{k}\mathbb{U}_{h}^{0}(\hat{K}_{1})}_{=\mathbb{W}_{h}^{2,V}}\oplus\underbrace{\iota(\mathbb{V}_{h}^{1}(\hat{K}_{2}))\otimes\mathbb{U}_{h}^{1}(\hat{K}_{1})}_{=\mathbb{W}_{h}^{2,H}},$
(39) $\displaystyle\mathbb{W}_{h}^{3}(\hat{K}_{3})$
$\displaystyle=\mathbb{V}_{h}^{2}(\hat{K}_{2})\otimes\mathbb{U}_{h}^{1}(\hat{K}_{1}).$
(40)

Here we write coordinates on $\hat{K}_{3}$ as $(x,z)$ with $x\in\hat{K}_{2}$
and $z\in\hat{K}_{1}$, defining the tensor product $V(K_{2})\otimes U(K_{1})$
as the span of function products $u(x)v(z)$ with $u\in V(K_{2})$ and $v\in
V(K_{1})$. Further, $\hat{k}$ is the unit upward pointing vector, $\iota$ is
the inclusion operator that maps 2D vectors into equivalent 3D vectors with
zero vertical part, and $\perp$ is the operator that rotates vectors by a
quarter of a rotation in the horizontal direction. These latter technicalities
involving $\iota$, $\hat{k}$ and $\perp$ can all be avoided in the unified
presentation of spaces of discrete differential forms in the finite element
exterior calculus, which also unifies many other aspects across dimensions,
numbering of the spaces, _etc._ See Arnold and Awanou (2014) for a full
presentation of tensor product discrete differential forms.

Again, in all of these three dimensional discrete de Rham complexes, we see
the reduction in interelement continuity properties moving across the discrete
de Rham complex: $\mathbb{W}_{h}^{0}$ contains only continuous functions,
$\mathbb{W}_{h}^{1}$ requires continuity of tangential components its vector
valued functions across cell facets (but not normal components),
$\mathbb{W}_{h}^{2}$ requires continuity of normal components its vector
valued functions across cell facets (but not tangential components), and
$\mathbb{W}_{h}^{3}$ has functions with no interelement continuity
constraints.

As indicated in (38-39). the $\mathbb{W}_{h}^{1}$ and $\mathbb{W}_{h}^{2}$
spaces can be split into vertical and horizontal parts, indicated by the $V$
and $H$ superscripts respectively. As suggested by this notation, the vertical
part contains vector fields that point in the vertical direction, whilst the
horizontal part contains vector fields that point in the horizontal direction.
After Piola transformation to mesh cells (discussed in Section 2.6), this
decomposition is not preserved in general. However, if mesh cells are arranged
into a global tensor product mesh (i.e, flat vertical layers), then this
decomposition is preserved. This will also occur in a spherical annulus meshed
by radially extruding a 2D surface mesh of the sphere (or an approximation of
one). If a terrain following mesh is used, so that the side walls of the mesh
are arranged vertically, but the horizontal layers move up and down to conform
to mountain ranges on the surface, the decomposition is only partially
preserved: $\mathbb{W}_{h}^{1,H}$ remains tangential to the “up” direction,
and $\mathbb{W}_{h}^{2,V}$ remains normal to the “up” direction, but
$\mathbb{W}_{v}^{1,V}$ will contain some horizontal component and
$\mathbb{W}_{h}^{2,V}$ will contain some vertical component.

In discretisations for geophysical fluid dynamics it is important that the
gravity term can be balanced by a vertical pressure gradient. To avoid
degeneracy in this balance, discussed in Section 3.8, it is necessary to use a
finite element space $\mathbb{W}_{\theta}$ for temperature variables (entropy,
temperature, potential temperature, _etc_.) that is adapted to the vertical
part of $\mathbb{W}_{h}^{2}$, which is used to represent the velocity in this
framework. We choose $\mathbb{W}^{\theta}(\hat{K})$ so that
$\theta\in\mathbb{W}^{\theta}(\hat{K})\implies\hat{k}\theta\in\mathbb{W}_{h}^{2,V}(\hat{K})$,
_i.e_.,

$\mathbb{W}_{\theta}(\hat{K})=\mathbb{V}_{h}^{0}(\hat{K}_{2})\otimes\mathbb{U}_{h}^{1}(\hat{K}_{1}).$
(41)

It is also useful to construct vertical slice models by making analogous
tensor product constructions combining the
$(\mathbb{U}_{h}^{0},\mathbb{U}_{h}^{1})$ de Rham complex with itself. Since
the only possible 1D meshes are intervals, this just leads to the usual tensor
product elements for quadrilaterals that we have already discussed above.

For further description of the construction and efficient implementation of
these tensor product elements within an automated system, see McRae et al.
(2016). The use of extruded meshes also has computational benefits that offset
the additional computational cost of using unstructured meshes in the
horizontal. If a semistructured data layout is used (i.e. a horizontal
unstructured index and a vertical structured one) then for a reasonable number
of vertical layers (20 is already enough in numerical experiments) the lookups
to find data in the unstructured grid data structure are negligible compared
to the computational work done on data loaded into memory. Thus there is no
significant performance penalty to using an unstructured data structure in the
horizontal (Bercea et al. 2016). This means that the benefits of the
flexibility and mesh invariance of the unstructured grid data structure in the
horizontal can be exploited in performant 3D geophysical fluid models.

###  2.6 Local to global mappings

Whilst there has been work on e.g. H(div) elements on quadrilaterals where the
polynomials are defined directly on the mesh elements (Arbogast and Correa
2016), here we mostly restrict discussion to discrete de Rham complexes that
are constructed on reference elements and mapped to mesh elements using Piola
maps. In three dimensions, this corresponds to the following set of relations,

$\displaystyle\psi\in\mathbb{W}_{h}^{0}(K)\\!:\,$ $\displaystyle\psi\circ
g_{K}=\hat{\psi},\,\hat{\psi}\in\mathbb{W}^{0}_{h}(\hat{K}),$ (42)
$\displaystyle u\in\mathbb{W}_{h}^{1}(K)\\!:\,$ $\displaystyle u\circ
g_{K}=J^{-T}\hat{u},\,\hat{u}\in\mathbb{W}^{1}_{h}(\hat{K}),$ (43)
$\displaystyle w\in\mathbb{W}_{h}^{2}(K)\\!:\,$ $\displaystyle u\circ
g_{K}=J\hat{w}/\det J,\,\hat{w}\in\mathbb{W}^{2}_{h}(\hat{K}),$ (44)
$\displaystyle\phi\in\mathbb{W}_{h}^{3}(K)\\!:\,$ $\displaystyle\phi\circ
g_{K}=\hat{\phi}/\det J,\,\hat{\phi}\in\mathbb{W}^{3}_{h}(\hat{K}),$ (45)

where $g_{K}$ is the map between reference cell $\hat{K}$ and mesh cell $K$,
with derivative $J$, $J^{-T}$ means the inverse of the transpose of $J$,and
$\circ$ indicates function composition _i.e._ $f\circ g$ is another function
with $(f\circ g)(x)=f(g(x))$. This ensures that the discrete de Rham complex
property is preserved under the mapping from $\hat{K}$ to $K$. To define
$\mathbb{W}_{\theta}(\Omega)$, we just use straightforward composition with
the reference to cell map as is done for $\mathbb{W}_{h}^{0}$. This leads to
some differences between $\mathbb{W}_{h}^{\theta}$ and the vertical compoent
of $\mathbb{W}_{h}^{2}(\Omega)$ when terrain following coordinates are used.

In two dimensions when using (18), this becomes

$\displaystyle\psi\in\mathbb{V}_{h}^{0}(K)\\!:\,$ $\displaystyle\psi\circ
g_{K}=\hat{\psi},\,\hat{\psi}\in\mathbb{V}^{0}_{h}(\hat{K}),$ (46)
$\displaystyle w\in\mathbb{V}_{h}^{1}(K)\\!:\,$ $\displaystyle u\circ
g_{K}=J\hat{w}/\det J,\,\hat{w}\in\mathbb{V}^{1}_{h}(\hat{K}),$ (47)
$\displaystyle\phi\in\mathbb{V}_{h}^{2}(K)\\!:\,$ $\displaystyle\phi\circ
g_{K}=\hat{\phi}/\det J,\,\hat{\phi}\in\mathbb{V}^{2}_{h}(\hat{K}).$ (48)

In geophysical fluid dynamics applications, being able to solve equations on
the surface of a sphere is important. In general, two dimensional complexes
can be extended to orientable manifolds embedded in three dimensions by
restricting vector fields to be tangential to the manifold at each point.
Then, for such a vector field $u$, we define $u^{\perp}=k\times u$,
$\nabla^{\perp}u=k\times\nabla u$ and $\nabla^{\perp}\cdot
u=k\cdot\nabla\times u$, where $\nabla$ is now the projection of the gradient
onto the tangent plane. In fact, these operations can all be given intrinsic
definitions on any two dimensional manifold without reference to an external
space $\mathbb{R}^{3}$ containing the manifold, best expressed using the
language of differential forms. However, we do not do this here (see Arnold et
al. (2006) for intrinsic constructions using differential forms).

The geometric factors in these formulae introduce complications. These are
related both to the approximation properties of the spaces and to their
computer implementation, the latter due to the resulting nonpolynomial
integrands in weak formulations. When $g_{K}$ is an affine transformation, $J$
is constant on each cell, and no alterations to the approximation properties
arise. However, when $g_{K}$ is nonaffine, $J$ is nonconstant. This means that
the transformed basis functions may not span the same polynomial spaces as
they do on the reference cell. (Arnold et al. 2005, Falk et al. 2011) showed
that this occurs for transformed H(div) and H(curl) elements on nonaffine
quadrilaterals and hexahedra. This interferes with the standard approximation
theory error estimates, because they apply the Bramble-Hilbert lemma
considering the largest polynomial space spanned by the basis. The degradation
of approximation theory was demonstrated in practice in those papers. This
presents a concern for the applicability of these spaces in geophysical fluid
dynamics, because we encounter nonaffine transformations when using
quadrilaterals or higher order triangular cells (i.e. triangles that have been
curved to better approximate the sphere) to approximate the surface of the
sphere. These spaces are also used when extruding the sphere radially to make
a spherical annulus, required for three dimensional atmosphere and ocean
models; this leads to nonaffine prismatic meshes. Nonaffine cells also arise
when terrain following meshes are used. These are meshes that slope layers up
and down to conform to mountain ranges at the Earth surface in the atmosphere
and ocean. Fortunately, these approximation issues can be avoided through the
framework of Holst and Stern (2012), who used Strang type estimates to
consider the “variational crime” of a sequence of meshes that only conform to
a manifold in the limit. Provided that the meshes can be obtained by piecewise
smooth mappings from an affine mesh (i.e. a mesh consisting of cells mapped to
the reference cell by affine transformations), then approximation error bounds
can be obtained that match those of the reference cell. For meshes of interest
in geophysical fluid dynamics, this covers the sphere meshes described above,
as well as terrain following meshes obtained by smooth transformation
(smoothing of topography is standard practice in atmosphere and ocean
modelling). For the case of meshes (nonaffine or otherwise) of the sphere
extruded into a spherical annulus, Natale et al. (2016) showed that these
meshes can be obtained via transformation from an affine mesh embedded in four
dimensions.

Additionally, on unstructured meshes one must take care that the degrees of
freedom are correctly matched up on facets, and with the correct sign (since
they involve normal and tangential components in general). A systematic
approach for this was set out by Rognes et al. (2010), which is now
implemented in a number of finite element software systems, such as FEnICS
(Logg et al. 2012) and Firedrake (Rathgeber et al. 2016). Such systems are
very useful as an aid to productivity when dealing with the complexities of
compatible finite elements.

###  2.7 Replacing $\nabla\cdot$ with $\operatorname{DIV}$

One practical modification to this framework on nonaffine meshes is to replace
the transformation for the $L^{2}$ space ($V=\mathbb{V}^{2}_{h}$ in 2D or
$V=\mathbb{W}^{3}_{h}$ in 3D) with straightforward composition, i.e.

$\displaystyle\phi\in{V}(K)\\!:\,$ $\displaystyle\phi\circ
g_{K}=\hat{\phi},\,\hat{\phi}\in{V}(\hat{K}).$ (49)

Then, the discrete de Rham complex can be restored by replacing $\nabla\cdot$
with $\operatorname{DIV}:=P_{V}\nabla\cdot$, where $P_{V}$ is the $L^{2}$
projection into $V$. This is always a local operation since $V$ is a
discontinuous space with no interelement coupling. The diagram still commutes
upon replacing $d^{N}$ (where $N$ is the dimension) with $\operatorname{DIV}$,
since if $u\in H(\operatorname{div};\Omega)$, then
$\pi_{N}\circ\nabla\cdot=\operatorname{DIV}$. This idea was originally
proposed in Bochev and Ridzal (2009), who defined $\operatorname{DIV}$
equivalently in the case of $RT_{0}$ spaces on quadrilateral grids using the
discrete exterior calculus framework by evaluating fluxes through edges, after
which the definition as projection of the divergence can be obtained using the
divergence theorem. The extension to more general H(div) finite element spaces
on nonaffine meshes was probably clear to those authors, but it was also
discussed from a practical perspective in Shipton et al. (2018).

###  2.8 Primal dual grids

There have been various attempts to build complementary spaces on dual grids,
in order to avoid the global mass solves that are required to compute
$\delta$. These spaces are linked by “discrete Hodge star” operators
$(\star_{0},\star_{1},\star_{2})$, with mappings indicated in the following
diagram,

$\begin{CD}{\mathbb{V}}^{0}_{h}@>{\operatorname{d}^{1}=\nabla^{\perp}}>{}>{\mathbb{V}}^{1}_{h}@>{\operatorname{d}^{2}=\nabla\cdot}>{}>\mathbb{V}^{2}_{h},\\\
@V{}V{\star_{0}}V@V{}V{\star_{1}}V@V{}V{\star_{2}}V\\\
\tilde{\mathbb{V}}^{2}_{h}<EMAIL_ADDRESS>(50)

It is important that the discrete Hodge star maps are invertible. The approach
is clearest for lowest order spaces, where there is at most one nodal degree
of freedom per edge, vertex or cell in the definition of each space, in
analogy with the discrete exterior calculus (Hirani 2003). Thuburn and Cotter
(2015), Melvin and Thuburn (2017) constructed complementary spaces on
triangulations and their duals by subdividing cells into triangles and placing
$RT_{0}$ elements in each subcell. A constraint is applied so that the $d_{2}$
operators produce constant functions over the primal and dual cells
respectively. This scheme produces an extension of the primal dual finite
difference C grid approach to the consistent finite element setting
(inconsistencies in the Coriolis term occur on dual icosahedral and cubed
sphere setups in the framework of Thuburn and Cotter (2012), Thuburn et al.
(2014a)). The framework of mimetic spectral elements extends this idea to
higher order polynomial spaces (Lee et al. 2018).

##  3 Wave propagation properties

In this section we review the properties of compatible finite element
discretisations applied to linearised geophysical fluid dynamics. We will see
that the favourable properties of the C grid finite difference approach to
finite element methods. It is these properties that underpin the Met Office’s
choice to use compatible finite element methods to build their “Gung Ho”
atmospheric dynamical core that lies at the centre of their next generation
LFRic modelling system, for example (Adams et al. 2019).

###  3.1 Compatible discretisation of the linear rotating shallow water
equations

For now, we consider compatible finite element methods applied to the
linearised rotating shallow water equations on the doubly periodic plane with
constant Coriolis parameter,

$u_{t}+\underbrace{fu^{\perp}}_{\mbox{Coriolis}}+\underbrace{g\nabla\eta}_{\mbox{pressure
gradient}}=0,\,\eta_{t}+H\nabla\cdot u=0,$ (51)

where $u$ is the horizontal velocity, $h=\eta+H$ is the layer height with $H$
constant and $\eta$ integrating to zero, $f$ is the (constant) Coriolis
parameter, and $g$ is the acceleration due to gravity. The $\perp$ operator
rotates vectors one quarter of a rotation, to obtain
$u^{\perp}=(-u_{2},u_{1})$. These equations are solved in two dimensions.

We use the Helmholtz decomposition $u=\nabla^{\perp}\psi+\bar{u}+\nabla\phi$,
where $\bar{u}$ is a spatially constant (but possibly time dependent) vector
field, as the harmonic vector fields in the doubly periodic plane are of this
form. Then, we may write

$\displaystyle\phi_{t}-f\psi+g\eta$ $\displaystyle=0,$ (52)
$\displaystyle\psi_{t}+f\phi$ $\displaystyle=0,$ (53)
$\displaystyle\bar{u}_{t}+f\bar{u}^{\perp}$ $\displaystyle=0,$ (54)
$\displaystyle\eta_{t}+H\nabla^{2}\phi$ $\displaystyle=0.$ (55)

First, we observe that the harmonic component $\bar{u}$ is decoupled and
rotates at frequency $f$. These are called inertial oscillations. Second, we
can find steady solutions with $\phi=0$, $\bar{u}=0$, $\psi=g\eta/f$. These
solutions correspond to states of “geostrophic balance”, where the velocity is
divergence free and the pressure gradient term cancels out the Coriolis term.
The remaining solutions are inertia-gravity waves with $\phi$ nonzero.
Applying time derivatives to (55) and (52), we eliminate $\psi$ and $\eta$ (by
also using (53), to obtain

$\phi_{tt}+f^{2}\phi-gH\nabla^{2}\phi=0,$ (56)

which is the inertia-gravity wave equation. When $f=0$, this becomes the wave
equation with wavespeed $\sqrt{gH}$. When $f\neq 0$ (a positive sign would be
used for Northern hemisphere dynamics) then the equation becomes a Klein-
Gordon equation, which is dispersive.

Numerically induced oscillations (physical and spurious) in shallow water
models have been extensively examined by Le Roux in a series of papers using
dispersion analysis (Le Roux et al. 2007, Roux and Pouliot 2008, Rostand and
Le Roux 2008, Le Roux 2012), including for discretisations in the compatible
finite element family. Much of our vocabulary we use in this section has been
taken from that work.

To discretise this equation in space using compatible finite elements, we pick
a two dimensional discrete de Rham complex
$(\mathbb{V}_{h}^{0},\mathbb{V}_{h}^{1},\mathbb{V}_{h}^{2})$, and choose
$u\in\mathbb{V}_{h}^{1}$, $\eta\in\mathbb{V}_{h}^{2}$. Then the usual
introduction of inner products with test functions and integration by parts
leads to the discrete formulation,

$\displaystyle\left\langle w,u_{t}\right\rangle+\left\langle
w,fu^{\perp}\right\rangle-\left\langle\nabla\cdot w,g\eta\right\rangle$
$\displaystyle=0,\quad\forall w\in\mathbb{V}_{h}^{1},$ (57)
$\displaystyle\left\langle\alpha,\eta_{t}+H\nabla\cdot u\right\rangle$
$\displaystyle=0,\quad\forall\alpha\in\mathbb{V}_{h}^{2}.$ (58)

Since $\eta_{t}+H\nabla\cdot u\in\mathbb{V}_{h}^{2}$, we may choose
$\alpha=\eta_{t}+H\nabla\cdot u$ in (58), concluding that

$\eta_{t}+H\nabla\cdot u=0\mbox{ in }L^{2}(\Omega).$ (59)

This is a useful property that we will use later.

###  3.2 Geostrophic balance

In large scale atmosphere and ocean applications, it is very important that
discretisations can preserve states of geostrophic balance well; this is
usually tested on the linear rotating shallow water equations as discussed
here. Cotter and Shipton (2012) proved that compatible finite element
discretisations will have exactly steady geostrophic balanced solutions. To be
precise, for any divergence-free velocity field $u$ with zero harmonic
component $\bar{u}$, there exists an $\eta$ such that $(u,\eta)$ form a steady
solution in geostrophic balance. To show this, we just find
$\psi\in\mathbb{V}_{h}^{0}$ such that $u=\nabla^{\perp}\psi$ (possible from
the discrete Helmholtz decomposition), and pick $\eta$ as the $L^{2}$
projection of $\psi$ into $\mathbb{V}_{2}$ before multiplying by $g/f$. Then,

$\displaystyle\left\langle\nabla\cdot
w,g\eta\right\rangle=\left\langle\nabla\cdot
w,f\psi\right\rangle=-\left\langle w,f\nabla\psi\right\rangle=\left\langle
w,fu^{\perp}\right\rangle,\quad\forall w\in\mathbb{V}_{h}^{1},$ (60)

where the first equality holds from the projection, since $\nabla\cdot
w\in\mathbb{V}_{h}^{2}$. The second equality holds by integration by parts,
which is exact because $w\in H(\operatorname{div})$ and $\psi\in H^{1}$. The
final equality follows from $u=\nabla^{\perp}\psi$.

What is not true is that for every $\eta$ there exists a $\psi$ giving a
steady state solution. This is because the $L^{2}$ projection from
$\mathbb{V}_{h}^{0}$ to $\mathbb{V}_{h}^{2}$ is not a bijection. However, this
does not hold in the C grid finite difference case either.

It is also important that discretisations correctly represent inertial
oscillations. In the linear rotating shallow water equations, the only
solutions with $\eta=0$ are the inertial oscillations with spatially constant
$u=\bar{u}$ rotating at frequency $f$. It is important that discretisations
are free of spurious additional inertial modes, _i.e._ solutions with $\eta=0$
but with spatially varying $u$. Le Roux (2012) examined spurious modes in
various finite element discretisations, showing that when they are present
they lead to degraded error convergence rates. They have also been observed to
lead to problems in practical ocean model simulations, where they can be
excited by nonlinearity in baroclinic jets. This results in the formation of
spurious gridscale oscillatory patterns that do not change the pressure/layer
depth (Danilov and Kutsenko 2019). In a closed bounded domain with boundary
condition $u\cdot n=0$, we do not expect inertial oscillations because the
space of harmonic vector fields only contains 0. However, it is possible for
spurious inertial oscillations to satisfy the boundary condition, leading to
their excitation.

###  3.3 Inertial oscillations

Natale et al. (2016) showed that compatible finite element discretisations
applied to the rotating shallow water equations in the periodic plane have the
following property: the only time-varying solutions of (57-58) with $\eta=0$
have spatially constant $u_{t}$, corresponding to inertial oscillations
oscillating with frequency $f$. Any time independent solutions are in the
kernel of the discrete Coriolis operator, i.e. $u\in\mathbb{V}_{h}^{1}$ such
that

$\left\langle w,u^{\perp}\right\rangle=0,\quad\forall w\in\mathbb{V}_{h}^{1}.$
(61)

In other words, these discretisations are free from inertial oscillations. To
see this, first note that if $\eta=0$, then (59) implies that $\nabla\cdot
u=0$, so that $u=k+\nabla^{\perp}\psi$ for $k\in\mathfrak{h}_{h}$ and
$\psi\in\mathbb{V}_{h}^{0}$, from the discrete Helmholtz decomposition. Using
$w=\nabla^{\perp}\gamma$ in (57), we get

$\left\langle\nabla^{\perp}\gamma,u_{t}\right\rangle=-f\left\langle\nabla^{\perp}\gamma,u^{\perp}\right\rangle=f\left\langle\gamma,\nabla\cdot
u\right\rangle=0,\quad\forall\gamma\in\mathbb{V}_{h}^{0}.$ (62)

Hence, $u_{t}$ is orthogonal to $B_{h}^{k}$, _i.e._ ,
$u_{t}\in\mathfrak{h}_{h}$. We know from the discrete Hodge-Helmholtz
decomposition that $\dim(\mathfrak{h}_{h})=\dim(\mathfrak{h})=2$ on the
periodic plane; $\mathfrak{h}$ are the constant vector fields. In fact, since
the constant vector fields are in $\mathbb{V}_{h}^{1}$, and they are
divergence free, and in the kernel of $\delta_{h}$. Time independent
divergence free solutions satisfy

$0=\left\langle w,u_{t}\right\rangle=-f\left\langle
w,u^{\perp}\right\rangle,\quad\forall w\in\mathbb{V}_{h}^{1},$ (63)

_i.e._ they are in the kernel of the discrete Coriolis operator. Vector fields
in this kernel are referred to as “Coriolis modes”. They are the main downside
of compatible discretisations, but the dimension of this kernel is always
found to be finite and resolution independent in analyses by Le Roux Rostand
and Le Roux (2008); in fact this number is typically very small. These
Coriolis modes are also always found in C grid finite difference
discretisations.

###  3.4 Inertia gravity waves

Cotter and Shipton (2012) also examined the discrete inertia-gravity waves
that correspond to solutions of the Klein Gordon equation above. Writing
$u=\nabla^{\perp}\psi+\delta\phi$ for $\psi\in\mathbb{V}_{0}$ and
$\phi\in\mathbb{V}_{2}$ (having already discarded the harmonic component since
it decouples), and choosing both $w=\nabla^{\perp}\gamma$ for
$\gamma\in\mathbb{V}_{h}^{0}$ and $w=\delta\alpha$ for
$\alpha\in\mathbb{V}_{h}^{2}$, we obtain

$\displaystyle\left\langle\delta\alpha,\delta\phi_{t}\right\rangle-f\left\langle\delta\alpha,\nabla\psi\right\rangle-g\left\langle\delta\alpha,\delta\eta\right\rangle$
$\displaystyle=0,\quad\forall\gamma\in\mathbb{V}_{2},$ (64)
$\displaystyle\left\langle\nabla\gamma,\nabla\psi_{t}\right\rangle+f\left\langle\nabla\gamma,\delta\phi\right\rangle$
$\displaystyle=0,\quad\forall\gamma\in\mathbb{V}_{h}^{0},$ (65)
$\displaystyle\left\langle\phi,\eta_{t}+H\nabla\cdot\delta\phi\right\rangle$
$\displaystyle=0,\quad\forall\phi\in\mathbb{V}_{h}^{2}.$ (66)

Using the definition of $\delta$, and the fact that the height equation holds
in $L^{2}$, we get

$\displaystyle-\left\langle\alpha,\nabla\cdot\delta\phi_{t}\right\rangle+f\left\langle\alpha,\nabla^{2}\psi\right\rangle+g\left\langle\alpha,\nabla\cdot\delta\eta\right\rangle$
$\displaystyle=0,\quad\forall\gamma\in\mathbb{V}_{2},$ (67)
$\displaystyle\left\langle\nabla\gamma,\nabla\psi_{t}\right\rangle-f\left\langle\gamma,\nabla\cdot\delta\phi\right\rangle$
$\displaystyle=0,\quad\forall\gamma\in\mathbb{V}_{h}^{0},$ (68)
$\displaystyle\eta_{t}+H\nabla\cdot\delta\phi=0.$ (69)

We can recognise the operator $\nabla\cdot\delta$ as the mixed approximation
$\tilde{\nabla}^{2}$ of the Laplacian defined by

$\displaystyle\left\langle\alpha,\tilde{\nabla}^{2}\phi\right\rangle-\left\langle\alpha,\nabla\cdot\sigma\right\rangle$
$\displaystyle=0,\quad\forall\alpha\in\mathbb{V}_{h}^{2}.$ (70)

After restricting to $\overline{\mathbb{V}}_{h}^{2}$, defined by

$\overline{\mathbb{V}}_{h}^{2}=\left\\{\phi\in\mathbb{V}_{h}^{2}:\int_{\Omega}\phi\operatorname{d}x=0\right\\},$
(71)

this discretisation is well known to be invertible, stable and convergent.
Hence, when $f=0$, we can deduce

$\phi_{tt}-gH\tilde{\nabla}^{2}\phi=0,$ (72)

_i.e._ the mixed approximation of the wave equation for
$\phi\in\overline{\mathbb{V}}_{h}^{2}$.

When $f\neq 0$, we have to introduce projection operators
$P_{0}:\overline{\mathbb{V}}_{h}^{2}\to\overline{\mathbb{V}}_{h}^{0}$,
$P_{2}:\overline{\mathbb{V}}_{h}^{0}\to\overline{\mathbb{V}}_{h}^{2}$, defined
by

$\displaystyle\left\langle\nabla\gamma,\nabla P_{0}\phi\right\rangle$
$\displaystyle=\left\langle\nabla\gamma,\delta\phi\right\rangle,\quad\forall\gamma\in\overline{\mathbb{V}}_{h}^{0},$
(73)
$\displaystyle\left\langle\alpha,-\tilde{\nabla}^{2}P_{2}\psi\right\rangle:=\left\langle\delta\alpha,\delta
P_{2}\psi\right\rangle$
$\displaystyle=\left\langle\delta\alpha,\nabla\phi\right\rangle,\quad\forall\gamma\in\overline{\mathbb{V}}_{h}^{2}.$
(74)

$P_{0}$ is well posed since it just involves solving the usual Galerkin
discretisation of the Laplacian on $\overline{\mathbb{V}}_{h}^{0}$, whilst
$P_{2}$ is well posed since it involves solving the mixed discretisation as we
have already discussed.

Using $P_{0}$ and $P_{2}$, we get

$\displaystyle\phi_{t}-fP_{2}\psi+gh$ $\displaystyle=0,$ (75)
$\displaystyle\psi_{t}+fP_{0}\phi$ $\displaystyle=0,$ (76)
$\displaystyle\eta_{t}+H\tilde{\nabla}^{2}\phi$ $\displaystyle=0,$ (77)

from which we can deduce the discrete Klein-Gordon equation

$\phi_{tt}+f^{2}P_{2}P_{0}\phi-gH\tilde{\nabla}^{2}\phi=0.$ (78)

The behaviour of the numerical dispersion relation depends on the kernel of
the composition $P_{2}P_{0}$, for which a lower bound is obtained by
considering the relative sizes of $\bar{\mathbb{V}}_{h}^{0}$ and
$\bar{\mathbb{V}}_{h}^{2}$. For $RT_{k}$ on quadrilaterals, these two spaces
have the same dimension in the periodic domain, so this projection is not too
harmful. For $RT_{k}$ spaces on triangles,
$\dim(\bar{\mathbb{V}}_{2})>\dim(\bar{\mathbb{V}}_{0})$, so $P_{2}P_{0}$ is
not surjective. In the case of the C grid finite difference method on
triangles, a similar issue arises, causing high and low frequency branches of
the numerical dispersion relation to intertangle, leading to numerical noise
when $f$ is sufficiently large (Danilov 2010). Further analysis is required to
really pin down these issues in the compatible finite element case. For
$BDM_{k}$ spaces on triangles, we have the opposite situation,
$\dim(\bar{\mathbb{V}}_{2})<\dim(\bar{\mathbb{V}}_{0})$, so there is at least
the chance for $P_{2}P_{0}$ to be surjective (although both projection
operators will have checkerboard modes in their kernel on structured meshes).
This suggests that $BDM$ spaces are more appropriate for geophysical fluid
dynamics using triangular meshes, but further analysis of these issues is
needed to make these statements more precise.

###  3.5 Spectral gaps and zero group velocity

Staniforth et al. (2013), Melvin et al. (2014) examined the numerical
dispersion relation for the cases of $RT_{0}$ and $RT_{1}$ on quadrilaterals,
motivated by building discretisations on the sphere using a cubed sphere grid.
This type of dispersion analysis allows to focus on group
($\partial\omega/\partial k$, where $\omega$ is the frequency and $k$ is the
wavenumber) and phase velocity ($\omega k/|k|^{2}$) for numerical
discretisations. For $RT_{0}$, the numerical dispersion relation is very
similar to the C grid finite difference numerical dispersion relation on
quadrilaterals, with no turning points for the group velocity except at
maximum wavenumbers. For $RT_{1}$, the numerical dispersion relation has two
roots for each wave number, corresponding to the resolution of higher
wavenumbers in the gridcell. When these are properly interpreted, the group
velocity is mostly well behaved except for a jump in the dispersion relation
at $\Delta x$ wavelengths where the group velocity goes to zero before and
after the jump. Remarkably, the jump is repairable by modifying the
coefficients in the mass matrix in such a way that the convergence rate is not
eroded. It is not really the jump itself which is the problem, but the repair
also makes the group velocity become nonzero through a L’Hôpital’s rule type
cancellation. This fix is independent of the value of $f$ and $H$, so is
useable in practice. In numerical experiments, Melvin et al. (2014) showed
that this modification leads to propagation of a wave packet with $\Delta x$
wavelength which otherwise stays in the same location, spuriously.

Eldred et al. (2019) introduced an alternative approach to avoiding these
dispersion relation spectral gaps, in which the nodal variables are the same
as $RT_{0}$ and $DG_{0}$ on quadrilaterals, but a higher order polynomial
expansion is constructed by using nodal variables from surrounding cells. This
can be seen as a form of spline, but does not increase the degree of
continuity, just the polynomial degree. The effect on the numerical dispersion
relation is that there is only one branch, so there cannot be jumps. The
downside of this approach is that the stencil of the operators is extended to
more cells, and the standard approach of finite element assembly becomes more
complicated. Eldred and Le Roux (2018, 2019) examined the spectral gaps in
$RT_{k}$ elements for larger k, and showed that this spline approach also
fixes the problem for higher k.

###  3.6 Rossby waves

Following the approach of Thuburn (2008), Cotter and Shipton (2012) also
examined the Rossby wave propagation properties of compatible finite element
spaces. Rossby waves occur in the situation where the Coriolis parameter $f$
is spatially varying (as is the case on the sphere under the “traditional
approximation”, where $f=2\Omega\cdot n$, with $\Omega$ being the rotational
velocity of the sphere, and $n$ being the normal to the sphere surface). To
perform the analysis, we consider solutions on an infinite plane, with
$f=f_{0}+\beta y$, with $f_{0}$ and $\beta$ being constants. Then, if the
Rossby number $\operatorname{Ro}=U/fL$ is small (where $U$ is a typical
velocity scale and $L$ is a typical spatial scale), and also $\beta
L/f_{0}=\mathcal{O}(\operatorname{Ro})$, then we neglect $u_{t}$ and $\beta
yu^{\perp}$ in the velocity equation. This gives the geostrophic balanced
states $(u,h)=(u_{g},\eta_{g})$ satisfying

$f_{0}u^{\perp}_{g}=-g\nabla\eta_{g},$ (79)

so that $u_{t}=0$, and $\eta_{t}=-H\nabla\cdot
u=H\nabla\cdot(f_{0}\nabla^{\perp}\eta)=0$, as we saw previously. This is the
main reason why it is important for the discretisation to represent these
geostrophic balanced states. To obtain dynamics, we consider
$\mathcal{O}(\operatorname{Ro})$ corrections to $(u,\eta)$, which we write as
$(u_{ag},\eta_{ag})$ (ageostrophic velocity and height). The equations at the
next order in $\mathcal{O}(\operatorname{Ro})$ then give

$\displaystyle(u_{g})_{t}+f_{0}u_{ag}^{\perp}+\beta
yu_{g}^{\perp}+g\nabla\eta_{ag}=0,(\eta_{g})_{t}+H\nabla\cdot u_{ag}=0.$ (80)

Applying $-\nabla^{\perp}\cdot$ to the first equation and using the second
gives

$-\nabla^{2}\psi_{t}+\frac{f_{0}}{H}(\eta_{g})_{t}-\beta u_{g}\cdot\hat{y}=0,$
(81)

where we used that $\nabla^{\perp}u_{g}=0$, so $u_{g}=\nabla^{\perp}\psi$
(having already eliminating inertial oscillations which are fast), and
$\hat{y}$ is the unit vector in the $y$-direction. Finally this becomes

$(\frac{f_{0}^{2}}{gH}-\nabla^{2})\psi_{t}-\beta\frac{\partial\psi}{\partial
x}=0,$ (82)

which is the Rossby wave equation, which exhibits waves propagating Westwards
when $\beta>0$ (_i.e._ , in the Northern Hemisphere).

Now we examine what happens with compatible finite element discretisations in
the low $\operatorname{Ro}$ limit. At leading order in Rossby number, we
obtain the geostrophic balance equation,

$\displaystyle\left\langle
w,fu^{\perp}_{g}\right\rangle-g\left\langle\nabla\cdot
w,\eta_{g}\right\rangle$ $\displaystyle=0,\quad\forall
w\in\mathbb{V}_{h}^{1},$ (83) $\displaystyle\left\langle\alpha,H\nabla\cdot
u_{g}\right\rangle$ $\displaystyle=0,\quad\forall\mathbb{V}_{h}^{2}.$ (84)

We have already seen that this has solutions $u_{g}=\nabla^{\perp}\psi$,
$h=P_{2}(f\psi)/h$ for $\psi\in\mathbb{V}_{h}^{0}$. At the next order we have

$\displaystyle\left\langle w,(u_{g})_{t}\right\rangle+\left\langle
w,f_{0}u_{ag}^{\perp}+\beta yu_{g}\right\rangle-\left\langle\nabla\cdot
w,g\eta_{ag}\right\rangle$ $\displaystyle=0,\quad\forall
w\in\mathbb{V}_{h}^{2},$ (85)
$\displaystyle\left\langle\alpha,(\eta_{g})_{t}+H\nabla\cdot(u_{g})\right\rangle$
$\displaystyle=0,\quad\forall\alpha\in\mathbb{V}_{h}^{1}.$ (86)

Choosing $w=\nabla^{\perp}\gamma$ with $\gamma\in\mathbb{V}_{h}^{0}$, and
integrating by parts in (85) (permissible since $\gamma\in H^{1}$ and
$u_{ag}\in H(\operatorname{div})$) gives

$\left\langle\nabla\gamma,\nabla\psi\right\rangle-\left\langle\gamma,f\nabla\cdot
u_{ag}\right\rangle-\left\langle\gamma,\beta\underbrace{u_{g}\cdot\hat{y}}_{=\frac{\partial\psi}{\partial
x}}\right\rangle=0,\quad\forall\gamma\in\mathbb{V}_{h}^{0}.$ (87)

(86) implies that

$f_{0}P_{2}(\psi)/g+H\nabla\cdot u_{g}=0,\mbox{ in }L^{2}.$ (88)

Combining this with (87) then gives the discrete Rossby wave equation,

$\left\langle\nabla\gamma,\nabla\psi_{t}\right\rangle+\left\langle\gamma,\frac{f_{0}^{2}}{gH}P_{2}\psi\right\rangle-\left\langle\gamma,\beta\frac{\partial\psi}{\partial
x}\right\rangle=0,\quad\forall\gamma\in\mathbb{V}_{h}^{0}.$ (89)

Without the $P_{2}$ projection operator, this would just be a regular $H^{1}$
finite element approximation of the Rossby wave equation. With it, there is
the possibility of some projection errors altering the numerical dispersion
relation, especially in the case of $BDM$ elements where
$\dim(\mathbb{V}_{h}^{0})>\dim(\mathbb{V}_{h}^{2})$. However, this will only
occur for high wavenumber waves, where the dynamics is already dominated by
Laplacian term, and the Rossby waves will have very slow phase and group
velocities in either case. Hence, we do not believe that this causes a problem
for discrete Rossby wave propagation. A similar argument was made in Thuburn
(2008) when considering C grid discretisations on hexagons.

Rostand and Le Roux (2008) examined the wave propagation properties of the
$RT_{0}$-$DG_{0}$ and $BDM_{1}$-$DG_{0}$ compatible finite element
discretisations (and $RT_{0}$-$CG_{1}$ and $BDM_{1}$-$CG_{1}$ discretisations
that are not compatible) using discrete dispersion relations computed through
Fourier analysis. These dispersion relations revealed the steady geostrophic
modes for constant $f$, and two branches of the dispersion relation for
$BDM_{1}$-$DG_{0}$, a primary one attached to the origin (zero frequency for
zero wavenumber) and a secondary one, which they described as spurious. It is
possible that this second branch can be interpreted as corresponding to higher
wavenumbers resolved in the cell, just as for $RT_{1}$ on quadrilaterals as
explored by Staniforth et al. (2013), but doing these calculations is
difficult on triangles and more work is needed to clarify this. In their
dispersion analysis, Rostand and Le Roux (2008) identified “CD modes” in the
$BDM_{1}$-$DG_{0}$ discretisation, which are modes in the intersection of the
Coriolis operator and the divergence operator. These precisely correspond to
the modes in the kernel of the $P_{2}$ operator appearing in (89). In
experiments with the linear rotating shallow water equations with balanced
initial data and $f=f_{0}+\beta y$, designed to examine Rossby wave
propagation, they observed accurate solutions on structured grids but very
noisy solutions on unstructured grids. The noise was attributed to
interactions with the CD modes. However, when we have repeated these
experiments using modern automated finite element systems such as FEniCS and
Firedrake, we have not observed this noisy behaviour with $BDM_{1}$ on
unstructured grids. It seems likely that Rostand and Le Roux (2008) had bugs
in their implementation related to the identification of the two nodal
variables between two cells across each edge (since no such problem arose with
$RT_{0}$, which only has one nodal variable per edge). This is understandable,
because the problem of how to systematically assemble $BDM_{1}$ and higher
order H(div) spaces on triangles was not solved until Rognes et al. (2010);
these things are very difficult to implement by hand. These results may have
discouraged the adoption of the $BDM$ family for atmosphere and ocean
modelling, but it seems like a good option for lowest order spaces to avoid
the issues with spurious inertia-gravity wave propagation with $RT_{0}$.

###  3.7 Consistent linear tidal response

Cotter and Kirby (2016) considered solutions of the linearised barotropic tide
equations, which take the form,

$u_{t}+fu^{\perp}+g\nabla(D+b)=-cu+F,\qquad D_{t}+H\nabla\cdot u=0.$ (90)

These are the linear rotating shallow water equations with additional
topography $b$, spatiotemporal lunar forcing $F$ and linear friction
coefficient $c$. They used the Helmholtz equation to obtain an exponentially
damping lower and upper bound on the energy in the absence of forcing. In the
presence of time dependent forcing (quasiperiodic forcing is appropriate for
tidal models), they proved that the solution converges at exponential rate to
a time dependent solution as $t\to\infty$, independent of the initial
condition. This is the solution that is of interest when predicting tides.
When compatible finite element methods are used to discretise the tide
equations, they proved that the continuous time discrete space solution also
converges exponentially to a time dependent numerical solution, independent of
the initial condition. Finally they showed that this discrete attracting
solution converges to the unapproximated attracting solution as the mesh is
refined. Cotter et al. (2018) extended this analysis to a nonlinear model with
$cu$ replaced by $c|u|u$, which is the more realistic damping model that is
actually used by oceanographers. This nonlinear case is surprisingly subtle
but they were able to prove long time stability of the system and obtain rates
of damping in the unforced case. These were used to prove error estimates for
the discrete solution obtained using compatible finite element methods. Kirby
and Kernell (2021) used the compatible finite element framework to design a
preconditioner for the implicit solver for the tidal equations, proving that
the convergence rates are independent of mesh resolution. Cotter et al. (2022)
extended this approach to the multiple layer version of this model.

###  3.8 Hydrostatic balance

Finally in this section, we discuss the discrete hydrostatic balance
properties of compatible finite element methods. Later, we shall introduce
three dimensional geophysical fluid dynamics models with gravity and pressure
gradient terms, so that the velocity equation takes the form

$\frac{\partial u}{\partial t}+\ldots+\underbrace{\nabla p}_{\mbox{pressure
gradient}}=\underbrace{b\hat{k}}_{\mbox{gravity}},$ (91)

in the case of the Boussinesq equations (typically used in ocean modelling),
where $p$ is the pressure, $b$ is the buoyancy, and $\hat{k}$ is the unit
normal vector in the “up” direction. In the case of the compressible Euler
equations (typically using in atmosphere modelling), we have

$\frac{\partial u}{\partial
t}+\ldots+c_{p}\underbrace{\theta\nabla\Pi}_{\mbox{pressure
gradient}}=\underbrace{-g\hat{k}}_{\mbox{gravity term}},$ (92)

where $\theta$ is the potential temperature, $\Pi$ is the Exner pressure, and
$g$ is the acceleration due to gravity. In both cases, we are concerned with
states of hydrostatic balance, which is when the vertical component of the
pressure gradient term balances the gravity term. Either we are considering
hydrostatic models, where this balance is enforced exactly in the model, or we
are considering nonhydrostatic models where it is important that these
hydrostatic states can be accurately represented. In the compatible finite
element case, to study the vertical part of the velocity equation, we restrict
the test function in the velocity equation to $\mathbb{W}_{h}^{2,V}$ (assuming
a tensor product discrete de Rham complex), the vertical part of the space
$\mathbb{W}_{h}^{2}$ containing the discretised velocity $u$.

In the case of the Boussinesq equations, the discrete hydrostatic balance is
written (after integration by parts) as

$-\left\langle\nabla\cdot w,p\right\rangle=\left\langle
w\cdot\hat{k},b\right\rangle,\quad\forall
w\in\mathring{\mathbb{W}}_{h}^{2,V},$ (93)

for pressure $p\in\mathbb{W}_{h}^{3}$ and buoyancy
$b\in\mathbb{W}_{h}^{\theta}$, where $\mathring{\mathbb{W}}_{h}^{2,V}$ is the
subspace of the vertical space $\mathbb{W}_{h}^{2,V}$ requiring the boundary
condition $u\cdot n=0$ at the top and the bottom. Despite appearances, this is
actually only defining the vertical part of the pressure gradient term, since
$w\in\mathring{\mathbb{W}}_{h}^{2,V}$ always points in the vertical direction.

Given $p$, there is a unique $b$ that satisfies this hydrostatic balance. To
see this, we note that if the layers of the mesh are flat, then $w\cdot
k\in\mathbb{W}_{h}^{\theta}$ for all $w\in\mathbb{W}_{h}^{2,V}$. If the layers
are not flat, _i.e._ for terrain following coordinates, then there exists
$0\leq\kappa\leq\infty$ such that $\kappa w\cdot k\in\mathbb{W}_{h}^{\theta}$
for all $w\in\mathbb{W}_{h}^{2,V}$. After replacing $w\cdot
k=\kappa^{-1}\gamma$ for $\gamma\in\mathbb{W}_{h}^{\theta}$, we recognise the
right hand side of (93) as a nondegenerate weighted $L^{2}$ inner product,
hence $b$ is unique.

To discuss the nature of the uniqueness of $p$, we consider an alternative
boundary condition with $u\cdot n=0$ on the bottom, and $p=p_{0}$ on the top
(for some chosen $p_{0}$ which may depend on the horizontal coordinate). The
equation after integration by parts and use of the top boundary condition
gives

$-\left\langle\nabla\cdot w,p\right\rangle+\left\llangle w\cdot
n,p_{0}\right\rrangle=\left\langle w\cdot\hat{k},b\right\rangle,\quad\forall
w\in\mathring{\mathbb{W}}_{h}^{2,V},$ (94)

where $\mathring{\mathbb{W}}_{h}^{2,V}$ is now the subspace with vanishing
normal component on the bottom only. To analyse this problem, Natale et al.
(2016) introduced the following formulation, defining
$(v,p)\in\mathbb{W}_{h}^{2,V}\times\mathbb{W}_{h}^{3}$ such that

$\displaystyle\left\langle w,v\right\rangle-\left\langle\nabla\cdot
w,p\right\rangle$ $\displaystyle=\left\langle
w\cdot\hat{k},b\right\rangle-\left\llangle w\cdot
n,p_{0}\right\rrangle,\quad\forall w\in\mathring{\mathbb{W}}_{h}^{2,V},$ (95)
$\displaystyle\left\langle\phi,\nabla\cdot v\right\rangle$
$\displaystyle=0,\qquad\forall\phi\in\mathbb{W}_{h}^{3},$ (96)

which we recognise as a mixed problem defined on
$\mathbb{W}_{h}^{2,V}\times\mathbb{W}_{h}^{3}$. Using the same arguments as we
have previously, at the solution we have $\nabla\cdot v=0$ in $L^{2}$. Since
$v\cdot n=0$ on the bottom surface, and $v$ points in the vertical direction,
we conclude that $v=0$, and therefore $p$ solves (94). Natale et al. (2016)
showed that this type of vertical mixed problem has a unique solution $(v,p)$.
Hence, there is a one to one correspondence between $p$ and $b$, as required.
If there are boundary conditions $u\cdot n=0$ on both top and bottom surfaces,
$p$ is only determined up to the value $p_{0}$ restricted to the upper
surface, also as required.

In the case of the compressible Euler equations, taking boundary conditions
$u\cdot n=0$ on the bottom surface and $\Pi=\Pi_{0}$ on the top surface, the
discrete hydrostatic balance is written (after integrating by parts) as

$-\left\langle\nabla\cdot(\theta w),\Pi\right\rangle=\left\langle
w\cdot\hat{k},g\right\rangle-\left\llangle w\cdot
n,\Pi_{0}\right\rrangle,\quad\forall w\in\mathring{\mathbb{W}}_{h}^{2,V},$
(97)

for $\Pi\in\mathbb{W}_{h}^{3}$ and $\theta\in\mathbb{W}_{h}^{\theta}$. Using
an extension of the techniques described for the Boussinesq equation, Natale
et al. (2016) proved similar results. This motivates the use of
$\mathbb{W}_{h}^{\theta}$ for temperature variables like $b$ or $\theta$.
Melvin et al. (2018) showed through linear dispersion analysis applied to the
compressible Boussinesq equations that this choice does indeed lead to an
absence of spurious hydrostatic modes that would appear if
$\theta\in\mathbb{W}_{h}^{3}$ or $\mathbb{W}_{h}^{0}$.

##  4 Transport and stabilisation

Hopefully it is clear from Section 3 that it might be interesting to consider
designing a numerical atmosphere or ocean model using compatible finite
element methods. One very important aspect of these models is the choice of
transport schemes, _i.e._ the discretisation of the advection operators. Since
different fields are restricted to different spaces from the discrete de Rham
complex (or $\mathbb{W}_{h}^{\theta}$) with different continuity constraints,
we need to consider a diverse range of transport schemes, some of which we
briefly survey in this section.

###  4.1 Transport of $H^{1}$ fields

For scalar fields in $\mathbb{W}_{0}$ or $\mathbb{V}_{0}$ (we will call it $V$
here), we consider the discretisation of the advection equation

$\frac{\partial q}{\partial t}+u\cdot\nabla q=0,$ (98)

for some specified $u\in\mathbb{W}_{h}^{2}$ or $\mathbb{V}_{h}^{1}$, where we
assume that $u\cdot n=0$ on exterior boundaries. Since $V$ is a continuous
finite element space, we simply take the $L^{2}$ inner product with a test
function and integrate by parts to obtain the standard continuous finite
element approximation, seeking $q\in V$ such that

$\left\langle\phi,\frac{\partial q}{\partial
t}\right\rangle-\left\langle\nabla\cdot(u\phi),q\right\rangle=0,\quad\forall\phi\in
V.$ (99)

As is well known, this discretisation tends to produce spurious oscillations
at regions of rapid changes in $q$. One way to suppress these oscillations is
to use the Streamline Upwind Petrov-Galerkin (SUPG) method (Brooks and Hughes
1982, Tezduyar et al. 1988, Tezduyar 1989). In this approach, the test
function $\phi$ is replaced by $\phi+\tau u\cdot\nabla\phi$, which biases it
in the upwind direction. $\tau$ is some chosen stabilisation parameter which
depends on the mesh, $u$ and other parameters. This leads to

$\left\langle\phi+\tau u\cdot\nabla\phi,\frac{\partial q}{\partial
t}\right\rangle-\left\langle\nabla\cdot(u\phi),q\right\rangle+{\left\langle\tau
u\cdot\nabla\phi,u\cdot\nabla q\right\rangle}=0,\quad\forall\phi\in V.$ (100)

The final term performs diffusion along streamlines of $u$, which tends to
reduce spurious oscillations. By applying this modification in the $\partial
q/\partial t$ term as well as the $u\cdot\nabla q$ term, we obtain a
consistent approximation (_i.e._ , substituting a smooth exact solution of the
unapproximated equation produces zero).

Another possibility is the edge stabilisation approach proposed for advection
equations and analysed (when combined with diffusion) in Burman and Hansbo
(2004), resulting in the formulation

$\left\langle\phi,\frac{\partial q}{\partial
t}\right\rangle-\left\langle\nabla\cdot(u\phi),q\right\rangle+\left\llangle\gamma
h^{2}\left[\\!\\!\left[\nabla\phi\right]\\!\\!\right],\left[\\!\\!\left[\nabla
q\right]\\!\\!\right]\right\rrangle_{\Gamma}=0,\quad\forall\phi\in V,$ (101)

where

$\left\llangle u,v\right\rrangle_{\Gamma}=\int_{\Gamma}u\cdot
v\operatorname{d}S,$ (102)

$\Gamma$ is the union of all interior facets $f$ in the mesh (_i.e._ , facets
joining two cells),
$\left[\\!\\!\left[v\right]\\!\\!\right]=v^{+}n^{+}+v^{-}n^{-}$ for vector
fields $v$, each facet $f$ has been arbitrarily assigned $+$ and $-$ labels to
its two sides, $a^{\pm}$ indicates the restriction of the discontinuous
function $a$ to the $\pm$ side of the facet, respectively, $h$ is a mesh edge
length parameter, and $\gamma$ is a (possibly $u$ and or $q$ dependent)
stabilisation parameter. This term has a diffusive effect across interior
facets, penalising jumps in $q$, without sacrificing consistency as $h\to 0$.

###  4.2 Transport of $L^{2}$ fields

The spaces $\mathbb{W}_{h}^{3}$ and $\mathbb{V}_{h}^{2}$ have no continuity
constraints, which allows for upwind stabilisation _via_ a discontinuous
Galerkin formulation. Here we consider the continuity equation,

$D_{t}+\nabla\cdot(uD)=0,$ (103)

where $u$ is as above. We introduce the discretisation by multiplying by a
test function and integrating over a single mesh cell $e$,

$\int_{e}\phi D_{t}-\nabla\phi\cdot uD\operatorname{d}x+\int_{\partial
e}\tilde{D}u\cdot n\phi\operatorname{d}S=0,\quad\forall\phi\in V(e),$ (104)

where $\partial e$ is the boundary of $e$ with outward pointing normal $n$,
$V$ is the chosen discontinuous space, and $\tilde{D}$ is the upwind value of
$D$, which must be defined in terms of the values of $D$ on the inside and the
outside of $e$. When $u\cdot n>0$, then $\tilde{D}$ is equal to $D$ from
inside $e$, otherwise the outside value is used. If we sum over all the cells
$e$ in the mesh, we obtain

$\left\langle\phi,D_{t}\right\rangle-\left\langle\nabla_{h}\phi,uD\right\rangle+\left\llangle\left[\\!\\!\left[\phi
u\right]\\!\\!\right],\tilde{D}\right\rrangle_{\Gamma}=0,\quad\forall\phi\in
V(\Omega),$ (105)

where $\nabla_{h}$ indicates the cellwise “broken” gradient,

$\nabla_{h}|_{e}D=\nabla|_{e}D,$ (106)

for each cell $e$ in the mesh. An analysis of how this choice of upwind
$\tilde{D}$ introduces stabilisation is provided in Brezzi et al. (2004). This
method is locally conservative.

###  4.3 Transport of $H(\operatorname{div})$ fields

In consideration of the $(u\cdot\nabla)u$ term in the velocity equation of
geophysical models, discretisations for vector advection equations of the form

$\frac{\partial v}{\partial t}+(u\cdot\nabla)v=0,$ (107)

for a vector field $v$ in $H(\operatorname{div})$ spaces
$\mathbb{W}_{h}^{2,V}$ or $\mathbb{V}_{h}^{1}$, and $u$ is again as above.
When the equation is solved on surface of the sphere, $v$ is constrained to be
tangential to the sphere (or the mesh approximating the sphere in the discrete
case). Then, the equation becomes

$\frac{\partial v}{\partial t}+\mathbb{P}_{S}\left((u\cdot\nabla)v\right)=0,$
(108)

where $\mathbb{P}_{S}$ is the Euclidean projection into the tangent plane to
the sphere.

Functions in $H(\operatorname{div})$ spaces are only partially continuous (in
the normal component across facets) so we need to start by considering an
upwind formulation on a single cell again,

$\int_{e}w\cdot v_{t}-\nabla\cdot(u\otimes w)\cdot
v\operatorname{d}x+\int_{\partial e}n\cdot
uw\cdot\tilde{v}\operatorname{d}S=0,\quad\forall w\in V(e),$ (109)

where $(a\otimes b)_{ij}=a_{i}b_{j}$ for vectors $a$ and $b$, $\tilde{v}$ is
the upwind value of $v$, $A:B=\sum_{ij}A_{ij}B_{ij}$ for two matrices $A$ and
$B$, and $V$ is the chosen $H(\operatorname{div})$ space. Summing up over all
of the cells in the mesh gives

$\int_{\Omega}w\cdot v_{t}-\nabla_{h}\cdot(u\otimes w)\cdot
v\operatorname{d}x+\int_{\partial\Gamma}\left[\\!\\!\left[u\otimes
w\right]\\!\\!\right]\cdot\tilde{v}\operatorname{d}S=0,\quad\forall w\in
V(\Omega),$ (110)

where $\left[\\!\\!\left[u\otimes w\right]\\!\\!\right]=(n^{+}\cdot
u^{+})w^{+}+(n^{-}\cdot u^{-})w^{-}$. Since $v\in V$ has continuous normal
components, $\tilde{v}$ only differs from $v$ in the tangential component.
Hence, the upwind stabilisation may be insufficient to adequately suppress
oscillations, depending on the shape of the mesh cells and the direction of
the velocity. On meshes approximating the sphere (and other manifolds), these
formulae require modification when $n^{+}\neq-n^{-}$ on an edge. The
modification rotates $\tilde{u}$ into the tangent plane of the cell $e$, as
described in Bernard et al. (2009). The projection of the advection equation
into the tangent to the sphere is naturally dealt with in (110), because
$w\in\mathbb{V}_{h}^{2}$ is always tangential to the surface mesh.

The vorticity form is an alternative form of the vector advection equation,
given by

$\frac{\partial v}{\partial t}+(\nabla\times v)\times
u+\frac{1}{2}\nabla(u\cdot v)+\frac{1}{2}\left((\nabla v)^{T}u-(\nabla
u)^{T}v\right)=0,$ (111)

in three dimensions, where $(\nabla v)^{T}_{ij}=\partial u_{j}/\partial
x_{i}$. In two dimensions this is written

$\frac{\partial v}{\partial t}+(\nabla^{\perp}\cdot v)\cdot
u^{\perp}+\frac{1}{2}\nabla(u\cdot v)+\frac{1}{2}\left((\nabla v)^{T}u-(\nabla
u)^{T}v\right)=0,$ (112)

where $w^{\perp}=(-w_{2},w_{1})$, $\omega=\nabla^{\perp}\cdot
w:=-\frac{\partial w_{2}}{\partial x_{1}}+\frac{\partial w_{1}}{\partial
x_{2}}$ for a vector field $w$ in planary geometry. On the sphere, with
outward pointing normal $\hat{k}=x/|x|$ (and $x$ is the three dimensional
coordinate with origin at the centre of the sphere), we have
$w^{\perp}=\hat{k}\times w$ and $\nabla^{\perp}\cdot
w=\hat{k}\cdot\nabla\times w$, where $\nabla$ is now the projection of the
gradient into the tangent to the mesh surface (see Rognes et al. (2013) for
implementation details).

When $v=u$, we have the “vector invariant form”,

$\frac{\partial u}{\partial t}+\underbrace{(\nabla\times u)\times u}_{\mbox{or
}(\nabla^{\perp}\cdot u)u^{\perp}}+\frac{1}{2}\nabla|u|^{2}=0,$ (113)

which is particularly useful on the surface of the sphere because it avoids
the need to rotate upwinded vectors. To see this, we multiply by a test
function $w\in V$ and integrate over one cell, integrating by parts to get

$\int_{e}w\cdot\frac{\partial u}{\partial t}-\nabla^{\perp}(w\cdot
u^{\perp})\cdot u-\nabla\cdot
w\frac{1}{2}|u|^{2}\operatorname{d}x+\int_{\partial e}w\cdot
u^{\perp}n^{\perp}\cdot\tilde{u}\operatorname{d}S=0,\quad\forall w\in V.$
(114)

Here, no rotation is required because the tangent to the edge between cells
agrees on both sides, it is just the facet normal (the normal to the cell edge
that is tangential to the cell surface) that can change on manifold meshes.
Summing over all of the cells in the mesh gives

$\int_{\Omega}w\cdot\frac{\partial u}{\partial t}-\nabla^{\perp}(w\cdot
u^{\perp})\cdot u-\nabla\cdot
w\frac{1}{2}|u|^{2}\operatorname{d}x-\int_{\Gamma}\left[\\!\\!\left[w\cdot
u^{\perp}\right]\\!\\!\right]\cdot\tilde{u}^{\perp}\operatorname{d}S=0,$ (115)

where for scalars $\phi$,
$\left[\\!\\!\left[\phi\right]\\!\\!\right]=\phi^{+}n^{+}+\phi^{-}n^{-}$. This
upwinded vector invariant form for $H(\operatorname{div})$ spaces first
appeared in Natale and Cotter (2018) for the incompressible Euler equations
and was used for the rotating shallow water equations on the sphere in Gibson
et al. (2019).

The polynomial spaces for lowest order RT elements do not span all linear
vector fields. This means that if we use the above scheme then it will only be
first order accurate. Bendall and Wimmer (2023) looked at using an auxiliary
field $q\in\mathbb{V}^{0}_{h}$ with $q=-\delta_{0}v$, _i.e._ $q$ approximates
$\nabla^{\perp}\cdot v$ if the solution domain $\Omega$ has no boundary. Then,
we can use this in an approximation of (112),

$\left\langle w,\frac{\partial v}{\partial t}\right\rangle+\left\langle
w,qu^{\perp}\right\rangle-\frac{1}{2}\left\langle\nabla\cdot w,u\cdot
v\right\rangle+G^{\prime}(v;w),\quad\forall w\in\mathbb{V}_{h}^{1}.$ (116)

This was inspired by the energy enstrophy conserving schemes that we discuss
in Section 7. Here, $G^{\prime}$ represents the discretisation of the last two
terms on the left hand side of (112), which we do not go into here (standard
upwind discontinuus Galerkin approaches were used, similar to those above). To
obtain the dynamics for $q$, we can select $w=-\nabla^{\perp}\gamma$ for
$\gamma\in\mathbb{V}_{h}^{0}$, and substitute into (116) to obtain

$\left\langle\gamma,\frac{\partial q}{\partial
t}\right\rangle-\left\langle\nabla\gamma,qu\right\rangle-G^{\prime}(v;\nabla^{\perp}\gamma),\quad\forall\gamma\in\mathbb{V}_{h}^{0}.$
(117)

We recognise the first two terms as the standard continuous finite element
approximation of $\partial q/\partial t+u\cdot\nabla q$. As we discussed
above, some method of stabilisation is usually needed to suppress oscillations
in this approximation. Modifying the test function according to the SUPG
approach is ungainly here, because of the surface terms in $G^{\prime}$.
Instead, Bendall and Wimmer (2023) used a residual based approach, writing

$\left\langle\gamma,\frac{\partial q}{\partial
t}\right\rangle-\left\langle\nabla\gamma,q^{*}u\right\rangle-G^{\prime}(v;\nabla^{\perp}\gamma),\quad\forall\gamma\in\mathbb{V}_{h}^{0},$
(118)

where

$q^{*}=q-\tau\left(\frac{\partial q}{\partial
t}+\nabla\cdot(uq)+\frac{1}{2}\nabla_{h}^{\perp}\left((\nabla_{h}v)^{T}u-(\nabla_{h}u)^{T}v\right)\right),$
(119)

where $\tau$ is a stabilisation parameter. We note that this change preserves
the consistency of the discretisation since the quantity in the brackets is
just the curl ($\nabla^{\perp}\cdot$) of (112). Then, the equation for $v$
becomes

$\left\langle w,\frac{\partial v}{\partial t}\right\rangle+\left\langle
w,q^{*}u^{\perp}\right\rangle-\frac{1}{2}\left\langle\nabla\cdot w,u\cdot
v\right\rangle+G^{\prime}(v;w),\quad\forall w\in\mathbb{V}_{h}^{1}.$ (120)

Bendall and Wimmer (2023) showed in numerical experiments (including a
nonlinear rotating shallow water equation test case on the sphere) that this
discretisation produces second order accurate solutions using lowest order RT
quadrilateral elements on a cubed sphere grid, whilst (115) only produces
first order accurate solutions. A predecessor of this scheme was considered in
Kent et al. (2022), but without the consistent definition of $q$ and $u$.

Transport schemes for $\mathbb{W}_{h}^{1}$ are discussed in Wimmer and Tang
(2022) in the context of magnetohydrodynamics, using similar ideas to those
discussed in this section for $\mathbb{W}_{h}^{2}$.

###  4.4 Temperature space transport schemes

The temperature space $\mathbb{W}^{\theta}$ produces similar challenges, since
it is continuous in the vertical and discontinuous in the horizontal. Yamazaki
et al. (2017) combined an upwind discontinuous Galerkin discretisation in the
horizontal with an SUPG discretisation in the vertical, with the modification
of test functions $\gamma\mapsto\gamma+\tau\hat{k}\cdot
u\frac{\partial\tau}{\partial z}$, producing a second order scheme when the
$RT_{1}$ discrete de Rham complex is used, so that temperature is continuous
quadratic in the vertical and discontinuous linear in the horizontal.

In staggered grid weather models, it is standard practice to collocate
thermodynamic tracers such as moisture, _etc._ , with temperature. This makes
it easier to localise thermodynamic processes that alter, and depend on, the
temperature (see Bush et al. (2020, for example)). Hence, we need to use
$\mathbb{W}_{h}^{\theta}$ transport schemes for those tracers as well. For
many of these tracers (moisture, chemical species, _etc._) it is important to
avoid numerical over- and undershoots leading to negative humidity, for
example. Hence, it is important to be able to incorporate limiters into
$\mathbb{W}_{h}^{\theta}$ transport schemes. Cotter and Kuzmin (2016) proposed
such a scheme for the vertically quadratic, horizontally linear
$\mathbb{W}_{h}^{\theta}$ also considered by Yamazaki et al. (2017). In that
scheme, at the beginning of the timestep, the vertical continuity conditions
are relaxed, and an upwind discontinuous Galerkin transport scheme is applied
over one timestep in $\hat{\mathbb{W}}_{h}^{\theta}$, the corresponding
discontinuous space. A slope limiter, such as the one in Kuzmin (2013), can
then be used to avoid under- and overshoots in this step. Then, an element
based flux corrected remapping is used to transform $\theta$ back to the
vertically continuous space $\mathbb{W}_{h}^{\theta}$. The flux correction
switches between a high order and low order mapping into order to maximise the
use of the high order solution unless under- or overshoots would otherwise
occur.

###  4.5 Recovered space schemes

The lowest order RT de Rham complex on hexahedra is attractive because it
allows storage of field values at the same grid locations as for the C grid
finite difference approximation. This is why this de Rham complex is being
used for the Met Office “Gung Ho” dynamical core (Melvin et al. 2019).
However, as we have already mentioned above, standard upwind finite element
schemes on these spaces are only first order accurate, because only
$\mathbb{W}_{h}^{0}$ has element shape functions that span the complete linear
polynomial space; $\mathbb{W}_{h}^{i}$, $i=1,2,3$, and
$\mathbb{W}_{h}^{\theta}$ do not. Bendall et al. (2019) chose to address this
by using recovery operators to construct higher order approximations of the
solution based on averaging cell values around vertices. If the original low
order solution is obtained by interpolating a smooth function, this recovery
step produces a higher order continuous finite element solution (Georgoulis
and Pryer 2018). Following Cotter and Kuzmin (2016), they then relaxed the
continuity of the higher order recovered solution and applied a discontinuous
Galerkin transport scheme step before remapping back to the original low order
finite element spaces. This was demonstrated in numerical experiments to
produce second order convergence of solutions. Bendall and Wimmer (2023)
introduced modifications to extend this approach to the sphere.

This recovery process also allows the introduction of limiters to prevent
over- and undershoots. This was done in Bendall et al. (2020), applied to
compressible Euler solutions with moisture, where limiters are critical for
stability. This produced the first atmospheric simulations using compatible
finite elements with moist physics. Bendall et al. (2022) then showed how to
achieve this framework whilst conserving mass and total moisture.

##  5 Example discretisations and iterative solution strategies

In this section, we survey some compatible finite element discretisations of
geophysical fluid dynamics models, concentrating on approaches that can be
considered as evolutions of existing approaches using more “traditional”
discretisations. More advanced structure preserving discretisations are
discussed in Sections 6-9.

###  5.1 Rotating shallow water equations on the sphere

We start with the rotating shallow water equations on the sphere, written

$\displaystyle\frac{\partial u}{\partial
t}+\mathbb{P}_{S}\left((u\cdot\nabla)u\right)+fu^{\perp}+g\nabla(D+b)$
$\displaystyle=0,$ (121) $\displaystyle\frac{\partial D}{\partial
t}+\nabla\cdot(Du)$ $\displaystyle=0,$ (122)

where $u$ is the horizontal velocity tangential to the sphere, $D$ is the
depth of the fluid layer, and $b$ is the height of the bottom surface.

Gibson et al. (2019) introduced a spatial discretisation built around the
vector invariant formulation (115) for velocity advection $u$, and the
discontinuous Galerkin formulation (105) for depth $D$. This leads a spatial
discretisation seeking $(u,D)\in\mathbb{V}_{h}^{1}\times\mathbb{V}_{h}^{2}$
such that

$\displaystyle\left\langle w,\frac{\partial u}{\partial
t}\right\rangle-\left\langle\nabla^{\perp}_{h}(w\cdot
u^{\perp}),u\right\rangle+\left\llangle\left[\\!\\!\left[w\cdot
u^{\perp}\right]\\!\\!\right],\tilde{u}^{\perp}\right\rrangle_{\Gamma}$
$\displaystyle\qquad-\nabla\cdot
w\left(\frac{1}{2}|u|^{2}+g(D+b)\right)\operatorname{d}x$
$\displaystyle=0,\quad\forall w\in\mathbb{V}_{h}^{1},$ (123)
$\displaystyle\left\langle\phi,D_{t}\right\rangle-\left\langle\nabla_{h}\phi,uD\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}.$ (124)

They used a semi implicit timestepping scheme, which is best described as some
form of iteration towards the fully implicit timestepping scheme

$\displaystyle\left\langle w,u^{n+1}-u^{n}\right\rangle-\Delta
t\left\langle\nabla^{\perp}_{h}(w\cdot\bar{u}^{\perp}),(u^{n+1/2})^{\perp}\right\rangle$
$\displaystyle\qquad+\Delta
t\left\llangle\left[\\!\\!\left[w\cdot\bar{u}^{\perp}\right]\\!\\!\right],(\tilde{u}^{n+1/2})^{\perp}\right\rrangle_{\Gamma}$
$\displaystyle\qquad-\Delta t\nabla\cdot
w\left(\frac{1}{2}|\bar{u}|^{2}+g(\bar{D}+b)\right)\operatorname{d}x$
$\displaystyle=0,\quad\forall w\in\mathbb{V}_{h}^{1},$ (125)
$\displaystyle\left\langle\phi,D_{t}\right\rangle-\Delta
t\left\langle\nabla_{h}\phi,\bar{u}D\right\rangle+\Delta
t\left\llangle\left[\\!\\!\left[\phi\bar{u}\right]\\!\\!\right],\tilde{D}\right\rrangle_{\Gamma}$
$\displaystyle=0,\quad\forall\phi\in\mathbb{V}_{h}^{2},$ (126)

where $u^{n+1/2}=\bar{u}=(u^{n+1}+u^{n})/2$. This mixture of two symbols for
the same thing is introduced to describe an iteration based around this
implicit discretisation. We write $v^{0}=u^{n},v^{1},v^{2},\ldots$ and
$D^{0}=D^{n},D^{1},D^{2},\ldots$ as a sequence of iterative approximations to
$u^{n+1}$ and $u^{n}$ respectively. For each iteration $k$, we write
$\bar{u}=(u^{n}+v^{k})/2$, $\bar{D}=(D^{n}+D^{k})/2$, solving (125-126) for
$u^{n+1}$ and $D^{n+1}$ using those values of $\bar{u}$ and $\bar{D}$. Then,
we use the linearisation about the state of rest to compute iterative
corrections $(\Delta u,\Delta D)$, according to

$\displaystyle\left\langle w,\Delta u\right\rangle+\frac{\Delta
t}{2}\left\langle w,f\Delta u^{\perp}\right\rangle$
$\displaystyle\qquad-\frac{\Delta t}{2}\left\langle\nabla\cdot w,g\Delta
D\right\rangle$ $\displaystyle=-R_{u}[w]:=-\left\langle
w,u^{n+1}-u^{n}\right\rangle,\quad\forall w\in\mathbb{V}_{h}^{1},$ (127)
$\displaystyle\left\langle\phi,\Delta D+\frac{H\Delta t}{2}\nabla\cdot\Delta
u\right\rangle$
$\displaystyle=-R_{D}[\phi]:=-\left\langle\phi,D^{n+1}-D^{n}\right\rangle,\quad\forall\phi\in\mathbb{V}_{h}^{2}.$
(128)

We will discuss the solution of this linear system later. The time integration
scheme applies a fixed number of iterations of this type (typically $2\leq
k_{\max}\leq 4$). We then take $(u^{n+1},D^{n+1})=(v^{k_{\max}},D^{k_{\max}})$
before moving to the next timestep. It is not intended to converge to the
solution of the implicit midpoint rule but just to produce a second order semi
implicit scheme that is stable conditional on the advective Courant number
$|u|\Delta t/\Delta x$, but unconditionally in the wave Courant number
$\sqrt{gH}\Delta t/\Delta x$. A probably more stable approach is to use the
form (112) for the velocity equation, substituting $\bar{u}$ in for $u$, and
$u$ in for $v$. This was done using quadrilateral $RT_{0}$ elements in Bendall
and Wimmer (2023). Another approach is to replace the implicit midpoint rule
update for $u^{n+1}$ and $D^{n+1}$ given $\bar{u}$, instead using an explicit
transport step (or several substeps). This can facilitate more sophisticated
transport schemes with limiters that are hard to implement in implicit
schemes. This was also done using quadrilateral $RT_{0}$ elements in Bendall
and Wimmer (2023).

###  5.2 Rotating incompressible Boussinesq equations

For ocean models, the most common setting is the incompressible Boussinesq
equations,

$\displaystyle\frac{\partial u}{\partial t}+(u\cdot\nabla)u+2\Omega\times
u+\nabla p-b\hat{k}$ $\displaystyle=0,$ (129) $\displaystyle\frac{\partial
b}{\partial t}+(u\cdot\nabla)b$ $\displaystyle=0,$ (130)
$\displaystyle\nabla\cdot u$ $\displaystyle=0,$ (131)

where $\Omega$ is the Earth’s rotation rate, $p$ is the pressure, and $b$ is
the buoyancy. Here we limit discussion to the rigid lid approximation with
boundary conditions $u\cdot n=0$ on all boundaries. In full ocean models there
are generally two thermodynamic tracers, the potential temperature $\theta$
and the salinity $S$, which are both transported by advection equations as $b$
is above, and $b$ is then a specified function of $\theta$ and $S$ _via_ an
equation of state. However, here we keep things to the simple formulation
above. Further, many models make the hydrostatic approximation, but we do not
discuss that here. Finally, these equations are extended in operational models
to include mixing parameterisations and representations of other physical
processes; we do not discuss those either.

Yamazaki et al. (2017) introduced a hybrid approach for velocity advection in
the vertical slice setting (the velocity is three dimensional but all fields
are independent of $y$, so the equations can be solved on a two dimensional
mesh in the $x-z$ plane), with (115) as the transport scheme for the $x-z$
components of velocity, and a standard upwind discontinuous Galerkin scheme
for the $y$ component. For bouyancy $b$ (represented in $\mathbb{W}^{\theta}$)
the hybrid scheme with upwind discontinuous Galerkin in the horizontal and
SUPG in the vertical was used. This spatial discretisation was combined with a
similar timestepping scheme to the one above, resulting in a linear system of
the form,

$\displaystyle\left\langle w,\Delta u\right\rangle+\frac{\Delta
t}{2}\left\langle w,2\Omega\times\Delta u\right\rangle$
$\displaystyle\qquad-\frac{\Delta t}{2}\left\langle\nabla\cdot w,\Delta
p\right\rangle-\left\langle w,b\hat{k}\right\rangle$
$\displaystyle=-R_{u}[w],\quad\forall w\in\mathbb{W}_{h}^{2},$ (132)
$\displaystyle\left\langle\gamma,\Delta b\right\rangle+\frac{\Delta
t}{2}\left\langle\gamma,\Delta u\cdot\hat{k}B_{z}\right\rangle$
$\displaystyle=-R_{b}[\gamma],\quad\forall\gamma\in\mathbb{W}_{h}^{\theta},$
(133) $\displaystyle\left\langle\phi,\nabla\cdot\Delta u\right\rangle$
$\displaystyle=0,\quad\forall\phi\in\mathbb{W}_{h}^{3},$ (134)

where $B_{z}$ is the vertical derivative of a reference buoyancy profile, to
compute the iterative linear corrections to $u^{n+1}$, $p^{n+1}$ (actually an
approximation to pressure at time level $t^{n+1/2}$) and $b^{n+1}$,
analogously to (127\- 128).

Equations (132-134) were solved by eliminating $\Delta b$. This is possible
without introducing errors on an extruded mesh with flat layers when $B_{z}$
is constant, since $\Delta b\hat{k}\in\mathbb{W}_{h}^{2,V}$. When that is not
the case, the errors from the approximate elimination can be removed by using
the elimination as a preconditioner for a Krylov method on (132-134). At the
time, the equations were modified by setting $\Omega$ to zero on the left hand
side (which prevented the iterative solver from being robust to large values
of $\Omega$), and using GMRES applied to the coupled system,

$\displaystyle\left\langle w,\Delta u+\frac{\Delta
t^{2}}{4}\hat{k}\hat{k}\Delta uB_{z}\right\rangle+\frac{\Delta
t}{2}\left\langle w,2\Omega\times\Delta u\right\rangle$
$\displaystyle\qquad-\frac{\Delta t}{2}\left\langle\nabla\cdot w,\Delta
p\right\rangle-\left\langle w,b\hat{k}\right\rangle$
$\displaystyle=-\tilde{R}_{u}[w],\quad\forall w\in\mathbb{W}_{h}^{2},$ (135)
$\displaystyle\nabla\cdot\Delta u$ $\displaystyle=0.$ (136)

This was preconditioned by an $H(\operatorname{div})$ block preconditioner,
which we do not describe here. However, it is much better to precondition
(5.2-136) using a hybridised solver, which we describe later in this section.

Yamazaki et al. (2017) showed that this suite of discretisation and solver
choices produces a scheme that can resolve fronts in the Eady vertical slice
frontogenesis problem to a similar degree as C grid finite difference methods
used previously. After this paper was published, the authors experimented with
replacing the hybrid advection scheme for velocity with a full vector
invariant form in all 3 components of velocity. In this case, oscillations in
the velocity field quickly emerge when the front sharpens. These oscillations
do not appear when (110) is used instead. It is possible that this is related
to the Hollingsworth instability associated to the vector invariant form used
with C grid finite difference methods (Hollingsworth et al. 1983), but this
requires further investigation.

###  5.3 Rotating compressible Euler equations

For global atmosphere models, a standard approach is to solve the rotating
compressible Euler equations, given by

$\displaystyle u_{t}+(u\cdot\nabla)u+2\Omega\times
u+c_{p}\theta\nabla\Pi+g\hat{k}$ $\displaystyle=0,$ (137)
$\displaystyle\theta_{t}+u\cdot\nabla\theta$ $\displaystyle=0,$ (138)
$\displaystyle D_{t}+\nabla\cdot(uD)$ $\displaystyle=0,$ (139)
$\displaystyle\Pi$ $\displaystyle=E(\theta,D),$ (140)

where $\theta$ is the potential temperature, $\Pi$ is the Exner pressure, $D$
is now the density for 3D models, $c_{p}$ is the specific heat at constant
pressure (a constant parameter in the ideal gas law), and $E$ is a prescribed
function describing the thermal equation of state relating $\Pi$, $D$ and
$\theta$. This form of the equations is known as the “theta-Pi” formulation.
Other thermodynamic formulations make use of pressure and temperature
directly, or other combinations of variables, but we do not discuss them here.
For simplicity we consider boundary conditions $u\cdot n=0$ on the bottom and
top boundaries of the domain, although representations of the top of the
atmosphere can be rather more complicated (since the real atmosphere has a
density that decreases with height until it is so low that assumptions
underlying the fluid dynamics model do not hold).

A compatible finite element formulation of the compressible Euler equations
uses $u\in\mathbb{W}_{h}^{2}$, $D\in\mathbb{W}_{h}^{3}$, and
$\theta\in\mathbb{W}_{h}^{\theta}$. Either $\Pi$ is solved as an independent
variable in $\mathbb{W}_{h}^{3}$, with the equation of state being projected
into $\mathbb{W}_{h}^{3}$, or $\Pi$ can just be replaced by $E(\theta,D)$ in
the velocity equation, leading to a system for $u$, $\theta$ and $D$.

Natale et al. (2016) proposed a formulation using (115) as the transport
scheme for velocity, upwind discontinuous Galerkin for density transport and
the hybrid scheme for potential temperature used in Yamazaki et al. (2017).
The main additional challenge is the discretisation of the pressure gradient
term $-\theta\nabla\Pi$ (which is scaled by $c_{p}$). In Yamazaki et al.
(2017), the pressure gradient $\nabla p$ was integrated by parts, but this is
more complicated in the compressible Euler case because of the presence of
$\theta\in\mathbb{W}_{h}^{\theta}$, which can have discontinuities in the
horizontal direction across vertical facets. We need to integrate by parts
because $\Pi$ is discontinuous, whether it is an independent variable in
$\mathbb{W}_{h}^{3}$ or the evaluation of $E(D,\theta)$ (since
$D\in\mathbb{W}_{h}^{3}$ is discontinuous). Natale et al. (2016) proposed to
apply integration by parts separately in each cell $e$ for this term,
obtaining

$\int_{e}\nabla\cdot(w\theta)\Pi\operatorname{d}x-\int_{\partial e}\theta
w\cdot n\\{\Pi\\}\operatorname{d}S,\quad\forall w\in\mathbb{W}_{h}^{2}(e),$
(141)

where $\\{\Pi\\}$ is the average value of $\Pi$ between the inside and the
outside of $e$ (since $\Pi$ takes two values on $\partial e$). Summing this
over all cells gives

$\left\langle\nabla_{h}\cdot(w\theta),\Pi\right\rangle-\left\llangle\left[\\!\\!\left[w\theta\right]\\!\\!\right],\\{\Pi\\}\right\rrangle_{\Gamma},\quad\forall
w\in\mathbb{W}_{h}^{2}(\Omega),$ (142)

where $\\{\Pi\\}$ is now defined as $(\Pi^{+}+\Pi^{-})/2$. Again, analogously
to (127\- 128), a semi implicit timestepping scheme can be used, this time
built around a linearisation about a state of rest with $u=0$,
$\theta=\bar{\theta}$ and $\Pi=\bar{\Pi}$. These reference profiles vary in
the vertical only (or following the ENDGame approach (Wood et al. 2014), the
values of $\bar{\theta}$ and $\bar{\Pi}$ can be used from $\theta$ and $\Pi$
at the previous timestep, but whilst neglecting their horizontal derivatives
in the linearisation to facilitate efficient solution). Finally, we neglect
horizontal derivatives of $\Delta\theta$ appearing in the
$\Delta\theta\nabla\bar{\Pi}$ and term, for the same reason. This results in
the following linear iteration, presented in Gibson (2019) (and used for the
numerical results in Natale et al. (2016)),

$\displaystyle\left\langle w,\Delta u\right\rangle+\frac{\Delta
t}{2}\left\langle w,2\Omega\times\Delta u\right\rangle$
$\displaystyle-\frac{c_{p}\Delta
t}{2}\left\langle\nabla_{h}\cdot(\bar{\theta}w),\Delta\Pi\right\rangle+\frac{c_{p}\Delta
t}{2}\left\llangle\left[\\!\\!\left[\bar{\theta}w\right]\\!\\!\right],\\{\Delta\Pi\\}\right\rrangle_{\Gamma}$
$\displaystyle\qquad-\frac{c_{p}\Delta
t}{2}\left\langle\nabla\cdot(\hat{k}\Delta\theta\theta
w\cdot\hat{k}),\bar{\Pi}\right\rangle$ $\displaystyle=-R_{u}[w],\quad\forall
w\in\mathbb{W}_{h}^{2},$ (143)
$\displaystyle\left\langle\gamma,\Delta\theta\right\rangle+\frac{\Delta
t}{2}\left\langle\gamma,\frac{\partial\bar{\theta}}{\partial z}\Delta
u\cdot\hat{k}\right\rangle$
$\displaystyle=-R_{\theta}[\gamma],\quad\forall\gamma\in\mathbb{W}_{h}^{\theta},$
(144) $\displaystyle\left\langle\phi,\Delta D\right\rangle-\frac{\Delta
t}{2}\left\langle\nabla_{h}\phi,\bar{D}\Delta u\right\rangle+\frac{\Delta
t}{2}\left\llangle\left[\\!\\!\left[\phi\Delta
u\right]\\!\\!\right],\\{\bar{D}\\}\right\rrangle_{\Gamma}$
$\displaystyle=-R_{D}[\phi],\quad\forall\phi\in\mathbb{W}_{h}^{3},$ (145)
$\displaystyle\Delta\Pi$ $\displaystyle=\frac{\partial
E}{\partial\theta}\Delta\theta+\frac{\partial E}{\partial D}\Delta D,$ (146)

for the iterative updates $(\Delta u,\Delta\theta,\Delta
D)\in\mathbb{W}_{h}^{2}\times\mathbb{W}_{h}^{\theta}\times\mathbb{W}_{h}^{3}$
to $(u^{n+1},\theta^{n+1},D^{n+1})$. Alternatively, as was done in Melvin et
al. (2019), we can add $\Delta\Pi\in\mathbb{W}^{3}_{h}$ to the list of
independent variables and replace (146) with

$\left\langle\alpha,\Delta\Pi\right\rangle=\left\langle\alpha,\frac{\partial
E}{\partial\theta}\Delta\theta+\frac{\partial E}{\partial D}\Delta
D\right\rangle,\quad\forall\alpha\in\mathbb{W}_{h}^{3}.$ (147)

Melvin et al. (2019) adopted a hybrid approach, using finite volume methods to
approximate the transport terms in a discretisation otherwise built using
compatible finite element methods.

###  5.4 Iterative solver strategies

Now we focus on the iterative solver strategies for these linear implicit
systems. In all of the strategies we discuss here, the temperature
$\Delta\theta$ is first (approximately) eliminated following our description
of the approach to incompressible Boussinesq equations discussed above. This
leads to the following system

$\displaystyle\left\langle w,\Delta u\right\rangle+\frac{\Delta
t}{2}\left\langle w,2\Omega\times\Delta u\right\rangle$
$\displaystyle-\frac{c_{p}\Delta
t}{2}\left\langle\nabla\cdot(\bar{\theta}w),\Delta\Pi\right\rangle+\frac{c_{p}\Delta
t}{2}\left\llangle\left[\\!\\!\left[\bar{\theta}w\right]\\!\\!\right],\\{\Delta\Pi\\}\right\rrangle_{\Gamma}$
$\displaystyle\qquad-\frac{c_{p}\Delta
t}{2}\left\langle\nabla\cdot(\hat{k}\Delta\theta\theta
w\cdot\hat{k}),\bar{\Pi}\right\rangle$ $\displaystyle=-R_{u}[w],\quad\forall
w\in\mathbb{W}_{h}^{2},$ (148) $\displaystyle\left\langle\phi,\Delta
D\right\rangle-\frac{\Delta t}{2}\left\langle\nabla_{h}\phi,\bar{D}\Delta
u\right\rangle+\frac{\Delta t}{2}\left\llangle\left[\\!\\!\left[\phi\Delta
u\right]\\!\\!\right],\\{\bar{D}\\}\right\rrangle_{\Gamma}$
$\displaystyle=-R_{D}[\phi],\quad\forall\phi\in\mathbb{W}_{h}^{3},$ (149)
$\displaystyle\Delta\theta$ $\displaystyle=-\frac{\Delta
t}{2}\bar{\theta}\Delta u\cdot\hat{k}+r_{\theta},$ (150)
$\displaystyle\Delta\Pi$ $\displaystyle=\frac{\partial
E}{\partial\theta}\Delta\theta+\frac{\partial E}{\partial D}\Delta D,$ (151)

for $(\Delta u,\Delta D)\in\mathbb{W}_{h}^{2}\times\mathbb{W}_{h}^{3}$, where
$r_{\theta}\in\mathbb{W}_{h}^{\theta}$ such that

$\left\langle
r_{\theta},\gamma\right\rangle=R_{\theta}[\gamma],\quad\forall\gamma\in\mathbb{W}_{h}^{\theta}.$
(152)

In other words, $r_{\theta}$ is the $L^{2}$ Riesz representer of $R_{\theta}$.
Mitchell and Müller (2016) proposed to solve this reduced system using GMRES
with a Schur complement preconditioner, using an approximate Schur complement
formed from the lumped velocity mass matrix (and setting $\Omega=0$). This was
incorporated into a horizontal multigrid scheme (coarsening the mesh in the
horizontal but not the vertical) using line smoothers (direct solves
neglecting horizontal coupling between columns) for the approximate Schur
complement on the levels. This combination of horizontal multigrid and
vertical line smoothers is necessary because of the small aspect ratio of the
atmosphere, and scalable parallel performance was observed over a large range
of resolutions. This solve approach was successfully implemented in the Met
Office system using the discretisation approach of Melvin et al. (2019) by
Maynard et al. (2020).

There is an alternative solution approach that has been applied to compatible
discretisations of elliptic problems since the mid 20th Century:
hybridisation. In hybridisation, the continuity constraints of the
$H(\operatorname{div})$ space are relaxed, and are enforced through Lagrange
multipliers as part of the solution formulation. The Lagrange multiplier space
$\operatorname{Tr}$, known as the trace space since it is supported only on
cell facets, is chosen to match the $H(\operatorname{div})$ space when
restricted to a facet and dotted with the normal component. This means that it
is discontinuous between facets that meet at a vertex in two dimensions, or at
an edge or vertex in three dimensions. For example, the hybridisable
formulation of (127-128) seeks $(\Delta u,\Delta
D,\lambda)\in\hat{\mathbb{V}}_{h}^{1}\times\mathbb{V}_{h}^{2}\times\operatorname{Tr}(\mathbb{V}_{h}^{1})$
such that

$\displaystyle\left\langle w,\Delta u\right\rangle+\frac{\Delta
t}{2}\left\langle w,f\Delta u^{\perp}\right\rangle$
$\displaystyle\quad-\frac{\Delta t}{2}\left\langle\nabla\cdot w,g\Delta
u\right\rangle+\left\llangle\left[\\!\\!\left[w\right]\\!\\!\right],\lambda\right\rrangle_{\Gamma}$
$\displaystyle=-\tilde{R}_{u}[w],\quad\forall w\in\hat{\mathbb{V}}_{h}^{1},$
(153) $\displaystyle\left\langle\phi,\Delta D+\frac{H\Delta
t}{2}\nabla\cdot\Delta u\right\rangle$
$\displaystyle=-R_{D}[\phi],\quad\forall\phi\in\mathbb{V}_{h}^{2},$ (154)
$\displaystyle\left\llangle\gamma,\left[\\!\\!\left[\Delta
u\right]\\!\\!\right]\right\rrangle_{\Gamma}$
$\displaystyle=0,\quad\forall\gamma\in\operatorname{Tr}(\mathbb{V}_{h}^{1}),$
(155)

where

$\tilde{R}_{u}[w]=R_{u}[w],\quad\forall w\in\mathbb{V}_{h}^{1}.$ (156)

To see that this is an equivalent formulation to (127-128), note that
$\hat{\mathbb{V}}_{h}^{1}\subset\mathbb{V}_{h}^{1}$, so we may choose
$w\in\mathbb{V}_{h}^{1}$ in (153). In that case, the $\lambda$ term vanishes
because $\left[\\!\\!\left[w\right]\\!\\!\right]=0$, and we recover 127.
Further, (155) ensures that $u\in\mathbb{V}_{h}^{1}$ at the solution, since
taking $\gamma=\left[\\!\\!\left[\Delta u\right]\\!\\!\right]$ implies that
$\left[\\!\\!\left[\Delta u\right]\\!\\!\right]=0$ in $L^{2}(\Gamma)$; we note
that

$\mathbb{V}_{h}^{1}(\Omega)=\\{u\in\hat{\mathbb{V}}_{h}^{1}(\Omega):\|\left[\\!\\!\left[u\right]\\!\\!\right]\|_{L^{2}(\Gamma)}=0\\}.$
(157)

The advantage of this formulation is that $\Delta u$ and $\Delta D$ can now
both be eliminated elementwise, leading to a sparse system for $\lambda$; the
reduced system is referred to as the hybridised system. This is possible
because we can take $w$ and $\phi$ supported in only one cell, and then given
$\lambda$, we can solve for $\Delta u$ and $\Delta D$ independently in each
cell (this is referred to as the local solver). Material discussing the well
posedness of the hybridised system is surveyed in Boffi et al. (2013), along
with postprocessing techniques for obtaining improved approximations using
$\lambda$. A nonrigorous intuitive explanation for this is that $\lambda$
gives an approximation of $\Delta D$ evaluated on facets, and hence the
hybridised equation has properties of an approximation to the Helmholtz
equation satisfied by $\Delta D$ after eliminating $\Delta u$ from the linear
PDE. This idea is built upon in Cockburn and Gopalakrishnan (2004), which
provides an explicit weak form characterisation of the hybridised method
(including the nonsymmetric term here containing $\Omega\times u$ is a
straightforward extension of that work). This idea was used to demonstrate
smoothing properties of standard iterative methods in Gopalakrishnan (2003)
and a convergent multigrid scheme in Gopalakrishnan and Tan (2009). Gibson et
al. (2020) applied the hybridisation technique to the linear compressible
Boussinesq equations, which are a minor modification of the linearisation of
the incompressible Boussinesq equations described above, incorporating linear
acoustic waves.

This strategy cannot be applied to (148 -151), because the averaging
$\\{\Delta\Pi\\}$ of $\Pi$ on facets couples the values of $\Delta\Pi$ between
cells (and so there is no local solver). Gibson (2019) proposed a modification
to address this, in which we seek $(\Delta u,\Delta
D,\lambda)\in\hat{\mathbb{W}}_{h}^{2}\times\mathbb{W}_{h}^{3}\times\operatorname{Tr}(\mathbb{W}_{h}^{2})$
such that

$\displaystyle\left\langle w,\Delta u\right\rangle+\frac{\Delta
t}{2}\left\langle w,2\Omega\times\Delta u\right\rangle$
$\displaystyle-\frac{c_{p}\Delta
t}{2}\left\langle\nabla_{h}\cdot(\bar{\theta}w),\Delta\Pi\right\rangle\qquad+\frac{c_{p}\Delta
t}{2}\left\llangle\left[\\!\\!\left[\bar{\theta}w\right]\\!\\!\right],\lambda\right\rrangle_{\Gamma}$
$\displaystyle\qquad-\frac{c_{p}\Delta
t}{2}\left\langle\nabla\cdot(\hat{k}\Delta\theta\theta
w\cdot\hat{k}),\bar{\Pi}\right\rangle$
$\displaystyle\qquad+\underbrace{\frac{c_{p}\Delta
t}{2}\left\llangle\left[\\!\\!\left[\hat{k}\Delta\theta\theta
w\cdot\hat{k}\right]\\!\\!\right],\\{\Pi\\}\right\rrangle_{\Gamma_{H}}}_{\star}$
$\displaystyle=-R_{u}[w],\quad\forall w\in\hat{\mathbb{W}}_{h}^{2},$ (158)
$\displaystyle\left\langle\phi,\Delta D\right\rangle-\frac{\Delta
t}{2}\left\langle\nabla\phi,\bar{D}\Delta u\right\rangle+\frac{\Delta
t}{2}\left\llangle\left[\\!\\!\left[\phi\Delta
u\right]\\!\\!\right],\\{\bar{D}\\}\right\rrangle_{\Gamma}$
$\displaystyle=-R_{D}[\phi],\quad\forall\phi\in\mathbb{W}_{h}^{3},$ (159)
$\displaystyle\left\llangle\gamma,\left[\\!\\!\left[u\right]\\!\\!\right]\right\rrangle_{\Gamma}$
$\displaystyle=0,\quad\forall\gamma\in\operatorname{Tr}(\mathbb{W}_{h}^{2}),$
(160) $\displaystyle\Delta\theta$ $\displaystyle=-\frac{\Delta
t}{2}\bar{\theta}\Delta u\cdot\hat{k},$ (161) $\displaystyle\Delta\Pi$
$\displaystyle=\frac{\partial E}{\partial\theta}\Delta\theta+\frac{\partial
E}{\partial D}\Delta D,$ (162)

where $\Gamma_{H}$ is the set of horizontal faces between cells in the same
vertical column. Here, the idea is that $\lambda$ is an approximation of
$c_{p}\Delta t\Delta\Pi/2$ on mesh facets. This system is not equivalent to
(148 -151), though. Note the addition of the term indicated with $\star$. This
term vanishes when $w\in\mathbb{W}_{h}^{2}$, so it does not change the
solution, but it was found that without it, iterative solvers do not perform
well; it appears to be required for the coercivity of the solution. This
system has not been analysed yet, but it was demonstrated to produce
comparible results to standard test cases when applied to the full nonlinear
compressible Euler equations. Scalable multigrid behaviour for the hybridised
system was also demonstrated. Bendall et al. (2020) extended this solver
approach to the compressible Euler equations with moisture, where the mass
lumping approach was found not to work well when $RT_{1}$ spaces were used.
Betteridge et al. (2022) showed that this hybridisation approach produces
scalable results for the Met Office formulation of Melvin et al. (2019).

###  5.5 Computing hydrostatic balanced states

The hybridisation approach also provides a useful way to solve for hydrostatic
balanced states, satisfying

$-\left\langle\nabla\cdot(\theta w),\Pi\right\rangle=\left\langle
w\cdot\hat{k},g\right\rangle-\left\llangle w\cdot
n,\Pi_{0}\right\rrangle\quad\forall w\in\mathring{\mathbb{W}}_{h}^{2,V}.$
(163)

Following the technique of Natale et al. (2016) previously discussed, for
given $\theta\in\mathbb{W}_{h}^{\theta}$, we seek
$(v,D)\in\mathbb{W}_{h}^{2}\times\mathbb{W}_{h}^{3}$ such that

$\displaystyle\left\langle w,v\right\rangle-\left\langle\nabla\cdot(\theta
w),\Pi\right\rangle$ $\displaystyle=\left\langle
w\cdot\hat{k},g\right\rangle-\left\llangle w\cdot
n,\Pi_{0}\right\rrangle\quad\forall w\in\mathring{\mathbb{W}}_{h}^{2,V},$
(164) $\displaystyle\left\langle\phi,\nabla\cdot v\right\rangle$
$\displaystyle=0,\quad\forall\phi\in\mathbb{W}_{h}^{3},$ (165)
$\displaystyle\Pi$ $\displaystyle=E(\theta,D),$ (166)

which is independent between columns. An equivalent hybridisable formulation
seeks $(v,D)\in\mathbb{W}_{h}^{2}\times\mathbb{W}_{h}^{3}$ such that

$\displaystyle\left\langle w,v\right\rangle-\left\langle\nabla\cdot(\theta
w),\Pi\right\rangle$
$\displaystyle\qquad+\left\llangle\left[\\!\\!\left[w\right]\\!\\!\right],\lambda\right\rrangle_{\Gamma}$
$\displaystyle=\left\langle w\cdot\hat{k},g\right\rangle-\left\llangle w\cdot
n,\Pi_{0}\right\rrangle\quad\forall w\in\hat{\mathbb{W}}_{h}^{2,V},$ (167)
$\displaystyle\left\langle\phi,\nabla\cdot v\right\rangle$
$\displaystyle=0,\quad\forall\phi\in\mathbb{W}_{h}^{3},$ (168)
$\displaystyle\left\llangle\gamma,\left[\\!\\!\left[v\right]\\!\\!\right]\right\rrangle_{\Gamma}$
$\displaystyle=0,\quad\forall\gamma\in\operatorname{Tr}(\mathbb{W}_{h}^{2,V}),$
(169) $\displaystyle\Pi$ $\displaystyle=E(\theta,D),$ (170)

where now $\Gamma$ includes the bottom boundary (but not the top), and
$\operatorname{Tr}(\mathbb{W}_{h}^{2,V})$ is only supported on horizontal
facets between neighbouring cells in the same column. This can be solved using
Newton’s method, with solution of the Jacobian system _via_ the hybridised
system for $\delta\lambda$.

###  5.6 Monolithic solvers

One more recent solver approach has been investigated in numerical experiments
in Cotter and Shipton (2022), who used a similar suite of discretisations to
that of Natale et al. (2016), with the exception of using edge stabilisation
for potential temperature instead of vertical SUPG. In this work, the fully
nonlinear implicit midpoint rule is solved using Newton’s method, and GMRES is
applied to the coupled system for $\Delta u,\Delta\theta,\Delta D$ without
elimination. The system is preconditioned by an additive Schwarz method,
computing exact solutions of the Jacobian system restricted to overlapping
columnar patches; each patch consists of the cells surrounding one vertical
edge (excluding degrees of freedom attached to the vertical facets on the side
boundaries of the patch). This scheme also requires further analysis, but was
shown to produce mesh independent iteration counts in numerical experiments.

##  6 Variational discretisations

Variational discretisations are discretisations that are derived from a
discrete Hamilton’s principle. They were originally introduced in the setting
of ordinary differential equations (ODEs), taking the name “variational
integrators”, surveyed in Marsden and West (2001). The idea behind variational
discretisations is that rather than discretising the equations directly, we
instead discretise the action functional from which the equations are derived.
In the case of ODEs, this means replacing the time integral with a discrete
quadrature rule involving the solution at discrete points in time. The
discretisation of the equations is then obtained by finding stationary points
of the discretised Lagrangian.

The advantage of variational discretisations is that if the discretised action
has symmetries, then these symmetries give rise to corresponding conserved
quantities via (the discrete) Noether’s theorem. In the context of variational
integrators for mechanical systems, this yields discrete conservation of
momentum and angular momentum, for example. Further, after Legendre
transformation the discrete timestepping map is symplectic, leading to the
conservation (up to exponentially small terms in $\Delta t$) of a modified
energy/Hamiltonian obtained through backward error analysis (Hairer et al.
2006, 2003, Leimkuhler and Reich 2004, Sanz-Serna 1992).

In principle, the variational discretisation approach can be extended directly
to partial differential equations by simply discretising the action functional
in space as well as in time (or one may consider spatial semidiscretisation by
discretising in space only, as we shall mostly do here). For fluid dynamics,
the situation is more challenging, because the underlying variational
principle is defined in terms of the Lagrangian flow map, rather than the
Eulerian quantities. This is discussed in the following subsection.

###  6.1 Hamilton’s Principle for fluid dynamics: continuous theory

The Lagrangian flow map is treated formally (we avoid discussions of
smoothness _etc._ here) as a diffeomorphism
$\chi(\cdot,t):\Omega_{0}\to\Omega$, mapping labels in a configuration space
$\Omega_{0}$ to fluid particle locations at time $t$ in the physical domain
$\Omega$. We have a time-$t$ family of maps222Here we use the notation
$\chi(\cdot,t)$ to indicate the function $x\mapsto\chi(x,t)$ for given $t$..
The variational formulation then follows by writing an action as an integral
over $\Omega_{0}$, and variations in $\chi$ are considered subject to the
usual endpoint conditions in time plus the requirement that $\chi$ be a
diffeomorphism. Working with these Lagrangian flow maps is difficult both for
theory and numerical computation; this was addressed by Arnold’s geometric
formulation in Eulerian variables for the incompressible Euler equations
(Arnold 1966). For each $(x,t)\in\Omega\times[0,T]$, where $T$ is the time
interval over which the equations are being solved, $u(x,t)$ is a vector
tangent to $x$ at $\Omega$333We have $u(x,t)\in\mathbb{R}^{N}$ for the case
where $\Omega$ is an $N$ dimensional subset of $\mathbb{R}^{N}$, but in
geophysical fluid dynamics we are also interested in solving problems where
$\Omega$ is the surface of a sphere. For each $t$, we say that
$u(\cdot,t)\in\mathfrak{X}(\Omega)$, the space of vector fields on $\Omega$.

The formulation stems from the observation that the Eulerian velocity $u$
satisfies

$\frac{\partial}{\partial t}\chi=u\circ\chi,$ (171)

where for functions of $(x,t)$ we write $(u\circ\chi)(x,t)=u(\chi(x,t),t)$. To
construct Hamilton’s principle, we need to consider perturbations to $\chi$
that are still diffeomorphisms for each $t$. If we consider such a continuous
2-parameter family of perturbed maps $\tilde{\chi}(x,t,\epsilon)$, with
$\tilde{\chi}(x,t,\epsilon=0)=\chi(x,t)$, then there exists a continuous
family of vector fields $\tilde{w}(\cdot,t,\epsilon)\in\mathfrak{X}(\Omega)$
such that

$\frac{\partial}{\partial\epsilon}\tilde{\chi}=\tilde{w}\circ\tilde{\chi},$
(172)

where we now extend the $\circ$ notation to the case of two parameters,
$\epsilon$ and $t$, writing
$(\tilde{w}\circ\tilde{\chi})(x,t,\epsilon)=\tilde{w}(\tilde{\chi}(x,t,\epsilon),t,\epsilon)$.
Differentiating with respect to $\epsilon$ and evaluating at $\epsilon=0$, we
get the infinitesimally perturbed $\chi$,

$\delta\chi:=\frac{\partial}{\partial\epsilon}|_{\epsilon=0}\tilde{\chi}=w\circ\chi,$
(173)

where $w(x,t):=\tilde{w}(x,t,\epsilon=0)$. Here $\delta\chi$ depends on the
direction of perturbation in the parameter $\epsilon$, so we always consider
$\delta\chi$ being defined with respect to a particular choice of $w$. By
differentiating (171) with respect to $\epsilon$, we obtain

$\displaystyle\delta\frac{\partial}{\partial t}\chi$ $\displaystyle=(\nabla
u)\circ\chi\cdot\delta{\chi}+\delta{u}\circ\chi,$ (174) $\displaystyle=(\nabla
u)\circ\chi\cdot w\circ\chi+\delta{u}\circ\chi,$ (175)

where

$\delta
u:=\frac{\partial}{\partial\epsilon}|_{\epsilon=0}\frac{\partial}{\partial
t}{\chi}\circ\chi^{-1}.$ (176)

Similarly, by differentiating (173) with respect to $t$, we obtain

$\displaystyle\frac{\partial}{\partial t}\delta\chi$ $\displaystyle=(\nabla
w)\circ\chi\cdot\frac{\partial}{\partial t}{\chi}+\frac{\partial}{\partial
t}{w}\circ\chi,$ (177) $\displaystyle=(\nabla w)\circ\chi\cdot
u\circ\chi+\frac{\partial}{\partial t}{w}\circ\chi.$ (178)

By subtracting Equations (175) and (178), noting the symmetry of second
derivatives and composing with $\chi^{-1}$, we obtain

$\delta u=\frac{\partial}{\partial t}{w}+\underbrace{u\cdot\nabla
w-w\cdot\nabla u}_{=[u,w]}.$ (179)

Hence, we have successfully related infinitesimal perturbations (“variations”)
in $u$ to variations in $\chi$ entirely in terms of the Eulerian quantities
$u$ and $w$. Arnold used this calculation to derive the incompressible Euler
equations by seeking $u\in\mathfrak{X}_{\operatorname{vol}}(\Omega)$ (the
subspace of $\mathfrak{X}(\Omega)$ containing divergence-free vector fields),
that extremise the reduced action functional

$S=\int_{0}^{T}\int_{\Omega}\frac{1}{2}|u|^{2}\operatorname{d}x\operatorname{d}t,$
(180)

subject to endpoint conditions $u|_{t=0}=u_{0}$, $u|_{t=T}=u_{T}$ (for some
chosen $u_{0}$ and $u_{T}$), and the boundary conditions $u\cdot n=0$ on the
boundary $\partial\Omega$ of the domain $\Omega$, where $n$ is the outward
pointing normal to $\partial\Omega$. If we perturb $u$ subject to these
conditions, this implies that $w_{t=0}=w_{t=T}=0$.

Taking variations gives

$\displaystyle 0=\delta S$ $\displaystyle=\int_{0}^{T}\int_{\Omega}\delta
u\cdot u\operatorname{d}x\operatorname{d}t,$ (181) $\displaystyle=$
$\displaystyle\int_{0}^{T}\int_{\Omega}(\frac{\partial}{\partial
t}{w}+[u,w])\cdot u\operatorname{d}x\operatorname{d}t,$ (182) $\displaystyle=$
$\displaystyle\int_{0}^{T}\int_{\Omega}w\cdot(-\frac{\partial}{\partial
t}{u}-\nabla\cdot(u\otimes u)-\underbrace{(\nabla
u)^{T}u}_{=\frac{1}{2}|u|^{2}})\operatorname{d}x\operatorname{d}t$
$\displaystyle\quad+\int_{0}^{T}\int_{\partial\Omega}\underbrace{u\cdot
n}_{=0}w\cdot u\operatorname{d}S\operatorname{d}t$
$\displaystyle\qquad+\left[\int_{\Omega}u\cdot\underbrace{w}_{=0}\operatorname{d}x\right]_{t=0}^{t=T},\quad\forall
w\in\mathfrak{X}_{\operatorname{vol}}(\Omega),$ (183)

where $(u\otimes u)_{ij}=u_{i}u_{j}$, and $(\nabla u)_{ij}=\partial
u_{i}/\partial x_{j}$. This formally implies that

$\mathbb{P}\left(\frac{\partial}{\partial t}{u}+u\cdot\nabla u\right)=0,$
(184)

where $\mathbb{P}$ is the $L^{2}$ projection onto divergence-free vector
fields. In other words, there exists a pressure $p$ such that

$\frac{\partial}{\partial t}{u}+u\cdot\nabla u+\nabla p=0,$ (185)

which is the familiar form of the incompressible Euler equation (together with
the divergence free constraint $\nabla\cdot u=0$).

Besides the usual translational and rotational symmetries (assuming
appropriate boundary conditions), the variational formulation has an
additional particle relabelling symmetry, meaning that we can replace $\chi$
by $\chi\circ\psi$ for any diffeomorphism $\psi:\Omega_{0}\to\Omega_{0}$ and
the action remains invariant (in fact, the value of $u$ does not change).
Physically this corresponds to fact that the fluid physics is independent of
the choice of label for a fluid particle (the label for a fluid particle at
time $t$ and position $x_{0}$ being $\chi^{-1}(x)$, _i.e._ $\chi(x_{0},t)=x$).
As discussed in many places (Morrison 1982, Salmon 1998, Shepherd 1990, for
example), this symmetry leads to conservation of circulation

$\frac{\operatorname{d}}{\operatorname{d}t}\int_{C(t)}u\cdot dx=0,$ (186)

for closed loops $C(t)$ that are being transported by the fluid velocity $u$.

The reduction from the flow map $\chi$ to the Eulerian velocity by Arnold was