2. The divergence and curl operators
a. The grid and coordinate system
b. Vector operators
3. Governing equations
4. Conservation principles
a. Conservation of the domain-mean mass and tracer
b. Conservation of the tracer variance
c. Conservation of total energy
1) Conservation of kinetic energy under advection
2) The energy conversion term
5. Properties of the finite-difference operators
6. Numerical tests
7. Discussion and conclusions
Comparison of Numerical Schemes
Mass equation
OS_Vor-Div scheme
NS_Vor_Div scheme
Vorticity equation
Divergence equation
Export References
Arakawa, A., and V. R. Lamb, 1977: Computational design of the basic dynamical process of the UCLA general circulation model. Methods in Computational Physics. J. Chang, Ed., Vol. 17. Academic Press, 173–265.
Arakawa, A., and V. R. Lamb, 1977: Computational design of the basic dynamical process of the UCLA general circulation model. Methods in Computational Physics. J. Chang, Ed., Vol. 17. Academic Press, 173–265.)| false
Arakawa, A., and V. R. Lamb, . 1981: A potential enstrophy and energy conserving scheme for the shallow water equations. Mon. Wea. Rev., 109 , 18–36.
Arakawa, A., and V. R. Lamb, . 1981: A potential enstrophy and energy conserving scheme for the shallow water equations. Mon. Wea. Rev., 109, 18–36.)| false
Baumgardner, J. R., and P. O. Frederickson, 1985: Icosahedral discretization of the two-sphere. SIAM J. Numer. Anal., 22 , 1107–1115.
Baumgardner, J. R., and P. O. Frederickson, 1985: Icosahedral discretization of the two-sphere. SIAM J. Numer. Anal., 22, 1107–1115.)| false
Giraldo, F. X., 2000: Lagrange–Galerkin methods on spherical geodesic grids: The shallow water equations. J. Comput. Phys., 160 , 336–368.
Giraldo, F. X., 2000: Lagrange–Galerkin methods on spherical geodesic grids: The shallow water equations. J. Comput. Phys., 160, 336–368.)| false
Heikes, R., 1993: The shallow water equations on a spherical geodesic grid. Department of Atmospheric Science Tech. Rep. 524, Colorado State University, 175 pp.
Heikes, R., 1993: The shallow water equations on a spherical geodesic grid. Department of Atmospheric Science Tech. Rep. 524, Colorado State University, 175 pp.)| false
Heikes, R., and D. A. Randall, 1995a: Numerical integration of the shallow-water equations on a twisted icosahedral grid. Part I: Basic design and results of tests. Mon. Wea. Rev., 123 , 1862–1880.
Heikes, R., and D. A. Randall, 1995a: Numerical integration of the shallow-water equations on a twisted icosahedral grid. Part I: Basic design and results of tests. Mon. Wea. Rev., 123, 1862–1880.)| false
Heikes, R., and D. A. Randall, . 1995b: Numerical integration of the shallow-water equations on a twisted icosahedral grid. Part II: A detailed description of the grid and an analysis of numerical accuracy. Mon. Wea. Rev., 123 , 1881–1887.
Heikes, R., and D. A. Randall, . 1995b: Numerical integration of the shallow-water equations on a twisted icosahedral grid. Part II: A detailed description of the grid and an analysis of numerical accuracy. Mon. Wea. Rev., 123, 1881–1887.)| false
Masuda, Y., and H. Ohnishi, 1987: An integration scheme of the primitive equations model with an icosahedral-hexagonal grid system and its application to the shallow water equations. Short- and Medium-Range Numerical Weather Prediction, T. Matsuno, Ed., Universal Academy Press, 317–326.
Masuda, Y., and H. Ohnishi, 1987: An integration scheme of the primitive equations model with an icosahedral-hexagonal grid system and its application to the shallow water equations. Short- and Medium-Range Numerical Weather Prediction, T. Matsuno, Ed., Universal Academy Press, 317–326.)| false
Pedlosky, J., 1987: Geophysical Fluid Dynamics. 2d ed. Springer-Verlag, 710 pp.
Pedlosky, J., 1987: Geophysical Fluid Dynamics. 2d ed. Springer-Verlag, 710 pp.)| false
Randall, D. A., 1994: Geostrophic adjustment and the finite-difference shallow water equations. Mon. Wea. Rev., 122 , 1371–1377.
Randall, D. A., 1994: Geostrophic adjustment and the finite-difference shallow water equations. Mon. Wea. Rev., 122, 1371–1377.)| false
Randall, D. A., R. Heikes, and T. D. Ringler, 2000: Global atmospheric modeling using a geodesic grid with an isentropic vertical coordinate. General Circulation Model Development, D. A. Randall, Ed., Academic Press, 509–538.
Randall, D. A., R. Heikes, and T. D. Ringler, 2000: Global atmospheric modeling using a geodesic grid with an isentropic vertical coordinate. General Circulation Model Development, D. A. Randall, Ed., Academic Press, 509–538.)| false
Ringler, T. D., and D. A. Randall, 2002: The ZM grid: An alternative to the Z grid. Mon. Wea. Rev., 130 , 1411–1422.
Ringler, T. D., and D. A. Randall, 2002: The ZM grid: An alternative to the Z grid. Mon. Wea. Rev., 130, 1411–1422.)| false
Ringler, T. D., R. P. Heikes, and D. A. Randall, 2000: Modeling the atmospheric general circulation using a spherical geodesic grid: A new class of dynamical cores. Mon. Wea. Rev., 128 , 2471–2490.
Ringler, T. D., R. P. Heikes, and D. A. Randall, 2000: Modeling the atmospheric general circulation using a spherical geodesic grid: A new class of dynamical cores. Mon. Wea. Rev., 128, 2471–2490.)| false
Sadourny, R., and P. Morel, 1969: A finite-difference approximation of the primitive equations for a hexagonal grid on a plane. Mon. Wea. Rev., 97 , 439–445.
Sadourny, R., and P. Morel, 1969: A finite-difference approximation of the primitive equations for a hexagonal grid on a plane. Mon. Wea. Rev., 97, 439–445.)| false
Sadourny, R., A. Arakawa, and Y. Mintz, 1968: Integration of the nondivergent barotropic vorticity equation with an icosahedal-hexagonal grid for the sphere. Mon. Wea. Rev., 96 , 351–356.
Sadourny, R., A. Arakawa, and Y. Mintz, 1968: Integration of the nondivergent barotropic vorticity equation with an icosahedal-hexagonal grid for the sphere. Mon. Wea. Rev., 96, 351–356.)| false
Salmon, R., 1998: Lectures on Geophysical Fluid Dynamics. Oxford University Press, 375 pp.
Salmon, R., 1998: Lectures on Geophysical Fluid Dynamics. Oxford University Press, 375 pp.)| false
Stuhne, G. R., and W. R. Peltier, 1996: Vortex erosion and amalgamation in a new model of large scale flow on a sphere. J. Comput. Phys., 128 , 58–81.
Stuhne, G. R., and W. R. Peltier, 1996: Vortex erosion and amalgamation in a new model of large scale flow on a sphere. J. Comput. Phys., 128, 58–81.)| false
Stuhne, G. R., and W. R. Peltier, . 1999: New icosahedral grid-point discretizations of the shallow water equations on the sphere. J. Comput. Phys., 144 , 23–58.
Stuhne, G. R., and W. R. Peltier, . 1999: New icosahedral grid-point discretizations of the shallow water equations on the sphere. J. Comput. Phys., 144, 23–58.)| false
Thuburn, J., 1997: A PV-based shallow-water model on a hexagonal-icosahedral grid. Mon. Wea. Rev., 125 , 2328–2347.
Thuburn, J., 1997: A PV-based shallow-water model on a hexagonal-icosahedral grid. Mon. Wea. Rev., 125, 2328–2347.)| false
Williamson, D. L., 1968: Integration of the barotropic vorticity equation on a spherical geodesic grid. Tellus, 20 , 642–653.
Williamson, D. L., 1968: Integration of the barotropic vorticity equation on a spherical geodesic grid. Tellus, 20, 642–653.)| false
Williamson, D. L., . 1969: Numerical integration of fluid flow over triangular grids. Mon. Wea. Rev., 97 , 885–895.
Williamson, D. L., . 1969: Numerical integration of fluid flow over triangular grids. Mon. Wea. Rev., 97, 885–895.)| false
Winninghoff, F. J., 1968: On the adjustment toward a geostrophic balance in a simple primitive equation model with application to the problems of initialization and objective analysis. Ph.D. thesis, University of California, Los Angeles, 161 pp.
Winninghoff, F. J., 1968: On the adjustment toward a geostrophic balance in a simple primitive equation model with application to the problems of initialization and objective analysis. Ph.D. thesis, University of California, Los Angeles, 161 pp.)| false
Export Figures
Grid points, shown as circles, can be connected to form a triangular grid shown by the dashed lines. Alternatively, cell walls can be positioned halfway between grid points to form a hexagonal grid shown by the solid lines. The hexagonal grid is the dual of the triangular grid
The grid is composed of hexagons (and possibly pentagons) with scalars defined at the cell centers and vectors defined at the cell corners. Each vector is described in the (e1, e2) coordinate system that is placed at each cell corner. The normal vectors to the cell walls are nc+ and nc−, where the c+ direction is in the counterclockwise direction from c and c− is in the clockwise direction. The area associated with each corner, Ac, is defined by the line segments connecting the grid-cell centers. R0c is the portion of Ac owned by cell 0. The area associated with the grid-cell center, Ai, is defined by the perimeter of the hexagon
The symbol hi for 0 ≤ i ≤ 2 represents an averaging of mass to the cell wall shared by grid cell i and i + 1. The area of the triangle is divided into three parts, Ric. The unit vectors normal to each cell-wall segment are given as ni and point from cell center i to cell center i + 1
Log base 10 of (top) the fractional change in total energy and (bottom) the mass-weighted potential enstrophy for each numerical simulation
(top) Energy power spectrum and (bottom) enstrophy power spectrum at day 40 for each numerical simulation, plotted as functions of total wavenumber. The open circles show the initial energy and enstrophy power spectra
(top) Total energy power spectrum and (bottom) enstrophy power spectrum averaged over days 950–1050. The open circles show the initial energy and enstrophy power spectra. Also shown in each panel, by a dashed line, is the theoretical prediction of the slope of the respective spectrum. The solid line denotes the equivalent wavenumber of the Rossby radius of deformation
The evolution of the relative vorticity field over the first 1000 days of integration. The relative vorticity scales range from approximately ±2.0 × 10−5, ±1.5 × 10−5, ±1.0 × 10−5, and ±3.0 × 10−6 s−1 at t = 0, t = 10, t = 100, and t = 1000 days, respectively.
View raw image
A Potential Enstrophy and Energy Conserving Numerical Scheme for Solution of the Shallow-Water Equations on a Geodesic Grid
Todd D. Ringler
Todd D. RinglerDepartment of Atmospheric Science, Colorado State University, Fort Collins, Colorado
Search for other papers by Todd D. Ringler in
David A. Randall
David A. RandallDepartment of Atmospheric Science, Colorado State University, Fort Collins, Colorado
Search for other papers by David A. Randall in
https://doi.org/10.1175/1520-0493(2002)130<1397:APEAEC>2.0.CO;2
Using the shallow water equations, a numerical framework on a spherical geodesic grid that conserves domain-integrated mass, potential vorticity, potential enstrophy, and total energy is developed. The numerical scheme is equally applicable to hexagonal grids on a plane and to spherical geodesic grids. This new numerical scheme is compared to its predecessor and it is shown that the new scheme does considerably better in conserving potential enstrophy and energy. Furthermore, in a simulation of geostrophic turbulence, the new numerical scheme produces energy and enstrophy spectra with slopes of approximately K−3 and K−1, respectively, where K is the total wavenumber. These slopes are in agreement with theoretical predictions. This work also exhibits a discrete momentum equation that is compatible with the Z-grid vorticity-divergence equation.
Corresponding author address: Todd Ringler, Department of Atmospheric Science, Colorado State University, Fort Collins, CO 80523-1371. Email: todd@atmos.colostate.edu
A spherical geodesic grid is a tessellation of the sphere that is generated by using the icosahedron as a starting point (e.g., Heikes and Randall 1995a). The potential applicability of spherical geodesic grids to the simulation of the atmospheric general circulation has been recognized since the 1960s (Sadourny et al. 1968; Williamson 1968, 1969; Sadourny and Morel 1969). Recently there has been renewed interest in this idea (Baumgardner and Frederickson 1985; Masuda and Ohnishi 1987; Heikes 1993; Heikes and Randall 1995a,b; Stuhne and Peltier 1996, 1999; Thuburn 1997; Giraldo 2000; Randall et al. 2000; Ringler et al. 2000).
Sadourny and Morel (1969) proposed that an optimal grid structure on the sphere should be as uniform as possible while preserving the isotropy of the spherical geometry. Spherical geodesic grids come close to meeting these ideal specifications. In terms of uniformity, a spherical geodesic grid can be constructed such that grid cell areas vary by less than 5% over the entire sphere (Heikes 1993). In terms of isotropy, relative to any given grid cell center all cell neighbors are nearly equidistant, and all the neighbors lie across cell walls (Randall et al. 2000). While Sadourny and Morel (1969) chose the spherical geodesic grid because of its uniformity, Williamson (1969) chose a spherical geodesic grid because the grid could be gradually distorted in space to produce increased resolution in critical regions such as the Gulf Stream.
In contrast to spherical geodesic grids, conventional latitude–longitude grids are neither uniform nor isotropic. The grid-pole singularities of latitude–longitude grids result in grid-cell areas that vary by O(1) over the globe. Furthermore, since latitude–longitude grids have cell neighbors that lie across both cell walls and cell corners, the grids are not isotropic.
While the merits of the spherical geodesic grid may have been appreciated decades ago, only recently have these merits been realized in simulations with full atmospheric general circulation models (AGCMs). Attempts to use spherical geodesic grids to solve the nondivergent shallow-water equations were successful (Sadourny et al. 1968; Williamson 1968). The difficulties involved in modeling the divergent shallow-water equations on the spherical geodesic grid were not overcome, however, until Masuda and Ohnishi (1987) proposed using the vorticity-divergence form of the equations. The vorticity-divergence formulation requires inverting elliptic equations at every time step to obtain the velocity field. The computational overhead required to do this with finite-difference models was generally considered prohibitive. Heikes and Randall (1995a,b, hereafter HR95) overcame this significant difficulty by implementing a multigrid method to invert the elliptic equations in a computationally efficient manner. Ringler et al. (2000) extended this framework from the shallow-water equations to the full 3D primitive equations and incorporated the physical parameterizations required for climate simulations.
Williamson (1969) chose to tile the sphere with triangles by connecting the grid points of the spherical geodesic grid (see Fig. 1). Alternatively, Sadourny and Morel (1969) used the inverse, or dual, of the triangular grid, which results in hexagonal grid cells also shown in Fig. 1. The lineage of work leading to the creation of a full AGCM based on a spherical geodesic grid in Ringler et al. (2000) emphasized the use of the hexagonal grid cells. An alternative line of research that emphasized the use of the triangular grid cells originated with both Williamson (1969) and Baumgardner and Frederickson (1985). The more recent works by Stuhne and Peltier (1996, 1999) and Giraldo (2000) use triangular finite elements to construct their numerical schemes. The present study makes extensive use of both the hexagonal and triangular grids.
As discussed by Arakawa and Lamb (1977), the problem of designing a numerical scheme for a general circulation model can be conceptually separated into the design of the linear properties of the scheme, and the design of its nonlinear properties. The numerical simulation of linear wave phenomena in geophysical fluid dynamics can be analyzed in terms of the geostrophic adjustment process (Winninghoff 1968; Arakawa and Lamb 1977; Randall 1994). The atmosphere is continually adjusting to forcings, such as diabatic heating, by the radiation of inertia–gravity waves, which leave behind an "adjusted" state that is close to geostrophic balance. In order to produce realistic results, a numerical scheme must faithfully reproduce this adjustment process. Following the lead of Masuda and Ohnishi (1987), HR95 used the "Z grid," on which vorticity, divergence, and mass are prognosed, in part because, as shown by Randall (1994), the Z grid simulates geostrophic adjustment very well; that is, schemes based on the Z grid have good linear properties. HR95 used a finite-volume method which guarantees conservation of mass, including tracer mass, as well as potential vorticity. HR95 did not discuss the conservation of nonlinear functions of the prognostic fields, such as the kinetic energy, total energy, and potential enstrophy.
The present study is aimed at modifying the scheme of HR95 so as to guarantee conservation of kinetic energy, total energy, and potential enstrophy under frictionless processes. Such conservation properties are particularly important in long-term simulations such as those that are needed for the study of climate. On long timescales, sources and sinks of energy and other quantities are fundamental to the circulation. Small but systematic spurious sources and sinks of fundamental quantities, such as energy, can lead to unrealistic circulation regimes.
The continuous two-dimensional nonlinear shallow-water equations are highly constrained in terms of their energy and enstrophy cascades. The net effect of advection is to transfer energy from shorter to longer spatial scales, while transferring enstrophy from longer to shorter spatial scales. Arakawa and Lamb (1977, 1981) point out that we have little hope of accurately simulating such energy cascades and the associated energy spectra if the numerical schemes do not conserve energy and potential enstrophy. They show that conserving basic quantities, such as mass, potential vorticity, potential enstrophy, and energy, can dramatically increase the overall accuracy of long numerical simulations, even with no change in the local "order of accuracy" of the schemes. While schemes of higher-order accuracy are generally to be preferred, all other things being equal, a formal increase in order of accuracy should not be chosen at the expense of the conservation principles. Arakawa and Lamb (1977) provide an example of a lower-order conservative numerical scheme that produces a more realistic simulation than a higher-order nonconserving scheme.
The scheme presented in this paper is equally applicable to planar hexagonal grids (see Fig. 1) and a spherical geodesic grid. The spherical geodesic grid is quite similar to the hexagonal grid; all the grid cells are hexagons, with the exception of 12 grid cells that are pentagons (Sadourny and Morel 1969). In section 2 we construct some of the building blocks of the numerical scheme, namely the divergence, and curl operators, from their respective definitions. We then use these basic operators, in section 3, to obtain a discrete form of the shallow-water equations. In section 4 we show that this discrete system conserves a number of important quantities: mass, potential vorticity, potential enstropy, and total energy. In section 5 we demonstrate some properties of our finite-difference operators. This numerical scheme has been implemented in a shallow water model on a doubly periodic plane. Results are shown in section 6, and these results are compared and contrasted to those produced with the numerical scheme currently used in the Colorado State University (CSU) AGCM (Ringler et al. 2000). Section 7 gives a discussion that places our numerical scheme into context.
The numerical scheme developed below is defined on a grid of polygons that are either hexagonal or pentagonal in shape. When discretizing the surface of a plane, all the polygons are hexagons. When discretizing the sphere, all of the polygons are hexagons with the exception of 12 polygons that are pentagons.
All scalars, such as the fluid depth, are defined at the centers of the grid cells and are referenced with the subscript i, while all vectors, such as the velocity vector, will be defined at the cell corners and referenced with the subscript c. When we wish to analyze a specific corner, we will set c = γ as shown in Fig. 2 and use the data from the surrounding numbered cell centers. In addition to cell centers and cell corners, we will need to reference segments of the cell walls defined by the symbol d. We will use the subscripts c+ and c− to denote the cell wall segments on d in the counterclockwise and clockwise directions from corner c, respectively. The unit vectors nc+ and nc− denote the directions perpendicular to dc+ and dc−, respectively, and point outward relative to cell center 0.
At each corner we have placed an arbitrary orthogonal coordinate system defined by the unit vectors (e1, e2). The area fluxes across the dc+ and dc− wall segments are Vc · nc+dc+, and Vc · nc−dc−, respectively. We use the symbols Ai and Ac to represent the total areas associated with the grid cell centers and corners, respectively. The area of the shaded region in Fig. 2, denoted by Ric, is the portion of Ac that is associated with cell center i.
With the coordinate system defined in Fig. 2, we can derive numerical approximations to the divergence and curl operators from their respective analytic definitions.1 The analytic form of the divergence operator is
The divergence operator yields a scalar, which will be defined at cell centers, by summing over the dot product of vector data defined at cell corners. Referring to Fig. 2, we can approximate the divergence operator at the cell centers as
where the summation is over all cell corners associated with grid cell i. The variable nc is equal to 6 for hexagons and 5 for pentagons. If we define
then we can rewrite (2) in the more compact form
The definition of the curl operator is
where we have explicitly selected the "vertical" component of the curl by dotting with k. Since the curl is a scalar quantity, it will be defined at grid-cell centers. We can approximate (5) as
where, as with the divergence operator, the summation is over all cell corners associated with the cell center i.
The shallow-water equations can be written as
where h is the fluid depth, hs is the height of the surface, V is the vector velocity, K is the kinetic energy, and η is the absolute vorticity. Equation (9) describes the evolution of an arbitrary tracer τ, where τ is a mixing ratio. Alternatively, we can take the curl and divergence of (8) to generate equations for the vorticity and divergence:
Here η ≡ f + k · ∇ × V is the absolute vorticity and δ ≡ ∇·V is the divergence. In the continuous equations, the vector momentum formulation and the vorticity-divergence formulation are equivalent. Given η and δ we can determine V in two steps. First, we solve the elliptic equations
for the streamfunction, ψ, and velocity potential, χ. We then compute the velocity using the relation
Vkψχ. (14)
Given the discrete analogs to the divergence and curl operators derived in section 2, we can write the discrete forms of (7), (8), and (9) as
where (15) and (17) describe the evolution of mass and mass-weighted tracer within a grid cell. Equation (16) describes the evolution of the velocity at a cell corner. For clarity we reiterate our notation: terms with a subscript i are defined at cell centers, terms with a subscript c are defined at cell corners, and terms with a subscript c+ or c− are defined at the cell walls. The terms hc+, τc+, hc−, and τc− should be interpreted as averages of mass and tracer from the cell centers to the cell wall segment c+ and c−. The term ηc is the average of absolute vorticity from the cell centers to the cell corners. All of these symbols are undefined at this point. The form of the kinetic energy, Ki, is undefined, as is the discrete form of the gradient operator.
The discrete forms of the vorticity and divergence equations can be written down by analogy with the analytic forms shown in (10) and (11):
Since we have not yet specified the form of the discrete gradient operator, the discrete Laplacian operator is unspecified. In section 5 we will show that (18) and (19) can, alternatively, be derived from the discrete form of the momentum equation. As a result, the discrete momentum formulation will be entirely consistent with the discrete vorticity-divergence formulation.
In the shallow-water equations the potential vorticity, q, is equal to the absolute vorticity divided by the layer thickness, h, so we can rewrite (18) as
which is identical in form to the tracer equation.
While the continuous equations allow an infinite number of quantities to be conserved, this is not possible within the discrete system. In this section we will show that with the appropriate choices of hc+, τc+, hc−, τc−, ηc, Ki, and the gradient operator, the discrete shallow-water equations shown above can be implemented such that the following quantities are conserved in the domain mean: mass, mass-weighted tracer, and mass-weighted potential vorticity, mass-weighted tracer variance, mass-weighted potential enstrophy, and total energy.
The domain-integrated mass, mass-weighted tracer, and mass-weighted potential vorticity are conserved simply by virtue of their flux forms, written in (15), (17), and (20), respectively. The only stipulation is that hc+, τc+, hc−, τc−, qc+, and qc− have the same values when referenced by either of the cell centers that share a given cell wall. This is equivalent to the requirement that the flux out of one cell is identical to the flux into its neighbor. This will become more explicit later when we more precisely define these quantities.
In order to show how the mass-weighted tracer variance can be conserved, we must start with the advective form of the tracer equation. The advective form is obtained by taking (17) − τi × (15) to yield
Following Arakawa and Lamb (1977) we introduce a function G = G(τ) that depends only on τ. We can construct an equation for G by taking (15) × G(τ) and adding it to (21) × dG/dτ. This results in
If we choose c = γ, as shown in Fig. 2, and recognize that Fc+ represents the flux between cell 0 and cell 2, while Fc− represents the flux between cell 0 and cell 1, we can see that G will be conserved if
Solving (23) for τc+ and (24) for τc−, we obtain
As an example, if we choose G(τ) = τ2 we can solve (25) and (26) for τc+ and τc− to obtain
Note that no additional constraints on hc+ or hc− are required.
Since the potential vorticity equation, (20), is identical in form to the tracer equation, (17), we can use the results of (27) and (28) to guarantee the conservation of potential enstrophy in the discrete system. We can write the vorticity equation by using the lhs of (18) and the rhs of (20) to yield
Instead of averaging vorticity to the cell walls, we average potential vorticity to the cell walls using (27) and (28) and then multiply by the averaged mass.
The total energy in the discrete shallow-water equations is conserved. We show this in two steps: first, that kinetic energy is conserved under the process of advection, and second, that the energy conversion term neither creates nor destroys total energy.
In this section we will find the form of the gradient operator that guarantees that the process of advection neither creates nor destroys kinetic energy. First, we must form the kinetic energy equation by taking the scalar product of Vc and the momentum equation given in (16):
This equation is valid for every cell corner. The first term on the rhs of (30) is identically zero at every cell corner because Vc is perpendicular to k × Vc, even in the discrete case. Recall that while the momentum equation is defined at cell corners, the kinetic energy is a scalar and will, therefore, be defined at cell centers. To move (30) from the corners to the centers, we weight each momentum point by the area, Ric, as defined in Fig. 2 and sum over all corners associated with grid cell i to yield
where we have multiplied by the cell mass, hi, after summing over the corners. Equation (31) holds at every cell center.
We now multiply the continuity equation, (15), by Ki to obtain
where we have moved the cell area, Ai, to the lhs of the equation. Equations (31) and (32) are both defined at cell centers, so we can add them together to yield
where we have used (3) to express (33) completely in terms of the velocity. The last term on the rhs of (33) is an energy conversion term. At this point we will drop the last term on the rhs of (33). We return to this term in the next section.
Summing both sides of (33) over the entire domain yields
where it is understood that the outside summation is over all the grid cell centers and the inside summation is over the corners associated with a given cell center i.
We will now show that the rhs of (34) sums to zero at the corner labeled c = γ in Fig. 2. Symmetry then implies that the sum is zero at every corner. Note that in Fig. 2 there is no way to tell whether the grid cells that share the corner c = γ are hexagons or pentagons. Both types of polygons are accounted for in this numerical scheme without exception.
At this point it is convenient to adopt a new naming convention as shown in Fig. 3. Since we are choosing a specific corner, we will drop the subscript c. The area of the triangle is partitioned into three subareas, Ri, for 0 ≤ i ≤ 2. The vectors normal to each cell-wall segment, ni, point away from cell center i toward cell center i + 1, where i + 1 is cyclic. The averaging of mass to cell walls, hi, denotes a yet-unspecified averaging from cell centers to the cell-wall segment shared by cell center i and i + 1.
If we sum (34) over all grid cells but keep only the terms in the summation that have c = γ and expand, we are left with
where we require that the expression will sum to zero. Now let
for all i between 0 and 2, and let S = R0 + R1 + R2. We can then rewrite (35) in the compact form of
where there is an implied sum over the index i and i − 1 is cyclic. We can factor out the velocity and rewrite (37) as
We want (38) to be true for an arbitrary vector velocity field, so (38) reduces to
If we take the dot product of (39) with the unit vectors that define the local coordinate system and solve for the components of ∇K, we find
Equations (40) and (41) serve as the definition of the gradient operator. Note that the form of K is still unspecified; the rhs of (34) sums to zero for arbitrary scalar fields. We will return to the gradient operator in the next section.
We can now write (34) as
If we interpret (42) as two terms in a chain rule expansion, we have
This will be true if we define the kinetic energy as
Equations (36), (39), and (44) are the constraints that allow kinetic energy to be conserved under advection.
In addition to showing that advection does not create or destroy kinetic energy, we must show that the terms that represent the conversion of energy between its kinetic and potential forms do not create or destroy energy. In order to do this we must derive the potential energy equation by multiplying gAi(hs + hi) through (15) to obtain
In an approach similar to that used above, we will show that the rhs of (45) and the last term in (33) cancel at every cell corner. Combining the rhs of (45) and the last term on the rhs of (33) and summing over the domain, we want to ensure that
If we compare (46) to the rhs of (34), we see that the two expressions are identical, except that K in (34) is replaced by g(hs + h) in (46). Since we have already shown that the rhs of (34) sums to zero for an arbitrary scalar, (46) is satisfied provided that we use an expression analogous to (39) to compute the gradient of g(hs + h).
The form of total energy that is conserved in this discrete system is
which is consistent with the continuous shallow-water equations.
We have shown that we can conserve total energy by choosing (39) as the discrete gradient operator. In this section we will show that (39) has an alternative geometric interpretation as the slope of the plane fit through the surrounding three data points. For clarity, we will assume a regular hexagonal grid in this section, but note that we have verified our findings on both distorted hexagonal grids and spherical geodesic grids.
Referring to Fig. 2 and assuming a regular hexagonal grid, if we take the dot product of (39) with the unit vector nc+, we find that
where Δl is the distance between grid cells 2 and 0. Equation (48) can be interpreted as the slope of the plane fit through (K0, K1, K2) in the nc+ direction. Taking the dot product of (39) with nc− yields
Equations (48) and (49) are sufficient to prove that (39) is the slope of the plane fit through (K0, K1, K2).
Once we interpret (39) as the slope of the plane fit through the surrounding data, the vector identity,
follows immediately in the discrete system. The physical interpretation of the curl operator shown in (5) and (6) is a line integral of the component of a vector field along the path of integration. In Fig. 2, if we start at γ and follow ∇K along the perimeter of grid cell 0, we return to where we started. Hence, (50) holds in the discrete system. While this is a rough argument, we have proven (50) rigorously.
With the divergence and gradient operators specified, we can determine the form of the discrete Laplacian operator as ∇2K = ∇·∇K, where K is a scalar field. If we define the vector V in (2) as ∇K and use (39) as the form of the gradient operator, we find that
where Δt is the length of the line segment connecting grid cell corners. Equation (51) can be simplified to
which is exactly the same form as derived by HR95.
Along with the vector identity shown in (50) and the definition shown in (52), we require two more equations to hold in the discrete system:
We simply note here that (53) and (54) hold in the discrete case solely because we assume that the vector V is perpendicular to the vector k.
Equations (50), (52), (53), and (54) suffice to derive the discrete vorticity and divergence equations [(18) and (19)] from the discrete momentum equation (16). Thus we conclude that the discrete momentum formulation is consistent with the discrete vorticity-divergence formulation.
We will show results from three numerical models. Two of the models are based on the framework developed above. The first uses the discrete vorticity-divergence formulation [Eqs. (15), (18), and (19)] and we will refer to this model as NS_Vor-Div. The second uses the discrete vector momentum formulation [Eqs. (15) and (16)] and will be referred to as NS_Momentum. The third model is based on the scheme developed by Masuda and Ohnishi (1987) and HR95. This scheme uses the vorticity-divergence formulation of the shallow-water equations and will be referred to as OS_Vor-Div. A fairly complete comparison of OS_Vor-Div and NS_Vor-Div is contained in the appendix.
The models are situated on a doubly periodic f plane with f0 = 1.4 × 10−4 s−1. The plane is discretized using a regular hexagonal grid of 128 × 128 (see Fig. 7). The distance between cell centers, Δn, is 100 km. The time-stepping scheme for all three models is third-order Adams–Bashforth with a time step of 100 s. All three models are started from the same initial condition, of the form
The first values on the rhs of (55), (56), and (57) are the respective domain mean values and the second values are the maximum perturbations from the means. The perturbations were generated by a random number generator. The initial conditions are highly unbalanced. Given the absolute vorticity and divergence fields, the discrete analogs of (12) and (13) are used to solve for the streamfunction and velocity potential, respectively. For the NS_Momentum model, we use the discrete gradient operator to differentiate the streamfunction and velocity potential to obtain the initial velocity field. The time rates of change of all prognostic variables are zero at t = 0. The surface topography is random and has the form hs = 0 ± 20 m. With gravity set to 9.81 m s−2, we obtain a Rossby radius of deformation of approximately 700 km. With Δn = 100 km, we marginally resolve eddies on this scale. All three models are inviscid in the sense that no explicit dissipation is introduced.
Figure 4a measures the error in domain-integrated total energy over a 40-day integration. The y axis is log base 10 of the absolute value of the fractional change in domain-integrated total energy. All of the simulations show substantial variability during the first few days of integration. During this period the state of the system is rapidly changing while it moves toward geostrophic balance. The time truncation error is relatively large during this period and we have verified that it accounts for the variations of the total energy over the first several days of integration. After 40 days, the NS_Vor-Div and NS_Momentum simulations have a total energy that is within 0.5% of the initial value. The total energy in the OS_Vor-Div simulation continues to drift and after 40 days the total energy has doubled.
Figure 4b shows the fractional change in domain-integrated mass-weighted potential enstrophy. Both the NS_Vor-Div and NS_Momentum simulations conserve mass-weighted potential enstrophy to within 0.05% after 40 days of integration. The NS_Vor-Div simulation displays a slow modulation of this quantity that is due to truncation error in the elliptic solver [(12) and (13)]. If we make the convergence threshold in the elliptic solver sufficiently stringent, we can eliminate this modulation. The OS_Vor-Div simulation shows a steady drift in potential enstrophy and after 40 days of integration the potential enstrophy has also doubled.
Figure 5a shows the wavenumber spectrum of total energy at the end of the 40-day integration for each simulation. Figure 5a also shows the energy spectrum at t = 0, which is the same of all simulations. Consistent with the white noise initial condition, the initial energy spectrum is also white. Figure 5b shows the enstrophy spectrum of each simulation after 40 days and also shows the initial enstrophy spectrum. The purpose of Fig. 5 is to show whether there is anomalous buildup of energy or enstrophy at any given wavenumber. In terms of energy, both the NS_Vor-Div simulation and NS_Momentum simulation show an appropriate upscale transport of energy with no apparent buildup of energy at the grid scale. The OS_Vor-Div simulation also shows the upscale cascade, but the lack of conservation of total energy is readily apparent. In terms of enstrophy, both the NS_Vor-Div simulation and NS_Momentum simulation show no buildup of enstrophy at the smallest scales. The OS_Vor-Div simulation shows a buildup of enstrophy at both the smallest and largest scales.
Now we wish to determine whether the new numerical scheme produces the appropriate energy and enstrophy cascades. In order to test this, we must include a sink of enstrophy at the smallest resolved scales in order to generate a downscale cascade of enstrophy. Using the NS_Momentum model, we do this by including a ∇6 diffusion in the momentum equation.2 The coefficient on the diffusion is 1 × 1024 m6 s−1; scale analysis suggests that this value is sufficient to dissipate a vorticity or divergence perturbation of size 1 × 10−4 s−1, with a spatial scale of 2Δn, in several hours. The model is initialized with the same initial condition shown above, but there is no topography. Figure 6 shows the energy and enstrophy spectra averaged over days 950–1050. Dimensional analysis suggests that, within the inertial range, the energy spectrum should decay as K−3, where K is the total wavenumber, and the enstrophy spectrum should decay as K−1 (Pedlosky 1987; Salmon 1998). Figure 6 indicates that the numerical model is accurately capturing both the upscale cascade of energy and the downscale cascade of enstrophy. After 1000 days of integration, the total energy has changed by only 0.5% (not shown). So while the ∇6 diffusion operator continues to destroy enstrophy, it destroys very little energy. Figure 7 shows the relative vorticity field at four stages of the integration: t = 0, t = 10, t = 100, and t = 1000 days. This figure qualitatively confirms the spectral analysis shown in Fig. 6. The flow continually evolves into larger structures. Also consistent with Fig. 6, we do not see any spurious buildup of energy or enstrophy at the grid scale.
We have constructed a shallow-water model, for use on a spherical geodesic grid, that conserves mass, potential vorticity, potential enstrophy, and total energy. We have demonstrated these conservation properties in a numerical simulation of geostrophic turbulence on an f plane. In addition to conserving energy and potential enstrophy to within the time truncation error, our numerical scheme was able to accurately capture the spectral distributions of energy and enstrophy in a simulation of geostrophic turbulence.
In the introduction we stated that conserving these basic quantities is only one necessary component of a robust numerical scheme; faithful simulation of the geostrophic adjustment process and numerical convergence are the other two necessary components. Looking at these other two necessary components will provide some context for the work we have completed here.
In a study of the geostrophic adjustment process, Randall (1994) discretizes the vorticity-divergence form of the shallow-water equations on the Z grid. The Z grid entails the use of the vorticity-divergence formulation with all scalar quantities defined at the grid-cell centers. He compares the discrete dispersion relation obtained using the Z grid to the dispersion relations obtained using the A, B, C, D, and E grids (Arakawa and Lamb 1977). The A–E grids entail the use of the momentum formulation. Randall (1994) clearly shows that the Z grid more accurately simulates the geostrophic adjustment process than any of the A–E-grid systems. On the Z grid, the discrete equations involve no spatial averaging and only a discrete Laplacian operator is required [see Randall (1994), Eqs. (4)–(6)]. As shown in (52), the Laplacian operator derived here is identical to that used by HR95. Furthermore, this Laplacian operator is consistent with the formulation used in Randall (1994). We conclude that the framework developed here not only conserves the quantities listed above, but is also consistent with the Z-grid discretization and, thus, does a better job of simulating the geostrophic adjustment process than the A–E-grid systems. This topic is more fully explored in Ringler and Randall (2002).
The order of accuracy of this numerical scheme is determined by three aspects: the accuracy of the interpolations [Eqs. (27), (28), and (36)], the accuracy of the vector operators [Eqs. (2), (6), and (39)], and the properties of the grid. On a regular hexagonal grid, all of these approximations are formally second-order accurate.
We interpreted the approximation to the gradient operator in two ways: first, the gradient ensures conservation of kinetic energy under the process of advection, and second, the gradient operator measures the slope of the plane fit through the surrounding three data points.
This work has shed light on the relationship between the Z grid and its discrete analog in the momentum equations. We have exhibited a momentum equation that is compatible with the Z grid. In the continuous system, we can move between the momentum formulation and vorticity-divergence formulation in a seamless manner using basic vector operators. Both forms contain the same information and yield the same result. In this work we were able to manipulate the discrete form of the vector momentum equation using discrete analogs of the gradient, divergence, and curl operators to derive a discrete form of vorticity-divergence equations. Since the derived vorticity-divergence equations are the Z-grid vorticity-divergence equations, we are guaranteed that the discrete momentum equation is consistent with the Z grid. The momentum analog to the Z grid (call it the ZM grid) solves the full momentum equation at every grid-cell corner.
On the hexagonal grid there are twice as many grid-cell corners as there are grid-cell centers. An easy way to see this is as follows. Each cell center is associated with an area, Ai, and each cell corner is associated with an area, Ac. On a regular hexagonal grid Ai = 2Ac, so in order to cover the same area there must be twice as many cell corners as there are cell centers. Since the ZM grid solves the momentum equation at every cell corner and the continuity equation at every cell center, there are twice as many momentum points as there are mass points. The redundancy of the momentum equation allows for the existence of computational modes. We define computational modes as solutions to the discrete equations that have no analogs in the continuous system. We are addressing the redundancy of the ZM grid in separate work by analyzing the linear geostrophic adjustment process on a hexagonal grid (Ringler and Randall 2002).
We have implemented the numerical scheme outlined in this paper on a spherical geodesic grid and we are currently conducting the standard suite of shallow-water test cases as well as several other tests.
We would like to thank Ross Heikes and Alistair Adcroft for useful discussions and for comments on an earlier draft of this paper. This work was supported by the U.S. Department of Energy's Climate Change Prediction Program under Grant DE-FG03-98ER62611 to Colorado State University.
In this appendix we compare the OS_Vor-Div scheme developed by Masuda and Ohnishi (1987) and HR95 to our new scheme, NS_Vor-Div. For this comparison we will assume a regular hexagonal grid where the distance between grid-cell centers is Δn and the distance between grid-cell corners is Δt.
Both schemes express the mass equation as fluxes across the cell walls and we will denote these fluxes as Fj, where the flux points from cell 0 to cell j. For both schemes, the mass equation for the fluid thickness at cell center 0 can be written as
The schemes will differ in the formulation of the flux, Fj.
The flux across cell walls can be express as
The first term on the rhs of (A2) represents the flux of mass across the cell wall shared by grid cells 0 and j due to divergent motion. The second term on the rhs of (A2) represents the flux across the cell wall due to rotational motion. Substituting (A2) into (A1) yields
The flux across each cell wall, Fj, is broken into two parts: one part associated with the corner j − 1, and one part associated with the corner j. The mass equation can be written just like (A1), but the fluxes are expressed as
Substituting (A4) into (A1) and doing some algebra yields
Comparison with (A3) shows that the rotational motion is modeled exactly as before, but that a slight modification to the flux due to divergent motion is introduced.
Similar to the mass equation, the vorticity equation for both schemes can be written as
where, similar to Fj in the mass equation [(A1)], Hj is the flux of absolute vorticity across the cell walls. In fact, the vorticity equation is identical to the mass equation, expect h is replaced by η. Furthermore, in both schemes the averaging of absolute vorticity is done in exactly the same manner as the averaging of mass. Instead of duplicating the same manipulations we just completed, we can simply state that the OS_Vor-Div scheme integrates the vorticity equation as shown in (A3) and the NS_Vor-Div scheme integrates the vorticity equation as shown in (A5).
For clarity, we reiterate the continuous form of the divergence equations:
which we can write symbolically as
In section 5 we show that both the OS_Vor-Div and NS-Vor-Div schemes use the same discrete Laplacian operator; therefore, regarding term 2, the only differences between the two schemes could be the form of kinetic energy, K. In comparing the two schemes we need to identify any differences in the treatment of term 1 or the formulation of K.
This scheme approximates term 1 as
Masuda and Ohnishi (1987) express the kinetic energy in continuous form as
Equation (A10) is derived in Heikes (1993). If we use the formulation developed by Masuda and Ohnishi (1987) we can write the discrete form of (A10) at cell center 0 as
In the process of manipulating the discrete analog of (A10) into (A11), we combined the first two terms on the rhs of (A10) into the first term on the rhs of (A11). The third and fourth terms on the rhs of (A10) were combined into the second term on the rhs of (A11). And finally, the last term on the rhs of (A10) was expanded to form the last two terms on the rhs of (A11). We can now write (A11) in terms of the velocity components as
where R and D are the rotational and divergent velocities, respectively. The subscripts N and T are the normal and tangent directions relative to each cell wall, respectively. For example, RT is the rotational part of the velocity vector in the tangent direction. There are several points to note. First, (A12) is not positive definite; the last two terms can be negative. Second, (A12) is biased in the sense that certain components of the velocity vector, such as RT and DN, more strongly control the value of K0 than other components of the velocity, such as RN and DT.
In comparing this to OS_Vor-Div we see that the part of term 1 due to divergent motion is modeled exactly the same in the two schemes, but the part due to rotational motion is different.
The form of kinetic energy used in NS_Vor-Div is shown in (44). When we use the vorticity-divergence formulation, we can rewrite (44) using (14) as
In practice, we use the gradient operator to compute (∇ψ)j and (∇χ)j at each corner, then we construct the velocity vector as each corner. We compute the kinetic energy as the sum of (1/2)(Vj · Vj) over the cell corners. The important point to note is that our formulation uses the full, unbiased, velocity vector to compute the kinetic energy.
Citation: Monthly Weather Review 130, 5; 10.1175/1520-0493(2002)130<1397:APEAEC>2.0.CO;2
We could also derive an approximation to the gradient operator from its analytical definition. Instead, we will use the constraint of conservation of total energy (section 4) to define the form of the gradient operator.
If we use the NS_Vor-Div model with a ∇6 diffusion operator of vorticity and divergence, our findings are the same.
https://journals.ametsoc.org/view/journals/mwre/130/5/1520-0493_2002_130_1397_apeaec_2.0.co_2.xml