Conservation and stability in a discontinuous Galerkin method for the vector invariant spherical shallow water equations

We use the vector invariant form of the RSWE on a 2D manifold $\Omega$ embedded in $R^{3}$ with periodic boundary conditions. We choose periodic boundary conditions as they simplify the analysis and the domain of most practical interest is the sphere. The RSWE are

$$\mathbf{u}_{t}+\omega\mathbf{k}\crossproduct\mathbf{u}+\nabla G=0\text{,}$$

(1)

$$D_{t}+\nabla\cdot\mathbf{F}=0\text{,}$$

(2)

with

$$\omega=\mathbf{k}\cdot\nabla\crossproduct\mathbf{u}+f\text{,}\quad\mathbf{F}=D\mathbf{u}\text{,}\quad G=\tfrac{1}{2}\mathbf{u}\cdot\mathbf{u}+gD\text{,}$$

(3)

where $t\geq 0$ denotes the time variable, the unknowns are the flow velocity vector $\mathbf{u}=[u,v,w]^{T}$ and the fluid depth $D$. Note that here and in the rest of the paper all differential operators are taken to be constrained to the manifold. Furthermore, $\omega$ is the absolute vorticity, $f$ is the Coriolis frequency, $\mathbf{F}$ is the mass flux, $G$ is the potential, and $\mathbf{k}$ is normal of the manifold. At $t=0$, we augment the RSWE (1)–(2) with the initial conditions

$\displaystyle\mathbf{u}=\mathbf{u}_{0}(\mathbf{x}),\quad{D}={D}_{0}(\mathbf{x}),\quad\mathbf{x}\in\Omega,$

(4)

such that the flow is constrained to the manifold by ensuring that $\mathbf{u}_{0}\cdot\mathbf{k}=0$. The constraint of the differential operators and the initial constraint $\mathbf{u}_{0}\cdot\mathbf{k}=0$ is sufficient to ensure that the $\mathbf{u}(\mathbf{x},t)\cdot\mathbf{k}=0$ holds $\forall\mathbf{x}\in\Omega,\;t\geq 0$.

We will show that the total mass, absolute vorticity and energy are conserved. Integrating the continuity equation (2) over a periodic domain yields conservation of total mass, that is

$\displaystyle\frac{d}{dt}\int_{\Omega}Dd\Omega=0.$

(5)

We can derive a continuity equation for absolute vorticity $\omega$ by applying $\mathbf{k}\cdot\nabla\crossproduct$ to the momentum equation (1) giving

$$\omega_{t}+\nabla\cdot(\omega\mathbf{u})=0\text{.}$$

(6)

Similarly, integrating equation (6) over a periodic domain yields conservation of total absolute vorticity,

$\displaystyle\frac{d}{dt}\int_{\Omega}\omega d\Omega=0.$

(7)

Additionally, all quantities of the form $Db\left(\frac{\omega}{D}\right)$, for arbitrary function $b$ of $\frac{\omega}{D}$, are also locally conserved by the RSWE (Vallis, 2017). The potential vorticity $\frac{\omega}{D}$ is transported by the flow, therefore $b\left(\frac{\omega}{D}\right)$ is also transported by the flow, and hence $Db\left(\frac{\omega}{D}\right)$ is locally conserved. The potential enstrophy is one such quantity, and while we do not target its conservation here several other numerical schemes have (McRae and Cotter, 2014; Lee et al., 2018).

Let the elemental energy $E$ be denoted by

$$E=\tfrac{1}{2}D\mathbf{u}\cdot\mathbf{u}+\tfrac{1}{2}gD^{2}.$$

(8)

The elemental energy $E$ also satisfies a continuity equation. To show this we take the time derivative of the elemental energy $E$ giving

$$E_{t}=\mathbf{F}\cdot\mathbf{u}_{t}+(\tfrac{1}{2}\mathbf{u}\cdot\mathbf{u}+gD)D_{t}=-\mathbf{F}\cdot\omega\mathbf{k}\crossproduct\mathbf{u}-\mathbf{F}\cdot\nabla G-G\nabla\cdot\mathbf{F}.$$

(9)

We have the continuity equation for the energy $E$

$$E_{t}+\nabla\cdot\big{(}G\mathbf{F}\big{)}=0,$$

(10)

where we have used $\mathbf{F}\cdot\omega\mathbf{k}\crossproduct\mathbf{u}=\omega D\mathbf{u}\cdot\mathbf{k}\crossproduct\mathbf{u}=0$ and the product rule for derivatives. As before, integrating equation (10) over a periodic domain yields conservation of total energy,

$\displaystyle\frac{d}{dt}\int_{\Omega}Ed\Omega=0.$

(11)

On domains with constant Coriolis parameter $f$, in addition to a set of inertio-gravity wave solutions, the linear rotating shallow water equations also exhibit steady geostrophic modes (Cotter and Shipton, 2012; Vallis, 2017). The RSWE (1)–(2) linearised around a state of constant mean height $H$ and zero mean flow are

$$\mathbf{u}_{t}+f{\mathbf{k}\crossproduct}\mathbf{u}+g\nabla D=0\text{,}$$

(12)

$$D_{t}+H\nabla\cdot\mathbf{u}=0\text{.}$$

(13)

The steady solutions to these equations can be found by choosing an arbitrary stream function $\psi$ and setting $D=-\tfrac{f}{g}\psi,\quad\mathbf{u}=\nabla{\crossproduct\psi\mathbf{k}}$. Then we have

$$D_{t}=-{H}\nabla\cdot\mathbf{u}=-{H\nabla\cdot\nabla\crossproduct\psi\mathbf{k}}=0,$$

(14)

$$\mathbf{u}_{t}=-(f{\mathbf{k}\crossproduct}\mathbf{u}+g\nabla D)=-(f\nabla\psi-f\nabla\psi{)=0}.$$

(15)

An entropy function is any convex combination of prognostic variables that also satisfies a conservation law (LeVeque, 1992). The importance of entropy functions is that they are conserved in smooth regions and are dissipated across discontinuities in physically relevant weak solutions. DG methods are continuous within each element but can be discontinuous across element boundaries, this motivates that the volume terms of DG methods should conserve entropy while numerical fluxes should dissipate entropy. We call this local entropy stability and in practice many numerical methods have found that this is sufficient for numerical stability (Giraldo, 2020).

It is well known that energy is an entropy function of the RSWE in conservative form (Gassner et al., 2016), here we show that energy is also an entropy for the RSWE in vector invariant form for subcritical flow $|\mathbf{u}|<\sqrt{gD}$. Atmospheric flows are subcritical and so this restriction does not cause any difficulty in practice, for the target applications. The analysis in section 2.3 shows that energy is conserved, all that remains is to show that energy is also a convex combination of $u$, $v$, $w$, and $D$. In the following analysis we non-dimensionalise all quantities by their mean values to avoid issues with conflicting units. The Hessian of the energy with respect to the prognostic variables is

$$H_{E}=\begin{pmatrix}D&0&0&u\\ 0&D&0&v\\ 0&0&D&w\\ u&v&w&g\\ \end{pmatrix}\text{,}$$

(16)

and its characteristic equation is

$$(D-\lambda)^{2}\big{(}\lambda^{2}-(D+g)\lambda+Dg-|\mathbf{u}|^{2}\big{)}\text{.}$$

(17)

As $D,g>0$ the eigenvalues corresponding to the left bracket are always positive, and as long as $Dg-|\mathbf{u}|^{2}>0$, which implies $|\mathbf{u}|<\sqrt{gD}$, the eigenvalues corresponding to the quadratic in the second bracket are also positive. Therefore for subcritical flow the energy is an entropy function of the vector invariant RSWE.

An effective numerical method should as far as possible conserve mass, absolute vorticity, energy, entropy and linear geostrophic balance.

For each element $\Omega^{q}$ we define an invertible map $r(\mathbf{x};q)$ to the reference square $[-1,1]^{2}$, and approximate solutions by polynomials of order $m$ in this reference square. We use a computationally efficient tensor product Lagrange polynomial basis with interpolation points collocated with Gauss-Lobatto-Legendre (GLL) quadrature points. Using this basis the polynomials of order $m$ on the reference square can be defined as

$$P_{m}:=\text{span}_{i,j=1}^{m+1}l_{i}(\xi)l_{j}(\eta)\text{,}$$

(18)

where $(\xi,\eta)\in[-1,1]^{2}$ and $l_{i}$ is the 1D Lagrange polynomial which interpolates the $i^{th}$ GLL node. We define a discontinuous scalar space $S$ which contains the depth $D^{h}$ as

$$S=\{\phi\in L^{2}(\Omega):r(\phi;q)\in P_{m},\forall q\}\text{.}$$

(19)

Note that within each element functions $\phi\in S$ can be expressed as

$$\phi(\mathbf{x})=\sum_{i,j=1}^{m+1}{\phi_{ijq}l_{i}(\xi)l_{j}(\eta)},\forall\mathbf{x}\in\Omega^{q}\text{,}$$

(20)

where $\xi,\eta=r(\mathbf{x};q)$, $\phi_{ijq}=\phi(r^{-1}(\xi_{i},\eta_{j};q))$.

We define a discontinuous vector space which contains the velocity $\mathbf{u}$. Let $\mathbf{v}_{1}(x,y,z)$ and $\mathbf{v}_{2}(x,y,z)$ be any independent set of vectors which span the tangent space of the manifold at each point $(x,y,z)$, then the discontinuous vector space $V$ containing the velocity $\mathbf{u}$ can be defined as

$$V=\{\mathbf{w}:\mathbf{w}\cdot\mathbf{v}_{i}\in S,i=1,2\}\text{.}$$

(21)

We note that $V$ is independent of the particular choice of $\mathbf{v}_{i}$ however two convenient choices are the contravariant and covariant basis vectors.

Here we provide a brief summary of the differential geometry tools we use in our discrete method. This summary closely follows (Taylor and Fournier, 2010) and these results can also be found in standard textbooks (Giraldo, 2020; Heinbockel, 2001). The covariant base vectors $\mathbf{g}_{1}$ and $\mathbf{g}_{2}$ are defined as

$$\mathbf{g}_{1}=\dfrac{\partial\mathbf{x}}{\partial\xi},\quad\mathbf{g}_{2}=\dfrac{\partial\mathbf{x}}{\partial\eta}\text{,}$$

(22)

with these the determinant of the Jacobian of $r(\mathbf{x};n)$ can be expressed as

$$J=|\mathbf{g}_{1}\times\mathbf{g}_{2}|\text{.}$$

(23)

Further, we assume that the map $r(\mathbf{x};q)$ is invertible and the Jacobian is strictly positive $J>0$.

The contravariant base vectors are defined as

$$\mathbf{g}^{1}=\nabla\xi,\quad{\mathbf{g}^{2}}=\nabla\eta\text{}$$

(24)

to fulfill the property

$$\mathbf{g}^{\alpha}\cdot\mathbf{g}_{\beta}=\delta{{}^{\alpha}_{\beta}}\text{,}$$

(25)

enabling vectors to be expressed as

$$\mathbf{w}=\sum_{\alpha=1,2}{w^{\alpha}\mathbf{g}_{\alpha}}=\sum_{\alpha=1,2}{w_{\alpha}\mathbf{g}^{\alpha}}\text{,}$$

(26)

where $w^{\alpha}=\mathbf{w}\cdot\mathbf{g}^{\alpha}$, $w_{\beta}=\mathbf{w}\cdot\mathbf{g}_{\beta}$.

For flows such that $\mathbf{w}\cdot\mathbf{k}=0$ this enables the divergence, gradient, and curl to be expressed as

$$\nabla\cdot\mathbf{w}=\dfrac{1}{J}\bigg{(}\dfrac{\partial Jw^{1}}{\partial\xi}+\dfrac{\partial Jw^{2}}{\partial\eta}\bigg{)}\text{,}$$

(27)

$$\nabla\phi=\dfrac{\partial\phi}{\partial\xi}\mathbf{g}^{1}+\dfrac{\partial\phi}{\partial\eta}\mathbf{g}^{2}\text{,}$$

(28)

$$\nabla\crossproduct\mathbf{w}=\dfrac{1}{J}\bigg{(}\dfrac{\partial w_{2}}{\partial\eta}-\dfrac{\partial w_{1}}{\partial\xi}\bigg{)}\mathbf{k}\text{,}$$

(29)

$$\nabla\crossproduct\phi\mathbf{k}=\dfrac{1}{J}\dfrac{\partial\phi}{\partial\eta}\mathbf{g}_{1}-\dfrac{1}{J}\dfrac{\partial\phi}{\partial\xi}\mathbf{g}_{2}\text{,}$$

(30)

where we have only used $\nabla\crossproduct$ on quantities either parallel or orthogonal to $\mathbf{k}$ to simplify the exposition.

We approximate integrals over the elements by using GLL quadrature. This defines a discrete element inner product

$$\langle f,g\rangle_{\Omega^{q}}:=\sum_{i,j=1}^{m+1}{w_{i}w_{j}J_{ij}f_{ij}g_{ij}}\text{,}$$

(31)

which approximates the integral

$$\int_{\Omega^{q}}{fgd\Omega^{q}}\approx\langle f,g\rangle_{\Omega^{q}}\text{.}$$

(32)

The discrete global inner product is defined simply as $\langle f,g\rangle_{\Omega}=\sum_{q}{\langle f,g\rangle_{\Omega^{q}}}$. Similarly for vectors

$$\langle\mathbf{f},\mathbf{g}\rangle_{\Omega^{q}}:=\sum_{i,j=1}^{m+1}{w_{i}w_{j}J_{ij}\mathbf{f}_{ij}\cdot\mathbf{g}_{ij}}\text{.}$$

(33)

We also approximate element boundary integrals using GLL quadrature

$$\int_{\partial\Omega^{q}}{\mathbf{f}\cdot\mathbf{n}dl}\approx\langle\mathbf{f},\mathbf{n}\rangle_{\partial\Omega^{q}}\text{,}$$

(34)

where

$$\begin{split}\langle\mathbf{f},\mathbf{n}\rangle_{\Omega^{q}}:=\sum_{j}{w_{j}\big{(}|g_{1}|_{1j}\mathbf{f}_{1j}\cdot\mathbf{n}_{1j}+|g_{1}|_{(m+1)j}\mathbf{f}_{(m+1)j}\cdot\mathbf{n}_{(m+1)j}}\\ {+|g_{2}|_{j1}\mathbf{f}_{j1}\cdot\mathbf{n}_{j1}+|g_{2}|_{j(m+1)}\mathbf{f}_{j(m+1)}\cdot\mathbf{n}_{j(m+1)}\big{)}}\text{,}\end{split}$$

(35)

and $\mathbf{n}$ is the outward facing unit normal of the element boundary $\partial\Omega^{q}$ given by

$$\mathbf{n}=\begin{cases}\xi\dfrac{\mathbf{g}^{1}}{|\mathbf{g}^{1}|},\;\text{for}\;\xi=-1,1\\ \eta\dfrac{\mathbf{g}^{2}}{|\mathbf{g}^{2}|},\;\text{for}\;\eta=-1,1\text{.}\end{cases}$$

(36)

The unit tangent vector $\mathbf{t}$ is similarly defined as

$$\mathbf{t}=\begin{cases}\xi\dfrac{\mathbf{g}_{2}}{|\mathbf{g}_{2}|},\;\text{for}\;\xi=-1,1\\ -\eta\dfrac{\mathbf{g}_{1}}{|\mathbf{g}_{1}|},\;\text{for}\;\eta=-1,1\end{cases}\text{.}$$

(37)

With the discrete inner products defined we now construct discrete $L^{2}$ projection operators $\Pi_{S}(\cdot)$ and $\Pi_{V}(\cdot)$ which preserve discrete $L^{2}$ inner products in $S$ and $V$ respectively. For example $\Pi_{S}(\cdot)$ is defined as the unique projection that satisfies

$$\langle\gamma,b\rangle_{\Omega}=\langle\gamma,\Pi_{S}(b)\rangle_{\Omega}\text{,}\;\forall\gamma\in S\text{,}$$

(38)

for all functions $b$. $\Pi(\cdot)_{V}$ is defined similarly. In practice the discontinuous projection operators $\Pi_{S}$ and $\Pi_{V}$ simply interpolate functions within each element

$$\Pi_{S}(f)(\mathbf{x})=\sum_{ij}{f_{ijq}l_{i}(\xi)l_{j}(\eta)},\forall\mathbf{x}\in\Omega^{q}\text{,}$$

(39)

$$\Pi_{V}(\mathbf{f})(\mathbf{x})=\sum_{ij}{\mathbf{f}_{ijq}l_{i}(\xi)l_{j}(\eta)},\forall\mathbf{x}\in\Omega^{q}\text{,}$$

(40)

where $(\xi,\eta)=r(\mathbf{x};q)$.

Following (Taylor and Fournier, 2010) we construct discrete operators by discretising equations (27)–(30)

$$\nabla_{d}\cdot\mathbf{w}=\Pi_{S}\Bigg{(}\dfrac{1}{J}\bigg{(}\dfrac{\partial}{\partial\xi}\Pi_{S}(Jw^{1})+\dfrac{\partial}{\partial\eta}\Pi_{S}(Jw^{2})\bigg{)}\Bigg{)}\text{.}$$

(41)

$$\nabla_{d}\phi=\dfrac{\partial\phi}{\partial\xi}\mathbf{g}^{1}+\dfrac{\partial\phi}{\partial\eta}\mathbf{g}^{2}\text{,}$$

(42)

$$\nabla_{d}\crossproduct\mathbf{w}=\Pi_{S}\Bigg{(}\dfrac{1}{J}\bigg{(}\dfrac{\partial w_{2}}{\partial\xi}-\dfrac{\partial w_{1}}{\partial\eta}\bigg{)}\Bigg{)}\mathbf{k}\text{,}$$

(43)

$$\nabla_{d}\crossproduct\phi\mathbf{k}=\Pi_{V}\Bigg{(}\dfrac{1}{J}\bigg{(}\dfrac{\partial\phi}{\partial\eta}\mathbf{g}_{1}-\dfrac{\partial\phi}{\partial\xi}\mathbf{g}_{2}\bigg{)}\Bigg{)}\text{,}$$

(44)

Here we present the properties of the function spaces and discrete operators necessary for the numerical stability and discrete conservation proofs performed later in this paper.

First we introduce notation for jumps $[[\cdot]]$ and averages $\{\{\cdot\}\}$ across element boundaries as

$$\{\{a\}\}=\tfrac{1}{2}\big{(}a^{+}+a^{-}\big{)}\text{,}$$

(53)

$$[[a]]=a^{+}-a^{-}\text{,}$$

(54)

where $\cdot^{+}$ and $\cdot^{-}$ denote points on opposites sides of the boundary element boundary, and $a$ can be either a scalar or vector. We now define our spatial discretisation as

$$\begin{split}\langle\mathbf{w}^{h},\mathbf{u}_{t}^{h}\rangle_{\Omega^{q}}+\langle\mathbf{w}^{h},\omega^{h}\mathbf{k}\crossproduct\mathbf{u}^{h}\rangle_{\Omega^{q}}+\langle\mathbf{w}^{h},\nabla_{d}G^{h}\rangle_{\Omega^{q}}\\ +\langle\mathbf{w}^{h},(\hat{G}-G)\mathbf{n}\rangle_{\partial\Omega^{q}}=0\text{,}\;\forall\mathbf{w}^{h}\in V\text{,}\end{split}$$

(55)

$$\langle\phi^{h},D_{t}^{h}\rangle_{\Omega^{q}}+\langle\phi^{h},\nabla_{d}\cdot\mathbf{F}^{h}\rangle_{\Omega^{q}}+\langle\phi^{h},\big{(}\hat{\mathbf{F}}-\mathbf{F}^{h}\big{)}\cdot\mathbf{n}\rangle_{\partial\Omega^{q}}=0\text{,}\;\forall\phi^{h}\in S\text{,}$$

(56)

where the numerical fluxes are defined for the scalars $\alpha,\beta\geq 0$ as

$$\hat{G}=\{\{G\}\}+\alpha[[\mathbf{F}^{h}]]\cdot{\mathbf{n}^{+}}\text{.}$$

(57)

$$\hat{\mathbf{F}}\cdot\mathbf{n}^{+}=\{\{\mathbf{F}^{h}\}\}\cdot\mathbf{n}^{+}+\beta[[G^{h}]]\text{,}$$

(58)

and the discrete absolute vorticity is

$$\langle\phi^{h},\omega^{h}\rangle_{\Omega^{q}}=\langle-\nabla_{d}\crossproduct\phi^{h}\mathbf{k},\mathbf{u}^{h}\rangle_{\Omega^{q}}+\langle\phi^{h},f\rangle_{\Omega^{q}}+\langle\phi^{h},\{\{\mathbf{u}^{h}\}\}\cdot\mathbf{t}\rangle_{\partial\Omega^{q}}\;\forall\phi^{h}\in S\text{,}$$

(59)

where $\mathbf{t}$ is the unit tangent vector along element boundaries.

As will be shown below, our scheme is discretely entropy stable so long as $\alpha,\beta\geq 0$, and is locally semi-discretely energy/entropy conserving for $\alpha=\beta=0$ and locally discretely energy/entropy dissipating for $\alpha,\beta>0$. We use a centred mass flux $\beta=0$ to ensure that potential energy is not dissipated, and either use centred fluxes $\alpha=0$ or

$$\alpha=\frac{1}{2}\max\left(\frac{c^{+}}{D^{+}},\frac{c^{-}}{D^{-}}\right)\text{,}$$

(60)

where $c^{\pm}=|\mathbf{u}^{\pm}|+\sqrt{gD^{\pm}}$ is the fastest wave speed on either side of the boundary. This ensures that $\hat{G}$ has the correct units and reduces to a Rusanov flux (Rusanov, 1970) in the case of the linearised RSWE in section 2.5. We refer to the parameter choices $\alpha=\beta=0$ and $\alpha=\frac{1}{2}\max\left(\frac{c^{+}}{D^{+}},\frac{c^{-}}{D^{-}}\right),\;\beta=0$ as the energy conserving and energy dissipating methods.

Element local conservation of total mass can be shown using the element local SBP property to transform the mass equation (2) into weak form

$$\langle\phi^{h},D_{t}^{h}\rangle_{\Omega^{q}}-\langle\nabla_{d}\phi^{h},\mathbf{F}^{h}\rangle_{\Omega^{q}}+\langle\phi,\hat{\mathbf{F}}\cdot\mathbf{n}\rangle_{\partial\Omega^{q}}=0\text{,}\;\forall\phi^{h}\in S\text{,}$$

(61)

and then choosing the test function $\phi^{h}=1$ to yield the rate of change of the total mass of element $k$

$$\langle 1,D_{t}^{h}\rangle_{\Omega^{q}}+\langle 1,\hat{\mathbf{F}}\cdot\mathbf{n}\rangle_{\partial\Omega^{q}}=0\text{.}$$

(62)

By construction the numerical fluxes are continuous across element boundaries and so the mass flowing through each point along the element boundary is equal and opposite to that of the adjacent element, therefore mass is locally conserved. Global conservation can be trivially obtained from local conservation in a periodic domain by noting that mass flux cancels for adjacent element and so $\sum_{q}{\langle 1,D_{t}^{h}\rangle_{\Omega^{q}}}=-\sum_{q}{\langle 1,\hat{\mathbf{F}}\cdot\mathbf{n}\rangle_{\partial\Omega^{q}}}=0$.

We show that the potential flux $G$ does not spuriously force the absolute vorticity, and that the absolute vorticity is also locally conserved. This is a direct consequence of our choice of the numerical flux $\hat{G}$ defined in (57), numerical approximation of the absolute vorticity $\omega^{h}$ giving in (59), and the element boundary SBP property discussed in subsection 3.5.5.

The evolution of discrete absolute vorticity follows from taking the time derivative of $\omega^{h}$ in (59), giving

$$\langle\phi^{h},\omega^{h}_{t}\rangle_{\Omega^{q}}=-\langle\nabla_{d}\crossproduct\phi^{h}\mathbf{k},\mathbf{u}^{h}_{t}\rangle_{\Omega^{q}}+\langle\phi^{h},\{\{\mathbf{u}^{h}_{t}\}\}\cdot\mathbf{t}\rangle_{\partial\Omega^{q}}\text{,}\;\forall\phi^{h}\in S\text{.}$$

(63)

First we expand the volume term by using the weak form of the velocity equation (55) and the curl gradient identity (50), that is $\mathbf{n}\cdot\nabla_{d}\crossproduct\phi^{h}\mathbf{k}=\nabla_{d}\phi^{h}\cdot\mathbf{t}$, and $\nabla_{d}\cdot\nabla_{d}\crossproduct\phi^{h}\mathbf{k}=0$ yielding

$$\langle\nabla_{d}\crossproduct\phi^{h}\mathbf{k},\mathbf{u}^{h}_{t}\rangle_{\Omega^{q}}=\langle\nabla_{d}\phi^{h},\omega^{h}\mathbf{u}^{h}\rangle_{\Omega^{q}}-\langle\hat{G}\nabla\phi^{h},\mathbf{t}\rangle_{\partial\Omega^{q}}\text{,}\;\forall\phi^{h}\in S\text{.}$$

(64)

Next we evaluate the boundary term, as this is a GLL collocated method

$$\langle l_{i}(\xi)l_{j}(\eta)\mathbf{t}_{ij},\mathbf{u}^{h}_{t}\rangle_{\Omega^{q}}=w_{i}w_{j}(\mathbf{u}^{h}_{ij})_{t}\cdot\mathbf{t}_{ij}\text{,}$$

(65)

substituting the test function $l_{i}(\xi)l_{j}(\eta)\mathbf{t}_{ij}$ into (55) yields

$$\langle l_{i}(\xi)l_{j}(\eta)\mathbf{t}_{ij},\mathbf{u}^{h}_{t}\rangle_{\Omega^{q}}=-w_{i}w_{j}\omega^{h}\mathbf{k}\crossproduct\mathbf{u}^{h}\cdot\mathbf{t}-w_{i}w_{j}\nabla_{d}G\cdot\mathbf{t}\text{}$$

(66)

$$=-w_{i}w_{j}\omega^{h}\mathbf{u}^{h}\cdot\mathbf{n}-w_{i}w_{j}\nabla_{d}G\cdot\mathbf{t}\text{,}$$

(67)

and therefore

$$\{\{\mathbf{u}_{t}\}\}\cdot\mathbf{t}=-\{\{\omega^{h}\mathbf{u}^{h}\}\}\cdot\mathbf{n}-\nabla_{d}\{\{G\}\}\cdot\mathbf{t}\text{.}$$

(68)

With these the absolute vorticity evolution becomes

$$\begin{split}\langle\phi^{h},\omega^{h}_{t}\rangle_{\Omega^{q}}&=\langle\nabla_{d}\phi^{h},\omega^{h}\mathbf{u}^{h}\rangle_{\Omega^{q}}-\langle\hat{G}\nabla\phi^{h}+\phi^{h}\{\{\nabla_{d}G\}\},\mathbf{t}\rangle_{\partial\Omega^{q}}\\ &-\langle\phi^{h}\{\{\omega^{h}\mathbf{u}^{h}\}\},\mathbf{n}\rangle_{\partial\Omega^{q}}\\ &=\langle\nabla_{d}\phi^{h},\omega^{h}\mathbf{u}^{h}\rangle_{\Omega^{q}}-\langle\{\{G\}\}\nabla\phi^{h}+\phi^{h}\{\{\nabla_{d}G\}\},\mathbf{t}\rangle_{\partial\Omega^{q}}\\ &-\langle\phi^{h}\{\{\omega^{h}\mathbf{u}^{h}\}\},\mathbf{n}\rangle_{\partial\Omega^{q}}-\langle\nabla\phi^{h}\cdot\mathbf{t},\alpha[[\mathbf{F}]]\cdot\mathbf{n}^{+}\rangle_{\partial\Omega^{q}}\text{,}\end{split}$$

(69)

where we have used the definition of $\hat{G}$. Using the tangent boundary SBP property in section 3.5.5 this becomes

$$\begin{split}\langle\phi^{h},\omega^{h}_{t}\rangle_{\Omega^{q}}-\langle\nabla_{d}\phi^{h},\omega^{h}\mathbf{u}^{h}\rangle_{\Omega^{q}}+\langle\phi^{h}\{\{\omega^{h}\mathbf{u}^{h}\}\},\mathbf{n}\rangle_{\partial\Omega^{q}}\\ +\langle\nabla\phi^{h}\cdot\mathbf{t},\alpha[[\mathbf{F}]]\cdot\mathbf{n}^{+}\rangle_{\partial\Omega^{q}}=0\text{,}\forall\phi^{h}\in S\text{,}\end{split}$$

(70)

demonstrating that the potential $G$ does not spuriously force the discrete absolute vorticity.

To show local absolute vorticity conservation we choose $\phi=1$ to yield

$$\langle 1,\omega^{h}_{t}\rangle_{\Omega^{q}}={-}\langle\{\{\omega^{h}\mathbf{u}^{h}\}\},\mathbf{n}\rangle_{\partial\Omega^{q}}\text{.}$$

(71)

By construction $\{\{\omega^{h}\mathbf{u}^{h}\}\}$ is continuous across element boundaries and therefore absolute vorticity is locally conserved. Global conservation in a periodic domain follows from cancellation of $\{\{\omega^{h}\mathbf{u}^{h}\}\}\cdot\mathbf{n}$ across element boundaries.

Here we show that, for centred flux with $\alpha=\beta=0$, the semi-discrete energy is conserved and for $\alpha,\beta>0$ the semi-discrete energy is dissipated. We will obtain fully discrete energy stability by using a strong stability preserving time integrator (Shu and Osher, 1988).

Consider the element ${\Omega^{q}}$, we define the discrete energy $\mathcal{E}^{q}$ within the element as

$$\mathcal{E}^{q}=\tfrac{1}{2}\langle D^{h}\mathbf{u}^{h},\mathbf{u}^{h}\rangle_{\Omega^{q}}+\tfrac{1}{2}\langle gD^{h},D^{h}\rangle_{\Omega^{q}}\text{,}$$

(72)

with the total energy given by the sum over all elements $\mathcal{E}=\sum_{q}\mathcal{E}^{q}$. Following the continuous form of the energy conservation analysis in section 2.2 the semi-discrete time evolution of $\mathcal{E}^{q}$ is given by

$$\begin{split}\mathcal{E}^{q}_{t}&=\langle D^{h}\mathbf{u}^{h},\mathbf{u}^{h}_{t}\rangle_{\Omega^{q}}+\langle\tfrac{1}{2}\mathbf{u}^{h}\cdot\mathbf{u}^{h},D^{h}_{t}\rangle_{\Omega^{q}}+\langle gD^{h},D^{h}_{t}\rangle_{\Omega^{q}}\\ &=\langle\mathbf{F}^{h},\mathbf{u}^{h}_{t}\rangle_{\Omega^{q}}+\langle G^{h},D^{h}_{t}\rangle_{\Omega^{q}}.\end{split}$$

(73)

Using the weak forms of the mass and velocity equations (55)–(56) gives

$$\begin{split}\mathcal{E}^{q}_{t}&=-\langle\mathbf{F}^{h},\nabla_{d}G^{h}\rangle_{\Omega^{q}}+\langle\nabla_{d}G^{h},\mathbf{F}^{h}\rangle_{\Omega^{q}}-\langle\hat{G},\mathbf{F}^{h}\cdot\mathbf{n}\rangle_{\partial\Omega^{q}}\\ &-\langle G^{h},(\hat{\mathbf{F}}^{h}-\mathbf{F}^{h})\cdot\mathbf{n}\rangle_{\partial\Omega^{q}}\\ &=-\langle\hat{G},\mathbf{F}^{h}\cdot\mathbf{n}\rangle_{\partial\Omega^{q}}-\langle G^{h},(\hat{\mathbf{F}}^{h}-\mathbf{F}^{h})\cdot\mathbf{n}\rangle_{\partial\Omega^{q}}\end{split}$$

(74)

where we have used the fact that $\mathbf{F}^{h}\cdot\mathbf{k}\crossproduct\mathbf{u}^{h}=0$. Note that $\mathbf{F}^{h}\cdot\mathbf{k}\crossproduct\mathbf{u}^{h}=D^{h}\mathbf{u}^{h}\cdot\mathbf{k}\crossproduct\mathbf{u}^{h}=0$ only holds for the special case of diagonal mass matrices (Giraldo, 2020).

To prove energy stability, that is $\mathcal{E}_{t}=\sum_{q}\mathcal{E}^{q}_{t}\leq 0$, it suffices to show that the jump of the numerical energy flux across element boundaries is negative semi-definite, that is

$$-\Big{[}\Big{[}\mathbf{F}^{h}\hat{G}+G^{h}\big{(}\hat{\mathbf{F}}^{h}-\mathbf{F}^{h}\big{)}\Big{]}\Big{]}\cdot\mathbf{n}^{+}\leq 0,$$

(75)

for all points on the element interface. Using the expansion $[[ab]]=\{\{a\}\}[[b]]+[[a]]\{\{b\}\}$ the energy jump condition becomes

$$\begin{split}-\Big{[}\Big{[}\mathbf{F}^{h}\hat{G}&+G^{h}\big{(}\hat{\mathbf{F}}^{h}-\mathbf{F}^{h}\big{)}\Big{]}\Big{]}\cdot\mathbf{n}^{+}=\\ &\big{(}\{\{G^{h}\}\}-\hat{G}\big{)}[[\mathbf{F}^{h}]]\cdot\mathbf{n}^{+}+[[G^{h}]]\big{(}\{\{\mathbf{F}^{h}\}\}-\hat{\mathbf{F}}\big{)}\cdot\mathbf{n}^{+}\leq 0\text{.}\end{split}$$

(76)

For centred fluxes with $\alpha=0$, $\beta=0$, we have $\hat{G}=\{\{G^{h}\}\}$ and $\hat{\mathbf{F}}=\{\{\mathbf{F}\}\}$ the jump is $0$ and therefore the energy is conserved,

$$-\Big{[}\Big{[}\mathbf{F}^{h}\hat{G}+G^{h}\big{(}\hat{\mathbf{F}}^{h}-\mathbf{F}^{h}\big{)}\Big{]}\Big{]}\cdot\mathbf{n}^{+}=0\implies\mathcal{E}_{t}=0.$$

When $\alpha\geq 0$, $\beta\geq 0$, we have $\hat{G}=\{\{G^{h}\}\}+\alpha[[\mathbf{F}^{h}]]\cdot\mathbf{n}^{+}$, $\hat{\mathbf{F}}\cdot\mathbf{n}^{+}=\{\{\mathbf{F}^{h}\}\}\cdot\mathbf{n}^{+}+\beta[[G^{h}]]$ with

$$-\Big{[}\Big{[}\mathbf{F}^{h}\hat{G}+G^{h}\big{(}\hat{\mathbf{F}}^{h}-\mathbf{F}^{h}\big{)}\Big{]}\Big{]}\cdot\mathbf{n}^{+}=-\alpha\big{(}[[\mathbf{F}^{h}]]\cdot\mathbf{n}^{+}\big{)}^{2}-\beta[[G^{h}]]^{2}\leq 0\text{.}$$

(77)

The energy/entropy is dissipated $\mathcal{E}_{t}\leq 0$.

For the linearised RSWE with constant $f$ detailed in section 2.3 we can preserve geostrophic balance by ensuring that the projection of any continuous stream function $\psi$ into the discrete space $S$ is in $S\cap H^{1}(\Omega)$, and that the resulting initial conditions satisfy $D^{h}\in S\cap H^{1}(\Omega)$ and $\mathbf{u}^{h}\in V\cap H(div)(\Omega)$ and remain in these subspsaces for the discrete linear system. For the linearised RSWE our method becomes

$$\langle\mathbf{w}^{h},\mathbf{u}^{h}_{t}\rangle_{\Omega^{q}}+\langle\mathbf{w}^{h},f\mathbf{k}\times\mathbf{u}^{h}+g\nabla D^{h}\rangle+\langle g(\hat{D}-D^{h})\mathbf{w},\mathbf{n}\rangle_{\partial\Omega^{q}}=0\text{,}$$

(78)

$$\langle\phi^{h},D^{h}_{t}\rangle_{\Omega^{q}}+\langle\phi^{h},H\nabla\cdot\mathbf{u}^{h}\rangle_{\Omega^{q}}+\langle H(\hat{\mathbf{u}}-\mathbf{u}^{h}),\mathbf{n}\rangle_{\partial\Omega^{q}}=0\text{.}$$

(79)

Analogous to the continuous case for any continuous stream function $\psi$ the discrete initial condition $D^{h}=-\tfrac{f}{g}\psi^{h}$, $\mathbf{u}^{h}=\nabla_{d}\crossproduct\psi^{h}\mathbf{k}$ where $\psi^{h}=\Pi_{S}(\psi)$ is an exact discrete steady state. To prove this we show that for these particular initial conditions both the interior and boundary terms vanish.

To see that the interior terms vanish we apply the div-curl cancellation property yielding

$$\nabla_{d}\cdot\mathbf{u}^{h}=\nabla_{d}\cdot\nabla_{d}\crossproduct\psi^{h}\mathbf{k}=0\text{,}$$

(80)

and by construction

$$\begin{split}\langle\mathbf{w}^{h},f\mathbf{k}\crossproduct\mathbf{u}^{h}+g\nabla_{d}D^{h}\rangle_{\Omega^{q}}&=\langle\mathbf{w}^{h},f\mathbf{k}\crossproduct\nabla_{d}\crossproduct\psi^{h}\mathbf{k}-f\nabla_{d}\psi^{h}\rangle_{\Omega^{q}}\\ &=\langle\mathbf{w}^{h},f\nabla_{d}\psi^{h}-f\nabla_{d}\psi^{h}\rangle_{\Omega^{q}}=0\text{.}\end{split}$$

(81)

To show that the boundary terms are $0$ we use the continuity property described in section 3.5.3 to show that $D^{h}$ and the normal components of $\mathbf{u}^{h}$ are continuous across element boundaries, and therefore $\hat{D}-D^{h}=0$ and $(\hat{\mathbf{u}}-\mathbf{u}^{h})\cdot\mathbf{n}=0$. To see that $D^{h}$ is continuous across element boundaries note that continuity of $\psi$ implies continuity of $\psi^{h}=\Pi_{S}(\psi)$, and therefore

$$D=-\dfrac{f}{g}\psi^{h}=-\dfrac{f}{g}\Pi_{S}(\psi)\in S\cap H^{1}(\Omega)\text{.}$$

(82)

is continuous as well. Next we show that $\mathbf{u}$ is continuous in the normal direction. As $\psi^{h}$ is continuous it’s derivative in the tangent direction $\big{(}\nabla^{d}\psi^{h}\big{)}\cdot\mathbf{t}$ is also continuous across element boundaries. The initial flow is given by

$$\mathbf{u}=\nabla_{d}\crossproduct\psi^{h}\mathbf{k}=\mathbf{k}\crossproduct\nabla_{d}\psi^{h}\text{,}$$

(83)

the flow in the normal direction is then

$$\mathbf{u}\cdot\mathbf{n}=\mathbf{k}\crossproduct\nabla_{d}\psi^{h}\cdot\mathbf{n}=\nabla_{d}\psi^{h}\cdot(\mathbf{n}\crossproduct\mathbf{k})=-\nabla_{d}\psi^{h}\cdot\mathbf{t}\text{,}$$

(84)

and is therefore continuous across element boundaries. Therefore both the interior and element boundary terms are zero, and hence for any geostrophic mode there is a corresponding discrete solution which is exactly steady and will be preserved up to machine precision. We note that this result is dependent on the use of a collocated tensor product basis.