Development and technical paper 22 Nov 2018
Development and technical paper  22 Nov 2018
Symmetric equations on the surface of a sphere as used by model GISS:IB
 ^{1}NASA Goddard Institute for Space Studies, 2880 Broadway, New York, NY 10025, USA
 ^{2}Center for Climate Research, Columbia University, New York, NY 10025, USA
 ^{1}NASA Goddard Institute for Space Studies, 2880 Broadway, New York, NY 10025, USA
 ^{2}Center for Climate Research, Columbia University, New York, NY 10025, USA
Correspondence: Gary L. Russell (gary.l.russell@nasa.gov)
Hide author detailsCorrespondence: Gary L. Russell (gary.l.russell@nasa.gov)
Standard vector calculus formulas of Cartesian three space are projected onto the surface of a sphere. This produces symmetric equations with three nonindependent horizontal velocity components. Each orthogonal axis has a velocity component that rotates around its axis (eastward velocity rotates around the north–south axis) and a specific angular momentum component that is the product of the velocity component multiplied by the cosine of axis' latitude. Angular momentum components align with the fixed axes and simplify several formulas, whereas the rotating velocity components are not orthogonal and vary with location. Three symmetric coordinates allow vector resolution and calculus operations continuously over the whole spherical surface, which is not possible with only two coordinates. The symmetric equations are applied to onelayer shallow water models on cubedsphere and icosahedral grids, the latter being computationally simple and applicable to an ocean domain. Model results are presented for three different initial conditions and five different resolutions.
According to the “hairy ball theorem” of Poincaré (proved by Brouwer), every continuous horizontal vector field on the surface of a sphere must have a 0; a unit vector field must have a discontinuity. A differentiable coordinate on the surface of a sphere, for example, latitude or longitude, will have a gradient unit vector that is tangent to the sphere; such a coordinate will also have a discontinuity. A horizontal vector component has a magnitude and a direction based on an underlying coordinate (e.g., eastward velocity and longitude). The coordinate will have a discontinuity somewhere, so if the component is to be continuous over the whole sphere, its magnitude must be 0 where the coordinate is discontinuous. Scalar quantities have no associated direction and may be continuous over the whole sphere. Spatial derivatives acting on scalar quantities may use local coordinates, and coordinates of a point need not be continuous with respect to those of an adjacent point. Spatial derivatives acting on vector components require the components to be continuous. Coordinate discontinuities occur at the poles on a latitude–longitude (lat–long) grid or at the face edges on a cubedsphere grid if coordinates are switched.
Many numerical schemes for solving the fluid dynamic equations on the surface of a sphere use two independent horizontal velocity components aligned with underlying coordinates. On a lat–long grid, polar singularities occur as well as other problems discussed in the introduction of Heikes and Randall (1995). Components are orthogonal on a lat–long grid, but if they are not orthogonal on another grid, then greater obtuseness of the components' angle decreases stability and precision of the results. In addition, formulas involving the two coordinates are often not symmetric; this is certainly the case for standard lat–long schemes (Williamson et al., 1992) or spectral harmonic schemes (Temperton, 1991; Swarztrauber, 1993).
To use the lat–long grid, but to avoid its polar deficiencies, the HYbrid Coordinate Ocean Model (HYCOM) (Sun and Bleck, 2005) uses a Mercator lat–long grid south of 58^{∘} N, and a different polar projection north of it, introducing the necessity to match variables along the boundary. The Yin–Yang grid (like the cover of a baseball) (Qaddouri, 2011) uses two orthogonal lat–long grids, one with a north–south axis and another with an equatorial axis; meeting of the domains is complicated, but the results are reasonable. Weller et al. (2012) tested a “skipped lat–long” grid that uses fewer cells in longitude poleward of 66^{∘} in a onelayer model. Other researchers have explored cubedsphere grids: Adcroft et al. (2004) and Putman and Lin (2007). The advantages of using a cubedsphere grid over a lat–long grid are that the two polar singularities on the lat–long grid are replaced by eight illbehaved corners on the cubedsphere grid and the extreme aspect ratio of grid cells near the poles on a lat–long grid is eliminated as well as the polar filter used to increase the time step. A more thorough discussion of advantages and disadvantages of a cubedsphere grid over a lat–long grid occurs in Putman and Lin (2007). A problem with the abovementioned grids is that the treatment and behavior at grid edges differs significantly from that away from the edges. Choices that must be made in the pursuit of consistency have the potential for inducing edge errors. This will not be the case for the approach presented here.
Recent research on onelayer models has been directed at icosahedral grids: Heikes and Randall (1995), Stuhne and Peltier (1999), Läuter et al. (2008), Lee and MacDonald (2009), Ringler et al. (2010), Weller et al. (2012), and others. The grid cells are usually pentagons and slightly irregular hexagons. Although irregularities are distributed all over the sphere, the hope is that icosahedral grid errors are less concentrated than edge errors of other grids and also less severe.
If twocomponent velocity is the prognostic transport variable that is advected in flux form, then spatial derivatives of vector components will cause discontinuities to occur. This is a principal reason why researchers developed forms of the shallow water equations wherein scalar quantities such as potential vorticity, specific kinetic energy, and divergence are continuous everywhere. Computations are performed on local spherical coordinates or on the local tangent plane after which the horizontal velocity components are resurrected or time integrated by manipulating spatial derivatives of scalar quantities; spatial derivatives of vector components are not needed. Such forms include vectorinvariant (Ringler et al., 2010), vorticity–divergence (Williamson et al., 1992), or stream function and velocity potential (Masuda and Ohnishi, 1986), but the equations and programming can be complex.
The approach here uses three symmetric coordinates on the surface of a sphere to represent twodimensional flow. When one coordinate reaches a singularity, it is ignored and the other two coordinates become perpendicular on the spherical surface. Symmetric equations are used to develop onelayer shallow water equation models: one for a gnomonic cubedsphere grid (CSK), one for an icosahedral B grid with momentum defined at the primary grid cell corners (IB), and one for an icosahedral grid with a Voronoi tessellation (IK).
Several formulas of the symmetric equations are simplified by using relative specific angular momentum on the unit sphere, A, as opposed to using the three velocity components. A is continuous everywhere; each component converges to 0 at its respective poles. A has been used by Ringler et al. (2010) and other researchers under the name u^{⟂}, but it was not recognized as the specific angular momentum vector on the surface of a sphere. The north–south axis component of A is identical to eastward ucos ϕ in the spectral model of Swarztrauber (1993) which also uses northward vcos ϕ. Both components are 0 and are continuous at the poles, but the other components of angular momentum are absent. Swarztrauber (1996) presents a spectral transform threedimensional Cartesian method to solve the shallow water equations on the sphere written in vorticity–divergence form. Putting aside the spectral transform method and vorticity–divergence form, there are some similarities (he again uses ucos ϕ) and major differences. His equations exist in R^{3} and are later restricted to the spherical surface, whereas symmetric equations are compressed to the surface from the beginning, and his equations use velocity instead of specific angular momentum.
The shallow water equations based on A are simpler than those using velocity or those using vectorinvariant methods. Components of A are symmetric; polar problems are absent; and the metric term disappears from the momentum equation that exists when using eastward and northward velocity. Conservation of A by advection in flux form is precise without time truncation errors. The symmetric equations are performed on the spherical surface without relying on tangent plane computations. Horizontal velocity and related A must be tangent to the spherical surface that allows only two degrees of freedom. Thus, three momentum components are not independent; there is a required alignment. If some process disturbs this alignment (e.g., advection or nonhorizontal acceleration), causing momentum to no longer be horizontal, then a simple algorithm brings the three components back into alignment. This problem was recognized by Coté (1988), who added a “Lagrange multiplier” term to the velocity equation of Cartesian R^{3} that was restricted to the spherical surface.
Computers require grid representations of differential equations; this causes numerical errors that relate to grid imprinting, mass variations, time integration, etc. Grid imprinting is easily recognized when integrating the solid body rotation Test Case 2 of Williamson et al. (1992), which lacks bottom topography. Errors due to grid imprinting should and do decrease with finer resolution.
Mass variations cause advection errors in numerical models and cause gridmatched alternating patterns. Mass is usually conserved by programming advection to use flux form, but when mass is needed at different locations, it is specific mass or concentration that is interpolated. Tracers that follow mass advection include linear momentum and velocity, angular momentum and specific angular momentum, kinetic energy and specific kinetic energy, or absolute vorticity and potential vorticity. Russell and Lerner (1981) investigated mass variations and stated, “Tracer concentration is defined relative to air mass and not relative to space. In fact, with nonuniform mass, second and fourthorder schemes become firstorder schemes …” The mass coordinate at a point, in a spatial onedimensional model, is the number of kilograms between the origin and the point. Secondorder, fourthorder, and spectral schemes, in one dimension, usually compute their polynomial or sine wave coefficients from equally distributed points in space, but the points are not equally distributed in the mass coordinate, and consequently the derived coefficients are erroneous.
If tracer concentration is a linear function of mass over several grid cells in one dimension, then “the linear upstream scheme” of Leer (1977) will perform the advection over those cells perfectly, even with arbitrary mass variations. Similarly, if concentration is a quadratic function of mass in one dimension, then “the quadratic upstream scheme” of Prather (1986) performs advection perfectly. Each of these schemes uses mean tracer values and prognostic tracer gradients (and secondorder moments) inside each grid cell; these schemes are less sensitive to mass variations than are nonupstream schemes. The onelayer icosahedral model IB, discussed in detail later in this paper, uses a combination of linear upstream and secondorder schemes for momentum advection but does not carry prognostic gradients. According to Weller et al. (2012), “an upwindbased interpolation of the potential vorticity controls the computational Rossby modes” in some onelayer models.
To be applied to the Earth, models should be tested with mass variations comparable to those on Earth. In some mountainous regions, surface pressure gradients and mass variations increase with finer resolution and so do their errors. Test Case 5 of Williamson et al. (1992) is insufficient to judge model quality considering the topography variations that occur on Earth. Surface pressure gradients in Williamson do not vary with resolution and are 2 orders of magnitude smaller than the gradient at some locations on Earth, particularly the Andes. When a onelayer model is tested with various resolutions, it should include tests where such increased gradients exist. If a onelayer model requires special filtering to produce good results with an Earth topography simulation, then that filtering should be used when the model is applied to the less rigorous tests of Williamson.
Section 2 explains symmetric mathematics, including coordinates and variables, symmetric calculus operators, and the differential solution to the shallow water equations in terms of symmetric coordinates. Section 3 presents the discrete implementation of alignment, the pressure gradient force, advection, the Coriolis force, and other aspects of icosahedral models including grid arrangements and time steps. Section 4 applies three test cases to lat–long, cubedsphere, and icosahedral onelayer models with various resolutions. Section 5 contains discussion and conclusions.
Unlike computer languages, division has lower precedence than does multiplication in this paper. Vector quantities are indicated by bold capital letters; when displayed by three coordinates, e.g., $\mathbf{\text{A}}=(a,b,c)$, they are the Cartesian coordinates of R^{3} on the fixed Earth and not a local coordinate system.
A sphere of radius 1 is centered at the origin in threedimensional Cartesian space with axis unit vectors $\mathbf{\text{X}}=(\mathrm{1},\mathrm{0},\mathrm{0})$, $\mathbf{\text{Y}}=(\mathrm{0},\mathrm{1},\mathrm{0})$, and $\mathbf{\text{Z}}=(\mathrm{0},\mathrm{0},\mathrm{1})$. A concentric sphere represents the fixed Earth with its radius and its geometric labels such as latitude, longitude, Equator, and poles. The Z axis is aligned with the Earth's north–south axis, but this is only necessary for simplifying the sine of latitude for the Coriolis force. Horizontal velocity is resolved by the three components, one for each of the mutually orthogonal axes. One component is eastward velocity and its coordinate is longitude, which rotates around the Z axis. The other two components rotate around the two equatorial axes. Symmetric versions of the shallow water equations and vector calculus formulas are presented, in which each coordinate has symmetric representation in the equations and formulas. Each coordinate is defined continuously except at its two poles. But, as a coordinate approaches one of its poles, its numerical importance in the formulas decreases, having no importance at the pole, and the other two coordinates become more nearly perpendicular. Consequently, vector resolution and calculus operations are continuous over the whole sphere including the poles.
2.1 Coordinates and variables
Cartesian coordinates on the surface of the unit sphere are labeled $\mathbf{\text{P}}=(p,q,r)$, where ${p}^{\mathrm{2}}+{q}^{\mathrm{2}}+{r}^{\mathrm{2}}=\mathrm{1}$; this label is also used for points on the Earth's surface, but they mean the projection onto the unit sphere. Note that r and $\sqrt{{p}^{\mathrm{2}}+{q}^{\mathrm{2}}}$ are the sine and cosine of the Earth's latitude.
Three horizontal velocity components, u, v, and w (m s^{−1}), defined everywhere on the Earth's surface, except at their poles, rotate around each respective axis. Velocity unit vectors of the three symmetric components and the northward velocity component n, at point P, are
S (m s^{−1}) is horizontal velocity at point P and is equal to u of Cartesian R^{3} in Ringler et al. (2010). Note that $u=\mathbf{\text{S}}\cdot \mathbf{\text{U}}$, $v=\mathbf{\text{S}}\cdot \mathbf{\text{V}}$, and $w=\mathbf{\text{S}}\cdot \mathbf{\text{W}}$. Although there are three components, there are only two degrees of freedom, since $\mathbf{\text{S}}\cdot \mathbf{\text{P}}=\mathrm{0}$.
Relative specific angular momentum on the unit sphere, A (m s^{−1}), uses the variables a, b, and c:
The component c pointing toward the Z axis is equal to the eastward velocity multiplied by the cosine of the Earth's latitude. If A and S are properly aligned, i.e., horizontal or tangent to the spherical surface at P, then they have the relationships $\mathbf{\text{S}}\cdot \mathbf{\text{P}}=\mathbf{\text{A}}\cdot \mathbf{\text{P}}=\mathrm{0}$, $\mathbf{\text{A}}=\mathbf{\text{P}}\times \mathbf{\text{S}}$, and $\mathbf{\text{S}}=\mathbf{\text{A}}\times \mathbf{\text{P}}$. Thus, on the surface of the sphere, A is at right angle to S and both are tangent to the surface. The components of A are mutually orthogonal, being aligned with the X, Y, and Z axes that are fixed with respect to the Earth, whereas the components of horizontal velocity aligned with the unit vectors U, V, and W are not orthogonal and change with location. Also, a, b, and c are continuous everywhere, whereas u, v, and w are discontinuous at their respective poles. S, written in terms of u, v, and w, can be simplified using the components of A and is most readily derived from A×P.
The northern velocity component rotated 90^{∘} from W is $n=\mathbf{\text{S}}\cdot \mathbf{\text{N}}=(qapb)/\sqrt{{p}^{\mathrm{2}}+{q}^{\mathrm{2}}}$. The velocity component rotated 90^{∘} from U is $\mathbf{\text{S}}\cdot \mathbf{\text{P}}\times \mathbf{\text{U}}=(rbqc)/\sqrt{{q}^{\mathrm{2}}+{r}^{\mathrm{2}}}$ and that rotated from V is $\mathbf{\text{S}}\cdot \mathbf{\text{P}}\times \mathbf{\text{V}}=(pcra)/\sqrt{{r}^{\mathrm{2}}+{p}^{\mathrm{2}}}$, which are used for the Coriolis force. Velocity squared is $\mathbf{\text{S}}\cdot \mathbf{\text{S}}=\mathbf{\text{A}}\cdot \mathbf{\text{A}}={a}^{\mathrm{2}}+{b}^{\mathrm{2}}+{c}^{\mathrm{2}}$.
Spherical angular rotation coordinates, measured in radians, that rotate around the X, Y, and Z axes are μ, ν, and λ, respectively. Angular rotation coordinates measured from pole to pole for each axis are δ, ϵ, and ϕ, respectively. λ and ϕ are the Earth's longitude and latitude. A point on the sphere can be designated by any of four different coordinate systems:
At any point, the gradients of μ, ν, λ, and ϕ on the spherical surface are parallel to the unit vectors U, V, W, and N. Partial derivatives of the angular rotation coordinates with respect to one another are needed to derive new and old forms of various terms. A change in Δϕ causes a change in Δμcos δ in the ratio of U⋅N. In the limit, $\partial \mathit{\mu}/\partial \mathit{\varphi}=\mathbf{\text{U}}\cdot \mathbf{\text{N}}/\mathrm{cos}\mathit{\delta}$. A few useful derivatives are
Prognostic variables for the shallow water equations on the sphere are the height field above the surface topography, h (m), and A. Because density (kg m^{−3}) is uniform and is set to 1, h and mass per unit area are used interchangeably. The surface topography, h_{S} (m), is specified. R (m), the Earth's radius, g (m s^{−2}), the downward vertical acceleration due to gravity, and Ω (s^{−1}), the Earth's angular rotation rate, are assumed to be uniform. The field top geopotential Φ (m^{2} s^{−2}) is g(h+h_{S}).
The new symmetric equations to be presented here are applicable to many grid arrangements; one is the gnomonic cubedsphere grid. A symmetric tessellation of the surface of a cube [$\mathrm{1}:\mathrm{1}$, $\mathrm{1}:\mathrm{1}$, $\mathrm{1}:\mathrm{1}$] by line segments of length 2 parallel to an edge projects a rectangular grid onto great circle arcs on the surface of the sphere. Faces of the cube are numbered 1 to 6: z=1, $y=\mathrm{1}$, x=1, $z=\mathrm{1}$, y=1, and $x=\mathrm{1}$; see Fig. 1. On the sphere, grid cells are shaped like parallelograms and at least one of the unit vectors, U, V, or W, will be perpendicular to any grid edge. It was this property that led to the exploration of symmetric equations. At corners of the cube, the unit vectors are separated by 60 or 120^{∘}. A proprietary onelayer shallow water equation model, labeled CSK, may be published later, but some results are presented in Sect. 4.
Two new shallow water equation models were developed on icosahedral grids. The centers of primary grid cells are the vertices of a triangular lattice covering the sphere; the cells themselves are pentagons and irregular hexagons. The Bgrid icosahedral model, labeled IB, is discussed in detail in this paper. Its momentum cells (described as the dual mesh by Ringler et al., 2010) are spherical triangles whose vertices are three primary cell centers and whose momentum centers are primary cell corners. The second icosahedral model, IK, uses proprietary techniques similar to CSK. All three models use the symmetric equations, but IB and IK perform less quantity averaging than CSK does with its unusual grid cell shape. Figure 2 shows a fourtriangle wedge of an icosahedral grid.
2.2 The symmetric ∇_{S} operator on the surface of a sphere
The symmetric del operator on the surface of a sphere ∇_{S}, or simply ∇ in this paper, is
which is equivalent to the common twodimensional del operator on the surface
of a sphere,
${\mathrm{\nabla}}_{\text{R}}=(\mathbf{\text{W}}\partial /\partial \mathit{\lambda}+\mathbf{\text{N}}\mathrm{cos}\mathit{\varphi}\partial /\partial \mathit{\varphi})/R\mathrm{cos}\mathit{\varphi}$ (Williamson et al., 1992,
Eq. 3). ∇_{S} or ∇_{R} can be applied to three space, and be equivalent
to the threedimensional Cartesian operator ${\mathrm{\nabla}}_{C}=(\partial /\partial x,\partial /\partial y,\partial /\partial z)$, by restricting quantities to
which they are applied to be radially constant or by adding the term
$\mathbf{\text{P}}\partial /\partial \mathit{\rho}$ to the del operator, where ρ (m) is
the radial coordinate. Computations of $\partial Q/\partial \mathit{\lambda}$ (or
$\partial Q/\partial \mathit{\varphi}$), an arbitrary scalar Q, are identical in
formulas that use ∇_{S} or ∇_{R}, but the factors that multiply
those partial derivatives in the formulas are different.
The gradient of a scalar h on the Earth's surface is
The gradient vector is tangent to the sphere at point P, and, as P approaches the North or South Pole, r approaches ±1; cos δ and cos ϵ approach 1; p, q, and cos ϕ approach 0; $\partial h/\partial \mathit{\varphi}$ and $\partial h/\partial \mathit{\lambda}$ become less important in Eq. (2.13); and U and V become more nearly perpendicular. Using perpendicular unit vectors W and N, the gradient is commonly written as
This common form is not valid near the poles because of cos ϕ in the denominator, while the new symmetric form, Eq. (2.13), is valid everywhere, a benefit of using the symmetric equations. The equivalence between the common form and the new form of ∇h is shown in Appendix A. If ∇h is treated as a velocity, then the specific angular momentum with which it is associated is
The divergence of a horizontal vector $\mathbf{\text{D}}=(d,e,f)$ at point P is
which is equivalent to the common form on the surface of a sphere: $[\partial (\mathbf{\text{D}}\cdot \mathbf{\text{W}})/\partial \mathit{\lambda}+\partial (\mathbf{\text{D}}\cdot \mathbf{\text{N}}\mathrm{cos}\mathit{\varphi})/\partial \mathit{\varphi}]/R\mathrm{cos}\mathit{\varphi}$. The first two forms of Eq. (2.16) are equivalent because $\mathbf{\text{D}}\cdot \mathbf{\text{P}}=\mathrm{0}$ and
Similar reasoning shows the equivalence between the third and fourth forms of Eq. (2.16). If D is horizontal velocity, D = S, then the divergence of S can also be written as
Also,
Noting that ∇h is perpendicular to P and reasoning similar to Eq. (2.17), the Laplacian is
The Laplacian of Eq. (2.20) is equivalent to the common form on the surface of
a sphere:
$[{\partial}^{\mathrm{2}}h/\partial {\mathit{\lambda}}^{\mathrm{2}}+\mathrm{cos}\mathit{\varphi}\partial (\mathrm{cos}\mathit{\varphi}\partial h/\partial \mathit{\varphi})/\partial \mathit{\varphi}]/{R}^{\mathrm{2}}{\mathrm{cos}}^{\mathrm{2}}\mathit{\varphi}$.
The curl of a horizontal vector $\mathbf{\text{D}}=(d,e,f)$ is
If ∇×D were completely vertical, then $(\mathrm{\nabla}\times \mathbf{\text{D}})\times \mathbf{\text{P}}$ would be 0. It is not; in fact, $(\mathrm{\nabla}\times \mathbf{\text{D}})\times \mathbf{\text{P}}=\mathbf{\text{D}}/R$. The upward vertical component of the curl of D is
which is equivalent to the common form on the surface of a
sphere:
$[\partial (\mathbf{\text{D}}\cdot \mathbf{\text{N}})/\partial \mathit{\lambda}\partial (\mathbf{\text{D}}\cdot \mathbf{\text{W}}\mathrm{cos}\mathit{\varphi})/\partial \mathit{\varphi}]/R\mathrm{cos}\mathit{\varphi}$.
If
D = S, then the vertical component of the curl of
S can also be written as
The upward vertical component of relative vorticity is
which is equivalent to the common form on the surface of a sphere: $[\partial n/\partial \mathit{\lambda}\partial (w\mathrm{cos}\mathit{\varphi})/\partial \mathit{\varphi}]/R\mathrm{cos}\mathit{\varphi}$.
If G and H are differentiable vectors in three space, then a wellknown identity using the Cartesian del operator is
If G is radially aligned and of constant magnitude (P for example), then ${\mathrm{\nabla}}_{C}\times \mathbf{\text{G}}=\mathbf{\text{0}}$. If G is a radially aligned unit vector and H is perpendicular to G, then
The above relationships apply to the present ∇ (or ∇_{S}) operator and are shown in the relationships of Eqs. (2.22) and (2.23).
Many of the new symmetric forms have no varying quantities outside their derivatives and can be integrated using Green's theorem. Proofs of several of the equivalences used above are available at https://aom.giss.nasa.gov (last access: 15 November 2018).
2.3 Differential form of the shallow water equations
The differential form for conservation of mass, using Eq. (2.18), is applied to mass per unit area:
This symmetric form is equivalent to the common equation:
The threecomponent advective form for specific angular momentum is
where f (s^{−1}) is the Coriolis parameter and Φ (m^{2} s^{−2}) is the fluid top geopotential:
Replacing h with hA in Eqs. (2.27) and (2.28) shows the differential form for momentum advection in flux form:
The last line of Eq. (2.32) evaluates to zero because of conservation of mass. Inside the square brackets, the penultimate line of Eq. (2.32) shows the form that may be used for advection of A at a grid edge when the edge is perpendicular to the unit vector W. Replacing A in the square brackets in that penultimate line with c=wcos ϕ and dividing by cos ϕ yields
which are the common shallow water forms for the time derivative, the advective terms, and the metric term of eastward velocity. Application of the momentum conservation form to A thus includes the metric term. Although the common shallow water form for $\partial w/\partial t$ can incorporate the metric term into the advective terms, $\partial n/\partial t$ cannot conveniently do so. In the present formulation, all three components act like $\partial c/\partial t$ with the metric term included into the advective terms.
Symmetric versions of the shallow water equations using vorticity and divergence or vectorinvariant form are shown at https://aom.giss.nasa.gov. Model IB does not use this form nor Eq. (2.33).
3.1 Alignment of velocity or specific angular momentum
Alignment of the three momentum components at P means that $\mathbf{\text{A}}=(u\mathrm{cos}\mathit{\delta},v\mathrm{cos}\mathit{\u03f5},w\mathrm{cos}\mathit{\varphi})$ and S are perpendicular to P ($\mathbf{\text{A}}\cdot \mathbf{\text{P}}=\mathbf{\text{S}}\cdot \mathbf{\text{P}}=\mathrm{0}$). Application of a horizontal acceleration vector to the components maintains alignment. Thus, the pressure gradient force vector and the Coriolis force maintain alignment. The pressure gradient force obtained via Green's theorem and advection may distort alignment. The following procedure brings distorted momentum components back into alignment.
Given unaligned velocity components u, v, and w at point P, determine the least square fit velocity vector S_{NEW} that is horizontal (${\mathbf{\text{S}}}_{\text{NEW}}\cdot \mathbf{\text{P}}=\mathrm{0}$) and best matches the components weighted by the square of the distance to their axes.
If u, v, and w were already aligned, then t=0, where
When s is minimized, alignment of distorted components u, v, and w produces
Performing this analysis with specific angular momentum components yields a=ucos δ, $t=pa+qb+rc=\mathbf{\text{P}}\cdot \mathbf{\text{A}}$, and
The minimization technique applied above is equivalent to projecting an unaligned S or A onto the tangent plane of the spherical surface at P. When P is close to an axis pole, say r is close to ±1, then cos ϕ is close to 0, w or c is most strongly modified, and u and v are modified weakly. Alignment has the same purpose as the “Lagrange multiplier” term of Coté (1988).
The benefit of using three aligned components instead of two for advection on a particular cubedsphere grid is shown at https://aom.giss.nasa.gov.
3.2 Pressure gradient force for model IB
Change of velocity by the pressure gradient force is proportional to the gradient of the field top geopotential Φ:
where Δt is the time step. Change of specific angular momentum, using Eq. (2.15), is
Application of the pressure gradient force to velocity averaged over an arc usually involves interpolating Φ to the corners of the arc and to Φ on either side of the arc. Application to velocity averaged over a cell is conveniently performed using Green's theorem. This computation is discussed first for cubedsphere models and later for icosahedral models, like IB.
Acceleration of velocity averaged over a primary cell of a cubedsphere model starts by knowing Φ averaged over the cell's edges from interpolation of primary cell values. For each cell mean velocity component, edge Φ multiplied by the cosine of the angle between the component's unit vector and a unit vector outwardly perpendicular to the edge is integrated with respect to distance around the cell's edges. This integral, multiplied by the time step and divided by the cell's area, accelerates the mean velocity component of the cell. On Face 1 (Fig. 1), where the outward perpendicular directions of the right and left edges are V and −V and those of the top and bottom edges are −U and U, and where dχ_{W}=Rcos ϕdλ is the spatial differential in the direction of W,
where L is the arc length of an edge and K is the primary cell area [$I\mathrm{1}:I,J\mathrm{1}:J$] in the computer implementation. For the right edge [$I,J\mathrm{1}:J$], LW⋅Vcos ϕ is evaluated as an integral over arc length as
Similarly, the top edge [$I:I\mathrm{1},J$] is integrated over arc length as
Integrating counterclockwise around the grid cell, similar formulas for different components yield
If all four edge values of Φ are the same, then cancellation of the P values causes cell mean ΔA to be 0.
Relationships like Eqs. (3.10) or (3.11) apply to any great circle arc, not just cubedsphere edges. To check this, any arc from P_{1} to P_{2} is a subset of a great circle that intersects the Equator at two points $\mathbf{\text{I}}=(i,j,\mathrm{0})$ and −I; the longitudinal angular rotation coordinate around this axis is ξ and the latitudinal coordinate from −I to I is η. P is computed as $(i\mathrm{sin}\mathit{\eta}j\mathrm{cos}\mathit{\xi}\mathrm{cos}\mathit{\eta},i\mathrm{cos}\mathit{\xi}\mathrm{cos}\mathit{\eta}+j\mathrm{sin}\mathit{\eta},\mathrm{sin}\mathit{\xi}\mathrm{cos}\mathit{\eta})$, and the horizontal unit vector perpendicular to the arc at point P is $\mathbf{\text{F}}=\mathbf{\text{I}}\times \mathbf{\text{P}}/\mathbf{\text{I}}\times \mathbf{\text{P}}$. $\mathbf{\text{W}}\cdot \mathbf{\text{F}}=\mathrm{sin}\mathit{\xi}\mathrm{sin}\mathit{\eta}/\mathrm{cos}\mathit{\varphi}$.
Renaming the vertices of a spherical polygon with N arcs P_{0}, P_{1}, ... P_{N}=P_{0} counted counterclockwise around the cell, and generalizing Eq. (3.12), yields
Appendix B shows that if Φ is a linear function of η in the arc, then integrals that include Φ in Eq. (3.13) can be computed in closed form knowing Φ at each end of the arc, and Eq. (3.14) can be modified to use Φ_{n} and Φ_{n−1} instead of ${\mathrm{\Phi}}_{n\mathrm{1}/\mathrm{2}}$, as shown in Eq. (B4). This is applicable to model IB where primary mass and Φ are centered at the momentum cell's vertices and Φ may be interpolated along the edges. For the Sect. 4 tests performed on model IB, Eqs. (3.14) and (B4) produce very similar results. Equation (3.14) is used because it is simpler.
3.3 Advection for model IB
Change in primary cell mass by advection, in flux form, is derived from Green's theorem and Eq. (2.27).
where M_{P} (kg s^{−1}) is the outwardly transported mass at primary cell edges. M_{P} at each edge is the product of the mass per unit area at the edge, the perpendicular velocity component, and the length of each edge. Using the same notation of variables I, η, and F as in the prior subsection, M_{P} along an edge from η_{1} to η_{2} is computed as
where E is the unit vector parallel to the edge. The final form of Eq. (3.16) shows the elegance of using specific angular momentum as the transport variable; integration around a cell's boundary follows the position vector without a perpendicular velocity computation. As with Eq. (3.14) for the pressure gradient force, ${h}_{n\mathrm{1}/\mathrm{2}}{\mathbf{\text{A}}}_{n\mathrm{1}/\mathrm{2}}$ can be the average value over the counterclockwise arcs, P_{n−1} to P_{n}, so that
The inherent symmetry between variables of the shallow water equations is evidenced by Eqs. (3.14) and (3.17). As shown in Appendix B, variables A or hA may also be linearly interpolated with respect to η, so that the formula for Δh is closer to Eq. (B4). This is applicable to model IB because A is located at the primary cell corners. For Sect. 4 tests performed on IB, Eq. (3.17) is used; the more complex version with linearly interpolated A produces similar results.
The change in grid cell mean relative angular momentum by advection, in flux form, is equal to the summation of the angular momentum fluxes V_{T} (kg m s^{−2}) at the edges multiplied by Δt. V_{T} is the product of A multiplied by M_{T}, the mass flux of momentum cell edges designated by subscript T.
If ${h}_{n\mathrm{1}/\mathrm{2}}$ and ${\mathbf{\text{A}}}_{n\mathrm{1}/\mathrm{2}}$ are average values over the arcs surrounding the momentum cells, then the change in angular momentum averaged over a cell is
For stability purposes, if A is constant throughout several cells in a neighborhood, even though M_{T} may be irregular, ΔA should be 0 throughout the neighborhood. If principal mass cells and momentum cells are identical, M_{P}=M_{T}, this principle is obeyed automatically, assuming that V_{T}=AM_{T} and A is the constant.
Primary and momentum cells are not identical for model IB; triangular momentum cells are centered at the corners of primary cells. The solution for momentum cell mass presented here is different from that of Ringler et al. (2010), where M_{P} and M_{T} are labeled F_{e} and ${F}_{e}^{\u27c2}$, respectively. The area of a momentum cell is the summation of quadrilateral areas inside three touching primary cells, and the momentum cell's mass is the summation of the three quadrilaterals' masses (Figs. 2 or 3). Mass is always uniformly distributed with respect to area in primary cells, and the change in mass in a momentum cell during an advective time step is determined by the change in mass of primary cells. Mass fluxes M_{T} must be chosen so that their mass changes into a momentum cell match the change caused by primary mass changes. M_{T} is initially computed along an arc between two primary cell centers, but M_{T} is then modified as explained in the caption to Fig. 3. M_{T} halves are minimally adjusted so that the change in mass per unit area of each quadrilateral area of a primary cell is identical during the time step. Adjusted M_{T} halves adjust full M_{T}. This technique satisfies the stability principle of the prior paragraph.
3.4 Coriolis force
For computer implementations of symmetric equations created so far, momentum components have not been defined on staggered locations; all three components reside at the same locations. For each velocity component, the other two components determine the velocity that is perpendicular to the first component. Thus, $(pcra)/\mathrm{cos}\mathit{\u03f5}$ is the velocity component that aligns with P×V. The Coriolis acceleration (m s^{−2}) acting on v is $\mathrm{2}\mathrm{\Omega}\mathrm{sin}\mathit{\varphi}(pcra)/\mathrm{cos}\mathit{\u03f5}$ and the acceleration acting on b is 2Ωsin ϕ(pc−ra). Each specific angular momentum component is accelerated by its local components:
3.5 Icosahedral models
Starting with 12 vertices and 20 triangles of the icosahedron, new vertices for the raw grid are formed from the center points of triangular edges. After m iterations of this process, the triangular lattice contains $\mathrm{2}+\mathrm{10}\cdot {\mathrm{4}}^{m}$ vertices as described by Stuhne and Peltier (1999). m is called the grid level. These vertices are the centers of primary cells where scalar quantities such as h are defined. Primary cells are mainly irregular hexagons, although there are 12 regular pentagons centered on the vertices of the original icosahedron. For model IB, each three nearest primary cell centers form the vertices of a spherical triangle where momentum is defined and which contains the common corner of the three primary cells; these corners designate momentum cell centers. There are 20⋅4^{m} triangular momentum cells that cover the whole sphere; it is called the “dual” grid by Ringler et al. (2010) and others; see Fig. 2.
Starting from the raw grid, Heikes et al. (2013) reposition and optimize the location of the cell centers producing a tweaked grid. For both the raw and tweaked grids, arc edges of primary cells are the perpendicular bisectors of nearby primary cell centers. The centroid grid's primary cell centers and corners are also repositioned from raw grid locations, but primary cell edges are neither perpendicular nor bisectors of primary center arcs. For the raw and tweaked grids, cell centers do not coincide with the cell centroid for neither primary nor momentum cells. For the centroid grid, centers and centroids coincide for both primary and momentum cells (see Table 1). The centroid grid here is different from the centroidal Voronoi tessellation grid of Du et al. (2003) which had only primary cell coincidence and different from the spherical centroidal barycenter method of Miura and Kimoto (2005) which had only momentum cell coincidence. Each icosahedral model, IB or IK, may use any grid; a capital letter (R indicates raw, T indicates tweaked, and C indicates centroid) is affixed to a model's label that show the grid. IKC works poorly.
Table 1 shows various properties of the three grids for grid levels 4 through 8. A significant property is “smallest arc length from primary cell corner to D divided by half of ArcA length”; ArcA is between two primary cell corners and D is the intersection of ArcA and the arc between the primary centers. The raw grid value in Table 1 stays at 81%, while the tweaked grid converges to unity with increasing resolution. For this reason, IKT is superior to IKR. IB performs edge computations using only values of the two primary cells, ignoring the fact that D is not the center of ArcA. IK, CSK, and the icosahedral model of Lee and MacDonald (2009) use additional primary cell data to compute cell edge values, but this makes things more difficult when applying such models to an ocean domain.
Considering the most extreme momentum cells in Table 1, the raw grid triangular momentum cells of model IB are more equilateral than are those of the tweaked grid. Momentum cells are more important for model IB, and consequently IBR is as good as IBT. An expanded version of Table 1 with additional parameters is available at https://aom.giss.nasa.gov.
The time scheme is leapfrog initialized every 8 to 10 time steps by forward–backward steps. Alignment of A (Sect. 3.1) is performed after every step. The geometry subroutine is complex, but the computational subroutines are simple because angular momentum components are treated identically and use the same lines of code. M_{P} is computed using twopoint secondorder differencing. Some momentum computations use upstream differencing. This causes the flow to be stable but slightly diffusive. However, no other smoothing or filtering subroutines are needed. Computations are performed as locally as possible, and consequently the icosahedral models are suitable for an ocean domain or stepmountain atmosphere.
For model IB, the unadjusted M_{T} is computed using secondorder differencing, but implementing the stability adjustment (Sect. 3.3) adds complexity to the code. When computing V_{T}=AM_{T}, A is a mixture of half secondorder and half linear upstream advection. Unlike IK and the unfiltered lat–long Bgrid scheme of Arakawa and University of California (1972), model IB seldom generates alternating patterns in the height field. The Fortran source code for IB is available at https://aom.giss.nasa.gov or on Zenodo at https://doi.org/10.5281/zenodo.1313736.
Model IK uses h at primary grid cell corners interpolated from its value of the three primary cells touching the corner. Different interpolation formulas were tested, but Eq. (3.5a) of Heikes et al. (2013) was deemed best.
The following two lat–long climate models, reduced to one layer for the shallow water equations, were applied to the test cases below. Arakawa's secondorder Bgrid (velocity components defined at primary cell corners) lat–long model with the NASA Goddard Institute for Space Studies (GISS) ModelE's pole modifications and filters on mass and velocity (Schmidt et al., 2006) is labeled LLB. Arakawa's secondorder Cgrid (velocity components perpendicular to primary cell edges) lat–long model as modified by Russell (2007) is labeled LLC, which also applies a filter to velocity components. The symmetric equation models (CSK, IB, and IK) do not apply any external filters to mass or momentum other than alignment. CSK performs more twopoint interpolations of variables than do the other models because of its parallelogramshaped cells; this probably explains CSK's poor performance for some of the tests. IB, the simplest model, has few choices to tune, the main ones being the mixture of secondorder versus linear upstream for momentum advection (0.5 and 0.5) and the choice of grid (raw, tweaked, or centroid). The proprietary IK model is more complicated, and the tweaked grid is superior for this model. Table 2 shows the abbreviations and grids of the models used below.
Symmetric equation models IBR, IBT, IBC, IKT, and CSK are represented; IKR and IKC are not, being worse than IKT. Each model was tested with different initial conditions and different horizontal resolutions that approximate 4, 2, 1, and 0.5^{∘}; the symmetric models were tested at 0.25^{∘} in addition. Approximately 1^{∘} uses 64 primary cells along the triangular edge of an icosahedron, 88 cells along the edge of a cube face, and 180 latitude bands for LL models; for other resolutions, these numbers are doubled or halved.
(Schmidt et al., 2006)(Russell, 2007)4.1 Solid body rotation without bottom topography (SBR)
Using the parameters of Test Case 2 of Williamson et al. (1992), a perfect model would maintain the initial mass and velocity fields indefinitely. For a secondorder model, the mass and velocity errors should decrease by a factor of 4 when doubling the horizontal resolution. Figure 4 shows the areaweighted root mean square l_{2} norm of the mass field after 5 days of integration. This figure may be compared to Fig. 2 of Stuhne and Peltier (1999) and to Fig. 3 of Lee and MacDonald (2009), whose data are averaged over the first 5 days. Both older icosahedral models, LLB and LLC, show error reduction factors of about 4 for doubling resolution. Models IBR and IBC show factors close to 3 that decrease with finer resolution but increase with less upstream advection of momentum. IKT and CSK show impressive reduction factors of 6 and 5, respectively. As mentioned in the third paragraph of Sect. 3.5, with finer resolution the intersection point D gets closer to the center of primary cell edge arc for model IKT; this increases the precision of this model.
Figure 5 shows the l_{1} and l_{∞} norms as a function of time for the first 5 days of the 2^{∘} resolution models IBR, IKT, CSK, and LLC. For all resolutions of the symmetric equation models, the mass field error increases linearly with time; the error after 50 days is nearly 10 times that after 5 days. For most resolutions of LLC, the mass field error after 80 days is the same as the error after 5 days. This is also the case for coarsest resolutions, 4 and 2^{∘}, of LLB, but 1^{∘} resolution LLB diverges after 78 days for any reasonable time step and 0.5^{∘} resolution produces polar instabilities although it survives for 100 days and beyond. The present symmetric models of Fig. 5 can be compared to Fig. 8 of Heikes and Randall (1995) that shows smaller errors, Fig. 2 of Stuhne and Peltier (1999) that shows larger errors, and Fig. 2 of Lee and MacDonald (2009) that is comparable to IBR. The older models show less error growth with increasing time. Both the present 2^{∘} LLC and the Cgrid model of Heikes and Randall (1995) at 4^{∘} resolution show occasional instabilities that do not affect the longterm error growth.
In general, the mass field error of IBR is less than that of IBT, often exceeding more than 10 % for the coarser resolutions, but closer to 1 % for 0.25^{∘} resolution. All of the symmetric equation models use some amount of upstream advection of momentum. This causes a continual loss of total energy and of structure to the prognostic variables, and is the reason for the continual degradation of the l_{2} norms in time. The Arakawa Bgrid and Cgrid schemes are designed to approximately conserve total energy.
4.2 Rossby–Haurwitz wave 3 (RH3)
For cubedsphere models that initialize from repetitive wave 4 initial conditions, each of the four wavelengths lies above the same grid arrangement and copies remain identical when the wave's axis passes through the center of a cube face as is the case for CSK. The same statement applies to icosahedral models with wave 5 initial conditions when the axis passes through a vertex of the icosahedron as is the case for IB and IK. The initial conditions used here are Rossby–Haurwitz wave 3, 3 being relatively prime to both 4 and 5. With wave 3, separate wavelengths diverge among themselves for models CSK, IB, and IK, and there is more variety in the errors that occur. Plots (not displayed) show that CSK for Rossby–Haurwitz wave 4 and IB or IK for Rossby–Haurwitz wave 5 maintain their wavelength shape, but other situations do not. Some results of CSK and IBR for Rossby–Haurwitz waves 3, 4, and 5 are shown at https://aom.giss.nasa.gov under “RHn”.
Figure 6 shows horizontal plots of h (h_{S}=0) after 45 days of integration at 2^{∘} resolution for six models: IBR, IBC, IKT, CSK, LLB, and LLC. Not displayed is model IBT, which no longer has three peaks at 2^{∘} resolution but is comparable to IBR for finer resolutions. Although IBC has destroyed the pattern for 2^{∘} resolution, it actually looks better than IBR for 0.5^{∘} after 45 days. CSK performed well for SBR initial conditions, but it is the worst model for RH3, being unable to maintain the RH3 pattern for 45 days except for 0.25^{∘} resolution. Its reduced quality is partially caused by the grid edge lines not being perpendicular. LLB at 2^{∘} resolution is beginning to develop alternating patterns in the orange color of Fig. 6, but LLB is excellent for finer resolutions.
Table 3 shows the change in spectral specific kinetic energy for wave numbers 0 and 3 and the change in total energy for models IBR, IBT, IBC, IKT, CSK, LLB, LLC, and National Center for Atmospheric Research spectral transform model (STSWM) at T42 resolution (Hack and Jakob, 1992) after 40 days of integration. Icosahedral and cubedsphere output was interpolated to a highresolution lat–long grid before performing spectral decomposition. For each model's finest resolution, the reduction of wave 3 energy is between 0 and 5.6 %. Courser resolutions show greater reductions in wave 3 energy and CSK is less similar to the highresolution results than are the icosahedral models. All of the symmetric models have some amount of upstream advection of momentum which reduces (kinetic and) total energy, as shown in Table 3. Model IB has more total energy loss than other models, but it is smoother and more stable.
Except for LLC, which maintains its pattern well for 100 days, the other 1^{∘} resolution models deviate from the wave 3 pattern after different numbers of days: IBR – 45, IBT – 45, IBC – 37, IKT – 50, CSK – 26, and LLB – 63. Both LLB and LLC 4^{∘} models diverge in the first day of integration, but all other models and resolutions survive for at least 100 days; the final pattern is smooth but may be unrecognizable.
4.3 Initial solid body rotation but with Earth's bottom topography (SBRZ)
The initial velocity is eastward, 50 m s^{−1} at the Equator multiplied by cosine of latitude. The initial global mean height field top is 10 000 m with the latitudinal distribution to maintain itself were there no bottom topography. But h is reduced by the Earth's bottom topography which has peaks as high as 5600 m. SBRZ is a more realistic and severe than Test Case 5 of Williamson et al. (1992) because it has faster initial velocity, higher mountains, and much larger topography gradients that increase with finer resolution, particularly at the Andes. In Williamson, topography gradients are independent of resolution. Greater accelerations by the pressure gradient force partially explain why the time step is not inversely proportional to the linear horizontal resolution.
Although the lat–long models use Fourier polar filters on east–west mass flux and pressure gradient force and other filters on prognostic variables, the symmetric equation models do not. The stability of IB, IK, and CSK is maintained by using the proper amount of linear upstream advection of momentum.
For the videos discussed later, the leapfrog time step of each model is interrupted every eight time steps, and 2Δt is an integral division of 900 s, the video time step. Table 4 shows, for each model and resolution, the largest time step for which the model survives for 50 days without major instabilities or after what day the model diverged. Time steps were limited to be greater than or equal to onethird that of IBR for the same approximate resolution. All IB models are stable. IBC requires smaller time steps than IBR or IBT because IBC has some smaller grid cells; this discrepancy increases with resolution, as shown in Table 4. IK models are not as stable as IB and require time steps that are half that of IB for finer resolutions. Except for 4^{∘}, LLC diverges for any time step during the 50 days. LLB diverges for resolutions 4 and 0.5^{∘}, and requires very short time steps for 1^{∘} resolution.
Table 5 shows the change in kinetic energy after 50 days of integration for different models and resolutions. For resolutions 1^{∘} or finer, the kinetic energy for all models that survive is within 5 % of the original value. Kinetic energy is not a conserved quantity, but for coarseresolution models, its numerical reduction coincides with the washing out of highs and lows of the height field, as shown in the videos discussed subsequently.
Figure 7 shows horizontal plots of h+h_{S} for resolutions 2^{∘} and 0.5^{∘}, and models IBR, IKT, and CSK after 34 days of integration. IBT and IBC are not displayed, but their plots are similar to IBR. Models IKT and CSK display alternating patterns aligned with grid lines in South America. This occurs for both 2 and 0.5^{∘} resolutions but is more difficult to see for 0.5^{∘} because of the small size of the printed page. For 2^{∘} models shown in Fig. 7, IBR is most similar to the 0.5^{∘} models. There is discrepancy among the nondisplayed 0.25^{∘} symmetric models after 50 days; IBR and IKT are more similar to each other than to CSK. Nondisplayed LLB diverges at 0.5^{∘} resolution, but at 1^{∘}, it is similar to IBR.
Videos of 50 simulated days for models IBR, IKT, CSK, LLB, and LLC are displayable of limited quality on YouTube for all resolutions (https://www.youtube.com/user/NASAGISStv/videos, last access: 15 October 2018). The label for each video contains the model acronym, grid edge resolution, and “Z”, e.g., IBR64Z, CSK88Z, or LLB180Z. IKT models show frequent multicell long linear alternating data that may last for several hours. IBR models also show the alternating patterns, but they are much less extensive in occurrence, duration, and spatial distance. CSK is somewhere between IBR and IKT in terms of these alternating patterns. LLB infrequently displays a checkerboard pattern for three or four adjacent cells. LLC, before running into polar problems, is smooth. At resolutions 1^{∘} and finer and up to 35 days of integration, all of the models that do not diverge are similar to each other.
This paper presents symmetric calculus operators on the surface of a sphere that are projected from threedimensional Cartesian formulas. Symmetric equations are simplified by using specific angular momentum on the unit sphere, A, instead of velocity. Components of A align with the fixed orthogonal axes, whereas the velocity unit vectors U, V, and W are not orthogonal and vary with location. A is continuous everywhere, whereas u, v, and w are discontinuous at their respective poles. A simplifies the equation for horizontal velocity S (Eq. 2.6) and several formulas in Sect. 2.2. Advection of A does not use the metric term, a correction term needed when advection is applied to velocity or linear momentum on the sphere. Applications of Green's theorem invoke the elegant formulas: Eq. (3.14) for the pressure gradient force, and Eqs. (3.17) and (3.19) for advection. All components of relative angular momentum, hA, are conserved without time truncation errors by the flux form advection (except for alignment) used by the symmetric models; LLC conserves the north–south axis component of relative angular momentum by advection, but other components and models lack conservation.
Summarizing the results from the test cases, no one model is clearly superior in all tests for all resolutions and times of integration.

Solid body rotation without bottom topography (SBR). For most resolutions, icosahedral models have the lowest relative mass error at day 5 of the simulation (Fig. 4). The Arakawa Cgrid model, LLC, maintains the same mass field error after 80 days that it had after 5 days, whereas the symmetric equation models all have mass field errors that increase linearly with time. LLC must be considered best for this test.

Rossby–Haurwitz wave 3 (RH3). Large losses of wave number 3 kinetic energy after 40 days (Table 3) show model deficiencies of which CSK is the most egregious. LLC maintains the wave 3 shape longer than the other models (Fig. 6) and again is the best. For Rossby–Haurwitz wave 4, CSK maintains the wavelength shape better than IBR does, while for Rossby–Haurwitz wave 5, IBR is superior to CSK for longer integrations (https://aom.giss.nasa.gov under “RHn”).

Initial solid body rotation with Earth's topography (SBRZ). Polar problems cause the lat–long models, LLB and LLC, to diverge for different resolutions including 0.5^{∘}. IKT and CSK to a lesser extent produce linear alternating patterns that are much diminished in IBR (Fig. 7 and the videos). IBR and IBT use a longer time step (Table 4) than do other models, and, after 34 days, 2^{∘} IBR is more similar to the higherresolution 0.5^{∘} models than are 2^{∘} IKT and CSK (Fig. 7). The icosahedral Bgrid models, IBR and IBK, are best.
The SBRZ test is the most difficult but the most realistic test. Given that these onelayer models will be the basis for multilayer climate models with realistic topography, IB models must be considered the best overall. Parallelogramshaped grid cells in CSK with nonperpendicular grid line edges lasting over large swaths of the globe cause a systematic error in numerical flow, most noticeably evident in the RH3 test case. Weller et al. (2012) also conclude that “the hexagonal icosahedron gives the most accurate results” for several test cases. Although IKT's reduction factor of 6 for doubling resolution with SBR initial conditions is impressive, it generates frequent alternating linear patterns for SBRZ conditions, as shown in the YouTube videos. The Williamson et al. (1992) test cases are inadequate for onelayer models that will be expanded to climate models simulating Earth; Williamson's Test Case 5 is much less demanding than using SBRZ. Grid imprinting errors in IB and IK are comparable to errors reported by other researchers (Stuhne and Peltier, 1999; Lee and MacDonald, 2009); the errors are not worse than would occur using two horizontal velocity components.
The symmetric equation models presented here use the smallest sensible grid cell stencil needed for a computation. Enlarging the stencil, as used by Lee and MacDonald (2009), may improve the results for tests SBR and RH3 by using fourthorder differencing or other methods, but in ocean domains or stepmountain atmospheres, large stencils require many different formulas for flow near ocean coastlines, based on their shape, and each formula is different from that used in the interior. CSK, with its parallelogramshaped cells, has more difficulty in conforming to an ocean domain than do the icosahedral models.
As noted in Sect. 3.5, IB uses only two adjacent primary cell centers when performing computations on their common edge, even though point D is not the center of the primary cell edge. Except for the stability requirement of mass fluxes entering momentum cells (Sect. 3.3 and Fig. 3), the computational subroutines of IB are extremely simple. There are no separate lines of code for angular momentum; all components use the same lines. IK is slightly more complicated, but CSK is worse, requiring frequent interpolation to position variables.
Flux form velocity, represented by two horizontal components, has significant problems where the coordinates become discontinuous. An improvement has been to use forms of the shallow water equations where scalar quantities such as potential vorticity, specific kinetic energy, divergence, stream function, and velocity potential are continuous over the whole sphere and from which the local horizontal velocity can be resurrected, or integrated using the manipulated scalar quantities. Deficiencies of these methods are complexity of understanding and computer coding. This paper presents another method: vector angular momentum is continuous over the whole sphere and its application via the symmetric equations is simpler than using velocity. Each component of relative angular momentum is conserved by flux form advection without discontinuities. Further work is needed to determine the practical advantages that one scheme may have over others.
Fortran source code for model GISS:IB is available at https://aom.giss.nasa.gov or on Zenodo at https://doi.org/10.5281/zenodo.1313736 (Russel et al., 2018).
Using the notation for F, η, and ξ of Sect. 3.2,
Assume that Φ was not constant throughout the arc but instead was a linearly function of η from Φ_{1} to Φ_{2}. Then, the arc integral of Eq. (3.13) with Φ included is
where r_{C} is the value of r at the center of the arc from η_{1} to η_{2}. As the resolution of the grid becomes finer, r_{C} approaches $({r}_{\mathrm{1}}+{r}_{\mathrm{2}})/\mathrm{2}$ and $\mathrm{sin}\left[\right({\mathit{\eta}}_{\mathrm{2}}{\mathit{\eta}}_{\mathrm{1}})/\mathrm{2}]$ approaches $({\mathit{\eta}}_{\mathrm{2}}{\mathit{\eta}}_{\mathrm{1}})/\mathrm{2}$, and the integral of Eq. (B3) approaches $R({r}_{\mathrm{2}}{r}_{\mathrm{1}})({\mathrm{\Phi}}_{\mathrm{1}}+{\mathrm{\Phi}}_{\mathrm{2}})/\mathrm{2}$, which is the result of Eq. (3.13) multiplied by Φ_{1.5}.
Using Eq. (B3) instead of Eq. (3.13), the formula of Eq. (3.14) is replaced with
Most of the formulas and equations of Sect. 2 have probably been known for at least two centuries, but were rederived from scratch by the authors. Use of three horizontal components, vector angular momentum instead of velocity, and the combination of both is new. The Fortran programs used for this paper are original and were written from scratch.
The authors declare that they have no conflict of interest.
The origins of this paper came from discussions with Maxwell Kelley in 2012 and 2013. He has made contributions to the present paper, and models IK and CSK are based on his ideas that may be published in the future. The authors thank Ross Heikes of Colorado State University for the tweaked grid locations. Rainer Bleck was helpful with knowledge of other icosahedral models. Robert Schmunk assisted in creating the YouTube video files. The paper was significantly improved by the editorial process of Monthly Weather Review.
Climate modeling at GISS is supported by the NASA Modeling, Analysis,
and Prediction program, and resources supporting this work were provided
by the NASA HighEnd Computing Program through the NASA Center for Climate
Simulation at Goddard Space Flight Center.
Edited by: Paul Ullrich
Reviewed by: two anonymous referees
Adcroft, A., Campin, J.M., Hill, C., and Marshall, J.: Implementation of an AtmosphereOcean General Circulation Model on the Expanded Spherical Cube, Mon. Weather Rev., 132, 2845–2863, https://doi.org/10.1175/mwr2823.1, 2004. a
Arakawa, A. and University of California, L. A. D. o. M.: Design of the UCLA General Circulation Model, Numerical simulation of weather and climate: Technical report, Department of Meteorology, University of California, available at: https://books.google.com/books?id=nzEESwAACAAJ (last access: 15 November 2018), 1972. a
Coté, J.: A Lagrange multiplier approach for the metric terms of semiLagrangian models on the sphere, Q. J. Roy. Meteor. Soc., 114, 1347–1352, https://doi.org/10.1002/qj.49711448310, 1988. a, b
Du, Q., Gunzburger, M. D., and Ju, L.: Constrained Centroidal Voronoi Tessellations for Surfaces, SIAM J. Sci. Comput., 24, 1488–1506, https://doi.org/10.1137/s1064827501391576, 2003. a
Hack, J. and Jakob, R.: Description of a Global Shallow Water Model Based on the Spectral Transform Method, NCAR Technical Note NCAR/TN343+STR, https://doi.org/10.5065/d64b2z73, 1992. a
Heikes, R. and Randall, D. A.: Numerical Integration of the ShallowWater Equations on a Twisted Icosahedral Grid. Part I: Basic Design and Results of Tests, Mon. Weather Rev., 123, 1862–1880, https://doi.org/10.1175/15200493(1995)123<1862:niotsw>2.0.co;2, 1995. a, b, c, d
Heikes, R. P., Randall, D. A., and Konor, C. S.: Optimized Icosahedral Grids: Performance of FiniteDifference Operators and Multigrid Solver, Mon. Weather Rev., 141, 4450–4469, https://doi.org/10.1175/mwrd1200236.1, 2013. a, b
Läuter, M., Giraldo, F. X., Handorf, D., and Dethloff, K.: A discontinuous Galerkin method for the shallow water equations in spherical triangular coordinates, J. Comput. Phys., 227, 10226–10242, https://doi.org/10.1016/j.jcp.2008.08.019, 2008. a
Lee, J.L. and MacDonald, A. E.: A FiniteVolume Icosahedral ShallowWater Model on a Local Coordinate, Mon. Weather Rev., 137, 1422–1437, https://doi.org/10.1175/2008mwr2639.1, 2009. a, b, c, d, e, f
Leer, B. V.: Towards the ultimate conservative difference scheme. IV. A new approach to numerical convection, J. Comput. Phys., 23, 276–299, https://doi.org/10.1016/00219991(77)90095x, 1977. a
Masuda, Y. and Ohnishi, H.: An Integration Scheme of the Primitive Equation Model with an IcosahedralHexagonal Grid System and its Application to the Shallow Water Equations, J. Meteorol. Soc. Japan. Ser. II, 64A, 317–326, https://doi.org/10.2151/jmsj1965.64a.0_317, 1986. a
Miura, H. and Kimoto, M.: A Comparison of Grid Quality of Optimized Spherical Hexagonal–Pentagonal Geodesic Grids, Mon. Weather Rev., 133, 2817–2833, https://doi.org/10.1175/mwr2991.1, 2005. a
Prather, M. J.: Numerical advection by conservation of secondorder moments, J. Geophys. Res., 91, 6671, https://doi.org/10.1029/jd091id06p06671, 1986. a
Putman, W. M. and Lin, S.J.: Finitevolume transport on various cubedsphere grids, J. Comput. Phys., 227, 55–78, https://doi.org/10.1016/j.jcp.2007.07.022, 2007. a, b
Qaddouri, A.: Nonlinear shallowwater equations on the YinYang grid, Q. J. Roy. Meteor. Soc., 137, 810–818, https://doi.org/10.1002/qj.792, 2011. a
Ringler, T., Thuburn, J., Klemp, J., and Skamarock, W.: A unified approach to energy conservation and potential vorticity dynamics for arbitrarilystructured Cgrids, J. Comput. Phys., 229, 3065–3090, https://doi.org/10.1016/j.jcp.2009.12.007, 2010. a, b, c, d, e, f
Russell, G. L.: StepMountain Technique Applied to an Atmospheric CGrid Model, or How to Improve Precipitation near Mountains, Mon. Weather Rev., 135, 4060–4076, https://doi.org/10.1175/2007mwr2048.1, 2007. a, b
Russell, G. L. and Lerner, J. A.: A New FiniteDifferencing Scheme for the Tracer Transport Equation, J. Appl. Meteorol., 20, 1483–1498, https://doi.org/10.1175/15200450(1981)020<1483:anfdsf>2.0.co;2, 1981. a
Russell, G. L., Rind, D. H., and Jonas, J.: Symmetric Equations on the Surface of a Sphere as Used by Model GISS:IB, Zenodo, https://doi.org/10.5281/zenodo.1313736, 2018. a
Schmidt, G. A., Ruedy, R., Hansen, J. E., Aleinov, I., Bell, N., Bauer, M., Bauer, S., Cairns, B., Canuto, V., Cheng, Y., Genio, A. D., Faluvegi, G., Friend, A. D., Hall, T. M., Hu, Y., Kelley, M., Kiang, N. Y., Koch, D., Lacis, A. A., Lerner, J., Lo, K. K., Miller, R. L., Nazarenko, L., Oinas, V., Perlwitz, J., Perlwitz, J., Rind, D., Romanou, A., Russell, G. L., Sato, M., Shindell, D. T., Stone, P. H., Sun, S., Tausnev, N., Thresher, D., and Yao, M.S.: PresentDay Atmospheric Simulations Using GISS ModelE: Comparison to In Situ, Satellite, and Reanalysis Data, J. Climate, 19, 153–192, https://doi.org/10.1175/jcli3612.1, 2006. a, b
Stuhne, G. and Peltier, W.: New Icosahedral GridPoint Discretizations of the Shallow Water Equations on the Sphere, J. Comput. Phys., 148, 23–58, https://doi.org/10.1006/jcph.1998.6119, 1999. a, b, c, d, e
Sun, S. and Bleck, R.: Multicentury simulations with the coupled GISS–HYCOM climate model: control experiments, Clim. Dynam., 26, 407–428, https://doi.org/10.1007/s0038200500917, 2005. a
Swarztrauber, P. N.: The Vector Harmonic Transform Method for Solving Partial Differential Equations in Spherical Geometry, Mon. Weather Rev., 121, 3415–3437, https://doi.org/10.1175/15200493(1993)121<3415:tvhtmf>2.0.co;2, 1993. a, b
Swarztrauber, P. N.: Spectral Transform Methods for Solving the ShallowWater Equations on the Sphere, Mon. Weather Rev., 124, 730–744, https://doi.org/10.1175/15200493(1996)124<0730:stmfst>2.0.co;2, 1996. a
Temperton, C.: On Scalar and Vector Transform Methods for Global Spectral Models, Mon. Weather Rev., 119, 1303–1307, https://doi.org/10.1175/1520049311951303.1, 1991. a
Weller, H., Thuburn, J., and Cotter, C. J.: Computational Modes and Grid Imprinting on Five QuasiUniform Spherical C Grids, Mon. Weather Rev., 140, 2734–2755, https://doi.org/10.1175/mwrd1100193.1, 2012. a, b, c, d
Williamson, D. L., Drake, J. B., Hack, J. J., Jakob, R., and Swarztrauber, P. N.: A standard test set for numerical approximations to the shallow water equations in spherical geometry, J. Comput. Phys., 102, 211–224, https://doi.org/10.1016/s00219991(05)800166, 1992. a, b, c, d, e, f, g
 Abstract
 Introduction
 Symmetric mathematics
 Discrete implementation of symmetric equations
 Test case results
 Discussion and conclusions
 Code availability
 Appendix A: Equivalence between new and old forms for ∇h
 Appendix B: Closed form integration along an arc
 Author contributions
 Competing interests
 Acknowledgements
 References
 Abstract
 Introduction
 Symmetric mathematics
 Discrete implementation of symmetric equations
 Test case results
 Discussion and conclusions
 Code availability
 Appendix A: Equivalence between new and old forms for ∇h
 Appendix B: Closed form integration along an arc
 Author contributions
 Competing interests
 Acknowledgements
 References