A mimetic, semi-implicit, forward-in-time, finite volume shallow water model: comparison of hexagonal–icosahedral and cubed-sphere grids

. A new algorithm is presented for the solution of the shallow water equations on quasi-uniform spherical grids. It combines a mimetic ﬁnite volume spatial discretization with a Crank–Nicolson time discretization of fast waves and an accurate and conservative forward-in-time advection scheme for mass and potential vorticity (PV). The algorithm is implemented and tested on two families of grids: hexagonal–icosahedral Voronoi grids, and modiﬁed equiangular cubed-sphere grids. Results of a variety of tests are presented, including convergence of the discrete scalar Laplacian and Coriolis operators, advection, solid body rotation, ﬂow over an isolated mountain, and a barotropically unstable jet. The results conﬁrm a number of desirable properties for which the scheme was designed: exact mass conservation, very good available energy and potential enstrophy conservation, consistent mass, PV and tracer transport, and good preservation of balance including vanishing ∇ × ∇ , steady geostrophic modes, and accurate PV advection. The scheme is stable for large wave Courant numbers and advective Courant numbers up to about 1.


Introduction
In order to achieve the parallel scalability needed to exploit future generations of supercomputers, weather and climate prediction models will need to use quasi-uniform spherical grids. A significant challenge is to develop schemes that can achieve comparable accuracy to current state-ofthe-art longitude-latitude grid models, without excessive cost (Staniforth and Thuburn, 2012). With this motivation, Thuburn and Cotter (2012) recently presented a framework for the construction of finite volume schemes for the solution of the rotating shallow water equations. In this framework, key physical properties related to conservation, balance, and potential vorticity (PV) dynamics are obtained by ensuring that the numerical scheme mimics certain mathematical properties of the continuous governing equations. The framework makes use of a primal polygonal grid with a C-grid placement of variables and the corresponding dual grid, along with a set of linear operators with certain symmetry properties for mapping between the two. It extends the work of Thuburn et al. (2009) and Ringler et al. (2010) to the case in which dual grid edges are not necessarily orthogonal to primal grid edges, thus making it applicable to quasiuniform cubed-sphere grids, for example, as well as Voronoi grids. However, Thuburn and Cotter (2012) did not provide any specific examples of the required linear operators for non-orthogonal grids. Here we present a set of operators suitable for a particular class of cubed-sphere grid, and compare the resulting model with one using the Ringler et al. (2010) operators on a hexagonal-icosahedral Voronoi grid.
Published by Copernicus Publications on behalf of the European Geosciences Union.

Governing equations
The vector invariant form of the shallow water equations is used: Here, φ is the geopotential given by the fluid depth times the gravitational acceleration, φ T = φ + φ orog is the total geopotential at the fluid's upper surface including the contribution from orography, v is the velocity, and k = |v| 2 /2. A superscript ⊥ indicates that the vector in question is rotated through π/2 in the positive (anticlockwise) direction: v ⊥ = k ×v where k is the unit vertical vector. Finally, vφq is the flux that appears in the conservation law for PV (derived from Eqs. 1 and 2): where q = ζ /φ is the PV, ζ = f + ξ is the absolute vorticity, and ξ = k · ∇ × v is the relative vorticity.

Grids
The shallow water code we have developed is formulated for an arbitrary unstructured grid. However, a considerable amount of grid-related information, including the operators listed in Table 3 below and the multigrid restriction and prolongation operators (Sect. 4), needs to be generated and pre-tabulated for any given grid. So far, grid generators and operators for two families of grids have been developed: hexagonal-icosahedral Voronoi grids and cubedsphere grids. The hexagonal-icosahedral grid is essentially that proposed by Heikes and Randall (1995a, b) (but without the twist). The primal grid comprises hexagonal and pentagonal Voronoi cells, while its dual is the corresponding Delaunay triangulation. The grid generation code iteratively adjusts the primal grid cell centres so as to minimize a cost function J = α HR J HR + α C J C , where J HR is the cost function used by Heikes and Randall (1995b), which penalizes failure of primal and dual edges to cross at their midpoints, and J C penalizes departures of the primal cell "centres" (i.e. dual grid vertices) from primal cell centroids. Setting (α HR , α C ) = (1, 0) gives the Heikes-Randall grid; setting (α HR , α C ) = (0, 1) gives a centroidal Voronoi grid as described by Du et al. (1999) and used by Ringler et al. (2010), and Skamarock et al. (2012). All of the results shown below use 40 iterations of a simple cell-by-cell minimization algorithm to minimize J with (α HR , α C ) = (1, 0). (A few experiments using (α HR , α C ) = (0, 1) suggest that there is only weak sensitivity of results to this choice.) The default orientation of the grid, used in all tests below, has a pentagon at each pole. For the cubed-sphere grid, the vertices are first positioned as on the equiangular cubed-sphere (e.g. Ronchi et al., 1996), then the dual vertices are positioned at the barycentres of the surrounding primal vertices, and finally the primal vertices are relocated to the barycentres of the surrounding dual vertices. The last step is needed in order to use the H operator described in Sect. 2 below and the Appendix. It has the effect of smoothing the grid somewhat across the cube edges, which is probably beneficial. (Iterating the last two steps leads to further smoothing and resolution clustering at the cube vertices; after many iterations the grid resembles the conformal cube (Rancic et al., 1996). For this reason the last two steps are not iterated for the results shown below.) The default orientation of the grid, used in all tests below, has the cube corners at latitudes ±π/4. Figure 1 shows coarse resolution versions of the two grids. Some characteristics of the grids at different resolutions are given in Table 1. The resolutions on the cubed-sphere have been chosen to give approximately the same number of degrees of freedom as one of the hexagonal-icosahedral grids, allowing a fair comparison between the two grid types.

Summary of the framework, and specific operators
The framework of Thuburn and Cotter (2012) is expressed in terms of variables integrated over relevant grid elements (cells or edges) or sampled at vertices. See Table 2. This has the advantage that many geometrical factors such as lengths and areas are absorbed into the field variables, helping to make both the mathematical formulation and the computer code simpler and clearer. We will use subscripts on variables (e.g. φ i or V e ) to refer to specific grid elements (geopotential at dual vertex i or circulation integrated along dual edge e), and omit subscripts (e.g. φ or V ) to refer to the entire vector of these variables at all relevant grid elements; this allows the use of a compact matrix-vector notation. In terms of these integrated variables, the vector calculus operations gradient, divergence and curl then take particularly simple forms: -Gradient integrated along a dual edge -Gradient integrated along a primal edge -Divergence integrated over a primal cell -Curl integrated over a dual cell where the sparse matrices D 1 , D 1 , D 2 , and D 2 are determined by the grid topology and their non-zero entries are all equal to +1 or −1. The meanings of these and the other mimetic operators are briefly summarized in Table 3 and Fig. 2. See Thuburn and Cotter (2012) for detailed mathematical definitions and discussion. Discretizing Eqs.
(1) and (2) in space giveṡ Here the prognostic variables are i , the integral of φ over primal cell i, and V e , the circulation along dual edge e. A time derivative is indicated by( ). F e is the mass flux across primal edge e; it is computed from the and U fields using the primal grid advection scheme -see Sect. 5. Q ⊥ e is the PV flux across dual edge e. It is computed from the C ⊥ and J −1 q fields using the dual grid advection scheme, where C ⊥ e = (WF ) e is the mass flux across dual edge e, and Finally, K i is the kinetic energy per unit mass integrated over primal cell i. Motivated by the fact that the discrete approximation to the global integral of kinetic energy is given Absolute vorticity integrated over dual cell v However, although this is first-order accurate on the hexagonal Voronoi grid, it is not on the cubed-sphere grid, and convergence of the maximum φ error was observed to stall for the solid body rotation test case (Sect. 6.5). Therefore, an alternative scheme is used for all results presented below. A constant vector velocity u i is constructed for each cell i that gives a least squares best fit to the V e at the edges of that cell. The kinetic energy in cell i is then approximated by |u i | 2 /2 times the area of cell i. Provided the operators satisfy certain symmetry conditions, this framework ensures a number of desirable properties for the scheme.
-The placement of degrees of freedom is that of a polygonal C-grid, which helps to ensure an accurate representation of the geostrophic adjustment process (Arakawa and Lamb, 1977).
-Equation (8) is manifestly in conservative form, ensuring conservation of mass.
-It is a topological identity that D 2 ≡ −D T 1 ; this, together with the condition that the matrices I and H should be symmetric, provides a discrete analogue of the property that ∇ is minus the adjoint of ∇ · ( ), ensuring that the geopotential gradient term is energy conserving for the linearized shallow water equations.
-The condition that W should be antisymmetric ensures that the Coriolis term is energy conserving for the linearized shallow water equations.
-Another topological identity is that D 2 D 1 ≡ 0, which leads to a discrete analogue of the identity ∇ × ∇ ≡ 0; this ensures that the geopotential gradient term cannot act as a source of vorticity or potential vorticity.
-The R operator must be local and conservative, so that the global integral of any variable over the primal grid is equal to the global integral of (v) = R over the dual grid. This property enables the construction of a unique W that satisfies −D 2 W = RD 2 (as well as antisymmetry), as described by Thuburn et al. (2009). Provided W and R are related in this way, any nondivergent velocity field can be balanced by some field, and will produce no vorticity via the Coriolis term; in other words, the scheme can support geostrophically balanced flows.
-A corollory is that the linearized potential vorticity   Lin and Rood (1997), in the following sense: we can define a PV flux Q e by discretizing the PV conservation law where Q = vφq; then, by substituting this Q e in Eq. (9), the PV can be made to evolve exactly as if we were to integrate Eq. (11) directly, even though we actually integrate Eqs. (8) and (9). This means we can use any desired advection scheme to compute the PV flux, giving a high degree of control over the PV evolution.
Here we have emphasized conservation of energy only for the linearized equations. If desired, a scheme can be constructed that conserves energy for the full nonlinear system (Thuburn and Cotter, 2012). However, this requires a somewhat artificial construction of the PV flux. The philosophy adopted here is that we prefer to have the PV evolve according to a chosen advection scheme, maximizing our control of the PV. Upwind advection schemes of the sort described below are inherently damping on small scales, leading to dissipation of potential enstrophy and, to some extent, energy. Kent et al. (2012) and Thuburn et al. (2014) discuss the extent to which this approach can provide an implicit subgrid model capturing cascades of potential enstrophy and energy for two-dimensional vortical flow. Apart from this inherent dissipation due to the advection scheme, no other dissipation mechanism (either explicit or inherent in the numerics) is needed in the model to maintain stability or to control dispersion errors or other numerical sources of noise.
All of the above properties are quite general. It remains to define specific instances of the I, H, J, and R operators for the two grids used here.
For the hexagonal-icosahedral grid, the I, H and J operators implemented are all diagonal, i.e. their stencil is a single cell or edge. This is just the translation into the present framework of the operators used by Thuburn et al. (2009) and Ringler et al. (2010).
where A i is the area of primal cell i. This I operator will be first-order accurate provided the dual vertex i lies within primal cell i. (It would be second-order accurate on a centroidal Voronoi grid.) H e e = l e /d e , e = e, 0, where l e is the length of primal edge e and d e is the length of dual edge e. On an orthogonal grid this diagonal H will be first-order accurate provided the primal and dual edges do intersect each other, and will be second-order accurate if the grid is constructed so that the intersection point approaches the midpoint of both the primal edge and the dual edge as resolution is increased (Heikes and Randall, 1995b). Finally, where v is the area of dual cell v. This J operator will be first-order accurate provided the primal vertex v lies within dual cell v. The R operator is defined as in Ringler et al. (2010), in which the mapping weights are proportional to the overlap area between primal and dual cells. This operator is first-order accurate. The W operator is constructed from R following Thuburn et al. (2009).
For the cubed-sphere grid, I and J are again the diagonal operators defined by Eqs. (12) and (14). However, because the primal and dual edges are not orthogonal to each other, a diagonal H operator would be inconsistent (i.e. not even first-order accurate). Instead we use the following: www.geosci-model-dev.net/7/909/2014/ Geosci. Model Dev., 7, 909-929, 2014 The stencil S comprises the five edges nearest to edge e, including edge e itself. d e is a vector of magnitude d e in the direction of dual edge e. Subscript c refers to the corner formed by the edges e and e . s c = 4 when dual edges e and e are edges of the same quadrilateral, and s c = 6 when dual edges e and e are edges of the same triangle. This operator exactly converts V to U for a constant velocity field on a plane, and is therefore first-order accurate, provided the primal grid vertices are located at the barycentres of the surrounding dual vertices. See the Appendix for more details.
The construction works when the dual cells are either quadrilaterals or triangles. Again, the R operator is defined as in Ringler et al. (2010) and W is constructed from R following Thuburn et al. (2009). Note that the above defines the I, H, J, and R operators everywhere on the two grids; no special handling is used near pentagons on the hexagonal-icosahedral grid or near cube edges or corners on the cubed-sphere grid.

Time integration scheme
The time integration scheme is motivated by the observation that a Crank-Nicolson scheme gives a robustly stable treatment of fast waves, with sufficient accuracy to capture geostrophic adjustment, while a forward-in-time finite volume advection scheme can provide exact mass conservation and stability up to Courant number ≈ 1 while accurately capturing Lagrangian conservation. Therefore, consider the following scheme, obtained by integrating Eqs. (8) and (9) over a time interval t: Here, indicates a (possibly off-centred) trapezoidal approximation to the Eulerian time integral for any variable ψ, F is the time integral of the mass flux, given by the primal grid advection scheme using the time-integrated fluxes U t , and Q ⊥ is the time integral of the potential vorticity flux, given by the dual grid advection scheme using the mass fluxes C ⊥ = W F . All of the results presented below, except in Sect. 6.8, use a centred approximation to the time integral: α = β = 0.5.

Iterative solution and Helmholtz problem
The K, F and Q ⊥ terms in Eqs. (16) and (17) depend nonlinearly on the values of the predicted variables at step n + 1. Hence we must solve a coupled nonlinear system of equations at each step. An incremental iterative approach is used. Let (l) and V (l) be our estimates for n+1 and V n+1 after l iterations, and let R and R V be the associated residuals in the and V equations: where (.) t quantities, including the velocities used for advection, are evaluated using the latest available estimates. Seek increments , V that will reduce the residuals towards zero: Here φ * is a reference geopotential field defined at cell edges.
In the current implementation it is updated at each time step based on φ at step n. This approach resembles Newton's method with an approximate Jacobian. Eliminating V leaves a Helmholtz problem for : A variety of methods are possible for solving the Helmholtz problem. The current implementation uses a single sweep of a full multigrid method (e.g. Fulton et al., 1986), which provides more than sufficient accuracy in the context of the iterative nonlinear solver. (There is growing evidence that the solution of elliptic subproblems need not be a barrier to scalability on massively parallel computers, e.g. Heikes et al. (2013); Müller and Scheichl (2014).) Once is found, back-substitution gives V , and then the latest estimates for and V are updated: (l+1) The first guess is given by the value at the previous step: = V n . In the tests described below, the residuals R and R V typically decrease by an order of magnitude at each iteration (by up to two orders of magnitude per iteration at very high resolution with short time step). Four iterations are sufficient to achieve stable results, and it might be feasible to use fewer iterations operationally.

Advection scheme
The advection scheme is a so-called forward-in-time scheme, a kind of finite volume scheme. (The approach is also referred to as "swept area" or "incremental remapping".) The time integral of the flux across a cell edge is replaced by a spatial integral of the advected field over the area swept across the edge during one time step. The idea is an extension to two dimensions and more general grids of the schemes described by Crowley (1968), Tremback et al. (1987), Leonard (1979) and Leonard et al. (1993Leonard et al. ( , 1995. Similar schemes have been described by Thuburn (1997), Lashley (2002), Lipscomb and Ringler (2005), Miura (2007), and Skamarock and Menchaca (2010). The swept area integral is computed by approximating the subgrid distribution of the advected field as a polynomial in terms of local x and y coordinates. The code has been implemented for an arbitrary degree polynomial (though with some approximations); the results below focus on the case of a quadratic fit. Some details are described in the following subsections.

Construction of stencils
A polynomial subgrid distribution of the advected field is constructed for each grid cell. A polynomial of degree d in two dimensions has (d + 1)(d + 2)/2 coefficients, and so requires at least this many pieces of information in order to determine the coefficients, i.e. we need a stencil of at least this size. The stencil should be as isotropic and symmetrical as possible to ensure that the construction of the fit is well conditioned, and because the same subgrid reconstruction is used for all the downstream edges of the given cell. The stencil is grown iteratively, as follows: -First sweep: the stencil comprises only the cell in question.
-Subsequent sweeps: if there are sufficient cells in the stencil then stop. If not then make a list of all cells that are not yet in the stencil but are neighbours of cells in the stencil. If any cells in the list are neighbours of two or more stencil cells then add these to the stencil and finish this sweep. Otherwise, add all the cells in the list (which are neighbours of only one stencil cell) to the stencil and finish this sweep. Figure 3 illustrates how the stencil is grown on three different grids in order to fit quartic polynomials, which need at least 15 stencil cells. The numbers indicate the number of the sweep in which the cell is added to the stencil. (The triangular case is relevant to PV advection on the dual of the hexagonal grid.) Note that the details of the stencil vary near anomalous regions of the grid such as pentagons or cube corners; these cases are left as an exercise for the interested reader!

Constructing the polynomial fit
We wish to construct a polynomial fit for cell i so as to fit the given data, which are grid cell area integrals j in each of the corresponding stencil cells. However, the stencil constructed using the above scheme will almost always contain more cells than needed to determine the polynomial coefficients; we have an overdetermined problem. An obvious solution is to construct a least squares fit to the data. However, Lashley (2002) found that this gave unstable results. He obtained stable results by demanding that the central cell be fitted exactly, with a least squares fit to the rest of the data.
Here we generalize this idea by demanding an exact fit to the data in some substencil containing the central cell, with a least squares fit to the rest of the data.
For any cell i, let B k (x) be a set of basis functions for the subgrid fit, so that where the coefficients a k are to be determined. (Later the basis functions will be chosen to be monomials, 1, x, y, x 2 , . . . etc, but for now we keep the theory general.) It is then a straightforward exercise in linear algebra to show that for some matrix C. For each cell i, the matrix C can be computed provided we can evaluate the integral of B k (x) over each stencil cell. Moreover, the C's do not change in time, so they can be evaluated just once at the start of the integration and stored for later use. One matrix C needs to be stored for each cell. Typical sizes of C per cell are 6 × 7 for a quadratic fit on a hexagonal grid (up to 7 stencil cells to determine 6 coefficients) and 15 × 19 for a quartic fit on a hexagonal grid (up to 19 stencil cells to determine 15 coefficients).
There is some freedom in the choice of substencil to be fitted exactly. All of the results shown below fit only the central cell exactly.

Local coordinate system
In the current scheme the basis functions B k (x) are monomials in local coordinates x and y. We choose a local coordinate system that approaches Cartesian as the ratio of grid length to Earth's radius tends to zero, and for which the results are independent of the choice of direction of the x axis.
The origin of the coordinate system is taken to be the centre of the cell in question, x 0 , say. The direction of the x axis is defined to be the direction from x 0 to the centre of an arbitrarily chosen neighbouring cell. This choice is simply a matter of convention; identical results should be obtained for other choices. Then for any point x in the neighbourhood of x 0 , it is straightforward to compute the spherical distance s between x 0 and x, and the angle θ between the x axis and the line joining x 0 to x. Finally, the local coordinates are given by x = s cos θ; y = s sin θ. (26)

Integrals of monomials
The basis functions B k are chosen to be the set of monomials in x and y up to some chosen degree d, for example, {1, x, y, x 2 , xy, y 2 } for d = 2. (It may be verified that the space of functions spanned by the basis is independent of the choice of x axis.) For each cell i, the construction of the subgrid fit requires the integral L j k of the kth basis function over the j th stencil cell, for all k and j . These are computed as follows.
The stencil cell is subdivided into subtriangles by joining its centre to each vertex. For each subtriangle, three approximate Gauss points are found by computing x g = (4x i +x j + x k )/6, where {i, j, k} is a cyclic permutation of the indices of the three subtriangle vertices, and then projecting x g back onto the unit sphere. The corresponding Gaussian weights are given by 1/3 of the true spherical triangle area. The integrals of the monomials over the subtriangle are thus evaluated by (approximate) Gaussian integration. The integrals over the stencil cell are obtained by summing over the corresponding subtriangles.
Three point Gaussian integration would be exact for integration of polynomials up to degree 2 for planar triangles. For spherical triangles it remains very accurate. Three point Gaussian integration is not exact for higher degree polynomials, even for planar triangles. However, because the stencil cell is subdivided into several triangles, the approximate integrals are still very accurate for d = 4.

Construction of swept area
The area swept across an edge during a time step is approximated as a parallelogram in the local x-y coordinate system of the upwind cell. The displacement in the normal direction is given by u t and in the tangential direction by v t , where u and v are the normal and tangential velocity components, making appropriate allowance for sign. On the primal grid, u is obtained from U and v is obtained from U ⊥ = HV ⊥ where V ⊥ = WU . On the dual grid, u is obtained from V ⊥ and v is obtained from V .

Integral over swept area
The swept area integral is evaluated by Gaussian integration over the parallelogram area. Sufficient Gauss points are used to evaluate the swept area integral exactly on a plane: 2 × 2 for a quadratic subgrid fit, 3 × 3 for a quartic subgrid fit. (The affine transformation that transforms a rectangle into a parallelogram also transforms the rectangle's Gauss points into the parallelogram's Gauss points.) On the primal grid, we wish to ensure that a constant φ field remains constant in a non-divergent flow. This requires that the subgrid reconstruction of a constant φ field should be constant, and also that the swept area integrals should be proportional to the velocity fluxes, so that the mass fluxes are non-divergent. The subgrid reconstruction described above preserves a constant. To ensure non-divergent mass fluxes for non-divergent velocity the Gaussian weights are normalized by the swept area: A sw e = U t e .
One dimensional experiments and stability analysis (Sect. 6.1) revealed the need for an important modification to this swept area calculation: where (ID 2 U n ) up is the divergence at time level n in the cell upwind of edge e. This modification is essential for stability when the advection of mass is coupled to the momentum equation. It may be interpreted as allowing for the timeintegrated effect of the divergence on the swept area.

Advection of tracers: primal grid
We require two important properties of the advection scheme: i. A constant φ field should remain constant in nondivergent flow.
ii. A constant tracer mixing ratio should remain constant.
On the primal grid, the first of these properties is guaranteed by the construction of the swept area integral described above.
The prognostic variables for primal grid tracers are assumed to be "concentrations", e.g. = φγ where γ is the "mixing ratio", and they are stored as area integrals (analogous to ). Provided the same subgrid reconstruction scheme is applied to the area integrals of as is applied to , the scheme will preserve a constant mixing ratio, because the subgrid distributions of and φ will be proportional to each Geosci. Model Dev., 7, 909-929, 2014 www.geosci-model-dev.net/7/909/2014/ other. (A flux limiter could also be used to ensure preservation of a constant mixing ratio, e.g. Thuburn, 1996; 2013, but has not been implemented for the results shown below.)

Dual grid advection
The spatial discretization implies that there exists a dual grid mass field (v) = R that satisfies its own mass continuity equatioṅ where C ⊥ = WF . In other words, predicting the (v) directly using Eq. (29) would give the same result as predicting then diagnosing (v) using Eq. (28).
To ensure that this property continues to hold for finite time steps, and in particular that properties (i) and (ii) hold on the dual grid, we must ensure that the dual grid advection is based on the time integrated mass fluxes C ⊥ = W F . Note that C ⊥ must be built in this way using the W operator; constructing a dual cell subgrid fit from the (v) data and a swept area based on V t ⊥ would give a different result for which Eqs. (28) and (29) are not consistent. Also note that there is no need to modify the swept mass to allow for divergence, as in Eq. (27); this has already been taken into account in the calculation of F .
The mass flux C ⊥ determines a swept mass M sw e = C ⊥ e rather than a swept area. Therefore, in order to ensure preservation of a constant tracer mixing ratio, the subgrid reconstruction must be for the mixing ratio rather than the concentration, so that the swept integral is evaluated as γ dM (30) rather than dA.
The rest of the calculation is the same as for primal grid advection, except that the Gaussian weights are normalized by the swept mass rather than the swept area. Numerical tests (see interactive discussion) confirm that properties (i) and (ii) do indeed hold in practice on both the primal and dual grids. The importance of such mass-tracer consistency for modelling atmospheric transport has been discussed, for example, by Zhang et al. (2008).

Results
A variety of tests have been applied to test the stability and different aspects of accuracy of the scheme on the hexagonal-icosahedral and cubed-sphere grids. A particular focus is on the properties listed in Sect. 2.

Stability
The stability of the scheme was investigated by applying it, with appropriate simplifications, to the one-dimensional shallow water equations. Specifically, for small perturbations to simple basic states with constant φ and u, the system matrix was generated numerically in order to compute the frequencies and structures of the eigenmodes. Provided the effect of divergence on swept area is included, as in Sect. 5.6, then, with α = β = 0.5, the scheme was found to be stable for large gravity wave Courant numbers and advective Courant numbers up to about 0.75, with only very small instability growth rate for advective Courant numbers between 0.75 and 1.0. A very small off-centring, α = 0.502, was enough to obtain stability for advective Courant numbers up to 1.
All testing of the two-dimensional shallow water model presented below has used α = β = 0.5. In practice the model is found to be stable for large gravity wave Courant numbers and advective Courant numbers less than about 1.

Convergence of Laplacian
Using the operators listed in Table 3, a discrete Laplacian operator can be built for scalars defined at the centres of primal cells (ID 2 HD 1 ) and for scalars defined at the centres of dual cells (−JD 2 H −1 D 1 ). Examining the convergence of the discrete Laplacian provides a basic test of the accuracy of some of the discrete operators. Moreover, the primal grid Laplacian arises when the discrete linearized mass and momentum equations are combined to obtain a discrete wave equation, and also in the vorticity form of the expression for geostrophic balance. Thus, the accuracy of the primal grid Laplacian will influence the accuracy with which gravity wave propagation and geostrophic balance are captured.
The code was used to test the convergence with increasing resolution of the primal grid discrete scalar Laplacian applied to the function cos ϕ sin λ (here ϕ is latitude, λ is longitude). The results are given in Table 4. For the hexagonal grid, the L 2 error is almost second order while the L ∞ error is first order or a little better. For the cubed-sphere grid the L 2 error is almost first order, while the L ∞ error does not decrease with increasing resolution.
The test was repeated for the dual grid discrete scalar Laplacian. (Again see Table 4.) For the hexagonal grid the L 2 error is approximately first order while the L ∞ error appears to be close to first order until the finest resolution where convergence almost stalls. For the cubed-sphere grid the L 2 error is worse than first order, and the L ∞ error does not decrease as resolution is refined.
A vector Laplacian can also be built from the operators in Table 3 (see Thuburn and Cotter, 2012). Its convergence was found to be similar to that of the dual grid scalar Laplacian on both the hexagonal and cubed-sphere grids (not shown).

Convergence of Coriolis operator
Although the R operator is at least first-order accurate, this does not imply any guarantee of accuracy of the W operator constructed from it. On a regular hexagonal or square grid on a plane, both R and W would be second-order accurate. However, for the distorted polygons of the quasi-uniform spherical grids the convergence rate must be checked empirically. A stream function equal to cos ϕ sin λ for a rotational flow was defined and sampled at both primal cell centres (ψ) and primal vertices (ψ (v) ). Dual edge normal fluxes were then calculated both directly from the stream function (V ⊥ = D 1 ψ) and by applying the W operator to primal edge normal fluxes (V ⊥ approx = WU = −WD 1 ψ (v) ). The difference between V ⊥ approx and V ⊥ , after dividing by the lengths of the corresponding dual edges, gives a measure of the accuracy of W. The variations of the L ∞ and L 2 errors with resolution are listed in Table 5. A similar calculation was carried out for a divergent flow by defining a velocity potential equal to cos ϕ sin λ, sampled at primal cell centres (χ) and at primal vertices (χ (v) ). The L ∞ and L 2 differences D 1 χ (v) − HWHD 1 χ (after dividing by primal edge lengths) are also shown in Table 5.
On both grids the L ∞ error fails to converge to zero, while the L 2 error converges very slowly, roughly proportional to the square root of the grid spacing, or O(N −1/4 ) on a grid with N cells. This is consistent with the observation that O(1) errors are found along the lines joining the pentagons on the hexagonal-icosahedral grid and along the cube edges on the cubed-sphere grid, thus affecting O(N 1/2 ) edges. The errors on the cubed sphere are significantly larger than those on the hexagonal-icosahedral grid.

Advection
Test case 1 of Williamson et al. (1992) tests the advection scheme in isolation from the rest of the dynamics. A cosine bell profile tracer is advected once around the sphere by a solid body rotation flow. We have carried out this test on both the hexagonal-icosahedral and cubed-sphere grids, for tracers stored on both primal and dual cells, at a range of grid resolutions and for flow at different angles relative to the grid. Figure 4 shows a sample of results. Generally the advection scheme shows -accurate phase speed, with weak dispersion error; -some erosion of the maximum but with rather isotropic error field (the cubed-sphere case in Fig. 4 is rather unusual in that the flow is aligned with the grid, leading to more elongation than broadening of the tracer profile); -weak undershoots provided the tracer is well resolved; -little sensitivity to the grid (hexagonal or cubed sphere, primal or dual); -little grid imprinting, as measured by the evolution of errors as the cosine bell crosses grid features such as pentagons or cube corners; -a convergence rate close to second order or better, depending on the error norm, at the resolutions tested.
The advection scheme has also been tested using a quartic subgrid fit d = 4. This option is considerably more expensive than the quadratic fit d = 2. The qualitative features of the results are rather similar to the d = 2 case. Generally the errors are reduced in magnitude by a factor in the range 3-5, though the convergence rate remains close to second order. Undershoots are a little worse on poorly resolved data. All of the results shown in later sections use d = 2.
Several factors could reduce the convergence rate below the third order and fifth order rates expected in the ideal case for quadratic and quartic subgrid fits, respectively. These factors include the approximation of the swept area as a parallelogram (Ullrich et al., 2013), the use of approximate quadrature in the quartic case, and the lack of perfect smoothness of the initial data (Holdaway et al., 2008). The parallelogram swept area should be an excellent approximation for the solid body flow of this test case. The observed convergence rates are consistent with the results of Holdaway et al. (2008) (case n = 2 in their Table I), suggesting that smoothness of the initial data is the dominant factor limiting the convergence rate.
We have also experimented with fitting a larger substencil exactly. In most cases this leads to a reduction in the errors. However, for advection on the dual of the hexagonal grid with d = 4 the scheme became unstable, for reasons we do not yet understand. All of the results shown below use exact fitting only to the central stencil cell.  triggers the radiation of fast inertio-gravity waves and slow Rossby waves.

760
Maps of the surface height field produced by the model at day 15 appear very similar to those published in the literature for other models. To obtain more detailed information we therefore examine the errors in the height field. Test case 5 has no analytical formula for the true solution, so 765 a high-resolution reference solution was generated using the ENDGame shallow water model (Zerroukat et al., 2009), which uses a well tested semi-implicit semi-Lagrangian solution scheme on a longitude-latitude C-grid. The reference solution is generated at a grid resolution 1024×512, and a time 770 step ∆t = 225 s. Standard semi-Lagrangian advection of φ (rather than SLICE) was used. The high-resolution solution for the surface height h was interpolated to all of the desired test grids and resolutions using cubic Lagrange interpolation (the same scheme as used for semi-Lagrangian advection).

775
Each test model computed a height error field by reading in the corresponding reference solution and subtracting it from the test model solution.
As well as the mimetic finite volume scheme on the hexagonal and cubed sphere grids, errors were also calculated for 780 the ENDGame shallow water model at similar resolutions, for comparison.
All three models were found to run stably with the time steps given in Table 11. However, all three models produced

Solid body rotation
The solid body rotation test case, test case 2 of Williamson et al. (1992), tests the ability of the scheme to maintain a steady, balanced, large-scale flow. Since the flow is steady, the exact solution is known, allowing the calculation of errors and convergence rates. Table 6 shows L 2 and L ∞ errors for φ and v at day 5 vs. resolution for the two grid types. The time steps used are given in Table 1. At the resolutions tested, v converges at close to second order on both grids, while φ converges at somewhere between first and second order. At any given resolution the errors are considerably smaller on the hexagonalicosahedral grid than on the cubed-sphere grid.
Maps of the error in φ on both grids (Fig. 5) show small scale and large scale components. The small-scale errors are concentrated along the lines joining the pentagons on the hexagonal grid and along the cube edges on the cubedsphere grid. They are stationary, and are present at their full amplitude after just a few time steps. The large scale error reflects the symmetry of the grid: zonal wave number 5 and antisymmetric about the equator for the hexagonal-icosahedral grid, zonal wave number 4 and symmetric about the equator for the cubed-sphere. This component is also stationary; it gradually emerges over about three days of integration.

Flow over an isolated mountain
A more complex flow field is produced in test case 5 of Williamson et al. (1992). An initial geostrophically balanced solid body rotation velocity field impinges on an isolated conical mountain at northern mid-latitudes. The mountain triggers the radiation of fast inertio-gravity waves and slow Rossby waves.
Maps of the surface height field produced by the model at day 15 appear very similar to those published in the literature for other models. To obtain more detailed information we therefore examine the errors in the height field. Test www.geosci-model-dev.net/7/909/2014/ Geosci. Model Dev., 7, 909-929, 2014  case 5 has no analytical formula for the true solution, so a high-resolution reference solution was generated using the ENDGame shallow water model (Zerroukat et al., 2009), which uses a well tested semi-implicit semi-Lagrangian solution scheme on a longitude-latitude C-grid. The reference solution is generated at a grid resolution 1024×512, and a time step t = 225 s. Standard semi-Lagrangian advection of φ (rather than SLICE) was used. The high-resolution solution for the surface height h was interpolated to all of the desired test grids and resolutions using cubic Lagrange interpolation (the same scheme as used for semi-Lagrangian advection). Each test model computed a height error field by reading in the corresponding reference solution and subtracting it from the test model solution. Table 6. Geopotential errors (m 2 s −2 ) and velocity errors (m s −1 ) at day 5 for the solid body rotation test case. As well as the mimetic finite volume scheme on the hexagonal and cubed-sphere grids, errors were also calculated for the ENDGame shallow water model at similar resolutions, for comparison.

Grid & Cells
All three models were found to run stably with the time steps given in Table 1. However, all three models produced almost identical height error fields (not shown), despite the use of very different grids and numerics. Moreover, convergence of the errors towards zero stalled between the two finest resolutions tested, for all models and all error norms. The common feature of the numerics of the three models is the semi-implicit treatment of gravity waves, which will artificially slow the highest frequency waves. Such waves are present at large amplitude because of the "impulsive" start to the test case. The evidence suggests that this slowing of gravity waves is the dominant source of error in all three models.
Geosci. Model Dev., 7, 909-929, 2014 www.geosci-model-dev.net/7/909/2014/ It is an encouraging result that the mimetic finite volume model produces such similar errors to ENDGame when using time steps of the size that would be used in practice at these resolutions. However, it is also of interest to try to assess the errors arising from the spatial discretization and the grid. To do this, the test case was repeated for all three models with t reduced by a factor of 4 on each grid in order to reduce the time truncation errors. The resulting error norms for height are shown in Table 7. All three models appear to be converging at a rate somewhere between first and second order, depending on the norm chosen. (Note that, because of the lack of smoothness of the forcing mountain, we cannot expect better than first-order convergence of L ∞ (h), even with a high-order scheme.) At any given resolution, the hexagonal and cubed-sphere grids have rather similar errors, while those from ENDGame are typically (though not always) slightly smaller.
The height error fields for the three models at the second highest resolution are shown in Fig. 6. Even with the reduced time step used here, there are noticeable similarities among the three models in both the length scales and the details of the error patterns, which are suggestive of wave trains radiating globally from the forcing mountain.
In order to examine some aspects of behaviour that depend on the mimetic properties of the scheme, the test case was run to 50 days at the highest resolution given in Table 7 (with t = 900 s). After about 20 days the PV dynamics becomes strongly nonlinear, leading to the production of thin PV filaments which stretch and wrap up, localized sharp PV gradients, and a cascade of potential enstrophy to small scales. Figure 7 shows the PV at day 50 on the cubed-sphere grid. It shows the production of PV filaments and sharp gradients, with no noise or other unphysical behaviour apparent. Results on the hexagonal grid are similar. Figure 8 shows several diagnostics of behaviour related to the mimetic properties, again for the cubed-sphere case. (The hexagonal grid case is very similar.) The top left panel shows the relative change in total mass, and confirms that this is at the level of round-off error. Equation (29) implies that advecting a dual grid tracer initialized with the dual grid geopotential should give the same result as diagnosing the dual grid geopotential from the predicted primal grid geopotential at each time step. The continuous curve in the top right panel shows the maximum absolute difference between these two dual grid φ fields, normalized by the maximum value of the field. The normalized difference is of order 10 −5 . Although the difference is tiny, it is significantly larger than round-off error. This is due to incomplete convergence of the nonlinear iterative solver; increasing the number of iterations from 4 to 8 reduces the discrepancy by a factor of 10 3 . The construction of the Q ⊥ Coriolis term ensures that the PV diagnosed from the predicted and V fields evolves as if it were a passive tracer advected by the mass fluxes F ⊥ . The dashed curve in the top right panel shows the normalized maximum absolute difference between PV diagnosed at each time step and an advected tracer initialized with the PV field. The differences are at the level of round-off error. The available potential energy is given by where φ T is the global mean of φ T . It gives an upper bound on the amount of potential energy that could be converted to kinetic energy, and is typically much smaller than the total potential energy. The bottom left panel of Fig. 8 shows the available potential energy, kinetic energy and their sum (total available energy). There is a significant conversion of available potential energy to kinetic energy over the 50 days, but their sum is well conserved. The bottom right panel shows the relative change in total available energy (continuous curve) and also the relative change in total potential enstrophy (dashed curve). As discussed in Sect. 2, the numerical methods are designed not to conserve these quantities exactly but to allow a transfer to unresolved scales when there are significant nonlinear cascades. In fact, during the first 15 days of this test case the flow is only weakly nonlinear and any downscale cascades are rather weak. During this time the losses of available energy and potential enstrophy are very small, of order 1 part per thousand. Very similar results are obtained for ENDGame (see interactive discussion).

Barotropically unstable jet
The test case described by Galewsky et al. (2004) produces a rapidly growing barotropic instability from a strong, narrow mid-latitude zonal jet. The true solution at day 6 has the instability localized over a certain range of longitudes with part of the jet remaining almost quiescent. Schemes with significant grid imprinting tend to excite the instability all along the jet, possibly with an incorrect zonal wave number. The test case is dominated by strongly nonlinear PV advection with rapid generation of small scales by straining and vortex roll-up. A good scheme should produce no small-scale noise or unphysical rippling in the PV or vorticity field. Since the features of interest can easily be identified from maps of vorticity, full fields rather than errors are presented. Finer resolution is needed for this test case than the earlier ones. Also, since the maximum velocity is around 80 m s −1 , the time step is halved for all models and grids compared to the values given in Table 1. Figure 9 shows the vorticity field at day 6 on the hexagonal grid with 10 242 and 163 842 cells, the cubed-sphere grid with 13 824 and 221 184 cells, and, for reference, from ENDGame with 640 × 320 cells. The ENDGame solution is similar to other published high resolution solutions (e.g. Li and Xiao, 2010; Salehipour et al., 2013). All of the results show a clean vorticity field, free from noise and spurious ripples. However, at the coarser resolutions shown, the finite volume model shows strong grid imprinting, with the wave number of the instability determined by the grid structure. At the finer resolutions shown, the finite volume model solutions are more similar to the reference solution, but both show significantly stronger development of the instability over the longitude range [π/2, π ] than the reference solution, implying that truncation errors are still significant even at these fine resolutions. Similar grid-triggering of instability is seen in the Rossby-Haurwitz wave test case (Williamson et al., 1992) (not shown). The Rossby-Haurwitz wave is a quasi-steadily propagating wave number 4 pattern that is actually unstable, but the time taken for the pattern to break down depends on how strongly numerical truncation or roundoff errors project onto the unstable modes, which involve wave numbers 1, 3 and 5 (Thuburn and Li, 2000). The five-fold symmetry of the hexagonal grid clearly has potential to trigger the instability, and in practice we found that on the hexagonal grid with 10 242 cells the Rossby-Haurwitz wave pattern breaks down after about 20 days. Thus, the mimetic properties of the finite volume scheme do not help to reduce this kind of gridtriggering of instability in delicately balanced initial states.

Balance
One of the main motivations for the development of the mimetic numerical method is the requirement for the numerical model to respect balance. In the regimes of small Rossby or Froude number, the generation of imbalance in the form of fast waves from an initially balanced flow should be very weak (e.g. Ford et al., 2000;Cullen, 2000). A numerical model should not generate excessive imbalance, and, ideally, should not require artificial damping mechanisms to control imbalance. A thorough investigation of this issue requires an examination of how the model performs in the asymptotic limits of small Rossby or Froude number (Cullen, 2007(Cullen, , 2008; this is the subject of ongoing work. For present purposes we examine some simple diagnostics of imbalance in the barotropic instability test case, applied to the mimetic finite volume model and, for comparison, to ENDGame. For each model a series of three integrations was run. a. The model was integrated with a centred semi-implicit time integration scheme (α = 0.5). The initial perturbation applied to the jet in the barotropic instability test case is unbalanced. With α = 0.5 it is found that the fast waves resulting from the initial perturbation dominate the divergence field throughout the integration to day 6. b. The model was run again using a fully off-centred scheme (α = 1.0). This is found to suppress the fast waves within about one day. The divergence pattern is now very different: it is collocated with the vortex rolls of the evolving barotropic instability and grows along with them, confirming that it is slaved to the balanced dynamics rather associated with fast waves. However, the use of α = 1.0 could be suppressing any numerically generated imbalance.
c. Therefore, for the third integration we set α = 1.0 for the first 500 steps (a little over 31 h) and α = 0.5 thereafter. This removes the imbalance associated with the initial perturbation, but would avoid artificially suppressing any imbalance (physical or numerical) generated subsequently as the instability grows. Figure 10 shows the divergence field at day 6 from the mimetic finite volume model on the hexagonal grid for the three integrations (a), (b), and (c). It confirms the behaviour described above, and shows that the results of runs (b) and (c) are almost identical. Very similar behaviour is found on the cubed-sphere grid and for ENDGame. Figure 11 shows time series of the root-mean-square divergence from the three models. It shows that the evolution of the divergence is very similar in runs (b) and (c) for all three models. These results confirm that any generation of imbalance by the mimetic finite volume scheme is extremely weak and is comparable to that in ENDGame.

Computational modes
Dispersion analysis for the regular hexagonal C-grid on a plane (Thuburn, 2008) shows that it supports an extra family of Rossby modes. These are characterised by small-scale vorticity and PV structures. Although they correctly have zero frequency on an f plane (in zero background flow), on a β plane their frequencies are anomalously small and strongly sensitive to the detailed formulation of the Coriolis operator. These unphysical aspects of their behaviour lead to concerns that these extra modes might adversely affect solutions, for example through the appearance of noise or through an incorrect response to forcing.
The following argument suggests that the extra modes should be harmless provided PV advection is well handled (see also Weller, 2012). Small scale Rossby waves have very small intrinsic frequency. In the presence of a background flow their absolute frequency is dominated by advection and Doppler shifting. Although these extra Rossby modes have excessively small intrinsic frequency, their absolute frequency will be quite accurate provided PV advection is adequately captured. Moreover, in order to avoid dispersion errors, an advection scheme must be inherently damping on small scales, so advection by a background wind should tend to damp the extra Rossby modes. The mimetic finite volume model discussed here is designed to have good PV advection properties, so we expect the above argument to hold.
To test this argument, test case 5 of Williamson et al. (1992) was run to day 15 on a hexagonal grid of 10 242 cells, then two patches of grid scale noise in the vorticity field were superposed on the solution. The noise patches were generated by starting with a zero stream Geosci. Model Dev., 7, 909-929, 2014 www.geosci-model-dev.net/7/909/2014/ function, introducing a "seed" delta function in the desired regions, applying a number of iterations of antidiffusion but with extrema limited to some maximum value to grow the patch, constructing a rotational velocity perturbation from the stream function, normalizing the maximum velocity perturbation to 1 m s −1 , and adding to the model's velocity field. The pattern generated in the vorticity field this way has structure very similar to that of the smallest scale extra Rossby modes. One patch of noise was introduced on the equator near (π, 0) and the other near the north pole. The evolution of the vorticity field over the next few time steps is shown in Fig. 12. The equatorial patch is located in a region of relatively strong wind, and it is found to be rapidly damped, almost completely disappearing within 2 h. The polar patch is located in a region of relatively weak wind. It is damped more slowly, and preferentially on the side that experiences slightly stronger wind. Nevertheless, after 24 h the noise has been almost completely removed.

Conclusions
A new finite volume shallow water model on the sphere has been described. The formulation allows the use of arbitrary polygonal grids, motivated by the need for quasi-uniform spherical grids to enable parallel scalability. Results of several test cases have been presented for two particular grids: a hexagonal-icosahedral Voronoi grid and a modified equiangular cubed-sphere grid.
The model uses a new time integration scheme (Sect. 3). It combines a semi-implicit treatment of fast waves with a forward-in-time advection scheme for mass and for the potential vorticity fluxes that appear in the velocity equation. Provided the advective swept areas are modified to allow for the effect of divergence (Eq. 27), this scheme is found to be stable and robust for advective Courant numbers up to about one and for wave Courant numbers greater than one, without the need for off-centring of the semi-implicit part of the scheme or other additional damping mechanisms.
The scheme is built around a framework that allows it to mimic key mathematical properties of the continuous equations that underpin important physical properties such as mass conservation, linear energy conservation, an accurate representation of balance and potential vorticity dynamics, and consistent advection of mass, PV and tracers. We have presented diagnostics (e.g. Sects. 6.6 and 6.8) confirming that these properties are obtained in practice.
The gradient, divergence and curl operators D 1 , D 2 , D 1 , and D 2 are exact within this finite volume framework, while the other operators I, J, H, and R are at least first-order accurate. However, a limitation of the scheme presented here is that the Coriolis operator W, which is crucial to some of the mimetic properties, is not, in fact numerically consistent (Sect. 6.3). Nevertheless, the model solution does appear to be converging at the resolutions tested in the idealized solid body rotation test case (Sect. 6.5), though the error patterns clearly reflect the grid structure. For the more complex flow of the isolated mountain test case the errors are comparable to those of a state-of-the-art semi-implicit semi-Lagrangian model on a longitude-latitude grid. For the barotropic instability test case (Sect. 6.7) truncation errors spuriously trigger the development of the instability even at quite fine resolution, though the solution does appear to be approaching the reference solution as resolution is increased. Thus, although the use of an inconsistent operator is unappealing, it is not clear that it will prevent convergence in flows of realistic complexity. The errors for both the Laplacian operator and the Coriolis operator are significantly larger on the cubed-sphere than on the hexagonal-icosahedral grid. The solution errors in the solid body rotation test are also significantly larger on the cubed-sphere than on the hexagonal-icosahedral grid. However, for the isolated mountain test case the errors on the two grid types are comparable. In this case the errors appear to be dominated by time truncation errors associated with the implicit treatment of gravity waves, rather than the spatial discretization. Finally, both grids appear to show a similar degree of spurious development in the barotropic instability test case, though some features are slightly better captured on the hexagonal grid.
In the test cases carried out, there is no evidence for any damaging effects of the extra Rossby modes supported on the hexagonal C-grid. When grid scale vorticity features, which should project strongly onto the extra Rossby modes, are forced into the solution, the model remains robust and the numerics are able to remove the noise on the advective timescale.
Because of the requirement of symmetry for the I, H and J operators, the construction of higher-order versions appears to be difficult. The construction of a consistent W operator appears to be even more difficult. The construction of W from R with the desired mimetic properties is only known for the case in which the stencil for R is the set of primal cells overlapped by each dual cell; and in this case the resulting W is uniquely determined . Motivated by these apparent limitations on the mimetic finite volume scheme, we are currently investigating the use of a finite element approach (Cotter and Shipton, 2012), which can give the same mimetic properties as the scheme presented here but with more accurate basic operators. This will help to determine whether the accuracy of the basic operators is indeed a limiting factor for solution accuracy.
The Supplement related to this article is available online at doi:10.5194/gmd-7-909-2014-supplement.