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

We develop a novel and efficient discontinuous Galerkin spectral element method (DG-SEM) for the spherical rotating shallow water equations in vector invariant form. We prove that the DG-SEM is energy stable, and discretely conserves mass, vorticity, and linear geostrophic balance on general curvlinear meshes. These theoretical results are possible due to our novel entropy stable numerical DG fluxes for the shallow water equations in vector invariant form. We experimentally verify these results on a cubed sphere mesh. Additionally, we show that our method is robust, that is can be run stably without any dissipation. The entropy stable fluxes are sufficient to control the grid scale noise generated by geostrophic turbulence without the need for artificial stabilisation.


Introduction
Variational methods such as mixed finite element (MFE) methods, and continuous Galerkin and discontinuous Galerkin spectral element methods (CG-SEM and DG-SEM respectively), have become an increasingly appealing choice for atmospheric modelling, owing to their scalable performance on modern supercomputers and ability to discretely preserve conservation properties [McRae and Cotter, 2014, Lee et al., 2018, Taylor and Fournier, 2010, Gassner et al., 2016, Waruszewski et al., 2022].In this paper we present an entropy stable DG-SEM method for the rotating shallow water equations (RSWE) in vector invariant form that we argue combines many of the strengths of these three approaches.
Ideally numerical methods should mimic the properties of the continuous system they approximate, but continuous systems have infinite invariants of which only a finite number can be preserved by any discretisation [Thuburn, 2008].This necessitates choosing a subset of the invariants, based on the particular modelling application, for the discrete model to conserve.For the RSWE in the context of atmospheric modelling we target the following: • Local conservation of mass.Arguably the most important of these properties, it helps to limit the bias of long climate simulations [Thuburn, 2008].
• Local conservation of vorticity.Mid-latitude weather is dominated by vorticity dynamics and its conservation assists in accurately modelling these dynamics.[Lee et al., 2018] Additionally this prevents the gravitational potential, by far the largest component of atmospheric energy, from spuriously forcing the vorticity which may destroy the meteorological signal [Staniforth and Thuburn, 2012].
• Local energy conservation and stability.The vast majority of energy dissipation within the atmosphere is due to three-dimensional effects not included in the RSWE.However there is evidence of a small amount of energy dissipation from two-dimensional layer wise effects on the order or 0.1W m −2 [Thuburn, 2008].Additionally the RSWE conserve energy and so discretisations of the RSWE should have minimal or no energy dissipation.Energy is also the entropy of the RSWE, and we follow the approach of using entropy stable fluxes to ensure non-linear stability [Gassner et al., 2016][Waruszewski et al., 2022][Carpenter et al., 2014].
• Exact steady discrete geostrophic balance.The atmosphere is often characterised by slow moving variations around a linear state, consisting of a collection of slowly evolving geostrophic modes.To accurately represent this process a discretisation of the RSWE applied to the linear RSWE should contain exactly steady discrete geostrophic modes.[Cotter andShipton, 2012] [Lee et al., 2018] The MFE approach is able to locally conserve mass, has steady discrete geostrophic modes, and globally conserves vorticity and energy [Cotter and Shipton, 2012, McRae and Cotter, 2014, Lee et al., 2018].MFE methods can also conserve potential enstrophy on affine geometries subject to exact integration of the chain rule, but this has not yet been extended to spherical geometry [Lee et al., 2018, McRae andCotter, 2014].Impressively, these schemes can be run stably without any dissipation [Lee et al., 2022].To simultaneously achieve all these properties MFE methods use a non-diagonal mass matrix making them more computationally expensive than alternative spectral element methods [Lee et al., 2018].
In parallel to this several groups have developed highly scalable collocated CG-SEM and DG-SEM methods for the RSWE [Taylor andFournier, 2010] [Gassner et al., 2016].CG-SEM applied to the vector invariant form of the RSWE has been shown to locally conserve mass, vorticity, and energy but has not been shown to possess steady discrete geostrophic modes.Due to dispersion errors CG-SEM requires filtering or artificial viscosity which may introduce undesirable numerical artefacts, excessively dissipate energy, break vorticity conservation, and adversely effect the accuracy of the solution [Melvin et al., 2012, Dennis et al., 2012].
To date DG-SEM has only been applied to the conservative form of the RSWE and therefore has only been shown to locally conserve mass and energy [Gassner et al., 2016].To the best of the authors' knowledge there is no DG-SEM which has steady discrete geostrophic modes or conserves vorticity, this is a direct result of solving the RSWE in conservative form as opposed to the vector invariant form.Similar to CG-SEM, DG-SEM has dispersion errors and therefore requires that high frequency waves be damped.In contrast to CG-SEM, DG-SEM can achieve this through entropy stable fluxes which avoids the pitfalls of artificial viscosity.
In this work we present a DG-SEM for the RSWE in vector invariant form that is arbitrarily high order in space, locally conserves mass and vorticity, locally semi-discretely conserves energy, and preserves linear geostrophic balance.We then show how to consistently modify this scheme to dissipate energy/entropy through the use of dissipative fluxes.
The rest of the article is as follows.Section 2 presents the continuous equations and derives the properties we target with our numerical method.Section 3 introduces our spatial discretisation and presents proofs for the discrete properties.Section 4 presents the results of numerical experiments that verify our theoretical analysis.

Continuous shallow water equations
In this section we introduce the RSWE and derive the continuous properties that our numerical scheme should discretely mimic.

Vector invariant equations
We use the vector invariant form of the rotating shallow water equations on a 2D manifold Ω 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 where t ≥ 0 denotes the time variable, the unknowns are the flow velocity vector u = [u, v, w] T and the fluid depth D. Furthermore, ω is the absolute vorticity, f is the Coriolis frequency, F is the mass flux, G is the potential, and k is normal of the manifold.At t = 0, we augment the SWE (1)-( 2) with the initial conditions such that the flow is constrained to the manifold u 0 • k = 0.

Mass and vorticity conservation
Integrating the mass equation (2) over a periodic domain yields conservation of mass.Similarly, we can derive a continuity equation for vorticity by applying k • ∇× to the momentum equation (1) giving Integrating this over a periodic domain yields conservation of vorticity.

Energy conservation
The energy E = 1 2 Du • u + 1 2 gD 2 also satisfies a continuity equation.It's time derivative is ) where we have used F • ωk × u = ωDu • k × u = 0. Therefore the energy evolves according to the continuity equation (11) and global conservation again follows from integration over a periodic domain.

Linear geostrophic balance
In addition to a set of inerto-gravity wave solutions, the linear rotating shallow water equations also exhibit steady geostrophic modes [Cotter andShipton, 2012, Vallis, 2017].The shallow water equations linearised around a state of constant mean height H and zero mean flow are The steady solutions to these equations can be found by choosing an arbitrary stream function ψ and setting

Entropy
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 |u| < √ gD.Atmospheric flows are subcritical and so this restriction does not cause any difficulty in practice.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 h.The Hessian of the energy with respect to the prognostic variables is and its characteristic equation is As h, g > 0 the eigenvalues corresponding to the left bracket are always positive, and as long as hg − |u| 2 > 0, which implies |u| < √ 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.

Spatial discretisation
In this section we describe the spatial discretisation of our method and prove several discrete stability and conservation properties.Our method decomposes the domain into quadrilateral elements Ω k and approximates the solution in each quadrilateral by a polynomial.

Function spaces
We construct our function spaces in typical DG-SEM fashion by mapping each element Ω k to the reference square [−1, 1] 2 with map r(x; k), and restricting functions to those mapped to polynomials of order m by r(x; k).For computational efficiency we use a tensor product Lagrange polynomial basis with interpolation points collocated with Gauss-Lobatto-Legendre (GLL) quadrature points.Using this basis the polynomials of up to order m on the reference square can be defined as where ξ, η ∈ [−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 Note that within each element functions φ ∈ S can be expressed as where ξ, η = r(x; k), φ ijk = φ(r −1 (ξ i , η j ; k)).
We define a discontinuous vector space which contains the velocity u.Let v 1 (x, y, z) and v 2 (x, y, z) be any independent set of vectors which span the tangent space of the manifold at each point, then the discontinuous vector space V containing the velocity u can be defined as We note that V is independent of the particular choice of v i however two convenient choices are the contravariant and covariant basis vectors.

Covariant and contravariant 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 vectors g 1 and g 2 are defined as with these the determinant of the Jacobian of r(x; k) can be expressed as The contravariant vectors are given by g 1 = ∇ξ, g 2 = ∇η, (22) which have the useful property that g α • g β = δ αβ , (23) enabling vectors to to expressed as w = α=1,2 where For flows such that w • k = 0 this enables the divergence, gradient, and curl to be expressed as where we have only used ∇× on quantities either parallel or orthogonal to k to simplify the exposition.

Discrete inner product and boundary integrals
We approximate integrals over the elements by using GLL quadrature, this defines a discrete element inner product which approximates the integral The discrete global inner product is defined simply as f, g We also approximate element boundary integrals using GLL quadrature where and n is the unit normal of the element boundary ∂Ω k given by The unit tangent vector t is similarly defined as With the discrete inner products defined we now construct discrete projection operators Π(•) S and Π(•) V which preserve discrete L 2 inner products in S and V respectively.For example Π(•) S is defined as the unique projection that satisfies for all functions b.Π(•) V is defined similarly.In practice the discontinuous projection operators Π S and Π V simply interpolate functions within each element.

Properties of function spaces
Here we present the properties of the function spaces and discrete operators that we use in our stability and conservation proofs.

Element local summation by parts
The discrete gradient and divergence operators satisfy an element local summation by parts (SBP) property This allows us to transform between the weak and strong form of our DG prognostic equations which we use in several proofs.

Divergence curl cancellation
The discrete curl operator satisfies a discrete grad − curl cancellation which we use for vorticity conservation.

Continuity of Interpolants
For tensor product bases with boundary interpolation points, the element local interpolation of continuous functions ψ h = Π S (ψ) are continuous across element boundaries.To see this consider two elements k 1 and k 2 that share a boundary at ξ = 1 and ξ = −1 respectively.The value of the interpolant ψ h in element k 1 at the shared boundary is and in element k 2 at that same boundary The mappings r(•; k) agree along element boundaries, therefore r −1 (1, η; k 1 ) = r −1 (−1, η; k 2 ) and implying that ψ h is continuous across element boundaries.

Curl gradient identity
We also note that which we use in our discrete linear geostrophic balance proof.
(47) This arises from a 1D SBP property for integrals along each side of the quadrilateral.Consider the ξ = 1 boundary, and therefore by the exact integration of 2m − 1 polynomials with GLL quadrature (48)

Discrete formulation
First we introduce notation for jumps [[•]] and averages {{•}} across element boundaries as {{a}} = 1 2 a + + a − , ( 49) where • + and • − 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 where the numerical fluxes are defined for the scalars α, β ≥ 0 as (54) and the discrete vorticity is ) where t is the unit tangent vector along element boundaries.
Our scheme is discretely entropy stable so long as α, β ≥ 0, and is locally semi-discretely energy/entropy conserving for α = β = 0 and locally discretely energy/entropy dissipating for α, β ≥ 0. We use a centred mass flux β = 0 to ensure that potential energy is not dissipated, and either use centred fluxes α = 0 or α = g 2c and β = 0 where is the fastest wave speed on either side of the boundary.This ensures that Ĝ has the correct units and reduces to a Rusanov flux [Rusanov, 1970] in the case of the linearised SWE in section 2.5.

Local conservation of mass
Local conservation of mass can be shown using the element local SBP property to transform the mass equation ( 2 57) and then choosing the test function φ h = 1 to yield the rate of change of the total mass of element k (58) By construction the numerical fluxes are continuous across element boundaries and so the mass flowing through each point along the 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

Local conservation of vorticity
We show that the potential G does not spuriously force the vorticity, and that the vorticity is also locally conserved.This is the result of our choice of Ĝ and ω h , and the element boundary SBP property discussed in section 3.5.5.
The evolution of the vorticity follows from taking the time derivative of ( 55) First we expand the volume term by using the weak form of ( 51), the curl gradient identity (46 Next we evaluate the boundary term, as this is a GLL collocated method substituting the test function l i (ξ)l j (η)t ij into (1) yields and therefore With these the vorticity evolution becomes where we have used the definition of Ĝ.Using the tangent boundary SBP property in section 3.5.5 this becomes demonstrating that the potential G does not spuriously force the discrete vorticity.
To show local vorticity conservation we choose φ = 1 to yield By construction {{ω h u h }} is continuous across element boundaries and therefore vorticity is locally conserved.Global conservation in a periodic domain follows from cancellation of {{ω h u h }} • n across element boundaries.

Local energy conservation and stability
Here we show that energy is semi-discretely locally conserved for centred flux α = β = 0, and that energy is semidiscretely dissipated for α, β > 0. We obtain fully discrete energy stability by using a strong stability preserving (SSP) time integrator.
We define the discrete energy of a single element E k to be with the total energy given by the sum over all elements E = k E k .Following the continuous form of the energy conservation analysis in section 2.3 the semi-discrete time evolution of E k is given by where we have used the weak form of the mass equation and that 0 only holds for the special case of diagonal mass matrices.
To show E t = k E k t ≤ 0 it suffices to show that the jump across element boundaries for all points.Using the expansion ]{{b}} the energy jump condition becomes For centred fluxes Ĝ = {{G h }} and F = {{F}} the jump condition is 0 and therefore the method is semi-discretely locally energy conserving.The jump condition is ≤ 0 and dissipates energy/entropy for Ĝ where α, β ≥ 0. To see this we expand (75) yielding

Discrete steady linear modes
We can preserve geostrophic balance by ensuring that the projection of any continuous stream function ψ into our discrete space S is in S ∩ H 1 (Ω), and that the resulting initial conditions satisfy and remain in these subspsaces for the discrete linear system.For the linearised SWE detailed in section 2.5 our method becomes Analogous to the continuous case for any continuous stream function ψ the discrete initial condition D h = − f g ψ h , u h = ∇ d × ψ h k where ψ h = Π S (ψ) is an exact discrete steady state.To prove this we show that for these particular initial conditions both the interior and boundary terms are 0.
To see that the interior terms are 0 we apply the div-curl cancellation property yielding and by construction 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 u h are continuous across element boundaries, and therefore D − D h = 0 and (û − u h ) • n = 0.
To see that D h is continuous across element boundaries note that continuity of ψ implies continuity of ψ h = Π S (ψ), and therefore is continuous as well.Next we show that u is continuous in the normal direction.As ψ h is continuous it's derivative in the tangent direction ∇ d ψ h • t is also continuous across element boundaries.The initial flow is given by the flow in the normal direction is then 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 heavily dependant on the use of a collocated tensor product basis.

Numerical results
In this section we present the results of numerical experiments which verify our theoretical results.For all experiments we use an equi-angular cubed sphere mesh [Rančić et al., 1996], an SSP-RK3 time integrator [Shu and Osher, 1988], and third-order polynomials.Except for the energy convergence test we use an adaptive time step with a fixed CFL of 0.8, ∆t = 4 5 CF L = 4∆x 5c(2p+1) where p is the polynomial order, c is the fastest wave speed, and ∆x is the element spacing.The code for these experiments can found at https://github.com/kieranricardo/dg-swe-code.

Williamson test case 2
We apply both the energy conserving and energy dissipating methods to the steady state Williamson test case 2 with angle α = 0 [Williamson et al., 1992].We use a resolution of 6 × n 2 elements for increasing n = [3,5,10,15,30].
Figure 1 shows that both methods are converging at greater than 3rd order validating the consistency of our scheme.The energy conservative method is converging at approximately 3.4 order, while the dissipative method is converging at 3.8 order.This difference in convergence is most likely due to the energy dissipating method damping spurious high frequency waves.

Galewsky test case
In order to validate the model in a more physically challenging turbulent regime we also simulate the Galewsky test case of a geostrophically balanced zonal jet destabilised by a Gaussian perturbation in the depth field [Galewsky et al., 2004].

Semi-discrete energy conservation
To test the semi-discrete energy conservation of our energy conserving scheme we ran the Galewsky test case at a low resolution of 6 × 5 × 5 elements for decreasing time steps ∆t = [50s, 40s, 30s, 20s, 10s].As expected for a 3rd order time integrator figure 2 shows that the energy conservation error decreases at 3rd order, verifying that the energy conserving scheme is semi-discretely energy conserving.

High resolution simulation
To verify the accuracy and robustness of our method we ran the Galewsky test case for 20 days at a resolution of 6 × 64 × 64 elements.Figure 3 shows that both energy conserving and energy dissipating methods are energy stable, conserve mass and vorticity, and appear to reasonably control potential enstrophy.Notably even our energy conserving method, which contains no dissipation of any kind, runs stably for the 20 day period.However, potential enstrophy is steadily increasing and the simulation may eventually become unstable.Figure 4 shows the relative vorticity at day 7 where the characteristic Kelvin-Helmholtz billows can be observed.At day 7 some noise can be seen in the energy conserving solution, most likely due to spurious high frequency waves.This is typical for SEM without dissipation but our entropy stable fluxes are able to damp these spurious numerical modes.To gain an indication as to whether the energy dissipating method is over damping we scale the final relative energy loss by the total energy in the Earth's atmosphere to obtain an adjusted average energy dissipation rate of 0.01W m −2 , an order of magnitude lower than the energy loss of the real atmosphere at the modelled length scales [Thuburn, 2008].

Discrete steady states
To verify that our scheme supports discrete steady modes we use the stream function ψ = 0.1 cos(λ) cos(θ), where λ and θ are latitude and longitude respectively, with 6 × 5 × 5 elements, a spatially constant Coriolis parameter f = 8, H = 0.2, and g = 8.Following the analysis in section 3.10 we use the initial condition u h = ∇ × Π S (ψ)k and D h = − f g Π S (ψ) and integrate forward in time.Figure 5 shows the time L 2 difference from the initial conditions, demonstrating that the initial state is preserved up to machine precision.

Conclusion
In this paper we presented a DG-SEM for the rotating shallow water equations in vector invariant form that possess many of the conservation properties and exact discrete linear geostrophic balance of recent mixed finite element methods.This enables our method to have the computational efficiency of SEM while having many of the mimetic properties of mixed finite element methods.We proved that on a general curvilinear grid our method locally conserves mass, vorticity, and energy, does not spuriously project the potential onto vorticity, and satisfies a discrete linear geostrophic balance.We then experimentally verified these properties and show that geostrophic turbulence can be well simulated with our entropy stable fluxes and no additional stabilisation.
While the issue of computational modes was not addressed, for the linear SWE our method reduces to a DG method with Rusanov fluxes which is well known to have good dispersion properties.Future work will focus on analysing the numerical modes of our method, investigating alternative numerical fluxes, and extending this method to the thermal shallow water equations.

Figure 1 :
Figure 1: Convergence plots for the energy dissipating and energy conserving methods applied to Williamson test case 2.

Figure 2 :
Figure 2: Energy conservation error for the centred flux scheme at low resolution for the Galewsky test case

Figure 3 :
Figure 3: Conservation errors for the Galewsky test case at high resolution.

Figure 4 :
Figure 4: Relative vorticity for the Galewsky test case at day 7 for the energy conserving (left), and energy dissipating methods (right).

Figure 5 :
Figure 5: L2 difference from discrete initial conditions for the linear geostrophic balance test case.