Symmetric equations on the surface of a sphere as used by model GISS:IB

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 one-layer shallow water models on cubed-sphere 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.

1 Introduction

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 cubed-sphere 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 one-layer model. Other researchers have explored cubed-sphere grids: Adcroft et al. (2004) and Putman and Lin (2007). The advantages of using a cubed-sphere grid over a lat–long grid are that the two polar singularities on the lat–long grid are replaced by eight ill-behaved corners on the cubed-sphere 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 cubed-sphere 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 one-layer 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 two-component 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 vector-invariant (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 two-dimensional 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 one-layer shallow water equation models: one for a gnomonic cubed-sphere 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 three-dimensional 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³ 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 vector-invariant 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 non-horizontal 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³ 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 grid-matched 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 fourth-order schemes become first-order schemes …” The mass coordinate at a point, in a spatial one-dimensional model, is the number of kilograms between the origin and the point. Second-order, fourth-order, 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 second-order moments) inside each grid cell; these schemes are less sensitive to mass variations than are non-upstream schemes. The one-layer icosahedral model IB, discussed in detail later in this paper, uses a combination of linear upstream and second-order schemes for momentum advection but does not carry prognostic gradients. According to Weller et al. (2012), “an upwind-based interpolation of the potential vorticity controls the computational Rossby modes” in some one-layer 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 one-layer model is tested with various resolutions, it should include tests where such increased gradients exist. If a one-layer 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, cubed-sphere, and icosahedral one-layer 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³ on the fixed Earth and not a local coordinate system.

2 Symmetric mathematics

A sphere of radius 1 is centered at the origin in three-dimensional 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⁻¹), 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

$$\begin{array}{ll}\text{(2.1)}& {\displaystyle }\mathbf{\text{U}}& {\displaystyle }={\displaystyle \frac{\mathbf{\text{X}}\times \mathbf{\text{P}}}{|\mathbf{\text{X}}\times \mathbf{\text{P}}|}}={\displaystyle \frac{(\mathrm{0},-r,q)}{\sqrt{q^{\mathrm{2}}+r^{\mathrm{2}}}}};\text{(2.2)}& {\displaystyle }\mathbf{\text{V}}& {\displaystyle }={\displaystyle \frac{\mathbf{\text{Y}}\times \mathbf{\text{P}}}{|\mathbf{\text{Y}}\times \mathbf{\text{P}}|}}={\displaystyle \frac{(r,\mathrm{0},-p)}{\sqrt{r^{\mathrm{2}}+p^{\mathrm{2}}}}};\text{(2.3)}& {\displaystyle }\mathbf{\text{W}}& {\displaystyle }={\displaystyle \frac{\mathbf{\text{Z}}\times \mathbf{\text{P}}}{|\mathbf{\text{Z}}\times \mathbf{\text{P}}|}}={\displaystyle \frac{(-q,p,\mathrm{0})}{\sqrt{p^{\mathrm{2}}+q^{\mathrm{2}}}}},\mathrm{\text{which points eastward}};{\displaystyle }\mathbf{\text{N}}& {\displaystyle }=\mathbf{\text{P}}\times \mathbf{\text{W}}={\displaystyle \frac{(-rp,-qr,p^{\mathrm{2}}+q^{\mathrm{2}})}{\sqrt{p^{\mathrm{2}}+q^{\mathrm{2}}}}},\\ \text{(2.4)}& {\displaystyle }& {\displaystyle }\mathrm{\text{which points northward}}.\end{array}$$

S (m s⁻¹) is horizontal velocity at point P and is equal to u of Cartesian R³ 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⁻¹), uses the variables a, b, and c:

$$\begin{array}{ll}{\displaystyle }\mathbf{\text{A}}& {\displaystyle }=\left(u\sqrt{q^{\mathrm{2}}+r^{\mathrm{2}}},v\sqrt{r^{\mathrm{2}}+p^{\mathrm{2}}},w\sqrt{p^{\mathrm{2}}+q^{\mathrm{2}}}\right)=\\ \text{(2.5)}& {\displaystyle }& {\displaystyle }=(a,b,c).\end{array}$$

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.

$$\begin{array}{}\text{(2.6)}& \begin{array}{rl}\mathbf{\text{S}}& =(rv\sqrt{r^{\mathrm{2}}+p^{\mathrm{2}}}-qw\sqrt{p^{\mathrm{2}}+q^{\mathrm{2}}},pw\sqrt{p^{\mathrm{2}}+q^{\mathrm{2}}}\\ & -ru\sqrt{q^{\mathrm{2}}+r^{\mathrm{2}}},qu\sqrt{q^{\mathrm{2}}+r^{\mathrm{2}}}-pv\sqrt{r^{\mathrm{2}}+p^{\mathrm{2}}})=\\ & =(rb-qc,pc-ra,qa-pb).\end{array}\end{array}$$

The northern velocity component rotated 90^∘ from W is $n=\mathbf{\text{S}}\cdot \mathbf{\text{N}}=(qa-pb)/\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}}=(rb-qc)/\sqrt{q^{\mathrm{2}}+r^{\mathrm{2}}}$ and that rotated from V is $\mathbf{\text{S}}\cdot \mathbf{\text{P}}\times \mathbf{\text{V}}=(pc-ra)/\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:

$$\begin{array}{}\text{(2.7)}& \begin{array}{rl}(p,q,r)& =(\sin \mathit{\delta },\cos \mathit{\mu }\cos \mathit{\delta },\sin \mathit{\mu }\cos \mathit{\delta })=\\ & =(\sin \mathit{\nu }\cos \mathit{ϵ},\sin \mathit{ϵ},\cos \mathit{\nu }\cos \mathit{ϵ})=\\ & =(\cos \mathit{\lambda }\cos \mathit{\varphi },\sin \mathit{\lambda }\cos \mathit{\varphi },\sin \mathit{\varphi }).\end{array}\end{array}$$

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}}/\cos \mathit{\delta }$. A few useful derivatives are

$$\begin{array}{}\text{(2.8)}& {\displaystyle \frac{\partial \mathit{\lambda }}{\partial \mathit{\mu }}}& {\displaystyle }={\displaystyle \frac{\mathbf{\text{W}}\cdot \mathbf{\text{U}}\cos \mathit{\delta }}{\cos \mathit{\varphi }}}={\displaystyle \frac{-rp}{p^{\mathrm{2}}+q^{\mathrm{2}}}}=-\cos \mathit{\lambda }\tan \mathit{\varphi };\text{(2.9)}& {\displaystyle \frac{\partial \mathit{\varphi }}{\partial \mathit{\mu }}}& {\displaystyle }=\mathbf{\text{N}}\cdot \mathbf{\text{U}}\cos \mathit{\delta }={\displaystyle \frac{q}{\sqrt{p^{\mathrm{2}}+q^{\mathrm{2}}}}}=\sin \mathit{\lambda };\text{(2.10)}& {\displaystyle \frac{\partial \mathit{\mu }}{\partial \mathit{\lambda }}}& {\displaystyle }={\displaystyle \frac{\mathbf{\text{U}}\cdot \mathbf{\text{W}}\cos \mathit{\varphi }}{\cos \mathit{\delta }}}={\displaystyle \frac{-rp}{q^{\mathrm{2}}+r^{\mathrm{2}}}};\text{(2.11)}& {\displaystyle \frac{\partial \mathit{\mu }}{\partial \mathit{\varphi }}}& {\displaystyle }={\displaystyle \frac{\mathbf{\text{U}}\cdot \mathbf{\text{N}}}{\cos \mathit{\delta }}}={\displaystyle \frac{q}{(q^{\mathrm{2}}+r^{\mathrm{2}})\sqrt{p^{\mathrm{2}}+q^{\mathrm{2}}}}}.\end{array}$$

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⁻³) 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⁻²), the downward vertical acceleration due to gravity, and Ω (s⁻¹), the Earth's angular rotation rate, are assumed to be uniform. The field top geopotential Φ (m² s⁻²) is g(h+h_S).

The new symmetric equations to be presented here are applicable to many grid arrangements; one is the gnomonic cubed-sphere 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 one-layer shallow water equation model, labeled CSK, may be published later, but some results are presented in Sect. 4.

Figure 1 Upper diagram shows arrangement of velocity components on Face 1 of cubed-sphere grid centered around the North Pole. Lower diagram shows arrangement of velocity components at the intersection of Faces 1, 2, and 3. Although U, V, and W are defined everywhere except at their poles, only the component perpendicular to an edge is displayed.

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 B-grid 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 four-triangle wedge of an icosahedral grid.

Figure 2 Four triangles of an icosahedron are surrounded by thick, bold lines. Vertices of the bold triangles are centers of incomplete pentagonal primary cells. Edges of primary cells, for grid level 2, are indicated by thin solid lines. Dotted lines, between primary cell centers, indicate edges of triangular momentum cells. Primary cell corners are momentum cell centers.

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

$$\begin{array}{}\text{(2.12)}& \begin{array}{rl}\mathrm{\nabla }& =\left(\cos \mathit{\delta }{\displaystyle \frac{\partial }{\partial \mathit{\delta }}},\cos \mathit{ϵ}{\displaystyle \frac{\partial }{\partial \mathit{ϵ}}},\cos \mathit{\varphi }{\displaystyle \frac{\partial }{\partial \mathit{\varphi }}}\right)/R=\\ & =\left(r{\displaystyle \frac{\partial }{\partial \mathit{\nu }}}-q{\displaystyle \frac{\partial }{\partial \mathit{\lambda }}},p{\displaystyle \frac{\partial }{\partial \mathit{\lambda }}}-r{\displaystyle \frac{\partial }{\partial \mathit{\mu }}},q{\displaystyle \frac{\partial }{\partial \mathit{\mu }}}-p{\displaystyle \frac{\partial }{\partial \mathit{\nu }}}\right)/R\end{array},\end{array}$$

which is equivalent to the common two-dimensional del operator on the surface of a sphere,
${\mathrm{\nabla }}_{\text{R}}=(\mathbf{\text{W}}\partial /\partial \mathit{\lambda }+\mathbf{\text{N}}\cos \mathit{\varphi }\partial /\partial \mathit{\varphi })/R\cos \mathit{\varphi }$ (Williamson et al., 1992, Eq. 3). ∇_S or ∇_R can be applied to three space, and be equivalent to the three-dimensional 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

$$\begin{array}{ll}{\displaystyle }\mathrm{\nabla }h& {\displaystyle }=\left(\cos \mathit{\delta }{\displaystyle \frac{\partial h}{\partial \mathit{\delta }}},\cos \mathit{ϵ}{\displaystyle \frac{\partial h}{\partial \mathit{ϵ}}},\cos \mathit{\varphi }{\displaystyle \frac{\partial h}{\partial \mathit{\varphi }}}\right)/R=\\ {\displaystyle }& {\displaystyle }=\left(r{\displaystyle \frac{\partial h}{\partial \mathit{\nu }}}-q{\displaystyle \frac{\partial h}{\partial \mathit{\lambda }}},p{\displaystyle \frac{\partial h}{\partial \mathit{\lambda }}}-r{\displaystyle \frac{\partial h}{\partial \mathit{\mu }}},q{\displaystyle \frac{\partial h}{\partial \mathit{\mu }}}-p{\displaystyle \frac{\partial h}{\partial \mathit{\nu }}}\right)/R=\\ \text{(2.13)}& {\displaystyle }& {\displaystyle }=\left(\mathbf{\text{U}}\cos \mathit{\delta }{\displaystyle \frac{\partial h}{\partial \mathit{\mu }}}+\mathbf{\text{V}}\cos \mathit{ϵ}{\displaystyle \frac{\partial h}{\partial \mathit{\nu }}}+\mathbf{\text{W}}\cos \mathit{\varphi }{\displaystyle \frac{\partial h}{\partial \mathit{\lambda }}}\right)/R.\end{array}$$

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

$$\begin{array}{}\text{(2.14)}& \mathrm{\nabla }h=(\mathbf{\text{W}}\partial h/\partial \mathit{\lambda }+\mathbf{\text{N}}\cos \mathit{\varphi }\partial h/\partial \mathit{\varphi })/R\cos \mathit{\varphi }.\end{array}$$

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

$$\begin{array}{}\text{(2.15)}& \mathbf{\text{A}}=\mathbf{\text{P}}\times \mathrm{\nabla }h=\left({\displaystyle \frac{\partial h}{\partial \mathit{\mu }}},{\displaystyle \frac{\partial h}{\partial \mathit{\nu }}},{\displaystyle \frac{\partial h}{\partial \mathit{\lambda }}}\right)/R.\end{array}$$

The divergence of a horizontal vector $\mathbf{\text{D}}=(d,e,f)$ at point P is

$$\begin{array}{ll}\text{(2.16)}& {\displaystyle }\mathrm{\nabla }\cdot \mathbf{\text{D}}& {\displaystyle }=\left(\cos \mathit{\delta }{\displaystyle \frac{\partial d}{\partial \mathit{\delta }}}+\cos \mathit{ϵ}{\displaystyle \frac{\partial e}{\partial \mathit{ϵ}}}+\cos \mathit{\varphi }{\displaystyle \frac{\partial f}{\partial \mathit{\varphi }}}\right)/R={\displaystyle }& {\displaystyle }=\left[{\displaystyle \frac{\partial \left(d\cos \mathit{\delta }\right)}{\partial \mathit{\delta }}}+{\displaystyle \frac{\partial \left(e\cos \mathit{ϵ}\right)}{\partial \mathit{ϵ}}}+{\displaystyle \frac{\partial \left(f\cos \mathit{\varphi }\right)}{\partial \mathit{\varphi }}}\right]/R=\\ {\displaystyle }& {\displaystyle }=\left(r{\displaystyle \frac{\partial d}{\partial \mathit{\nu }}}-q{\displaystyle \frac{\partial d}{\partial \mathit{\lambda }}}+p{\displaystyle \frac{\partial e}{\partial \mathit{\lambda }}}-r{\displaystyle \frac{\partial e}{\partial \mathit{\mu }}}+q{\displaystyle \frac{\partial f}{\partial \mathit{\mu }}}-p{\displaystyle \frac{\partial f}{\partial \mathit{\nu }}}\right)/R=\\ {\displaystyle }& {\displaystyle }=\left[{\displaystyle \frac{\partial (qf-re)}{\partial \mathit{\mu }}}+{\displaystyle \frac{\partial (rd-pf)}{\partial \mathit{\nu }}}+{\displaystyle \frac{\partial (pe-qd)}{\partial \mathit{\lambda }}}\right]/R,\end{array}$$

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}}\cos \mathit{\varphi })/\partial \mathit{\varphi }]/R\cos \mathit{\varphi }$. The first two forms of Eq. (2.16) are equivalent because $\mathbf{\text{D}}\cdot \mathbf{\text{P}}=\mathrm{0}$ and

$$\begin{array}{ll}\text{(2.17)}& {\displaystyle }& {\displaystyle \frac{\partial \left(d\cos \mathit{\delta }\right)}{\partial \mathit{\delta }}}+{\displaystyle \frac{\partial \left(e\cos \mathit{ϵ}\right)}{\partial \mathit{ϵ}}}+{\displaystyle \frac{\partial \left(f\cos \mathit{\varphi }\right)}{\partial \mathit{\varphi }}}={\displaystyle }& {\displaystyle }=\cos \mathit{\delta }{\displaystyle \frac{\partial d}{\partial \mathit{\delta }}}-d\sin \mathit{\delta }+\cos \mathit{ϵ}{\displaystyle \frac{\partial e}{\partial \mathit{ϵ}}}-e\sin \mathit{ϵ}+\cos \mathit{\varphi }{\displaystyle \frac{\partial f}{\partial \mathit{\varphi }}}\\ {\displaystyle }& {\displaystyle }-f\sin \mathit{\varphi }=\cos \mathit{\delta }{\displaystyle \frac{\partial d}{\partial \mathit{\delta }}}+\cos \mathit{ϵ}{\displaystyle \frac{\partial e}{\partial \mathit{ϵ}}}+\cos \mathit{\varphi }{\displaystyle \frac{\partial f}{\partial \mathit{\varphi }}}-\mathbf{\text{D}}\cdot \mathbf{\text{P}}.\end{array}$$

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

$$\begin{array}{}\text{(2.18)}& \mathrm{\nabla }\cdot \mathbf{\text{S}}=\mathrm{\nabla }\cdot \mathbf{\text{A}}\times \mathbf{\text{P}}=\left({\displaystyle \frac{\partial a}{\partial \mathit{\mu }}}+{\displaystyle \frac{\partial b}{\partial \mathit{\nu }}}+{\displaystyle \frac{\partial c}{\partial \mathit{\lambda }}}\right)/R.\end{array}$$

Also,

$$\begin{array}{ll}{\displaystyle }\mathrm{\nabla }\cdot \mathbf{\text{A}}& {\displaystyle }=\mathrm{\nabla }\cdot \mathbf{\text{P}}\times \mathbf{\text{S}}=-\mathrm{\nabla }\cdot \mathbf{\text{S}}\times \mathbf{\text{P}}=\\ \text{(2.19)}& {\displaystyle }& {\displaystyle }-\mathrm{\nabla }\cdot \mathbf{\text{D}}\times \mathbf{\text{P}}=-\left({\displaystyle \frac{\partial d}{\partial \mathit{\mu }}}+{\displaystyle \frac{\partial e}{\partial \mathit{\nu }}}+{\displaystyle \frac{\partial f}{\partial \mathit{\lambda }}}\right)/R.\end{array}$$

Noting that ∇h is perpendicular to P and reasoning similar to Eq. (2.17), the Laplacian is

$$\begin{array}{ll}{\displaystyle }{\mathrm{\nabla }}^{\mathrm{2}}h& {\displaystyle }=\left({\cos }^{\mathrm{2}}\mathit{\delta }{\displaystyle \frac{{\partial }^{\mathrm{2}}h}{\partial {\mathit{\delta }}^{\mathrm{2}}}}+{\cos }^{\mathrm{2}}\mathit{ϵ}{\displaystyle \frac{{\partial }^{\mathrm{2}}h}{\partial {\mathit{ϵ}}^{\mathrm{2}}}}+{\cos }^{\mathrm{2}}\mathit{\varphi }{\displaystyle \frac{{\partial }^{\mathrm{2}}h}{\partial {\mathit{\varphi }}^{\mathrm{2}}}}\right)/R^{\mathrm{2}}=\\ \text{(2.20)}& {\displaystyle }& {\displaystyle }=\left({\displaystyle \frac{{\partial }^{\mathrm{2}}h}{\partial {\mathit{\mu }}^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}h}{\partial {\mathit{\nu }}^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}h}{\partial {\mathit{\lambda }}^{\mathrm{2}}}}\right)/R^{\mathrm{2}}.\end{array}$$

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}}+\cos \mathit{\varphi }\partial (\cos \mathit{\varphi }\partial h/\partial \mathit{\varphi })/\partial \mathit{\varphi }]/R^{\mathrm{2}}{\cos }^{\mathrm{2}}\mathit{\varphi }$.

The curl of a horizontal vector $\mathbf{\text{D}}=(d,e,f)$ is

$$\begin{array}{ll}{\displaystyle }\mathrm{\nabla }\times \mathbf{\text{D}}=& {\displaystyle }\left(\cos \mathit{ϵ}{\displaystyle \frac{\partial f}{\partial \mathit{ϵ}}}-\cos \mathit{\varphi }{\displaystyle \frac{\partial e}{\partial \mathit{\varphi }}},\cos \mathit{\varphi }{\displaystyle \frac{\partial d}{\partial \mathit{\varphi }}}-\cos \mathit{\delta }{\displaystyle \frac{\partial f}{\partial \mathit{\delta }}},\right.\\ \text{(2.21)}& {\displaystyle }& {\displaystyle }\left.\cos \mathit{\delta }{\displaystyle \frac{\partial e}{\partial \mathit{\delta }}}-\cos \mathit{ϵ}{\displaystyle \frac{\partial d}{\partial \mathit{ϵ}}}\right)/R.\end{array}$$

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

$$\begin{array}{}\text{(2.22)}& \begin{array}{rl}\mathbf{\text{P}}\cdot \mathrm{\nabla }\times \mathbf{\text{D}}& =\left[p\left(\cos \mathit{ϵ}{\displaystyle \frac{\partial f}{\partial \mathit{ϵ}}}-\cos \mathit{\varphi }{\displaystyle \frac{\partial e}{\partial \mathit{\varphi }}}\right)\right.\\ & +q\left(\cos \mathit{\varphi }{\displaystyle \frac{\partial d}{\partial \mathit{\varphi }}}-\cos \mathit{\delta }{\displaystyle \frac{\partial f}{\partial \mathit{\delta }}}\right)\\ & \left.+r\left(\cos \mathit{\delta }{\displaystyle \frac{\partial e}{\partial \mathit{\delta }}}-\cos \mathit{ϵ}{\displaystyle \frac{\partial d}{\partial \mathit{ϵ}}}\right)\right]/R\\ & =\left({\displaystyle \frac{\partial d}{\partial \mathit{\mu }}}+{\displaystyle \frac{\partial e}{\partial \mathit{\nu }}}+{\displaystyle \frac{\partial f}{\partial \mathit{\lambda }}}\right)/R=-\mathrm{\nabla }\cdot \mathbf{\text{P}}\times \mathbf{\text{D}},\end{array}\end{array}$$

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}}\cos \mathit{\varphi })/\partial \mathit{\varphi }]/R\cos \mathit{\varphi }$.
If D = S, then the vertical component of the curl of S can also be written as

$$\begin{array}{ll}{\displaystyle }\mathbf{\text{P}}& {\displaystyle }\cdot \mathrm{\nabla }\times \mathbf{\text{S}}=\mathbf{\text{P}}\cdot \mathrm{\nabla }\times (\mathbf{\text{A}}\times \mathbf{\text{P}})=-\mathrm{\nabla }\cdot \mathbf{\text{A}}=\\ \text{(2.23)}& {\displaystyle }& {\displaystyle }=-\left[{\displaystyle \frac{\partial \left(a\cos \mathit{\delta }\right)}{\partial \mathit{\delta }}}+{\displaystyle \frac{\partial \left(b\cos \mathit{ϵ}\right)}{\partial \mathit{ϵ}}}+{\displaystyle \frac{\partial \left(c\cos \mathit{\varphi }\right)}{\partial \mathit{\varphi }}}\right]/R.\end{array}$$

The upward vertical component of relative vorticity is

$$\begin{array}{}\text{(2.24)}& \mathit{\zeta }=\mathbf{\text{P}}\cdot \mathrm{\nabla }\times \mathbf{\text{S}}=\mathbf{\text{P}}\cdot \mathrm{\nabla }\times (\mathbf{\text{A}}\times \mathbf{\text{P}}),\end{array}$$

which is equivalent to the common form on the surface of a sphere: $[\partial n/\partial \mathit{\lambda }-\partial (w\cos \mathit{\varphi })/\partial \mathit{\varphi }]/R\cos \mathit{\varphi }$.

If G and H are differentiable vectors in three space, then a well-known identity using the Cartesian del operator is

$$\begin{array}{}\text{(2.25)}& \mathbf{\text{G}}\cdot {\mathrm{\nabla }}_C\times \mathbf{\text{H}}=\mathbf{\text{H}}\cdot {\mathrm{\nabla }}_C\times \mathbf{\text{G}}-{\mathrm{\nabla }}_C\cdot \mathbf{\text{G}}\times \mathbf{\text{H}}.\end{array}$$

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

$$\begin{array}{}\text{(2.26)}& \mathbf{\text{G}}\cdot {\mathrm{\nabla }}_C\times (\mathbf{\text{H}}\times \mathbf{\text{G}})=-{\mathrm{\nabla }}_C\cdot \mathbf{\text{G}}\times (\mathbf{\text{H}}\times \mathbf{\text{G}})=-{\mathrm{\nabla }}_C\cdot \mathbf{\text{H}}.\end{array}$$

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:

$$\begin{array}{ll}\text{(2.27)}& {\displaystyle \frac{\partial h}{\partial t}}& {\displaystyle }=-\mathrm{\nabla }\cdot \left(h\mathbf{\text{S}}\right)=-\mathrm{\nabla }\cdot (h\mathbf{\text{A}}\times \mathbf{\text{P}})={\displaystyle }& {\displaystyle }=-\left[{\displaystyle \frac{\partial \left(ha\right)}{\partial \mathit{\mu }}}+{\displaystyle \frac{\partial \left(hb\right)}{\partial \mathit{\nu }}}+{\displaystyle \frac{\partial \left(hc\right)}{\partial \mathit{\lambda }}}\right]/R=\\ {\displaystyle }& {\displaystyle }=-\left[\cos \mathit{\delta }{\displaystyle \frac{\partial \left(hu\right)}{\partial \mathit{\mu }}}+\cos \mathit{ϵ}{\displaystyle \frac{\partial \left(hv\right)}{\partial \mathit{\nu }}}+\cos \mathit{\varphi }{\displaystyle \frac{\partial \left(hw\right)}{\partial \mathit{\lambda }}}\right]/R.\end{array}$$

This symmetric form is equivalent to the common equation:

$$\begin{array}{}\text{(2.28)}& {\displaystyle \frac{\partial h}{\partial t}}+\left[{\displaystyle \frac{\partial \left(hw\right)}{\partial \mathit{\lambda }}}+{\displaystyle \frac{\partial \left(hn\cos \mathit{\varphi }\right)}{\partial \mathit{\varphi }}}\right]/R\cos \mathit{\varphi }=\mathrm{0}.\end{array}$$

The three-component advective form for specific angular momentum is

$$\begin{array}{ll}{\displaystyle \frac{\partial \mathbf{\text{A}}}{\partial t}}& {\displaystyle }=-\left(a{\displaystyle \frac{\partial \mathbf{\text{A}}}{\partial \mathit{\mu }}}+b{\displaystyle \frac{\partial \mathbf{\text{A}}}{\partial \mathit{\nu }}}+c{\displaystyle \frac{\partial \mathbf{\text{A}}}{\partial \mathit{\lambda }}}\right)/R+\\ \text{(2.29)}& {\displaystyle }& {\displaystyle }+f\mathbf{\text{A}}\times \mathbf{\text{P}}-\mathbf{\text{P}}\times \mathrm{\nabla }\mathrm{\Phi },\end{array}$$

where f (s⁻¹) is the Coriolis parameter and Φ (m² s⁻²) is the fluid top geopotential:

$$\begin{array}{}\text{(2.30)}& {\displaystyle }& {\displaystyle }f=\mathrm{2}\mathrm{\Omega }\sin \mathit{\varphi };\text{(2.31)}& {\displaystyle }& {\displaystyle }\mathrm{\Phi }=g(h+h_{\text{S}}).\end{array}$$

Replacing h with hA in Eqs. (2.27) and (2.28) shows the differential form for momentum advection in flux form:

$$\begin{array}{}\text{(2.32)}& \begin{array}{rl}& {\displaystyle \frac{\partial \left(h\mathbf{\text{A}}\right)}{\partial t}}+\left[{\displaystyle \frac{\partial \left(ha\mathbf{\text{A}}\right)}{\partial \mathit{\mu }}}+{\displaystyle \frac{\partial \left(hb\mathbf{\text{A}}\right)}{\partial \mathit{\nu }}}+{\displaystyle \frac{\partial \left(hc\mathbf{\text{A}}\right)}{\partial \mathit{\lambda }}}\right]/R=\\ & ={\displaystyle \frac{\partial \left(h\mathbf{\text{A}}\right)}{\partial t}}+\left[{\displaystyle \frac{\partial \left(hw\mathbf{\text{A}}\right)}{\partial \mathit{\lambda }}}+{\displaystyle \frac{\partial \left(hn\mathbf{\text{A}}\cos \mathit{\varphi }\right)}{\partial \mathit{\varphi }}}\right]/R\cos \mathit{\varphi }=\\ & =h\left[{\displaystyle \frac{\partial \mathbf{\text{A}}}{\partial t}}+\left(w{\displaystyle \frac{\partial \mathbf{\text{A}}}{\partial \mathit{\lambda }}}+n\cos \mathit{\varphi }{\displaystyle \frac{\partial \mathbf{\text{A}}}{\partial \mathit{\varphi }}}\right)/R\cos \mathit{\varphi }\right]+\\ & +\mathbf{\text{A}}\left\{{\displaystyle \frac{\partial h}{\partial t}}+\left[{\displaystyle \frac{\partial \left(hw\right)}{\partial \mathit{\lambda }}}+{\displaystyle \frac{\partial \left(hn\cos \mathit{\varphi }\right)}{\partial \mathit{\varphi }}}\right]/R\cos \mathit{\varphi }\right\}.\end{array}\end{array}$$

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

$$\begin{array}{ll}\text{(2.33)}& {\displaystyle }& {\displaystyle }\left\{{\displaystyle \frac{\partial \left(w\cos \mathit{\varphi }\right)}{\partial t}}+\right.{\displaystyle }& {\displaystyle }\left.+\left[w{\displaystyle \frac{\partial \left(w\cos \mathit{\varphi }\right)}{\partial \mathit{\lambda }}}+n\cos \mathit{\varphi }{\displaystyle \frac{\partial \left(w\cos \mathit{\varphi }\right)}{\partial \mathit{\varphi }}}\right]/R\cos \mathit{\varphi }\right\}/\cos \mathit{\varphi }=\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\partial w}{\partial t}}+\left(w{\displaystyle \frac{\partial w}{\partial \mathit{\lambda }}}+n\cos \mathit{\varphi }{\displaystyle \frac{\partial w}{\partial \mathit{\varphi }}}\right)/R\cos \mathit{\varphi }-nw\tan \mathit{\varphi }/R\end{array}$$

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 vector-invariant form are shown at https://aom.giss.nasa.gov. Model IB does not use this form nor Eq. (2.33).

3 Discrete implementation of symmetric equations

3.1 Alignment of velocity or specific angular momentum

Alignment of the three momentum components at P means that $\mathbf{\text{A}}=(u\cos \mathit{\delta },v\cos \mathit{ϵ},w\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.

$$\begin{array}{ll}{\displaystyle }s& {\displaystyle }={\cos }^{\mathrm{2}}\mathit{\delta }({\mathbf{\text{S}}}_{\text{NEW}}\cdot \mathbf{\text{U}}-u{)}^{\mathrm{2}}+{\cos }^{\mathrm{2}}\mathit{ϵ}({\mathbf{\text{S}}}_{\text{NEW}}\cdot \mathbf{\text{V}}-v{)}^{\mathrm{2}}\\ \text{(3.1)}& {\displaystyle }& {\displaystyle }+{\cos }^{\mathrm{2}}\mathit{\varphi }({\mathbf{\text{S}}}_{\text{NEW}}\cdot \mathbf{\text{W}}-w{)}^{\mathrm{2}}\end{array}$$

If u, v, and w were already aligned, then t=0, where

$$\begin{array}{}\text{(3.2)}& t=pu\cos \mathit{\delta }+qv\cos \mathit{ϵ}+rw\cos \mathit{\varphi }.\end{array}$$

When s is minimized, alignment of distorted components u, v, and w produces

$$\begin{array}{}\text{(3.3)}& {\displaystyle }{u}_{\text{NEW}}& {\displaystyle }={\mathbf{\text{S}}}_{\text{NEW}}\cdot \mathbf{\text{U}}=u-pt/\cos \mathit{\delta };\text{(3.4)}& {\displaystyle }{v}_{\text{NEW}}& {\displaystyle }={\mathbf{\text{S}}}_{\text{NEW}}\cdot \mathbf{\text{V}}=v-qt/\cos \mathit{ϵ};\text{(3.5)}& {\displaystyle }{w}_{\text{NEW}}& {\displaystyle }={\mathbf{\text{S}}}_{\text{NEW}}\cdot \mathbf{\text{W}}=w-rt/\cos \mathit{\varphi }.\end{array}$$

Performing this analysis with specific angular momentum components yields a=ucos δ, $t=pa+qb+rc=\mathbf{\text{P}}\cdot \mathbf{\text{A}}$, and

$$\begin{array}{}\text{(3.6)}& {\mathbf{\text{A}}}_{\text{NEW}}=\mathbf{\text{A}}-\mathbf{\text{P}}(\mathbf{\text{P}}\cdot \mathbf{\text{A}}).\end{array}$$

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 cubed-sphere 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 Φ:

$$\begin{array}{}\text{(3.7)}& \mathrm{\Delta }\mathbf{\text{S}}=-\mathrm{\Delta }t\mathrm{\nabla }\mathrm{\Phi }=-g\mathrm{\Delta }t\mathrm{\nabla }(h+h_{\text{S}}),\end{array}$$

where Δt is the time step. Change of specific angular momentum, using Eq. (2.15), is

$$\begin{array}{}\text{(3.8)}& \mathrm{\Delta }\mathbf{\text{A}}=-\mathrm{\Delta }t\mathbf{\text{P}}\times \mathrm{\nabla }\mathrm{\Phi }=-\mathrm{\Delta }t\left({\displaystyle \frac{\partial \mathrm{\Phi }}{\partial \mathit{\mu }}},{\displaystyle \frac{\partial \mathrm{\Phi }}{\partial \mathit{\nu }}},{\displaystyle \frac{\partial \mathrm{\Phi }}{\partial \mathit{\lambda }}}\right)/R.\end{array}$$

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 cubed-sphere models and later for icosahedral models, like IB.

Acceleration of velocity averaged over a primary cell of a cubed-sphere 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,

$$\begin{array}{ll}{\displaystyle }\mathrm{\Delta }c=& {\displaystyle }\mathrm{\Delta }w\cos \mathit{\varphi }=-\mathrm{\Delta }t\cos \mathit{\varphi }{\displaystyle \frac{\partial \mathrm{\Phi }}{\partial {\mathit{\chi }}_W}}=-\mathrm{\Delta }t{\displaystyle \frac{\partial \left(\mathrm{\Phi }\cos \mathit{\varphi }\right)}{\partial {\mathit{\chi }}_W}}=\\ {\displaystyle }=& {\displaystyle }-\mathrm{\Delta }t\left[\right(\mathrm{\Phi }L\mathbf{\text{W}}\cdot \mathbf{\text{V}}\cos \mathit{\varphi }{)}_{\mathrm{RHT}}-(\mathrm{\Phi }L\mathbf{\text{W}}\cdot \mathbf{\text{V}}\cos \mathit{\varphi }{)}_{\mathrm{LEF}}-\\ \text{(3.9)}& {\displaystyle }& {\displaystyle }-(\mathrm{\Phi }L\mathbf{\text{W}}\cdot \mathbf{\text{U}}\cos \mathit{\varphi }{)}_{\mathrm{TOP}}+(\mathrm{\Phi }L\mathbf{\text{W}}\cdot \mathbf{\text{U}}\cos \mathit{\varphi }{)}_{\mathrm{BOT}}]/K,\end{array}$$

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

$$\begin{array}{ll}{\displaystyle }\underset{{\mathit{ϵ}}_{I,J-\mathrm{1}}}{\overset{{\mathit{ϵ}}_{I,J}}{\int }}R\mathbf{\text{W}}\cdot \mathbf{\text{V}}\cos \mathit{\varphi }d\mathit{ϵ}& {\displaystyle }=-R\cos \mathit{\nu }\underset{{\mathit{ϵ}}_{I,J-\mathrm{1}}}{\overset{{\mathit{ϵ}}_{I,J}}{\int }}\sin \mathit{ϵ}d\mathit{ϵ}=\\ \text{(3.10)}& {\displaystyle }& {\displaystyle }=R(r_{I,J}-r_{I,J-\mathrm{1}}).\end{array}$$

Similarly, the top edge [$I:I-\mathrm{1},J$] is integrated over arc length as

$$\begin{array}{ll}{\displaystyle }\underset{{\mathit{\delta }}_{I,J}}{\overset{{\mathit{\delta }}_{I-\mathrm{1},J}}{\int }}R\mathbf{\text{W}}\cdot \mathbf{\text{U}}\cos \mathit{\varphi }d\mathit{\delta }& {\displaystyle }=-R\sin \mathit{\mu }\underset{{\mathit{\delta }}_{I,J}}{\overset{{\mathit{\delta }}_{I-\mathrm{1},J}}{\int }}\sin \mathit{\delta }d\mathit{\delta }=\\ \text{(3.11)}& {\displaystyle }& {\displaystyle }=R(r_{I-\mathrm{1},J}-r_{I,J}).\end{array}$$

Integrating counterclockwise around the grid cell, similar formulas for different components yield

$$\begin{array}{}\text{(3.12)}& \begin{array}{rl}\mathrm{\Delta }\mathbf{\text{A}}=-R\mathrm{\Delta }t[& ({\mathbf{\text{P}}}_{I,J-\mathrm{1}}-{\mathbf{\text{P}}}_{I-\mathrm{1},J-\mathrm{1}}){\mathrm{\Phi }}_{\mathrm{BOT}}+\\ +& ({\mathbf{\text{P}}}_{I,J}-{\mathbf{\text{P}}}_{I,J-\mathrm{1}}){\mathrm{\Phi }}_{\mathrm{RHT}}+\\ +& ({\mathbf{\text{P}}}_{I-\mathrm{1},J}-{\mathbf{\text{P}}}_{I,J}){\mathrm{\Phi }}_{\mathrm{TOP}}+\\ +& ({\mathbf{\text{P}}}_{I-\mathrm{1},J-\mathrm{1}}-{\mathbf{\text{P}}}_{I-\mathrm{1},J}){\mathrm{\Phi }}_{\mathrm{LEF}}]/K.\end{array}\end{array}$$

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 cubed-sphere edges. To check this, any arc from P₁ to P₂ 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\sin \mathit{\eta }-j\cos \mathit{\xi }\cos \mathit{\eta },i\cos \mathit{\xi }\cos \mathit{\eta }+j\sin \mathit{\eta },\sin \mathit{\xi }\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}}=-\sin \mathit{\xi }\sin \mathit{\eta }/\cos \mathit{\varphi }$.

$$\begin{array}{}\text{(3.13)}& \underset{{\mathit{\eta }}_{\mathrm{1}}}{\overset{{\mathit{\eta }}_{\mathrm{2}}}{\int }}R\mathbf{\text{W}}\cdot \mathbf{\text{F}}\cos \mathit{\varphi }d\mathit{\eta }=-R\sin \mathit{\xi }\underset{{\mathit{\eta }}_{\mathrm{1}}}{\overset{{\mathit{\eta }}_{\mathrm{2}}}{\int }}\sin \mathit{\eta }d\mathit{\eta }=R(r_{\mathrm{2}}-r_{\mathrm{1}})\end{array}$$

Renaming the vertices of a spherical polygon with N arcs P₀, P₁, ... P_N=P₀ counted counterclockwise around the cell, and generalizing Eq. (3.12), yields

$$\begin{array}{}\text{(3.14)}& \mathrm{\Delta }\mathbf{\text{A}}=R\mathrm{\Delta }t\sum {\mathrm{\Phi }}_{n-\mathrm{1}/\mathrm{2}}({\mathbf{\text{P}}}_{n-\mathrm{1}}-{\mathbf{\text{P}}}_n)/K.\end{array}$$

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).

$$\begin{array}{}\text{(3.15)}& K\mathrm{\Delta }h=-\mathrm{\Delta }t\sum M_{\text{P}},\end{array}$$

where M_P (kg s⁻¹) 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 η₁ to η₂ is computed as

$$\begin{array}{}\text{(3.16)}& \begin{array}{rl}{M}_{\text{P}}& =\underset{{\mathit{\eta }}_{\mathrm{1}}}{\overset{{\mathit{\eta }}_{\mathrm{2}}}{\int }}Rh\mathbf{\text{S}}\cdot \mathbf{\text{F}}d\mathit{\eta }=R\underset{{\mathit{\eta }}_{\mathrm{1}}}{\overset{{\mathit{\eta }}_{\mathrm{2}}}{\int }}h(\mathbf{\text{P}}\times \mathbf{\text{S}})\cdot (\mathbf{\text{P}}\times \mathbf{\text{F}})d\mathit{\eta }=\\ & =R\underset{{\mathit{\eta }}_{\mathrm{1}}}{\overset{{\mathit{\eta }}_{\mathrm{2}}}{\int }}h\mathbf{\text{A}}\cdot \mathbf{\text{E}}d\mathit{\eta }=R\underset{{\mathit{\eta }}_{\mathrm{1}}}{\overset{{\mathit{\eta }}_{\mathrm{2}}}{\int }}h\mathbf{\text{A}}\cdot d\mathbf{\text{P}},\end{array}\end{array}$$

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

$$\begin{array}{}\text{(3.17)}& \mathrm{\Delta }h=R\mathrm{\Delta }t\sum h_{n-\mathrm{1}/\mathrm{2}}{\mathbf{\text{A}}}_{n-\mathrm{1}/\mathrm{2}}\cdot ({\mathbf{\text{P}}}_{n-\mathrm{1}}-{\mathbf{\text{P}}}_n)/K.\end{array}$$

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⁻²) 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.

$$\begin{array}{ll}{\displaystyle }{\mathbf{\text{V}}}_T& {\displaystyle }=\int Rh\mathbf{\text{A}}(\mathbf{\text{S}}\cdot \mathbf{\text{F}})d\mathit{\eta }=R\int h\mathbf{\text{A}}\left[\right(\mathbf{\text{P}}\times \mathbf{\text{S}})\cdot (\mathbf{\text{P}}\times \mathbf{\text{F}}\left)\right]d\mathit{\eta }=\\ {\displaystyle }& {\displaystyle }=R\int h\mathbf{\text{A}}(\mathbf{\text{A}}\cdot \mathbf{\text{E}})d\mathit{\eta }=R\int h\mathbf{\text{A}}(\mathbf{\text{A}}\cdot d\mathbf{\text{P}})=\\ \text{(3.18)}& {\displaystyle }& {\displaystyle }=\mathbf{\text{A}}{M}_T\end{array}$$

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

$$\begin{array}{ll}{\displaystyle }\mathrm{\Delta }\left(h\mathbf{\text{A}}\right)=R\mathrm{\Delta }t\sum & {\displaystyle }{h}_{n-\mathrm{1}/\mathrm{2}}{\mathbf{\text{A}}}_{n-\mathrm{1}/\mathrm{2}}\\ \text{(3.19)}& {\displaystyle }& {\displaystyle }[{\mathbf{\text{A}}}_{n-\mathrm{1}/\mathrm{2}}\cdot ({\mathbf{\text{P}}}_{n-\mathrm{1}}-{\mathbf{\text{P}}}_n\left)\right]/K.\end{array}$$

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^{⟂}$, 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.

Figure 3 C points are primary cell centers or triangular momentum cell corners. A and B points are primary cell corners or momentum cell centers. Quadrilaterals are the intersection area between primary and momentum cells. Arc intersection D points are overwritten by M_P in the diagram. D is halfway between C points but not halfway between A and B. Mass flux M_P along a primary cell edge is computed and not changed; it is partitioned to each quadrilateral edge proportional to arc length. Change in mass of a primary cell by advection is computed by summing M_P along the cell's edges. This total change is partitioned into each quadrilateral of the cell proportional to area. The initial value of mass flux M_T is computed along each momentum cell edge; half of M_T is attributed to each primary cell through which it passes. In each primary cell, the half M_T values are adjusted minimally so that the sum of the mass fluxes into each quadrilateral matches the quadrilateral's expected partitioned change. Each final M_T is equal to the sum of the two adjusted “half” M_T values.

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, $(pc-ra)/\cos \mathit{ϵ}$ is the velocity component that aligns with P×V. The Coriolis acceleration (m s⁻²) acting on v is $\mathrm{2}\mathrm{\Omega }\sin \mathit{\varphi }(pc-ra)/\cos \mathit{ϵ}$ and the acceleration acting on b is 2Ωsin ϕ(pc−ra). Each specific angular momentum component is accelerated by its local components:

$$\begin{array}{}\text{(3.20)}& \mathrm{\Delta }\mathbf{\text{A}}=\mathrm{2}\mathrm{\Omega }\mathrm{\Delta }t\mathbf{\text{A}}\times \mathbf{\text{P}}\sin \mathit{\varphi }.\end{array}$$

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.

Table 1 Properties of the raw, tweaked, and centroid grids; numerical columns are for grid levels 4 through 8. ArcA is the arc between adjacent primary cell corners or momentum cell centers. ArcC is the arc between nearby primary cell centers or momentum cell corners. Point D is intersection of ArcA and ArcC. The ideal number is 1 for first three properties and 0 for last two properties. Radius is the square root of cell area divided by π.

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 leap-frog 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 two-point second-order 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 step-mountain atmosphere.

For model IB, the unadjusted M_T is computed using second-order 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 second-order and half linear upstream advection. Unlike IK and the unfiltered lat–long B-grid 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.

4 Test case results

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 second-order B-grid (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 second-order C-grid (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 two-point interpolations of variables than do the other models because of its parallelogram-shaped 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 second-order 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)

Table 2 Abbreviations, number of primary cells at 1^∘ resolution, and model description.

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 second-order model, the mass and velocity errors should decrease by a factor of 4 when doubling the horizontal resolution. Figure 4 shows the area-weighted root mean square l₂ 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 4 Mass field error l₂ as function of grid level for solid body rotation initial conditions without bottom topography on day 5. Black line shows ideal mass error: factor of 4 error reduction for doubling horizontal resolution. Models represented are IBR, IBC, IKT, CSK, LLB, and LLC. For last three models, grid level is computed as ${\log }_{\mathrm{4}}\left[\right(N-\mathrm{2})/\mathrm{10}]$, where N is the number of primary grid cells.

Figure 5 shows the l₁ 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 C-grid model of Heikes and Randall (1995) at 4^∘ resolution show occasional instabilities that do not affect the long-term 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₂ norms in time. The Arakawa B-grid and C-grid schemes are designed to approximately conserve total energy.

Figure 5 Mass field errors l₁ and l_∞ as function of time for solid body rotation initial conditions without bottom topography. Models represented are IBR, IKT, CSK, and LLC at 2^∘ resolution.

Figure 6 Horizontal plots of h for models IBR, IBC, IKT, CSK, LLB, and LLC at 2^∘ resolution after 45 days of integration starting from Rossby–Haurwitz wave 3 conditions.

4.2 Rossby–Haurwitz wave 3 (RH3)

For cubed-sphere 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 Relative change (%) of specific kinetic energy for wave numbers 0 and 3 and total energy for various resolutions and models IBR, IBT, IBC, IKT, CSK, LLB, LLC, and National Center for Atmospheric Research spectral transform model at T42 resolution (Hack and Jakob, 1992) after 40 days of integration. Initial conditions are Rossby–Haurwitz wave 3. “Diverge” means the simulation diverged.

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 cubed-sphere output was interpolated to a high-resolution 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 high-resolution 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⁻¹ 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 leap-frog 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 one-third 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 coarse-resolution models, its numerical reduction coincides with the washing out of highs and lows of the height field, as shown in the videos discussed subsequently.

Table 4 For each model and approximate horizontal resolution, the largest dynamical time step, half of leap-frog time step, in seconds for which the model survives for 50 days without major instabilities, or after what day the model diverged (in parentheses), for solid body rotation initial conditions with Earth's bottom topography.

Table 5 Percentage change in kinetic energy for different models and resolutions after 50 days of integration from solid body rotation initial conditions with Earth's bottom topography. “Diverge” means the simulation diverged.

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 non-displayed 0.25^∘ symmetric models after 50 days; IBR and IKT are more similar to each other than to CSK. Non-displayed LLB diverges at 0.5^∘ resolution, but at 1^∘, it is similar to IBR.

Figure 7 Horizontal plots of h+h_S for models IBR, IKT, and CSK and resolutions 2 and 0.5^∘ after 34 days of integration starting from 50 m s⁻¹ eastward velocity multiplied by cosine of latitude, global mean height field top of 10 000 m but latitudinal distribution, and Earth's fixed bottom topography (Test Case 3 SBRZ).

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 multi-cell 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 checker-board 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.

5 Discussion and conclusions

This paper presents symmetric calculus operators on the surface of a sphere that are projected from three-dimensional 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.

1. 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 C-grid 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.

2. 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”).

3. 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 higher-resolution 0.5^∘ models than are 2^∘ IKT and CSK (Fig. 7). The icosahedral B-grid models, IBR and IBK, are best.

The SBRZ test is the most difficult but the most realistic test. Given that these one-layer models will be the basis for multi-layer climate models with realistic topography, IB models must be considered the best overall. Parallelogram-shaped grid cells in CSK with non-perpendicular 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 one-layer 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 fourth-order differencing or other methods, but in ocean domains or step-mountain 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 parallelogram-shaped 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.

Appendix A: Equivalence between new and old forms for ∇h

$$\begin{array}{}\text{(A1)}& {\displaystyle }\mathbf{\text{W}}& {\displaystyle }=(-q,p,\mathrm{0})/\sqrt{p^{\mathrm{2}}+q^{\mathrm{2}}}=(-\sin \mathit{\lambda },\cos \mathit{\lambda },\mathrm{0})\text{(A2)}& {\displaystyle }\mathbf{\text{N}}& {\displaystyle }=(-rp,-qr,p^{\mathrm{2}}+q^{\mathrm{2}})/\sqrt{p^{\mathrm{2}}+q^{\mathrm{2}}}=\text{(A3)}& {\displaystyle }& {\displaystyle }=(-\cos \mathit{\lambda }\sin \mathit{\varphi },-\sin \mathit{\lambda }\sin \mathit{\varphi },\cos \mathit{\varphi })\text{(A4)}& {\displaystyle }\partial \mathit{\lambda }/\partial \mathit{\delta }& {\displaystyle }=\mathbf{\text{W}}\cdot \mathbf{\text{P}}\times \mathbf{\text{U}}/\cos \mathit{\varphi }=-\sin \mathit{\lambda }/\cos \mathit{\varphi }\cos \mathit{\delta }\text{(A5)}& {\displaystyle }\partial \mathit{\lambda }/\partial \mathit{ϵ}& {\displaystyle }=\mathbf{\text{W}}\cdot \mathbf{\text{P}}\times \mathbf{\text{V}}/\cos \mathit{\varphi }=\cos \mathit{\lambda }/\cos \mathit{\varphi }\cos \mathit{ϵ}\text{(A6)}& {\displaystyle }\partial \mathit{\varphi }/\partial \mathit{\delta }& {\displaystyle }=\mathbf{\text{N}}\cdot \mathbf{\text{P}}\times \mathbf{\text{U}}=\mathbf{\text{W}}\cdot \mathbf{\text{U}}=-\cos \mathit{\lambda }\sin \mathit{\varphi }/\cos \mathit{\delta }\text{(A7)}& {\displaystyle }\partial \mathit{\varphi }/\partial \mathit{ϵ}& {\displaystyle }=\mathbf{\text{N}}\cdot \mathbf{\text{P}}\times \mathbf{\text{V}}=\mathbf{\text{W}}\cdot \mathbf{\text{V}}=-\sin \mathit{\lambda }\sin \mathit{\varphi }/\cos \mathit{ϵ}\end{array}$$

$$\begin{array}{}\text{(A8)}& \begin{array}{rl}\mathrm{\nabla }h=& \left(\cos \mathit{\delta }{\displaystyle \frac{\partial h}{\partial \mathit{\delta }}},\cos \mathit{ϵ}{\displaystyle \frac{\partial h}{\partial \mathit{ϵ}}},\cos \mathit{\varphi }{\displaystyle \frac{\partial h}{\partial \mathit{\varphi }}}\right)/R=\\ =& \left[\cos \mathit{\delta }\left({\displaystyle \frac{\partial \mathit{\lambda }}{\partial \mathit{\delta }}}{\displaystyle \frac{\partial h}{\partial \mathit{\lambda }}}+{\displaystyle \frac{\partial \mathit{\varphi }}{\partial \mathit{\delta }}}{\displaystyle \frac{\partial h}{\partial \mathit{\varphi }}}\right),\right.\\ & \left.\cos \mathit{ϵ}\left({\displaystyle \frac{\partial \mathit{\lambda }}{\partial \mathit{ϵ}}}{\displaystyle \frac{\partial h}{\partial \mathit{\lambda }}}+{\displaystyle \frac{\partial \mathit{\varphi }}{\partial \mathit{ϵ}}}{\displaystyle \frac{\partial h}{\partial \mathit{\varphi }}}\right),\cos \mathit{\varphi }{\displaystyle \frac{\partial h}{\partial \mathit{\varphi }}}\right]/R=\\ =& \left(-\sin \mathit{\lambda }{\displaystyle \frac{\partial h}{\partial \mathit{\lambda }}}/\cos \mathit{\varphi }-\cos \mathit{\lambda }\sin \mathit{\varphi }{\displaystyle \frac{\partial h}{\partial \mathit{\varphi }}},\right.\\ & \left.\cos \mathit{\lambda }{\displaystyle \frac{\partial h}{\partial \mathit{\lambda }}}/\cos \mathit{\varphi }-\sin \mathit{\lambda }\sin \mathit{\varphi }{\displaystyle \frac{\partial h}{\partial \mathit{\varphi }}},\cos \mathit{\varphi }{\displaystyle \frac{\partial h}{\partial \mathit{\varphi }}}\right)/R=\\ =& \left(\mathbf{\text{W}}{\displaystyle \frac{\partial h}{\partial \mathit{\lambda }}}/\cos \mathit{\varphi }+\mathbf{\text{N}}{\displaystyle \frac{\partial h}{\partial \mathit{\varphi }}}\right)/R=\\ =& \left(\mathbf{\text{W}}{\displaystyle \frac{\partial h}{\partial \mathit{\lambda }}}+\mathbf{\text{N}}\cos \mathit{\varphi }{\displaystyle \frac{\partial h}{\partial \mathit{\varphi }}}\right)/R\cos \mathit{\varphi }\end{array}\end{array}$$

Appendix B: Closed form integration along an arc

Using the notation for F, η, and ξ of Sect. 3.2,

$$\begin{array}{}\text{(B1)}& {\displaystyle }\mathbf{\text{W}}\cdot \mathbf{\text{F}}& {\displaystyle }=-\sin \mathit{\xi }\sin \mathit{\eta }/\cos \mathit{\varphi };\text{(B2)}& {\displaystyle }r& {\displaystyle }=\sin \mathit{\xi }\cos \mathit{\eta }.\end{array}$$

Assume that Φ was not constant throughout the arc but instead was a linearly function of η from Φ₁ to Φ₂. Then, the arc integral of Eq. (3.13) with Φ included is

$$\begin{array}{}\text{(B3)}& {\displaystyle }\begin{array}{rl}\int & R\mathrm{\Phi }\mathbf{\text{W}}\cdot \mathbf{\text{F}}\cos \mathit{\varphi }d\mathit{\eta }=\\ =& \underset{{\mathit{\eta }}_{\mathrm{1}}}{\overset{{\mathit{\eta }}_{\mathrm{2}}}{\int }}R\left({\displaystyle \frac{{\mathrm{\Phi }}_{\mathrm{1}}{\mathit{\eta }}_{\mathrm{2}}-{\mathrm{\Phi }}_{\mathrm{2}}{\mathit{\eta }}_{\mathrm{1}}}{{\mathit{\eta }}_{\mathrm{2}}-{\mathit{\eta }}_{\mathrm{1}}}}+\mathit{\eta }{\displaystyle \frac{{\mathrm{\Phi }}_{\mathrm{2}}-{\mathrm{\Phi }}_{\mathrm{1}}}{{\mathit{\eta }}_{\mathrm{2}}-{\mathit{\eta }}_{\mathrm{1}}}}\right)\mathbf{\text{W}}\cdot \mathbf{\text{F}}\cos \mathit{\varphi }d\mathit{\eta }=\\ =& -R\sin \mathit{\xi }\underset{{\mathit{\eta }}_{\mathrm{1}}}{\overset{{\mathit{\eta }}_{\mathrm{2}}}{\int }}\left({\displaystyle \frac{{\mathrm{\Phi }}_{\mathrm{1}}{\mathit{\eta }}_{\mathrm{2}}-{\mathrm{\Phi }}_{\mathrm{2}}{\mathit{\eta }}_{\mathrm{1}}}{{\mathit{\eta }}_{\mathrm{2}}-{\mathit{\eta }}_{\mathrm{1}}}}+\mathit{\eta }{\displaystyle \frac{{\mathrm{\Phi }}_{\mathrm{2}}-{\mathrm{\Phi }}_{\mathrm{1}}}{{\mathit{\eta }}_{\mathrm{2}}-{\mathit{\eta }}_{\mathrm{1}}}}\right)\sin \mathit{\eta }d\mathit{\eta }=\\ =& R({\mathrm{\Phi }}_{\mathrm{2}}{r}_{\mathrm{2}}-{\mathrm{\Phi }}_{\mathrm{1}}{r}_{\mathrm{1}})-R{\displaystyle \frac{{\mathrm{\Phi }}_{\mathrm{2}}-{\mathrm{\Phi }}_{\mathrm{1}}}{{\mathit{\eta }}_{\mathrm{2}}-{\mathit{\eta }}_{\mathrm{1}}}}\sin \mathit{\xi }(\sin {\mathit{\eta }}_{\mathrm{2}}-\sin {\mathit{\eta }}_{\mathrm{1}})=\\ =& R({\mathrm{\Phi }}_{\mathrm{2}}{r}_{\mathrm{2}}-{\mathrm{\Phi }}_{\mathrm{1}}{r}_{\mathrm{1}})-\\ & -R{\displaystyle \frac{{\mathrm{\Phi }}_{\mathrm{2}}-{\mathrm{\Phi }}_{\mathrm{1}}}{{\mathit{\eta }}_{\mathrm{2}}-{\mathit{\eta }}_{\mathrm{1}}}}\mathrm{2}\sin \mathit{\xi }\cos \left({\displaystyle \frac{{\mathit{\eta }}_{\mathrm{1}}+{\mathit{\eta }}_{\mathrm{2}}}{\mathrm{2}}}\right)\sin \left({\displaystyle \frac{{\mathit{\eta }}_{\mathrm{2}}-{\mathit{\eta }}_{\mathrm{1}}}{\mathrm{2}}}\right)=\\ =& R({\mathrm{\Phi }}_{\mathrm{2}}{r}_{\mathrm{2}}-{\mathrm{\Phi }}_{\mathrm{1}}{r}_{\mathrm{1}})-R{\displaystyle \frac{{\mathrm{\Phi }}_{\mathrm{2}}-{\mathrm{\Phi }}_{\mathrm{1}}}{{\mathit{\eta }}_{\mathrm{2}}-{\mathit{\eta }}_{\mathrm{1}}}}\mathrm{2}{r}_C\sin \left({\displaystyle \frac{{\mathit{\eta }}_{\mathrm{2}}-{\mathit{\eta }}_{\mathrm{1}}}{\mathrm{2}}}\right),\end{array}\end{array}$$

where r_C is the value of r at the center of the arc from η₁ to η₂. As the resolution of the grid becomes finer, r_C approaches $(r_{\mathrm{1}}+r_{\mathrm{2}})/\mathrm{2}$ and $\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

$$\begin{array}{ll}{\displaystyle }\mathrm{\Delta }\mathbf{\text{A}}& {\displaystyle }=R\mathrm{\Delta }t\sum \left[{\mathrm{\Phi }}_{n-\mathrm{1}}{\mathbf{\text{P}}}_{n-\mathrm{1}}-{\mathrm{\Phi }}_n{\mathbf{\text{P}}}_n+\mathrm{2}{\displaystyle \frac{{\mathrm{\Phi }}_n-{\mathrm{\Phi }}_{n-\mathrm{1}}}{{\mathit{\eta }}_n-{\mathit{\eta }}_{n-\mathrm{1}}}}\right.\\ \text{(B4)}& {\displaystyle }& {\displaystyle }\left.\sin \left({\displaystyle \frac{{\mathit{\eta }}_n-{\mathit{\eta }}_{n-\mathrm{1}}}{\mathrm{2}}}\right){\mathbf{\text{P}}}_{n-.5}\right]/K.\end{array}$$