# An interpretation of TRiSK-type schemes from a discrete exterior calculus
perspective

Christopher Eldred111Center for Computing Research, Sandia National
Laboratories<EMAIL_ADDRESS>Werner Bauer222University of Surrey,
Department of Mathematics<EMAIL_ADDRESS>
(October 2022)

###### Abstract

TRiSK-type numerical schemes are widely used in both atmospheric and oceanic
dynamical cores, due to their discrete analogues of important properties such
as energy conservation and steady geostrophic modes. In this work, we show
that these numerical methods are best understood as a discrete exterior
calculus (DEC) scheme applied to a Hamiltonian formulation of the rotating
shallow water equations based on split exterior calculus. This comprehensive
description of the differential geometric underpinnings of TRiSK-type schemes
completes the one started in [65, 23], and provides a new understanding of
certain operators in TRiSK-type schemes as discrete wedge products and
topological pairings from split exterior calculus. All known TRiSK-type
schemes in the literature are shown to fit inside this general framework, by
identifying the (implicit) choices made for various DEC operators by the
different schemes. In doing so, unexplored choices and combinations are
identified that might offer the possibility of fixing known issues with TRiSK-
type schemes such as operator accuracy and Hollingsworth instability.

###### Contents

  1. 1 Introduction
  2. 2 Split Exterior Calculus for $n=2$
    1. 2.1 Manifolds, Orientations and (Twisted) Forms
    2. 2.2 (Twisted) Hodge Star
    3. 2.3 Vector Proxies
    4. 2.4 Topological Pairing and Topological Functional Derivatives
  3. 3 Continuous RSWE
    1. 3.1 Vector Calculus
      1. 3.1.1 Hamiltonian Form
    2. 3.2 Split Exterior Calculus
      1. 3.2.1 Hamiltonian Form
  4. 4 A Discrete Exterior Calculus in 2D
    1. 4.1 Grids
      1. 4.1.1 Grid Geometry
    2. 4.2 Discrete Differential Forms
    3. 4.3 Operators
      1. 4.3.1 Exterior derivative
      2. 4.3.2 Wedge product
      3. 4.3.3 Hodge star
    4. 4.4 Topological Pairing, Functional Derivatives and Inner Products
    5. 4.5 Operator Adjoints
    6. 4.6 Operator Properties
      1. 4.6.1 Exterior Derivative and Topological Pairing
      2. 4.6.2 Wedge Product
      3. 4.6.3 Hodge Star
  5. 5 Discretization of RSWE with DEC
    1. 5.1 Discrete Hamiltonian and Functional Derivatives
    2. 5.2 Discrete Poisson Brackets
    3. 5.3 Semi-discrete Equations of Motion
    4. 5.4 Linearized Equations
  6. 6 Scheme Properties
    1. 6.1 Topological Properties
      1. 6.1.1 Total Energy Conservation
      2. 6.1.2 Total Mass Conservation
      3. 6.1.3 Total Circulation Conservation
      4. 6.1.4 Total Potential Enstrophy Conservation
      5. 6.1.5 No Spurious Vorticity Production
      6. 6.1.6 PV Compatibility/Steady Geostrophic Modes
    2. 6.2 Metric Properties
    3. 6.3 Summary of Scheme Properties
  7. 7 Specific Operator Choices
    1. 7.1 TRiSK Notation
    2. 7.2 Choices for Hodge Stars: $\operatorname{\mathbf{I}}$, $\operatorname{\mathbf{J}}$ and $\operatorname{\mathbf{H}}$
      1. 7.2.1 Alternative Hodge Stars
    3. 7.3 Choices for PV Wedge Products $\operatorname{\mathbf{R}}$, $\operatorname{\mathbf{W}}$ and $\operatorname{\mathbf{Q}}$
      1. 7.3.1 Choice of $\operatorname{\mathbf{R}}$/$\operatorname{\mathbf{W}}$
      2. 7.3.2 Choice of $\operatorname{\mathbf{Q}}$
    4. 7.4 Choices for KE Wedge Product
    5. 7.5 Choice of grid
    6. 7.6 Grid Optimization
    7. 7.7 Effects of choices on properties
  8. 8 Schemes from the literature
    1. 8.1 Degree of freedom scaling
    2. 8.2 Specific schemes
    3. 8.3 Unexplored approaches
  9. 9 Conclusions
  10. 10 Acknowledgements
  11. A Proof of $\mathbf{D}_{k}$/$\mathbf{\bar{D}}_{k}$/$\left<\left<{},{}\right>\right>$ Properties
  12. B Relationship between PV wedge product Leibniz rules
    1. B.1 Partial Leibniz rule (96)
    2. B.2 Partial Leibniz rule (97)
  13. C Derivation of $\operatorname{\mathcal{H}}$ functional derivatives
    1. C.1 Simplified form for $\frac{\tilde{\delta}\\!\operatorname{\mathcal{H}}}{\delta{v}^{1}}$
  14. D 2D Split Exterior $\leftrightarrow$ Vector Calculus Identities

##  1 Introduction

The rotating shallow water equations (RSWE) are a useful simplified model in
the development process of geophysical fluid dynamics models, as they possess
many of the same properties: conservation laws such as total energy and
potential enstrophy, a (quasi-)balanced state and waves that mediate
departures from this state (inertia-gravity and Rossby waves). However, their
2D nature permits much quicker testing and development of numerical schemes,
since the expense of 3D modelling can be avoided. Additionally, the numerical
methods used in a shallow-water model can often be partially or fully adopted
for use in a 3D model.

A prominent family of numerical methods used for the numerical modelling of
the RSWE are TRiSK-type schemes. These schemes are characterized by a number
of features: Arakawa C-grid staggering333Fluid height at cells centers and
normal component of velocity at cell edges., no spurious vorticity production,
potential vorticity (PV) compatibility/steady geostrophic modes and various
conservation laws (total mass, total circulation, total energy and/or
potential enstrophy). These schemes form the basis for the numerics of both
atmospheric (DYNAMICO [16], LMD-Z [72], MPAS-A [58], Wavetrisk [39], ICON-IAP
[26, 27]) and oceanic (MPAS-O [51]) models. TRiSK-type schemes originated near
the dawn of atmospheric model development, starting with the celebrated
Arakawa-Lamb scheme (AL81) [2] that conserves both total energy and potential
enstrophy. Similar efforts occurred in [54], which also conserves total energy
but only conserves enstrophy in the non-divergent limit, instead of potential
enstrophy. AL81 was extended to incorporate higher-order advection for the
vorticity term in [62]. These schemes are all restricted to uniform
rectangular or lat-lon grids. This shortcoming was addressed in [68, 52],
which generalized the approach of AL81 to arbitrary orthogonal grids444and
gave rise to the name TRiSK, based on the initials of the authors of [52].. A
partial understanding of TRiSK-type schemes in terms of discrete exterior
calculus (DEC)555Discrete exterior calculus [32] is a type of structure-
preserving numerical method for arbitrary polygonal grids, based on the double
deRham complex [41, 30] (implemented through a straight-twisted grid
structure) and split exterior calculus. appeared in [65], which provided the
general framework for further work [75, 77, 76, 67] extending TRiSK-type
schemes to arbitrary non-orthogonal grids and exploring some of the
shortcomings of these schemes (such as operator accuracy). However, this work
also did not recognize the presence of discrete versions of the wedge product
and the topological pairing, which turn out to be key to a complete
understanding of TRiSK-type schemes as a DEC.

In a different direction, the Hamiltonian formulation of the continuous
equations was exploited in [55] to develop a quasi-Hamiltonian discrete model
that generalized both [2] and [62]. This quasi-Hamiltonian formulation was
used to develop a multilayer TRiSK-type shallow water model in [61]. A direct
extension of AL81 to the case of logically square non-orthogonal grids, not
based on either quasi-Hamiltonian formulations or DEC, is found in [73]. An
important step towards a deeper understanding of TRiSK-type schemes occurred
in [19, 23], which unified the quasi-Hamiltonian and DEC viewpoints but still
lacked the recognition of the wedge product and topological pairing.

This paper aims to fill that gap, and will present a complete understanding of
TRiSK-type schemes as a type of DEC applied to a Hamiltonian formulation of
the RSWE based on split exterior calculus, building on prior work in [65, 19,
23]. Split exterior calculus is the extension of standard exterior calculus to
take into account the orientation of the underlying manifold through the
introduction of straight and twisted differential forms, where the latter
change sign under a change of orientation. The use of split exterior calculus
allows one to actually “split” the RSWE into topological and metric parts [4].
When viewed from a split exterior calculus perspective there is a 1-1 mapping
between quantities and operators in the continuous and discrete settings, and
TRiSK-type schemes are quite elegant. This is in contrast to the vector
calculus perspective, where TRiSK-type schemes appear relatively complex, with
many distinct operators.

Despite their popularity, TRiSK-type schemes suffer from several significant
shortcomings. In particular, they have issues with operator accuracy on some
grids [67, 23, 49, 9, 74], especially in a grid refinement context and for
quasi-uniform spherical grids commonly used in atmospheric and oceanic models.
Additionally, TRiSK can suffer from Hollingsworth instability [35, 42, 48, 7]
(also known as Symmetric Instability of the Computational Kind, SICK), a non-
cancellation error in the momentum equations which results in an unphysical
loss of neutral stability of internal inertia-gravity waves. This instability
can cause spurious energy exchanges in balanced motions [17] and errors in the
representation of deep water renewal [59] for ocean models, along with
computational instability in atmospheric models [48, 7]. A general
prescription for avoiding Hollingsworth instability has yet to be found, but
there are several known remedies that can alleviate it [48, 7, 27]. Efforts to
fix the accuracy issues in TRiSK are detailed in [49], but the proposed fixes
lose the desirable properties such as energy conservation and PV
compatibility. Tantalisingly, the scheme in [73] does not suffer from any
accuracy issues despite being utilized on highly distorted non-orthogonal
grids (but still topologically square and periodic). It is hoped that the more
complete understanding of the differential geometric underpinnings of TRiSK-
type schemes presented in this paper will allow the resolution of these issues
in a way that does not break the desirable properties, for arbitrary grid
topologies and geometries. It is anticipated that this will heavily utilize
the separation of the discrete equations into topological parts (associated
with the Poisson bracket) and metric parts (associated with the Hamiltonian),
and is closely related to the ideas expressed in [4, 70, 71]. By categorizing
existing TRiSK-type schemes in terms of the new general DEC, unexplored
operator choices are identified that do not expand the stencil beyond nearest-
neighbor and offer the possibility of achieving better accuracy or alleviating
Hollingsworth instability; while still keeping the desirable properties.
Similar efforts in the context of split finite element methods rather than DEC
are detailed in [5, 6].

The rest of this paper is structured as follows: Section 2 provides a brief
review of split exterior calculus, focusing on new aspects compared to
standard exterior calculus; Section 3 reviews the continuous RSWE in
Hamiltonian form, using both vector calculus and split exterior calculus;
Section 4 introduces a discrete exterior calculus (DEC) in 2D; Section 5
discretizes the RSWE using this DEC; Section 6 discusses properties of the
scheme and it’s dependence on operator properties from DEC; Section 7
identifies specific DEC operator choices with the required properties; Section
8 explores the relationship between TRiSK-type schemes in the literature and
the DEC-based scheme presented here by identifying the choices these schemes
(implicitly) made for various operators in the DEC context, and includes a
discussion in Section 8.3 of unexplored possibilities in choices for these
operators; and finally Section 9 offers some conclusions. Appendix A proves
some useful properties of the discrete exterior derivative and topological
pairing, Appendix B derives some relationships between various desired
properties of the PV wedge product, Appendix C derives the functional
derivatives of the discrete Hamiltonian of the RSWE, and Appendix D explores
the relationships between quantities and operators for split exterior calculus
and vector calculus in 2D.

##  2 Split Exterior Calculus for $n=2$

Split exterior calculus [4] is standard exterior calculus extended with an
explicit notion of orientation for the underlying manifold. It is assumed that
the reader has a basic knowledge of standard exterior calculus and
differential geometry, for reference see standard textbooks on differential
geometry [25, 1] or the review paper at [4]. In this section only the new
concepts specific to split exterior calculus are introduced: straight and
twisted differential forms, the twisted Hodge star and the topological
pairing.

###  2.1 Manifolds, Orientations and (Twisted) Forms

Consider a Riemannian manifold $\mathbb{M}$ with boundary $\delta\mathbb{M}$
and choice of ambient orientation $Or$, on which the space of straight
$k$-forms ${\omega}^{k}$ is denoted with $\Lambda^{k}$. The opposite
orientation is denoted with $-Or$. A classical example of orientation on a
manifold is the choice between a right-hand or left-hand rule for the cross-
product in $\mathbb{R}^{3}$.

Then the space of twisted differential $k$-forms ${\tilde{\omega}}^{k}$ is
defined similarly to straight differential forms, except that twisted
differential forms change sign under a change of orientation from $Or$ to
$-Or$, while straight forms do not.666More concretely, twisted differential
forms are a type of (vector)-bundle valued differential forms (BVDFs), where
the bundle is the pseudoscalar bundle $\Psi$. In contrast, straight
differential forms are BVDFs where the bundle is $\mathbb{R}$. Since $\Psi$
has the canonical trivializations $\Psi\otimes\Psi\rightarrow\mathbb{R}$ and
$\Psi\otimes\mathbb{R}\rightarrow\Psi$, the wedge product of two twisted forms
produces a straight form and the wedge product of a straight form and a
twisted form produces a twisted form. This change of sign is analogous to the
situation with pseudovectors and vectors or pseudoscalars and scalars in
vector calculus.

Using both straight and twisted differential forms to describe a physical
system allows a formulation that is manifestly independent of the choice of
ambient orientation [4]. Since the choice of orientation is (for classical
physics at least) arbitrary, this is a desirable feature. These ideas are
closely related to the work of Tonti [71, 70], who characterized a wide range
of variables for physical theories in terms of their association with oriented
geometric entities: source (twisted differential forms) and configuration
(straight differential forms).

###  2.2 (Twisted) Hodge Star

Rather than the standard orientation-dependent Hodge star $\star$, an
orientation independent Hodge star ($\operatorname{\tilde{\star}}$, the
twisted Hodge star) is defined through

${a}^{k}\wedge\operatorname{\tilde{\star}}{b}^{k}=<{a}^{k},{b}^{k}>\tilde{\mu}^{n}$
(1)

where $<>$ is the pointwise $L_{2}$ inner product that produces a straight
$0$-form and $\tilde{\mu}^{n}$ is the twisted volume form, with a similar
equivalent definition for $\tilde{a}^{k}$ and $\tilde{b}^{k}$. The definition
(1) is very similar to the one used for the standard Hodge star, but with the
twisted volume form $\tilde{\mu}^{n}$ instead of the straight volume form. It
is clear from this definition that the twisted Hodge star
$\operatorname{\tilde{\star}}$ must map between straight $k$-forms and twisted
$(n-k)$-forms or twisted $k$-forms and straight $(n-k)$-forms, since
${a}^{k}\wedge\operatorname{\tilde{\star}}{b}^{k}$ must produce a twisted
$n$-form. Additionally, it has the property
$\operatorname{\tilde{\star}}\operatorname{\tilde{\star}}=(-1)^{k(n-k)}$, just
like $\star$.

###  2.3 Vector Proxies

Consider a vector $\operatorname{\mathbf{v}}$ and a pseudovector
$\tilde{\operatorname{\mathbf{v}}}$. Using the flat operator $\flat$, the
interior product (contraction) $\mathrm{i}$ and the twisted volume form
$\tilde{\mu}^{n}$, there are four vector proxies that can be constructed from
$\operatorname{\mathbf{v}}$ and $\tilde{\operatorname{\mathbf{v}}}$:

$\displaystyle{v}^{1}$ $\displaystyle=$
$\displaystyle\operatorname{\mathbf{v}}^{\flat},$ (2)
$\displaystyle\tilde{v}^{1}$ $\displaystyle=$
$\displaystyle\tilde{\operatorname{\mathbf{v}}}^{\flat},$ (3)
$\displaystyle{v}^{n-1}$ $\displaystyle=$
$\displaystyle\mathrm{i}_{\tilde{\operatorname{\mathbf{v}}}}\tilde{\mu}^{n}=\operatorname{\tilde{\star}}\tilde{v}^{1},$
(4) $\displaystyle\tilde{v}^{n-1}$ $\displaystyle=$
$\displaystyle\mathrm{i}_{\operatorname{\mathbf{v}}}\tilde{dV}^{n}=\operatorname{\tilde{\star}}{v}^{1},$
(5)

where in the last two equations Hirani’s formula
$\mathrm{i}_{\mathbf{x}}\alpha=(-1)^{k(n-k)}\operatorname{\tilde{\star}}(\operatorname{\tilde{\star}}\alpha\wedge\mathbf{x}^{\flat})$
[32] for an arbitrary $k$-form $\alpha$ has been used. From these proxies it
is clear that $\operatorname{\tilde{\star}}$ transforms (straight) twisted
circulations ($1$-forms) into (twisted) straight fluxes ($(n-1)$-forms), or
vice versa (sometimes with a minus sign since
$\operatorname{\tilde{\star}}\operatorname{\tilde{\star}}=(-1)^{k(n-k)}$). A
direct transformation of circulations into circulations or fluxes into fluxes
(i.e. between straight and twisted $k$-forms) is accomplished through the use
of the wedge product $\tilde{I}^{0}\wedge$, for example
${x}^{n-1}=\tilde{I}^{0}\wedge\tilde{x}^{n-1}$.

###### Remark 2.1

Note here the distinction between $1$-forms (circulations) and $(n-1)$-forms
(fluxes). In 2D ($n=2$) they appear to be the same object since $n-1=1$.
However, this obscures their true nature, and leads to convoluted notation
(see for example [12]) that attempts to distinguish between them. In this work
we will keep a clear distinction between the two.

###  2.4 Topological Pairing and Topological Functional Derivatives

Here we introduce new concepts that extend the existing split exterior
calculus as introduced in [4], namely the topological pairing and a
corresponding topological functional derivative, both crucial for this work.
The first new concept in split exterior calculus is the _topological pairing_
, which is defined by:

$\left<\left<{{x}^{k}},{\tilde{y}^{n-k}}\right>\right>:=\int_{\mathbb{M}}{x}^{k}\wedge\tilde{y}^{n-k}=(-1)^{k(n-k)}\int_{\mathbb{M}}\tilde{y}^{n-k}\wedge{x}^{k}=(-1)^{k(n-k)}\left<\left<{\tilde{y}^{n-k}},{{x}^{k}}\right>\right>.$
(6)

This pairing is purely topological, since both the wedge product and
integration of a differential form do not require reference to a metric. The
topological pairing can be connected to the inner product (which does require
a metric, for the Hodge star $\operatorname{\tilde{\star}}$) through:

$\displaystyle\left<\left<{{a}^{k}},{\tilde{b}^{n-k}}\right>\right>$
$\displaystyle=$ $\displaystyle\left<{{a}^{k}},{{b}^{k}}\right>,$ (7)
$\displaystyle\left<\left<{\tilde{a}^{k}},{{b}^{n-k}}\right>\right>$
$\displaystyle=$
$\displaystyle(-1)^{k(n-k)}\left<{\tilde{a}^{k}},{\tilde{b}^{k}}\right>,$ (8)

recalling the inner product of two straight $k$-forms is defined as
$\left<{{a}^{k}},{{b}^{k}}\right>:=\int_{\mathbb{M}}{a}^{k}\wedge\tilde{\star}{b}^{k}$
(a definition that also holds for twisted forms).

Additionally, an integration by parts rule holds:

$\displaystyle\left<\left<{\operatorname{d}{x}^{k-1}},{\tilde{y}^{n-k}}\right>\right>+(-1)^{k-1}\left<\left<{{x}^{k-1}},{\operatorname{d}\tilde{y}^{n-k}}\right>\right>=\int_{\delta\mathbb{M}}tr({x}^{k-1}\wedge\tilde{y}^{n-k}),$
(9)
$\displaystyle\left<\left<{\operatorname{d}\tilde{x}^{k-1}},{{y}^{n-k}}\right>\right>+(-1)^{k-1}\left<\left<{\tilde{x}^{k-1}},{\operatorname{d}{y}^{n-k}}\right>\right>=\int_{\delta\mathbb{M}}tr(\tilde{x}^{k-1}\wedge{y}^{n-k}),$
(10)

where $tr$ is the trace operator. In the case of a domain without boundaries
(such as the ones we consider here), the boundary term on the right-hand side
drops out.

Secondly, we define topological functional derivatives relative to this
pairing, through

###### Definition 2.2

The _topological_ functional derivatives of
$\operatorname{\mathcal{F}}[{x}^{k},\tilde{y}^{k}]:\Lambda^{k}\times\tilde{\Lambda}^{k}\rightarrow\mathbb{R}$
with respect to the straight $k$-form ${x}^{k}$ and twisted $k$-form
$\tilde{x}^{k}$ and with respect to the topological pairing
$\left<\left<{},{}\right>\right>$ are defined by

$\begin{split}&\delta\operatorname{\mathcal{F}}:=\lim_{\epsilon\rightarrow
0}\frac{1}{\epsilon}\big{(}\operatorname{\mathcal{F}}[{x}^{k}+\epsilon\,{\omega}^{k}]-\operatorname{\mathcal{F}}[{x}^{k}]\big{)}=:\left<\left<{{\omega}^{k}},{\frac{\tilde{\delta}\\!\operatorname{\mathcal{F}}}{\delta{x}^{k}}}\right>\right>=\int_{\mathcal{M}}{\omega}^{k}\wedge\frac{\tilde{\delta}\\!\operatorname{\mathcal{F}}}{\delta{x}^{k}}\quad\forall{\omega}^{k}\in\Lambda^{k},\\\
\text{resp.}\ &\delta\operatorname{\mathcal{F}}:=\lim_{\epsilon\rightarrow
0}\frac{1}{\epsilon}\big{(}\operatorname{\mathcal{F}}[\tilde{y}^{k}+\epsilon\,\tilde{\omega}^{k}]-\operatorname{\mathcal{F}}[\tilde{y}^{k}]\big{)}=:\left<\left<{\tilde{\omega}^{k}},{\frac{\tilde{\delta}\\!\operatorname{\mathcal{F}}}{\delta\tilde{y}^{k}}}\right>\right>=\int_{\mathcal{M}}\tilde{\omega}^{k}\wedge\frac{\tilde{\delta}\\!\operatorname{\mathcal{F}}}{\delta\tilde{y}^{k}}\quad\forall\tilde{\omega}^{k}\in\tilde{\Lambda}^{k},\end{split}$
(11)

for arbitrary test functions ${\omega}^{k}$ (resp. $\tilde{\omega}^{k}$). In
particular for ${x}^{k}\in\Lambda^{k}$ the topological functional derivative
$\frac{\tilde{\delta}\\!\operatorname{\mathcal{F}}}{\delta{x}^{k}}\in\tilde{\Lambda}^{n-k}$
is a twisted $(n-k)$-form, while for $\tilde{y}^{k}\in\tilde{\Lambda}^{k}$ the
topological functional derivative
$\frac{\tilde{\delta}\\!\operatorname{\mathcal{F}}}{\delta\tilde{y}^{k}}\in\Lambda^{(n-k)}$
is a straight $(n-k)$-form.

Note that counting the tildes reveals easily the nature of the functional
derivatives: an odd number of tildes give twisted differential forms, while an
even number gives straight differential forms.

##  3 Continuous RSWE

We will now write the rotating shallow water equations (RSWE) in Hamiltonian
form using both vector calculus and split exterior calculus. The former is
mostly a review of material from [55, 19]. The two subsections are structured
identically to help facilitate understanding, and enable quick comparison
between these two approaches. These equations apply to an arbitrary 2D
manifold without boundaries (for example, the sphere or doubly periodic
plane), taking into account the fact that the operators can all be defined
intrinsically and the vector fields are tangent to the manifold. The
coordinate system is assumed to be an arbitrary rotating Eulerian coordinate
system, with rotational velocity $\operatorname{\mathbf{R}}$ and Coriolis
parameter $f=\nabla^{\perp}\cdot\operatorname{\mathbf{R}}$. This rotation is
almost always chosen to correspond with the rotation of the underlying
geophysical system (for example, the rotation of the Earth).

In the vector calculus form, the physical quantities are written in terms of
scalars ($a$ and $b$) and vectors ($\mathbf{x}$ and $\mathbf{y}$), and the key
operators are: the differential operators gradient $\nabla a$, skew-gradient
$\nabla^{\perp}a$, divergence $\nabla\cdot\mathbf{x}$ and curl
$\nabla^{\perp}\cdot\mathbf{x}$; the ”product” operators vector transpose
$\mathbf{x}^{\perp}$, vector dot product $\mathbf{x}\cdot\mathbf{y}$, scalar
product $ab$ and mixed product $a\mathbf{y}$ or $a\mathbf{y}^{\perp}$; and the
inner product $\left<{},{}\right>$, which is defined as
$\left<{a},{b}\right>=\int_{\Omega}ab$ for scalars and
$\left<{\mathbf{x}},{\mathbf{y}}\right>=\int_{\Omega}\mathbf{x}\cdot\mathbf{y}$
for vectors. Functional derivatives are defined relative to the appropriate
inner product.

###### Remark 3.1

For a 2D manifold embedded in $\mathbb{R}^{3}$ with unit normal vector
$\operatorname{\mathbf{\hat{k}}}$ (in $\mathbb{R}^{3}$), if we consider the
vector $\mathbf{x}$ to lie in $\mathbb{R}^{3}$ (recalling it is still tangent
to the manifold), some of these operators take simpler forms in terms of
common operators (gradient $\nabla_{3}$ and curl $\nabla_{3}\times$, along
with cross product $\times$ and dot product $\cdot$) from vector calculus in
$\mathbb{R}^{3}$:
$\mathbf{x}^{\perp}=\operatorname{\mathbf{\hat{k}}}\times\mathbf{x}$,
$\nabla^{\perp}x=\operatorname{\mathbf{\hat{k}}}\times\nabla_{3}x$,
$\nabla^{\perp}\cdot\mathbf{x}=\operatorname{\mathbf{\hat{k}}}\cdot\nabla_{3}\times\mathbf{x}$.
However, all of the 2D operators can also be defined intrinsically without
requiring that the manifold is embedded in $\mathbb{R}^{3}$ and without
reference to a $\operatorname{\mathbf{\hat{k}}}$.

In the split exterior calculus form the physical quantities are written in
terms of straight $k$-forms ${x}^{k}$ and twisted $k$-forms $\tilde{x}^{k}$;
and the key operators are the exterior derivative $\operatorname{d}$, the
wedge product $\wedge$, the twisted Hodge star $\operatorname{\tilde{\star}}$,
the metric inner product $\left<{},{}\right>$ and the topological pairing
$\left<\left<{},{}\right>\right>$. Compared to vector calculus, the number of
operators is greatly reduced. As discussed above, for split exterior calculus
we use the topological functional derivative defined relative to the
topological pairing, since the topological functional derivative gives a
natural splitting between topological parts of the equations (the Poisson
brackets) and the metric parts of the equations (the Hamiltonian and it’s
functional derivatives). This splitting is explored in more detail below.

See Appendix D for identities relating the quantities and operators of vector
calculus and split exterior calculus in 2D. In particular, scalars are either
$0$-forms or $2$-forms and vectors are either $1$-forms or $(n-1)$-forms
(recalling that we are in 2D and therefore $n-1=1$), the differential
operators are composed of the exterior derivative $\operatorname{d}$ plus the
Hodge star $\operatorname{\tilde{\star}}$, the product operators are composed
of the wedge product $\wedge$ plus the Hodge star
$\operatorname{\tilde{\star}}$, and the inner products for scalars and vectors
become either the inner product for differential forms or the topological
pairing.

###  3.1 Vector Calculus

Consider the RSWE with fluid height $h$, relative velocity
$\operatorname{\mathbf{u}}$, absolute velocity
$\operatorname{\mathbf{v}}=\operatorname{\mathbf{u}}+\operatorname{\mathbf{R}}$,
bottom topography $h_{s}$, Coriolis parameter $f$, relative vorticity
$\zeta=\nabla^{\perp}\cdot\operatorname{\mathbf{u}}$, absolute vorticity
$\eta=\nabla^{\perp}\cdot\operatorname{\mathbf{v}}=\zeta+f$, and potential
vorticity $q$ (defined through $hq=\eta$). The equations of motion for $h$ and
$\operatorname{\mathbf{v}}$ are

$\displaystyle\frac{\partial\operatorname{\mathbf{v}}}{\partial
t}+f\operatorname{\mathbf{u}}^{\perp}+\zeta\operatorname{\mathbf{u}}^{\perp}+\nabla(gh+gh_{s}+\frac{\operatorname{\mathbf{u}}\cdot\operatorname{\mathbf{u}}}{2})$
$\displaystyle=$ $\displaystyle 0,$ (12) $\displaystyle\frac{\partial
h}{\partial t}+\nabla\cdot(h\operatorname{\mathbf{u}})$ $\displaystyle=$
$\displaystyle 0.$ (13)

Their linearized form (around $h=H$, $\operatorname{\mathbf{u}}=0$) is given
by

$\displaystyle\frac{\partial\operatorname{\mathbf{v}}}{\partial
t}+f\operatorname{\mathbf{u}}^{\perp}+\nabla(gh+gh_{s})$ $\displaystyle=$
$\displaystyle 0,$ (14) $\displaystyle\frac{\partial h}{\partial
t}+\nabla\cdot(H\operatorname{\mathbf{u}})$ $\displaystyle=$ $\displaystyle
0,$ (15)

where primes have been dropped for convenience. The linear equations will be
useful when discretizing, to understand the behaviour of the scheme for steady
geostrophic modes.

####  3.1.1 Hamiltonian Form

These equations can be put into Hamiltonian form (see [56]) using the
Hamiltonian $\operatorname{\mathcal{H}}[\operatorname{\mathbf{v}},h]$:

$\mathcal{H}[\operatorname{\mathbf{v}},h]=\frac{g}{2}\left<{h},{h}\right>+g\left<{h},{h_{s}}\right>+\left<{h},{\frac{\operatorname{\mathbf{u}}\cdot\operatorname{\mathbf{u}}}{2}}\right>,$
(16)

and Poisson brackets
$\\{\operatorname{\mathcal{A}},\operatorname{\mathcal{B}}\\}$ for arbitrary
functionals $\operatorname{\mathcal{A}},\operatorname{\mathcal{B}}$:

$\\{\operatorname{\mathcal{A}},\operatorname{\mathcal{B}}\\}=-\left<\frac{\delta\\!\operatorname{\mathcal{A}}}{\delta\operatorname{\mathbf{v}}},q\frac{\delta\\!\operatorname{\mathcal{B}}}{\delta\operatorname{\mathbf{v}}}^{\perp}\right>-\left<\frac{\delta\\!\operatorname{\mathcal{A}}}{\delta\operatorname{\mathbf{v}}},\nabla\frac{\delta\\!\operatorname{\mathcal{B}}}{\delta
h}\right>-\left<\frac{\delta\\!\operatorname{\mathcal{A}}}{\delta
h},\nabla\cdot\frac{\delta\\!\operatorname{\mathcal{B}}}{\delta\operatorname{\mathbf{v}}}\right>,$
(17)

recalling that $q$ is defined through
$hq=\nabla^{\perp}\cdot\operatorname{\mathbf{v}}=\eta$. The last term in
$\operatorname{\mathcal{H}}[\operatorname{\mathbf{v}},h]$ (the kinetic energy)
can also be written as
$\frac{1}{2}\left<{h\operatorname{\mathbf{u}}},{\operatorname{\mathbf{u}}}\right>$.
The functional derivatives of
$\operatorname{\mathcal{H}}[\operatorname{\mathbf{v}},h]$ are

$\displaystyle\operatorname{\mathbf{F}}:=\frac{\delta\\!\operatorname{\mathcal{H}}}{\delta\operatorname{\mathbf{v}}}=h\operatorname{\mathbf{u}},\quad\quad\quad
B:=\frac{\delta\\!\operatorname{\mathcal{H}}}{\delta
h}=gh+gh_{s}+\frac{\operatorname{\mathbf{u}}\cdot\operatorname{\mathbf{u}}}{2}.$
(18)

Inserting the functional derivatives (18) into the Poisson brackets (17) and
using $\operatorname{\mathcal{F}}=\left<\hat{x},x\right>$ for
$x\in(h,\operatorname{\mathbf{v}})$ with arbitrary test function
$\hat{x}\in(\hat{h},\hat{\operatorname{\mathbf{v}}})$ gives the equations of
motion:

$\displaystyle\frac{\partial\operatorname{\mathbf{v}}}{\partial
t}+q\operatorname{\mathbf{F}}^{\perp}+\nabla B$ $\displaystyle=$
$\displaystyle 0,$ (19) $\displaystyle\frac{\partial h}{\partial
t}+\nabla\cdot\operatorname{\mathbf{F}}$ $\displaystyle=$ $\displaystyle 0,$
(20)

which upon substitution of the specific values of the functional derivatives
reduce to (12) - (13).

##### Casimirs

The general expression for the RSWE Casimirs
$\operatorname{\mathcal{C}}[\operatorname{\mathbf{v}},h]$ is given by

$\operatorname{\mathcal{C}}[\operatorname{\mathbf{v}},h]=\int hf(q),$ (21)

with arbitrary function $f(q)$, which has functional derivatives

$\displaystyle\frac{\delta\\!\operatorname{\mathcal{C}}}{\delta\operatorname{\mathbf{v}}}=-\nabla^{\perp}f^{\prime}(q),\quad\quad\quad\frac{\delta\\!\operatorname{\mathcal{C}}}{\delta
h}=f(q)-qf^{\prime}(q),$ (22)

where we have used $q=\frac{\nabla^{\perp}\cdot\operatorname{\mathbf{v}}}{h}$
and also integrated by parts (and dropped boundary terms since we are in a
domain without boundaries). Insertion of (22) into the Poisson brackets (17)
confirms (after some algebra) that

$\\{\operatorname{\mathcal{C}},\operatorname{\mathcal{A}}\\}=0$ (23)

for any $\operatorname{\mathcal{A}}$. Important examples of RSWE Casimirs are
total mass ($f=1$), total mass-weighted potential vorticity (total
circulation) ($f=q$) and total potential enstrophy ($f=\frac{q^{2}}{2}$).

##### Vorticity Dynamics

The evolution equation for absolute vorticity $\eta$ is obtained by taking
$\nabla^{\perp}\cdot$ of (19) yielding

$\frac{\partial\eta}{\partial t}+\nabla\cdot(q\operatorname{\mathbf{F}})=0.$
(24)

The last term in (24) can also be written as
$\nabla^{\perp}\cdot(q\operatorname{\mathbf{F}}^{\perp})$. Using (24) and (20)
we can get an evolution equation for $q$ as

$\frac{\partial q}{\partial t}+\frac{\operatorname{\mathbf{F}}}{h}\cdot\nabla
q=0.$ (25)

The last term in (25) can also be written as
$\frac{\operatorname{\mathbf{F}}^{\perp}}{h}\cdot\nabla^{\perp}q$. From (25)
we see that $q$ is materially conserved, since
$\frac{\operatorname{\mathbf{F}}}{h}=\operatorname{\mathbf{u}}$ is the
material velocity. Vorticity dynamics will be useful for understanding the PV
compatibility behaviour of the discrete scheme.

##### Linearized Equations

Following standard procedures [57] the Hamiltonian form can be linearized
around $h=H$, $\operatorname{\mathbf{u}}=0$ to yield for arbitrary
$\operatorname{\mathcal{A}},\operatorname{\mathcal{B}}$

$\\{\operatorname{\mathcal{A}},\operatorname{\mathcal{B}}\\}_{lin}=-\left<\frac{\delta\\!\operatorname{\mathcal{A}}}{\delta\operatorname{\mathbf{v}}},\frac{f}{H}\frac{\delta\\!\operatorname{\mathcal{B}}}{\delta\operatorname{\mathbf{v}}}^{\perp}\right>-\left<\frac{\delta\\!\operatorname{\mathcal{A}}}{\delta\operatorname{\mathbf{v}}},\nabla\frac{\delta\\!\operatorname{\mathcal{B}}}{\delta
h}\right>-\left<\frac{\delta\\!\operatorname{\mathcal{A}}}{\delta
h},\nabla\cdot\frac{\delta\\!\operatorname{\mathcal{B}}}{\delta\operatorname{\mathbf{v}}}\right>,$
(26)

and

$\mathcal{H}_{lin}[\operatorname{\mathbf{v}},h]=\int
g\frac{h^{2}}{2}+ghh_{s}+H\frac{\operatorname{\mathbf{u}}\cdot\operatorname{\mathbf{u}}}{2}.$
(27)

The functional derivatives are thus

$\displaystyle\operatorname{\mathbf{F}}_{lin}:=\frac{\delta\\!\operatorname{\mathcal{H}}_{lin}}{\delta\operatorname{\mathbf{v}}}=H\operatorname{\mathbf{u}},\quad\quad\quad
B_{lin}:=\frac{\delta\\!\operatorname{\mathcal{H}}_{lin}}{\delta
h}=gh+gh_{s},$ (28)

with equations of motion

$\displaystyle\frac{\partial\operatorname{\mathbf{v}}}{\partial
t}+\frac{f}{H}\operatorname{\mathbf{F}}_{lin}^{\perp}+\nabla B_{lin}$
$\displaystyle=$ $\displaystyle 0,$ (29) $\displaystyle\frac{\partial
h}{\partial t}+\nabla\cdot\operatorname{\mathbf{F}}_{lin}$ $\displaystyle=$
$\displaystyle 0.$ (30)

Insertion of (28) into (29) - (30) gives (14) - (15). The linearized dynamics
will be useful for understanding the linear modes behaviour of the discrete
scheme.

###  3.2 Split Exterior Calculus

We will now present the same development using split exterior calculus instead
of vector calculus. The two sections are structured the same, and it is highly
instructive to compare them.

Consider the RSWE with fluid height twisted 2-form $\tilde{h}^{2}$, relative
velocity straight 1-form ${u}^{1}$, rotational velocity straight 1-form
${R}^{1}$, absolute velocity straight 1-form ${v}^{1}={u}^{1}+{R}^{1}$, bottom
topography twisted 2-form $\tilde{h_{s}}^{2}$, Coriolis parameter straight
2-form ${f}^{2}=\operatorname{d}{R}^{1}$, relative vorticity straight 2-form
${\zeta}^{2}=\operatorname{d}{u}^{1}$, absolute vorticity straight 2-form
${\eta}^{2}=\operatorname{d}{v}^{1}={\zeta}^{2}+{f}^{2}$, and potential
vorticity twisted 0-form $\tilde{q}^{0}$ (defined through
$\tilde{h}^{2}\wedge\tilde{q}^{0}={\eta}^{2}$). The correspondence between
these objects and their vector calculus analogues is given in Table 1. These
choices of types and degrees of forms for the various physical quantities are
motivated by the work of Tonti and collaborators [71, 70], who conclusively
demonstrated the association of physical quantities with oriented geometric
entities, which are nothing more than straight and twisted differential forms.
The same choices were made in [4, 23, 65].

Vector Calculus | Split Exterior Calculus
---|---
$h$ |  $\tilde{h}^{2}$ or ${h}^{0}$
$\operatorname{\mathbf{v}}$ | ${v}^{1}$
$\operatorname{\mathbf{u}}$ |  $\tilde{u}^{n-1}$ or ${u}^{1}$
$\operatorname{\mathbf{R}}$ | ${R}^{1}$
$h_{s}$ |  $\tilde{h_{s}}^{2}$ or ${h_{s}}^{0}$
$f$ |  ${f}^{2}$ or $\tilde{f}^{0}$
$\eta$ |  ${\eta}^{2}$ or $\tilde{\eta}^{0}$
$\zeta$ |  ${\zeta}^{2}$ or $\tilde{\zeta}^{0}$
$q$ | $\tilde{q}^{0}$
$\operatorname{\mathbf{F}}$ | $\tilde{F}^{n-1}$
$B$ | ${B}^{0}$
Table 1: The relationships between quantities in vector calculus and split
exterior calculus forms of the RSWE. Note that many of the scalar quantities
such as $h$ appear as both $0$-forms and $2$-forms. This corresponds to the
fact that a standard vector calculus representation does not take into account
whether a ”scalar” is actually a density, and also whether it is a normal or
pseudo quantity. A similar thing happens with $\operatorname{\mathbf{u}}$,
which can be represented as either a circulation straight $1$-form or a flux
twisted $(n-1)$-form.

Introduce also the auxiliary quantities fluid height straight 0-form
${h}^{0}=\operatorname{\tilde{\star}}\tilde{h}^{2}$, bottom topography
straight 0-form ${h_{s}}^{0}=\operatorname{\tilde{\star}}\tilde{h_{s}}^{2}$,
relative velocity twisted $(n\\!-\\!1)$-form
$\tilde{u}^{n-1}=\operatorname{\tilde{\star}}{u}^{1}$, relative vorticity
twisted 0-form $\tilde{\zeta}^{0}=\operatorname{\tilde{\star}}{\zeta}^{2}$ and
Coriolis parameter twisted 0-form
$\tilde{f}^{0}=\operatorname{\tilde{\star}}{f}^{2}$, which are related to the
corresponding primary quantities through the Hodge star
$\operatorname{\tilde{\star}}$.

The equations of motion for $\tilde{h}^{2}$ and ${v}^{1}$ are given by:

$\displaystyle\frac{\partial{v}^{1}}{\partial
t}+\tilde{f}^{0}\wedge\tilde{u}^{n-1}+\tilde{\zeta}^{0}\wedge\tilde{u}^{n-1}+\operatorname{d}(g{h}^{0}+g{h_{s}}^{0}+\operatorname{\tilde{\star}}\frac{({u}^{1}\wedge\tilde{u}^{n-1})}{2})=0,$
(31) $\displaystyle\frac{\partial\tilde{h}^{2}}{\partial
t}+\operatorname{d}({h}^{0}\wedge\tilde{u}^{n-1})=0.$ (32)

Their linearized form (around ${h}^{0}=H$ and ${u}^{1}=0$) is

$\displaystyle\frac{\partial{v}^{1}}{\partial
t}+\tilde{f}^{0}\wedge\tilde{u}^{n-1}+\operatorname{d}(g{h}^{0}+g{h_{s}}^{0})=0,$
(33) $\displaystyle\frac{\partial\tilde{h}^{2}}{\partial
t}+\operatorname{d}(H\tilde{u}^{n-1})=0,$ (34)

where primes have been dropped for convenience.

The equations (31) - (34) can be obtained by applying the relationships in
Appendix D to (12) - (15).777A more fundamental derivation [36] is possible by
starting with an Euler-Poincaré variational formulation based on semi-direct
product theory and written in terms of split exterior calculus; and then
applying a Legendre transform to get a Lie-Poisson Hamiltonian formulation,
and finally making a change of variables from momentum to velocity to get a
curl-form Hamiltonian formulation as used here.

####  3.2.1 Hamiltonian Form

Equations (31) - (32) can be put into Hamiltonian form using the Hamiltonian
$\operatorname{\mathcal{H}}[{v}^{1},\tilde{h}^{2}]$:

$\operatorname{\mathcal{H}}[{v}^{1},\tilde{h}^{2}]=\frac{g}{2}\left<{\tilde{h}^{2}},{\tilde{h}^{2}}\right>+g\left<{\tilde{h}^{2}},{\tilde{h_{s}}^{2}}\right>+\left<{\tilde{h}^{2}},{\frac{{u}^{1}\wedge\tilde{u}^{n-1}}{2}}\right>,$
(35)

(the ordering ${u}^{1}\wedge\tilde{u}^{n-1}$ is important due to anti-symmetry
of the wedge product) and Poisson bracket for arbitrary
$\operatorname{\mathcal{A}},\operatorname{\mathcal{B}}$:

$\\{\operatorname{\mathcal{A}},\operatorname{\mathcal{B}}\\}=-\left<\left<{\frac{\tilde{\delta}\\!\operatorname{\mathcal{A}}}{\delta{v}^{1}}},{\tilde{q}^{0}\wedge\frac{\tilde{\delta}\\!\operatorname{\mathcal{B}}}{\delta{v}^{1}}}\right>\right>-\left<\left<{\frac{\tilde{\delta}\\!\operatorname{\mathcal{A}}}{\delta{v}^{1}}},{\operatorname{d}\frac{\tilde{\delta}\\!\operatorname{\mathcal{B}}}{\delta\tilde{h}^{2}}}\right>\right>-\left<\left<{\frac{\tilde{\delta}\\!\operatorname{\mathcal{A}}}{\delta\tilde{h}^{2}}},{\operatorname{d}\frac{\tilde{\delta}\\!\operatorname{\mathcal{B}}}{\delta{v}^{1}}}\right>\right>,$
(36)

where $\left<\left<{},{}\right>\right>$ is the topological pairing defined in
Eqn. (6). The last term in $\operatorname{\mathcal{H}}[{v}^{1},\tilde{h}^{2}]$
can also be written as
$\frac{1}{2}\left<{{u}^{1}},{{h}^{0}\wedge{u}^{1}}\right>=\frac{1}{2}\left<{\tilde{u}^{n-1}},{{h}^{0}\wedge\tilde{u}^{n-1}}\right>$
and we recall that ${v}^{1}={u}^{1}+{R}^{1}$.

A key point is that the Poisson brackets are purely topological: they involve
only the wedge product, exterior derivative and topological pairing, all of
which can be defined without using a metric. All of the metric information is
present in the Hamiltonian, through the Hodge star that is part of the
definition of the inner product.

The functional derivatives of
$\operatorname{\mathcal{H}}[{v}^{1},\tilde{h}^{2}]$ are

$\tilde{F}^{n-1}:=\frac{\tilde{\delta}\\!\operatorname{\mathcal{H}}}{\delta{v}^{1}}=\frac{1}{2}\left[{h}^{0}\wedge\tilde{u}^{n-1}+\operatorname{\tilde{\star}}({h}^{0}\wedge{u}^{1})\right],\quad\quad{B}^{0}:=\frac{\tilde{\delta}\\!\operatorname{\mathcal{H}}}{\delta\tilde{h}^{2}}=g{h}^{0}+g{h_{s}}^{0}+\operatorname{\tilde{\star}}\frac{({u}^{1}\wedge\tilde{u}^{n-1})}{2}.$
(37)

In taking the first functional derivative
($\frac{\tilde{\delta}\\!\operatorname{\mathcal{H}}}{\delta{v}^{1}}$), we have
used the linearity of the Hodge star to write
$\delta\tilde{u}^{n-1}=\delta\operatorname{\tilde{\star}}{u}^{1}=\operatorname{\tilde{\star}}\delta{u}^{1}$.
In the continuous case, it possible to simplify $\tilde{F}^{n-1}$ as
${h}^{0}\wedge\tilde{u}^{n-1}$. However, this simplification will be valid
only for certain choices of operators in the discrete case.

Inserting the functional derivatives (37) into the Poisson brackets (36) and
using $\operatorname{\mathcal{F}}=\left<\hat{x},x\right>$ for
$x\in(\tilde{h}^{2},{v}^{1})$ with arbitrary test function
$\hat{x}\in(\widehat{\tilde{h}^{2}},\widehat{{v}^{1}})$ gives the equations of
motion:

$\displaystyle\frac{\partial{v}^{1}}{\partial
t}+\tilde{q}^{0}\wedge\tilde{F}^{n-1}+\operatorname{d}{B}^{0}=0,$ (38)
$\displaystyle\frac{\partial\tilde{h}^{2}}{\partial
t}+\operatorname{d}\tilde{F}^{n-1}=0,$ (39)

which upon substitution of the specific values of the functional derivatives
reduce to (31) - (32).

##### Casimirs

The general expression for RSWE Casimirs
$\operatorname{\mathcal{C}}[{v}^{1},\tilde{h}^{2}]$ is given by

$\operatorname{\mathcal{C}}[{v}^{1},\tilde{h}^{2}]=\left<{{h}^{0}},{F({q}^{0})}\right>,$
(40)

where $F({q}^{0})$ is an arbitrary function of potential vorticity straight
$0$-form ${q}^{0}=\tilde{I}^{0}\wedge\tilde{q}^{0}$. The functional
derivatives of $\operatorname{\mathcal{C}}[{v}^{1},\tilde{h}^{2}]$ are

$\displaystyle\frac{\tilde{\delta}\\!\operatorname{\mathcal{C}}}{\delta{v}^{1}}=\operatorname{d}(\tilde{I}^{0}\wedge
F^{\prime}),\quad\quad\frac{\tilde{\delta}\\!\operatorname{\mathcal{C}}}{\delta\tilde{h}^{2}}=F-{q}^{0}\wedge
F^{\prime},$ (41)

where $F^{\prime}=\frac{dF}{d{q}^{0}}$. Important cases are $F=1$ (total
mass), $F={q}^{0}$ (total circulation) and $F=\frac{{q}^{0}\wedge{q}^{0}}{2}$
(total potential enstrophy).

##### Vorticity Dynamics

The evolution equation for absolute vorticity ${\eta}^{2}$ is obtained by
taking $\operatorname{d}$ of (38) yielding

$\frac{\partial{\eta}^{2}}{\partial
t}=\frac{\partial(\tilde{q}^{0}\wedge\tilde{h}^{2})}{\partial
t}=\frac{\partial({q}^{0}\wedge{h}^{2})}{\partial
t}=-\operatorname{d}(\tilde{q}^{0}\wedge\tilde{F}^{n-1}),\\\ $ (42)

and for ${h}^{2}=\tilde{I}^{0}\wedge\tilde{h}^{2}$ by taking
$\tilde{I}^{0}\wedge$ of (39) resulting in

$\frac{\partial{h}^{2}}{\partial
t}+\operatorname{d}(\tilde{I}^{0}\wedge\tilde{F}^{n-1})=0.$ (43)

These can be combined with (39) to yield evolution equations for
$\tilde{q}^{0}$ and ${q}^{0}$:

$\frac{\partial\tilde{q}^{0}}{\partial
t}+\operatorname{\tilde{\star}}(\frac{1}{{h}^{0}}\wedge\tilde{F}^{n-1}\wedge\operatorname{d}\tilde{q}^{0})=0,$
(44) $\frac{\partial{q}^{0}}{\partial
t}+\operatorname{\tilde{\star}}(\frac{1}{{h}^{0}}\wedge\tilde{F}^{n-1}\wedge\operatorname{d}{q}^{0})=0.$
(45)

Note that the right-hand sides of (44) and (45) are just Lie derviatives,
since
$\mathcal{L}_{\operatorname{\mathbf{u}}}\tilde{q}^{0}=\operatorname{\tilde{\star}}(\frac{1}{{h}^{0}}\wedge\tilde{F}^{n-1}\wedge\operatorname{d}\tilde{q}^{0})$,
$\mathcal{L}_{\operatorname{\mathbf{u}}}{q}^{0}=\operatorname{\tilde{\star}}(\frac{1}{{h}^{0}}\wedge\tilde{F}^{n-1}\wedge\operatorname{d}{q}^{0})$
and
$\mathcal{L}_{\operatorname{\mathbf{u}}}{\eta}^{2}=\operatorname{d}(\tilde{I}^{0}\wedge\tilde{F}^{n-1})$.
Thus we see that $\tilde{q}^{0}$ and ${q}^{0}$ are materially conserved.

##### Linearized Equations

Following standard procedures [57] the Hamiltonian form can be linearized
around ${h}^{0}=H$, ${u}^{1}=0$ to yield for arbitrary
$\operatorname{\mathcal{A}},\operatorname{\mathcal{B}}$:

$\\{\operatorname{\mathcal{A}},\operatorname{\mathcal{B}}\\}_{lin}=-\left<\left<{\frac{\tilde{\delta}\\!\operatorname{\mathcal{A}}}{\delta{v}^{1}}},{\frac{\tilde{f}^{0}}{H}\wedge\frac{\tilde{\delta}\\!\operatorname{\mathcal{B}}}{\delta{v}^{1}}}\right>\right>-\left<\left<{\frac{\tilde{\delta}\\!\operatorname{\mathcal{A}}}{\delta{v}^{1}}},{\operatorname{d}\frac{\tilde{\delta}\\!\operatorname{\mathcal{B}}}{\delta\tilde{h}^{2}}}\right>\right>-\left<\left<{\frac{\tilde{\delta}\\!\operatorname{\mathcal{A}}}{\delta\tilde{h}^{2}}},{\operatorname{d}\frac{\tilde{\delta}\\!\operatorname{\mathcal{B}}}{\delta{v}^{1}}}\right>\right>,$
(46)

and

$\mathcal{H}_{lin}[{v}^{1},\tilde{h}^{2}]=\frac{g}{2}\left<{\tilde{h}^{2}},{\tilde{h}^{2}}\right>+g\left<{\tilde{h}^{2}},{\tilde{h_{s}}^{2}}\right>+\frac{H}{2}\left<{{u}^{1}},{{u}^{1}}\right>.$
(47)

The functional derivatives are thus

$\displaystyle\tilde{F}^{n-1}_{lin}:=\frac{\tilde{\delta}\\!\operatorname{\mathcal{H}}_{lin}}{\delta{v}^{1}}=H\tilde{u}^{n-1},\quad\quad\quad
B_{lin}:=\frac{\tilde{\delta}\\!\operatorname{\mathcal{H}}_{lin}}{\delta\tilde{h}^{2}}=g{h}^{0}+g{h_{s}}^{0},$
(48)

with equations of motion

$\displaystyle\frac{\partial{v}^{1}}{\partial
t}+\frac{\tilde{f}^{0}}{H}\wedge\tilde{F}^{n-1}_{lin}+\operatorname{d}{B}^{0}_{lin}$
$\displaystyle=$ $\displaystyle 0,$ (49)
$\displaystyle\frac{\partial\tilde{h}^{2}}{\partial
t}+\operatorname{d}\tilde{F}^{n-1}_{lin}$ $\displaystyle=$ $\displaystyle 0.$
(50)

Insertion of (48) into (49) - (50) gives (33) - (34).

##  4 A Discrete Exterior Calculus in 2D

In order to discretize the split exterior calculus form of the RSWE, we now
introduce a discrete exterior calculus (DEC) in 2D without boundaries. Much of
this material is based on [65] and [23], with some extensions. Fundamentally,
DEC is composed of discrete versions of the key elements of split exterior
calculus: straight $k$-forms ${x}^{k}$, twisted $k$-forms $\tilde{x}^{k}$,
exterior derivative $\operatorname{d}$, wedge product $\wedge$, Hodge star
$\operatorname{\tilde{\star}}$, inner product $\left<{},{}\right>$ and the
topological pairing $\left<\left<{},{}\right>\right>$ (along with functional
derviatives associated to the topological pairing). As shown in Section 8,
these operators are in fact those that underlie TRiSK-type schemes. This
connection was partially explored in [65, 23], but the importance of the wedge
product and topological pairing was not recognized. In particular, what was
missing was the formulation of the continuous RSWE in terms of split exterior
calculus, from which the full correspondence between TRiSK-type schemes and
DEC can be deduced.

In this paper we will (other than in Section 7 and 8, to facilitate comparison
with the more common TRiSK notation) denote the discrete exterior derivative
with $\mathbf{D}_{k}$ and $\mathbf{\bar{D}}_{k}$; the discrete Hodge star with
$\mathbf{H}_{k}$ and $\mathbf{\bar{H}}_{k}$ and the discrete wedge product
with ${}\boldsymbol{\wedge}_{h}{}$ (all explained further below). The discrete
inner product $\left<{},{}\right>$, topological pairing
$\left<\left<{},{}\right>\right>$ and topological functional derivative
$\frac{\tilde{\delta}\\!}{\delta}$ will share the same notation as their
continuous counterparts.

###  4.1 Grids

For DEC, a pair of related grids are used: a straight grid and a twisted grid,
each composed of $k$-cells ($k=0,1,2$). A $k$-cell is an oriented geometric
entity of dimension $k$: a $0$-cell is a point ($v$ or $\tilde{v}$), a
$1$-cell is an edge ($e$ or $\tilde{e}$) and a $2$-cell is a cell ($c$ or
$\tilde{c}$). A collection of $k$-cells for a given $k$ is a $k$-chain.

###### Remark 4.1

In the DEC literature, the terminology primal (=straight) and dual (=twisted)
are usually used. However, in the TRiSK-type scheme literature, this
convention is flipped and primal=twisted and dual=straight. In this paper we
avoid this confusion and will consistently use straight and twisted grid.

By oriented we mean that there is a value in $\\{1,-1\\}$ assigned to each of
the $(k-1)$-cells that make up the boundary of each $k$-cell, that denotes
whether this boundary $(k-1)$-cell is oriented ”towards” or ”away” from the
$k$-cell. To each $1$-cell (edge) we associate unit normal vectors
$\operatorname{\mathbf{\hat{n}_{s}}}$ (straight) and
$\operatorname{\mathbf{\hat{n}_{t}}}$ (twisted); and unit tangent vectors
$\operatorname{\mathbf{\hat{t}_{s}}}$ (straight),
$\operatorname{\mathbf{\hat{t}_{t}}}$ (twisted) such that
$\operatorname{\mathbf{\hat{n}_{s}}}\cdot\operatorname{\mathbf{\hat{t}_{s}}}=0$
and
$\operatorname{\mathbf{\hat{n}_{t}}}\cdot\operatorname{\mathbf{\hat{t}_{t}}}=0$.
These vectors can vary across the edge, for example in the case of a curved
edge. An orthogonal pair of grids will have
$\operatorname{\mathbf{\hat{n}_{s}}}=-\operatorname{\mathbf{\hat{t}_{t}}}$ and
$\operatorname{\mathbf{\hat{t}_{s}}}=\operatorname{\mathbf{\hat{n}_{t}}}$,
while a non-orthogonal pair of grids will not. An excellent description of the
grid topology and orientation in DEC can be found in [71].

From these unit vectors we can define the orientation elements $n_{ec}$,
$t_{ev}$ $\tilde{n}_{\tilde{e}\tilde{c}}$ and
$\tilde{t}_{\tilde{e}\tilde{v}}$, which all take values in $\\{-1,1\\}$.
Following exterior calculus conventions, the straight grid will be inner-
oriented, and the twisted grid will be outer-oriented. An inner-oriented grid
requires only the elements of the grid to orient itself, while an outer-
oriented grid requires a notion of embedding in a larger space; or of duality
with another grid. In fact, given an orientation for the straight grid, the
twisted grid orientation is induced by simply setting
$\tilde{t}_{\tilde{e}\tilde{v}}=n_{ec}$ and
$\tilde{n}_{\tilde{e}\tilde{c}}=-t_{ev}$ (where duality between straight
$k$-cells and twisted $(n-k)$-cells has been used, see below). Using straight
grid orientation to set twisted grid orientation ensures that the straight and
twisted grid are consistently oriented.

These oriented $k$-cells are arranged into a chain complex (also known as a CW
complex from algebraic topology) that can be described by a directed acyclic
graph (DAG) [71]. The two grids are related through topological duality: their
respective DAG’s are dual to each other in the graph theoretic sense. In
particular, duality means that there is a 1-1 relationship between straight
$k$-cells and twisted $(n-k)$-cells. This relationship is a fundamental part
of DEC, and is used extensively (for example, in the discrete Hodge star
operator). In what follows, we use $\tilde{}$ to denote the twisted grid
because it will be used to represent twisted $k$-forms, and as noted before
the $k$-cells themselves are denoted with $v,e,c$ (straight) or
$\tilde{v},\tilde{e},\tilde{c}$ (twisted) for $k=0,1,2$, respectively.
Topological duality means that there is duality between $v$ and $\tilde{c}$,
$e$ and $\tilde{e}$ and $c$ and $\tilde{v}$. The topology (including
orientation) for an example planar grid is shown in Figure 1.

Figure 1: An example straight-twisted grid without boundaries, with the
straight grid in solid blue and the twisted grid in dashed red. The unit
normal vectors $\operatorname{\mathbf{\hat{n}_{s}}}$ (straight) and
$\operatorname{\mathbf{\hat{n}_{t}}}$ (twisted); and unit tangent vectors
$\operatorname{\mathbf{\hat{t}_{s}}}$ (straight) and
$\operatorname{\mathbf{\hat{t}_{t}}}$ (twisted) for edges $e$ and $\tilde{e}$
(which are dual to each other) are also shown. The 1-1 relationship between
$k$-cells and $(n-k)$-cells is clearly visible. Some $k$-cells on both the
straight and twisted grid are labeled, which correspond to the discrete
differential forms in Figure 3.

###### Remark 4.2

In the case that $\mathbb{M}\subset\mathbb{R}^{3}$, it can be useful to define
a unit outward normal $\operatorname{\mathbf{\hat{k}}}$ to $\mathbb{M}$. This
can then be used to define $\operatorname{\mathbf{\hat{t}_{s}}}$ and
$\operatorname{\mathbf{\hat{t}_{t}}}$ using a right-hand rule as
$\operatorname{\mathbf{\hat{t}_{s}}}=\operatorname{\mathbf{\hat{k}}}\times\operatorname{\mathbf{\hat{n}_{s}}}$
and
$\operatorname{\mathbf{\hat{t}_{t}}}=\operatorname{\mathbf{\hat{k}}}\times\operatorname{\mathbf{\hat{n}_{t}}}$.
In fact, this is what is typically done for TRiSK-type schemes. This
constitutes a choice of orientation for the underlying grid. However, this is
not necessary and the discrete grids can be oriented without a
$\operatorname{\mathbf{\hat{k}}}$, including the case when
$\mathbb{M}\not\subset\mathbb{R}^{3}$.

On each grid, we can define a set of stencils that encompass the lower and
higher-dimensional nearest-neighbor $l$-cells that surround a given $k$-cell.
Specifically, for vertices $v$ we have $CV(v)$ for the nearest-neighbor cells
$c$ and $EV(v)$ for the nearest-neighbor edges; for edges $e$ we have $CE(e)$
and $VE(e)$ and for cells $c$ we have $EC(c)$ and $VC(c)$. It will also be
useful to introduce the stencil $ECP(e)$, which denotes the set of all edges
$e^{\prime}$ that are in $EC(c)$ for cells $c$ in $CE(e)$; in other words it
is the composition $EC(CE(e))$. Examples of these stencils can be found in
[65].

####  4.1.1 Grid Geometry

So far only the topological aspects of the straight and twisted grids have
been discussed. To complete their description, geometric information must be
introduced. Each geometric entity has a size associated with it: lengths
$A_{e}$ and $A_{\tilde{e}}$ for edges, areas $A_{c}$ and $A_{\tilde{c}}$ for
cells, and sizes for vertices $A_{v}$ and $A_{\tilde{v}}$ (defined to be $1$).
It is also useful to introduce the extended area $\tilde{A}_{\tilde{e}}$ and
the overlap areas $A_{c,\tilde{c}}$ and $A_{\tilde{c},\tilde{e}}$.
$\tilde{A}_{\tilde{e}}$ is the extended area associated with edge $\tilde{e}$;
while $A_{c,\tilde{c}}$ is the overlap area between straight cell $c$ and
twisted cell $\tilde{c}$ and $A_{\tilde{c},\tilde{e}}$ is the overlap area
between twisted cell $\tilde{c}$ and extended edge area
$\tilde{A}_{\tilde{e}}$. Other overlap areas and extended areas can be
similarly defined, but they are not needed for the DEC so we omit them. These
sizes and overlap areas are shown in Figure 2. In the same way that the
topology of the twisted grid can be obtained from the topology of the straight
grid through a notion of topological duality, the geometry of the twisted grid
can be obtained from the straight grid through a notion of geometric duality.
A commonly used approach is circumcentric/Voronoi duality, which places
twisted grid vertices at the circumcenters of the straight grid cells. If the
grid is optimized such that the circumcenter is also the centroid, this is a
centroidal Voronoi tessellation (CVT) [43, 38, 15]. Other approaches to
geometric duality and grid optimization include generalized power grids [24],
spring-dynamics [69, 37], barycentric [32], Hodge-optimized triangulations
[46] and tweaking [28, 29]. The choice of geometric duality is often
intimately connected to the choice of discrete Hodge star operator, for
example the Voronoi Hodge star is associated with circumcentric duality and
the barycentric Hodge star with barycentric duality. Certain choices of
geometric duality lead to orthogonality between the straight and dual grids.
It is also possible to use a grid with curved instead of straight edges: the
topology does not change but the geometry does. Curved edges are useful, for
example, to treat complex boundaries more accurately or to define a grid as a
continuous deformation applied to a regular polygonal grid, as done in [73].

Figure 2: Example geometric objects for the same straight-twisted grid as in
Figure 1. Specifically, the overlap areas $A_{\tilde{c},c}$ (stripes) and
$A_{\tilde{c}_{1},\tilde{e}_{1}}$ (waves) are shown, along with the extended
edge area $\tilde{A}_{\tilde{e}}$ (checkerboard). These quantities are used in
the definition of metric PV and KE wedge products (see definitions below). Not
illustrated are the straight $k$-cells areas $A_{v}$ (defined to be 1),
$A_{e}$ (length of edge $e$) and $A_{c}$ (area of cell $c$) and the twisted
$k$-cell areas $A_{\tilde{v}}$ (defined to be 1), $A_{\tilde{e}}$ (length of
edge $\tilde{e})$ and $A_{\tilde{c}}$ (area of cell $\tilde{c}$), which are
self-explanatory, and are used in the Voronoi Hodge star.

###  4.2 Discrete Differential Forms

A discrete differential $k$-form (${x}^{k}$ or $\tilde{x}^{k}$)888Note that we
use the same symbols for continuous and discrete quantities when there is no
danger of confusion. is represented by its action through assigning a real
number to each element of a $k$-chain: this is a cochain. In the following we
will use only the notation discrete $k$-form, not $k$-cochains. We will use
the straight grid for straight (inner-oriented) forms and the twisted grid for
twisted (outer-oriented) forms. For notation, we will use ${x}^{k}$ (resp.
$\tilde{x}^{k}$) to represent the column vector of discrete straight (resp.
twisted) $k$-forms, and ${x}^{k}_{p}$ (resp. $\tilde{x}^{k}_{p}$) to represent
the element of this vector at a specific geometric entity (where
$p\in\\{v,e,c,\tilde{v},\tilde{e},\tilde{c}\\}$ is the appropriate $k$-cell).
As needed, we will denote the space of discrete straight $k$-forms with
$\Lambda^{k}$ and discrete twisted $k$-forms with $\tilde{\Lambda}^{k}$, just
as in the continuous case. Some example discrete $k$-forms are shown in Figure
3.

Figure 3: Example discrete differential forms for the same straight-twisted
grid as in Figure 1, with the straight variables in solid blue and the twisted
variables in dashed red. Note the 1-1 duality between straight $k$-cells and
twisted $(n-k)$-cells, and therefore between straight $k$-forms and twisted
$(n-k)$-forms. The placement of quantities integrated over vertices ($0$-forms
${x}^{0}$ and $\tilde{x}^{0}$), edges ($1$-forms ${x}^{1}$ and $\tilde{x}^{1}$
and $(n-1)$-forms ${x}^{n-1}$ and $\tilde{x}^{n-1}$) and cells ($2$-forms
${x}^{2}$ and $\tilde{x}^{2}$) is illustrated, for the same $k$-cells labeled
in Figure 1.

More specifically, a discrete $k$-form can be thought of as the integration of
the associated continuous quantity (either scalar field x or vector field
$\mathbf{x}$) from vector calculus over the relevant $k$-cell:

$\begin{array}[]{rclrcl}{x}^{0}_{v}&=&x(v),&\tilde{x}^{0}_{\tilde{v}}&=&x(\tilde{v}),\\\
{x}^{1}_{e}&=&\int_{e}(\mathbf{x}\cdot\operatorname{\mathbf{\hat{t}_{s}}}dL),&\tilde{x}^{1}_{\tilde{e}}&=&\int_{\tilde{e}}(\mathbf{x}\cdot\operatorname{\mathbf{\hat{t}_{t}}}dL),\\\
{x}^{n-1}_{e}&=&\int_{e}(\mathbf{x}\cdot\operatorname{\mathbf{\hat{n}_{s}}}dL),&\tilde{x}^{n-1}_{\tilde{e}}&=&\int_{\tilde{e}}(\mathbf{x}\cdot\operatorname{\mathbf{\hat{n}_{t}}}dL),\\\
{x}^{2}_{c}&=&\int_{c}(xdA),&\tilde{x}^{2}_{\tilde{c}}&=&\int_{\tilde{c}}(xdA).\end{array}$
(51)

where $dL$ and $dA$ are the differential line and area elements of edges $e$
and $\tilde{e}$ and cells $c$ and $\tilde{c}$, respectively.

As discussed in Section 2, although we are in $n=2$ and therefore $n-1=1$,
using the notation ${x}^{1}$ and $\tilde{x}^{1}$ to represent both $1$-forms
and $(n-1)$-forms is confusing and therefore we retain the notation
${x}^{n-1}$ and $\tilde{x}^{n-1}$ to facilitate this distinction. Using $1$
and $(n-1)$ makes clear the difference between a circulation $1$-form along an
edge, and a flux $(n-1)$-form along an edge. This distinction is also seen in
the unit vectors used for integration: $\operatorname{\mathbf{\hat{t}_{s}}}$
and $\operatorname{\mathbf{\hat{t}_{t}}}$ for $1$-forms and
$\operatorname{\mathbf{\hat{n}_{s}}}$ and
$\operatorname{\mathbf{\hat{n}_{t}}}$ for $(n-1)$-forms.

###  4.3 Operators

The discrete operators are the discrete exterior derivative (denoted with
$\mathbf{D}_{k}$ and $\mathbf{\bar{D}}_{k}$); the discrete Hodge star (denote
with $\mathbf{H}_{k}$ and $\mathbf{\bar{H}}_{k}$) and the discrete wedge
product (denoted with ${}\boldsymbol{\wedge}_{h}{}$). For the first two we use
matrix notation to indicate the representation of these operators as (sparse)
matrices acting on a discrete $k$-form to produce another discrete form. In
the unary operators $\mathbf{\bar{D}}_{k}$ and $\mathbf{\bar{H}}_{k}$ the
overline indicates that the operator acts on twisted grid quantities. The
wedge product ${}\boldsymbol{\wedge}_{h}{}$ can be represented as a (sparse)
3-tensor, since it is a binary operation that takes a discrete $k$-form and a
discrete $l$-form to produce a discrete $(k+l)$-form.

Note that the scheme presented here is very general, and to close it specific
choices for the sparse coefficients and stencils appearing in the discrete
Hodge stars $\mathbf{H}_{k}$ and $\mathbf{\bar{H}}_{k}$ and wedge products
${}\boldsymbol{\wedge}_{h}{}$ must be made, ideally such that the properties
in Section 4.6 are obtained. These choices are usually specific to a given
grid topology and geometry. Some possible choices with the desired properties
are discussed in Section 7, and the specific (implicit) choices that were made
by TRiSK-type schemes in the literature are identified in Section 8.

####  4.3.1 Exterior derivative

We will use the notation

$\mathbf{D}_{k}:\Lambda^{k-1}\rightarrow\Lambda^{k},\quad\quad\mathbf{\bar{D}}_{k}:\tilde{\Lambda}^{k-1}\rightarrow\tilde{\Lambda}^{k},$
(52)

for the discrete exterior derivative. $\mathbf{D}_{k}$ and
$\mathbf{\bar{D}}_{k}$ can be represented as (sparse) matrices acting on
$(k-1)$-forms and producing a $k$-forms. These forms always belong to the same
grid. Note that this operator exists only for $k>0$. More specifically, the
discrete exterior derivatives are defined using the co-boundary operator [71].
This definition is a combinatorial discretization: it depends only on the
topology and orientation of the underlying cell complex. Since
$\operatorname{d}$ is a topological operator (one that does not require a
metric), a combinatorial definition makes sense. In terms of grid orientations
$n_{ec}$, $t_{ve}$, $\tilde{n}_{\tilde{e}\tilde{c}}$ and
$\tilde{t}_{\tilde{e}\tilde{v}}$, the exterior derivatives are written as

$\displaystyle(\mathbf{D}_{1}{x}^{0})_{e}$ $\displaystyle=$
$\displaystyle\sum_{v\in VE(e)}t_{ve}{x}^{0}_{v},$ (53)
$\displaystyle(\mathbf{D}_{2}{x}^{1})_{c}$ $\displaystyle=$
$\displaystyle\sum_{e\in EC(c)}n_{ec}{x}^{1}_{e},$ (54)
$\displaystyle(\mathbf{\bar{D}}_{1}\tilde{x}^{0})_{\tilde{e}}$
$\displaystyle=$ $\displaystyle\sum_{\tilde{v}\in
VE(\tilde{e})}\tilde{t}_{\tilde{v}\tilde{e}}\tilde{x}^{0}_{\tilde{v}},$ (55)
$\displaystyle(\mathbf{\bar{D}}_{2}\tilde{x}^{1})_{\tilde{c}}$
$\displaystyle=$ $\displaystyle\sum_{\tilde{e}\in
EC(\tilde{c})}\tilde{n}_{\tilde{e}\tilde{c}}\tilde{x}^{1}_{\tilde{e}}$ (56)

More generally, the discrete exterior derivative for a discrete $k$-form is
simply the weighted sum of nearest-neighbor discrete $(k-1)$-forms, with the
weights $w\in\\{-1,1\\}$ given by the orientation. This definition is
dimension and form degree independent.

####  4.3.2 Wedge product

We will use the notation

$\displaystyle{{x}^{k}}\boldsymbol{\wedge}_{h}{{y}^{l}}:\Lambda^{k}\times\Lambda^{l}\rightarrow\Lambda^{k+l},$
(57)
$\displaystyle{\tilde{x}^{k}}\boldsymbol{\wedge}_{h}{{y}^{l}}:\tilde{\Lambda}^{k}\times\Lambda^{l}\rightarrow\tilde{\Lambda}^{k+l},$
(58)
$\displaystyle{{x}^{k}}\boldsymbol{\wedge}_{h}{\tilde{y}^{l}}:\Lambda^{k}\times\tilde{\Lambda}^{l}\rightarrow\tilde{\Lambda}^{k+l},$
(59)
$\displaystyle{\tilde{x}^{k}}\boldsymbol{\wedge}_{h}{\tilde{y}^{l}}:\tilde{\Lambda}^{k}\times\tilde{\Lambda}^{l}\rightarrow\Lambda^{k+l},$
(60)

to represent the discrete wedge product. This can be represented as a (sparse)
3-tensor acting on $k$-forms and $l$-forms to produce a $(k+l)$-forms. These
forms can belong to the same grid, or different grids, depending on the
specific wedge product. Note that this operator exists only when $k+l\leq n$.
In algebraic topology, this operator is known as the cup product [40]. As an
example, consider
${z}^{1}={\tilde{x}^{0}}\boldsymbol{\wedge}_{h}{\tilde{y}^{1}}$. An explicit
form for this wedge product is

${z}^{1}_{e}=\sum_{\tilde{v}}\sum_{\tilde{e}}W_{e\tilde{v}\tilde{e}}\tilde{x}^{0}_{\tilde{v}}\tilde{y}^{1}_{\tilde{e}},$
(61)

for some set of arbitrary (sparse) coefficients $W_{e\tilde{v}\tilde{e}}$,
where the target is in the first slot, and the two sources arguments are in
the second and third slots.

For TRiSK-type schemes, there are in fact only a few wedge products that must
be considered. Specifically, to compute the nonlinear PV flux term and PV
itself requires, respectively:

${z}^{1}={\tilde{x}^{0}}\boldsymbol{\wedge}_{h}{\tilde{y}^{1}},\quad\quad{z}^{2}={\tilde{x}^{0}}\boldsymbol{\wedge}_{h}{\tilde{y}^{2}}.$
(62)

To compute the kinetic energy part of the Hamiltonian and the associated mass
flux and kinetic energy functional derivatives requires

$\tilde{z}^{2}={{x}^{1}}\boldsymbol{\wedge}_{h}{\tilde{y}^{1}},\quad\quad\tilde{z}^{1}={{x}^{0}}\boldsymbol{\wedge}_{h}{\tilde{y}^{1}},\quad\quad{z}^{1}={{x}^{0}}\boldsymbol{\wedge}_{h}{{y}^{1}},$
(63)

These cannot be chosen independently, instead one chooses one of three
possible forms for the kinetic energy part of the Hamiltonian (each of which
requires one of these wedge products), and then in the process of taking
functional derivatives the other two wedge products arise in terms of the
adjoint (see Section 4.5 and Appendix C).

####  4.3.3 Hodge star

We will use the notation

$\mathbf{H}_{k}:\Lambda^{k}\rightarrow\tilde{\Lambda}^{n-k},\quad\quad\mathbf{\bar{H}}_{k}:\tilde{\Lambda}^{k}\rightarrow\Lambda^{n-k},$
(64)

to represent the discrete Hodge star. It can be represented as a (sparse)
matrix acting on a $k$-form from one grid and producing a $(n-k)$-form on the
other grid. By the duality between straight (twisted) $k$-cells and twisted
(straight) $(n-k)$-cells, these forms have the same number of degrees of
freedom. As an example, consider $\tilde{y}^{n-1}=\mathbf{H}_{1}{x}^{1}$. This
can be written as

$\tilde{y}^{n-1}_{\tilde{e}}=\sum_{e}H_{\tilde{e}e}{x}^{1}_{e}$ (65)

for some set of arbitrary (sparse) coefficients $H_{\tilde{e}e}$, where the
target is in the first slot and the source is in the second slot.

For TRiSK-type schemes, only three discrete Hodge stars are needed:

$\tilde{y}^{1}=\mathbf{H}_{1}{x}^{1},\quad\quad{y}^{0}=\mathbf{\bar{H}}_{2}\tilde{x}^{2},\quad\quad\tilde{y}^{0}=\mathbf{H}_{2}{x}^{2}.$
(66)

The remaining Hodge star operators can be implicitly defined in terms of these
operators as

$\displaystyle\mathbf{H}_{0}:=\mathbf{\bar{H}}_{2}^{-1},\quad\quad\mathbf{\bar{H}}_{0}:=\mathbf{H}_{2}^{-1},\quad\quad\mathbf{\bar{H}}_{1}:=-\mathbf{H}_{1}^{-1}.$
(67)

This definition is chosen so that

$\mathbf{\bar{H}}_{k}\mathbf{H}_{n-k}=(-1)^{k(n-k)}\mathbf{I}$ (68)

which is a discrete analogue of
$\operatorname{\tilde{\star}}\operatorname{\tilde{\star}}=(-1)^{k(n-k)}$. The
sign in the last definition is a slight difference from [65]. This definition
requires that $\mathbf{H}_{1}$, $\mathbf{\bar{H}}_{2}$ and $\mathbf{H}_{2}$
are invertible, and for computational efficiency should also be sparse.
Although $\mathbf{H}_{0}$, $\mathbf{\bar{H}}_{0}$ and $\mathbf{\bar{H}}_{1}$
will be dense (unless $\mathbf{H}_{1}$, $\mathbf{\bar{H}}_{2}$ and
$\mathbf{H}_{2}$ are diagonal, such as in the case of the Voronoi Hodge star),
this is not an issue since they are not needed in the general DEC-based TRiSK-
type scheme. We could of course have started with a different set of primary
Hodge stars if the scheme required it, the only requirement is that
$\mathbf{\bar{H}}_{k}\mathbf{H}_{n-k}=(-1)^{k(n-k)}\mathbf{I}$ holds.

###  4.4 Topological Pairing, Functional Derivatives and Inner Products

It remains to define the inner product $\left<{},{}\right>$, topological
pairing $\left<\left<{},{}\right>\right>$ and functional derivatives (with
respect to the topological pairing) $\frac{\tilde{\delta}\\!}{\delta}$. These
definitions are done by using the discrete Hodge stars $\mathbf{H}_{k}$ and
$\mathbf{\bar{H}}_{k}$, in a way that mimics
$\left<{{a}^{k}},{{b}^{k}}\right>=\int{a}^{k}\wedge\operatorname{\tilde{\star}}{b}^{k}$
and
$\left<{\tilde{a}^{k}},{\tilde{b}^{k}}\right>=\int\tilde{a}^{k}\wedge\operatorname{\tilde{\star}}\tilde{b}^{k}$.
For all three operators, we will use the same notation for the discrete
objects as their continuous counterparts.

Start with the inner product $\left<{},{}\right>$ for straight and twisted
forms:

$\left<{{x}^{k}},{{y}^{k}}\right>:=({x}^{k})^{T}\mathbf{H}_{k}{y}^{k},\quad\quad\left<{\tilde{x}^{k}},{\tilde{y}^{k}}\right>:=(-1)^{k(n-k)}(\tilde{x}^{k})^{T}\mathbf{\bar{H}}_{k}\tilde{y}^{k}.$
(69)

If $\mathbf{H}_{k}$ and $\mathbf{\bar{H}}_{k}$ are positive-definite, this
inner product will be positive-definitive. Some care is required with
$\mathbf{\bar{H}}_{1}$, as it is actually negative-definite (assuming
$\mathbf{H}_{1}$ is positive-definite) since
$-\mathbf{\bar{H}}_{1}=\mathbf{H}_{1}^{-1}$. However, the minus sign for
twisted 1-forms in the above definition ensure that the resulting inner
product is still positive-definite.

Additionally, if $\mathbf{H}_{k}$ and $\mathbf{\bar{H}}_{k}$ are symmetric,
then

$\displaystyle\left<{{x}^{k}},{{y}^{k}}\right>$ $\displaystyle=$
$\displaystyle\left<{{y}^{k}},{{x}^{k}}\right>$ (70)
$\displaystyle\left<{\tilde{x}^{k}},{\tilde{y}^{k}}\right>$ $\displaystyle=$
$\displaystyle\left<{\tilde{y}^{k}},{\tilde{x}^{k}}\right>$ (71)

just as in the continuous case.

Now we can define the topological pairing using these inner products
(similarly to (7) - (8)) as

$\displaystyle\left<\left<{{a}^{k}},{\tilde{b}^{n-k}}\right>\right>$
$\displaystyle:=\left<{{a}^{k}},{{b}^{k}}\right>=({a}^{k})^{T}\mathbf{H}_{k}{b}^{k}=({a}^{k})^{T}\tilde{b}^{n-k},$
(72) $\displaystyle\left<\left<{\tilde{a}^{k}},{{b}^{n-k}}\right>\right>\\!$
$\displaystyle:=\\!(-1)^{k(n-k)}\left<{\tilde{a}^{k}},{\tilde{b}^{k}}\right>=(\tilde{a}^{k})^{T}\mathbf{\bar{H}}_{k}\tilde{b}^{k}=(\tilde{a}^{k})^{T}\mathbf{\bar{H}}_{k}\mathbf{H}_{n-k}{b}^{n-k}=(-1)^{k(n-k)}(\tilde{a}^{k})^{T}{b}^{n-k},$
(73)

where $\tilde{b}^{n-k}=\mathbf{H}_{k}{b}^{k}$ and
$\tilde{b}^{k}=\mathbf{H}_{n-k}{b}^{n-k}$. The last definition here relied on
$\mathbf{\bar{H}}_{k}\mathbf{H}_{n-k}=(-1)^{k(n-k)}\mathbf{I}$.

This definition relies on the 1-1 relationship between straight and twisted
quantities. In particular, these definitions ensure that the discrete
analogues of the continuous relationships

$\displaystyle\left<\left<{{a}^{k}},{\tilde{b}^{n-k}}\right>\right>$
$\displaystyle\stackrel{{\scriptstyle\eqref{discrete-topo-
pairing-1}\&\eqref{discrete-topo-pairing-2}}}{{=}}$
$\displaystyle(-1)^{k(n-k)}\left<\left<{\tilde{b}^{n-k}},{{a}^{k}}\right>\right>,$
(74) $\displaystyle\left<\left<{{a}^{k}},{\tilde{b}^{n-k}}\right>\right>$
$\displaystyle\stackrel{{\scriptstyle\eqref{discrete-topo-
pairing-1}\&\eqref{stf-inner-symmetric}}}{{=}}$
$\displaystyle\left<\left<{{b}^{k}},{\tilde{a}^{n-k}}\right>\right>,$ (75)
$\displaystyle\left<\left<{{a}^{k}},{\tilde{b}^{n-k}}\right>\right>$
$\displaystyle\stackrel{{\scriptstyle\eqref{discrete-topo-pairing-1}}}{{=}}$
$\displaystyle\left<{{a}^{k}},{{b}^{k}}\right>,$ (76)
$\displaystyle\left<\left<{\tilde{a}^{k}},{{b}^{n-k}}\right>\right>$
$\displaystyle\stackrel{{\scriptstyle\eqref{discrete-topo-pairing-2}}}{{=}}$
$\displaystyle(-1)^{k(n-k)}\left<{\tilde{a}^{k}},{\tilde{b}^{k}}\right>,$ (77)

hold. Additionally, the topological pairing is a topological operator, so the
definition used above (which is combinatorial) also makes sense.

Finally, using this discrete topological pairing, discrete topological
functional derivatives can be defined using

$\displaystyle\delta\operatorname{\mathcal{F}}$ $\displaystyle:=$
$\displaystyle\left<\left<{{\omega}^{k}},{\frac{\tilde{\delta}\\!\operatorname{\mathcal{F}}}{\delta{x}^{k}}}\right>\right>=({\omega}^{k})^{T}\frac{\tilde{\delta}\\!\operatorname{\mathcal{F}}}{\delta{x}^{k}}\quad\forall{\omega}^{k}\in\Lambda^{k},$
(78) $\displaystyle\delta\operatorname{\mathcal{F}}$ $\displaystyle:=$
$\displaystyle\left<\left<{\tilde{\omega}^{k}},{\frac{\tilde{\delta}\\!\operatorname{\mathcal{F}}}{\delta\tilde{y}^{k}}}\right>\right>=(-1)^{k(n-k)}(\tilde{\omega}^{k})^{T}\frac{\tilde{\delta}\\!\operatorname{\mathcal{F}}}{\delta\tilde{y}^{k}}\quad\forall\tilde{\omega}^{k}\in\Lambda^{k}.$
(79)

As in the continuous case, a key point is that the functional derivatives lie
in the dual space (defined through the Hodge star) to the variable they are
taken with respect to, i.e.

$\frac{\tilde{\delta}\\!\operatorname{\mathcal{A}}}{\delta{x}^{k}}\in\tilde{\Lambda}^{n-k},\quad\quad\frac{\tilde{\delta}\\!\operatorname{\mathcal{A}}}{\delta\tilde{y}^{k}}\in\Lambda^{n-k}.$
(80)

Appendix C gives an example of how this definition is used to derive
functional derivatives with respect to the topological pairing.

###  4.5 Operator Adjoints

It is also useful to consider the adjoints with respect to the topological
pairing of the discrete operators defined above. The adjoints will be denoted
with a superscript $\star$: $\mathbf{D}_{k}^{\star}$,
$\mathbf{\bar{D}}_{k}^{\star}$, ${}\boldsymbol{\wedge}_{h}^{\star}{}$,
$\mathbf{H}_{k}^{\star}$ and $\mathbf{\bar{H}}_{k}^{\star}$. Note that for the
wedge product, there are actually two adjoints to consider since it is a
binary operator. These adjoints are topological adjoints, in that they can be
defined independently of a metric since they are with respect to the
topological pairing.

In the continuous setting, the topological adjoint of a unary operator
$\mathbf{X}$ is defined through

$\left<\left<{\mathbf{X}^{*}{x}^{k}},{{y}^{l}}\right>\right>:=\left<\left<{{x}^{k}},{\mathbf{X}{y}^{l}}\right>\right>.$
(81)

where $k$ and $l$ are defined such that the resulting topological pairing
makes sense (${y}^{l}$ or ${x}^{k}$ might even need to become twisted forms).
In the discrete setting, this is written as

$\left<\left<{{x}^{k}},{\mathbf{X}{y}^{l}}\right>\right>=({x}^{k})^{T}\mathbf{X}{y}^{l}=(\mathbf{X}^{T}{x}^{k})^{T}{y}^{l}=\left<\left<{\mathbf{X}^{T}{x}^{k}},{{y}^{l}}\right>\right>=\left<\left<{\mathbf{X}^{*}{x}^{k}},{{y}^{l}}\right>\right>,$
(82)

and thus we see that $\mathbf{X}^{*}=\mathbf{X}^{T}$, or in other words the
coefficients that define the discrete adjoint are just the transpose of those
that define the original operator. Given $\mathbf{H}_{1}$, for example,
$\mathbf{H}_{1}^{*}$ is defined by

$H_{e\tilde{e}}^{*}=H_{\tilde{e}e}.$ (83)

For a binary operator such as the wedge product, there are in fact two
adjoints depending on which argument we are taking the adjoint with respect
to. For example, consider the adjoint with respect to the second argument:

$\left<\left<{\tilde{x}^{m}},{{{y}^{k}}\boldsymbol{\wedge}_{h}{{z}^{l}}}\right>\right>=\left<\left<{{z}^{l}},{{{y}^{k}}\boldsymbol{\wedge}_{h}^{\star}{\tilde{x}^{m}}}\right>\right>,$
(84)

for $m=n-(k+l)$. In the discrete setting this is

$\left<\left<{\tilde{x}^{m}},{{{y}^{k}}\boldsymbol{\wedge}_{h}{{z}^{l}}}\right>\right>=\sum_{m^{\prime}}\tilde{x}^{m}_{m^{\prime}}\sum_{k^{\prime}}\sum_{l^{\prime}}W_{m^{\prime}k^{\prime}l^{\prime}}{y}^{k}_{k^{\prime}}{z}^{l}_{l^{\prime}}=\sum_{l^{\prime}}{z}^{l}_{l^{\prime}}\sum_{m^{\prime}}\sum_{k^{\prime}}W_{l^{\prime}k^{\prime}m^{\prime}}^{*}{y}^{k}_{k^{\prime}}\tilde{x}^{m}_{m^{\prime}}=\left<\left<{{z}^{l}},{{{y}^{k}}\boldsymbol{\wedge}_{h}^{\star}{\tilde{x}^{m}}}\right>\right>,$
(85)

and thus we see that $W_{lkm}^{*}=W_{mkl}$, i.e. a transposition of the
relevant indices of the 3-tensor.

###  4.6 Operator Properties

It is important for the discrete operators to have (some of) the same
properties as their continuous counterparts. These properties are required in
order for the discretization of the RSWE in Section 5 to have the desirable
properties discussed in Section 6.

####  4.6.1 Exterior Derivative and Topological Pairing

The discrete exterior derivatives $\mathbf{D}_{k}$ and $\mathbf{\bar{D}}_{k}$
for $k>0$ (defined using the co-boundary operator) and the topological pairing
$\left<\left<{},{}\right>\right>$ have the following set of properties:

  * •

Annihilation:

$\mathbf{D}_{k}\mathbf{D}_{k-1}=0,\quad\quad\quad\mathbf{\bar{D}}_{k}\mathbf{\bar{D}}_{k-1}=0;$
(86)

Equation (86) is the discrete analogue of
$\operatorname{d}\operatorname{d}=0$.

  * •

Integration by parts (IBP):

$\displaystyle\left<\left<{\mathbf{D}_{k}{x}^{k-1}},{\tilde{y}^{n-k}}\right>\right>+(-1)^{k-1}\left<\left<{{x}^{k-1}},{\mathbf{\bar{D}}_{n-k+1}\tilde{y}^{n-k}}\right>\right>=0,$
(87)
$\displaystyle\left<\left<{\mathbf{\bar{D}}_{k}\tilde{x}^{k-1}},{{y}^{n-k}}\right>\right>+(-1)^{k-1}\left<\left<{\tilde{x}^{k-1}},{\mathbf{D}_{n-k+1}{y}^{n-k}}\right>\right>=0.$
(88)

which is equivalent (for the case of $n=2$) to

$\mathbf{\bar{D}}_{2}=-\mathbf{D}_{1}^{T},\quad\quad\quad\mathbf{D}_{2}=\mathbf{\bar{D}}_{1}^{T};$
(89)

Equations (87) and (88) are the discrete analogue of

$\displaystyle\left<\left<{\operatorname{d}{x}^{k-1}},{\tilde{y}^{n-k}}\right>\right>+(-1)^{k-1}\left<\left<{{x}^{k-1}},{\operatorname{d}\tilde{y}^{n-k}}\right>\right>=0,$
(90)
$\displaystyle\left<\left<{\operatorname{d}\tilde{x}^{k-1}},{{y}^{n-k}}\right>\right>+(-1)^{k-1}\left<\left<{\tilde{x}^{k-1}},{\operatorname{d}{y}^{n-k}}\right>\right>=0.$
(91)

  * •

Zero exterior derivative for constants:

$\mathbf{D}_{1}{c}^{0}=0,\quad\quad\quad\mathbf{\bar{D}}_{1}\tilde{c}^{0}=0,$
(92)

where ${c}^{0}$ and $\tilde{c}^{0}$ are constants. Equation (92) is the
discrete analogue of $\operatorname{d}{c}^{0}=0$ and
$\operatorname{d}\tilde{c}^{0}=0$. Note these last two properties (zero
exterior derivative for constants and integration by parts) imply:

$(\tilde{I}^{0})^{T}\mathbf{D}_{2}{x}^{1}=0,\quad\quad({I}^{0})^{T}\mathbf{\bar{D}}_{2}\tilde{x}^{1}=0,$
(93)

which is the discrete analogue of Stokes theorem
$\int_{\Omega}\operatorname{d}{x}^{k}=\int_{\delta\Omega}\text{tr }{x}^{k}=0$.

A proof of these properties for $\mathbf{D}_{k}$ and $\mathbf{\bar{D}}_{k}$
based on the coboundary operator is found in Appendix A.

####  4.6.2 Wedge Product

The continuous wedge product has three key properties:

  * •

Leibniz rule:
$\operatorname{d}({x}^{k}\wedge{y}^{l})=\operatorname{d}{x}^{k}\wedge{y}^{l}+(-1)^{k}{x}^{k}\wedge\operatorname{d}{y}^{l}$
,

  * •

Anti-symmetry: ${x}^{k}\wedge{y}^{l}=(-1)^{kl}{y}^{l}\wedge{x}^{k}$ ,

  * •

Associative:
${x}^{m}\wedge({y}^{k}\wedge{z}^{l})=({x}^{m}\wedge{y}^{k})\wedge{z}^{l}$ .

Unfortunately, it has been proved using algebraic topology that it is
impossible to define a discrete wedge product (i.e. a cup product) with all
three properties [40]. However, since only the first two (anti-symmetry and
the Leibniz rule) are actually needed for TRiSK-type schemes, this is not an
impediment to practical usage of DEC in this case.

In fact, only the wedge products used to compute the nonlinear PV flux term
and PV itself, namely

${z}^{1}={\tilde{x}^{0}}\boldsymbol{\wedge}_{h}{\tilde{y}^{1}},\quad\quad{z}^{2}={\tilde{x}^{0}}\boldsymbol{\wedge}_{h}{\tilde{y}^{2}},$
(94)

respectively, are required to have these properties. Specifically, we would
like

  * •

Anti-symmetry for ${\tilde{x}^{0}}\boldsymbol{\wedge}_{h}{\tilde{y}^{1}}$:
Namely, that

$\left<\left<{\tilde{z}^{1}},{{\tilde{x}^{0}}\boldsymbol{\wedge}_{h}{\tilde{y}^{1}}}\right>\right>=-\left<\left<{\tilde{y}^{1}},{{\tilde{x}^{0}}\boldsymbol{\wedge}_{h}{\tilde{z}^{1}}}\right>\right>.$
(95)

In other words, the coefficients that define
${\tilde{x}^{0}}\boldsymbol{\wedge}_{h}{\tilde{y}^{1}}$ must be antisymmetric
in the arguments corresponding to $\tilde{z}^{1}$ and $\tilde{y}^{1}$. (95)
ensures that the (nonlinear) PV flux term conserves energy (as shown later).
The vector calculus analogue of this is that
$\mathbf{x}\cdot\mathbf{x}^{\perp}=0$.

  * •

(Partial) Leibniz rule I for $\operatorname{d}$ and $\wedge$: Namely, that

$\mathbf{D}_{2}({\tilde{I}^{0}}\boldsymbol{\wedge}_{h}{\tilde{y}^{1}})={\tilde{I}^{0}}\boldsymbol{\wedge}_{h}{\mathbf{\bar{D}}_{2}\tilde{y}^{1}}.$
(96)

Equation (96) ensures PV compatibility and steady geostrophic modes (also
shown below) The vector calculus analogue of this is that
$\nabla\cdot\mathbf{x}=\nabla^{\perp}\cdot\mathbf{x}^{\perp}$.

  * •

(Partial) Leibniz rule II for $\operatorname{d}$ and $\wedge$ II: Namely, that

${\tilde{x}^{0}}\boldsymbol{\wedge}_{h}{\mathbf{\bar{D}}_{1}\tilde{x}^{0}}-\frac{1}{2}\mathbf{D}_{1}({\tilde{x}^{0}}\boldsymbol{\wedge}_{h}^{\star}{\tilde{x}^{0}})=0.$
(97)

Equation (97) ensures potential enstrophy conservation (shown below). The
vector calculus analogue of this is that
$q(\nabla^{\perp}q)^{\perp}+\nabla\frac{q^{2}}{2}=0$.

  * •

The wedge product
${z}^{2}={\tilde{x}^{0}}\boldsymbol{\wedge}_{h}{\tilde{y}^{2}}$ has a
solvability restriction on the coefficients: they must be chosen such that
given ${z}^{2}$ and $\tilde{y}^{2}$ an explicit formula for $\tilde{x}^{0}$ is
available, otherwise a linear system must be solved and the resulting scheme
will be inefficient. An implicit formula for $\tilde{x}^{0}$ will not affect
any of it’s desirable properties from Section 6, however.

An interesting question that arises is whether a connection exists between the
two partial Leibniz rules (96) and (97), motivated by the fact that in [23] a
discrete wedge product that satisfied (97) and (95) was found to also
automatically satisfy (96). This is explored in Appendix B.

####  4.6.3 Hodge Star

The discrete Hodge stars $\mathbf{\bar{H}}_{k}$ and $\mathbf{H}_{k}$ should
satisfy $\mathbf{\bar{H}}_{n-k}\mathbf{H}_{k}=(-1)^{k(n-k)}\mathbf{I}$, which
is ensured by the construction above in Section 4.3.3. Additionally,
$\mathbf{H}_{k}$ and $\mathbf{\bar{H}}_{k}$ must be symmetric positive-
definite for the inner product $\left<{},{}\right>$ to have the appropriate
properties to make it an inner product, i.e. $\left<{a},{b}\right>\geq 0$ with
equality only if $a=0$ or $b=0$ along with
$\left<{a},{b}\right>=\left<{b},{a}\right>$. This is clearly seen in (69).

##  5 Discretization of RSWE with DEC

Equipped with the 2D discrete exterior calculus from Section 4, it is
straightforward to spatially discretize the split exterior calculus
formulation of the RSWE from Section 2. Exactly as in the continuous case, the
predicted discrete variables are the absolute velocity ${v}^{1}$ and fluid
height $\tilde{h}^{2}$. The discretization simply consists in replacing $d$
with $\mathbf{D}_{k}$ or $\mathbf{\bar{D}}_{k}$, $\wedge$ with
${}\boldsymbol{\wedge}_{h}{}$ or ${}\boldsymbol{\wedge}_{h}^{\star}{}$,
$\operatorname{\tilde{\star}}$ with $\mathbf{H}_{k}$ or $\mathbf{\bar{H}}_{k}$
and using discrete versions of the inner product $\left<{},{}\right>$ and
topological pairing $\left<\left<{},{}\right>\right>$. The properties of these
operators will ensure that the discretization inherits most of the key
properties of the continuous equations.

###  5.1 Discrete Hamiltonian and Functional Derivatives

The discrete Hamiltonian $\operatorname{\mathcal{H}}[{v}^{1},\tilde{h}^{2}]$
is

$\operatorname{\mathcal{H}}[{v}^{1},\tilde{h}^{2}]=\frac{g}{2}\left<{\tilde{h}^{2}},{\tilde{h}^{2}}\right>+g\left<{\tilde{h}^{2}},{\tilde{h_{s}}^{2}}\right>+\left<{\tilde{h}^{2}},{\frac{{{u}^{1}}\boldsymbol{\wedge}_{h}{\tilde{u}^{n-1}}}{2}}\right>,$
(98)

which has functional derivatives relative to the topological pairing (6)

$\tilde{F}^{n-1}:=\frac{\tilde{\delta}\\!\operatorname{\mathcal{H}}}{\delta{v}^{1}}=\frac{1}{2}{{h}^{0}}\boldsymbol{\wedge}_{h}^{\star}{\mathbf{H}_{1}{u}^{1}}+\frac{1}{2}\mathbf{H}_{1}({{h}^{0}}\boldsymbol{\wedge}_{h}^{\star}{{u}^{1}}),\quad\quad{B}^{0}:=\frac{\tilde{\delta}\\!\operatorname{\mathcal{H}}}{\delta\tilde{h}^{2}}=\frac{1}{2}\mathbf{\bar{H}}_{2}({{u}^{1}}\boldsymbol{\wedge}_{h}{\tilde{u}^{n-1}})+g({h}^{0}+{h_{s}}^{0}),$
(99)

with topography $\tilde{h_{s}}^{2}$, recalling
${h_{s}}^{0}=\mathbf{\bar{H}}_{2}\tilde{h_{s}}^{2}$,
${h}^{0}=\mathbf{\bar{H}}_{2}\tilde{h}^{2}$ and
$\tilde{u}^{n-1}=\mathbf{H}_{1}{u}^{1}$. The derivation of these can be found
in Appendix C.

Thus we see that the discrete kinetic energy is
$\tilde{K}^{2}=\frac{1}{2}({{u}^{1}}\boldsymbol{\wedge}_{h}{\tilde{u}^{n-1}})$,
just as in the continuous case. The mass flux is equivalent to the continuous
form
$\frac{1}{2}{h}^{0}\wedge\tilde{u}^{n-1}+\frac{1}{2}\operatorname{\tilde{\star}}({h}^{0}\wedge{u}^{1})={h}^{0}\wedge\tilde{u}^{n-1}$.
For certain choices of $\mathbf{H}_{1}$ and ${}\boldsymbol{\wedge}_{h}{}$ this
simplification is also true in the discrete setting, i.e.
$\mathbf{H}_{1}({{h}^{0}}\boldsymbol{\wedge}_{h}{{u}^{1}})={{h}^{0}}\boldsymbol{\wedge}_{h}{\mathbf{H}_{1}{u}^{1}}$
and then
$\tilde{F}^{n-1}={{h}^{0}}\boldsymbol{\wedge}_{h}{\mathbf{H}_{1}{u}^{1}}={{h}^{0}}\boldsymbol{\wedge}_{h}{\tilde{u}^{n-1}}$.
This is the case for the Voronoi Hodge star and the metric or combinatorial
wedge products from Section 6, as shown in Appendix C.1. However, it will not
be true in general.

Just as in the continuous case, the last term in
$\operatorname{\mathcal{H}}[{v}^{1},\tilde{h}^{2}]$ (the kinetic energy) can
also be written as
$\frac{1}{2}\left<{{u}^{1}},{{{h}^{0}}\boldsymbol{\wedge}_{h}{{u}^{1}}}\right>$
or
$\frac{1}{2}\left<{\tilde{u}^{n-1}},{{{h}^{0}}\boldsymbol{\wedge}_{h}{\tilde{u}^{n-1}}}\right>$,
leading respectively to slightly different functional derivatives

$\tilde{F}^{n-1}:=\frac{\tilde{\delta}\\!\operatorname{\mathcal{H}}}{\delta{v}^{1}}=\frac{1}{2}{{h}^{0}}\boldsymbol{\wedge}_{h}^{\star}{\mathbf{H}_{1}{u}^{n-1}}+\frac{1}{2}\mathbf{H}_{1}({{h}^{0}}\boldsymbol{\wedge}_{h}{{u}^{1}}),\quad\quad{B}^{0}:=\frac{\tilde{\delta}\\!\operatorname{\mathcal{H}}}{\delta\tilde{h}^{2}}=\frac{1}{2}\mathbf{\bar{H}}_{2}({{u}^{1}}\boldsymbol{\wedge}_{h}^{\star}{\tilde{u}^{n-1}})+g({h}^{0}+{h_{s}}^{0}),$
(100)

and

$\tilde{F}^{n-1}:=\frac{\tilde{\delta}\\!\operatorname{\mathcal{H}}}{\delta{v}^{1}}=\frac{1}{2}{{h}^{0}}\boldsymbol{\wedge}_{h}{\mathbf{H}_{1}{u}^{n-1}}+\frac{1}{2}\mathbf{H}_{1}({{h}^{0}}\boldsymbol{\wedge}_{h}^{\star}{{u}^{1}}),\quad\quad{B}^{0}:=\frac{\tilde{\delta}\\!\operatorname{\mathcal{H}}}{\delta\tilde{h}^{2}}=\frac{1}{2}\mathbf{\bar{H}}_{2}({{u}^{1}}\boldsymbol{\wedge}_{h}^{\star}{\tilde{u}^{n-1}})+g({h}^{0}+{h_{s}}^{0}).$
(101)

These functional derivatives have the same form as (99), the only difference
is which wedge products are adjoints (${}\boldsymbol{\wedge}_{h}^{\star}{}$)
and which ones are primary (${}\boldsymbol{\wedge}_{h}{}$).

###  5.2 Discrete Poisson Brackets

The discrete Poisson brackets are

$\\{\operatorname{\mathcal{A}},\operatorname{\mathcal{B}}\\}=-\left<\left<{\frac{\tilde{\delta}\\!\operatorname{\mathcal{A}}}{\delta{v}^{1}}},{{\tilde{q}^{0}}\boldsymbol{\wedge}_{h}{\frac{\tilde{\delta}\\!\operatorname{\mathcal{B}}}{\delta{v}^{1}}}}\right>\right>-\left<\left<{\frac{\tilde{\delta}\\!\operatorname{\mathcal{A}}}{\delta{v}^{1}}},{\mathbf{D}_{1}\frac{\tilde{\delta}\\!\operatorname{\mathcal{B}}}{\delta\tilde{h}^{2}}}\right>\right>-\left<\left<{\frac{\tilde{\delta}\\!\operatorname{\mathcal{A}}}{\delta\tilde{h}^{2}}},{\mathbf{\bar{D}}_{2}\frac{\tilde{\delta}\\!\operatorname{\mathcal{B}}}{\delta{v}^{1}}}\right>\right>,$
(102)

with potential vorticity $\tilde{q}^{0}$ defined through

${\tilde{q}^{0}}\boldsymbol{\wedge}_{h}{\tilde{h}^{2}}={\eta}^{2}=\mathbf{D}_{2}{v}^{1}={\zeta}^{2}+{f}^{2}=\mathbf{D}_{2}{u}^{1}+{f}^{2}.$
(103)

The definition of $\tilde{q}^{0}$ highlights the solvability restriction on
the coefficients for the ${}\boldsymbol{\wedge}_{h}{}$ used to compute
$\tilde{q}^{0}$, since $\tilde{q}^{0}$ should be computable without requiring
a linear solve in order for the scheme to be efficient.

###  5.3 Semi-discrete Equations of Motion

The semi-discrete equations of motion that arise from these brackets and
Hamiltonian are therefore

$\displaystyle\frac{\partial\tilde{h}^{2}}{\partial
t}+\mathbf{\bar{D}}_{2}\tilde{F}^{n-1}=0,$ (104)
$\displaystyle\frac{\partial{v}^{1}}{\partial
t}+{\tilde{q}^{0}}\boldsymbol{\wedge}_{h}{\tilde{F}^{n-1}}+\mathbf{D}_{1}{B}^{0}=0.$
(105)

These are obtained by using
$\frac{\operatorname{d}\\!\operatorname{\mathcal{F}}}{\operatorname{d}\\!t}=\\{\operatorname{\mathcal{F}},\operatorname{\mathcal{H}}\\}$
for a functional $\operatorname{\mathcal{F}}=x_{p}$, since in this case
$\frac{\delta\\!\operatorname{\mathcal{F}}}{\delta x_{p}}=1$ and
$\frac{\delta\\!\operatorname{\mathcal{F}}}{\delta x}=0$ for all other $x$;
and repeating this process for all
$x_{p}\in(\tilde{h}^{2}_{\tilde{c}},{v}^{1}_{e})$. This is the discrete
analogue of letting $\operatorname{\mathcal{F}}=\left<\hat{x},x\right>$ for
$x\in(\tilde{h}^{2},{v}^{1})$ with arbitrary test function
$\hat{x}\in(\widehat{\tilde{h}^{2}},\widehat{{v}^{1}})$. Standard time
integration schemes can then be applied to obtain the fully discrete scheme,
for example explicit Runge-Kutta schemes or an (implicit) energy-conserving
Poisson integrator [11] that preserves properties such as energy conservation
in the fully discrete case.

###  5.4 Linearized Equations

Similar to Eqns.(38)–(39), we linearize the above system around ${h}^{0}=H$
and ${u}^{1}=0$ to get

$\displaystyle\frac{\partial{v}^{1}}{\partial
t}+{\frac{\tilde{f}^{0}}{H}}\boldsymbol{\wedge}_{h}{\tilde{F}^{n-1}_{lin}}+\mathbf{D}_{1}{B}^{0}_{lin}$
$\displaystyle=$ $\displaystyle 0,$ (106)
$\displaystyle\frac{\partial\tilde{h}^{2}}{\partial
t}+\mathbf{\bar{D}}_{2}\tilde{F}^{n-1}_{lin}$ $\displaystyle=$ $\displaystyle
0,$ (107)

where the discrete functional derivatives read

$\displaystyle\tilde{F}^{n-1}_{lin}:=\frac{\tilde{\delta}\\!\operatorname{\mathcal{H}}_{lin}}{\delta{v}^{1}}=H\tilde{u}^{n-1},\quad\quad\quad
B_{lin}:=\frac{\tilde{\delta}\\!\operatorname{\mathcal{H}}_{lin}}{\delta\tilde{h}^{2}}=g({h}^{0}+{h_{s}}^{0}),$
(108)

for the discrete linearized Hamiltonian

$\mathcal{H}_{lin}[{v}^{1},\tilde{h}^{2}]=\frac{g}{2}\left<{\tilde{h}^{2}},{\tilde{h}^{2}}\right>+g\left<{\tilde{h}^{2}},{\tilde{h_{s}}^{2}}\right>+\frac{H}{2}\left<{{u}^{1}},{{u}^{1}}\right>.$
(109)

As expected, almost all the wedge products (which give the nonlinearities)
have disappeared, other than the linearized PV wedge product
${\frac{\tilde{f}^{0}}{H}}\boldsymbol{\wedge}_{h}{\tilde{F}^{n-1}_{lin}}$.

##  6 Scheme Properties

The scheme presented in Section 5 has, by construction, several desirable
properties (see [60, 21] for more discussion of desirable properties for
schemes used in atmospheric and oceanic models). These properties come from
the operator properties discussed in Section 4 and are detailed below. They
can broadly be divided into topological properties and metric properties
(similarly to the split finite element approach taken in [6]). Topological
properties are those that depend only on the properties of the topological
operators (the exterior derivative, topological pairing and wedge product), or
in other words the Poisson bracket. The metric properties involve the Hodge
star and/or inner product, i.e. the Hamiltonian.

###  6.1 Topological Properties

Topological properties include structural properties such as no spurious
vorticity production and PV compatibility/steady geostrophic modes; and
conservation laws (total energy, mass, circulation and potential enstrophy)
which arise because the discrete Poisson bracket (102) is anti-symmetric and
has a subset of the continuous Casimirs. Fundamentally, all of these rely on
the discrete operators having discrete analogues of the key split exterior
calculus identities: integration by parts for
$\operatorname{d}$/$\left<\left<{},{}\right>\right>$, annihilation for
$\operatorname{d}$, anti-symmetry for $\wedge$ and a (partial) Leibniz rule
for $\operatorname{d}$/$\wedge$.

####  6.1.1 Total Energy Conservation

Consider $\\{\operatorname{\mathcal{B}},\operatorname{\mathcal{A}}\\}$:

$\\{\operatorname{\mathcal{B}},\operatorname{\mathcal{A}}\\}=-\left<\left<{\frac{\tilde{\delta}\\!\operatorname{\mathcal{B}}}{\delta{v}^{1}}},{{\tilde{q}^{0}}\boldsymbol{\wedge}_{h}{\frac{\tilde{\delta}\\!\operatorname{\mathcal{A}}}{\delta{v}^{1}}}}\right>\right>-\left<\left<{\frac{\tilde{\delta}\\!\operatorname{\mathcal{B}}}{\delta{v}^{1}}},{\mathbf{D}_{1}\frac{\tilde{\delta}\\!\operatorname{\mathcal{A}}}{\delta\tilde{h}^{2}}}\right>\right>-\left<\left<{\frac{\tilde{\delta}\\!\operatorname{\mathcal{B}}}{\delta\tilde{h}^{2}}},{\mathbf{\bar{D}}_{2}\frac{\tilde{\delta}\\!\operatorname{\mathcal{A}}}{\delta{v}^{1}}}\right>\right>.$
(110)

Using the anti-symmetry property (95) for ${}\boldsymbol{\wedge}_{h}{}$ the
first term can be written as

$\left<\left<{\frac{\tilde{\delta}\\!\operatorname{\mathcal{A}}}{\delta{v}^{1}}},{{\tilde{q}^{0}}\boldsymbol{\wedge}_{h}{\frac{\tilde{\delta}\\!\operatorname{\mathcal{B}}}{\delta{v}^{1}}}}\right>\right>.$
(111)

Using the integration by parts property (87) for $\mathbf{D}_{1}$ and
$\mathbf{\bar{D}}_{2}$ the last two terms can be written as

$\left<\left<{\frac{\tilde{\delta}\\!\operatorname{\mathcal{A}}}{\delta{v}^{1}}},{\mathbf{D}_{1}\frac{\tilde{\delta}\\!\operatorname{\mathcal{B}}}{\delta\tilde{h}^{2}}}\right>\right>+\left<\left<{\frac{\tilde{\delta}\\!\operatorname{\mathcal{A}}}{\delta\tilde{h}^{2}}},{\mathbf{\bar{D}}_{2}\frac{\tilde{\delta}\\!\operatorname{\mathcal{B}}}{\delta{v}^{1}}}\right>\right>.$
(112)

Therefore, combining (111) and (112) it is clear that the discrete Poisson
bracket (102) is anti-symmetric:

$\\{\operatorname{\mathcal{A}},\operatorname{\mathcal{B}}\\}=-\\{\operatorname{\mathcal{B}},\operatorname{\mathcal{A}}\\},$
(113)

and total energy is conserved since

$\displaystyle\frac{\operatorname{d}\\!\operatorname{\mathcal{H}}}{\operatorname{d}\\!t}=\\{\operatorname{\mathcal{H}},\operatorname{\mathcal{H}}\\}=-\\{\operatorname{\mathcal{H}},\operatorname{\mathcal{H}}\\}=0.$
(114)

####  6.1.2 Total Mass Conservation

Total mass is given by

$\operatorname{\mathcal{M}}=\left<{\tilde{I}^{2}},{\tilde{h}^{2}}\right>,$
(115)

which has functional derivatives

$\frac{\tilde{\delta}\\!\operatorname{\mathcal{M}}}{\delta{v}^{1}}=0,\quad\quad\quad\frac{\tilde{\delta}\\!\operatorname{\mathcal{M}}}{\delta\tilde{h}^{2}}={I}^{0}.$
(116)

Inserting these into (102) gives

$\\{\operatorname{\mathcal{M}},\operatorname{\mathcal{A}}\\}=0\quad\quad\forall\operatorname{\mathcal{A}},$
(117)

provided that $\mathbf{D}_{1}{I}^{0}=0$ i.e. (92), and thus

$\displaystyle\frac{\operatorname{d}\\!\operatorname{\mathcal{M}}}{\operatorname{d}\\!t}=\\{\operatorname{\mathcal{M}},\operatorname{\mathcal{H}}\\}=0.$
(118)

####  6.1.3 Total Circulation Conservation

Total mass-weighted potential vorticity (i.e. total circulation) is given by

$\operatorname{\mathcal{PV}}=\left<{{I}^{2}},{{\eta}^{2}}\right>\stackrel{{\scriptstyle\eqref{discrete-
topo-
pairing-1}}}{{=}}\left<\left<{\tilde{I}^{0}},{\mathbf{D}_{2}{v}^{1}}\right>\right>\stackrel{{\scriptstyle
IBP}}{{=}}-\left<\left<{\mathbf{\bar{D}}_{1}\tilde{I}^{0}},{{v}^{1}}\right>\right>\stackrel{{\scriptstyle\eqref{D-const}}}{{=}}0,$
(119)

which has functional derivatives

$\frac{\tilde{\delta}\\!\operatorname{\mathcal{PV}}}{\delta{v}^{1}}=0,\quad\quad\quad\frac{\tilde{\delta}\\!\operatorname{\mathcal{PV}}}{\delta\tilde{h}^{2}}=0.$
(120)

The zero functional derivatives are just a discrete analogue of the well-known
fact that the total circulation for a manifold without boundaries is zero, or
that the conservation of PV occurs independently of the equations of motion
(it is off-shell). Therefore

$\\{\operatorname{\mathcal{PV}},\operatorname{\mathcal{A}}\\}=0\quad\quad\forall\operatorname{\mathcal{A}},$
(121)

and

$\displaystyle\frac{\operatorname{d}\\!\operatorname{\mathcal{PV}}}{\operatorname{d}\\!t}=\\{\operatorname{\mathcal{PV}},\operatorname{\mathcal{H}}\\}=0.$
(122)

Another way to see this is to simply take $\mathbf{D}_{2}$ of the ${v}^{1}$
equation and then use the discrete Stokes theorem (93):

$\displaystyle(\tilde{I}^{0})^{T}\frac{\partial(\mathbf{D}_{2}\operatorname{\mathbf{v}})}{\partial
t}=0$ (123)

####  6.1.4 Total Potential Enstrophy Conservation

Total potential enstrophy is given by

$\operatorname{\mathcal{PE}}=\frac{1}{2}\left<{{q}^{2}},{{\eta}^{2}}\right>=\frac{1}{2}\left<{\mathbf{\bar{H}}_{0}\tilde{q}^{0}},{{\tilde{q}^{0}}\boldsymbol{\wedge}_{h}{\tilde{h}^{2}}}\right>\stackrel{{\scriptstyle\eqref{discrete-
topo-
pairing-1}}}{{=}}\frac{1}{2}\left<\left<{\tilde{q}^{0}},{{\tilde{q}^{0}}\boldsymbol{\wedge}_{h}{\tilde{h}^{2}}}\right>\right>\stackrel{{\scriptstyle\eqref{Q-antisymmetry}}}{{=}}-\frac{1}{2}\left<\left<{\tilde{h}^{2}},{{\tilde{q}^{0}}\boldsymbol{\wedge}_{h}{\tilde{q}^{0}}}\right>\right>,$
(124)

where ${q}^{2}=\mathbf{\bar{H}}_{0}\tilde{q}^{0}$, which has functional
derivatives

$\displaystyle\frac{\tilde{\delta}\\!\operatorname{\mathcal{PE}}}{\delta{v}^{1}}=-\mathbf{\bar{D}}_{1}\tilde{q}^{0},\quad\quad\frac{\tilde{\delta}\\!\operatorname{\mathcal{PE}}}{\delta\tilde{h}^{2}}=-\frac{1}{2}({\tilde{q}^{0}}\boldsymbol{\wedge}_{h}{\tilde{q}^{0}}).$
(125)

Inserting these into (102) and using the (partial) Leibniz rule II (97) along
with annihilation (86) gives

$\\{\operatorname{\mathcal{PE}},\operatorname{\mathcal{A}}\\}=0\quad\quad\forall\operatorname{\mathcal{A}},$
(126)

and therefore

$\displaystyle\frac{\operatorname{d}\\!\operatorname{\mathcal{PE}}}{\operatorname{d}\\!t}=\\{\operatorname{\mathcal{PE}},\operatorname{\mathcal{H}}\\}=0.$
(127)

It is rare that TRiSK-type schemes choose a nonlinear PV flux operator
${\tilde{q}^{0}}\boldsymbol{\wedge}_{h}{\tilde{F}^{n-1}}$ that satisfies this.
Instead, they usually opt for more sophisticated forms of PV advection with
better behaviour for sharp gradients, such as upwinding, WENO or CLUST.

####  6.1.5 No Spurious Vorticity Production

The evolution equation for absolute vorticity ${\eta}^{2}$ (obtained by
$\mathbf{D}_{2}$ applied to (105)) is:

$\frac{\partial{\eta}^{2}}{\partial
t}+\mathbf{D}_{2}({\tilde{q}^{0}}\boldsymbol{\wedge}_{h}{\tilde{F}^{1}})+\cancel{\mathbf{D}_{2}\mathbf{D}_{1}{B}^{0}}=0.$
(128)

The last term is equal to zero since $\mathbf{D}_{2}\mathbf{D}_{1}=0$, i.e.
(86), which means there is no spurious vorticity production.

####  6.1.6 PV Compatibility/Steady Geostrophic Modes

If $\tilde{q}^{0}=\tilde{I}^{0}$ (or in fact any constant), then the
${\eta}^{2}$ evolution equation becomes

$\frac{\partial({\tilde{I}^{0}}\boldsymbol{\wedge}_{h}{\tilde{h}^{2}})}{\partial
t}+\mathbf{D}_{2}({\tilde{I}^{0}}\boldsymbol{\wedge}_{h}{\tilde{F}^{1}})=0.$
(129)

Now compare this to taking ${\tilde{I}^{0}}\boldsymbol{\wedge}_{h}{\square}$
of the $\tilde{h}^{2}$ equation

$\frac{\partial({\tilde{I}^{0}}\boldsymbol{\wedge}_{h}{\tilde{h}^{2}})}{\partial
t}+{\tilde{I}^{0}}\boldsymbol{\wedge}_{h}{\mathbf{\bar{D}}_{2}\tilde{F}^{1}}=0.$
(130)

If the (partial) Leibniz rule I (96) holds, then
$\mathbf{D}_{2}({\tilde{I}^{0}}\boldsymbol{\wedge}_{h}{\tilde{F}^{1}})={\tilde{I}^{0}}\boldsymbol{\wedge}_{h}{\mathbf{\bar{D}}_{2}\tilde{F}^{1}}$
for arbitrary $\tilde{F}^{1}$, and (129) and (130) are the same equation. In
that case, we say that the scheme has PV compatibility (a uniform PV field
remains uniform). In the linearized equations, the same requirements arise for
the presence of steady geostrophic modes (see [65] for more details).

###  6.2 Metric Properties

The metric properties involve a combination of the topological operators along
with the metric operators. They are

  * •

Linear modes: spurious stationary modes are avoided through the choice of grid
staggering in the scheme (equivalent to an Arakawa C-grid) except for the
Coriolis mode [53]. The absence or presence of spurious branches of the
dispersion relationship will depend on the grid topology. For example, a grid
of quadrilaterals will have no spurious branches, a grid of triangles will
have spurious inertia-gravity waves and a grid of hexagons will have spurious
Rossby waves [50, 47, 64]. Commonly used quasi-uniform spherical grids such as
the cubed-sphere and hexagonal-pentagonal icosahedron have similar
considerations, for details see [19, 68]. The actual numerical values for the
dispersion relationship itself will be a function of the (metric dependent)
Hodge star operators and the choice of
${\tilde{I}^{0}}\boldsymbol{\wedge}_{h}{\tilde{x}^{1}}$ (the linearized
version of the operator used to compute the PV flux term).

  * •

Accuracy: the accuracy of the operators are determined by the exterior
derivative and Hodge stars for the differential operators, and by the wedge
products and Hodge stars for the product operators. This last point is
especially important, since it is usually these product operators that suffer
from insufficient accuracy (see Section 8 and [67, 23]).

  * •

Hollingsworth instability [35, 42, 48, 7]: Hollingsworth instability occurs
due to a combination of wedge products, Hodge stars and exterior derivative
operators. A general prescription for avoiding it has yet to be found, but
there are several known remedies that can alleviate it in the literature, all
based on modifying either the kinetic energy part of the Hamiltonian [35, 26]
or the nonlinear PV flux term [27]. Such modifications can be done in a way
that preserves the other desirable properties, except potential enstrophy
conservation.

###  6.3 Summary of Scheme Properties

The key features of TRiSK-type discretizations are conservation laws, no
spurious vorticity production and PV compatibility/steady geostrophic modes.
They all arise due to the operator properties of the topological operators
(exterior derivatives, topological pairing and wedge products) and are
independent of the choice of Hodge stars and grid topology and geometry. The
remaining desirable properties are a good representation of linear modes,
accuracy and (avoidance of) Hollingsworth instability, which interestingly are
also the areas where TRiSK-type schemes can have issues. These remaining
properties are largely a function of the choice of Hodge stars, specific wedge
product coefficients and the choice of grid topology and geometry. This
highlights the possibility of changing these choices to improve the desirable
metric properties while keeping the desirable topological properties intact.
This is explored more in Section 8.3.

##  7 Specific Operator Choices

To close the general DEC presented in Section 4, choices must be made for the
Hodge star and wedge product operators. Specifically, there are three choices
to be made here: a choice of Hodge star, a choice of PV wedge product, and a
choice of KE wedge product. These choices will determine the supported grid
topologies and geometries. In this section we identify possible choices from
the TRiSK, DEC and split FE literature, along with simple extensions of
existing approaches to new grid topologies or geometries.

###  7.1 TRiSK Notation

To start, we will first establish the correspondence between the notation used
in [19, 23] and the one used in this paper, to enable easier comparison with
schemes in the literature in Section 8. The former notation will be used only
in this section and the next one.

For the Hodge star, the correspondence is

$\displaystyle\operatorname{\mathbf{I}}=\mathbf{\bar{H}}_{2},\quad\quad\operatorname{\mathbf{H}}=\mathbf{H}_{1},\quad\quad\operatorname{\mathbf{J}}=\mathbf{H}_{2}.$
(131)

For the PV wedge products, the correspondence is

$\operatorname{\mathbf{R}}\tilde{x}^{2}={\tilde{I}^{0}}\boldsymbol{\wedge}_{h}{\tilde{x}^{2}},\quad\quad\operatorname{\mathbf{W}}\tilde{x}^{1}={\tilde{I}^{0}}\boldsymbol{\wedge}_{h}{\tilde{x}^{1}},\quad\quad\operatorname{\mathbf{Q}}\tilde{x}^{1}={\tilde{q}^{0}}\boldsymbol{\wedge}_{h}{\tilde{x}^{1}},$
(132)

This form of $\operatorname{\mathbf{R}}$ is written such that an explicit
formula for $\tilde{x}^{0}$ always exists, which is true for all the wedge
products discussed here. The conditions on the PV wedge product required for
energy conservation (i.e. antisymmetry (95)) can be written as

$\operatorname{\mathbf{W}}=-\operatorname{\mathbf{W}}^{T},\quad\quad\quad\operatorname{\mathbf{Q}}=-\operatorname{\mathbf{Q}}^{T},$
(133)

those for PV compatibility (partial Leibnitz rule I (96)) as

$\operatorname{\mathbf{R}}\mathbf{\bar{D}}_{2}=\mathbf{D}_{2}\operatorname{\mathbf{W}},$
(134)

and those for potential enstrophy conservation (partial Leibnitz rule II (97))
as

$\operatorname{\mathbf{Q}}\mathbf{\bar{D}}_{1}\tilde{q}^{0}-\mathbf{D}_{1}\operatorname{\mathbf{R}}^{T}\frac{\tilde{q}^{0}}{2}=0.$
(135)

In the last equation,
$\operatorname{\mathbf{R}}^{T}\frac{\tilde{q}^{0}}{2}={\tilde{q}^{0}}\boldsymbol{\wedge}_{h}^{\star}{\tilde{q}^{0}}$,
where the adjoint is from
${\tilde{x}^{0}}\boldsymbol{\wedge}_{h}{\tilde{x}^{2}}$ i.e.
$\operatorname{\mathbf{R}}$. These are the same formulas given in [65, 19,
23], up to different sign conventions.

For the kinetic energy, we will use the approach of starting with a definition
for $\tilde{K}^{2}={{x}^{1}}\boldsymbol{\wedge}_{h}{\tilde{y}^{1}}$, and then
determining the wedge products needed for the mass flux as adjoint operators.
The operator ${{x}^{1}}\boldsymbol{\wedge}_{h}{\tilde{y}^{1}}$ is closely
related to the $T$ operator in [65], but acts on ${x}^{1}$ and
$\tilde{y}^{1}$, while $T$ acts on ${x}^{1}$ and ${y}^{1}$. We do not use the
$\operatorname{\boldsymbol{\phi}^{T}}$ and $\operatorname{\boldsymbol{\phi}}$
notation from [23], since it is not general enough to encompass the two
different mass flux wedge products that arise in the case of non-diagonal
$\operatorname{\mathbf{H}}$.

###  7.2 Choices for Hodge Stars: $\operatorname{\mathbf{I}}$,
$\operatorname{\mathbf{J}}$ and $\operatorname{\mathbf{H}}$

A wide range of Hodge stars exist in the literature. All of the Hodge stars
discussed here are low-order (between 1st and 2nd order) accurate and
geometrically inflexible: they impose some sort of restriction on the topology
and/or geometry of the grids. Generally speaking, they can be categorized into
three types:

  * •

Voronoi Hodge star999also known as the circumcentric or diagonal Hodge star.
[33, 34]. This Hodge star is based on the 1-1 relationship between $k$-cells
and $(n-k)$-cells on opposite grids, and has (diagonal) stencil coefficients
equal to the ratio of areas for these two geometric entities:

$I_{\tilde{c},v}=\frac{1}{A_{\tilde{c}}},\quad\quad
J_{c,\tilde{v}}=\frac{1}{A_{c}},\quad\quad
H_{\tilde{e},e}=\frac{A_{e}}{A_{\tilde{e}}}.$ (136)

More generally, it is the ratio of measures between the target $(n-k)$-cell
and the source $k$-cell (areas for $2$-cells, lengths for $1$-cells) for
relevant geometric entities, where the measures of $0$-cells $v$ and
$\tilde{v}$ are $1$ by definition. This Hodge star requires grids with
circumcentric duality, which implies orthogonality. However, there has been
some work to modify $\operatorname{\mathbf{H}}$ to support non-orthogonal
grids, given in [76, 67]. Unfortunately, this operator is inconsistent on many
commonly used quasi-uniform spherical grids such as the cubed sphere.

  * •

Barycentric Hodge star: Using a simplicial straight grid along with a
barycentric twisted grid, a non-diagonal Hodge star can be defined based on
nearest-neighbors. Examples of this can be found in [32, 3].

  * •

Galerkin Hodge star [32]: Using the mass matrices induced by a set of basis
functions for discrete differential forms, a Hodge star operator can be
defined. This is usually done with Whitney forms [32, 63], although
alternative basis functions have been developed [10] based on energetics
definitions. These approaches require a simplicial straight grid along with a
barycentric twisted grid. An interesting extension of this idea to arbitrary
complexes is provided in [66], based on primal-dual finite elements.
Unfortunately, this Hodge star requires the solution of a linear system of
equations, and is therefore computationally challenging for TRiSK-type
schemes.

A comparison between the three types of Hodge stars for 2D surfaces can be
found in [45, 44], while connections between various Hodge star definitions
and how they make an inner product are explored in [31].

####  7.2.1 Alternative Hodge Stars

In addition to the Voronoi, barycentric and Galerkin Hodge stars, there are
some additional alternative Hodge stars that have appeared in the literature.
The $G_{ij}$ Hodge star [73] is based on the metric tensor of the underlying
coordinate system used to define the computational grid. It is defined around
a topologically square, possibly non-orthogonal grid. In the case of an
orthogonal grid it will reduce to the Voronoi Hodge star. The extension of
this Hodge star to arbitrary topologies and geometries seems possible,
although $\operatorname{\mathbf{H}}$ will likely be quite complicated. A
higher-order Hodge star on uniform rectangular grids based on coefficient
fitting can be found in [14]. The extension of this idea to arbitrary
topologies and geometries seems to be quite intricate, however, and it is
unclear if it can be generalized. Finally, a diagonal, positive-definite Hodge
star for non-orthogonal grids based on constitutive relationships can be found
in [18].

###  7.3 Choices for PV Wedge Products $\operatorname{\mathbf{R}}$,
$\operatorname{\mathbf{W}}$ and $\operatorname{\mathbf{Q}}$

The major innovation of [68, 52], and indeed the basis for all TRiSK-type
schemes, is an explicit formula for $\operatorname{\mathbf{W}}$ coefficients
given a choice of $\operatorname{\mathbf{R}}$ coefficients with a nearest-
neighbor stencil, such that all the desirable properties for
$\operatorname{\mathbf{W}}$ are obtained
($\operatorname{\mathbf{W}}=-\operatorname{\mathbf{W}}^{T}$ and
$\operatorname{\mathbf{R}}\mathbf{\bar{D}}_{2}=\mathbf{D}_{2}\operatorname{\mathbf{W}}$).
This is hereafter refered to as the TRSK2010 approach. Therefore, the choice
of PV wedge product really amounts to a choice of $\operatorname{\mathbf{R}}$,
along with a choice of how to construct $\operatorname{\mathbf{Q}}$ given
$\operatorname{\mathbf{R}}$ and $\operatorname{\mathbf{W}}$.

####  7.3.1 Choice of $\operatorname{\mathbf{R}}$/$\operatorname{\mathbf{W}}$

General nearest-neighbor stencil definitions of $\operatorname{\mathbf{R}}$
and $\operatorname{\mathbf{W}}$ are given by [32]

$\displaystyle(\operatorname{\mathbf{R}}\tilde{x}^{2})_{c}$
$\displaystyle=\sum_{\tilde{c}\in
VC(c)}R_{\tilde{c},c}\tilde{x}^{2}_{\tilde{c}},$ (137)
$\displaystyle(\operatorname{\mathbf{W}}\tilde{x}^{1})_{e}$
$\displaystyle=\sum_{\tilde{e}\in
ECP(e)}W_{\tilde{e},e}\tilde{x}^{1}_{\tilde{e}}.$ (138)

Using the TRSK2010 approach, the $\operatorname{\mathbf{W}}$ coefficients
($W_{\tilde{e},e}$) are uniquely determined from the
$\operatorname{\mathbf{R}}$ coefficients ($R_{\tilde{c},c}$). Additionally,
this choice of $\operatorname{\mathbf{R}}$ is explicitly solvable for
$\tilde{q}^{0}$. The specific choice of coefficients for
$\operatorname{\mathbf{R}}$ follows one of two approaches:

  * •

Metric [23, 65]: $R_{\tilde{c},c}$ is defined in terms of overlap areas
$A_{\tilde{c},c}$ between $c$ and $\tilde{c}$ and the twisted cell area
$A_{\tilde{c}}$ (see Figure 2) as

$\displaystyle R_{\tilde{c},c}$ $\displaystyle=$
$\displaystyle\frac{A_{\tilde{c},c}}{A_{\tilde{c}}}.$ (139)

This approach is valid for arbitrary grid geometries.

  * •

Combinatorical [32]: $R_{\tilde{c},c}$ is defined in terms of the size
$[CV(v)]$ (number of elements) of $CV(v)$ as

$\displaystyle R_{\tilde{c},c}$ $\displaystyle=$
$\displaystyle\frac{1}{[CV(v)]},$ (140)

where $v$ is the straight grid vertex dual to the dual grid cell $\tilde{c}$.
This approach is valid for arbitrary grid geometries as well, although it has
only been used, as far as we aware, on square topologies (with possible lat-
lon singularities).

These two approaches give the same results on a uniform grid. Another option
for $\operatorname{\mathbf{R}}$ and $\operatorname{\mathbf{W}}$ is to use the
primal-dual FE operators from [66], since they share the same stencil as the
DEC ones and have the same key properties.

####  7.3.2 Choice of $\operatorname{\mathbf{Q}}$

Given $\operatorname{\mathbf{W}}$ and $\operatorname{\mathbf{R}}$, there are
four (non-exhaustive) possible definitions of $\operatorname{\mathbf{Q}}$,
differing in whether they conserve total energy
($\operatorname{\mathbf{Q}}^{TE}$), potential enstrophy
($\operatorname{\mathbf{Q}}^{PE}$), both ($\operatorname{\mathbf{Q}}^{DBL}$)
or none ($\operatorname{\mathbf{Q}}^{accur}$):

  * •

$\operatorname{\mathbf{Q}}^{TE}$ satifies
$\operatorname{\mathbf{Q}}=-\operatorname{\mathbf{Q}}^{T}$, and therefore
conserves total energy (as shown in Section 6.1.1):

$\displaystyle(\operatorname{\mathbf{Q}}^{TE}\tilde{x}^{1})_{e}$
$\displaystyle=\sum_{\tilde{e}\in
ECP(e)}\frac{q^{e}+q^{\tilde{e}}}{2}W_{\tilde{e},e}\tilde{x}^{1}_{\tilde{e}},$
(141)

for an arbitrary dual edge reconstruction $q^{e}$. This definition is
equivalent to

$\operatorname{\mathbf{Q}}^{TE}=\frac{1}{2}(q^{e}\operatorname{\mathbf{W}}+\operatorname{\mathbf{W}}q^{e}),$

It is clear from this latter form of the definition that
$\operatorname{\mathbf{Q}}=-\operatorname{\mathbf{Q}}^{T}$, provided
$\operatorname{\mathbf{W}}=-\operatorname{\mathbf{W}}^{T}$. Oftentimes,
sophisticated advection schemes (APVM, LUST, CLUST, upwinding, WENO, etc.) are
used to compute $q^{e}$ [75, 77], especially on hexagonal grids to control the
spurious Rossby mode. A clever choice of $q^{e}$ was developed in [27] to help
alleviate Hollingsworth instability.

  * •

$\operatorname{\mathbf{Q}}^{PE}$ satisfies
$\operatorname{\mathbf{Q}}\mathbf{\bar{D}}_{1}\tilde{q}^{0}-\mathbf{D}_{1}\operatorname{\mathbf{R}}^{T}\frac{(\tilde{q}^{0})^{2}}{2}=0$,
and therefore conserves potential enstrophy (as shown in Section 6.1.4):

$\displaystyle(\operatorname{\mathbf{Q}}^{PE}\tilde{x}^{1})_{e}$
$\displaystyle=\sum_{\tilde{e}\in ECP(e)}\sum_{\tilde{v}\in
VE(\tilde{e})}\frac{\tilde{q}^{0}_{\tilde{v}}}{2}W_{\tilde{e},e}\tilde{x}^{1}_{\tilde{e}}.$
(142)

  * •

$\operatorname{\mathbf{Q}}^{DBL}$ satisfies both
$\operatorname{\mathbf{Q}}=-\operatorname{\mathbf{Q}}^{T}$ and
$\operatorname{\mathbf{Q}}\mathbf{\bar{D}}_{1}\tilde{q}^{0}-\mathbf{D}_{1}\operatorname{\mathbf{R}}^{T}\frac{(\tilde{q}^{0})^{2}}{2}=0$,
and therefore conserves total energy and potential enstrophy:

$\displaystyle(\operatorname{\mathbf{Q}}^{DBL}\tilde{x}^{1})_{e}$
$\displaystyle=\sum_{\tilde{e}\in ECP(e)}\sum_{\tilde{v}\in
VE(\tilde{e})}Q_{\tilde{e},e,\tilde{v}}\tilde{x}^{1}_{\tilde{e}}\tilde{q}^{0}_{\tilde{v}}.$
(143)

The $Q_{\tilde{e},e,\tilde{v}}$ coefficients in
$\operatorname{\mathbf{Q}}^{DBL}$ are uniquely determined by the choice of
$\operatorname{\mathbf{R}}$ coefficients. The first example of this appeared
in [2], with an extension to higher-order accuracy in the implied vorticity
advection equation in [62]. A general framework that combines [2] and [62]
appears in [55], which identifies other version of $\operatorname{\mathbf{Q}}$
with the same doubly conservative properties. All of these approaches are
restricted to logically square topologies. An extension of [55] to arbitrary
grid topologies and geometries is found in [23]. However, only metric
$\operatorname{\mathbf{R}}$ coefficients are evaluated in that work. A
$\operatorname{\mathbf{Q}}$ with partial double conservation is found in [54]:
it conserves total energy but only enstrophy in the non-divergent limit, not
potential enstrophy.

  * •

$\operatorname{\mathbf{Q}}^{accur}$ gives up conservation in exchange for
accurate advection of PV:

$\displaystyle(\operatorname{\mathbf{Q}}^{accur}\tilde{x}^{1})_{e}$
$\displaystyle=\sum_{\tilde{e}\in
ECP(e)}q^{e}W_{\tilde{e},e}\tilde{x}^{1}_{\tilde{e}}.$ (144)

recalling $q^{e}$ is an arbitrary dual edge reconstruction, usually obtained
from a sophisticated advection scheme such as APVM, LUST, CLUST, upwinding,
WENO, etc.

For all four definitions, if $\tilde{q}^{0}=\tilde{I}^{0}$, then
$\operatorname{\mathbf{Q}}$ becomes equal to the $\operatorname{\mathbf{W}}$
defined from $\operatorname{\mathbf{R}}$. This is clear for
$\operatorname{\mathbf{Q}}^{TE}$, $\operatorname{\mathbf{Q}}^{PE}$ and
$\operatorname{\mathbf{Q}}^{accur}$ since they are based on the
$\operatorname{\mathbf{W}}$ coefficients and $q^{e}=1$ when
$\tilde{q}^{0}=\tilde{I}^{0}$. For $\operatorname{\mathbf{Q}}^{DBL}$ see the
proof in Appendix B that enforcing both total energy and potential enstrophy
conservation automatically gives PV compatibility.

###  7.4 Choices for KE Wedge Product

Unlike the PV wedge product, the KE wedge product does not have to satisfy any
properties. Therefore, there is some additional freedom in how it is defined.
A general nearest-neighbor stencil definition of
${{x}^{1}}\boldsymbol{\wedge}_{h}{\tilde{y}^{1}}$ is given by

$({{x}^{1}}\boldsymbol{\wedge}_{h}{\tilde{y}^{1}})_{\tilde{c}}=\sum_{\tilde{e}\in
EC(\tilde{c})}\mathfrak{T}_{\tilde{c},e,\tilde{e}}{x}^{1}_{e}\tilde{y}^{1}_{\tilde{e}},$
(145)

where $e$ is the straight grid edge corresponding to the twisted grid edge
$\tilde{e}$. As for $\operatorname{\mathbf{R}}$, either metric or
combinatorial definitions can be used for
$\mathfrak{T}_{\tilde{c},e,\tilde{e}}$:

  * •

Metric [52, 67]: $\mathfrak{T}_{\tilde{c},e,\tilde{e}}$ is defined in terms of
overlap areas $A_{\tilde{c},\tilde{e}}$ and extended edge area
$\tilde{A}_{\tilde{e}}$ as

$\mathfrak{T}_{\tilde{c},e,\tilde{e}}=\frac{A_{\tilde{c},\tilde{e}}}{\tilde{A}_{\tilde{e}}}.$
(146)

This approach is valid for arbitrary grid geometries.

  * •

Combinatorical [32]: $\mathfrak{T}_{\tilde{c},e,\tilde{e}}$ is defined in
terms of the size $[CE(\tilde{e})]$ of $CE(\tilde{e})$, which is always
$\frac{1}{2}$ on a conforming grid without boundaries, i.e.

$\displaystyle\mathfrak{T}_{\tilde{c},e,\tilde{e}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}.$ (147)

This approach is valid for arbitrary grid geometries as well, although it has
only been used on square topologies (with possible lat-lon singularities).

From the coefficients $\mathfrak{T}_{\tilde{c},e,\tilde{e}}$, it is easy to
define the adjoint wedge products
${{x}^{0}}\boldsymbol{\wedge}_{h}{\tilde{y}^{1}}$ and
${{x}^{0}}\boldsymbol{\wedge}_{h}{{y}^{1}}$ needed to compute the mass flux
$\tilde{F}^{n-1}$:

$\displaystyle({{h}^{0}}\boldsymbol{\wedge}_{h}^{\star}{{u}^{1}})_{e}$
$\displaystyle=$ $\displaystyle\sum_{\tilde{c}\in
CE(e)}\mathfrak{T}_{\tilde{c},e,\tilde{e}}{h}^{0}_{v}{u}^{1}_{e},$ (148)
$\displaystyle({{h}^{0}}\boldsymbol{\wedge}_{h}^{\star}{\tilde{u}^{n-1}})_{\tilde{e}}$
$\displaystyle=$ $\displaystyle\sum_{\tilde{c}\in
CE(\tilde{e})}\mathfrak{T}_{\tilde{c},e,\tilde{e}}{h}^{0}_{v}\tilde{u}^{n-1}_{\tilde{e}}\
.$ (149)

In the case of a uniform grid, the metric and combinatorical approaches give
the same results.

##### Alternative Choices

Clever choices of KE wedge products with slightly different stencils were
developed in [35, 26] to help alleviate Hollingsworth instability. Another
alternative KE wedge product with improved accuracy was developed in [13],
based on a least square reconstruction, although it was used only to compute
kinetic energy and not the corresponding mass flux. As for
$\operatorname{\mathbf{R}}$ and $\operatorname{\mathbf{W}}$, it would also be
possible to use the primal-dual FE approach from [66] to define a KE wedge
product, although this has not been worked out yet. A final view of the wedge
product is through the connection with the interior product (contraction),
given by Hirani’s formula
$\mathrm{i}_{\mathbf{x}}\alpha=(-1)^{k(n-k)}\operatorname{\tilde{\star}}(\operatorname{\tilde{\star}}\alpha\wedge\mathbf{x}^{\flat})$.
A discrete version of this is explored in [8].

###  7.5 Choice of grid

As discussed in Section 4.1, grids have two aspects: topology and geometry.
Within TRiSK-type schemes, there is a distinction between schemes that are
capable of handling only topologically square grids (along with possibly a
lat-lon grid singularity), and those capable of handling arbitrary topologies.
Additionally, there is a distinction between schemes that can handle arbitrary
(non-orthogonal) grid geometries, and those that can handle only orthogonal
grid geometries. For quasi-uniform grids such as the icosahedral-pentagonal
and cubed sphere, the ability to handle general conforming topologies (for
both) and non-orthogonal geometries (for cubed sphere) is required. These
restrictions on grid topology/geometry arise due to the choices made for the
Hodge star and wedge product operators, as discussed in the previous
subsections.

###  7.6 Grid Optimization

In addition to these choices of operators, it is also possible to apply many
different optimization schemes to improve the geometry of the grid. This
includes tweaking [28, 29], centroidal voronoi tesselations [43, 38, 15],
spring-dynamics [69, 37], generalized power grids [24] and Hodge-optimized
triangulations [46]. Some of these optimization choices were compared in [19,
23, 24]; but only for certain choices of operators. Grid optimization is
primarily utilized to improve grid regularity, and therefore improve operator
and scheme accuracy, and to allow a larger explicit time step.

###  7.7 Effects of choices on properties

In general, the topological properties are obtained by any of the choices
discussed above, with the possible exception of total energy or potential
enstrophy conservation depending on the choice of $\operatorname{\mathbf{Q}}$.
The choice of operators does affect the metric properties, specifically:

  * •

The Hodge star choice (specifically $\mathbf{D}_{k}$/$\mathbf{\bar{D}}_{k}$
plus $\operatorname{\mathbf{I}}$, $\operatorname{\mathbf{J}}$ and
$\operatorname{\mathbf{H}}$) determines the accuracy of the differential
operators, and the combination of the Hodge star and the wedge product
(specifically
$\operatorname{\mathbf{R}}$/$\operatorname{\mathbf{W}}$/$\operatorname{\mathbf{Q}}$
or ${{x}^{1}}\boldsymbol{\wedge}_{h}{\tilde{y}^{1}}$ plus
$\operatorname{\mathbf{I}}$, $\operatorname{\mathbf{J}}$ and
$\operatorname{\mathbf{H}}$) determines the accuracy of the scalar/vector
product operators.

  * •

Hollingsworth instability is controlled by the choice of wedge products used
in the kinetic energy/mass flux
(${{x}^{1}}\boldsymbol{\wedge}_{h}{\tilde{y}^{1}}$) and/or the PV flux term
($\operatorname{\mathbf{Q}}$), as well as Hodge stars
$\operatorname{\mathbf{H}}$ and $\operatorname{\mathbf{I}}$.

  * •

The choice of grid topology determines the presence or absence of spurious
branches of dispersion relationship.

  * •

The choice of grid optimization can affect the accuracy of the differential
operators.

The main issues seen with TRiSK-type schemes are the inability to find (on
arbitrary grids) accurate and Hollingsworth free discrete versions of the
operators. However, as discussed in Section 8.3, only a small portion of the
possible choices have been explored. With careful grid optimization it seems
possible to, at least, get 1st order in $L_{\infty}$ for some tested quasi-
uniform spherical grids for differential operators, but scalar/vector products
are trickier. The combinatorial wedge products do not have accuracy issues,
including lat-lon variants and the extension [73] to topologically square non-
orthogonal grids. However, these wedge products have not been tested on
arbitrary grids. The metric wedge products
($\operatorname{\mathbf{R}}$/$\operatorname{\mathbf{W}}$/$\operatorname{\mathbf{Q}}$
and ${{x}^{1}}\boldsymbol{\wedge}_{h}{\tilde{y}^{1}}$) on arbitrary grids do
have accuracy issues in at least $L_{\infty}$. Hollingsworth instability
appears avoidable (or at least controllable) through careful specification of
either $q^{e}$ using in $\operatorname{\mathbf{Q}}^{TE}$ [27] and/or
${{x}^{1}}\boldsymbol{\wedge}_{h}{\tilde{y}^{1}}$ [35, 26].

##  8 Schemes from the literature

In this section, we identify the operator and grid optimization choices that
were made101010Usually implicitly, unaware of the DEC structure. in a variety
of TRiSK-type schemes from the literature. Specifically, we will look at the
schemes in [2] (AL81), [54] (SD75), [62] (TW84), [52] (TRSK2010), [67]
(Thuburn2014) , [75, 77, 76] (Weller2014), [19, 23] (Eldred2017) [55] (S04)
and [73] (Toy2017). Some of this connection was covered in [19, 23], but the
understanding of the discrete wedge product was lacking. A detailed summary of
the various choices made by TRiSK-type schemes in the literature is given in
Table 2.

Scheme | Grid |  $\operatorname{\mathbf{I}}$/$\operatorname{\mathbf{J}}$/$\operatorname{\mathbf{H}}$ |  $\operatorname{\mathbf{R}}$/$\operatorname{\mathbf{W}}$ | $\operatorname{\mathbf{Q}}$ | ${{x}^{1}}\boldsymbol{\wedge}_{h}{\tilde{y}^{1}}$ | Dof Scaling | Grid Opt
---|---|---|---|---|---|---|---
AL81 | S-O/LL | V | C | DBL | C | yes | N/A
SD75 | S-O/LL | V | C | DBL1 | C | yes | N/A
TW84 | S-O/LL | V | C | DBL3 | C | yes | N/A
TRSK2010 | US-O | V | M | TE/PE | M | yes | CVT
Thuburn2014 | US-NO | V2 | M | ACCUR4 | other4 | no | CVT
Weller2014 | US-NO | V2 | M | TE | other5 | partial5 | CVT
Eldred2017 | US-NO | V2 | M | DBL/TE/PE | M | no | CVT, SD, TW
S04 | S-O | V | C | DBL3 | C | yes | N/A
Toy2017 | S-NO | $G_{ij}$ | C | DBL | C | yes | N/A
Table 2: A summary of the TRiSK-type schemes in the literature: the grid
topologies and geometries they support; the choice of Hodge stars
$\operatorname{\mathbf{I}}$/$\operatorname{\mathbf{J}}$/$\operatorname{\mathbf{H}}$,
PV wedge products $\operatorname{\mathbf{R}}$/$\operatorname{\mathbf{W}}$ and
$\operatorname{\mathbf{Q}}$, KE wedge product
${{x}^{1}}\boldsymbol{\wedge}_{h}{\tilde{y}^{1}}$ they make; whether dof
scaling is required to put them into the DEC framework used in this paper; and
the choice of grid optimization. Here S-O indicates square topology with
orthogonal geometry, LL indicates lat-lon topology and geometry, US-O
indicates an unstructured topology with orthogonal geometry, US-NO indicates
an unstructured topology with non-orthogonal geometry and S-NO indicates a
square topology with non-orthogonal geometry. A V indicates a Voronoi Hodge
star, while $G_{ij}$ indicates the $G_{ij}$ based Hodge star. A C indicates a
combinatorial discretization for the PV or KE wedge product, while an M
indicates a metric one. For $\operatorname{\mathbf{Q}}$, DBL indicates the
doubly conservative version $\operatorname{\mathbf{Q}}^{DBL}$, TE indicates
the total energy version $\operatorname{\mathbf{Q}}^{TE}$, PE indicates the
potential enstrophy conserving version $\operatorname{\mathbf{Q}}^{PE}$ and
ACCUR indicates the accurate PV advection version
$\operatorname{\mathbf{Q}}^{accur}$. CVT indicates centroidal Voronoi
tessellation grid optimization, SD indicates spring dynamics grid optimization
and TW indicates the tweaking grid optimization.
1 The SD75 scheme only conserves total energy plus enstrophy in the divergent
limit, not potential enstrophy.
2 The Thuburn2014 and Eldred2017 schemes use a modified
$\operatorname{\mathbf{H}}$ for non-orthogonal grids, that reduces to the
diagonal Voronoi $\operatorname{\mathbf{H}}$ for an orthogonal grid.
Weller2014 explores the same modified $\operatorname{\mathbf{H}}$, but also
introduces non-symmetric $\operatorname{\mathbf{H}}$ choices that do not
conserve energy.
3 The TW84 and S04 schemes are 4th-order for the implied vorticity equation
(only optionally in the case of S04).
4 The Thuburn2014 scheme gives up energy conservation in exchange for accurate
nonlinear advection (through choice of $\operatorname{\mathbf{Q}}$ and mass
flux), and computes kinetic energy using a least-squares approach.