Total energy and potential enstrophy conserving schemes for the shallow water equations using Hamiltonian methods: Derivation and Properties (Part 1)

The shallow water equations provide a useful analogue of the fully compressible Euler equations since they have similar characteristics: conservation laws, inertia-gravity and Rossby waves and a (quasi-) balanced state. In order to obtain realistic simulation results, it is desirable that numerical models have discrete analogues of these properties. Two prototypical examples of such schemes are the 1981 Arakawa and Lamb (AL81) C-grid total energy and potential enstrophy conserving scheme, and the 2007 Salmon (S07) Z-grid total energy and potential enstrophy conserving scheme. Unfortunately, the AL81 scheme is restricted to logically square, orthogonal grids; and the S07 scheme is restricted to uniform square grids. The current work extends the AL81 scheme to arbitrary non-orthogonal polygonal grids and the S07 scheme to arbitrary orthogonal spherical polygonal grids in a manner that allows both total energy and potential enstrophy conservation, by combining Hamiltonian methods (work done by Salmon, Gassmann, Dubos and others) and Discrete Exterior Calculus (Thuburn, Cotter, Dubos, Ringler, Skamarock, Klemp and others). Detailed results of the schemes applied to standard test cases are deferred to Part 2 of this series of papers.


Introduction
Consider the motion of a (multi-component) fluid on a rotating spheroid under influence of gravity and radiation. This is the fundamental subject of inquiry for geophysical fluid dynamics, covering fields such as weather prediction, climate dynamics and planetary atmospheres. Central to our current understanding of these subjects is the use of numerical models to solve the otherwise intractable equations (such as the fully compressible Euler equations) that result. As a first step towards developing a numerical model for simulating geophysical fluid dynamics, schemes are usually developed for the rotating shallow water equations (RSWs). The RSWs provide a useful analogue of the fully compressible Euler equations since they have similar conservation laws, many of the same types of waves and a similar (quasi-) balanced state. It is desirable that a numerical model posses as least some these same properties (see Figure 1, and the discussion in [26]).
In fact, there exists some evidence ( [4]) that schemes without the appropriate conservation properties can fail to correctly capture long-term statistical behaviour, at least for simplified models without any dissipative effects. However, questions remain as to the relative importance of various conservation properties for a full atmospheric model, especially in the presence of forcing and dissipation ( [28]). This subject deserves further study, but a key first step is the development of a numerical scheme that posses the relevant conserved quantities; and is capable of being run at realistic resolutions on the types of grids that are used in operational weather and climate models.
A pioneering scheme developed over 30 years ago possesses many of these properties (including both total energy and potential enstrophy conservation): the 1981 Arakawa and Lamb scheme (AL81, [1]). Unfortunately, this scheme is restricted to logically square, orthogonal grids such as the lat-lon or conformal cubed-sphere grid. These grids are not quasi-uniform under refinement of resolution, and this leads to clustering at typically target resolutions for next generation weather and climate models (such as 2-3km for weather; and 10-15km for climate). Such clustering will introduce strong CFL limits, and in the case of the lat-lon grid requires polar filtering (which is not scalable on current computational architectures) in order to take realistic time steps. For these reasons, it is desirable to be able to use quasi-uniform grids such as the icosahedral grid (orthogonal but non-square) or gnomic cubed-sphere (square but non-orthogonal). In addition to the restriction to logically square, orthogonal grids, the AL81 scheme also suffers from poor wave dispersion properties when the Rossby radius is underresolved ( [17]). In fact, the unavoidable averaging required for the Coriolis term in a C grid scheme is expected to lead to poor wave dispersion properties for an underresolved Rossby radius regardless of the specific discretization employed.
Recently, there has been an effort to extend the AL81 scheme to more general grids, using tools from discrete exterior calculus (commonly referred to as the TRiSK scheme, [31], [18], [29], [35], [30]). This has lead to the development of a family of schemes on general non-orthogonal (spherical) polygonal meshes that posses all of the desirable properties of AL81 except for: extra modes branches on non quadrilateral meshes, which are unavoidable for C grid schemes; and lack of either total energy or potential enstrophy conservation. It is possible to obtain one or the other, but not both at the same time. Along different lines, Salmon ([20]) showed that AL81 and other doubly-conservative schemes (such as [27]) are all members of a another family of schemes on logically square orthogonal meshes. This was done using tools from Hamiltonian methods, which are an area of active research in atmospheric model development.
As an alternative to the AL81 scheme that preserves many of its valuable mimetic properties, but has good wave dispersion properties independent of Rossby radius, [17] introduced a scheme for uniform square grids based on the vorticity-divergence formulation (termed the Z grid) of the continuous equations. Subsequently, this approach was extended to arbitrary (spherical) orthogonal polygonal grids with a triangular dual in [12] and [13], which included the important case of an icosahedral-hexagonal grid. Although this scheme posses many of the desirable properties from AL81, it does not conserve total energy or potential enstrophy. However, a similar Z grid scheme based on a Helmholtz decomposition of the momentum instead of the wind that does conserve both total energy and potential enstrophy was developed by Salmon ([22], [21]) using techniques from Hamiltonian mechanics (specifically, Nambu brackets). The idea of using Hamiltonian mechanics to derive conservative models for atmospheric dynamical cores has seen a great deal of interest and progress in the past 10 years (see ([11], [10], [25], [16][6], [5], [34], [19], [24]).). With the recent development of Hamiltonian formulations for essentially all of the equation sets and vertical coordinates used in atmospheric dynamics, it seems likely that this approach will continue to be employed in the future. Unfortunately, the scheme in S07 is defined only for planar grids, and in the key case of general polygonal grids no expression for discrete Hamiltonian or Casimirs was given. This precludes its further development for implementation into an operational dynamical core.
This work combines the discrete exterior calculus approach from [29] and the Hamiltonian approach from [20] to extend AL81 to general non-orthgonal (spherical) polygonal grids in a manner that conserves both total energy and potential enstrophy; and to extend S07 to arbitrary (spherical) orthogonal polygonal grids. The extension of AL81 is done through the development of a new Q (the discretization of qkˆ, which is also known as the nonlinear potential vorticity flux) operator, using tools from Hamiltonian methods. S07 is extended by combining the Nambu bracket based approach from [22] with the discrete exterior calculus tools introduced in [29]. It should be noted that this work deals only with spatially conservative discretization. Conservation errors introduced due to time discretization are typically much smaller than those due to space discretization. However, the extension of this approach to fully conservative discretization would be a useful contribution.
The remainder of this paper is structured as follows: Section 2 introduces the rotating shallow water equations in both their familiar vector-invariant form and the less familiar Hamiltonian forms. Section 3 presents a family of C grid numerical schemes that posses many of the desirable properties, and discusses the specific member of this family introduced here. Section 4 introduces the new operator Q that enables the conservation of both total energy and potential enstrophy in the C grid scheme. Section 5 presents the Z grid scheme and discusses its key mimetic and conservation properties. Finally, some conclusions (Section 6) Figure 1: A diagram of some desirable model properties for the shallow water equations, organized thematically into groups. Similar considerations apply for the Euler, hydrostatic primitive and other equation sets used in atmospheric models. There is vigorous discussion in the literature and between model designers about the importance of various properties for different applications (such as weather forecasting or longterm climate prediction). The schemes presented here satisfy all of these properties, with the exception of accuracy. There are additional desirable model properties, such as consistent physics-dynamics coupling, compatible and accurate tracer advection, and tractable treatment of acoustic waves that are not presented. are drawn. The appendices discuss various ancillary topics such as the computational grid used (Appendix A), the specific discrete operators employed (Appendices B, C and D), and the discrete variables used in the C and Z grid schemes (Appendices E and F).

Rotating Shallow Water Equations
The rotating shallow water equations (RSWs) for both planar and spherical domains are presented below in several forms: the vector invariant formulation, the vorticity-divergence formulation, the symplectic Hamiltonian formulation based on the vector-invariant form and both Poisson bracket and Nambu bracket formulations based on the vorticity-divergence formulations. Although all of these formulations are equivalent in the continuous case, they lead to very different discretizations.

Vector Invariant Formulation
The mass continuity equation for the RSWs is expressed in vector invariant form as: where h is the fluid height and u is the fluid velocity. Similarly, the momentum equation is expressed as: where F " h u is the mass flux, q " η h is the potential vorticity, η " ζ`f is the absolute vorticity, ζ "k¨ ∇ˆ u is the relative vorticity, f is the Coriolis force, Φ " gh`K`gh s is the Bernoulli function, h s is the topography height, g is gravity and K " u¨ u 2 is the kinetic energy.

Poisson Bracket Formulation (Vector Invariant)
As discussed in [20], let the Hamiltonian H be given by and x " ph, uq. Then the time evolution of an arbitrary functional F can be written as where the Poisson bracket tF, Hu (which is a bilinear, antisymmetric operator that satifies the Jacobi identity) is It is useful to split this into two separate brackets as tF, Hu " tF, Hu R`t F, Hu Q where tF, Hu R " encompasses the gradient and divergence terms; and tF, Hu Q " encompasses the nonlinear PV flux term. The functional derivatives δH δ x of the Hamiltonian are given by This formulation is useful for development of a scheme that posses discrete conservation properties, as discussed below. A functional derivative of some functional Fr xs is defined as

Conserved Quantities
Since the rotating shallow water equations form a (non-canonical) Hamiltonian system, we know from Noether's theorem and other considerations (such as the singular nature of the symplectic operator) that there are at least two categories of conserved quantities: Hamiltonian and Casimirs.

Energy (Hamiltonian)
The first is simply the Hamiltonian itself. In this case, the Hamiltonian is the total energy of the system. Conservation of the Hamiltonian arises due to the skew-symmetric nature of the Poisson bracket. In particular, using (4) the evolution of H is given by dH dt " tH, Hu "´tH, Hu " 0 (11) since t, u is skew-symmetric. For the rotating shallow water equations, the Hamiltonian is the total energy of the system. The elegant derivation of energy conservation and its simplicity (relying ONLY on the skewsymmetry of t, u) motivates the use of the Hamiltonian formulation for development of numerical schemes that conserve energy.

Casimirs
The second category of conserved quantities consists of Casimir invariants. Since the rotating shallow water equations are a non-canonical Hamiltonian system, the Poisson bracket t, u is singular and thus it possesses Casimir invariants C that satisfy tF, Cu " 0 for any functional F. Note that from above, this implies that For the rotating shallow water equations, the Casimirs take the form where F pqq is an arbitrary function of the potential vorticity and Important cases include F " 1 (mass conservation), F " q (circulation or mass-weighted potential vorticity) and F " q 2 2 (potential enstrophy).

Vorticity-Divergence Formulation
By taking the divergence ( ∇¨q and curl ( ∇ K¨) of (2), we obtain the vorticity-divergence form of the equations: where µ " ∇¨ u is the divergence. The mass flux can then be split into rotational and divergent components (ie a Helmholtz decomposition) as: where ph uq div " ∇χ and ph uq rot " ∇ K ψ. The streamfunction ψ and velocity potential χ can be related to the vorticity and divergence as µ " ∇¨ph´1 ∇χq`Jpψ, h´1q (20) where Jpa, bq " ∇¨pa ∇ T bq " ∇ T¨p a ∇bq is the Jacobian operator. The Hemholtz decomposition connects the vorticity-divergence formulation and the vector invariant formulations. In the preceding, we have neglected the possibility of a harmonic component (a component A for which ∇ 2 A " 0), which works because the harmonic component on the sphere is zero. On the doubly periodic plane, it would be possible to have a constant harmonic component. Finally, (1) and (2) can be re-written in terms of χ and ψ directly as Bζ Bt " Jpq, ψq´ ∇¨pq ∇χq Bµ Bt " Jpq, χq` ∇¨pq ∇ψq´∇ 2 Φ (23)

27)
(this is the functional derivative of the Hamiltonian with respect to x). Also define a Poisson bracket (which is bilinear, anti-symmetric and satisfies the Jacobi identity) as tA, Bu " tA, Bu µµ`t A, Bu ζζ`t A, Bu µζh (28) where tA, Bu ζζ " tA, Bu µµ " for arbitrary functionals A and B. As before, the time evolution of an arbitrary functional A is then given by dA dt " tA, Hu It is easy to see that (21), (22) and (23) are recovered when A is set equal to h,ζ or µ, respectively. Note that each of the brackets (29), (30) and (31) are anti-symmetric, and that the Casimirs C " ş Ω hF pqqdΩ satisfy tA, Cu " 0 (where F is an arbitrary function and A is an arbitrary functional) independently for each bracket. The use of the Poisson (and Nambu) bracket formulation of the shallow water equations is motivated by the intimate connection between these formulations and the conserved quantities. As is well-known, the conservation of energy H rests solely on the anti-symmetry of the Poisson bracket, and a numerical scheme that retains this feature will automatically conserve energy. However, potential enstrophy is a Casimir, and therefore developing a numerical scheme using the Poisson formulation that conserves it requires that the discrete potential enstrophy lies in the null space of the resulting discrete bracket. This can be difficult, especially on arbitrary grids, and this motivates the use of a continuous formulation that does not contain a null space, which is discussed below.

Nambu Bracket Formulation (Vorticity-Divergence)
Fortunately, there is a closely related formulation of the shallow water equations in terms of Nambu brackets (see [22]): where cyc is a cyclic permutation, Z " ş Ω dΩh q 2 2 is the potential enstrophy, and the multipart dot product is simply the product of the individual components, summed over each basis (for example, in 2D doubly periodic flow the first term is These brackets are useful because they are triply anti-symmetric (which ensures the conservation of H and Z) and non-degenerate (they have no Casimirs). In fact, discrete conservation of both total energy and potential enstrophy requires only the triply anti-symmetric nature is retained. It is also possible to generalize these brackets to ANY Casimir (as shown in [21]), but since we are interested mostly in potential enstrophy conservation this is not necessary. These brackets will form the basis of the Z grid discretization method discussed below.

C Grid Scheme
Following [29], the prognostic variables for the C grid scheme are the mass primal 2-form m i and the wind dual 1-form u e . These are naturally staggered, since primal 2-forms are associated with primal grid cells and dual 1-forms are associated with dual grid edges. Letting x " pm i , u e q, the vector-invariant Poisson bracket can be discretized in a manner that preserves its anti-symmetric character (which ensures total energy conservation) and a subset of the Casimir invariants (specifically: mass, potential vorticity and potential enstrophy). Combined with a choice for the discrete Hamiltonian, this constitutes a complete discretization for the nonlinear rotating shallow water equations. Ideally, one would use a Nambu bracket formulation of the vector invariant shallow water equations rather than the Poisson bracket formulation in order to avoid the difficulties associated with developing a discretization that has the correct Casimirs, since in the Nambu bracket case only anti-symmetry must be enforced. Unfortunately, the only known Nambu bracket for the vector invariant shallow water equations possesses intractable singularities and is not suitable as the basis for developing a discretization ( [32]).
Specifically, the brackets 7 and 8 are discretized using the operators from Appendices C and B as: where the discrete functionals (such as A) are expressed as inner products using the Hodge stars. Note that these discrete brackets are only bilinear and anti-symmetric, they do not satisfy the Jacobi identity. In addition, they posses only a subset of the Casimirs of the continuous brackets. Therefore they should be properly be termed quasi-Poisson brackets. The brackets given in (37) and (38) are essentially a generalization of the brackets introduced in S04 from uniform square grids to arbitrary polygonal grids, using operators from discrete exterior calculus. The discrete function derivative with respect to a particular discrete form is the corresponding dual form. For example, consider F " pA i , B i q I , where A i and B i are primal 2-forms. Then δF δAi " IB i , which is a dual 0-form. The Hamiltonian H is discretized as: where g is the acceleration due to gravity, C e " m e u e and m e " φIm i . Taking functional derivatives yields where Φ i is the Bernoulli function dual 0-form and F e is the mass flux primal 1-form. Computing actual values yields: , where b i is the topographic height primal 2-form and K i is the kinetic energy primal 2-form; and F e " HC e . A detailed description of these discrete variables and their staggering on the computational grid can be found in Appendix E, and a diagram of their staggering is in Figure 2. The resulting discrete evolution equations are In fact, by making alternative choices for F e , Q and Φ i (along with the operators discussed below) it is possible to recover a wide range of C grid schemes present in the literature (such as [18], [30] and [35]), see THESIS for more details). The operators D 2 , D 1 ,D 1 ,D 2 , I, J, R, W and H are defined in Appendices B and C (and can also be found in a general form in [29]). The novelty of the current scheme is a new definition of Q, such that the properties of total energy conservation, potential enstrophy conservation and steady geostrophic modes hold simultaneously. This is the subject of Section 4.

Linearized Scheme
As is well-known, the linearized version of a Hamiltonian system about a steady state can be found by evaluating the brackets at that state and using the quadratic approximation to the associated psuedo-energy as the Hamiltonian ( [23]). Following this procedure and letting the Coriolis force f be a constant, b i " 0 and assuming a background state of x " pH, 0q, we obtain for the brackets (where W " Q qv"1 is the linearized version of Q) and for the Hamiltonian, which has associated functional derivatives of The resulting evolution equations are Bu e Bt´f WHu e`gD1 Im i " 0 (48)

Properties of Scheme
This scheme has many important properties, including: 1. Mass and potential vorticity conservation: Both mass m i and mass-weighted potential vorticity m v q v are conserved in both a local (flux-form) and global (integral) sense.
2. No spurious vorticity production: By construction, D 2 D 1 " 0 and there is no spurious production of vorticity due to the gradient term in the wind equation.
3. Linear stability (pressure gradient force and Coriolis force conserve energy): This due to the fact that I, J and H are all symmetric positive-definite; D T 2 "´D 1 ;D 2 T " D 1 and W "´W T .
4. Steady geostrophic modes: By construction,´RD 2 " WD 2 (noting that W is the same for all members of this family), which gives steady geostrophic modes.
5. PV Compatibility: again by construction´RD 2 " WD 2 with Q qv"c Ñ cW, and therefore the potential vorticity equation is compatible with the diagnostic mass equation (a constant PV field remains constant). Note that this is same as the condition required for steady geostrophic modes.
6. Other conservation properties: see below for a discussion on total energy and potential enstrophy conservation. Table 1 shows a summary of the required properties in order for the resulting scheme to have all of the mimetic and conservation properties discussed above. Table 1: Summary of required operator properties for obtaining the desirable mimetic properties along with total energy and potential enstrophy conservation. For example I is a discrete Hodge star that maps from primal 2-forms to dual 0-forms, and must be symmetric positive. The only operator that merits additional explanation here is φ-it is used to construct mass at edges for use in determining the mass flux, and its transpose φ T is used for kinetic energy calculations. This ensures that the scheme conserves energy, see [29] or THESIS for more details.
Operator Properties Notes see text see text see text

Total Energy Conservation
Following S04, total energy will be conserved for any choice of H if the discrete brackets retain their antisymmetric character. This requires that D T 2 "´D 1 , and that Q "´Q T . The first condition is satisfied by construction of the discrete exterior derivative operators D 2 andD 1 . The second condition is satisfied only for certain choices of Q. One example is Q " 1 2 q e W`1 2 Wq e (as used in [18]), where q e is any function that, given the set of q v at primal vertices, computes a unique q e at primal edges (such as q e " 1 2 ř vPV Epeq q v ). Flexibility in the choice of q e allows a wide variety of stabilization methods such as CLUST or APVM ( [36] and [37]). Unfortunately, this choice does not conserve potential enstrophy.

Potential Enstrophy Conservation
Following S04, potential enstrophy is a Casimir and therefore will be conserved when tZ, Au " 0 (49) holds for any choice of functional A. Note that is the potential enstrophy where q v is potential vorticity primal 0-form, Note that q v " ηv mv , where m v " Rm i is the mass dual 2-form and η v " ζ v`fv "D 2 u e`fv is the absolute vorticity dual 2-form. Its functional derivatives are δZ δ x "˜´R Using the chain rule for functional derivatives, it suffices to show that equation (49) holds for A " ř i m i and A " ř e u e . Therefore equation (49) reduces to which must hold for any choice of q v . The first of these is again satisfied by construction for D 2 and D 1 . The second is much trickier, and is the main subject of section 4. One example is Q " q e W (as used in [18]), where q e " 1 2 ř vPV Epeq q v . Unfortunately, this choice does not conserve total energy. It would be possible to explore alternative definitions of Z, but these would lead to different, less natural stencils for q v .

Arakawa and Lamb 1981
In the case of a uniform square grid, the C scheme grid above reduces to the well-known Arakawa and Lamb 1981 total energy and potential enstrophy scheme (modified to prognose m i and u e if their choice of Q is used. Unfortunately, the definition of Q presented in AL81 works only for logically square, orthogonal grids. For more general, non-orthogonal polygonal grids, a new operator Q must be found. This is the subject of the next section.

Hollingsworth Instability
Since this is an extension of Arakawa and Lamb 1981 scheme, it seems extremely likely that the proposed scheme will suffer from the Hollingsworth instability, especially if applied in a height coordinate framework using a Lorenz staggering in the vertical (as discussed in [2] and [15]). It also seems likely that proposed scheme will avoid the Hollingsworth instability when used with an isentropic or Lagrangian vertical coordinate, or when a Charney-Phillips staggering is used in the vertical. If the instability is encountered, it would be simple to modify the stencil of the kinetic energy in a consistent manner (to preserve total energy conservation, by simply modifying the Hamiltonian itself), which has been shown to be sufficient to prevent the instability ( [15]). Therefore, the possible presence of the instability is not expected to prevent use of this scheme in a full 3D model.

Operator Q
The principal novelty of the new C grid scheme is the specification of a Q operator that simultaneously conserves total energy and potential enstrophy, and also supports PV compatibility. Previous work found choices for Q that conserved either total energy or potential enstrophy, but not both. The key lies in S04, showing that the AL81 approach could be extended to more general stencils (although retaining a logically square, orthogonal grid). This work takes the Salmon 2004 approach in a different direction, keeping the same stencil as AL81 but considering a general polygonal grid.

Definition of Q
Loosely following S04, define Q as where i is the primal grid cell covered by both e and e 1 . A diagram of this operator is shown in Figure  3. An equivalent alternative form for Q given in terms of the Poisson bracket that closely mimics the one found in S04 can be found in the appendix. It is easy to see that in the case of a logically square orthogonal grid, this approach reduces to the same stencil considered by AL81. At this point, the coefficients α e,e 1 ,v are undetermined.

Linear System for α
It remains to determine the coefficients α e,e 1 ,v in a manner such that the resulting operator Q conserves both total energy and potential enstrophy, and satisfies PV consistency.

Requirements introduced by energy conservation
Following S04, in order for Q to be energy conserving then Q "´Q T . In terms of the coefficients, this implies that α e,e 1 ,v "´α e 1 ,e,v , or in other words, they are anti-symmetric under an interchange of e and e 1 .

Requirements introduced by potential enstrophy conservation
From (53), in order for Q to conserve potential enstrophy´D 1 R q 2 v 2`Q D 1 q v " 0 must hold for any choice of q v . Expanding this out yields ÿ e 1 PECP peq¨ÿ vPEV Cpe,e 1 q α e,e 1 ,v q v‚ for every e, which must hold for any choice of q v . For a given edge e, the vertices in question are v P CV Epeq (shown in Figure 4) where CV Epeq " V Epi1q Y V Epi2q and pi1, i2q " CEpeq. Both the left and right hand side of these equations are a quadratic form in this set of vertices, and for this to hold for arbitrary q v the coefficients in these two quadratic forms must be equal. These coefficients are linear combinations of the α's, and therefore the equality of these quadratic forms implies a set of linear equations for the α's. Specifically, for each grid cell i with n e edges and n v vertices (note that n e " n v for a polygonal grid cell, but it is useful to keep distinct notation to ease exposition), there are n e nvpnv`1q 2 equations (coefficients in the quadratic forms) and n v nepne´1q 2 unknowns (the coefficients α e,e 1 ,v ). This is therefore an overdetermined system, and the coefficient will be found through a least squares procedure. At least some of the additional freedom will be used to split the equations into independent subset for each grid cell (see below), which makes implementation practical for operational grids. The equations come from equating the coefficients in the two quadratic forms: there are nvpnv`1q 2 independent vertex pairs, and n e edges. The unknowns are the coefficients α e,e 1 ,v that are associated with the grid cell: there are nepne´1q 2 independent unique edge pairs, and n v vertices. Note that this has already taken into account the fact that α e,e 1 ,v "´α e 1 ,e,v (hence the wording unique edge pair) which reduces the number of independent coefficients in half. Letting v and v 1 loop over the vertices in the cell (they are the unique members of V CpiqˆV Cpiq), the equations are given by A v,v " ÿ e 1 PEV Epv,e,iq α e,e 1 ,v t e 1 ,v sgnpe, e 1 q (56) where the sum for B v,v occurs only when v P V Epeq; and where e loops over each edge in i and EV Epv, e, iq " ECpiq X EV pvq´e; and sgnpe, e 1 q " 1 "´sgnpe 1 , eq (which ensures that the scheme is also energy conservative). A diagram of EV Epv, e, iq is provided in Figure  5. Note that coefficients in one cell are coupled with adjacent cells when v P V Epeq or v 1 P V Epeq; that is to say, the equations involve coefficients that are associated with other grid cells. On a non-uniform mesh, this means that the entire set of coefficients must be solved for at the same time.
The solution procedure outlined above gives a large matrix system where each row in A represents an equation obtained by equating coefficients in the quadratic forms, and α is the vector of unknown coefficients. This system can be solved (via a least-squares approach) to yield a set of coefficients α such that Q conserves potential enstrophy. This procedure is essentially identical to the one employed in S04; when applied to a uniform square grid it reproduces AL81 and produces a total energy and potential enstrophy conserving scheme on a uniform hexagonal grid (not shown, verified numerically). In addition, the coefficients only have to be computed once, and then stored for later use. Unfortunately, the system that results from this procedure is impractical to solve for realistic non-uniform meshes: it is too large and ill-conditioned. For example, on an icosahedral-hexagonal mesh with O(1 million) grid cells, there will be O(90 million) coupled coefficients that need to be solved for.

Practical Solution
Instead, following [31], the coefficients can be uncoupled by defining when v P V Epeq or v 1 P V Epeq, where C "´1{6. On all meshes tested (including uniform square and uniform grid) there are enough degrees of freedom to do this, and the least-squares problem has a unique, exact solution. This has enabled the solution of the system for cubed-sphere meshes with up to 884736 grid cells and icosahedral-hexagonal meshes with up to 655363 grid cells in a few hours using an unoptimized, serial algorithm on a laptop computer. Furthermore, the uncoupled nature of the problem (one small independent least-squares problem per grid cell) would facilitate easy parallelism if needed for larger meshes (and again, the coefficients only need to be computed once). Figure 5: A diagram of the stencil EV Epv, e, iq " ECpiq X EV pvq´e. Consider the set pv, e, iq denoted in green: then EV Eppv, e, iq are the red edges. Now consider the set pv, e, iq denoted in blue: then EV Eppv, e, iq is the brown edge.

PV Compatibility
The astute reader will note that nothing has been said yet about enforcing PV compatibility (Q qv"c " cW. It was originally believed that PV compatibility would have to added as additional equations in the matrixvector system. However, it was found that enforcing potential enstrophy conservation (even using the cell split form) was sufficient to ensure that Q was PV compatible. This corresponds with the results of S04 ( [20]), who did not explicitly add PV compatibility, yet all of his schemes had this property. The reasons behind this result are not yet understood. If PV compatibility had to be added explicitly, it would simply mean that ÿ vPV Cpiq for every edge pair pe, e 1 q; which could be easily added to the independent system of equations solved in each grid cell.

Z Grid Scheme
Unlike the C grid scheme, the Z grid scheme starts with Nambu brackets rather than Poisson brackets. This greatly simplifies the derivation, since only the triply anti-symmetric nature of the brackets must be retained to ensure total energy and potential enstrophy conservation: there is no consideration of Casimirs. Start by defining a set of collocated discrete variables which are pointwise values of h, ζ and µ at primal grid centers. More details about the grid, discrete operators and discrete variables can be found in Appendices A,D and F.

Functional Derivatives
The functional derivative of a general functional F with respect to discrete variable x i is then defined as where A i is the area of primal grid cell i. The diagnostic variables Φ i , χ i , ψ i and q i are defined through the functional derivatives of the discrete Hamiltonian H and discrete Potential Enstrophy Z as: At this point the discrete Hamiltonian H and discrete Potential Enstrophy Z are left unspecified.

Discrete Nambu Brackets
Following [22], the general discretization starts from the Nambu brackets (33), (34) and (35) for the shallow water equations in vorticity-divergence form. As long as these brackets retain their triply anti-symmetric structure when discretized, total energy and potential enstrophy will be automatically conserved for any definition of the total energy and potential enstrophy (with one caveat explained below). In addition, the bracket structure ensures that this conservation is local as well as global. That is, the evolution of a conserved quantity can be written in flux-form for each grid cell, where cancellation of fluxes between adjacent cells leads to the global integral being invariant. This is in contrast to a method that conserves the global integral, but cannot be written in flux-form for each grid cell. In what follows below, we will consider only the case where Z is the potential enstrophy, although this approach could be easily generalized to arbitrary Casimirs (see [21] for an example of this on a uniform square grid).

Jacobian Brackets
Loosely following S07, the tF, H, Zu ζζζ bracket can be discretized as Note that this bracket is triply anti-symmetric (due to the cyclic permutation), as required. The tF, H, Zu µµζ bracket can be similarly discretized as tF, H, Zu µµζ " This bracket is only doubly anti-symmetric (in H and F due to the anti-symmetry of J), but it will conserve Z as well provided that δZ δµi " 0 (since JpA, Bq " 0 when either A " 0 or B " 0). These brackets are essentially those encountered when discretizing the Arakawa Jacobian, as detailed in [21].

Mixed Bracket
The mixed bracket is trickier since it contains an apparent singularity ( 1 ∇q q. On closer inspection, in the continuous case this singularity cancels out when combined with the functional derivative of the potential enstrophy. This is the caveat mentioned above-the discrete mixed bracket must be constructed such that the apparent singularity cancels out with the discrete functional derivative of the potential enstrophy. With this in mind, the general form of the discrete mixed bracket is chosen as: where, from before, q i " δZ δζi . This bracket is triply anti-symmetric (again due to the cyclic permutation), and the apparent singularity will cancel if Z is chosen with care.

Conservation
Since the tF, H, Zu ζζζ and tF, H, Zu µζh brackets are triply anti-symmetric, and the tF, H, Zu ζµµ bracket is doubly anti-symmetric, both total energy and potential enstrophy will be conserved for any choice of H and Z; provided that the caveats mentioned above are obeyed. Those are: 1. δZ δµi " 0 (ensures that the tF, H, Zu ζµµ bracket conserves potential enstrophy) 2. Z chosen such that the apparent singularity (D 1pZh q D1qi term + cycpF, H, Zq terms) in the tF, H, Zu µζh bracket cancels out These are fairly minimal requirements, and many reasonable choices for Z satisfy them.

Discrete Hamiltonian and Helmholtz Decomposition
The Hamiltonian H can be split into three parts: H F D , H J and H P E , where the first two are the kinetic energy due to flux-divergence terms and Jacobian terms, and the last is the potential energy. In the continuous system we have where These can be discretized as

Helmholtz Decompositions and Bernoulli Function
By taking variations of H we obtain After a lot of algebra, these can be grouped (half of each term involving δh i goes to Φ i and half to µ i /ζ i ) to obtain δH "´χ i δµ i`´ψi δζ i`Φi δh i where (using the definition of functional derivative) The latter two equations (85 and 86) are the discrete version of the Helmholtz decomposition, and form a pair of non-singular elliptic equations. They can be combined into a single equation as where, for example, FDχ i " 1 Ai hv qψ e . Note that (without the 1 Ai factors) FD is symmetric and JA is anti-symmetric, which means that A "´A T (ie A itself is skewsymmetric). Also note that when h i " H is a constant (and therefore h e " H), they reduce to where L " 1 Ai D 2 le deD 1 , which is the correct linearization behaviour.

Discrete Potential Enstrophy
A natural definition of the discrete potential enstrophy is where η i " ζ i`fi . Taking variations of this yields Then the natural definition for q i " ηi hi works, and the above simplifies to By plugging these back into the tF, H, Zu µζh bracket, it is seen that this choice of Z also ensures that the singularity cancels.

Independence between choices for H/Z and Nambu Brackets
As noted before, the mimetic and conservation properties of the discrete scheme are completely independent of the choice of discrete Hamiltonian H, provided the Hamiltonian is positive definite and produces invertible elliptic equations for the Helmholtz decomposition. If the resulting elliptic equations were singular, then the scheme would have a computational mode (as discussed in [22]). Additionally, the discrete Helmholtz decomposition should also simplify to a pair of uncoupled Poisson problems when linearized. The mimetic and conservation properties are also independent of the specific choice of Z, provided that the singularity in the mixed bracket cancels and Z δ " 0. The given choices of H and Z were selected to have these properties, and also correspond with those in S07 for the special cases of a uniform planar square grid and an orthogonal polygonal planar grid with a triangular dual.

Discrete Evolution Equations
By setting F " ph i , ζ i , µ i q in turn, the following evolution equations are obtained: where L is the Laplacian, FD is the Flux-Divergence and J is the Jacobian. Note that these operators on an icosahedral hexagonal-pentagonal grid are the same as those from [12]. The only difference is in the arguments (q i instead of η i , and different definitions for χ i and ψ i .)

Laplacian and Flux-Div Operators
The Laplacian and Flux-Divergence operators (which come from the mixed bracket) can be written as where α e " ř iPCEpeq αi 2 .

Jacobian Operators
The Jacobian operators (which come from the Jacobian brackets) can be written as Note that on a polygonal grid with a purely triangular dual (including the important case of an icosahedral grid), J δ " J ζ .

Linearized Version
Under the assumption of linear variations around a state of rest (h i " H, ζ i " µ i " 0, q i " f H ) on a f-plane, this scheme reduces to: where the Helmholtz equations given by (88) and (89) have been used to simplify the scheme (to the point that it no longer requires solving any elliptic equations). In the case of a uniform square grid (uniform hexagonal grid) this scheme is identical to the one studied in [17] ( [9]), and it shares the same excellent linear wave properties found for those schemes.

Relation to Salmon Schemes
For the cases of a uniform planar square grid and a general orthogonal planar polygonal grid with triangular dual, the general discretization scheme presented above reduces to the schemes given in S07. However, this discretization scheme is more general, and it also makes specific choices for the total energy H and potential enstrophy Z when using a general polygonal grid.

Properties of Scheme
The discrete scheme as outlined above posses the following (among others) key properties: 1. Linear stability (Coriolis and pressure gradient forces conserve energy): Provided that L " L T (which is satisfied for the L given above, and the majority of discrete Laplacians), the scheme will conserve energy in the linear case.
2. No spurious vorticity production: By construction, the pressure gradient term does not produce spurious vorticity since the curl is taken in the continuous system, prior to discretization.
3. Conservation: By construction, this scheme conserves mass, potential vorticity, total energy and potential enstrophy in both a local (flux-form) sense and global (integral) sense.
4. PV compatibility and consistency: By inspection, the mass-weighted potential vorticity equation is a flux-form equation that ensures both local and global conservation of mass-weighted potential vorticity. In addition, an initially uniform potential vorticity field will remain uniform. This rests on the fact that J ζ pq i , ψ i q " 0 and FDpq i , χ i q " cLχ i when q i " c is constant.
5. Steady geostrophic modes: Since the same divergence µ i appears in both the linearized vorticity and continuity equations, the scheme posses steady geostrophic modes.
6. Linear properties (dispersion relations, computational modes): As expected, the scheme possesses the same linear mode properties on uniform planar grids as those presented in [17] and [9]; and it does not have any computational modes. More details of the linear mode properties of the scheme on both uniform planar and quasi-uniform spherical grids can be found in a forthcoming paper [7].
7. Accuracy: Unfortunately, as shown in [14], the Jacobian operator as given is inconsistent on general grids. Even more unfortunately, the fix proposed in that paper breaks key properties of the Jacobian necessary to retain total energy and potential enstrophy conservation. Surprisingly, as shown in [8], the inconsistency of the Jacobian operator does not appear to cause issues in the test cases that were run. More details on possible fixes to the accuracy issue are discussed in [8].

Conclusions
This paper presents an extension of AL81 to arbitrary non-orthogonal (spherical) polygonal grids in a manner that preserves almost all of the desirable properties of that scheme (including both total energy and potential enstrophy conservation) through a new Q operator. Unfortunately, on non-quadrilateral grids such as the icosahedral grid there will be extra branches of the dispersion relationship due to a mismatch in the number of degrees of freedom in the wind and mass fields inherent to the C grid approach. Switching from a C grid type staggering (to an A grid staggering, for example) is undesirable for many reasons, foremost among them being the natural association of physical variables with geometric entities in a staggered grid as suggested by exterior calculus and differential geometry (see [33] and [3]). Fortunately, other than these extra mode branches on the icosahedral grid the proposed C grid scheme does not posses any additional computational modes. Furthermore, extensive testing has thus far been unable to show negative impacts from this extra mode branch, especially when running full-physics simulations with realistic topography and initial conditions (John Thuburn and Bill Skamarock, personal communication). This work has also presented an extension of the total energy and potential enstrophy conserving Z grid scheme in S07 from planar grids to arbitrary orthogonal (spherical) polygonal grids, using the same toolkit of Nambu brackets and Hamiltonian methods. The restriction to orthogonal grids (geodesic grids are the only orthogonal quasi-uniform spherical grid the author is aware of) rather than more general nonorthogonal grids is a drawback. However, the major motivations for using a cubed-sphere grid are the ability to properly balance degrees of freedom when using a staggered C grid methods (and therefore avoid spurious branches of the dispersion relationship), a tensor-product grid structure for spectral or finite element type methods (which ensures a diagonal mass matrix for spectral element methods and eases implementation of finite element methods) and higher-order finite volume methods (enabling easy dimension splitting), and an underlying piecewise continuous coordinate system for higher-order finite volume methods (allowing extended stencils). None of these considerations apply to a Z grid method, so the restriction to icosahedral grids is not anticipated to be a significant hurdle.
A detailed comparison of the two schemes, including an analysis of the accuracy of the operators used and results from a variety of test cases can be found in second part of this series [8]. In addition, an analysis of the linear mode properties of these two schemes on various quasi-uniform grids is undertaken in the third part of this paper series [7].

Code Availability
The schemes described in this manuscript have been implemented in a Python/Fortran mixed language code, and are freely available at https://bitbucket.org/chris_eldred/phd_thesis under a GNU Lesser General Public License Version 3.

A Discrete Grid
The schemes described above are designed to work on arbitrary (spherical) polygonal grids along with an associated dual grid. In the case of the C grid scheme, the grid can be either orthogonal or non-orthgonal, while the Z grid scheme is restricted to orthogonal grids. A description of the this grid framework is given in what follows.

A.1 General Non-Orthogonal Polygonal Grid
Consider a (primal) conformal grid constructed of polygons (or spherical polygons). A dual grid is constructed such that there is a unique one to one relationship between elements of the primal grid and element of the dual grid: primal grid cells are associated with dual grid vertices, primal grid edges are associated with dual grid edges and primal grid vertices are associated with dual grid cells. This grid configuration covers the majority of grids that are used in current and upcoming atmospheric dynamical cores, including cubed-sphere Figure 6: The geometric quantities on a planar grid. Primal grid edge lengths are denoted as de, dual grid edge lengths are denoted as le, the area associated with an edge by A e , the overlap between primal grid cell i and edge e by A ie and the overlap between dual grid cell v and edge e by A iv . Note that the same definitions can be used on a spherical grid, provided the appropriate measures are used (such as geodesic lengths for distances, and spherical polygonal areas for areas). See [35] for more details. and icosahedral grids (both hexagonal-pentagonal and triangular variants). Once the dual grid vertices have been placed, there are several important geometric quantities that are needed in order to construct the discrete operators (shown graphically in Figure 6). Specifically, we need the primal cell area A i , the dual cell area A v , the distance between primal grid centers le, the distance between dual grid centers de and the overlap areas A iv and A ie . On a planar grid, these are easily defined using the standard Euclidean metric and formulas. On a spherical grid, distances must be calculated using geodesic arcs; and areas are calculated by subdividing into spherical triangles as needed and then applying the relevant spherical area formulas. See the discussion in [35] for more details.

B Discrete Exterior Calculus Operators
Following [29], a set of discrete exterior derivative operators can be defined as: where n e,i is an indicator that is 1 when e is oriented out of a primal grid cell and -1 when e is oriented into a primal grid cell, and t e,v is an indicator that is 1 when e is oriented into a dual grid cell and -1 when e is oriented out of a dual grid cell. Note that by construction, these satisfy D 2 D 1 " 0,D 2D1 " 0, D T 2 "´D 1 andD 2 T " D 1 for arbitrary polygonal grids.

C Specific Choices for Various C Grid Operators
In order to close the C grid scheme presented in Section 3, specific choices must be made for I, J, H, R, φ and W. The ones used here (and in [18] and [30]) are: where H O is used on orthogonal grids such as the icosahedral grid, H N O is used on non-orthogonal grids such as the cubed-sphere grid (the details of the construction of this operator, including the stencil Speq and the weights H e,e 1 , can be found in [30]) and the weights W e,e 1 are chosen such that W "´W T and RD 2 "D 2 W (the details for this operator can be found in [31]). On an orthogonal grid, I, J, H correspond to the choice of a Voronoi hodge star from discrete exterior calculus.

D Specific Choices for Various Z Grid Operators
For the Z grid scheme, the following operators are needed: Note that JpA, Bq is anti-symmetric (JpA, Bq "´JpB, Aq) and satisfies JpA, 0q " JpB, 0q " JpA, Aq " 0. In addition, two different interpolations (from cell centers to vertices and to edges, respectively) are defined: where C is a constant given by 1 n , where n is the size of CV pvq (equal to 4 for quadrilateral dual grid cells and 3 for triangular dual grid cells). Table 2 gives the discrete variables used in the C grid scheme, their type (which indicates the staggering on the grid), and their diagnostic equation (where applicable). For the type, the first designator indicates the form type (primal or dual) and the second designator indicates the form degree (0,1 or 2). For example, C e is a dual 1-form. The only exceptions to this are the edge mass m e , which is used in constructing the dual mass flux C e ; and the edge PV q e , which is used in constructing Q for the variants that conserve only total energy or potential enstrophy. These quantities are not really physical, but instead are just used computationally to construct other, physical quantities or operators.

E Discrete Variables (C Grid Scheme)
Kinetic Energy  Table 3 gives the discrete variables used in the Z grid scheme and their type (either prognostic or diagnostic).