Extending Square Conservation to Arbitrarily Structured C-grids with Shallow Water Equations

The square conservation law is implemented in the 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 developed on quasi-uniform grids, such as arbitrarily structured Cgrids. 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, total absolute vorticity 15 and conserving total energy with time truncation error, the first 2 integral invariants are exactly conserved during integration, but the total energy dissipates when using the dissipative temporal integration schemes, i.e., Runge-Kutta. The method of strictly conserving the total energy simultaneously, which means conserving energy in round-off error during entire temporal integration period, uses both an anti-symmetrical spatial discrete operator and square conservative temporal integration scheme. In this study, we demonstrate that square conservation is equivalent to energy conservation in both a continuous shallow water 20 system and a discrete shallow water system of TRiSK, attempting 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 by Institute of Atmospheric Physics (IAP) transformation, the temporal derivates of new evolution variables can be expressed by a combination of temporal derivates of original evolution variables, which means the square conservative spatial discrete operator can be obtain by using original spatial discrete operator in TRiSK. 25 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.


Introduction
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, 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 5 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 the quadratic quantities. Toy and Nair (2017) developed an energy and potential enstrophy conserving scheme for 10 shallow water equations on generalized curvilinear coordinates, they mentioned conserving analogues of total energy and total potential enstrophy in numerical models are known to prevent a spurious cascade of energy toward small scales. For a shortterm 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, 1996). 15 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, 2013;Zeng, 1981). On a latitude-longitude grid, energy is able to 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 of which is better developed. Wang and Ji (1994a) discussed the square conservative scheme (SCS), the 20 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 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 25 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 then integrate the model with a square conservative temporal integration scheme, i.e., a modified predict-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 Runge-Kutta 30 schemes, hereafter NRK, 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 scheme, hereafter RK.
The SCS was implemented in the grid-point atmospheric model of IAP LASG (GAMIL, Wang et al. 2004, Wang andJi, 2006). GAMIL is widely used in climate simulation (Li et al., 2007(Li et al., , 2013Wu 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 antisymmetrical property on those grids.
In the most recent two decades, to avoid the polar singularity of the latitude-longitude grid, increasingly more atmosphere 5 models have been built on the quasi-uniform grid, i.e., spectral transform methods (Swarztrauber, 1996); the finite volume method (Lin, 2004;Putman and Lin, 2007;Chen and Xiao, 2008); and an extension on the finite difference method to the generalized curvilinear coordinates (Toy and Nair, 2017). Thuburn et al. (2009) andRingler 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 the total mass and total absolute vorticity, 10 and the total energy is conserved within time-truncation error. These important properties enable models using quasi-uniform Voronoi grids, the accuracies of which are similar to latitude-longitude grids . Based on Thuburn et al. (2009) andRingler 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(Skamarock et al., , 2018. 15 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, this is so-call 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 round-off error, but not just in time truncation error, which means that the variation of total energy should be in 20 round-off error during 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 25 conservation property of TRiSK to overcome the challenge of constructing an anti-symmetrical spatial operator on a quasiuniform grid. After using NRK as a temporal integration scheme, the square conservative constrains 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 section 2, we review the TRiSK framework which was described by Ringler et al. (2010). 30 Section 3 describes the method of extending square conservation to TRiSK in 3 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 section 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 section 3, we review a new type of Runge-Kutta with 4th-order precision, which was developed by Wang et al. (1996) as the square conservative temporal integration scheme. In section 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 2nd, 5th and 6th test cases mentioned by Williamson (1992).

Introduction of TRiSK 5
The shallow water equation set may be written in a vector-invariant form as follows: where, = + denotes the absolute vorticity; = ∇ × represents the relative vorticity; = 2 is the Coriolis parameter; = + (ℎ + ℎ ) = + + , = ℎ is the geopotential depth of the fluid; = ℎ is the geopotential 10 height of the underlying surface; = + is the free surface (top) of the fluid; = | | 2 2 is the kinetic energy; is the velocity vector; ℎ and ℎ are the fluid thickness and surface height, respectively; represents the latitude; and 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 ‖•‖ = √(•,•) denotes the 2-norm. The inner product (•,•) is defined as where Ω is the entire spherical surface.
where , are the normal velocity and geopotential height. The subscript signifies that the variable is defined on edge ; the subscript signifies that the variable is defined at the center of th cell. ⊥ is the absolute vorticity flux on the tangent 10 direction ⊥ of the edge , which is computed according to Eq. (49) in Ringler et al. (2010).
where , is an indicator function, defined as , = 1 when is an outward normal vector of cell , and , = −1 when is an inward normal vector of cell ; is the length of edge ; ∈ ( ) denotes the two cells that share edge ; and ∈ ( ) is the set of edges that define the boundary of cell . The potential vorticity on edge 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 5 Figure 1.

Extending the square conservation to TRiSK
As mentioned in section 1, to obtain the complete square conservation property, the spatial discrete operator must be antisymmetrical, 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. Subsection 3.1 reviews the concept of square conservation, 10 demonstrating the equivalent relationship between the square conservation and energy conservation. Subsection 3.2 constructs the anti-symmetrical spatial discrete operator. Subsection 3.3 introduces the square conservative temporal integration scheme by reviewing a class of Runge-Kutta (NRK), which was developed by Wang et al. (1996).

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. 15 (1996), considering the nonlinear evolution equation in operator form: Definition. Suppose H is a complete inner product space on R and ℒ is an → operator; if ℒ satisfies the following inner product equation 20 then ℒ is termed an anti-symmetrical operator.
The result of multiplying on both sides of (6) and integrating globally is the square conservation property: Next, we begin determine the relationship between energy conservation and square conservation. In the TRiSK framework, the evolution variables are and . 7 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 is replaced by the new evolution variable = √ , after completing IAP transformation.
The physical significance of √ is the phase speed of the external-gravity wave, and the shallow water equation set may be 5 rewritten as a vector format: As defined in section 2, = +

= +
The surface height is determined to be independent of time, 10 = 0 Therefore, Defining = ( ), according to (9) and (11), we have Accordingly, square conservation is equivalent to energy conservation in a continuous system.

Constructing the anti-symmetrical spatial discrete operator
In this subsection, we construct the anti-symmetrical spatial discrete operator by a specific combination of original operators in TRiSK. Firstly, we need to demonstrate the equivalent relationship between square conservative spatial discrete operator and energy conservative spatial discrete operator in continuous system, then prove that this relationship is also apply to discrete system. 5 Assuming a continuous-in-time system, the evolution equation of is able to be expressed as This formula is key to connecting square conservation and energy conservation; it is difficult to directly construct the antisymmetrical operator on quasi-uniform grids.
This theorem is proved in a continuous system, but the model is built in a discrete system; therefore, it is necessary to discuss the situation in a discrete system. 25 Following Ringler et al. (2010), we set the discrete functions = ( , ) and = ( , ) as: 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 5 ( , ) + ( , ) = 0. (Details in Ringler et al. (2010), section 3.7.2) 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: where , are the variables in the cell; , denote any variables on the edge; , are the areas for each cell and edge; = × , is the distance between the two cells' centers on edge ; is the length of edge ; denotes the total cell number; and is the total edge number.
As shown in the Appendix A, we have the discrete anti-symmetrical operator 20 The subscript 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 (21) and (22) as the evolution equation set.

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 4 th -order NRK scheme in TRiSK is adopted to obtain high integral 5 precision and a long-time step. To completely introduce the square conservative temporal integration scheme, we review some details in Wang et al. (1996).
The 4 th -order NRK may be expressed as where is an adjustable time step and is the integral time step of the model. 10 ‖ +1 ‖ 2 = ‖ ‖ 2 + 2 ( , ) + 2 ‖ ‖ 2 , We notice that although the spatial discrete operator is anti-symmetrical, the total energy at the + 1 time point remains 15 different from that at the time point. Energy is able to be completely conserved by satisfying the following equation: Considering the fitness when → 0, as described in Eqs. (17)-(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 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 = . 25 While the integral time step is modified from to , the precision order of NRK is the same as RK, when constructing NRK based on the th order RK, NRK has th order precision either, 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 , the total mass and total absolute vorticity are nevertheless conserved. In the following demonstrations, we notice that the mass conservation property and absolute vorticity conservation property are independent of temporal integration.

Mass Conservation 5
Considering the total mass, multiplying (5) by and summing all cells, Notice that the mass conservation property is independent of temporal integration.

Absolute Vorticity Conservation 10
According to Ringler et al. (2010) formula (23) is the unit vector, which points in the local vertical direction. See Figure 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 Ringler et al. 20 (2010) Eq. (33)

Multiplying the above equation by and summing all the vertices yields
Therefore, the relative vorticity is conserved during temporal integration. 25 Taking the partial derivative of the absolute vorticity with time yields

Numerical Tests
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 is replaced by , as we described above, and the temporal 5 integration scheme is selected from RK or NRK, both of which are in 4 th -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 (1992) 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 round-off error, meanwhile, it would be even better if NRK can bring us less simulation error. 15 Note, total energy is not merely conserved in time truncation error anymore, we need the change ratio of total energy to be limited in at least 10 −14 magnitude.

Error measure methods
Global invariants error measure: is the variable at the th time point on the ith cell and 0 is the variable at the initial time. The function is the change ratio of the invariants.
The total mass error measure: The total energy error measure: 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 4 th -order precision and RK4 represents the original Runge-Kutta scheme with 4 th -order precision. 10 The differences of the error norms between NRK4 and RK4 schemes by using the different ratios of 2 (L2DR) and ∞ where 2 4 and 2 4 are the 2 error norms of NRK4 and RK4, respectively, which is similar for ∞ 4 and ∞ 4 . 15 NRK4 has better performance than RK4 when the different ratios are negative; otherwise, NRK4 has worse performance than RK4.
In SWTC2, the true solution of , , 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 1 st to the 7 th year to compare the error norms of NRK4 and RK4. Figure 2 measures the 2 and ∞ error norms of geopotential height. In the first 4 years, the NRK4 and RK4 exhibit similar 5 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 2 and ∞ 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 in 10 −15 magnitude, and total absolute vorticity is oscillating around 10 −20 magnitude, which means these two invariants are strictly conserved. 10 The total energy of RK4 decreased approximately 0.5% in the final year, but NRK4 maintains strict energy conservation (in 10 −15 magnitude). Although the geopotential flow has been broken, NRK4 prevents an increasing rate of total potential enstrophy.

Zonal Flow Over an Isolated Mountain (SWTC5)
SWTC5 is the 5 th test case described by Williamson 1992; the wind and height fields are similar to SWTC2, but ℎ 0 = 5960 , 0 = 20 / and mountain height is determined according to the following equation: where ℎ 0 = 2000 ; = 9 ; = √ [ 2 , ( − ) 2 + ( − ) 2 ]; and and are the center longitude and latitude of 5 the mountain, respectively. Here, we set = 3 2 and = 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 4 / is selected as the coefficient for the ∇ 4 dissipation, and the test case is integrated for 50 days. (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. Figure 4 presents the different ratios of error norms. In the first 35 days, the 2 and ∞ error norms of NRK4 are considerably smaller than those of RK4. Compared with RK4, the 2 error norm of NRK4 decreases by 5 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 days; 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 (in 10 −15 magnitude) from the 10 beginning to the end, but the total energy of RK4 is dissipative. The CSCS exhibits no influence on the total potential enstrophy.

Rossby-Haurwitz Wave (SWTC6)
The classical 4 zonal wavenumber Rossby-Haurwitz wave was selected as the third test case. The initial condition follows Williamson (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 15 atmospheric model from GFDL. To limit the noise of the spectral model, we use 5 × 10 12 m 4 /s as the coefficient for the ∇ 4 dissipation. As presented by Williamson, 1992, the phase speed of the Rossby-Haurwitz wave is calculated as follows: where = 4 is the zonal wavenumber of the Rossby-Haurwitz wave; = 7.848 × 10 −6 −1 ; and Ω = 7.292 × 10 −5 −1 is the rotation rate of the earth; therefore, the 4 zonal wavenumber period of the Rossby-Haurwitz wave is approximately 29.52 days. We integrate the test case over one period (33 days). In both simulations of NRK4 and RK4, the Rossby-Haurwitz wave begins to distort at the 25 th day and then collapse a few days later. Figure 6 presents the error norm difference ratios. NRK4 has a smaller 2 error norm than RK4 in the first 20 days. With growth of the 2 error norm, the difference between NRK4 and RK4 trends toward zero. At the 4 th day, the 2 error norm of NRK4 is more than 0.11% less than that of RK4. NRK4 also has a smaller ∞ error norm a majority of the time. At the 6 th day, 10 the ∞ error norm of NRK4 is more than 0.08% less 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 (in 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 days. There is no influence of NRK4 to potential enstrophy in this case. 15 the blue lines present the property of conserving energy in time truncation error, and the red lines show conserving energy in round-off error, it's clear to see differences. As we mentioned above, conserving total energy in time truncation error leads to slightly energy dissipative, but the dissipate accumulates during integration, total energy may become zero after a long-term simulation, which is unreasonable for a pure dynamic core without any energy sink and source. On the other hand, complete 5 energy conservation scheme maintains strictly energy conservation in entire integration period, even though it is not able to prevent the collapse of SWTC2, the collapse time is delayed, meanwhile, the simulation errors in SWTC5 and SWTC6 are reduced even if in a short-term simulation.

Summary
In this paper, we extend the CSCS to arbitrarily structured C-grids with shallow water equations, and we estimate the 10 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 15 variables are unified, and the evolution variable is replaced with = √ . According to the derivative rule (19), the temporal trend of is expressed as a combination of the temporal trends of and , and we demonstrate that the spatial discrete operator of 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 20 continuous system and discrete system. In most of previous study, 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 quasiuniform 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 25 test cases, and the error norms decrease in varying degrees. The square conservation scheme improves the stability in SWTC2, and the error norms of NRK4 are evidently less 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 but also leads to less error norms than RK4.
In this appendix, we attempt to demonstrate that the spatial discrete operator is energy conservative. Our objective is to prove that the following equation is able to be satisfied by , Competing interests. The authors declare that they have no conflict of interest.