Energy-enstrophy conserving compatible finite element schemes for the rotating shallow water equations with slip boundary conditions

We describe an energy-enstrophy conserving discretisation for the rotating shallow water equations with slip boundary conditions. This relaxes the assumption of boundary-free domains (periodic solutions or the surface of a sphere, for example) in the energy-enstrophy conserving formulation of McRae and Cotter (2014). This discretisation requires extra prognostic vorticity variables on the boundary in addition to the prognostic velocity and layer depth variables. The energy-enstrophy conservation properties hold for any appropriate set of compatible finite element spaces defined on arbitrary meshes with arbitrary boundaries. We demonstrate the conservation properties of the scheme with numerical solutions on a rotating hemisphere.


Introduction
For large scale balanced flows, energy and enstrophy are important quantities for the rotating shallow water equations due to the cascade of energy to large scales whilst enstrophy cascades to small scales. At the level of numerical discretisations, energy conservation becomes important over long time integrations, whilst enstrophy conservation (or dissipation at the small scale) provides control of the regularity of the velocity field over long times.
Energy and enstrophy conserving schemes for the rotating shallow water equations have a long history that goes back to Arakawa and Lamb (1981); Sadourny (1975). The finite difference schemes in these papers were constructed from two important ingredients: (1) the vector-invariant form of the equations, and (2) the use of staggered grid finite difference methods built around discretisations of the div, grad and curl operators that preserve the vanishing div-curl and curl-grad identities at the discrete level. These discretisations form the foundations of several operational weather, ocean and climate models that are in current use. Another important practical aspect is that discretisations should preserve stationary geostrophic modes when applied to the f -plane linearisation of the shallow water equations. Ringler et al. (2010) addressed the issue of extending these properties to C-grid staggered finite difference discretisations on unstructured orthogonal grids, describing separate energy-conserving and enstrophy-conserving schemes; Thuburn and Cotter (2012) extended these ideas to non-orthogonal grids, making use of ideas from discrete exterior calculus (Hirani, 2003). Ringler et al. (2010) also considered enstrophy dissipation through the Anticipated Potential Vorticity method, following the structured rectangular grid formulation of Arakawa and Hsu (1990). There is still no known closed form for an energy-enstrophy conserving C-grid formulation on unstructured grids with an f -plane linearisation that preserves stationary geostrophic modes, but Eldred and Randall (2017) showed that such schemes can be obtained computationally through numerical optimisation.
In a series of papers, Salmon (2004Salmon ( , 2005Salmon ( , 2007, Salmon showed how to use Poisson and Nambu brackets to build conservation into numerical discretisations. For example, Stewart and Dellar (2016) provided a C-grid discretisation for the multi-layer shallow-water equations with complete Coriolis force. The variational formulation finite element method makes it easier to mimic the Poisson bracket structure of the vector-invariant shallow water equations at the discrete level, whilst compatible finite element spaces replicate the div-curl and curl-grad identities in the discrete setting. McRae and Cotter (2014) showed that this leads to a natural energy-enstrophy conserving compatible finite element scheme, with the bracket structure being exposed in the appendix. The finite element exterior calculus framework underpinning these properties was exposed by Cotter and Thuburn (2014). The same structure has been exploited to produce energy-enstrophy conserving discretisations using more exotic finite element spaces. Eldred et al. (2016) constructed compatible spaces from splines that allow higher-order approximations constructed around the low-order C-grid data structure, and Lee et al. (2017) used mimetic spectral elements. In the context of imcompressible two-dimensional turbulence, Natale and Cotter (2017) considered consistent energy-conserving/enstrophy-dissipating finite element schemes, including a formulation that extends to a consistent energy-conserving enstrophy-dissipating version of the McRae and Cotter (2014) scheme, and showed that these schemes have favourable turbulent backscatter properties.
One aspect missing from this framework is the treatment of lateral boundaries. These are necessary for ocean modelling in the presence of coastlines, and also in the extension to 3D vorticity conserving schemes in the atmosphere (since there is a boundary at the Earth's surface). In this paper, we are considering the inviscid equations with slip boundaries conditions. The addition of dissipative terms will introduce further boundary conditions. We do not discuss in-flow or out-flow boundary conditions. Arakawa and Hsu (1990) avoided the consideration of boundary conditions by considering layers that taper to zero depth, allowing wetting and drying near coastlines as well as allowing layers to outcrop the top surface in multilayer models. This approach has underpinned the formulation of isopycnal ocean models (Hallberg and Rhines, 1996). More recently, (Ketefian and Jacobson, 2009) produced an energy-enstrophy conserving discretisation using the C-grid approach on structured meshes with boundaries, by making the key observation that since vorticity cannot be diagnosed from velocity and height at boundary points, then vorticity on the boundary should be treated as a prognostic variable with its own conservation equation. Salmon (2009) used a similar idea but blended between a vorticity-divergence-depth and velocity-depth formulation. In this paper we show that the introduction of vorticity degrees of freedom on the boundary can also lead to an energy-enstrophy conserving formulation in the compatible finite element setting.
The rest of the paper is structured as follows. In Section 2 we review the compatible finite element energyenstrophy conserving formulation, to motivate the issues relating to boundaries. We then describe the new scheme that incorporates vorticity degrees of freedom on the boundary in order to recover energy-enstrophy conservation when boundaries are present. In Section 3 we demonstrate this scheme with numerical results. Finally in Section 4 we provide a summary and outlook.

Review of the energy-enstrophy conserving formulation on domains without boundaries
The energy-enstrophy conserving formulation in McRae and Cotter (2014) starts from the rotating shallow-water equations in vector-invariant form for velocity u and layer depth D, where where F is the mass flux, and q is the potential vorticity, defined by and where F ⊥ = k × F , k is the unit normal to the domain surface, and g is the acceleration due to gravity. If we consider a fluid flow where the fluid is moving in a container Ω such that now fluid can enter or leave, the relevant boundary conditions are u ⋅ n = F ⋅ n = 0, on ∂Ω.
There are then no boundary conditions for advected quantities such as the layer depth D, since there are no advective fluxes through ∂Ω. Inviscid fluid equations such as the rotating shallow water equations can be derived from Hamilton's principle considering a fluid flow with slip boundary conditions, leading to Poisson bracket formulations; for more details, see Holm et al. (1998).
In the absence of boundaries (periodic solutions or the surface of a sphere being the main relevant cases), we can introduce a weak formulation of Equations (1-4) as follows.
We introduce the compatible finite element spaces (V 0 , V 1 , V 2 ), that form a discrete de Rham sequence, with commuting, bounded, surjective projections (π 0 , π 1 , π 2 ). For the examples in this paper we have concentrated in the spaces V 0 = CG k , V 1 = BDM k−1 , V 2 = DG k−2 . These are then used to formulate a Galerkin finite element discretisation of Equations (6-9).
Definition 2. The compatible finite element discretisation of the weak formulation in Definition 1 seeks Remark 3. Since D t and ∇ ⋅ F δ are both in V 2 , Equation (12) is equivalent to the equation D t + ∇ ⋅ F δ = 0 holding in L 2 (Ω).

Remark 4.
It is important to note that q δ and F δ are merely diagnostic variables that can be computed from the prognostic variables u δ and D δ at any time.
These equations have an equivalent (almost 1 ) Poisson bracket formulation.
Definition 5 (Almost Poisson bracket formulation). Let H ∶ V 1 × V 2 → R be the Hamiltonian functional, defined by We define the bracket {⋅, ⋅} by Then the corresponding almost Poisson bracket formulation defines dynamics for any functional where F δ and q δ are considered to be functions of u δ and D δ defined by equations (12) and (14) respectively.
Proof. This follows immediately from the anti-symmetry of the bracket, Remark 8. The almost Poisson formulation has an equivalent strong form formulation, given by ∂ ∂t where J is an (u δ , D δ )-dependent skew-adjoint structure operator defined by and This will be useful for describing energy-conserving time integration methods later.
Equations (11-14) also conserve total potential vorticity Q and enstrophy Z, defined by This follows directly from the implied potential vorticity equation as shown in the next proposition.
Proposition 9. Let (u δ , D δ , F δ , q δ ) solve the compatible finite element discretisation in Definition 2. Then q δ satisfies the equation Proof. Taking the time derivative of (14), we obtain If w = −∇ ⊥ γ, then w ∈ V 1 , and we may use it in (11) to obtain Combining these two equations gives the result. (28) is the standard Galerkin discretisation of the conservation law
Since 1 ∈ V 0 , we may use γ = 1 in Equation (28), to geṫ Similarly, q δ ∈ V 0 , so we may use γ = q δ to geṫ Remark 12. Similar calculations show that Q and Z are Casimirs of the almost Poisson bracket, i.e.
for any functional F .

Energy-enstrophy conserving formulation on domains with boundaries
We now consider the case of slip boundary conditions on ∂Ω. This requires us to consider a modified de Rham complexH whereH and where Tr ∂Ω is the trace operator returning functions defined in L 2 (∂Ω). The presence of ∂Ω requires the modification of (14) to include a boundary integral, where ⟪⋅, ⋅⟫ defines the L 2 inner product on ∂Ω, If we just apply this modification to the discretisation in Definition 2, replacing V 1 withV 1 and V 2 withV 2 , then we still have an almost Poisson bracket formulation (so energy is still conserved). Proposition 9 does not hold though, because we can only take w = ∇ ⊥ γ in (30) if γ ∈V 0 , but we need to be able to take γ ∈ V 0 . In particular, the test functions 1 and q δ required to show conservation of total potential vorticity and enstrophy are not inV 0 (at least, not in general for q δ ). In numerical experiments with this formulation we do indeed see unbounded growth in enstrophy due to sources at the boundary which eventually pollute the solution throughout the domain.
To resolve this problem, we introduce a vorticity variable Z δ ∈ V 0 , such that The projectionZ 0 of Z δ intoV 0 , is thus a diagnostic variable that can be obtained from u δ according to since now ∇ ⊥ γ ∈V 1 . However, to compute Z δ we also need to know its projection Z ′ ontoV ⊥ 0 , defined as the before projecting toV ⊥ 0 to obtain Z ′ . After initialisation, Z ′ has its own dynamics as given below. Definition 13. The compatible finite element discretisation of the rotating shallow water equations seeks The dynamics of Z ′ are precisely chosen so as to recover the correct dynamics for potential vorticity q δ ∈ V 0 (and not just the projection of potential vorticity intoV 0 ), i.e. so that Proposition 17 will hold. This idea is analogous to Ketefian and Jacobson (2009), who introduced vorticity variables at boundary vertices that have dynamics that imply the PV equation at the boundary, which is otherwise undefined.
Since Z ′ is independent of u δ and D δ , we need to enlarge the phase space to obtain the bracket formulation for this extended system to include it.

Definition 14 (Extended almost Poisson bracket formulation). Let
We define the bracket {⋅, ⋅} by Then the corresponding almost Poisson bracket formulation defines dynamics for any functional whereZ, F δ and q δ are considered to be functions of u δ and D δ defined by the equations (54-56) respectively.
Proof. The proof is similar to the proof of Proposition 6, except that we need to also check the dynamics for Z ′ . Taking F = ⟨γ, Z ′ ⟩ for γ ∈V ⊥ 0 , so that δF δZ ′ = γ, we havė as required.
Proof. Follows directly from the almost Poisson bracket formulation.
Proof. The proof follows from the implied PV equation, as for the boundary-free case.
It is inpractical to deal withV ⊥ 0 as there is no local basis. However, an equivalent formulation exists.
We obtain this by time differentiating the definition ofZ as the projection of q δ D δ intoV 0 to obtain as required.
Definition 21 (Poisson brackets for equivalent extended formulation). Let H ∶V 1 × V 2 × V 0 → R be the Hamiltonian functional, defined by We define the bracket {⋅, ⋅} by where q δ is understood to be a function of Z δ and D δ given by and where F δ is considered to be a function of u δ and D δ defined by equation (55). Then the corresponding almost Poisson bracket formulation defines dynamics for any functional Proposition 22. Equations (67-70) imply the bracket formulation in Definition 21.
On the other hand, time differentiating (79) gives hence the result.

Energy-conserving time integration
In this section we demonstrate the numerical scheme, using the energy-preserving Poisson integrator of Cohen and Hairer (2011). Whilst the implicit midpoint rule conserves quadratic Hamiltonians exactly, the energypreserving Poisson integrator extends this property to higher-degree polynomials, such as the shallow-water Hamiltonian which is cubic. Given a Poisson system of the forṁ the Poisson integrator takes the form For our equivalent extended shallow water discretisation, we have ⟨φ, δH δ D δ ⟩ = ⟨φ, g Hence, we obtain the following time discretisation, This is equivalent to solving for (u δ n+1 , where F δ = u δ n+1 D δ n+1 3 + u δ n D δ n+1 6 + u δ n+1 D δ n 6 + u δ n D δ n 3, K = u δ n+1 2 3 + u δ n ⋅ u δ n+1 3 + u δ n 2 3. In this section, we demonstrate the energy-enstrophy conserving scheme through the adaptation of two of the popular shallow water sphere testcases (Williamson et al., 1992) to the hemisphere domain, so that the equator becomes a boundary. In all cases the mesh used is (half of) an octahedral mesh obtained by hierarchical refinement of the triangular faces of an octahedron, before mapping to the sphere so that lines of constant height are mapped to lines of constant latitude (see Figure 1). We use the compatible finite element spaces as follows: P 3 for vorticity, BDM 2 for velocity, P 1 DG for layer depth. All numerical calculations are performed using Firedrake, the automated code generation finite element library (Rathgeber et al., 2016). In all cases we use a sphere of radius R 0 = 6371220m, with rotation rate Ω = 7.292 × 10 −5 s −1 so that f = 2Ωz R 0 , and graviational acceleration g = 9.810616ms −2 .
First we consider a convergence test based on Test Case 2, which is a steady state solid rotation solution, which is also a solution when restricted to the Northern hemisphere, because the velocity is tangential to the equator, and the solution is zonally symmetric. The velocity and height are intialised according to the steady state solution where h 0 = 5960m and u 0 = 2πR 0 (12 days ). The equations are then integrated numerically for 15 days and the solution fields are compared against the initial conditions. We might expect any issues associated with consistency errors at the boundary to manifest themselves as loss of optimal convergence rates in the L 2 norm of these errors; the convergence rates (see Figure 2) match those of the scheme applied to the full sphere. In these calculations the relative energy is conserved to 7 decimal places, although we do observe some small fluctuations below that level of precision which are due to the tolerance settings of the Newton solver. As is expected for a steady state solution, we also observe conservation of relative enstrophy to 6 decimal places even though conservation is not guaranteed by the time integration method. Plots of energy and enstrophy for this experiment are shown in Figure  3.
To verify that the correct boundary condition is being enforced, we ran a Kelvin wave testcase in a disk with constant f . There are no analytic solutions for this problem, either for nonlinear or linear equations, but it is possible to observe that the Kelvin wave propagation mechanism is functioning correctly. In the limit of large radius compared to the Rossby deformation radius and the limit of low amplitude, the solution should approach that of the linear Kelvin wave along a straight boundary, which propagates at speed √ gH. We verified that this is approximately the case by considering a disk of radius 1, H = g = 1 so that the wave speed is 1, and f = 10 so that the deformation radius is 1/10, with initial data where a 0 = 0.01 is the wave amplitude. The results of this testcase are shown in Figure 4.

Energy-conserving, enstrophy-dissipating scheme
For turbulent large scale balanced dynamics, enstrophy conservation is not necessarily desirable since the enstrophy cascade to small scales will mean that the enstrophy will accumulate at the gridscale. Arakawa and Hsu (1990) introduced an anticipated potential vorticity approach, which modifies the equations so that they still conserve energy but have an upwinded/anticipated potential vorticity so that enstrophy is dissipated at small scales. Applying the modification to the energy-enstrophy conserving discretisation of Arakawa and Lamb (1981) preserves these properties at the discrete level. This approach has been followed by a number of authors extending the C-grid approach to unstructured grids and compatible finite element methods. In this paper we modify this approach using a streamline-upwind Petrov-Galerkin (SUPG) method for the implied potential vorticity equation. Unlike the anticipated potential vorticity method, this remains a consistent scheme to the equations without dissipation, but still introduces enstrophy dissipation at the gridscale. This approach was originally advocated in Cotter and Thuburn (2014);McRae (2015). The SUPG method, surveyed in Hughes (1995), replaces the test function in the weak formulation by a modified test function that is biased in the upwind direction. Since this test function is applied to the entire equation, this does not alter the residual formulation, only the nature of the test functions, and hence the scheme is expected to remain consistent at the appropriate order. SUPG was first applied to the Euler equations in streamfunction-vorticity formulation by Tezduyar et al. (1988); Tezduyar (1989), and the multiscale behaviour of the resulting scheme was examined by Natale and Cotter (2017).  The SUPG modification of the energy-enstrophy conserved shallow water scheme of this paper is obtained by replacing q δ in Equation (51) by q * given by where τ is a chosen time parameter. Equivalently, making use of the fact that D δ t + ∇ ⋅ F δ = 0 in L 2 , we have If q δ , F δ and D δ are replaced with the exact solutions q, F and D respectively, this extra term vanishes, because the residual of the equation is zero. We see that this is a consistent modification to the equations. Further, since this replacement is equivalent to replacing q δ with q * in the Poisson bracket from Definition 14, we see that the modified formulation still conserves energy.
After manipulations identical to the previous section, the equivalent PV equation is which rearranges to This is the SUPG discretisation of the potential vorticity conservation equation. Selecting γ = 1 implies that the total potential vorticity is still conserved. After using (101) to rewrite the SUPG formulation as we see that there is the possibility for enstrophy dissipation since the term is positive semi-definite, i.e. setting γ = q δ gives a non-negative number. Indeed, this is the term that arises in the anticipated potential vorticity method translated to compatible finite element methods in McRae and Cotter (2014). The other SUPG term ⟨τ F δ ⋅∇γ, q δ t ⟩ = 0, is sign indefinite, and so τ needs to be sufficiently large to guarantee monotonic decay in Z. The presence of this term means that the scheme is consistent. When the solution is well-resolved, there will be no dissipation of enstrophy, which only becomes significant once the solution becomes marginally-resolved. In this numerical example, we used a semi-implicit timestepping formulation for a more efficient implementation. Our semi-implicit scheme can be thought of as a fixed number of Picard iterations towards the energy-conserving time integrator given above. Writing Equations (88-89) in the form where any dependency on q δ n+1 or F δ n+1 2 is obtained by solving Equations (93) and (94) using u δ n+1 and D δ n+1 . Then to implement the timestep, we perform some number (4, in the case of this example) of fixed point iterations as follows. In the first iteration we take the initial guess u δ n+1 = u δ n , D δ n+1 = D δ n . Then we solve for ∆u ∈V 1 , ∆D ∈ V 2 such that Then we update u δ n+1 ← u δ n+1 +∆u, D δ n+1 ← D δ n+1 +∆D. This completes one fixed point iteration.
The system (108-109) can be hybridised and statically condensed to obtain a sparse discrete Helmholtz problem which we assembled using the Slate subpackage of Firedrake (Gibson et al., 2018). We used τ = ∆t, with ∆t =50s on mesh refinement level 5, where we observed monotonic decay in the enstrophy.
In this example, we used the initial condition of test case 5 from Williamson et al. (1992), in which a balanced flow in solid rotation is disturbed by the sudden appearance of a conical mountain at time 0. The difference in this example is that the equations are solved on the hemisphere, not the sphere. This means that a comparison cannot be made with other integrations, but we can check the stability of the scheme and observe the energy and enstrophy behaviour for long integration times once small scale vorticity filaments appear.
Snapshots of the solution are shown in Figure 5. These solutions show the formation of vortex filaments that roll up and stretch until they reach the grid scale, where they become more noisy (but do not pollute the whole domain with oscillations as they do in the enstrophy-conserving case). Plots of energy and enstrophy are shown in Figure 6. As expected we see approximate energy conservation over all times, and enstrophy eventually starts to decay when gridscale features appear. We also observed conservation of total potential vorticity up to round-off error.
To make a more careful study of the energy conservation, we performed a convergence experiment where we ran the same mountain test case for 15 days at mesh refinement level 6, to obtain a more interesting flow. We then used this as the initial condition for a finite time interval convergence test checking for energy conservation as ∆t goes to zero, verifying that the spatial discretisation is exactly energy-conserving. This is demonstrated in Figure  7.

Summary and outlook
In this paper we provided a compatible finite element scheme for the rotating shallow water equations in the presence of boundaries. The prognostic variables are velocity and height, plus vorticity on the boundary. The scheme has an equivalent formulation where one discards the potential vorticity -velocity relationship (which is then a consequence of the subsequent equations) in favour of a potential vorticity evolution equation everywhere in the domain. This becomes reminiscent of the energy-enstrophy mimetic spectral element conserving formulation of Palha and Gerritsma (2017), in which both vorticity and velocity are maintained as prognostic variables and a timestaggered discretisation leads to efficient implementation with energy and enstrophy conservation, although in that case the correspondence between velocity and vorticity is not guaranteed by the timestepping scheme. In our numerical experiments we used an implicit energy-conserving timestepping scheme, which preserves the correspondence between velocity and vorticity, but does not preserve enstrophy exactly. This scheme was implemented via Newton's method which does require the assembly of a Jacobian operator during each iteration. We observed optimal convergence rates when running the solid body rotation testcase. We then modified the scheme to an SUPG enstrophy dissipating scheme which still conserves energy. Using a more efficient time integration scheme that backs off from energy conservation, we demonstrated that this scheme conserves total potential vorticity and has good energy behaviour whilst dissipating enstrophy once fine vorticity filaments form, as is consistent with the enstrophy cascade for 2D geostrophic turbulence.
In future work we will explore the development of vorticity-based schemes for the three-dimensional compressible Euler equations in the context of numerical weather prediction.  : Plots showing convergence of energy error on a finite time interval to zero as ∆t is reduced. The mountain test was run for 15 days at refinement level 6, and then this was used as an initial condition for the convergence test. This shows that the spatial discretisation is exactly preserving energy and any errors are due to the fact that the Picard iteration scheme has not fully converged to the Poisson integrator.