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

Christopher Eldred¹¹1Center for Computing Research, Sandia National Laboratories, celdred@sandia.gov, Werner Bauer²²2University of Surrey, Department of Mathematics, w.bauer@surrey.ac.uk

(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.

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 staggering³³3Fluid 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 grids⁴⁴4and 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)⁵⁵5Discrete 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.⁶⁶6More 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).⁷⁷7A 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.

Refer to caption — 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].

4.2 Discrete Differential Forms

A discrete differential $k$ -form ( ${x}^{k}$ or $\tilde{x}^{k}$ )⁸⁸8Note 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.

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 star⁹⁹9also 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 made¹⁰¹⁰10Usually 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	DBL¹	C	yes	N/A
TW84	S-O/LL	V	C	DBL³	C	yes	N/A
TRSK2010	US-O	V	M	TE/PE	M	yes	CVT
Thuburn2014	US-NO	V²	M	ACCUR⁴	other⁴	no	CVT
Weller2014	US-NO	V²	M	TE	other⁵	partial⁵	CVT
Eldred2017	US-NO	V²	M	DBL/TE/PE	M	no	CVT, SD, TW
S04	S-O	V	C	DBL³	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.
¹ The SD75 scheme only conserves total energy plus enstrophy in the divergent limit, not potential enstrophy.
² 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.
³ The TW84 and S04 schemes are 4th-order for the implied vorticity equation (only optionally in the case of S04).
⁴ 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.
⁵ The Weller2014 scheme computes kinetic energy using a metric wedge product, but uses an accurate advection scheme for mass flux

\tilde{F}^{n-1}

; and requires scaling of fluid height but not velocity to fit into the DEC framework.

8.1 Degree of freedom scaling

Many TRiSK-type schemes in the literature work with point values, instead of the integral quantities $\tilde{h}^{2}$ and ${v}^{1}$ . Therefore, for these schemes to be put into the same form as the general DEC framework the degrees of freedom must be scaled. This involves multiplying the predicted pointwise height by the area of a twisted $2$ -cell $A_{\tilde{c}}$ to obtain $\tilde{h}^{2}$ , and the predicted pointwise velocity by the length of a straight $1$ -cell (edge) $A_{e}$ to obtain ${v}^{1}$ .

8.2 Specific schemes

AL81

The AL81 scheme [2] uses a Voronoi Hodge star, combinatorical PV and KE wedge products and $\operatorname{\mathbf{Q}}^{DBL}$ . It was developed for square, orthogonal grids, including the modifications needed to handle polar singularities on a latitude-longitude grid. A modification of the KE wedge product to suppress Hollingsworth instability can be found in [35, 42].

SD75

The SD75 scheme [54] uses a Voronoi Hodge star, combinatorical PV and KE wedge products and a slightly modified $\operatorname{\mathbf{Q}}^{DBL}$ that only conserves total energy and enstrophy in the non-divergent limit, not potential enstrophy. It was developed for square, orthogonal grids, including the modifications needed to handle polar singularities on a latitute-longitude grid.

TW84

The TW84 scheme [62] extended the AL81 scheme with a $\operatorname{\mathbf{Q}}^{DBL}$ that is 4th order for the implied vorticity equation. It is otherwise the same, including the choice of grids.

S04

The S04 scheme [55] generalizes AL81 and TW84 to a family of $\operatorname{\mathbf{Q}}^{DBL}$ (optionally extended to 4th order for the implied vorticity equation) operators, through a quasi-Hamiltonian approach. It therefore also uses a Voronoi Hodge star and combinatorical PV and KE wedge products. It was developed for square, orthogonal grids.

TRSK2010

The TRSK2010 scheme [52] uses a Voronoi Hodge star, metric PV and KE wedge products, and either $\operatorname{\mathbf{Q}}^{TE}$ or $\operatorname{\mathbf{Q}}^{PE}$ . It is designed for unstructured, orthogonal grids, primarily Delauney triangulations and their Voronoi duals such as the icosahedral-hexagonal grid. These grids are usually optimized using the centroidal Voronoi tesselation approach. TRSK2010 can be viewed as an extension of AL81 to unstructured grids with only partial conservation. For $\operatorname{\mathbf{Q}}^{TE}$ , a variety of choices for $q^{e}$ were explored, focusing mainly on the anticipated potential vorticity method (APVM). A modification of the KE wedge product to suppress Hollingsworth instability [26], and clever choice of $q^{e}$ to do the same for $\operatorname{\mathbf{Q}}^{TE}$ can be found in [27].

Thuburn2014

The Thuburn2014 scheme [67] uses a Voronoi Hodge star (with a possible modification of $\operatorname{\mathbf{H}}$ to support non-orthogonal grids) and metric PV wedge products. However, it gives up conservation of total energy and potential enstrophy for a more accurate computation of kinetic energy and better nonlinear advection properties. Specifically, the kinetic energy is computed using a least-squares approach; and accurate advection schemes are used to compute $q^{e}$ and the mass flux, using $\operatorname{\mathbf{Q}}^{accur}$ . It is designed for unstructured, non-orthogonal grids, including the icosahedral-hexagonal and cubed-sphere grids. The grid is optimized using the centroidal Voronoi tesselation approach or a related approach on cubed-sphere grids. This scheme can be thought of as an extension of TRSK2010 to non-orthogonal grids that gives up energy conservation for accurate kinetic energy and better advection.

Weller2014

The Weller2014 scheme [75, 77, 76] uses metric PV wedge products and $\operatorname{\mathbf{Q}}^{TE}$ . Voronoi Hodge stars are used for $\operatorname{\mathbf{I}}$ and $\operatorname{\mathbf{J}}$ , and several options are explored for $\operatorname{\mathbf{H}}$ : the modified $\operatorname{\mathbf{H}}$ from [67], and a non-symmetric $\operatorname{\mathbf{H}}$ that reduces to the Voronoi $\operatorname{\mathbf{H}}$ in the case of an orthogonal grid. A non-symmetric $\operatorname{\mathbf{H}}$ leads to a loss of energy conservation. The kinetic energy is computed using a metric wedge product, while the mass flux $\tilde{F}^{n-1}$ is computed using an accurate advection scheme. This mass flux also gives up energy conservation. In $\operatorname{\mathbf{Q}}^{TE}$ , $q^{e}$ is computed for a variety of different advection schemes. It is designed for unstructured, non-orthogonal grids, including the icosahedral-hexagonal and cubed-sphere grids. The grid is optimized used the centroidal Voronoi tesselation approach or a related approach on cubed-sphere grids. This scheme can be thought of as an extension of TRSK2010 to non-orthogonal grids that gives up energy conservation for better advection.

Eldred2017

The Eldred2017 scheme [19, 23] uses a Voronoi Hodge star (with a possible modification of $\operatorname{\mathbf{H}}$ to support non-orthogonal grids from [67]), metric PV and KE wedge products and any of $\operatorname{\mathbf{Q}}^{TE}$ , $\operatorname{\mathbf{Q}}^{PE}$ or $\operatorname{\mathbf{Q}}^{DBL}$ . It is essentially a merger of S04 and TRSK2010, with the main novelty being an extension of $\operatorname{\mathbf{Q}}^{DBL}$ to arbitrary grids. It is designed for unstructured, non-orthogonal grids, including the icosahedral-hexagonal and cubed-sphere grids. The grid is optimized using the centroidal Voronoi tesselation, spring dynamics or tweaking approaches.

Toy2017

The Toy2017 scheme [73] uses a $G_{ij}$ Hodge star, combinatorical PV and KE wedge products and $\operatorname{\mathbf{Q}}^{DBL}$ . It is designed for a square, non-orthogonal grid; and primarily tested for smooth deformations of a uniform square grid. The new Hodge star reduces to a Voronoi Hodge star for orthogonal grids, and therefore the scheme can be thought of as an extension of AL81 to non-orthogonal grids.

8.3 Unexplored approaches

This (brief) review of existing schemes immediately highlights some unexplored operator choices that fit within the DEC framework. Specifically, for the Hodge stars:

•

barycentric Hodge stars [10] (straightforward, already developed);
•

$G_{ij}$ Hodge star from [73] on arbitrary grids (tricky, not yet developed);
•

constitutive relationship based Hodge star from [18] (straightforward, already developed);

for the PV wedge product (specifically $\operatorname{\mathbf{R}}$ and $\operatorname{\mathbf{W}}$ ):

•

combinatorical wedge products on arbitrary grids (straightforward, already developed),
•

primal-dual FE wedge product operators from [66] (straightforward, already developed);

and for the KE wedge product:

•

combinatorical wedge products on arbitrary grids (straightforward, already developed),
•

primal-dual FE wedge product operators from [66] (straightforward, not yet developed),
•

least squares reconstruction wedge product from [13] (straightforward, already developed),
•

extrusion/contraction based wedge product from [8] (straightforward, already developed).

Most of these options have already been developed, just not studied in the context of TRiSK-type schemes. Only the primal-dual FE versions of the KE wedge product (straightforward) and the $G_{ij}$ Hodge star on arbitrary grids (tricky) would require further research. These alternative choices do not expand the stencil of the general scheme beyond nearest-neighbor (although the Hodge star is no longer diagonal) and keep the topological properties of the scheme.

These alternative choices of operators are interesting for several reasons:

•

Accuracy results suggest that the combinatorical approach to wedge products might be superior, although this is yet to be verified for quasi-uniform spherical grids. Specifically, [73] uses combinatorical wedge products and is accurate on highly-distorted non-orthogonal planar grids, while [23] showed that metric wedge products are inaccurate on quasi-uniform spherical grids. A combinatorical definition also fits with differential geometry: the wedge product is topological.
•

Barycentric Hodge stars expand the possible grid geometries, which is especially important for grids with complicated lateral boundaries or highly deformed grid cells.
•

The primal-dual FE scheme operators have proven accuracy [66], albeit when used in scheme with mass matrices.
•

The least squares based wedge product has proven accuracy [13], albeit when used in a non-mimetic scheme.

Therefore, it seems possible that some combination of these choices for operators and grid optimizations could obtain acceptable accuracy on arbitrary grids and hopefully avoid the Hollingsworth instability.

9 Conclusions

In this paper, we have shown that TRiSK-type schemes are best understood as a discrete exterior calculus applied to a Hamiltonian formulation based on split exterior calculus. The key to this was the wedge product and the topological pairing, with the latter newly introduced here in this manuscript. This new understanding gives a clear separation between the topological Poisson brackets, written entirely in terms of the wedge product, exterior derivative and topological pairing; and the metric Hamiltonian, which contains, in addition, the Hodge star and inner product. This separation is somewhat mimicked in the discrete setting, through the use of combinatorial operators for the topological parts and metric operators for the metric parts. However, the separation can be broken in the case of the wedge product, where metric operators are sometimes used. All known TRiSK-type schemes can be written in terms of this framework, differing only through their choice of Hodge stars and wedge products, along with the grids supported by these choices.

This new framework has several immediate possible applications. The first is a comprehensive investigation of accuracy and Hollingsworth instability for both explored and unexplored combinations from Section 8.3, on a variety of grids with various optimization choices. The main goal would be to identify a combination of grid choice/optimization, wedge products and Hodge stars that gives acceptable accuracy and avoids Hollingsworth instability. Such a scheme could serve as a revised basis for the operational atmospheric and oceanic models using TRiSK-type schemes that offers improved accuracy while keeping the desirable properties and nearest-neighbor stencil (and thus computational efficiency).

A second application is exploring TRiSK-type schemes with lateral boundaries, which are needed for ocean models and limited-area atmospheric models. Existing treatment of boundaries in DEC (and therefore TRiSK-type schemes) is somewhat ad-hoc, and lacks the generality to treat arbitrary boundary conditions in a consistent way. In recent work [20], a DEC with consistent treatment of boundaries and arbitrary boundary conditions was developed. The application of this DEC to TRiSK-type schemes will be reported in a future publication.

A final application is studying whether the (nonlinear) conservation properties of total energy and potential enstrophy can be kept while introducing accurate advection schemes for the fluid height and the potential vorticity, or in other words, can the loss of conservation properties in Thuburn2014 and Weller2014 be remedied in a way that still keeps accurate advection? A positive answer to this is given in [22], where a variant of TRiSK with WENO/FCT advection (that also has high-order Hodge stars) is developed for uniform rectangular grids. This does require expanding the stencil of the advection operators and the Hodge stars, but no more than required by the advection operators in the Weller2014 or Thuburn2014 schemes. Work is currently ongoing to extend this idea to unstructured grids.

10 Acknowledgements

This article has been co-authored by an employee of National Technology & Engineering Solutions of Sandia, LLC under Contract No. DE-NA0003525 with the U.S. Department of Energy (DOE). The employee owns right, title and interest in and to the article and is responsible for its contents. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this article or allow others to do so, for United States Government purposes. The DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan https://www.energy.gov/downloads/doe-public-access-plan.

References

[1] R. Abraham, J. E. Marsden, and T. Ratiu. Manifolds, tensor analysis, and applications, volume 75. Springer Science & Business Media, 2012.
[2] A. Arakawa and V. R. Lamb. A potential enstrophy and energy conserving scheme for the shallow water equations. Monthly Weather Review, 109(1):18–36, 1981.
[3] B. Auchmann and S. Kurz. A geometrically defined discrete hodge operator on simplicial cells. IEEE Transactions on Magnetics, 42(4):643–646, 2006.
[4] W. Bauer. A new hierarchically-structured n-dimensional covariant form of rotating equations of geophysical fluid dynamics. GEM - International Journal on Geomathematics, 7(1):31–101, Apr 2016.
[5] W. Bauer and J. Behrens. A structure-preserving split finite element discretization of the split wave equations. Applied Mathematics and Computation, 325:375 – 400, 2018.
[6] W. Bauer, J. Behrens, and C. J. Cotter. A structure-preserving approximation of the discrete split rotating shallow water equations. In F. J. Vermolen and C. Vuik, editors, Numerical Mathematics and Advanced Applications ENUMATH 2019, pages 103–113, Cham, 2021. Springer International Publishing.
[7] M. J. Bell, P. S. Peixoto, and J. Thuburn. Numerical instabilities of vector-invariant momentum equations on rectangular c-grids. Quarterly Journal of the Royal Meteorological Society, 143(702):563–581, 2017.
[8] A. Bossavit. Extrusion, contraction: their discretization via whitney forms. Compel-the International Journal for Computation and Mathematics in Electrical and Electronic Engineering, 22:470–480, 2003.
[9] S. Calandrini, D. Engwirda, and M. Petersen. Comparing numerical accuracy and stability for different horizontal discretizations in MPAS-Ocean. Ocean Modelling, 168:101908, 2021.
[10] Codecasa, Lorenzo, Specogna, Ruben, and Trevisan, Francesco. Geometrically defined basis functions for polyhedral elements with applications to computational electromagnetics. ESAIM: M2AN, 50(3):677–698, 2016.
[11] D. Cohen and E. Hairer. Linear energy-preserving integrators for poisson systems. BIT Numerical Mathematics, 51(1):91–101, Mar 2011.
[12] C. Cotter and J. Thuburn. A finite element exterior calculus framework for the rotating shallow-water equations. Journal of Computational Physics, 257(Part B):1506 – 1526, 2014. Physics-compatible numerical methods.
[13] Darren Engwirda. personal communication.
[14] C. Deimert, M. Potter, and M. Okoniewski. Collocated electrodynamic fdtd schemes using overlapping yee grids and higher-order hodge duals. Journal of Computational Physics, 326:629–649, 2016.
[15] Q. Du, M. Gunzburger, and L. Ju. Advances in studies and applications of centroidal voronoi tessellations. Numerical Mathematics: Theory, Methods and Applications, 3(2):119–142, 2010.
[16] T. Dubos, S. Dubey, M. Tort, R. Mittal, Y. Meurdesoif, and F. Hourdin. Dynamico-1.0, an icosahedral hydrostatic dynamical core designed for consistency and versatility. Geoscientific Model Development, 8(10):3131–3150, 2015.
[17] N. Ducousso, J. Le Sommer, J.-M. Molines, and M. Bell. Impact of the “symmetric instability of the computational kind” at mesoscale- and submesoscale-permitting resolutions. Ocean Modelling, 120:18–26, 2017.
[18] A. F. El Ouafdi, H. El Houari, and D. Ziou. Adaptive estimation of hodge star operator on simplicial surfaces. The Visual Computer, 37(6):1433–1445, 2021.
[19] C. Eldred. Linear and Nonlinear Properties of Numerical Methods for the Rotating Shallow Water Equations. PhD thesis, Colorado State University, 2015.
[20] C. Eldred. Structure-preserving numerical discretizations for domains with boundaries. Technical Report SAND2021-11517, Sandia National Laboratories, 2021.
[21] C. Eldred, T. Dubos, and E. Kritsikis. A quasi-Hamiltonian discretization of the thermal shallow water equations. Journal of Computational Physics, 2018.
[22] C. Eldred, M. Norman, and M. Taylor. PAM dynamical core- continuous formulation, discrete numerics and implementation. Technical Report WBS 2.2.3.05, Milestone ECP-AD-SE-15-1675, Sandia National Laboratories, 2021.
[23] C. Eldred and D. Randall. Total energy and potential enstrophy conserving schemes for the shallow water equations using hamiltonian methods – part 1: Derivation and properties. Geoscientific Model Development, 10(2):791–810, 2017.
[24] D. Engwirda. Generalised primal-dual grids for unstructured co-volume schemes. Journal of Computational Physics, 375:155–176, 2018.
[25] H. Flanders. Differential forms with applications to the physical sciences, volume 11. Courier Corporation, 1963.
[26] A. Gassmann. A global hexagonal C-grid non-hydrostatic dynamical core (ICON-IAP) designed for energetic consistency. Quarterly Journal of the Royal Meteorological Society, 139(670):152–175, 2013.
[27] A. Gassmann. Discretization of generalized Coriolis and friction terms on the deformed hexagonal C-grid. Quarterly Journal of the Royal Meteorological Society, 144(716):2038–2053, 2018.
[28] R. Heikes and D. A. Randall. Numerical integration of the shallow-water equations on a twisted icosahedral grid. part ii. a detailed description of the grid and an analysis of numerical accuracy. Monthly Weather Review, 123(6):1881 – 1887, 1995.
[29] R. P. Heikes, D. A. Randall, and C. S. Konor. Optimized icosahedral grids: Performance of finite-difference operators and multigrid solver. Monthly Weather Review, 141(12):4450 – 4469, 2013.
[30] R. Hiemstra, D. Toshniwal, R. Huijsmans, and M. Gerritsma. High order geometric methods with exact conservation properties. Journal of Computational Physics, 257(Part B):1444 – 1471, 2014. Physics-compatible numerical methods.
[31] R. Hiptmair. Discrete hodge-operators: an algebraic perspective - abstract. Journal of Electromagnetic Waves and Applications, 15(3):343–344, 2001.
[32] A. N. Hirani. Discrete Exterior Calculus. PhD thesis, California Institute of Technology, May 2003.
[33] A. N. Hirani, K. Kalyanaraman, and E. B. VanderZee. Delaunay hodge star. Computer-Aided Design, 45(2):540–544, 2013. Solid and Physical Modeling 2012.
[34] A. N. Hirani, K. Kalyanaraman, and E. B. VanderZee. Corrigendum to “delaunay hodge star” [comput. aided des. 45 (2013) 540–544]. Computer-Aided Design, 96:59–60, 2018.
[35] A. Hollingsworth, P. Kållberg, V. Renner, and D. M. Burridge. An internal symmetric computational instability. Quarterly Journal of the Royal Meteorological Society, 109(460):417–428, 1983.
[36] D. D. Holm, J. E. Marsden, and T. S. Ratiu. The Euler–Poincaré equations and semidirect products with applications to continuum theories. Advances in Mathematics, 137(1):1 – 81, 1998.
[37] S. ichi Iga and H. Tomita. Improved smoothness and homogeneity of icosahedral grids using the spring dynamics method. Journal of Computational Physics, 258:208–226, 2014.
[38] D. W. Jacobsen, M. Gunzburger, T. Ringler, J. Burkardt, and J. Peterson. Parallel algorithms for planar and spherical delaunay construction with an application to centroidal voronoi tessellations. Geoscientific Model Development, 6(4):1353–1365, 2013.
[39] N. K.-R. Kevlahan and T. Dubos. Wavetrisk-1.0: an adaptive wavelet hydrostatic dynamical core. Geoscientific Model Development, 12(11):4901–4921, 2019.
[40] P. R. Kotiuga. Theoretical limitations of discrete exterior calculus in the context of computational electromagnetics. IEEE Transactions on Magnetics, 44(6):1162–1165, 2008.
[41] J. Kreeft, A. Palha, and M. Gerritsma. Mimetic framework on curvilinear quadrilaterals of arbitrary order, 2011.
[42] L. Lazić, Z. Janjić, and F. Mesinger. “non-cancellation” instability in horizontal advection schemes for momentum equations. Meteorology and Atmospheric Physics, 35(1):49–52, Mar 1986.
[43] L. Lu, F. Sun, H. Pan, and W. Wang. Global optimization of centroidal voronoi tessellation with monte carlo approach. IEEE Transactions on Visualization and Computer Graphics, 18(11):1880–1890, 2012.
[44] M. S. Mohamed, A. N. Hirani, and R. Samtaney. Comparison of discrete hodge star operators for surfaces. Computer-Aided Design, 78:118–125, 2016. SPM 2016.
[45] M. S. Mohamed, A. N. Hirani, and R. Samtaney. Numerical convergence of discrete exterior calculus on arbitrary surface meshes. International Journal for Computational Methods in Engineering Science and Mechanics, 19(3):194–206, 2018.
[46] P. Mullen, P. Memari, F. de Goes, and M. Desbrun. Hot: Hodge-optimized triangulations. In ACM SIGGRAPH 2011 papers, pages 1–12. 2011.
[47] S. Ničković, M. B. Gavrilov, and I. A. Tošić. Geostrophic adjustment on hexagonal grids. Monthly Weather Review, 130(3):668–683, 2002.
[48] P. Peixoto, J. Thuburn, and M. Bell. Numerical instabilities of spherical shallow water models considering small equivalent depths. Quarterly Journal of the Royal Meteorological Society, pages n/a–n/a. QJ-17-0211.R1.
[49] P. S. Peixoto. Accuracy analysis of mimetic finite volume operators on geodesic grids and a consistent alternative. Journal of Computational Physics, 310(Supplement C):127 – 160, 2016.
[50] D. A. Randall. Geostrophic adjustment and the finite-difference shallow-water equations. Monthly Weather Review, 122(6):1371–1377, 1994.
[51] T. Ringler, M. Petersen, R. L. Higdon, D. Jacobsen, P. W. Jones, and M. Maltrud. A multi-resolution approach to global ocean modeling. Ocean Modelling, 69:211–232, 2013.
[52] T. Ringler, J. Thuburn, J. Klemp, and W. Skamarock. A unified approach to energy conservation and potential vorticity dynamics for arbitrarily-structured c-grids. Journal of Computational Physics, 229(9):3065 – 3090, 2010.
[53] D. Y. L. Roux. Spurious inertial oscillations in shallow-water models. Journal of Computational Physics, 231(24):7959 – 7987, 2012.
[54] R. Sadourny. The dynamics of finite-difference models of the shallow-water equations. Journal of the Atmospheric Sciences, 32(4):680–689, 1975.
[55] R. Salmon. Poisson-bracket approach to the construction of energy- and potential-enstrophy-conserving algorithms for the shallow-water equations. Journal of the Atmospheric Sciences, 61(16):2016–2036, 2004.
[56] T. G. Shepherd. Symmetries, conservation laws, and Hamiltonian structure in geophysical fluid dynamics. volume 32 of Advances in Geophysics, pages 287–338. Elsevier, 1990.
[57] T. G. Shepherd. A unified theory of available potential energy. Atmosphere-Ocean, 31(1):1–26, 1993.
[58] W. C. Skamarock, J. B. Klemp, M. G. Duda, L. D. Fowler, S.-H. Park, and T. D. Ringler. A multiscale nonhydrostatic atmospheric model using centroidal voronoi tesselations and c-grid staggering. Monthly Weather Review, 140(9):3090 – 3105, 2012.
[59] N. Soontiens and S. E. Allen. Modelling sensitivities to mixing and advection in a sill-basin estuarine system. Ocean Modelling, 112:17–32, 2017.
[60] A. Staniforth and J. Thuburn. Horizontal grids for global weather and climate prediction models: a review. Quarterly Journal of the Royal Meteorological Society, 138(662):1–26, 2012.
[61] A. L. Stewart and P. J. Dellar. An energy and potential enstrophy conserving numerical scheme for the multi-layer shallow water equations with complete coriolis force. Journal of Computational Physics, 313:99–120, 2016.
[62] K. Takano and M. G. Wurtele. A fourth order energy and potential enstrophy conserving difference scheme. California Univ., Los Angeles Report, June 1982.
[63] T. Tarhasaari, L. Kettunen, and A. Bossavit. Some realizations of a discrete hodge operator: a reinterpretation of finite element techniques [for em field analysis]. IEEE Transactions on Magnetics, 35(3):1494–1497, 1999.
[64] J. Thuburn. Numerical wave propagation on the hexagonal c-grid. Journal of Computational Physics, 227(11):5836 – 5858, 2008.
[65] J. Thuburn and C. J. Cotter. A framework for mimetic discretization of the rotating shallow-water equations on arbitrary polygonal grids. SIAM Journal on Scientific Computing, 34(3):B203–B225, 2012.
[66] J. Thuburn and C. J. Cotter. A primal–dual mimetic finite element scheme for the rotating shallow water equations on polygonal spherical meshes. Journal of Computational Physics, 290(Supplement C):274 – 297, 2015.
[67] J. Thuburn, C. J. Cotter, and T. Dubos. A mimetic, semi-implicit, forward-in-time, finite volume shallow water model: comparison of hexagonal-icosahedral and cubed-sphere grids. Geoscientific Model Development, 7(3):909–929, 2014.
[68] J. Thuburn, T. Ringler, W. Skamarock, and J. Klemp. Numerical representation of geostrophic modes on arbitrarily structured c-grids. Journal of Computational Physics, 228(22):8321 – 8335, 2009.
[69] H. Tomita, M. Satoh, and K. Goto. An optimization of the icosahedral grid modified by spring dynamics. Journal of Computational Physics, 183(1):307–331, 2002.
[70] E. Tonti. The mathematical structure of classical and relativistic physics. Springer, 2013.
[71] E. Tonti. Why starting from differential equations for computational physics? Journal of Computational Physics, 257(Part B):1260 – 1290, 2014. Physics-compatible numerical methods.
[72] M. Tort, T. Dubos, and T. Melvin. Energy-conserving finite-difference schemes for quasi-hydrostatic equations. Quarterly Journal of the Royal Meteorological Society, 141(693):3056–3075, 2015.
[73] M. D. Toy and R. D. Nair. A potential enstrophy and energy conserving scheme for the shallow-water equations extended to generalized curvilinear coordinates. Monthly Weather Review, 145(3):751–772, 2017.
[74] A. K. Turner, W. H. Lipscomb, E. C. Hunke, D. W. Jacobsen, N. Jeffery, D. Engwirda, T. D. Ringler, and J. D. Wolfe. Mpas-seaice (v1.0.0): Sea-ice dynamics on unstructured voronoi meshes. Geoscientific Model Development Discussions, 2021:1–46, 2021.
[75] H. Weller. Controlling the computational modes of the arbitrarily structured c grid. Monthly Weather Review, 140(10):3220–3234, 2012.
[76] H. Weller. Non-orthogonal version of the arbitrary polygonal c-grid and a new diamond grid. Geoscientific Model Development, 7(3):779–797, 2014.
[77] H. Weller, J. Thuburn, and C. J. Cotter. Computational modes and grid imprinting on five quasi-uniform spherical c grids. Monthly Weather Review, 140(8):2734–2755, 2012.

Appendix A Proof of $\mathbf{D}_{k}$ / $\mathbf{\bar{D}}_{k}$ / $\left<\left<{},{}\right>\right>$ Properties

Start by noting that the co-boundary operator versions of the discrete exterior derivatives $\mathbf{D}_{k}$ and $\mathbf{\bar{D}}_{k}$ satisfy by construction [32]:

\mathbf{\bar{D}}_{2}=-\mathbf{D}_{1}^{T},\quad\quad\quad\mathbf{D}_{2}=\mathbf{\bar{D}}_{1}^{T},

(150)

and

\mathbf{D}_{1}{c}^{0}=0,\quad\quad\quad\mathbf{\bar{D}}_{1}\tilde{c}^{0}=0,

(151)

for arbitrary constants ${c}^{0}$ and $\tilde{c}^{0}$ , which are just (92) and (89). Using (150) plus the definition of the topological pairing (72) - (73), it is easy to show that integration by parts (87) and (88) hold:

	$\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,$		(152)
	$\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.$		(153)

Finally, applying integration by parts for the case $k=2$ (recalling we are in $n=2$ ) plus (151) with $\tilde{y}^{n-k}=\tilde{I}^{0}$ and ${y}^{n-k}={I}^{0}$ gives the discrete Stokes theorems (93):

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

(154)

Appendix B Relationship between PV wedge product Leibniz rules

Consider a discrete wedge product ${\tilde{x}^{0}}\boldsymbol{\wedge}_{h}{\tilde{y}^{1}}$ that satisfies the full Leibniz rule

{\tilde{x}^{0}}\boldsymbol{\wedge}_{h}{\mathbf{\bar{D}}_{1}\tilde{y}^{0}}+{\tilde{y}^{0}}\boldsymbol{\wedge}_{h}{\mathbf{\bar{D}}_{1}\tilde{x}^{0}}=\mathbf{D}_{1}({\tilde{x}^{0}}\boldsymbol{\wedge}_{h}^{\star}{\tilde{y}^{0}})

(155)

along with anti-symmetry (95), where ${\tilde{x}^{0}}\boldsymbol{\wedge}_{h}^{\star}{\tilde{y}^{0}}$ is the adjoint of ${\tilde{x}^{0}}\boldsymbol{\wedge}_{h}{\tilde{y}^{2}}$ . Given coefficients for ${\tilde{x}^{0}}\boldsymbol{\wedge}_{h}^{\star}{\tilde{y}^{0}}$ (i.e. really coefficients for ${\tilde{x}^{0}}\boldsymbol{\wedge}_{h}{\tilde{y}^{2}}$ ), it seems likely that ${\tilde{x}^{0}}\boldsymbol{\wedge}_{h}{\tilde{y}^{1}}$ can be defined in such a way that (155) and (95) are both satisfied, following the same sort of general procedure as in [23]. In fact, it is hypothesized that the ${\tilde{x}^{0}}\boldsymbol{\wedge}_{h}{\tilde{y}^{1}}$ wedge product defined in [23] satisfies the more general conditions (155) rather than just (97). Since, as shown below, a wedge product that satisfies (155) also satisfies both (96) and (97), this would explain why PV compatibility (96) did not have to be enforced in [23], but was found to hold anyways.

B.1 Partial Leibniz rule (96)

Start by taking the topological pairing of (155) with an arbitrary twisted $1$ -form $\tilde{z}^{1}$ :

\left<\left<{\tilde{z}^{1}},{{\tilde{x}^{0}}\boldsymbol{\wedge}_{h}{\mathbf{\bar{D}}_{1}\tilde{y}^{0}}}\right>\right>+\left<\left<{\tilde{z}^{1}},{{\tilde{y}^{0}}\boldsymbol{\wedge}_{h}{\mathbf{\bar{D}}_{1}\tilde{x}^{0}}}\right>\right>=\left<\left<{\tilde{z}^{1}},{\mathbf{D}_{1}({\tilde{x}^{0}}\boldsymbol{\wedge}_{h}^{\star}{\tilde{y}^{0}})}\right>\right>.

(156)

Now let $\tilde{x}^{0}=\tilde{I}^{0}$ such that

\left<\left<{\tilde{z}^{1}},{{\tilde{I}^{0}}\boldsymbol{\wedge}_{h}{\mathbf{\bar{D}}_{1}\tilde{y}^{0}}}\right>\right>+\left<\left<{\tilde{z}^{1}},{{\tilde{y}^{0}}\boldsymbol{\wedge}_{h}{\mathbf{\bar{D}}_{1}\tilde{I}^{0}}}\right>\right>=\left<\left<{\tilde{z}^{1}},{\mathbf{D}_{1}({\tilde{I}^{0}}\boldsymbol{\wedge}_{h}^{\star}{\tilde{y}^{0}})}\right>\right>,

(157)

and use (151) to eliminate the middle term yielding

\left<\left<{\tilde{z}^{1}},{{\tilde{I}^{0}}\boldsymbol{\wedge}_{h}{\mathbf{\bar{D}}_{1}\tilde{y}^{0}}}\right>\right>=\left<\left<{\tilde{z}^{1}},{\mathbf{D}_{1}({\tilde{I}^{0}}\boldsymbol{\wedge}_{h}^{\star}{\tilde{y}^{0}})}\right>\right>.

(158)

Apply IBP (87) and (88) on the rhs which gives

\left<\left<{\tilde{z}^{1}},{{\tilde{I}^{0}}\boldsymbol{\wedge}_{h}{\mathbf{\bar{D}}_{1}\tilde{y}^{0}}}\right>\right>=\left<\left<{\mathbf{\bar{D}}_{2}\tilde{z}^{1}},{{\tilde{I}^{0}}\boldsymbol{\wedge}_{h}^{\star}{\tilde{y}^{0}}}\right>\right>.

(159)

Now use the definition of operator adjoints (82) or (84) on the right term and antisymmetry (95) on the lhs to get

-\left<\left<{\mathbf{\bar{D}}_{1}\tilde{y}^{0}},{{\tilde{I}^{0}}\boldsymbol{\wedge}_{h}^{\star}{\tilde{z}^{1}}}\right>\right>=\left<\left<{\tilde{y}^{0}},{{\tilde{I}^{0}}\boldsymbol{\wedge}_{h}{\mathbf{\bar{D}}_{2}\tilde{z}^{1}}}\right>\right>.

(160)

Finally, apply IBP (87) and (88) on the lhs to result in

\left<\left<{\tilde{y}^{0}},{\mathbf{D}_{2}({\tilde{I}^{0}}\boldsymbol{\wedge}_{h}^{\star}{\tilde{z}^{1}})}\right>\right>=\left<\left<{\tilde{y}^{0}},{{\tilde{I}^{0}}\boldsymbol{\wedge}_{h}{\mathbf{\bar{D}}_{2}\tilde{z}^{1}}}\right>\right>.

(161)

Since this must hold for arbitrary $\tilde{y}^{0}$ , we get

\mathbf{D}_{2}({\tilde{I}^{0}}\boldsymbol{\wedge}_{h}^{\star}{\tilde{z}^{1}})={\tilde{I}^{0}}\boldsymbol{\wedge}_{h}{\mathbf{\bar{D}}_{2}\tilde{z}^{1}},

(162)

which is just (96).

B.2 Partial Leibniz rule (97)

The partial Leibniz rule (97) is just a special case of (155) with $\tilde{x}^{0}=\tilde{y}^{0}=\tilde{q}^{0}$ :

{\tilde{q}^{0}}\boldsymbol{\wedge}_{h}{\mathbf{\bar{D}}_{1}\tilde{q}^{0}}+{\tilde{q}^{0}}\boldsymbol{\wedge}_{h}{\mathbf{\bar{D}}_{1}\tilde{q}^{0}}=\mathbf{D}_{1}({\tilde{q}^{0}}\boldsymbol{\wedge}_{h}^{\star}{\tilde{q}^{0}}),

(163)

which by inspection is equal to (97).

Appendix C Derivation of $\operatorname{\mathcal{H}}$ functional derivatives

Start by splitting the Hamiltonian into two parts: a kinetic energy and a potential energy as

\operatorname{\mathcal{H}}[{v}^{1},\tilde{h}^{2}]=\operatorname{\mathcal{H}}_{p}[\tilde{h}^{2}]+\operatorname{\mathcal{H}}_{kin}[{v}^{1},\tilde{h}^{2}],

(164)

where

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

(165)

recalling that ${v}^{1}={u}^{1}+{R}^{1}$ . The variations of $\operatorname{\mathcal{H}}_{p}$ are simple

\delta\operatorname{\mathcal{H}}_{p}=g\left<{\delta\tilde{h}^{2}},{\tilde{h}^{2}}\right>+g\left<{\delta\tilde{h}^{2}},{\tilde{h_{s}}^{2}}\right>,

(166)

and yield

{B}^{0}_{p}:=\frac{\tilde{\delta}\!\operatorname{\mathcal{H}}_{p}}{\delta\tilde{h}^{2}}=g({h}^{0}+{h_{s}}^{0}).

(167)

The variations of $\operatorname{\mathcal{H}}_{kin}$ are more complicated:

\delta\operatorname{\mathcal{H}}_{kin}=\left<{\delta\tilde{h}^{2}},{\frac{{{u}^{1}}\boldsymbol{\wedge}_{h}{\tilde{u}^{n-1}}}{2}}\right>+\left<{\tilde{h}^{2}},{\frac{{\delta{u}^{1}}\boldsymbol{\wedge}_{h}{\tilde{u}^{n-1}}}{2}}\right>+\left<{\tilde{h}^{2}},{\frac{{{u}^{1}}\boldsymbol{\wedge}_{h}{\delta\tilde{u}^{n-1}}}{2}}\right>.

(168)

The first term gives

{B}^{0}_{kin}:=\frac{\tilde{\delta}\!\operatorname{\mathcal{H}}_{kin}}{\delta\tilde{h}^{2}}=\frac{1}{2}\mathbf{\bar{H}}_{2}({{u}^{1}}\boldsymbol{\wedge}_{h}{\tilde{u}^{n-1}});

(169)

the remaining terms are

\left<{\tilde{h}^{2}},{\frac{{\delta{u}^{1}}\boldsymbol{\wedge}_{h}{\tilde{u}^{n-1}}}{2}}\right>+\left<{\tilde{h}^{2}},{\frac{{{u}^{1}}\boldsymbol{\wedge}_{h}{\delta\tilde{u}^{n-1}}}{2}}\right>.

(170)

We continue by writing these in terms of topological pairings as

\left<\left<{{h}^{0}},{\frac{{\delta{u}^{1}}\boldsymbol{\wedge}_{h}{\tilde{u}^{n-1}}}{2}}\right>\right>+\left<\left<{{h}^{0}},{\frac{{{u}^{1}}\boldsymbol{\wedge}_{h}{\delta\tilde{u}^{n-1}}}{2}}\right>\right>.

(171)

Now use the definition of operator adjoints (82) and (84) to get

\left<\left<{\delta{u}^{1}},{\frac{{{h}^{0}}\boldsymbol{\wedge}_{h}^{\star}{\tilde{u}^{n-1}}}{2}}\right>\right>+\left<\left<{\delta\tilde{u}^{n-1}},{\frac{{{u}^{1}}\boldsymbol{\wedge}_{h}^{\star}{{h}^{0}}}{2}}\right>\right>.

(172)

Since $\delta\tilde{u}^{n-1}=\mathbf{H}_{1}\delta{u}^{1}$ and $\tilde{u}^{n-1}=\mathbf{H}_{1}{u}^{1}$ , this gives

\tilde{F}^{n-1}_{kin}=\frac{\tilde{\delta}\!\operatorname{\mathcal{H}}_{kin}}{\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}}).

(173)

Finally, we can obtain $\tilde{F}^{n-1}$ and ${B}^{0}$ as

\tilde{F}^{n-1}:=\frac{\tilde{\delta}\!\operatorname{\mathcal{H}}}{\delta{v}^{1}}=\tilde{F}^{n-1}_{kin},

(174)

and

{B}^{0}:=\frac{\tilde{\delta}\!\operatorname{\mathcal{H}}}{\delta\tilde{h}^{2}}={B}^{0}_{p}+{B}^{0}_{kin},

(175)

which are the same as (99).

For alternative choice of $\operatorname{\mathcal{H}}_{kin}$ similar procedures would yield the same form for the functional derivatives, just with a change of which wedge product is primary ( ${}\boldsymbol{\wedge}_{h}{}$ ) and which ones are adjoints ( ${}\boldsymbol{\wedge}_{h}^{\star}{}$ ).

C.1 Simplified form for $\frac{\tilde{\delta}\!\operatorname{\mathcal{H}}}{\delta{v}^{1}}$

In the continuous RSWE, the mass flux $\frac{1}{2}\left[{h}^{0}\wedge\tilde{u}^{n-1}+\operatorname{\tilde{\star}}({h}^{0}\wedge{u}^{1})\right]$ can be simplified as ${h}^{0}\wedge\tilde{u}^{n-1}$ . A discrete version of this would require that

\mathbf{H}_{1}({{h}^{0}}\boldsymbol{\wedge}_{h}^{\star}{{u}^{1}})={{h}^{0}}\boldsymbol{\wedge}_{h}^{\star}{\mathbf{H}_{1}{u}^{1}},

(176)

and therefore $\tilde{F}^{n-1}={{h}^{0}}\boldsymbol{\wedge}_{h}^{\star}{\mathbf{H}_{1}{u}^{1}}={{h}^{0}}\boldsymbol{\wedge}_{h}^{\star}{\tilde{u}^{n-1}}$ . To determine the conditions on the operators required, recall the adjoint wedge product for the mass flux (148) - (149) defined from the nearest-neighbor stencil form of KE wedge product (145):

	$\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},$		(177)
	$\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}},$		(178)

where $v$ is the straight grid vertex corresponding to twisted grid cell $\tilde{c}$ . Inserting (177) - (178) plus the definition of $\mathbf{H}_{1}$ into (176) gives an equation

\sum_{e}H_{\tilde{e},e}\sum_{\tilde{c}\in CE(e)}{h}^{0}_{v}{u}^{1}_{e}=\sum_{\tilde{c}\in CE(\tilde{e})}{h}^{0}_{v}\sum_{e^{\prime}}H_{\tilde{e},e^{\prime}}{u}^{1}_{e^{\prime}}

(179)

that must hold at each twisted grid edge $\tilde{e}$ . The two sides of this equation will be equal only for a diagonal Hodge star (such as the Voronoi Hodge star), otherwise there will be unmatched ${h}^{0}_{v}$ terms on the left hand side. The conditions for a more general definition of KE wedge product could be derived in a similar way, but it seems quite unlikely that anything but a diagonal Hodge star would ensure (176) holds.

Appendix D 2D Split Exterior $\leftrightarrow$ Vector Calculus Identities

The relationships between 2D split exterior calculus and 2D vector calculus are collected here for reference. Start by identifying a scalar $f$ with a straight $0$ -form ${f}^{0}$ and a pseudoscalar $\tilde{f}$ with a twisted $0$ -form $\tilde{f}^{0}$ . As discussed in Section 2.3, given the twisted volume form $\tilde{\mu}^{n}$ , the vector $\mathbf{x}$ , and the pseudovector $\tilde{\mathbf{x}}$ , there are four vector proxies: the straight $1$ -form ${x}^{1}=\mathbf{x}^{\flat}$ and twisted $n-1$ form $\tilde{x}^{n-1}=\mathrm{i}_{\mathbf{x}}\tilde{\mu}^{n}=\operatorname{\tilde{\star}}{x}^{1}$ associated with $\mathbf{x}$ ; and the twisted $1$ -form $\tilde{x}^{1}=\mathbf{\tilde{x}}^{\flat}$ and straight $n-1$ form ${x}^{n-1}=\mathrm{i}_{\mathbf{\tilde{x}}}\tilde{\mu}^{n}=-\operatorname{\tilde{\star}}\tilde{x}^{1}$ associated with $\tilde{\mathbf{x}}$ . In 2D, $n-1=1$ and $1$ -form and $(n-1)$ -forms appear to be the same object, leading to much confusion. Therefore we retain the notation ${x}^{n-1}$ and $\tilde{x}^{n-1}$ to properly distinguish between $1$ and $(n-1)$ -forms. We can convert between twisted and straight forms using $\tilde{I}^{0}\wedge{x}^{k}=\tilde{x}^{k}$ and $\tilde{I}^{0}\wedge\tilde{x}^{k}={x}^{k}$ , and between straight 0-forms and twisted $n$ -forms using $\tilde{x}^{n}={x}^{0}\wedge\tilde{\mu}^{n}=\operatorname{\tilde{\star}}{x}^{0}$ .

In 2D (i.e. $n=2$ ) the following relationships hold (see [1] for proofs), where $f$ and $g$ are scalars, $\tilde{f}$ and $\tilde{g}$ are pseudoscalars, $\mathbf{x}$ and $\mathbf{y}$ are vectors and $\tilde{\mathbf{x}}$ and $\mathbf{\tilde{y}}$ are pseudovectors:

•

gradient, skew-gradient, divergence and curl in vector calculus $\leftrightarrow$ Hodge star and exterior derivative in split exterior calculus:

	$\displaystyle(\nabla f)^{\flat}=\operatorname{d}{f}^{0},\quad\quad(\nabla\tilde{f})^{\flat}=\operatorname{d}\tilde{f}^{0},\quad\quad(\nabla^{\perp}f)^{\flat}=\operatorname{\tilde{\star}}\operatorname{d}{f}^{0},\quad\quad(\nabla^{\perp}\tilde{f})^{\flat}=\operatorname{\tilde{\star}}\operatorname{d}\tilde{f}^{0},$		(180)
	$\displaystyle\nabla^{\perp}\cdot\mathbf{x}=\operatorname{\tilde{\star}}\operatorname{d}{x}^{1},\quad\quad\nabla^{\perp}\cdot\mathbf{\tilde{x}}=\operatorname{\tilde{\star}}\operatorname{d}\tilde{x}^{1},$		(181)
	$\displaystyle\nabla\cdot\mathbf{x}=\operatorname{\tilde{\star}}\operatorname{d}\operatorname{\tilde{\star}}{x}^{1}=\operatorname{\tilde{\star}}\operatorname{d}\tilde{x}^{n-1},\quad\quad\nabla\cdot\mathbf{\tilde{x}}=\operatorname{\tilde{\star}}\operatorname{d}\operatorname{\tilde{\star}}\tilde{x}^{1}=\operatorname{\tilde{\star}}\operatorname{d}{x}^{n-1};$		(182)

•

scalar product in vector calculus $\leftrightarrow$ Hodge star and wedge product in split exterior calculus:

$fg={f}^{0}\wedge{g}^{0},\quad\quad\tilde{f}g=\tilde{f}^{0}\wedge{g}^{0},\quad\quad\tilde{f}\tilde{g}=\tilde{f}^{0}\wedge\tilde{g}^{0};$ (183)

•

scalar-vector product in vector calculus $\leftrightarrow$ Hodge star and wedge product in split exterior calculus:

		$\displaystyle f\mathbf{x}={f}^{0}\wedge{x}^{1},\quad\quad\tilde{f}\mathbf{x}=\tilde{f}^{0}\wedge{x}^{1},\quad\quad f\mathbf{\tilde{x}}={f}^{0}\wedge\tilde{x}^{1},\quad\quad\tilde{f}\mathbf{\tilde{x}}=\tilde{f}^{0}\wedge\tilde{x}^{1},$		(184)
		$\displaystyle f\mathbf{x}^{\perp}={f}^{0}\wedge\operatorname{\tilde{\star}}{x}^{1}={f}^{0}\wedge\tilde{x}^{n-1},\quad\quad\tilde{f}\mathbf{x}^{\perp}=\tilde{f}^{0}\wedge\operatorname{\tilde{\star}}{x}^{1}=\tilde{f}^{0}\wedge\tilde{x}^{n-1},$		(185)
		$\displaystyle f\mathbf{\tilde{x}}^{\perp}={f}^{0}\wedge\operatorname{\tilde{\star}}\tilde{x}^{1}={f}^{0}\wedge{x}^{n-1},\quad\quad\tilde{f}\mathbf{\tilde{x}}^{\perp}=\tilde{f}^{0}\wedge\operatorname{\tilde{\star}}\tilde{x}^{1}=\tilde{f}^{0}\wedge{x}^{n-1};$		(186)

•

dot product in vector calculus $\leftrightarrow$ Hodge star and wedge product in split exterior calculus:

	$\displaystyle\mathbf{x}\cdot\mathbf{y}=\operatorname{\tilde{\star}}({x}^{1}\wedge\operatorname{\tilde{\star}}{y}^{1})=\operatorname{\tilde{\star}}({x}^{1}\wedge\tilde{y}^{n-1}),\quad\quad\mathbf{x}\cdot\mathbf{\tilde{y}}=\operatorname{\tilde{\star}}({x}^{1}\wedge\operatorname{\tilde{\star}}\tilde{y}^{1})=\operatorname{\tilde{\star}}({x}^{1}\wedge{y}^{n-1}),\quad\quad$		(187)
	$\displaystyle\mathbf{\tilde{x}}\cdot\mathbf{y}=\operatorname{\tilde{\star}}(\tilde{x}^{1}\wedge\operatorname{\tilde{\star}}{y}^{1})=\operatorname{\tilde{\star}}(\tilde{x}^{1}\wedge\tilde{y}^{n-1}),\quad\quad\mathbf{\tilde{x}}\cdot\mathbf{\tilde{y}}=\operatorname{\tilde{\star}}(\tilde{x}^{1}\wedge\operatorname{\tilde{\star}}\tilde{y}^{1})=\operatorname{\tilde{\star}}(\tilde{x}^{1}\wedge{y}^{n-1});$		(188)

•

inner product in vector calculus $\leftrightarrow$ inner product in split exterior calculus:

	$\displaystyle(f,g)=\left<{f}^{0},{g}^{0}\right>=\left<\tilde{f}^{n},\tilde{g}^{n}\right>,\quad\quad(\tilde{f},\tilde{g})=\left<\tilde{f}^{0},\tilde{g}^{0}\right>=\left<{f}^{n},{g}^{n}\right>,$		(189)
	$\displaystyle(\mathbf{x},\mathbf{y})=\left<{x}^{1},{y}^{1}\right>=\left<\tilde{x}^{n-1},\tilde{y}^{n-1}\right>,\quad\quad(\mathbf{\tilde{x}},\mathbf{\tilde{y}})=\left<\tilde{x}^{1},\tilde{y}^{1}\right>=\left<{x}^{n-1},{y}^{n-1}\right>.$		(190)

Using these identities is one way to derive the split exterior calculus form of equations from vector calculus.

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

Abstract

1 Introduction

2 Split Exterior Calculus for n=2n=2

2.1 Manifolds, Orientations and (Twisted) Forms

2.2 (Twisted) Hodge Star

2.3 Vector Proxies

Remark 2.1

2.4 Topological Pairing and Topological Functional Derivatives

Definition 2.2

3 Continuous RSWE

Remark 3.1

3.1 Vector Calculus

3.1.1 Hamiltonian Form

Casimirs

Vorticity Dynamics

Linearized Equations

3.2 Split Exterior Calculus

3.2.1 Hamiltonian Form

Casimirs

Vorticity Dynamics

Linearized Equations

4 A Discrete Exterior Calculus in 2D

4.1 Grids

Remark 4.1

Remark 4.2

4.1.1 Grid Geometry

4.2 Discrete Differential Forms

4.3 Operators

4.3.1 Exterior derivative

4.3.2 Wedge product

4.3.3 Hodge star

4.4 Topological Pairing, Functional Derivatives and Inner Products

4.5 Operator Adjoints

4.6 Operator Properties

4.6.1 Exterior Derivative and Topological Pairing

4.6.2 Wedge Product

4.6.3 Hodge Star

5 Discretization of RSWE with DEC

5.1 Discrete Hamiltonian and Functional Derivatives

5.2 Discrete Poisson Brackets

5.3 Semi-discrete Equations of Motion

5.4 Linearized Equations

6 Scheme Properties

6.1 Topological Properties

6.1.1 Total Energy Conservation

6.1.2 Total Mass Conservation

6.1.3 Total Circulation Conservation

6.1.4 Total Potential Enstrophy Conservation

6.1.5 No Spurious Vorticity Production

6.1.6 PV Compatibility/Steady Geostrophic Modes

6.2 Metric Properties

6.3 Summary of Scheme Properties

7 Specific Operator Choices

7.1 TRiSK Notation

7.2 Choices for Hodge Stars: 𝐈\operatorname{\mathbf{I}}, 𝐉\operatorname{\mathbf{J}} and 𝐇\operatorname{\mathbf{H}}

7.2.1 Alternative Hodge Stars

7.3 Choices for PV Wedge Products 𝐑\operatorname{\mathbf{R}}, 𝐖\operatorname{\mathbf{W}} and 𝐐\operatorname{\mathbf{Q}}

7.3.1 Choice of 𝐑\operatorname{\mathbf{R}}/𝐖\operatorname{\mathbf{W}}

7.3.2 Choice of 𝐐\operatorname{\mathbf{Q}}

7.4 Choices for KE Wedge Product

Alternative Choices

7.5 Choice of grid

7.6 Grid Optimization

7.7 Effects of choices on properties

8 Schemes from the literature

8.1 Degree of freedom scaling

8.2 Specific schemes

AL81

SD75

TW84

S04

TRSK2010

Thuburn2014

Weller2014

Eldred2017

Toy2017

8.3 Unexplored approaches

9 Conclusions

10 Acknowledgements

2 Split Exterior Calculus for $n=2$

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

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

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

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

Appendix A Proof of $\mathbf{D}_{k}$ / $\mathbf{\bar{D}}_{k}$ / $\left<\left<{},{}\right>\right>$ Properties

Appendix C Derivation of $\operatorname{\mathcal{H}}$ functional derivatives

C.1 Simplified form for $\frac{\tilde{\delta}\!\operatorname{\mathcal{H}}}{\delta{v}^{1}}$

Appendix D 2D Split Exterior $\leftrightarrow$ Vector Calculus Identities