Journal cover Journal topic
Geoscientific Model Development An interactive open-access journal of the European Geosciences Union
Journal topic
GMD | Articles | Volume 13, issue 2
Geosci. Model Dev., 13, 581–595, 2020
https://doi.org/10.5194/gmd-13-581-2020
© Author(s) 2020. This work is distributed under
the Creative Commons Attribution 4.0 License.
Geosci. Model Dev., 13, 581–595, 2020
https://doi.org/10.5194/gmd-13-581-2020
© Author(s) 2020. This work is distributed under
the Creative Commons Attribution 4.0 License.

Development and technical paper 14 Feb 2020

Development and technical paper | 14 Feb 2020

Extending square conservation to arbitrarily structured C-grids with shallow water equations

Extending square conservation to arbitrarily structured C-grids with shallow water equations
Lilong Zhou1,3, Jinming Feng2, Lijuan Hua4, and Linhao Zhong2 Lilong Zhou et al.
  • 1College of Earth and Planetary Sciences, University of Chinese Academy of Sciences, Beijing, China
  • 2Key Laboratory of Regional Climate-Environment for Temperate East Asia, Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing, China
  • 3Numerical Weather Prediction Center of China Meteorological Administration, Beijing, China
  • 4Key Laboratory of Computational Geodynamics, University of Chinese Academy of Sciences, Beijing, China

Correspondence: Jinming Feng (fengjm@tea.ac.cn)

Abstract
Back to toptop

The square conservation law is implemented in atmospheric dynamic cores on latitude–longitude grids, but it is rarely implemented on quasi-uniform grids, given the difficulty involved in constructing anti-symmetrical spatial discrete operators on these grids. Increasingly more models are being developed on quasi-uniform grids, such as arbitrarily structured C-grids. Thuburn–Ringler–Skamarock–Klemp (TRiSK) is a shallow water dynamic core on an arbitrarily structured C-grid. The spatial discrete operator of TRiSK is able to naturally maintain the conservation properties of total mass and total absolute vorticity and conserving total energy with time truncation error; the first two integral invariants are exactly conserved during integration, but the total energy dissipates when using the dissipative temporal integration schemes, i.e., Runge–Kutta (RK). The method of strictly conserving the total energy simultaneously, which means conserving energy in the round-off error over the entire temporal integration period, uses both an anti-symmetrical spatial discrete operator and a square conservative temporal integration scheme. In this study, we demonstrate that square conservation is equivalent to energy conservation in both a continuous shallow water system and a discrete shallow water system of TRiSK. After that, we attempt to extend the square conservation law to the TRiSK framework. To overcome the challenge of constructing an anti-symmetrical spatial discrete operator, we unify the unit of evolution variables of shallow water equations using the Institute of Atmospheric Physics (IAP) transformation, and the temporal derivatives of new evolution variables can be expressed by a combination of temporal derivatives of original evolution variables, which means the square conservative spatial discrete operator can be obtain by using original spatial discrete operators in TRiSK. Using the square conservative Runge–Kutta scheme, the total energy is completely conserved, and there is no influence on the properties of conserving total mass and total absolute vorticity. In the standard shallow water numerical test, the square conservative scheme not only helps maintain total conservation of the three integral invariants but also creates less simulation error norms.

1 Introduction
Back to toptop

In a statistical sense, the integral constraints make the physics of the discrete model more analogous to the physics of the continuous atmosphere, and also make the errors less systematic (Arakawa and Lam, 1977). Shallow water equation sets, without any outer sources and frictions, have five basic physical conservation properties, including total mass, total energy, total absolute vorticity (total potential vorticity), total potential enstrophy and total angular momentum. These conservation properties are important in an atmosphere model, especially with regard to long-term simulation; however, in a discrete system, some conservation properties cannot be maintained (Wang, 2008). If the square of a quantity is conserved with time when summed up over all the grid points in a domain, the quantity itself will be bounded, at every individual grid point, throughout the entire period of integration; this might be helpful for preventing nonlinear computational instability (Arakawa, 1966), and energy is one kind of quadratic quantity. Toy and Nair (2017) developed an energy and potential enstrophy conserving scheme for shallow water equations on generalized curvilinear coordinates; they mentioned conserving analogues of total energy, and total potential enstrophy in numerical models is known to prevent a spurious cascade of energy toward small scales. For a short-term simulation, the influences of slight energy dissipation are not obvious, but this dissipation accumulates in every time step, and finally, in a long-term simulation, leads to a quiet different result; i.e., an ellipse orbit tends to a single point after 108 steps (Wang et al., 1996).

A numerical scheme, with an energy conservation or energy dissipative property, is prerequisite to prevent nonlinear computational instability; however, an energy dissipative scheme will limit short waves, which are meaningful for the atmosphere (Shen et al., 2015; Zeng and Ji, 1981). On a latitude–longitude grid, energy can be entirely conserved by constructing a square conservative finite-difference scheme (Ji and Wang, 1991) or a multi-conservation finite-difference scheme (Wang and Ji, 2003), the former being better developed. Wang and Ji (1994a) discussed the square conservative scheme (SCS), the complete square conservative scheme (CSCS), the instantaneous square conservative scheme (ISCS) and the explicit complete square conservative scheme with adjustable time intervals (ECSCSA). The ISCS maintains square conservation only in the spatial discrete scheme and not in the temporal integration scheme, which implies that the spatial discrete operator of the model is a square conservative (i.e., an anti-symmetrical operator). However, the temporal integration scheme does not possess the square conservation property; therefore, the model is energy dissipative during integration. The CSCS maintains square conservation in both the spatial and temporal schemes. The model, which adopts CSCS, is able to maintain complete energy conservation during integration. The first step of applying the square conservation theory is to construct an anti-symmetrical spatial discrete operator, and the second step is to integrate the model with a square conservative temporal integration scheme, e.g., a modified predictor–corrector, modified leap-frog (Wang and Ji, 2006), harmonious dissipative operators (Wang and Ji, 1994b), etc.

To improve integration precision on the temporal direction of the square conservative scheme, a new class of new Runge–Kutta (NRK) schemes were developed by Wang et al. (1996). The NRK scheme maintains the complete square conservation property by adjusting the length of temporal integration steps and maintaining the same integral precision order as the original Runge–Kutta (RK) scheme.

The SCS was implemented in the grid-point atmospheric model of IAP LASG (GAMIL, Wang et al., 2004; Wang and Ji, 2006). GAMIL is widely used in climate simulation (Li et al., 2007, 2013; Wu and Li, 2008). The square conservation theory is rarely used on quasi-uniform grids or nonuniform grids because it is hard to construct a spatial discrete operator with an anti-symmetrical property on those grids.

In the past 2 decades, to avoid the polar singularity of the latitude–longitude grid, increasingly more atmosphere models have been built on the quasi-uniform grid, e.g., spectral transform methods (Swarztrauber, 1996), the finite-volume method (Lin, 2004; Putman and Lin, 2007; Chen and Xiao, 2008) and an extension of the finite-difference method to the generalized curvilinear coordinates (Toy and Nair, 2017).

Thuburn et al. (2009) and Ringler et al. (2010) provided a spatial discrete scheme based on arbitrarily structured C-grids, known as Thuburn–Ringler–Skamarock–Klemp (TRiSK). TRiSK is able to conserve total mass and total absolute vorticity, and total energy is conserved within time truncation error. These important properties allow for models using quasi-uniform Voronoi grids, the accuracies of which are similar to latitude–longitude grids (Weller et al., 2012). Based on Thuburn et al. (2009) and Ringler et al. (2010), a global/regional model, the Model for Prediction Across Scales (MPAS), was developed by the National Center for Atmospheric Research (NCAR) and Los Alamos National Laboratory (LANL; Skamarock et al., 2012, 2018).

Although the spatial discrete operator designed by Ringler et al. (2010) results in energy conservation, the total energy is still dissipative while using dissipative temporal integration schemes, i.e., Runge–Kutta; in other words, the conservation property of spatial discrete operator is not able to be maintain during temporal integration, which is referred to as conserving total energy in time truncation error. In this paper we attempt to construct a square conservative scheme for TRiSK, which is able to conserve total energy in the round-off error, but not just in time truncation error; this means that the variation of total energy should be in the round-off error during the entire temporal integration period. We call this complete energy conservation.

Total energy will be completely conserved only if the spatial discrete operator is anti-symmetrical and the temporal integration scheme is square conservative (Wang and Ji, 2006). The main obstacle of extending square conservation to the quasi-uniform grids is constructing the anti-symmetrical spatial discrete operator. Because many quasi-uniform grids are unstructured and the shapes of cells are not uniform, it is difficult to find the next or previous cell. In this study, we use the instantaneous energy conservation property of TRiSK to overcome the challenge of constructing an anti-symmetrical spatial operator on a quasi-uniform grid. After using NRK as a temporal integration scheme, the square conservative constraints are satisfied for both spatial and temporal directions, and the total energy, total mass and total absolute vorticity are completely conserved during the integration.

This paper is presented as follows: in Sect. 2, we review the TRiSK framework described by Ringler et al. (2010). Section 3 describes the method of extending square conservation to TRiSK in three parts. The first part presents a review of square conservation and a demonstration of the equivalent relationship between square conservation and energy conservation in a continuous shallow water system. In the second part of Sect. 3, we obtain the anti-symmetrical spatial discrete operator by using the derivative rule and the energy conservation property of TRiSK, a method that is key to extending square conservation to TRiSK. In the last part of Sect. 3, we review a new type of Runge–Kutta with fourth-order precision, which was developed by Wang et al. (1996), as the square conservative temporal integration scheme. In Sect. 4, by using the square conservation scheme, we demonstrate that the total mass and total absolute vorticity remain perfectly conservative. Section 5 exhibits the results of three different numerical tests, including the second, fifth and sixth test cases mentioned by Williamson et al. (1992).

2 Introduction of TRiSK
Back to toptop

The shallow water equation set may be written in a vector-invariant form as follows:

(1)ut-ξak×u+E=0,(2)ϕt+(ϕu)=0,

where ξa=ξ+f denotes the absolute vorticity; ξ=×u represents the relative vorticity; f=2Ωsin θ is the Coriolis parameter; E=K+g(h+hs)=K+ϕ+ϕs, ϕ=gh is the geopotential depth of the fluid; ϕs=ghs is the geopotential height of the underlying surface; ϕt=ϕ+ϕs is the free surface (top) of the fluid; K=|u|22 is the kinetic energy; u is the velocity vector; h and hs are the fluid thickness and surface height, respectively; θ represents the latitude; and g and Ω are acceleration of gravity and angular velocity of the earth.

In Ringler et al. (2010), the total energy is defined as

(3)ϵR10=hK+gh12h+hs.

To simplify the derivation in the following context, we define the total energy as

(4)Ωϵds=ΩgϵR10ds=ΩϕK+12ϕ2+ϕϕsds=ϕK+12ϕ2+ϕϕs,

where =(,) denotes the 2-norm. The inner product (,) is defined as

(5)(X,Y)=ΩXYds,

where Ω is the entire spherical surface.

Per the description provided in Ringler et al. (2010), velocity points are on the edges of each cell, the mass and kinetic energy points are in the center of each cell and vorticity points are on the vertices of each cell. The shallow water equation set may be expressed as a discrete form:

(6)uet-Qe+[E]e=0,(7)ϕit+[(ϕu)]i=0,

where ue, ϕi are the normal velocity and geopotential height. The subscript “e” signifies that the variable is defined on edge e; the subscript i signifies that the variable is defined at the center of ith cell. Qe is the absolute vorticity flux on the tangent direction  of the edge e, which is computed according to Eq. (49) in Ringler et al. (2010). Another two functions in Eqs. (6) and (7) are given by the introduction of Fig. 3 in Ringler et al. (2010):

(8)[E]e=1deiCE(e)-ne,iEi,(9)[(ϕu)]i=1AieEC(i)ne,ileϕeue,

where ne,i is an indicator function, defined as ne,i=1 when ne is an outward normal vector of cell i, and ne,i=-1 when ne is an inward normal vector of cell i; le is the length of edge e; i∈CE(e) denotes the two cells that share edge e; and e∈EC(i) is the set of edges that define the boundary of cell i. The potential vorticity on edge qe may be computed by the midpoint method (Ringler et al., 2010, Eq. 50) or the linear interpolation method (Weller, 2012, Eq. 5). The details are presented in Fig. 1.

https://www.geosci-model-dev.net/13/581/2020/gmd-13-581-2020-f01

Figure 1Definition of elements in a discrete system. Blue arrows represent the indicator function ne,i and red arrows are the indicator function te,v. The uniform grid is here for clearly introducing SCVT meshes, the situation on non-uniform grid is similar.

Download

3 Extending the square conservation to TRiSK
Back to toptop

As mentioned in Sect. 1, to obtain the complete square conservation property, the spatial discrete operator must be anti-symmetrical, and the temporal integration scheme is square conservative. Therefore, in this section, the method of extending the square conservation to TRiSK is introduced in three parts. Section 3.1 reviews the concept of square conservation, demonstrating the equivalent relationship between the square conservation and energy conservation. Section 3.2 constructs the anti-symmetrical spatial discrete operator. Section 3.3 introduces the square conservative temporal integration scheme by reviewing a class of NRK, developed by Wang et al. (1996).

3.1 Relationship between square conservation and energy conservation

First, we review the concept of anti-symmetrical operators and square conservation according to the study of Wang et al. (1996), considering the nonlinear evolution equation in operator form:

(10)Ft+L(F)=0.

Definition. Suppose H is a complete inner product space on R and  is an HH operator; if  satisfies the following inner product equation

(11)(L(F),F)=0,

then  is termed an anti-symmetrical operator.

The result of multiplying F on both sides of Eq. (10) and integrating globally is the square conservation property:

(12)ddtF2=0.

Next, we begin determining the relationship between energy conservation and square conservation. In the TRiSK framework, the evolution variables are u and ϕ.

The unified unit of evolution variables is the prerequisite of constructing the square conservation system. The evolution variables are unified by IAP transformation (Zeng and Zhang, 1987; Wang et al., 2004), and the original evolution variable u is replaced by the new evolution variable U=ϕu, after completing IAP transformation.

(13)F=Uϕ=ϕuϕ.

The physical significance of ϕ is the phase speed of the external gravity wave, and the shallow water equation set may be rewritten as a vector format:

(14)Ft+L(F)=tUϕ+LUϕ=0.

As defined in Sect. 2, ϕt=ϕ+ϕs, and

(15)ϕtt=ϕt+ϕst.

The surface height is determined to be independent of time,

(16)ϕst=0.

Therefore,

(17)ϕtt=ϕt.

Defining Ft=Uϕt, according to Eqs. (13) and (17), we have

(18)Ftt=Ft.

Multiplying Eq. (18) by Ft on both sides, and integrating globally, we have

(19)ddt12Ft2=Ft,Ft=ΩUUt+ϕ+ϕsϕtds=Ωt12|U|2+t12ϕ2+ϕϕsds=ddtϕK+12ϕ2+ϕϕs=dϵdt=0.

Accordingly, square conservation is equivalent to energy conservation in a continuous system.

3.2 Constructing the anti-symmetrical spatial discrete operator

In this subsection, we construct the anti-symmetrical spatial discrete operator using a specific combination of original operators in TRiSK. Firstly, we need to demonstrate the equivalent relationship between the square conservative spatial discrete operator and energy conservative spatial discrete operator in a continuous system; then we prove that this relationship is also applicable to discrete system.

Assuming a continuous-in-time system, the evolution equation of U is able to be expressed as

(20)Ut=ϕut+u2ϕϕt.

This formula is key to connecting square conservation and energy conservation; in the following theorem, we obtain an anti-symmetrical operator by using Eq. (20).

Theorem. The functions M=M(ϕ,u) and N=N(ϕ,u) satisfy

(21)ut+M=0ϕt+N=0,

and

(22)(M,ϕu)+(N,E)=0.

After IAP transformation Eq. (13), the evolution equation of U may be expressed as Eq. (20). Meanwhile, Eq. (21) may be rewritten as Eq. (14).

If the operator  satisfies Eq. (14), then  is an anti-symmetrical operator.

Proof.

(23)Ut=ϕut+u2ϕϕt=-ϕM-u2ϕN,

(24)(L(F),F)=-Ut,U-ϕt,ϕ=Ω-UUt-ϕϕtds=Ω-U-ϕM-u2ϕN+ϕNds=ΩϕuM+|u|22N+ϕNds=(M,ϕu)+(N,E)=0.

This theorem is proved in a continuous system, but the model is built as a discrete system; therefore, it is necessary to discuss the situation in a discrete space.

Following Ringler et al. (2010), we define the functions M=M(ϕ, u) and N=N(ϕ, u) in discrete space as follows:

(25)M=[E]e-Qe,(26)N=[(ϕu)]i.

And the semi-discrete shallow water equation set becomes

(27)ut+M=0,(28)ϕt+N=0.

Because the spatial discrete operator of TRiSK has an instantaneous energy conservation property, it is easy to prove (M,ϕu)+(N,E)=0 (details in Sect. 3.7.2 of Ringler et al., 2010).

There are cells, edges and vertices presented as three types of points on a spherical centroidal Voronoi tessellation (SCVT) grid, which is the mesh used by TRiSK. We define the inner product for different types of points as follows:

(29)(X,Y)i=i=1nCellsXiYiAi,(30)(X,Y)e=e=1nEdgesXeYeAe,

where XiYi are the variables in the cell; XeYe denote any variables on the edge; AiAe are the areas for each cell and edge; Ae=de×le, de is the distance between the two cells' centers on edge e; le is the length of edge e; “nCells” denotes the total cell number; and “nEdges” is the total edge number.

(31)(M,ϕu)e+(N,E)i=0.

Combining Eqs. (14) and (20), and rewriting it into a discrete system:

(32)Ft+L(F)=tUeϕi+LUeϕi=0,

where Ue=ϕeue, ϕe is the phase speed of external gravity wave on edge e. Note that we need to interpolate ϕi from cell center i to edge e; here we set ϕe=12iCE(e)ϕi.

(33)Uet=ϕeuet+ue2ϕeϕet=-ϕeM-ue2ϕeN.

As shown in Appendix A, we have the discrete anti-symmetrical operator L:

(34)(L(F),F)d=-U,Ute-ϕ,ϕti=0.

The subscript “d” represents that the inner product is computed in a discrete system.

Thus, the discrete evolution equation set becomes

(35)Uet+ϕeM+ue2ϕeN=0,(36)ϕit+N=0.

The model will be instantaneous square conservative by incorporating Eqs. (35) and (36) as the evolution equation set.

3.3 Constructing the temporal integration scheme with the square conservation property

The model is integrated in a discrete-in-time system; for the sake of guaranteeing complete square conservation, a square conservative temporal integration scheme is necessary. As NRK has the advantage of maintaining complete square conservation with a high order of integral precision, the fourth-order NRK scheme in TRiSK is adopted to obtain high integral precision and a long time step. To completely introduce the square conservative temporal integration scheme, we review some details of Wang et al. (1996).

The fourth-order NRK may be expressed as

(37)Fn+1=Fn+τnφFn,τ,

where τn is an adjustable time step and τ is the integral time step of the model.

(38)φFn,τ=τR1+2R2+2R3+R46.(39)R1=-LFnR2=-LFn+12τR1R3=-LFn+12τR2R4=-LFn+τR3.

Taking square operators on both sides of Eq. (37), delineating φn=φ(Fn,τ),

(40)Fn+12=Fn2+2τnφn,Fn+τn2φn2.

We note that although the spatial discrete operator L is anti-symmetrical, the total energy at the n+1 time point remains different from that at the n time point. Energy is able to be completely conserved by satisfying the following equation:

(41)Fn+12=Fn2.

Therefore,

(42)2τnφn,Fn+τn2φn2=0,(43)τn=-2φn,Fnφn2.

Considering the situation when τ→0, and according to Eqs. (17) and (18) in Wang et al. (1996),

(44)τn=τ3φn2R1,R2+R2,R3+R3,R4.

Once adopting the NRK scheme as the temporal integration scheme, the model will be completely square conservative, which implies that the total energy will be completely conserved from the beginning to the end of the integration. The NRK scheme is expected to perform better than RK in a numerical test. Moreover, NRK decays to RK by setting τn=τ.

While the integral time step is modified from τ to τn, the precision order of NRK is the same as RK; when constructing NRK based on the nth order RK, NRK has nth order precision as well, a conclusion proven by Theorem 1 in Wang et al. (1996).

4 Mass and absolute vorticity conservation
Back to toptop

In the CSCS introduced above, although the integral time step is modified from τ to τn, total mass and total absolute vorticity are nevertheless conserved. In the following demonstrations, we note that the mass conservation property and absolute vorticity conservation property are independent of temporal integration.

4.1 Mass conservation

Considering the total mass, multiplying Eq. (7) by Ai and summing all cells,

(45)i=1nCellsAiϕit=-i=1nCellsAi[(ϕu)]i=-i=1nCellseEC(i)ne,ileϕeue=-e=1nEdgesiCE(e)ne,ileϕeue=-e=1nEdgesleϕeue-leϕeue=0.

Note that the mass conservation property is independent of temporal integration.

4.2 Absolute vorticity conservation

According to Eq. (23) in Ringler et al. (2010), the relative vorticity is calculated according to the following diagnostic equation:

(46)ξ=1AveEV(v)te,vuede.

Multiplying by Av and summing all of the vertices yields

(47)v=1nVerticesAvξ=v=1nVerticeseEV(v)te,vuede=e=1nEdgesvVE(e)te,vuede=e=1nEdgesuede-uede=0,

where e∈EV(v) represents the set of edges that share the vertex v; v∈VE(e) are the two vertices on edge e. The indicator function te,v always points to the left side of ne,i. If k×ne,i is directed toward vertex v, then te,v=1; otherwise, te,v=-1. k is the unit vector, which points in the local vertical direction. See Fig. 1 for details. The total relative vorticity is shown to always be zero and independent of time.

Another method to compute the relative vorticity is to use the following prognostic equation, as described in Eq. (33) of Ringler et al. (2010):

(48)ξt+1AveEV(v)-te,vQede=0.

Multiplying the above equation by Av and summing all the vertices yields

(49)v=1nVerticesAvξt=v=1nVerticeseEV(v)te,vQede=e=1nEdgesvVE(e)te,vQede=e=1nEdgesQede-Qede=0.

Therefore, the relative vorticity is conserved during temporal integration.

Taking the partial derivative of the absolute vorticity with time yields

(50)ξat=ξt+ft.

The Coriolis parameter is independent of time, ft=0; thus

(51)v=1nVerticesAvξat=0.
5 Numerical tests
Back to toptop

To test the square conservation schemes using TRiSK, a new TRiSK-based shallow water dynamic core is developed, which is named TRiSK-based Multiple-Conservation dynamical cORE (TMCORE). The spatial discrete operators are the same as those introduced by Ringler et al. (2010), the evolution variable ue is replaced by Ue, as we described above, and the temporal integration scheme is selected from RK or NRK, both of which are of fourth-order precision.

We expected that the square conservation scheme will work on arbitrarily structured C-grids with a different initial field and mesh of a different resolution. In this section, we test the new scheme by using standard shallow water test cases 2, 5 and 6 (SWTC2, SWTC5, SWTC6) with two different meshes, as presented by Williamson et al. (1992). The first mesh has 2562 Voronoi cells (x1.2562), with an approximate resolution of 480 km, and the second mesh contains 40 962 Voronoi cells (x1.40962), with an approximate resolution of 120 km. The corresponding integral time steps to x1.2562 and x1.40962 are 900 and 360 s. Here, the midpoint scheme is selected as the method for interpolating the potential vorticity from vertices to edges for all tests.

In all of the test cases, we expect the complete energy conservation scheme (NRK) is able to conserve total mass, total absolute vorticity and total energy in the round-off error; meanwhile, it would be even better if NRK could provide us with a lower simulation error. Note, total energy is not merely conserved in time truncation error anymore; we need the change ratio of total energy to be limited to at least 10−14 magnitude.

5.1 Error measure methods

The global invariants error measure is as follows:

(52)IXn=SXin-SXi0SXi0,

where Xin is the variable at the nth time point on the ith cell and Xi0 is the variable at the initial time. The I function is the change ratio of the invariants.

The total mass error measure is as follows:

(53)Xin=hin,(54)Masschangeratio=Ihn.

The total energy error measure is as follows:

(55)Xin=ϵin,(56)Energychangeratio=Iϵn.

Measuring the fluid thickness error by L2 and L error norms is expressed as

(57)L2=Sfm(i)-fR(i)212SfR(i)212,(58)L=maxfm(i)-fR(i)maxfR(i),

where i denotes the index of each cell; fm(i) and fR(i), respectively, are the model solution and reference solution at the ith cell on the mesh; and the S function is the area-weighted accumulation of an arbitrary variable X.

(59)S(X)=i=1NX(i)A(i)i=1NA(i),

where A(i) is the area of the ith cell.

The reference solution should be an analytical solution or, when an analytical solution is not available, a high-resolution solution from the model with sufficient accuracy.

In the following context, NRK4 represents the NRK with fourth-order precision and RK4 represents the original Runge–Kutta scheme with fourth-order precision.

The differences of the error norms between NRK4 and RK4 schemes by using the different ratios of L2 (L2DR) and L (LInfDR) are expressed as follows:

(60)L2DR=L2CRK4-L2RK4L2RK4,(61)LInfDR=LCRK4-LRK4LRK4,

where L2NRK4 and L2RK4 are the L2 error norms of NRK4 and RK4, respectively, which is similar for LNRK4 and LRK4. NRK4 performs better than RK4 when the different ratios are negative; otherwise, NRK4 performs worse than RK4.

5.2 Global steady-state nonlinear zonal geostrophic flow (SWTC2)

For the global steady-state nonlinear zonal geostrophic flow test case, the initial velocity field has the following form

(62)u=u0cosθ,(63)v=0.

The geopotential height field is

(64)gh=gh0-aΩu0+u022sin2θ.

Here, we set Ω=7.292×10-5 s−1, g=9.80616 m s−2, a=6.37122×106 m, gh0=2.94×104 m2 s−2 and u0=2πa/(12 d), where h is fluid thickness and θ denotes latitude. In this test case, the exact solution is the initial state, and any difference between the numerical solution and the initial state is the simulation error.

In SWTC2, the true solution of ut, vt, ϕt is always zero; therefore, this test case can only represent the precision of spatial discrete operators but not the precision of temporal integration. This simulation is integrated for 10 years, but the shape of geostrophic flow breaks after 7 years. Therefore, we choose the simulation results from the first to the seventh year to compare the error norms of NRK4 and RK4.

Figure 2 measures the L2 and L error norms of geopotential height. In the first 4 years, the NRK4 and RK4 exhibit similar results, but in the last 3 years, the shape of geopotential flow tends to break. The error norms increase sharply after 6 years, and the differences between NRK4 and RK4 become more evident. Both the L2 and L error norms of NRK4 are evidently smaller than RK4, and the collapse of geopotential flow is delayed approximately 1 month when using NRK4.

https://www.geosci-model-dev.net/13/581/2020/gmd-13-581-2020-f02

Figure 2Geopotential height error norms of SWTC2. (a) L2 error norm; (b) L error norm. The results of RK4 and NRK4 are represented by blue and red lines, respectively. The model mesh is x1.2562.

Download

Figure 3 presents the variation of invariants as a function of time. The change ratio of total mass is limited to 10−15 magnitude, and total absolute vorticity oscillates around 10−20 magnitude, which means these two invariants are strictly conserved. The total energy of RK4 decreased approximately 0.5 % in the final year, but NRK4 maintains strict energy conservation (at 10−15 magnitude). Although the geopotential flow has been broken, NRK4 prevents an increasing rate of total potential enstrophy.

https://www.geosci-model-dev.net/13/581/2020/gmd-13-581-2020-f03

Figure 3The variation of integral invariants as a function of time of SWTC2. (a) Total mass change ratio; (b) total absolute vorticity; (c) total energy change ratio; (d) total potential enstrophy change ratio. The results of RK4 are represented by blue lines; the results of NRK4 are represented by red lines. The model mesh is x1.2562.

Download

5.3 Zonal flow over an isolated mountain (SWTC5)

SWTC5 is the fifth test case described by Williamson et al. (1992); the wind and height fields are similar to SWTC2, but h0=5960 m, u0=20 m s−1 and mountain height is determined according to the following equation:

(65)hs=hs01-rR,

where hs0=2000 m; R=π9; r=min[R2,(λ-λc)2+(θ-θc)2]; and λc and θc are the center longitude and latitude of the mountain, respectively. Here, we set λc=3π2 and θc=π6. As the analytical solution is not available, the reference solution is provided by a T511 idealized global spectral atmospheric model from GFDL, where 8×1012 m4 s−1 is selected as the coefficient for the 4 dissipation, and the test case is integrated for 50 d.

https://www.geosci-model-dev.net/13/581/2020/gmd-13-581-2020-f04

Figure 4The different ratios of fluid thickness error norms of SWTC5. (a) L2 error norm difference ratio; (b) L error norm difference ratio. The model mesh is x1.40962.

Download

Figure 4 presents the different ratios of error norms. In the first 35 d, the L2 and L error norms of NRK4 are considerably smaller than those of RK4. Compared with RK4, the L2 error norm of NRK4 decreases by approximately 2.5 % at the minimum point of L2DR, and the L error norm also decreases by approximately 3 % at the minimum point of LInfDR. The error norms increase very quickly after 35 d; therefore, the differences between error norm ratios for NRK4 and RK4 tend to be similar, along with time.

https://www.geosci-model-dev.net/13/581/2020/gmd-13-581-2020-f05

Figure 5Integral invariants of SWTC5. (a) Total mass change ratio; (b) total absolute vorticity; (c) total energy change ratio; (d) total potential enstrophy change ratio. The results of RK4 and NRK4 are represented by blue and red lines, respectively. The model mesh is x1.40962.

Download

Figure 5 presents the variation of the invariants as a function of time. The total mass and total absolute vorticity are completely conserved for both NRK4 and RK4. NRK4 is able to maintain strict energy conservation (at 10−15 magnitude) from the beginning to the end, but the total energy of RK4 is dissipative. The CSCS exhibits no influence on the total potential enstrophy.

https://www.geosci-model-dev.net/13/581/2020/gmd-13-581-2020-f06

Figure 6The different ratios of fluid thickness error norms of SWTC6. (a) L2 error norm difference ratio; (b) L error norm difference ratio. The model mesh is x1.40962.

Download

5.4 Rossby–Haurwitz wave (SWTC6)

The classical four-zonal wavenumber Rossby–Haurwitz wave was selected as the third test case. The initial condition follows Williamson et al. (1992). The initial state is the analytical solution of the nonlinear barotropic vorticity equation on the sphere but not the analytical solution of the shallow water equations. The reference field is provided by a T511 idealized global spectral atmospheric model from GFDL. To limit the noise of the spectral model, we use 5×1012 m4 s−1 as the coefficient for the 4 dissipation. As presented by Williamson et al. (1992), the phase speed of the Rossby–Haurwitz wave is calculated as follows:

(66)c=R(R+3)ω-2Ω(R+1)(R+2),

where R=4 is the zonal wavenumber of the Rossby–Haurwitz wave; ω=7.848×10-6 s−1; and Ω=7.292×10-5 s−1 is the rotation rate of the earth; therefore, the four-zonal wavenumber period of the Rossby–Haurwitz wave is approximately 29.52 d. We integrate the test case over one period (33 d).

In both simulations of NRK4 and RK4, the Rossby–Haurwitz wave begins to distort on the 25th day and then collapse a few days later.

Figure 6 presents the error norm difference ratios. NRK4 has a smaller L2 error norm than RK4 in the first 20 d. With growth of the L2 error norm, the difference between NRK4 and RK4 trends toward zero. On the fourth day, the L2 error norm of NRK4 is more than 0.11 % lower than that of RK4. NRK4 also has a smaller L error norm a majority of the time. On the sixth day, the L error norm of NRK4 is more than 0.08 % lower than that of RK4.

Figure 7 presents the variation of invariants as a function of time. The total mass and total absolute vorticity are strictly conserved for both NRK4 and RK4. As expected, NRK4 maintains strict energy conservation (at 10−15 magnitude), and RK4 cannot conserve the total energy during integration. With the Rossby–Haurwitz wave collapse, the total energy of RK4 rapidly dissipates after 25 d. There is no influence of NRK4 to potential enstrophy in this case.

The red lines and blue lines in Figs. 3c, 5c, and 7c are exactly the results of two kinds of energy conservation schemes; the blue lines present the property of conserving energy in time truncation error, and the red lines show conserving energy in the round-off error. The differences are clear to see. As mentioned above, conserving total energy in time truncation error leads to slight energy dissipation, and the dissipation accumulates during integration. Total energy may become zero after a long-term simulation, which is unreasonable for a pure dynamic core without any energy sinks or sources. On the other hand, the complete energy conservation scheme maintains strictly energy conservation in the entire integration period, even though it is not able to prevent the collapse of SWTC2, but the collapse time is delayed. Meanwhile, the simulation errors in SWTC5 and SWTC6 are reduced even in a short-term simulation.

6 Summary
Back to toptop

In this paper, we extend the CSCS to arbitrarily structured C-grids with shallow water equations, and we estimate the performance of the CSCS by using standard shallow water test cases.

There are two prerequisites for constructing CSCS: the anti-symmetrical spatial discrete operator and the square conservative temporal integration scheme. It is difficult to directly construct an anti-symmetrical spatial discrete operator on quasi-uniform grids; therefore, we take advantage of the instantaneous energy conservation property of the spatial discrete operators, as described by Ringler et al. (2010), to obtain the anti-symmetrical operator. After the IAP transformation, the units of evolution variables are unified, and the evolution variable ue is replaced with Ue=ϕeue. According to the derivative rule (Eq. 33), the temporal trend of Ue is expressed as a combination of the temporal trends of ue and ϕi, and we demonstrate that the spatial discrete operator of Ue is an anti-symmetrical operator. Then, we integrate the model with the square conservative temporal integration scheme NRK4, and the complete square conservation property is achieved.

https://www.geosci-model-dev.net/13/581/2020/gmd-13-581-2020-f07

Figure 7Integral invariants of SWTC6. (a) Total mass change ratio; (b) total absolute vorticity; (c) total energy change ratio; (d) total potential enstrophy change ratio. The results of RK4 and NRK4 are represented by blue and red lines, respectively. The model mesh is x1.40962.

Download

An important finding is the equivalency between the energy conservative operator and anti-symmetrical operator for both the continuous system and discrete system. In most previous studies, anti-symmetrical operators were constructed on uniform grids, especially longitude–latitude grids, and the equation's advection term was in the advection flux form. We extend the square conservation theory to a more general situation. The anti-symmetrical spatial discrete operator is constructed on quasi-uniform grids, and the equation is in the vector-invariant form.

The CSCS is able to maintain three integral invariants, including total mass, total absolute vorticity and total energy, in all the test cases, and the error norms decrease to varying degrees. The square conservation scheme improves the stability in SWTC2, and the error norms of NRK4 are evidently lower than RK4 after 4 years of simulation. For RK4, the total energy dissipates very quickly after the geostrophic flow collapse, but NRK4 maintains complete energy conservation for the entire period, and the increasing rate of the total potential enstrophy is also limited by the square conservation scheme. In both SWTC5 and SWTC6, NRK4 not only maintains strict conservation for three integral invariants; it also leads to less error norms than RK4.

Appendix A
Back to toptop

In this appendix, we attempt to demonstrate that the spatial discrete operator L is energy conservative. Our objective is to prove that the following equation is able to be satisfied by L:

(A1)U,Ute+ϕ,ϕti=0.

Consider the inner product on edge

(A2)U,Ute=e=1nEdgesUeUetAe.

Substituting Eq. (33) into the above formula,

(A3)U,Ute=e=1nEdgesUeCeuet+ue2CeϕetAe,

where Ce=ϕe is the phase speed of external gravity wave on edge e.

According to Eq. (52) in Ringler et al. (2010),

(A4)U,Ute=e=1nEdgesϕeueuet+ue24iCE(e)ϕitAe.

According to Eqs. (63) and (A8) in Ringler et al. (2010),

(A5)U,Ute=e=1nEdgesϕeueuetAe+i=1nCellsKiϕitAi.

Substituting Eq. (6) into the above formula,

(A6)U,Ute=e=1nEdgesϕeueAeQe+1deiCE(e)ne,iEi+i=1nCellsKiϕitAi.

According to Sect. 3.7.2 in Ringler et al. (2010),

(A7)e=1nEdgesϕeueAeQe=0.

Since Ae=lede,

(A8)U,Ute=e=1nEdgesϕeueleiCE(e)ne,iEi+i=1nCellsKiϕitAi

where Ei=Ki+ϕi.

According to Eq. (A4) in Ringler et al. (2010),

(A9)U,Ute=i=1nCellsEieEC(i)ne,iϕeuele+i=1nCellsKiϕitAi.

According to Eq. (7),

(A10)-Aiϕit=eEC(i)ne,ileϕeue.

Therefore,

(A11)U,Ute=-i=1nCellsEiϕitAi+i=1nCellsKiϕitAi.

Consider the inner product on cell

(A12)ϕ,ϕti=i=1nCellsϕiϕitAi.

Thus,

(A13)U,Ute+ϕ,ϕti=-i=1nCellsEiϕitAi+i=1nCellsKiϕitAi+i=1nCellsϕiϕitAi=0.
Code availability
Back to toptop
Code availability. 

Idealized Global Spectral Atmospheric Models (GFDL): https://www.gfdl.noaa.gov/idealized-spectral-models-quickstart/ (last access: 3 May 2019). TMCORE is available at https://github.com/TMCORE-Project/TMCORE (last access: 3 May 2019). The digital object identifier for Idealized Global Spectral Atmospheric Models (GFDL) with standard shallow water test cases is https://doi.org/10.5281/zenodo.3249878 (GFDL, 2007). The digital object identifier for TMCORE v1.0 is https://doi.org/10.5281/zenodo.3241647 (Zhou and Dong, 2019).

Author contributions
Back to toptop
Author contributions. 

LZ developed the theory, wrote the initial version of TMCORE, implemented the experiments and wrote the manuscript. JF, LH and LZ revised the context structure of the manuscript and gave some useful technical advice.

Competing interests
Back to toptop
Competing interests. 

The authors declare that they have no conflict of interest.

Acknowledgements
Back to toptop
Acknowledgements. 

This work was supported by the National Key R & D Program of China (2016YFA0600403) and the General Project of the National Natural Science Foundation of China (grant 41875134).

Financial support
Back to toptop
Financial support. 

This research has been supported by the National Key R & D Program of China (grant no. 2016YFA0600403) and the General Project of the National Natural Science Foundation of China (grant no. 41875134).

Review statement
Back to toptop
Review statement. 

This paper was edited by Simone Marras and reviewed by two anonymous referees.

References
Back to toptop

Arakawa, A.: Computational Design for Long-Term Numerical Integration of the Equations of Fluid Motion: Two-Dimensional Incompressible Flow. Part I, J. Comput. Phys., 1, 119–143, 1966. 

Arakawa, A. and Lam, V. R.: Computational Design of the Basic Dynamical Processes of the UCLA General Circulation Model, Meth. Comput. Phys., 17, 173–265, 1977. 

Chen, C., and Xiao, F.: Shallow water model on cubed-sphere by multi-moment finite volume method, J. Comput. Phys., 227, 5019–5044, https://doi.org/10.1016/j.jcp.2008.01.033, 2008. 

GFDL – Geophysical Fluid Dynamics Laboratory: Idealized Global Atmospheric Models with Spectral Dynamics, https://doi.org/10.5281/zenodo.3249878, 2007. 

Ji, Z. and Wang, B.: Further Discussion on the Construction and Application of Difference Scheme of Evolution Equations, Scient. Atmos. Sin., 15, 1–10, 1991. 

Li, L., Wang, B., Yuqing, W., and Hui, W.: Improvements in climate simulation with modifications to the Tiedtke convective parameterization in the grid-point atmospheric model of IAP LASG (GAMIL), Adv. Atmos. Sci., 24, 323–335, https://doi.org/10.1007/s00376-007-0323-3, 2007. 

Li, L., Wang, B., Dong, L., Liu, L., Shen, S., Hu, N., Sun, W., Wang, Y., Huang, W., Shi, X., Pu, Y., and Yang, G.: Evaluation of grid-point atmospheric model of IAP LASG version 2 (GAMIL2), Adv. Atmos. Sci., 30, 855–867, https://doi.org/10.1007/s00376-013-2157-5, 2013. 

Lin, S.-J.: A “Vertically Lagrangian” Finite-Volume Dynamical Core for Global Models, Mon. Weather Rev., 132, 2293–2307, 2004. 

Putman, W. M. and Lin, S.-J.: Finite-volume transport on various cubed-sphere grids, J. Comput. Phys., 227, 55–78, https://doi.org/10.1016/j.jcp.2007.07.022, 2007. 

Ringler, T. D., Thuburn, J., Klemp, J. B., and Skamarock, W. C.: A unified approach to energy conservation and potential vorticity dynamics for arbitrarily-structured C-grids, J. Comput. Phys., 229, 3065–3090, https://doi.org/10.1016/j.jcp.2009.12.007, 2010. 

Swarztrauber, P. N.: Spectral Transform Methods for Solving the Shallow-Water Equations on the Sphere, Mon. Weather Rev., 124, 730–744, 1996. 

Shen, T., Tian, Y., Chen, D., Shen, X., Sun, X., Liang, X., and Lu, W.: Numerical Weather Prediction, 2nd Edn., China Meteorological Press, Beijing, 487 pp., 2015. 

Skamarock, W. C., Klemp, J. B., Duda, M. G., Fowler, L. D., Park, S.-H., and Ringler, T. D.: A Multiscale Nonhydrostatic Atmospheric Model Using Centroidal Voronoi Tesselations and C-Grid Staggering, Mon. Weather Rev., 140, 3090–3105, https://doi.org/10.1175/mwr-d-11-00215.1, 2012. 

Skamarock, W. C., Duda, M. G., Ha, S., and Park, S.-H.: Limited-Area Atmospheric Modeling Using an Unstructured Mesh, Mon. Weather Rev., 146, 3445–3460, https://doi.org/10.1175/mwr-d-18-0155.1, 2018. 

Thuburn, J., Ringler, T. D., Skamarock, W. C., and Klemp, J. B.: Numerical representation of geostrophic modes on arbitrarily structured C-grids, J. Comput. Phys., 228, 8321–8335, https://doi.org/10.1016/j.jcp.2009.08.006, 2009. 

Toy, M. D. and Nair, R. D.: A Potential Enstrophy and Energy Conserving Scheme for the Shallow-Water Equations Extended to Generalized Curvilinear Coordinates, Mon. Weather Rev., 145, 751–772, https://doi.org/10.1175/mwr-d-16-0250.1, 2017. 

Wang, B.: Design of a new dynamical core for global atmospheric models based on some efficient numerical methods, Sci. China Ser. A, 47, 4–21, 2004. 

Wang, B.: An Explicit Multi-Conservation Finite-Difference Scheme for Shallow-Water Equation, J. Comput. Math., 26, 404–409, 2008. 

Wang, B. and Ji, Z.: An Explicit Complete Square Conservative Scheme with Adjustable Time Intervals, ACTA Meteorol. Sin., 8, 403–409, 1994a. 

Wang, B. and Ji, Z.: The Harmonious Dissipative Operators and the Completely Square Conservative Difference Scheme in an Explicit Way Science in China, Sci. China Ser. B, 37, 462–471, 1994b. 

Wang, B. and Ji, Z.: Construction and Numerical Test of the Multi-Conservation Difference Scheme, Chinese Sci. Bull., 48, 1016–1020, 2003. 

Wang, B. and Ji, Z.: The New Numerical Methods for Atmospheric Science with Applications, 1st Edn., Science Press, Beijing, 208 pp., 2006. 

Wang, B., Ji, Z., and Zeng, Q.: A Class of New Explicit Runge-Kutta Schemes, Prog. Nat. Sci., 6, 195–205, 1996. 

Weller, H.: Controlling the Computational Modes of the Arbitrarily Structured C Grid, Mon. Weather Rev., 140, 3220–3234, https://doi.org/10.1175/mwr-d-11-00221.1, 2012. 

Weller, H., Thuburn, J., and Cotter, C. J.: Computational Modes and Grid Imprinting on Five Quasi-Uniform Spherical C Grids, Mon. Weather Rev., 140, 2734–2755, https://doi.org/10.1175/mwr-d-11-00193.1, 2012. 

Williamson, D. L., Drake, J. B., Hack, J. J., Jakob, R., and Swarztrauber, P. N.: A Standard Test Set for Numerical Approximations of the Shallow Water Equations in Spherical Geometry, J. Comput. Phys., 102, 211–224, 1992.  

Wu, Z. and Li, J.: Prediction of the Asian-Australian Monsoon Interannual Variations with the Grid-Point Atmospheric Model of IAP LASG (GAMIL), Adv. Atmos. Sci., 25, 387–394, 2008. 

Zeng, Q. and Ji, Z.: On The Computational Stability of Evolution Equations, Math. Numer. Sin., 1, 79–86, 1981. 

Zeng, Q. and Zhang, X.: Available Energy Conserving Schemes for Spherical Baroclinic Primitive Equations, Scient. Atmos. Sin., 11, 113–127, 1987. 

Zhou, L. and Dong, L.: TRiSK-based Multiple-Conservation dynamical cORE, https://doi.org/10.5281/zenodo.3241647, 2019. 

Publications Copernicus
Download
Short summary
Atmospheric simulation should obey the physical conservation law: in a closed atmospheric system, some of the invariants should be conserved, i.e. the total mass, total energy, total absolute vorticity and so on. In this paper, we have improved the conservative properties of the atmospheric model which is based on an arbitrarily structured C-grid. The tests show that better conservative properties bring us less simulation error, and the stability of model is also improved.
Atmospheric simulation should obey the physical conservation law: in a closed atmospheric...
Citation