**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

Lilong Zhou Jinming Feng Lijuan Hua and Linhao Zhong

^{1,3},

^{2},

^{4},

^{2}

**Lilong Zhou et al.**Lilong Zhou Jinming Feng Lijuan Hua and Linhao Zhong

^{1,3},

^{2},

^{4},

^{2}

^{1}College of Earth and Planetary Sciences, University of Chinese Academy of Sciences, Beijing, China^{2}Key Laboratory of Regional Climate-Environment for Temperate East Asia, Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing, China^{3}Numerical Weather Prediction Center of China Meteorological Administration, Beijing, China^{4}Key Laboratory of Computational Geodynamics, University of Chinese Academy of Sciences, Beijing, China

^{1}College of Earth and Planetary Sciences, University of Chinese Academy of Sciences, Beijing, China^{2}Key Laboratory of Regional Climate-Environment for Temperate East Asia, Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing, China^{3}Numerical Weather Prediction Center of China Meteorological Administration, Beijing, China^{4}Key Laboratory of Computational Geodynamics, University of Chinese Academy of Sciences, Beijing, China

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

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

Received: 02 May 2019 – Discussion started: 21 Jun 2019 – Revised: 23 Sep 2019 – Accepted: 13 Jan 2020 – Published: 14 Feb 2020

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.

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 10^{8} 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).

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

where ${\mathit{\xi}}_{\mathrm{a}}=\mathit{\xi}+f$ denotes the absolute vorticity; $\mathit{\xi}=\mathrm{\nabla}\times \mathit{u}$ represents the relative
vorticity; *f*=2Ωsin *θ* is the Coriolis parameter; $E=K+g(h+{h}_{\mathrm{s}})=K+\mathit{\varphi}+{\mathit{\varphi}}_{\mathrm{s}}$, *ϕ*=*g**h* is the geopotential depth of the fluid; *ϕ*_{s}=*g**h*_{s} is the geopotential height of the underlying surface; ${\mathit{\varphi}}_{\mathrm{t}}=\mathit{\varphi}+{\mathit{\varphi}}_{\mathrm{s}}$ is the free surface (top) of the fluid; $K=\frac{|\mathit{u}{|}^{\mathrm{2}}}{\mathrm{2}}$ is the kinetic energy; ** u** is the velocity vector;

*h*and

*h*

_{s}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

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

where $\Vert \cdot \Vert =\sqrt{(\cdot ,\cdot )}$ denotes the 2-norm. The inner product ($\cdot ,\cdot $) is defined as

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:

where *u*_{e}, *ϕ*_{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 *i*th cell. ${Q}_{e}^{\perp}$ 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):

where *n*_{e,i} is an indicator function, defined as
${n}_{e,i}=\mathrm{1}$ when *n*_{e} is an outward normal vector of
cell *i*, and ${n}_{e,i}=-\mathrm{1}$ when *n*_{e} is an inward
normal vector of cell *i*; *l*_{e} 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 *q*_{e} 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.

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:

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

*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:

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 $\mathit{U}=\sqrt{\mathit{\varphi}}\mathit{u}$, after completing IAP transformation.

The physical significance of $\sqrt{\mathit{\varphi}}$ is the phase speed of the external gravity wave, and the shallow water equation set may be rewritten as a vector format:

As defined in Sect. 2, ${\mathit{\varphi}}_{\mathrm{t}}=\mathit{\varphi}+{\mathit{\varphi}}_{\mathrm{s}}$, and

The surface height is determined to be independent of time,

Therefore,

Defining ${\mathit{F}}_{\mathrm{t}}=\left(\begin{array}{l}\mathit{U}\\ {\mathit{\varphi}}_{\mathrm{t}}\end{array}\right)$, according to Eqs. (13) and (17), we have

Multiplying Eq. (18) by *F*_{t} on both sides, and integrating globally, we have

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

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* $\mathcal{M}=\mathcal{M}(\mathit{\varphi},u)$ *and* $\mathcal{N}=\mathcal{N}(\mathit{\varphi},u)$ *satisfy*

*and*

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

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*(

*ϕ*,

**) in discrete space as follows:**

*u*And the semi-discrete shallow water equation set becomes

Because the spatial discrete operator of TRiSK has an instantaneous energy conservation property, it is easy to prove $(M,\mathit{\varphi}u)+(N,E)=\mathrm{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:

where *X*_{i}*Y*_{i} are the variables in the cell; *X*_{e}*Y*_{e} denote any variables on the edge; *A*_{i}*A*_{e} are the areas for each cell and edge;
${A}_{e}={d}_{e}\times {l}_{e}$, *d*_{e} is the distance between the two cells' centers on edge *e*; *l*_{e} is the length of edge *e*; “nCells” denotes the total cell number; and “nEdges” is the total edge number.

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

where ${U}_{e}=\sqrt{{\mathit{\varphi}}_{e}}{u}_{e}$, $\sqrt{{\mathit{\varphi}}_{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 ${\mathit{\varphi}}_{e}=\frac{\mathrm{1}}{\mathrm{2}}\sum _{i\in \mathrm{CE}\left(e\right)}{\mathit{\varphi}}_{i}$.

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

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

Thus, the discrete evolution equation set becomes

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

where *τ*_{n} is an adjustable time step and *τ* is the integral time step of the model.

Taking square operators on both sides of Eq. (37), delineating ${\mathit{\phi}}^{n}=\mathit{\phi}({\mathit{F}}^{n},\mathit{\tau})$,

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:

Therefore,

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

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
*n*th order RK, NRK has *n*th order precision as well, a conclusion proven by Theorem 1 in Wang et al. (1996).

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 *A*_{i} and summing all
cells,

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:

Multiplying by *A*_{v} and summing all of the vertices yields

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 *t*_{e,v} always points to the left side of *n*_{e,i}. If $\mathit{k}\times {n}_{e,i}$ is directed toward vertex *v*, then ${t}_{e,v}=\mathrm{1}$; otherwise, ${t}_{e,v}=-\mathrm{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):

Multiplying the above equation by *A*_{v} and summing all the vertices
yields

Therefore, the relative vorticity is conserved during temporal integration.

Taking the partial derivative of the absolute vorticity with time yields

The Coriolis parameter is independent of time, $\frac{\partial f}{\partial t}=\mathrm{0}$; thus

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 *u*_{e} is replaced by *U*_{e}, 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:

where ${X}_{i}^{n}$ is the variable at the *n*th time point on the *i*th cell and ${X}_{i}^{\mathrm{0}}$ 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:

The total energy error measure is as follows:

Measuring the fluid thickness error by *L*_{2} and *L*_{∞} error norms
is expressed as

where *i* denotes the index of each cell; *f*_{m}(*i*) and *f*_{R}(*i*), respectively, are the model solution and reference solution at the *i*th cell on the mesh; and the *S* function is the area-weighted accumulation of an arbitrary variable *X*.

where *A*(*i*) is the area of the *i*th 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 *L*_{2} (L2DR) and *L*_{∞} (LInfDR) are expressed
as follows:

where ${L}_{{\mathrm{2}}_{\mathrm{NRK}\mathrm{4}}}$ and ${L}_{{\mathrm{2}}_{\mathrm{RK}\mathrm{4}}}$ are the *L*_{2} error norms of NRK4 and RK4, respectively, which is similar for ${L}_{{\mathrm{\infty}}_{\mathrm{NRK}\mathrm{4}}}$ and ${L}_{{\mathrm{\infty}}_{\mathrm{RK}\mathrm{4}}}$. 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

The geopotential height field is

Here, we set $\mathrm{\Omega}=\mathrm{7.292}\times {\mathrm{10}}^{-\mathrm{5}}$ s^{−1}, *g*=9.80616 m s^{−2},
$a=\mathrm{6.37122}\times {\mathrm{10}}^{\mathrm{6}}$ m, gh${}_{\mathrm{0}}=\mathrm{2.94}\times {\mathrm{10}}^{\mathrm{4}}$ m^{2} s^{−2} and ${u}_{\mathrm{0}}=\mathrm{2}\mathit{\pi}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 $\frac{\partial u}{\partial t}$, $\frac{\partial v}{\partial t}$, $\frac{\partial \mathit{\varphi}}{\partial 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 *L*_{2} 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 *L*_{2} 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.

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.

## 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 *h*_{0}=5960 m, *u*_{0}=20 m s^{−1} and mountain height is determined according to the following equation:

where *h*_{s0}=2000 m; $R=\frac{\mathit{\pi}}{\mathrm{9}}$; $r=\sqrt{\mathrm{min}[{R}^{\mathrm{2}},(\mathit{\lambda}-{\mathit{\lambda}}_{\mathrm{c}}{)}^{\mathrm{2}}+(\mathit{\theta}-{\mathit{\theta}}_{\mathrm{c}}{)}^{\mathrm{2}}]}$; and *λ*_{c} and *θ*_{c} are the center longitude and latitude of the mountain, respectively. Here, we set ${\mathit{\lambda}}_{\mathrm{c}}=\frac{\mathrm{3}\mathit{\pi}}{\mathrm{2}}$ and ${\mathit{\theta}}_{\mathrm{c}}=\frac{\mathit{\pi}}{\mathrm{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×10^{12} m^{4} s^{−1} is selected as the coefficient for the ∇^{4} dissipation, and the test case is integrated for 50 d.

Figure 4 presents the different ratios of error norms. In the first 35 d, the *L*_{2} and *L*_{∞} error norms of NRK4 are considerably smaller than those of RK4. Compared with RK4, the *L*_{2} 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.

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.

## 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×10^{12} m^{4} 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:

where *R*=4 is the zonal wavenumber of the Rossby–Haurwitz wave; $\mathit{\omega}=\mathrm{7.848}\times {\mathrm{10}}^{-\mathrm{6}}$ s^{−1}; and $\mathrm{\Omega}=\mathrm{7.292}\times {\mathrm{10}}^{-\mathrm{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 *L*_{2} error norm than RK4 in the first 20 d. With growth of the *L*_{2} error norm, the difference between NRK4 and RK4 trends toward zero. On the fourth day, the *L*_{2} 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.

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 *u*_{e} is
replaced with ${U}_{e}=\sqrt{{\mathit{\varphi}}_{e}}{u}_{e}$. According to the derivative rule (Eq. 33), the temporal trend of *U*_{e} is expressed as a combination of the temporal trends of *u*_{e} and *ϕ*_{i}, and we demonstrate that the spatial discrete operator of *U*_{e} 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.

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.

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*:

Consider the inner product on edge

Substituting Eq. (33) into the above formula,

where ${C}_{e}=\sqrt{{\mathit{\varphi}}_{e}}$ is the phase speed of external gravity wave on edge *e*.

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

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

Substituting Eq. (6) into the above formula,

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

Since *A*_{e}=*l*_{e}*d*_{e},

where ${E}_{i}={K}_{i}+{\mathit{\varphi}}_{i}$.

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

According to Eq. (7),

Therefore,

Consider the inner product on cell

Thus,

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

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.

The authors declare that they have no conflict of interest.

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

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

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

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.