Higher-order compatible finite element schemes for the nonlinear rotating shallow water equations on the sphere

We describe a compatible finite element discretisation for the shallow water equations on the rotating sphere, concentrating on integrating consistent upwind stabilisation into the framework. Although the prognostic variables are velocity and layer depth, the discretisation has a diagnostic potential vorticity that satisfies a stable upwinded advection equation through a Taylor-Galerkin scheme; this provides a mechanism for dissipating enstrophy at the gridscale whilst retaining optimal order consistency. We also use upwind discontinuous Galerkin schemes for the transport of layer depth. These transport schemes are incorporated into a semi-implicit formulation that is facilitated by a hybridisation method for solving the resulting mixed Helmholtz equation. We illustrate our discretisation with some standard rotating sphere test problems.


Introduction
The development of new numerical discretisations based on finite element methods is being driven by the need for more flexibility in mesh geometry.The scalability bottleneck arising from the latitude-longitude grid means that weather and climate model developers are searching for numerical discretisations that are stable and accurate on pseudo-uniform grids without sacrificing properties of conservation, balance and wave propagation that are important for accurate atmosphere modelling on the scales relevant to weather and climate [Staniforth and Thuburn, 2012].There is also ongoing interest in adaptively refined meshes as a way of seamlessly coupling global scale and local scale atmosphere simulations, as well as dynamic adaptivity or even moving meshes; using these meshes requires numerical methods that can remain stable and accurate on multiscale meshes.Further, there is an interest in using higher-order spaces to try to offset the inhomogeneity in the error due to using grids that break rotational symmetry.
Compatible finite element methods are a form of mixed finite element methods (meaning that different finite element spaces are used for different fields) that allow the exact representation of the standard vector calculus identities div-curl=0 and curl-grad=0.This necessitates the use of H(div) finite element spaces for velocity, such as Raviart-Thomas and Brezzi-Douglas-Marini, and discontinuous finite element spaces for pressure (stable pairing of velocity and pressure space relies on the existence of bounded commuting projections from continuous to discrete spaces, as detailed in Boffi et al. [2013], for example).The main reason for choosing compatible finite element spaces is that they have a discrete Helmholtz decomposition of the velocity space; this means that there is a clean separation between divergence-free and rotational velocity fields.Cotter and Shipton [2012] used this decomposition to demonstrate that compatible finite element discretisations for the linear shallow water equations on arbitrary grids satisfy the basic conservation, balance and wave propagation properties listed in Staniforth and Thuburn [2012].In particular, it was shown that the discretisation has a geostrophic balancing pressure for every velocity field in the divergence-free subspace of the H(div) finite element space.A survey of the stability and approximation properties of compatible finite element spaces is provided in Natale et al. [2016], including a proof of the absence of spurious inertial oscillations.
The challenge of building atmosphere models using compatible finite elements is that there is no freedom to select finite element spaces in order to ensure good representation of the nonlinear equations (such as conservation, or accurate advection, for example), because the choice has already been made to satisfy linear requirements.In the case of the rotating shallow water equations, the use of discontinuous finite element spaces for the layer depth field encourages us to use upwind discontinuous Galerkin methods, to solve the continuity equation describing layer depth transport.
The nonlinearity in the momentum/velocity equation is more challenging.In McRae and Cotter [2014], the energy-enstrophy conserving formulation of Arakawa and Lamb [1981] was extended to compatible finite element methods.This extension is closely related to C-grid methods for the shallow water equations on more general meshes in Ringler et al. [2010], Thuburn and Cotter [2012].Following these approaches, the compatible finite element formulation, which has velocity and height as prognostic variables, has a diagnostic potential vorticity that satisfies a conservation equation that is implied by the prognostic dynamics for velocity and height.A finite element exterior calculus structure in this formulation was exposed in Cotter and Thuburn [2014], which also provided an alternative formulation based around low-order finite element methods on dual grids.In Thuburn and Cotter [2015], the close relationship of the dual grid formulation to finite volume methods was exploited to obtain a stable discretisation of the nonlinear shallow water equations on the sphere where the finite element formulation of the wave dynamics was coupled with high-order finite volume methods for the layer depth and prognostic potential vorticity fields.The essential idea is to select a particular stable accurate finite volume scheme for the diagnostic potential vorticity, and to then find the update for the prognostic velocity which implies it.In this paper we address the issue of extending this idea to higher-order finite element spaces, for which there is no analogue of the dual grid spaces.This means that we must return to the formulation of McRae and Cotter [2014], where the potential vorticity is stored in a continuous finite element space.We then seek stable accurate higher-order discretisations of the potential vorticity equation using continuous finite element methods that make it possible to find the corresponding update for prognostic velocity.It turns out that this is indeed possible for advection methods from the SUPG/Taylor-Galerkin family of methods.
Finally, we show how these discretisations can be embedded within a semi-implicit time-integration scheme.We again follow the formulation in Thuburn and Cotter [2015], in which advection terms are obtained from explicit time integration methods applied using the (iterative) velocity at time level n + 1/2.The linear system solved during each nonlinear iteration for the corrections to the field values also requires attention.The standard approach of eliminating velocity to solve a Helmholtz problem for the correction to the layer depth is problematic because the inverse velocity mass matrix is dense.We instead use a hybridised formulation where one solves for the Lagrange multipliers that enforce normal continuity of the velocity field [Boffi et al., 2013, for example].
In section 2 we describe the finite element spaces that satisfy the properties outlined above and provide details of how to construct such spaces on the sphere.Section 3 contains information on the timestepping procedure, including advection schemes for both discontinuous and continuous fields as required.In section 4 we present the results of applying our scheme to some of the standard set of test cases for simulation of the rotating shallow water equations on the sphere as described in Williamson et al. [1992] and Galewsky et al. [2004].Section 5 provides a summary and brief outlook.

Shallow water equations
We begin with the vector invariant form of the nonlinear shallow water equations on a two dimensional surface Ω embedded in three dimensions, where u is the horizontal velocity, D is the layer depth, b is the height of the lower boundary, g is the gravitational acceleration, f is the Coriolis parameter and is the normal to the surface Ω, and where the ∇ and ∇• operators are defined intrinsically on the surface.These equations have the important property that the shallow water potential vorticity (PV) This can be seen by applying ∇ ⊥ • to equation (1).Equation (2) then implies that q is constant along characteristics moving with the flow velocity u, i.e., ∂q ∂t + (u • ∇) q = 0. (5) Numerical discretisations that preserve some aspects of these properties have been demonstrated to be very successful at obtaining long time integrations of the rotating shallow water equations on sphere in the quasigeostrophic flow regime.This is partially because they provide a way to avoid discretising the vector advection term (u • ∇)u in the velocity equation directly; instead one can choose a suitable stable, accurate and conservative scalar advection scheme for the potential vorticity (treated as a diagnostic variable with u and D being prognostic variables) and use it to diagnose a form of the (ζ + f )u ⊥ term in equation ( 1) that leads to stable advection of u.These ideas were introduced in the compatible finite element context in McRae and Cotter [2014] in order to obtain energy-enstrophy conserving discretisations; here we concentrate on stable, accurate and conservative advection of q, improving on the low-order APVM stabilisation suggested there.We also replace the centred discretisation of equation ( 2) with stable and accurate Discontinuous Galerkin (DG) advection schemes for D, and show how these can be incorporated into the PV conserving formulation.

Finite element spaces
In this section we shall summarise the properties we require from our finite element spaces and the operators between them.We start with the space H(div) of square integrable velocity fields, whose divergence in also square integrable.The condition that the discrete velocity belongs to the finite element subspace V 1 ⊂ H(div) means that the velocity must have continuous normal components across element edges.Having chosen V 1 , we select a finite element space V 2 ⊂ L 2 such that This necessarily requires that V 2 is a discontinuous space.We also define a space V 0 ⊂ H 1 consisting of continuous fields γ such that k × ∇ψ ∈ V 1 , where the curl k × ∇, henceforth written as ∇ ⊥ , maps from V 0 onto the kernel of The proof that the mixed finite element discretisation of the linear shallow water equations has steady geostrophic modes relies on the existence of a discrete Helmholtz decomposition for the velocity field [Cotter and Shipton, 2012].As described in Arnold et al. [2006], this decomposition exists if the following diagram commutes with bounded projections π 1 , π 2 , π 3 , that is, the result of applying an operator to the continuous field and projecting into the discrete space is the same as the result of first projecting the field into the discrete space and then applying the operator.Cotter and Shipton [2012] reviewed several sets of finite element spaces that satisfy these requirements, together with further requirements on the degree-of-freedom (DOF) ratios between V 1 and V 2 that are necessary to exclude the possibility of spurious mode branches in the dispersion relation for the linear shallow water equations.In this paper with shall present results from choosing V 1 to be the first order Brezzi-Douglas-Marini (BDM 2 ) function space on triangles which requires that V 0 is P 3 , the space of piecewise continuous cubic functions.V 2 is the space P DG 1 of piecewise linear functions that can be discontinuous at the element boundaries.

Constructing the finite element spaces on the sphere
In order to implement the finite element method, we need to expand our variables in terms of suitable basis functions and compute integrals of combinations of these basis functions over elements.This is done by defining a reference element on which these integrals can be calculated and mappings from each reference element to the physical element, which we describe in this section.A more detailed and general exposition, plus description of how these mappings are implemented in the FEniCS project, is provided in Rognes et al. [2013].In the rest of this section, hatted quantities refer to those defined on the reference element.We start by defining finite element spaces V i ( ê), i = 0, 1, 2, on the reference element ê.These are constructed from polynomials in the usual manner for the chosen spaces as described in Boffi et al. [2013], for example.Then, for each element e, we construct a polynomial mapping g e : ê → R 3 , such that where Ω is the computational domain, which is a piecewise polynomial approximation to the sphere domain Ω.
Then, we use g e to relate functions in our finite element spaces on Ω, restricted to each element e in Ω, to functions in V i ( ê).
For V 0 , V 2 , these functions are simply related by function composition, i.e., For V 1 , we use the Piola transformation where u| e is the restriction of u to element e, and is the Jacobian of g, to ensure that u is tangent to Ω.The Piola transformation has the crucial property that When the discrete space is made up of flat elements, as in Cotter and Shipton [2012] and McRae and Cotter [2014], det J e is constant.This guarantees that if ∇ However, for meshes made up of general quadrilaterals (i.e.cubed sphere meshes such as those in Putman and Lin [2007]) or high-order, curved triangles, the mapping g is no longer affine and the determinant of its Jacobian is no longer constant.This means that in general ∇ • u / ∈ V 1 .This clearly violates the commutative diagram property (7).There are 3 ways to remedy the situation.

Modify the mapping for
This option is the choice that is most consistent with the FEEC methodology.
2. Modify the mapping for V 1 so that the factor of det J e in ∇ • u| e • g e is replaced by the element average of det J e .This approach was described in Boffi and Gastaldi [2009] for the case of lowest order Raviart-Thomas elements; the extension to other H(div) elements is straightforward but the construction is quite complicated.

Replace the divergence operator
By defining a new divergence operator in this way we immediately recover the required commutation property, since π 2 appears in the commutative diagram (7) and is a projection.This is a generalisation of the "rehabilitation" technique described for lowest order Raviart-Thomas elements applied to mixed elliptic problems in Bochev and Ridzal [2008].In fact, as Bochev and Ridzal noticed, introduction of the ∇• operator does not require any changes to a code implementation for mixed elliptic operators since the ∇• operator only appears in an inner product with a test function from V 2 , and hence can be safely replaced by ∇• there.
We shall see that this continues to be the case in our nonlinear shallow water formulation.
In general, Option 1 is the most mathematically elegant choice.However, since our software implementation is based upon the Firedrake finite element library [Rathgeber et al., 2016, Luporini et al., 2016] which is already used to develop DG schemes that use the transformations in (9), Option 3 was much less pervasive through the code base.We will investigate the differences between Options 1 and 3 in future work.Arnold et al. [2014] showed that there is a potential problem with loss of consistency with compatible finite elements on quadrilateral and cubic curvilinear cells; it is likely that this problem also can be exhibited on triangular cells.In particular, for general constructions there may be loss of consistency for the V 1 spaces used here (the consistency problem for V 2 is avoided through rehabilitation).However, Holst and Stern [2012] demonstrated that consistent approximation can still be obtained if the computational domain can be obtained via a transformation from a mesh of affine elements onto the higher-order nodal interpolation of a C ∞ manifold (such as the sphere).This aspect is discussed further in the context of geophysical fluid dynamics in Natale et al. [2016].

Mixed finite element formulation
The mixed finite element discretisation of equations ( 1)-( 2) is formed by restricting u ∈ V 1 , D ∈ V 2 , multiplying the equations by appropriate test functions, w ∈ V 1 and φ ∈ V 2 , and integrating over the domain, giving where we have introduced the mass flux F ∈ V 1 (an approximation to uD), and the vorticity flux Q (an approximation to uDq).F and Q are defined by appropriate choice of advection schemes for D and the diagnostic potential vorticity q, which we shall define later.
We have integrated the gradient term in equation ( 15) by parts to avoid taking the gradient of the layer depth which is undefined since D ∈ V 2 is discontinuous.A similar problem is posed for the definition of the potential vorticity q since ζ = ∇ ⊥ • u, which appears in the definition of q, is similarly undefined.In order to fix this we integrate the ∇ ⊥ term by parts (with no surface term since we are solving our equations on the sphere), defining our discrete potential vorticity q ∈ V 0 as Provided D remains positive (this can be enforced using a slope limiter, although in the test cases used here the mean depth is sufficiently high that a slope limiter is not required), then this equation can be solved for q ∈ V 0 from known D and u.This provides a one-to-one mapping between q and the weak curl of u, which leads to a oneto-one mapping between q and the divergence-free part of u via the discrete Helmholtz decomposition (without needing to solve a global Poisson problem).Hence, control of q in the L 2 norm from a chosen stable advection scheme provides strong control over the divergence-free component of u without compromising the divergent part.This is one reason for the success of this type of potential vorticity conserving discretisation for the rotating shallow water equations.We differentiate this equation in time, and substitute for u t using equation ( 15) with w = −∇ ⊥ γ.Since ∇•∇ ⊥ ≡ 0 and we assume f t = 0, this gives which is the Galerkin projection of the PV conservation law into V 0 .If we select γ = 1, then we obtain conservation of total PV, On the other hand, if we are on the sphere, then this quantity is zero as can be computed directly from (17); this topological result stems from the fact that in this case w = −∇γ = 0.If we choose our advection scheme so that Q = qF if q is a constant, then we may integrate by parts without introducing error (since q ∈ H 1 and F ∈ H(div)), and we obtain and hence For flat elements, or if we had chosen Option 1, the right hand side is zero since D and ∇ • F are in the same finite element space and therefore equation ( 16) holds pointwise, i.e., Hence, we conclude that if q is constant, then q t = 0, and so q remains constant.As well as being a statement of first-order consistency for the advection scheme for q, this is also an important property of equation ( 5).
The formulation requires some adaptation for the finite element spaces used in this paper to obtain this result, since ∇ • F is not guaranteed to be in the same space as D, so equation ( 22) does not hold pointwise.To recover it, some further "rehabiliation" is required; we amend equation ( 17) by defining D ∈ V 0 such that Dt where Comparing the weak form of this equation with ( 16) we see that which is consistent with the definition hence we can solve for D (since τ is a positive quantity and just alters the metric).D is discontinuous and hence this equation can be solved separately in each element.Using this in equation ( 17) gives Differentiating, rearranging and assuming q is constant as before we obtain as required.
In fact, the finite element formulation allows us to go beyond the first order consistency result.For example, if we choose Q = F q, we have enough continuity to integrate by parts, and we obtain The left-hand side vanishes if q is an exact solution of the equation where D and F are the discrete mass and mass flux, indicating that the discretisation is consistent at the order of approximation of the finite element space V 0 .Unfortunately, this is not a good choice for Q since the implied discrete advection operator is not stable.An alternative is to choose where α is a stabilisation parameter.Substitution, rearrangement and integration by parts on the F q term leads to which is a streamline-upwind Petrov-Galerkin (SUPG) spatial discretisation of the PV conservation equation.This discretisation is stable, and also consistent, i.e. the equation vanishes when the exact solution to (30) is substituted.
In this paper we will use a Taylor-Galerkin discretisation for the potential vorticity equation; this discretisation achieves the same aims as the SUPG discretisation described above, but arises more naturally in the discrete time setting and hence we shall postpone our discussion of it until we have described the time-discrete formulation of the full shallow water system in the next section.

Time discretisation
We shall build a semi-implicit time discrete formulation.First we write We can now write equations ( 15) and ( 16) as where and where the time-averaged mass and vorticity fluxes F and Q are yet to be defined.The idea is that we choose F to be a time-independent function such that D n+1 is obtained from D n via a high-order stable time discretisation over one timestep for the equation i.e. with the advecting velocity frozen to the value of u * .Similarly, Q is chosen so that q n+1 is related to q n via a high-order stable time discretisation over one timestep for the equation i.e. with the advecting mass flux frozen to the time-averaged flux F .This means that for θ = 1/2 we obtain a scheme that is overall second-order in time, but that uses higher-order advection schemes for D and q.The rationale is that for near-linear waves, we would like the propagation to be as conservative as possible, so the semi-linear formulation should be based around a time-centred scheme.However, we would also like to obtain good quality solutions over long integrations when close to geostrophic balance, in which case the important quantity is q, and it is important that we transport D consistently with q to stay close to the balanced state.In that regime, it is thought to be important to use a high odd-order time integration scheme, since for odd-order schemes the error is dominated by diffusion rather than dispersion (the latter leads to oscillations near to near-discontinuous data).It is also thought that the use of a time-averaged velocity to transport q and D helps to preserve geostrophic balance.This was also the rationale behind choosing the 3rd order Forward-in-Time advection schemes used in Thuburn and Cotter [2015].
An implicit formulation such as the one above requires Newton or Picard to iterate to convergence.In practice, we perform a small fixed number of Picard iterations per timestep, since our aim is to obtain a stable timestepping method with accurate q and D transport, rather than iterating to convergence to obtain an exact implementation of the fully implicit scheme described above.These are "quasi" Newton iterations since we replace the Jacobian obtained from linearisation around the current state with the Jacobian linearised around the system at a state of rest.This means that it is possible to reduce the system as described in the next section, and that the same solver context can be reused during each Picard iteration and timestep.
Hence, each Picard cycle consists of the following steps.
2. Use the current values of ∆u and ∆D to compute corresponding values for u * , D * .
3. Use the chosen mass advection scheme with velocity u * to update from D n to D n+1 , and find F such that We will explain how to construct F in Section 3.2.The mass residual R D : V 2 → R is defined as 4. Diagnose the PV q n at time t n .
5. Use the chosen PV advection scheme with mass flux F to update from q n to q n+1 and compute the corresponding PV flux Q such that 7. The increments ∆u → ∆u + δu, ∆D → ∆D + δD (44) are then obtained by solving the coupled system where H 0 is the mean layer depth at rest.
8. If more iterations are required, apply these updates and return to 2.
In the following sections we describe the construction of F , Q and the solution of the coupled system in detail.

Solving the coupled linear system
In our formulation, we use hybridisation to solve equations (( 45)-( 46)).Hybridisation is a technique for efficiently solving mixed finite element problems that has been used since the 1960s; in the 1980s it was also discovered that the hybridised formulation could also be used to obtain more accurate approximations of the solution (see Boffi et al. [2013] for a general survey).Obtaining a hybridised formulation requires two steps.First, we introduce a finite element space Ṽ1 by relaxing the normal continuity constraints within V 1 .In other words, functions in Ṽ1 have the same local polynomial representation as functions in V 1 , but there are no requirements of continuity between edges.In particular, we note that V 1 ⊂ Ṽ1 .Second, we introduce a trace space Tr(V 1 ), defined on the element facet set Γ, such that functions λ ∈ Tr(V 1 ) are scalar functions which when restricted to a single element facet f , are from the same polynomial space as u • n restricted to that facet.Having relaxed the continuity requirements for δu ∈ Ṽ1 , we enforce them again by adding another equation, where we use the usual "jump" notation having arbitrarily labelled each side of each facet with + and −, so that n + points from the + side to the − side and vice versa.To avoid an over-determined system, we introduce Lagrange multipliers λ ∈ Tr(V 1 ) and rewrite equations (( 45)-( 46)) as together with equation ( 47).Note that the residual R u must now be evaluated with w ∈ Ṽ1 .All of the inter-element coupling in equations (( 49)-( 50)) takes place in the λ term.This means that if λ is known, then it is possible to obtain u and D independently in each element.To enable elimination of u and D, we define a lifting operator ), which gives the solution u for a given λ in the case when R D [φ] and R u [w ] are replaced by zero.We also define u 0 which is the solution obtained when λ is zero, but R D [φ] and R u [w ] are present.Then, the general solution of this equation given particular R D and R u is Substituting into equation ( 47) gives This reduced equation can be solved for λ before reconstructing u and D by solving equations (( 49)-( 50)) independently in each element.Since the value of Lλ in each element only depends on the values of λ on the facets of that element, equation ( 52) only couples together values of λ on facets that share an element.This means that the matrix-vector form of this equation is very sparse.In fact, the matrix can be assembled by visiting each element separately, performing inversion on element block systems.Further, this equation has the same spectral properties as the Helmholtz operator, and hence can be solved with Krylov methods and standard preconditioners such as SOR, algebraic multigrid; a geometric multigrid for general higher-order hybridised mixed finite element elliptic problems was provided in Gopalakrishnan and Tan [2009].Note that the Coriolis term can be included in the linear system in this approach without altering the sparsity of the reduced system.One important aspect is that if the solver for this system is not iterated to convergence, the resulting velocity field will not be exactly divergence-free.We address this in our implementation by simply projecting the velocity back into V 1 after reconstruction, which appears not to cause any problems with stability.A more sophisticated approach would use the hybridised solver as a preconditioner for the (∆u, ∆p) system on (V 1 , V 2 ); we will investigate this in further work.

Advection scheme for layer depth D
We now need to solve the mass continuity equation ( 16) for the update to D. As D ∈ V 2 and is therefore discontinuous, we can use standard upwind discontinuous Galerkin methods to obtain an approximation to D t , for each element e, where V 1 (e) is the space V 1 restricted to the element e, and D u is the value of D on the upwind side of the element boundary ∂e.We then use the standard 3-stage Strong Stability Preserving Runge-Kutta scheme [Gottlieb et al., 2001], Later we shall discuss a consistency property of the potential vorticity conserving discretisation of the velocity equation; this property requires that we find a time-integrated mass flux F ∈ V 1 that satisfies This is done by finding F given D t such that and then substituting into equations (( 54)-( 56)) to construct F .To find a flux F for each Runge-Kutta stage, we solve where V * 1 (e) is the right size to close the system and ∇φ ∈ V * 1 (e) for all φ ∈ V 2 (e).The left-hand sides of equations (( 59)-( 60)) are the same as the left-hand sides of the definition of the commuting operator π : H(div) → V 1 that features in stability proofs for mixed finite element methods (see Boffi et al. [2013], for example).This means that the above construction is well-posed.To check that equation ( 58) is satisfied, we integrate (57) by parts.Then, substituting the above relations ( 59)-( 60 We note, although we did not use it in this paper, that the use of a slope limiter can also be incorporated into the mass flux computation.If a slope limiter is used after a Runge-Kutta stage, then D is replaced by D , with then integration over a single element immediately tells us that this can be satisfied by i.e.F • n = 0 on the boundary ∂e.We then solve a local mixed problem for (F , p) given by subject to the above zero Dirichlet boundary conditions for F .

Potential vorticity conserving formulation
The advection scheme described in this section follows the following design strategy: find an advection scheme for q that is compatible with equation ( 35), and select the corresponding Q for insertion into that equation to compute R u [w ].Selecting w = −∇ ⊥ γ in equation ( 35), evaluating equation ( 17) at time levels n and n + 1 and substituting gives after noting that w is divergence-free.McRae and Cotter [2014] used Q = F q n+1/2 to obtain an implicit midpoint rule time discretisation for an energy-enstrophy conserving formulation.In this paper, we aim to use higher-order stabilised advection schemes in order to obtain accurate representation of potential vorticity transport.We note that Streamline Upwind Petrov Galerkin methods [Brooks and Hughes, 1982] and Taylor-Galerkin methods [Donea, 1984] all result in time-discrete formulations equivalent to equation (65).
It is desirable to use a higher order timestepping scheme for q, to obtain accurate transport of potential vorticity using the balanced flow.In particular, odd-ordered schemes are attractive since the leading order error is diffusive rather than dispersive.Hence, in this paper we make use of the two-step third-order unconditionally-stable Taylor-Galerkin scheme of Safjan and Oden [1993].
Taylor-Galerkin schemes are built by transforming time derivatives into space derivatives using the advection equation.
The general form of a multistage Taylor-Galerkin method is where η is a stabilisation parameter, the subscript i is the stage index, and the {µ} ij and {ν} ij are coefficients defined in Safjan and Oden [1993].After computing these stage variables, the value of Z k is copied into Z n+1 .
In our discretisation we use a formulation where Z = qD.Then we have recalling that F = D ū is considered to be time-independent over the advection step as part of the semi-implicit discretisation.Combination with Equation (66) leads to Note that this equation involves solving a Helmholtz-type equation with derivatives in the streamline direction for each stage variable.This equation is symmetric positive definite and well-conditioned for O(1) Courant numbers; the conjugate gradient method converges quickly with simple SOR preconditioning.
In this paper we took the following coefficients Safjan and Oden [1993] showed this scheme to be third order and unconditionally stable for η > 0.473.
Having solved for q n+1 , we take i = 2, and notice that equation (69) takes the form of equation ( 65), and hence Q can be extracted for insertion into equation (35).For the case of curved elements, D must be replaced by D/τ throughout.

Numerical Results
In this section we show numerical results from three standard test cases on icosahedral grids, using the spaces (P 3 , BDM 2 , DG 1 ), and a piecewise cubic approximation to the surface of the sphere.The code was implemented using the Firedrake library, with the GAMG algebraic multigrid preconditioner being used for the hybridised system in the implicit solver.
The first two test cases are numbers 2 (solid body rotation) and 5 (flow over a mountain) from Williamson et al. [1992] and the third is the barotropically unstable jet from Galewsky et al. [2004].Table 1 contains information on the properties of the 4 grids that we use for the Williamson convergence tests, along with the timesteps used.The timestep was chosen to give a constant Courant number across the different resolutions; 0.2 in the solid body rotation test and 0.06 in the flow over a mountain test.Figure 1 shows the lowest resolution icosahedral grid that we use.
Where shown, normalised errors in a field q are computed as in Williamson et al. [1992] as where q T is the true solution, specified either analytically (as in the solid body rotation test) or from a high resolution reference solution (as in the flow over a mountain test).1: Grid properties for the 4 grids used in the convergence tests, including the number of degrees of freedom (DOFs) for velocity u, depth D and potential vorticity q, along with the timestep used for the solid body rotation test (2) and flow over a mountain test (5).
Figure 1: Icosahedral sphere grid, viewed looking down on the North pole, corresponding to the lowest resolution described in table 1.

Solid body rotation (Williamson, test case 2)
This test case is initialised with depth and velocity fields that are in geostrophic balance: where R = 6.37122 × 10 6 m is the radius of the Earth, Ω = 7.292 × 10 −5 s −1 is the rotation rate of the Earth, (x, y, z) are the 3D Cartesian coordinates, g = 9.80616ms −2 is the gravitational acceleration, D 0 = 2.94 × 10 4 /g ≈ 2998m and u 0 = 2πR/(12 days) ≈ 38.6ms −1 .Since there is no topography or other forcing, the flow should remain in this steady, balanced state.As the flow is steady we have an analytic solution which allows us to compute errors and hence an order of convergence for our method.We present these results in table 2 and Figure 2 reveals that our method is converging at second order.The depth error at day 15, shown in Figures 3 and 4, is large scale and shows some evidence of grid imprinting.

Flow over a mountain (Williamson, test case 5)
This test case uses the same zonal flow initial conditions ( 73)-(74) as the previous test but with D 0 = 5960m and u 0 = 20m.An isolated, conical mountain, given by cells    ] is placed with its centre at latitude φ = π/6 and longitude λ = −π/2.As the zonal flow interacts with the mountain, it produces waves that travel around the globe.As there is no analytical solution for this problem, the model output at 15 days is compared to a high resolution reference solution generated using a semi-Lagrangian shallow water code provided by John Thuburn.We plot the depth errors at day 15 in Figure 6.We can see that the error is small scale and is not dominated by errors due to grid imprinting.
Up until this point the flow is only weakly nonlinear (see Figure 7) so, as in Thuburn et al. [2013], we continued the highest resolution simulation until day 50.By this time, fine scale structure has been generated as the flow becomes more nonlinear; the PV field develops sharp gradients and filaments that stretch out and roll up.These features can be seen in the potential vorticity field at 50 days, shown in Figure 8.

Barotropically unstable jet (Galewsky)
The details of this test case are specified in Galewsky et al. [2004].The initial condition consists of a strong midlatitude jet, with an added perturbation, which is barotropically unstable and evolves to produce vortices and small scale structure.It has become a particularly useful test for models on grids that are not aligned with latitude/longitude because these grids can induce early development of the instability and even lead to the final solution (at day 6) having the incorrect wavenumber.Again, there is no analytic solution so results are compared to those in the literature (see for example Thuburn et al. [2013] and Weller [2013]).An important feature to reproduce is the relatively straight path of the jet across approximately quarter of the globe -this is not seen in models where grid imprinting has resulted in the generation of instability along the length of the jet [Weller, 2013].
Figure 9 shows the vorticity field, computed on the highest resolution grid (see table 1) with timestep of 120s, between 10 • and 80 • latitude after 6 days.The maximum Courant number reached during this simulation is 0.28.We note that the jet has the correct wavenumber and the expected quiescent section.Lower resolution simulations, as others have found [Thuburn et al., 2013], did not have these requisite features.

Summary and Outlook
We have built upon the work of Cotter and Shipton [2012] and McRae and Cotter [2014] to produce a semiimplicit compatible finite element model for the nonlinear rotating shallow water equations on the sphere.The important features are that we introduce higher-order upwind advection schemes that maintain PV conservation, and consistency of PV advection with the mass conservation law.By applying the model to standard test cases we have demonstrated that this model has the expected second order convergence rate and can produce the required features of both large scale balanced flows and unstable turbulent flows.The developments of this paper inform our ongoing development of a three dimensional dynamical core on the sphere.

Figure 2 :Figure 3 :Figure 4 :
Figure 2: Solid body rotation test case: normalised depth, D, and velocity, u, errors at day 5 versus average mesh size ∆x.

Figure 5 :
Figure 5: Flow over a mountain test case: normalised depth errors at day 15 versus average mesh size ∆x.

Figure 6 :
Figure 6: Flow over a mountain test case: depth errors (in metres) at day 15.The range is min:−1.87m,max: 1.60m.

Figure 9 :
Figure 9: Vorticity at day 6 for the Galewsky jet test.

Table 2 :
Normalised depth, D, and velocity, u, errors at day 5 for the solid body rotation test case.