Generalized Z-Grid Model for Numerical Weather Prediction

Abstract

Z-grid finite volume models conserve all-scalar quantities as well as energy and potential enstrophy and yield better dispersion relations for shallow water equations than other finite volume models, such as C-grid and C-D grid models; however, they are more expensive to implement. During each time integration, a Z-grid model must solve Poisson equations to convert its vorticity and divergence to a stream function and velocity potential, respectively. To optimally utilize these conversions, we propose a model in which the stability and possibly accuracy on the sphere are improved by introducing more stencils, such that a generalized Z-grid model can utilize longer time-integration steps and reduce computing time. Further, we analyzed the proposed model’s dispersion relation and compared it to that of the original Z-grid model for a linearly rotating shallow water equation, an important property for numerical models solving primitive equations. The analysis results suggest a means of balancing stability and dispersion. Our numerical results also show that the proposed Z-grid model for a shallow water equation is more stable and efficient than the original Z-grid model, increasing the time steps by more than 1.4 times.

1. Introduction

The societal demand for methods that facilitate future accurate weather prediction is pushing resolvable scales of global weather forecasts down to the kilometer range. Compared to other schemes of spectral and finite element methods, the use of finite volume methods has several advantages in terms of grid nesting, handling non-hydrostatic phenomena, flexibility because of parallel computation, etc. Many meteorological operation centers are developing their next generation numerical prediction models and most of them are based on a finite volume scheme, for example, NOAA, Met Office, and ECMWF (European Centre for Medium-Range Weather Forecasts). They choose different global grid structures and forms of the finite element methods. These finite element methods can be classified into several types of models according to the grid arrangements of the model state variables: A-grid, C-grid, and C-D grid [1,2,3,4,5] as well as Z-grid (Heikes and Randall [6]). These models place mass fields at the centers of grid cells, for example, Voronoi or centroid centers. With the exception of Z-grid, all other methods use velocity vector equations in their prediction models. A-grid models place velocity vectors at the cell centers along with their mass fields; C-grid models place normal velocity vectors at the centers of cell edges; and C-D grid models place both tangential and normal velocity vectors at the centers of cell edges. These models in the velocity vector form are relatively easy to implement and run efficiently even though vectors are not easy to process on a sphere. For example, the six edges of a hexagonal cell of an icosahedral grid on a sphere may not be on the same plane and the velocity vectors at these edges may be on different planes. This issue affects some of the elegant theoretical conclusions developed for these models. Ringler et al. [4] showed that their C-grid model conserves kinetic energy for shallow water equations by properly interpolating the neighbors’ normal fluxes to obtain a tangential flux at an edge center. However, if these normal fluxes are not on the same plane as the interpolated tangential flux, their analyses might require further improvement.

The Z-grid model is unique in predicting all the scalar states of vorticity and divergence. The model can easily treat momentum on a sphere as scalars, provide the best dispersion relation among the finite volume schemes [7], and is free of computational modes [8]. It resolves the stationary geostrophic mode perfectly; this is difficult to achieve for all other grids, except the C-grid scheme, which can acquire this property by carefully constructing tangential velocity at an edge center from the normal velocities of neighbor edges [3]. Thuburn et al. [3] also assumed that these normal velocities are on the same plane as the tangential velocities. However, despite its elegance and perfection, a Z-grid model is relatively expensive to solve compared to other models in terms of the velocity vector. Its continuous model equations include the vorticity, divergence, stream function, and velocity potential. This implies that Poisson equations must be solved at each time-integration step to convert vorticity and divergence into the stream function and velocity potential, respectively. Given that a multigrid technique [9], which is vital in the implementation of the Z-grid model, is an iterative method, it can utilize the solutions of the stream function and velocity potential from the previous time step and obtain stream functions and velocity potentials in a few iterations. Further improvement of the Z-grid efficiency may be critical to its future implementation for operational applications. In this study, we investigated the possibility of modifying the Z-grid scheme to improve its stability such that a longer time step could be used in model integration.

Heikes and Randall [6] and Eldred and Randall [10] proposed a Z-grid model and placed all state variables at the Voronoi center of its cell. The normal derivative of a state variable at an edge’s center is approximated according to the finite difference between two state values at the two adjacent Voronoi centers of this edge, that is, a two-grid stencil is used. For tangential derivatives, the state values at each end of an edge are first calculated, and then a finite difference is applied between the interpolated values at these ends. Each end for a given edge of an icosahedral grid has three neighbors; hence, the tangential derivative calculation involves four stencils. Heikes and Randall [6] demonstrated their numerical results and showed the stability and conservation properties for numerous testing cases. For simplicity, we refer to this Z-grid scheme as the Voronoi Z-grid. When it is applied on a rectangular grid, six grid cells are used for calculating the tangential derivatives. To summarize, the Voronoi Z-grid model uses a two-grid stencil for normal derivatives and stencils with four or more grids for tangential derivatives.

To improve the efficiency and possible accuracy on spheres, we propose the generalization of the Voronoi Z-grid model by using state variables at the centroidal centers and the same grid stencil for calculating both the normal and tangential derivatives. The motivation of this generalization in this study was to seek a more stable scheme so that it allowed us to apply longer time steps than those applied in the Voronoi Z-grid. If such a scheme could be achieved, it would be more efficient and save computing time. In this paper, we refer to this generalization as the centroidal Z-grid and present this generalization scheme as a regional model for shallow water equations to demonstrate the improvement. While this generalization can possess the same properties as the Voronoi Z-grid, it is important to study how the dispersion relation is affected as it reflects how closely the discretized schemes represent the continuous equations. We analyzed the dispersion relation of this centroidal Z-grid by comparing it to the Voronoi Z-grid and developed some principles for balancing the stability and dispersion relation in the construction of the centroidal Z-grid.

To demonstrate the stability and efficiency of this scheme, we chose two typical test cases that possess nonlinear instability: the Rossby–Haurwitz case [11] and the unstable jet over an isolated terrain, the Galewsky case [12]. A rectangular domain is sufficient to evaluate the efficiency and we used the well-known Poisson equation solver to implement Z-grid models to simplify our discussions in this paper. Note we have modified these test cases from their spherical shallow water equations to a rectangular domain for the numerical experiments performed in this study. The Rossby–Haurwitz case is not an exact solution for shallow water equations. However, a perfect solution satisfying the continuous equations is more appropriate for comparison and analyzing the accuracy and stability. In this study, we modified the Rossby–Haurwitz test case and derived the forcing terms of these equations so that the it became an exact solution for the continuous and forced shallow water equations. This methodology is also recommended by Browning et al. [13]. For the Galewsky test case, we just modified the height derivation for a rectangular domain simulating a global latitudinal/longitudinal grid.

The remainder of this paper is organized as follows. In Section 2, we introduce the new Z-grid for shallow water equations and show its second-order accuracy. The analysis of the dispersion relation of the centroidal Z-grid is presented in Section 3. In Section 4, we discuss an improvement of the Rossby–Haurwitz test case for a more rigorous test of our model’s stability with respect to shallow water equations on a latitudinal/longitudinal plane. We also describe the height field derivation for the Galewsky test case. Section 5 presents the numerical experiments for the centroidal Z-grid and Voronoi Z-grid models. Finally, we provide our conclusions and remarks in Section 6.

2. Generalization of the Z-Grid Model

Shallow water equations are used in this paper to illustrate the generalization of the Z-grid model, analyze dispersion relation, and demonstrate numerical results. A Z-grid model solves the following system of shallow water equations:

$$\frac{\partial \eta }{\partial t}+\nabla ·\left(\eta \nabla \chi \right)-J\left(\eta ,\psi \right)=0,$$

(1)

$$\frac{\partial \delta }{\partial t}-\nabla ·\left(\eta \nabla \psi \right)-J\left(\eta ,\chi \right)+{\nabla }^2\left[K+g\left(h+h_s\right)\right]=0,$$

(2)

$$\frac{\partial h}{\partial t}+\nabla ·\left(h\nabla \chi \right)-J\left(h,\psi \right)=0,$$

(3)

where $\eta$ is the absolute vorticity, $\delta$ is the divergence, $\psi$ is the stream function, and $\chi$ is the velocity potential, which satisfies the following equations:

$${\nabla }^2\psi =\eta -f\mathrm{and}{\nabla }^2\chi =\delta ,$$

(4)

where $f$ is the Coriolis parameter, and the kinetic energy is

$$K=\frac{1}{2}\left[\nabla ·\left(\psi \nabla \psi \right)-\psi {\nabla }^2\psi +\nabla ·\left(\chi \nabla \chi \right)-\chi {\nabla }^2\chi \right]+J\left(\psi ,\chi \right).$$

(5)

In Equation (2), $g$ is the gravitational constant, and $h$ and $h_s$ are the fluid depth and height of the underlying surface topography, respectively. Note that the energy- and enstrophy-conservation Z-grid models solve a slightly different form of the abovementioned equations, as given by Eldred and Randall [10]. To evaluate the stability, we still used the abovementioned equations in this study.

Z-grid-based models have three important properties. First, they deal with all-scalar variables, which are relatively easy to solve (particularly on a sphere), thereby avoiding the complexity accrued when considering velocity vectors. Second, they separate vorticity from divergence. The fast mode waves are in the divergence equation given by Equation (2). For divergence-free modes, the vorticity evolves on its own without any high-frequency waves. In the linearized system, the geostrophic mode is perfectly stationary. Third, because all state variables are placed at the centers of cells, the number of state variables match the number of equations. Such a system does not have computational modes [8], with the only drawback being that Equations (4) must be solved at every time step because this system requires stream function $\psi$ and velocity potential $\chi$.

A multigrid technique [9] is commonly used to efficiently solve Equations (4) and is an iterative method. That is, the closer the initial value is to the true solution, the fewer are the required iterations. In the temporal integration of a Z-grid model, the $\psi$ and $\chi$ in the next time step are very close to their counterparts in the current time step. The multigrid iterative process can fully take advantage of this fact by using the $\psi$ and $\chi$ in the current time as the initial guesses for the solutions of $\psi$ and $\chi$ in the next time step. In this study, instead of determining ways to further improve the multigrid solver, we investigated whether we could construct a more stable Z-grid model so that a longer time step could be used to improve the efficiency.

A finite-volume Z-grid model constructs its scheme by integrating Equations (1)–(3) over each grid cell $c$ and using the following integral relations [6] of Jacobian, divergence, and Laplacian operators:

$${\int }_c^{}J\left(\sigma ,\beta \right)da={\oint }_{\partial c}^{}\alpha \frac{\partial \beta }{\partial \tau }dl,$$

(6)

$${\int }_c^{}\nabla ·\left(\sigma \nabla \beta \right)da={\oint }_{\partial c}^{}\alpha \frac{\partial \beta }{\partial n}dl,$$

(7)

$${\int }_C^{}{\nabla }^2\alpha da={\oint }_{\partial c}^{}\frac{\partial \alpha }{\partial n}dl,$$

(8)

where $\alpha$ and $\beta$ are arbitrary scalar functions, $c$ is any given horizontal grid cell, $\partial c$ is the perimeter of $c$, and $n$ and $\tau$ are the unit outward normal and tangential vectors, respectively, following the right-hand rule with the thumb pointing upward. In other words, Equations (1)–(3) can be solved by applying these relations repeatedly. For any given grid cell $c$, regardless of its shape, a Voronoi Z-grid scheme places all state variables at Voronoi center $c$, as illustrated in Figure 1.

By considering the Laplacian operator of Equation (8) as an example, we have the following approximation:

$${\oint }_{\partial c}^{}\frac{\partial \alpha }{\partial n}dl\cong {\displaystyle \sum _i}\frac{{\alpha }_2-{\alpha }_1}{d}{l}_i.$$

(9)

The continuous curve integral along $\partial c$ is approximated by a set of piece-wise line integrals along $\left\{l_i\right\}$ surrounding cell $c$, where $i$ runs through all edges of $\partial c$. These line integrals in (9) are then approximated by the directional derivatives multiplied by the edge lengths. For normal derivative $\frac{\partial \alpha }{\partial n}$, a two-stencil scheme is applied: $\frac{{\alpha }_2-{\alpha }_1}{d}$. If all of the cells are regular (i.e., the Voronoi and centroid centers are at the same location), then this approximation can achieve second-order accuracy. To simplify the discussion, the Z-grid models are discussed for planes and the cells are assumed to be regular.

A generalization of a Voronoi Z-grid model is achieved by using a grid stencil beyond the two Voronoi centers for normal derivatives of edges in (9). On a grid of a two-dimensional plane, the second-order accurate approximation of $\frac{\partial \alpha }{\partial n}$ can be ensured using a six-grid or larger stencil based on Taylor expansions, regardless of whether a Voronoi or centroid center is used. As a finite volume method defines the state variables at centroidal centers, we always choose a centroidal center in this generalized scheme regardless of whether it is used on a plane or sphere where Voronoi and centroid centers may differ. Thus, we call this scheme the centroidal Z-grid hereafter. To ensure that the generalized scheme conserves scalar quantities, the line integrals are constructed using the same stencils for two cells sharing common edges. Hence, the six-grid stencil must be symmetric on both sides of the edge. Hexagons of an icosahedral grid cannot be symmetric; thus, we divided these hexagons and pentagons into triangles. Figure 2 shows the symmetry for both triangular and rectangular grids. A convenient feature of this triangular grid on a sphere is that each cell has three edges, with each edge being associated with the six closest grid cells. Moreover, all the three edges of a cell are in the same plane. This homogeneity of the triangular grid makes parallel computing much easier, and in addition, the grid can be easily nested compared to a hexagonal/pentagonal grid.

The Voronoi Z-grid model interpolates function values at both the ends of an edge and uses a centered finite difference to approximate the tangential derivative at the edge center. For tangential derivatives, the new proposed Z-grid model applies the same stencil as the Voronoi Z-grid model for achieving second-order accuracy. For the sake of simplicity, in the current study, we used a rectangular grid for the model formulation and numerical experiments; we will discuss an icosahedral grid on a sphere in more detail in our future works.

Note that the subindex integers of 1,2,…,6 in Figure 2 serve a symbolic purpose only. For each edge, the indices are different. For a rectangular grid, the Voronoi Z-grid model approximates the line integrals of the operators of Equations (6)–(8) as follows:

$${\oint }_{\partial c}^{}\alpha \frac{\partial \beta }{\partial \tau }dl\cong {\displaystyle \sum _l}\frac{\left({\alpha }_2+{\alpha }_5\right)}{2}\frac{{\beta }_3+{\beta }_4-{\beta }_1-{\beta }_6}{4},$$

(10)

$${\oint }_{\partial c}^{}\alpha \frac{\partial \beta }{\partial n}dl\cong {\displaystyle \sum _l}\frac{\left({\alpha }_2+{\alpha }_5\right)}{2}\frac{{\beta }_5-{\beta }_2}{d}l,$$

(11)

$${\oint }_{\partial c}^{}\frac{\partial \alpha }{\partial n}dl\cong {\displaystyle \sum _l}\frac{{\alpha }_5-{\alpha }_2}{d}l.$$

(12)

For the tangential derivative, the model actually uses

$$\left(\frac{{\beta }_3+{\beta }_4+{\beta }_2+{\beta }_5}{4}-\frac{{\beta }_2+{\beta }_5+{\beta }_1+{\beta }_6}{4}\right)/l$$

(13)

and the grid values of ${\beta }_2$ and ${\beta }_5$ in (13) cancel out in Equation (10) for the rectangular grid. Thus, we showed that the model definitely uses a six-grid stencil. The proposed centroidal Z-grid uses the same stencil as that of the Voronoi Z-grid. However, the Voronoi Z-grid model uses different stencils for approximating tangential and normal derivatives. Our new generalization retains the tangential approximation and uses a six-grid stencil to approximate all the normal derivatives in the following approximations:

$${\oint }_{\partial c}^{}\alpha \frac{\partial \beta }{\partial n}dl\cong {\displaystyle \sum _l}\frac{\left({\alpha }_2+{\alpha }_5\right)}{2}\left(w_t\frac{{\beta }_4-{\beta }_3}{d}+w_m\frac{{\beta }_5-{\beta }_2}{d}+w_b\frac{{\beta }_6-{\beta }_1}{d}\right)l,$$

(14)

$${\oint }_{\partial c}^{}\frac{\partial \alpha }{\partial n}dl\cong {\displaystyle \sum _l}\left(w_t\frac{{\alpha }_4-{\alpha }_3}{d}+w_m\frac{{\alpha }_5-{\alpha }_2}{d}+w_b\frac{{\alpha }_6-{\alpha }_1}{d}\right)l,$$

(15)

where $w_t$, $w_m,$ and $w_b$ are weights satisfying $w_t\ge 0$, $w_m\ge 0$, and $w_b\ge 0$, respectively, and

$$w_t+w_m+w_b=1$$

(16)

Parameters $t$, $m$, and $b$ denote the top, middle, and bottom cells on the right column of Figure 2. When $w_t=w_b=0$ and $w_m=1$, the corresponding model reduces to the Voronoi Z-grid model shown previously. Thus, this proposed model with choices of parameters $w_t$, $w_m,$ and $w_b$ is a generalization of the Z-grid model. By applying Taylor expansion, we can prove that this new model is second-order accurate on a rectangular grid, possesses all conservation properties of the Voronoi Z-grid model, and provides flexibility in testing other choices of the weights. Moreover, we demonstrate that some selection of nonzero $w_t$ and $w_b$ provides more stability and can utilize longer time steps in the temporal integration of the discretized shallow water equations. Our numerical experiments showed that if $w_t=w_b<w_m$ satisfying (16) then the Z-grid model is stable, but a scheme of $w_t=w_b=w_m=\frac{1}{3}$ does not show stability.

3. Dispersion Relation Analysis

Before demonstrating the stability of the centroidal Z-grid, determining how this generalization affects the dispersion relation would be interesting, as it is an important feature approximating the continuous equations. We analyzed the centroidal Z-grid dispersion relation and compared it to the Voronoi Z-grid in order to have a better understanding of how the dispersion relation changes. Furthermore, the relation may provide some guidance to choosing weights $w_t$, $w_m$, and $w_b$.

In this study, a linearized rotating shallow-water-equation system was used in the dispersion relation analysis:

$$\frac{\partial u}{\partial t}-fv+g\frac{\partial h}{\partial x}=0,$$

(17)

$$\frac{\partial v}{\partial t}+fu+g\frac{\partial h}{\partial y}=0,$$

(18)

$$\frac{\partial h}{\partial t}+H\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)=0.$$

(19)

To convert these equations into the Z-grid equations, we need to derive the vorticity and divergence equations. First, we applied $\nabla ·$ in Equations (17) and (18) and obtained

$$\frac{\partial \delta }{\partial t}-f\zeta +g{\nabla }^{2}h=0,$$

(20)

where $\zeta =-u_y+v_x$ and $\delta =u_x+v_y$. Then, by applying $\nabla \times$ on these Equations, we obtained

$$\frac{\partial \zeta }{\partial t}+f\delta =0.$$

(21)

We applied $\frac{\partial }{\partial t}$ in Equation (20) to obtain

$$\frac{{\partial }^2\delta }{\partial t^2}-f\frac{\partial \zeta }{\partial t}+g\frac{\partial {\nabla }^{2}h}{\partial t}=0,$$

(22)

and by using (21), we have

$$\frac{{\partial }^2\delta }{\partial t^2}+f^2\delta +g\frac{\partial {\nabla }^{2}h}{\partial t}=0.$$

(23)

Then, Equation (19) can be rewritten as

$$\frac{\partial h}{\partial t}=-H\delta ,$$

(24)

and used to replace the divergence in Equation (23) so that

$$-\frac{{\partial }^{3}h}{\partial t^3}-f^2\frac{\partial h}{\partial t}+gH\frac{\partial {\nabla }^{2}h}{\partial t}=0.$$

(25)

Assume a wave solution of $h$ to have the form $h=h_0{e}^{i\left(kx+ly-\omega t\right)}$. By substituting this solution into Equation (25), we have

$$-i{\omega }^3+i{f}^2\omega +igH\omega \left(k^2+l^2\right)=0$$

(26)

or

$$\omega \left[{\omega }^2-f^2-gH\left(k^2+l^2\right)\right]=0.$$

Let $\lambda =\frac{\sqrt{gH}}{f}$ be the deformation radius of this shallow water Equation; then the dispersion relation satisfies

$$\omega \left[{\left(\frac{\omega }{f}\right)}^2-1-{\lambda }^2\left(k^2+l^2\right)\right]=0.$$

(27)

The three eigenvalues are represented as

$${\omega }_1=0$$

and

$${\left(\frac{{\omega }_{2,3}}{f}\right)}^2=1+{\lambda }^2\left(k^2+l^2\right)=1+{\left(\frac{\lambda }{d}\right)}^2\left[{\left(kd\right)}^2+{\left(ld\right)}^2\right],$$

(28)

which is the continuous dispersion relation and where $d$ is any constant. Any discrete numerical scheme should approach this relation as much as possible. The closer a scheme is to the continuous relation of Equation (28), the better is the scheme’s dispersion relation.

Now, we look into the discrete forms (10), (14), and (15) of a centroidal Z-grid model in solving Equation (20). Equation (20) mainly involves calculating a Laplacian operator. For a uniform grid with a grid spacing of $d$ in both $x$ and $y$ directions, a Z-grid model places the model states at the center of a square cell. The Z-grid model should approximate ${\nabla }^{2}h$ over a cell by using the Laplacian relation

$$\underset{c}{\overset{}{{\displaystyle \iint }}}{\nabla }^{2}hda={\displaystyle \underset{\partial c}{\overset{}{\oint }}}\frac{\partial h}{\partial n}dl={\displaystyle \sum _{l=1}^4\underset{l}{\overset{}{\int }}}\frac{\partial h}{\partial n}dl,$$

(29)

where $n$ is the outward normal direction. For a cell located at $\left(m,j\right)$, the centroidal Z-grid model uses six stencils to approximate the normal derivative. For the left edge, we have

$${\frac{\partial h}{\partial n}|}_{m-\frac{1}{2},j}\cong w_t\frac{h_{m-1,j+1}-h_{m,j+1}}{d}+w_m\frac{h_{m-1,j}-h_{m,j}}{d}+w_b\frac{h_{m-1,j-1}-h_{m,j-1}}{d},$$

(30)

where $w_t,w_m$ and $w_b$ satisfying (16). By assuming $w_t=w_b$ for the right, bottom, and top edges, we also have

$${\frac{\partial h}{\partial n}|}_{m+\frac{1}{2},j}\cong w_t\frac{h_{m+1,j+1}-h_{m,j+1}}{d}+w_m\frac{h_{m+1,j}-h_{m,j}}{d}+w_t\frac{h_{m+1,j-1}-h_{m,j-1}}{d}$$

$${\frac{\partial h}{\partial n}|}_{m,j-\frac{1}{2}}\cong w_t\frac{h_{m-1,j-1}-h_{m-1,j}}{d}+w_m\frac{h_{m,j-1}-h_{m,j}}{d}+w_t\frac{h_{m+1,j-1}-h_{m+1,j}}{d}$$

$${\frac{\partial h}{\partial n}|}_{m,j+\frac{1}{2}}\cong w_t\frac{h_{m-1,j+1}-h_{m-1,j}}{d}+w_m\frac{h_{m,j+1}-h_{m,j}}{d}+w_t\frac{h_{m+1,j+1}-h_{m+1,j}}{d}.$$

Let $h_{m,j}=h_0{e}^{i\left(kmd+ljd-\omega t\right)}$, then

$$h_{m+1,j}-h_{m,j}=h_0{e}^{i\left[k\left(m+\frac{1}{2}\right)d+ljd-\omega t\right]}\left(e^{i\frac{kd}{2}}-e^{-i\frac{kd}{2}}\right)$$

$$=h_0{e}^{i\left[k\left(m+\frac{1}{2}\right)d+ljd-\omega t\right]}2i\sin \left(\frac{kd}{2}\right).$$

Thus, Equation of (30) and the three Equations next it become the following forms,

$${\frac{\partial h}{\partial n}|}_{m-\frac{1}{2},j}\cong -2i{h}_0{e}^{i\left[k\left(m-\frac{1}{2}\right)d+ljd-\omega t\right]}\left[2w_t\cos \left(ld\right)+w_m\right]\sin \left(\frac{kd}{2}\right)/d$$

$${\frac{\partial h}{\partial n}|}_{m,j-\frac{1}{2}}\cong -2i{h}_0{e}^{i\left[kmd+l\left(j-\frac{1}{2}\right)d-\omega t\right]}\left[2w_t\cos \left(kd\right)+w_m\right]\sin \left(\frac{ld}{2}\right)/d$$

$${\frac{\partial h}{\partial n}|}_{m+\frac{1}{2},j}\cong 2i{h}_0{e}^{i\left[k\left(m+\frac{1}{2}\right)d+ljd-\omega t\right]}\left[2w_t\cos \left(ld\right)+w_m\right]\sin \left(\frac{kd}{2}\right)/d$$

$${\frac{\partial h}{\partial n}|}_{m,j+\frac{1}{2}}\cong 2i{h}_0{e}^{i\left[kmd+l\left(j+\frac{1}{2}\right)d-\omega t\right]}\left[2w_t\cos \left(kd\right)+w_m\right]\sin \left(\frac{ld}{2}\right)/d.$$

Therefore, the finite volume scheme gives

$$\frac{1}{d^2}{\displaystyle \underset{c}{\overset{}{\iint }}}{\nabla }^{2}hda=\frac{1}{d^2}\left({\frac{\partial h}{\partial n}|}_{m-\frac{1}{2},j}+{\frac{\partial h}{\partial n}|}_{m+\frac{1}{2},j}+{\frac{\partial h}{\partial n}|}_{m,j-\frac{1}{2}}+{\frac{\partial h}{\partial n}|}_{m,j+\frac{1}{2}}\right)d$$

$$=2i{h}_0{e}^{i\left[kmd+ljd-\omega t\right]}\left[2w_t\cos \left(ld\right)+w_m\right]{\sin }^2\left(\frac{kd}{2}\right)/d^2+$$

$$2i{h}_0{e}^{i\left[kmd+ljd-\omega t\right]}\left[2w_t\cos \left(kd\right)+w_m\right]{\sin }^2\left(\frac{ld}{2}\right)/d^2.$$

By assuming ${\delta }_{m,j}={\delta }_0{e}^{i\left(kmd+ljd-\omega t\right)}$, ${\zeta }_{m,j}={\zeta }_0{e}^{i\left(kmd+ljd-\omega t\right)}$ and replacing ${\nabla }^{2}h$ in Equation (20) with this discrete form, Equation (20) is rewritten as

$$-\omega {\delta }_0-f{\zeta }_0+4\frac{g{h}_0\left[w_t\cos \left(kd\right)+w_t\cos \left(ld\right)+w_m\right]\left[{\sin }^2\left(\frac{kd}{2}\right)+{\sin }^2\left(\frac{ld}{2}\right)\right]}{d^2}=0.$$

Further, Equation (21) gives

$$-\omega {\zeta }_0+f{\delta }_0=0,$$

and Equation (19) gives

$$-\omega h_0=-H{\delta }_0.$$

By replacing ${\delta }_0$ and ${\zeta }_0$ in the abovementioned spectral Equation, we have

$$-{\omega }^2\frac{h_0}{H}+f^2\frac{h_0}{H}+4\frac{g{h}_0\left[w_t\cos \left(kd\right)+w_t\cos \left(ld\right)+w_m\right]\left[{\sin }^2\left(\frac{kd}{2}\right)+{\sin }^2\left(\frac{ld}{2}\right)\right]}{d^2}=0$$

or

$${\left(\frac{\omega }{f}\right)}^2=1+4{\left(\frac{\lambda }{d}\right)}^2\left[w_t\cos \left(kd\right)+w_t\cos \left(ld\right)+w_m\right]\left[{\sin }^2\left(\frac{kd}{2}\right)+{\sin }^2\left(\frac{ld}{2}\right)\right].$$

(31)

For the given relation $w_m=1-2w_t$ and assuming that $\frac{\lambda }{d}=2$, we can find optimal weights that minimize the difference between Equations (28) and (31) in $l_2$ norm. Over the region $\left(kd,ld\right)\in \left(0,\pi \right)\times \left(0,\pi \right)$, the optimal weights are close to $w_m=2$ and $w_t=-0.5$, which possess an even better dispersion relation than the Voronoi Z-grid. However, we found that these optimal weights could result in an unstable model. Thus, we determined to seek weights that provided stability and the best possible dispersion relation.

To determine the difference between the continuous and discrete dispersion relations, we calculated the difference between the right hand sides of Equations (28) and (31) as functions of $\left(kd,ld\right)$ over three domains from the origin through 1) $\sqrt{{\left(kd\right)}^2+{\left(ld\right)}^2}\le \frac{\pi }{2}$; 2) $\sqrt{{\left(kd\right)}^2+{\left(ld\right)}^2}\le \pi$; and 3) $0\le kd\le \pi$ and $0\le ld\le \pi$. We compared three typical centroidal Z-grids with weights selected as ($w_m=0.5$, $w_t=0.25$), ($w_m=0.75$, $w_t=0.125$), and ($w_m=1.02$, $w_t=-0.01$) to the Voronoi Z-grid. The differences between Equations (28) and (31) over these three domains are shown in

Table 1. Dispersion relation difference between discrete and continuous dispersion relations.

Radius	Voronoi Z-Grid	Centroidal Z-Grid (0.5, 0.25)	Centroidal Z-Grid (0.75, 0.125)	Centroidal Z-Grid (1.02, −0.01)
$\pi /2$	0.113	0.307	0.207	0.106
$\pi$	0.887	1.957	1.358	0.854
$\left(\pi ,\pi \right)$	1.199	2.927	1.910	1.150

.

This table shows that there is a better choice of weights than those possessed by the Voronoi Z-grid as they give a better dispersion relation ($w_m=1.02$, $w_t=-0.01$). We also found in our numerical experiments that, with this choice of weights, the centroidal Z-grid is as stable as the Voronoi Z-grid. Thus, for the centroidal Z-grid models, we can choose weights such that their dispersion relation is as good as possible and provides more stability.

4. Modification of the Test Cases

To show the stability of the centroidal Z-grid, we chose two typical test cases for shallow water equations: the Rossby–Haurwitz [11] and Galewsky [12] test cases. We modified these cases to fit our numerical experiments over a rectangular grid simulating a global latitudinal/longitudinal grid.

The Rossby–Haurwitz test case is used for testing various models, including spectral and gridded models, because of its realism and nonlinear instability that closely reflect the real atmosphere. The stream function in this case provides a realistic wave pattern around a globe, and we extended it in this study to make the stream function move with time $t$ as follows:

$$\psi =-a^2\omega \sin \theta +a^2\kappa {\cos }^R\theta \sin \theta \cos R\left(\lambda +\upsilon t\right),$$

(32)

where $a$ is the radius of the Earth; $\theta \in \left(-\frac{\pi }{2},\frac{\pi }{2}\right)$ is the latitude; $\lambda \in \left(0,2\pi \right)$ is the longitude; $\omega$, $\kappa$ are constants and $R$ is an integer constant; and $t$ is the time in seconds. This stream function moves from west to east at an angular velocity given by the following Equation:

$$\upsilon =\frac{R\left(3+R\right)\omega -2\Omega }{\left(1+R\right)\left(2+R\right)},$$

(33)

where $\Omega$ is the rotational angular velocity of the Earth. This is a divergence-free test case, that is, $\delta =0$ and $\chi =0$. We adapted it to a latitudinal/longitudinal plane for testing the numerical accuracy and stability of the centroidal Z-grid scheme. However, this test case (Williamson et al. [11]) does not satisfy the continuous shallow water Equations (1)–(5). A common practice is to compare a given model output to the numerical output produced through a spectral model that solves the dissipative shallow water equations by using an appropriate initial condition (e.g., [4]). A ${\nabla }^4$ dissipation was used in a spectral model to control noise from high-frequency waves. Such a comparison may not serve our purpose, as we need to know the stability of our Z-grid models. A diffusive solution may not be able to demonstrate the differences of these Z-grid models. Hence, we attempted to search for a non-diffusive solution to check the stability of our model for the Rossby–Haurwitz case.

We selected a method suggested by Browning et al. [13]. This method adds forcing terms to a given system of equations to force a Rossby–Haurwitz case to become the exact solution for the non-diffusive shallow water Equations (1)–(5). For our test on the latitudinal/longitudinal plane, we used the Rossby–Haurwitz stream function in Equation (32) and chose the fluid depth $h$ such that the divergence Equation (2) was satisfied perfectly. In other words, by substituting $h_s=0$, $\delta =0$, and $\chi =0$ in Equations (2) and (5), we solved $h$, satisfying the following:

$$g{\nabla }^{2}h=\nabla ·\left(\eta \nabla \psi \right)-\frac{1}{2}{\nabla }^2\left(\nabla ·\left(\psi \nabla \psi \right)-\psi {\nabla }^2\psi \right).$$

(34)

By choosing $h$ in such a manner, Equation (2) is satisfied, and no forcing term is required for Equation (2). We then applied the stream function given by Equation (32), its corresponding vorticity $\eta$, and fluid depth $h$ that satisfies Equation (34) into Equations (1) and (3) and obtained two forcing functions, denoted as $f_{\eta }$ and $f_h$. Then, the functions $\left(\psi ,\chi ,h\right)$, $\chi =0$, and solution $h$ of Equation (34) satisfy the following non-diffusive Equations:

$$\frac{\partial \eta }{\partial t}+\nabla ·\left(\eta \nabla \chi \right)-J\left(\eta ,\psi \right)=f_{\eta },$$

(35)

$$\frac{\partial \delta }{\partial t}-\nabla ·\left(\eta \nabla \psi \right)-J\left(\eta ,\chi \right)+{\nabla }^2\left[K+g\left(h+h_s\right)\right]=0,$$

(36)

$$\frac{\partial h}{\partial t}+\nabla ·\left(h\nabla \chi \right)-J\left(h,\psi \right)=f_h.$$

(37)

Our numerical forecast results can be verified against this true solution at any given time $t$ to demonstrate the stability and efficiency of the centroidal and Voronoi Z-grid models.

This modification methodology for shallow-water-equation test cases can also be applied to the Rossby–Haurwitz test case on a sphere. The only difference is that the differential operators will be spherical and not Cartesian in the latitudinal and longitudinal directions. A similar process can also be applied to other test cases of shallow water equations. However, the Rossby–Haurwitz case is a very challenging case because of its nonlinearity. Hence, the comparison of the stability of both the models is sufficient.

Note that Equation (34) may be too complicated to solve for an analytic solution of $h$. Instead, for any given resolution and time $t$ in our numerical experiments, we solved Equation (34) numerically.

To test a regional Z-grid model, we must specify the boundary conditions. In general, a global or larger domain model provides Dirichlet boundary conditions to a regional Z-grid model. For the Rossby–Haurwitz test case, we chose a latitudinal and longitudinal domain $\left(-\frac{\pi }{2},\frac{\pi }{2}\right)\times \left(0,2\pi \right)$ and considered the east–west boundaries as periodic boundaries, thereby simulating a global latitudinal and longitudinal grid and using the Dirichlet boundary conditions at the north and south boundaries from the continuous solution of this test case.

Numerical noise usually appears near the Dirichlet boundaries because the boundary values do not match the predicted values of the model due to truncation errors, even though these errors are very small. A sponge boundary treatment [14] is usually used to minimize the spurious wave reflections. In our implementation of the Z-grid models, we used a simple sponge method to avoid the spurious reflections at the north and south boundaries. We used the true boundary values and forecast grid values at the second interior to interpolate the grid values at the first interior next to the boundaries (hereafter referred to as the averaging boundary conditions). This method removes boundary reflections without loss of model accuracy. Therefore, we applied this averaging boundary condition throughout our numerical experiments.

The Galewsky test case was also modified to fit our numerical experiments over a rectangular domain. We kept the zonal flow formulation obtained from [12],

$$u\left(\varphi \right)=\{\begin{array}{cc}0& \varphi \le {\varphi }_0\\ \frac{u_{max}}{e_n}exp\left[\frac{1}{\left(\varphi -{\varphi }_0\right)\left(\varphi -{\varphi }_1\right)}\right]& {\varphi }_0\le \varphi \le {\varphi }_1\\ 0& \varphi \ge {\varphi }_1\end{array}$$

(38)

and the perturbation field,

$$h^{\prime }\left(\lambda ,\varphi \right)=\widehat{h}\cos \varphi e^{-{\left(\lambda /\alpha \right)}^2}{e}^{-{\left[\left({\varphi }_2-\varphi \right)/\beta \right]}^2},$$

(39)

where $u_{max}=80$ ms⁻¹, ${\varphi }_0=\pi /7$, ${\varphi }_1=\frac{\pi }{2}-{\varphi }_0$, $e_n=exp\left[-4/{\left({\varphi }_1-{\varphi }_0\right)}^2\right]$, ${\varphi }_2=\pi /4$, $\alpha =1/3$, $\beta =1/15$, and $\widehat{h}=120$ m. Instead of using the height field in [12], we solved Equation (34) as well to obtain a balanced height field for the background flow given by Galewsky et al. [12]. We modified the derivative operators in (34) as a rectangular grid coordinate of $\left(a\lambda ,a\varphi \right)$ so that the Equation for the height field has the following form:

$$g{\nabla }^{2}h=-\frac{2\Omega \left(u\cos \varphi +u_{\varphi }\sin \varphi \right)}{a}.$$

(40)

A mean constant height, $h_0={10}^4$ m, was added to the solution of (40) so that the mean average height was around 10 km, as suggested by Galewsky et al. [12]. All the partial derivatives were treated as a Cartesian coordinate of $\left(\lambda ,\varphi \right)$ instead of a spherical coordinate.

5. Numerical Experiments

As stated earlier, a multigrid technique is usually applied in Z-grid models for solving the Poisson Equations in (4), converting vorticity and divergence to a stream function and velocity potential, respectively. For an accurate comparison between the centroidal and Voronoi Z-grid models, we used FISHPACK [15] in our regional Z-grid model experiments. This is not an iterative method but solves Poisson equations numerically with second-order accuracy. However, this solver is not sufficiently efficient to implement a Z-grid model in general. By using this solver, we avoid clarifying the implementation of a multigrid method and its convergence solving the Poisson equations numerically. Given that our goal is to show the improved stability and longer time step integration of the new Z-grid model, FISHPACK is appropriate.

The FISHPACK solver enables us to specify the Dirichlet boundary conditions at the north and south boundaries, and periodic boundary conditions at the east and west boundaries. Our domain involves the latitudinal and longitudinal grids of the whole globe, and these boundary conditions meet our needs. For all numerical experiments presented herein, the same solver was used for all Z-grid-model comparisons.

We implemented three Z-grid models to demonstrate the improvement of the new Z-grid model. For the Voronoi Z-grid model, we used the finite volume scheme defined by Equations (10), (14), and (15) by setting $w_t=w_b=0$ and $w_m=1$, denoted as VORO (Voronoi scheme) hereafter. For the new centroidal Z-grid model, we used $w_t=w_b=1/8$ and $w_m=3/4$; this scheme is denoted as CENT. We also chose the centroidal Z-grid with the “best” dispersion relation, with $w_t=w_b=-0.01$ and $w_m=1.02,$ and denoted it as BEST.

In our experiments, we chose parameters $\kappa =\omega =7.848\times {10}^{-6}\mathrm{and}R=4$ for the Rossby–Haurwitz test case. Other physical parameters used in the experiments of both cases were as follows: $g=9.8$, $a=6.37122\times {10}^6$, and $\Omega =7.2722\times {10}^{-5}$.

Figure 3 shows the stream function given by Equation (32), along with its absolute vorticity and fluid depth of the continuous Rossby–Haurwitz test case on the latitudinal/longitudinal plane in a one-degree resolution. For the Galewsky case, the background flow is given in Equations (38) and the solution of (40) and the perturbation (39) is shown in [12].

For the time integration of our numerical experiments, we used a Runge–Kutta fourth-order scheme, even though it is relatively expensive compared to the third-order Adams–Bashforth scheme [16]. Both Runge–Kutta and Adams–Bashforth schemes are classified as Eulerian-based time-integration (EBIT) schemes [17]. They are sufficient to support the numerical experiments in this study. The Runge–Kutta fourth-order scheme is free of a computational mode and is less diffusive. It serves the purpose of comparing the model stability well. An Adams–Bashforth scheme or some path-based time-integration (PBTI) schemes [17] will be implemented in practical applications.

All numerical experiments in this Rossby–Haurwitz case are free of any numerical dissipation. Even though Galewsky et al. [12] suggested a second-order dissipative for all equations, we applied a second-order dissipation to the vorticity and divergence Equations (1) and (2) only and no numerical diffusion for (3) for the Galewsky case. The dissipation in vorticity equation is sufficiently small to remove unresolvable noise. For each model, we tried several time step lengths until the longest one was found for stability. For the Rossby–Haurwitz case, the numerical forecast ran for 14 days and for the Galewsky case, it ran for 6 days, the same as in [12]. The time step used by the CENT scheme can be as much as 1.4 times that of the VORO and BEST schemes for the Rossby–Haurwitz case and 1.6 times for the Galewsky case. For example, at a one-degree resolution, the CENT scheme can take a time step of 840 s but VORO took just 600 s for the Rossby–Haurwitz case; CENT can take one of 480 s but VORO took just 300 s for the Galewsky case. These longer time-integration steps of CENT are indeed significant. The computer runtime confirms the efficiency improvement. For the Rossby–Haurwitz case, the runtime is shown in Figure 4 with three different resolutions. Computing time is saved significantly as the resolutions increase.

As we know the exact solution of the forced shallow-water Equations for the Rossby–Haurwitz case, we calculated the errors against the true solution of $\left(\psi ,\chi ,h\right)$ at the end of the 14th day in two error norms: $l_2$ and $l_{\infty }$. In the error table (

Table 2. Forecast errors on the 14th day of the two Z-grid models of CENT, VORO and BEST.

Errors	$\mathit{\eta }$ (${\mathit{E}}_{\mathbf{2}}\mathbf{/}{\mathit{E}}_{\mathbf{\infty }}$)	$\mathit{\delta }$ (${\mathit{E}}_{\mathbf{2}}\mathbf{/}{\mathit{E}}_{\mathbf{\infty }}$)	$\mathit{h}$ (${\mathit{E}}_{\mathbf{2}}\mathbf{/}{\mathit{E}}_{\mathbf{\infty }}$)
CENT (2o)	0.1884 × 10−1/	0.1233 × 106/	0.2518 × 10−2/
	0.1661 × 10−5	0.4392 × 10−6	0.1067 × 102
VORO (2o)	0.1818 × 10−1/	0.7939 × 10−7/	0.2373 × 10−2/
	0.1628 × 10−5	0.4392 × 10−6	0.1051 × 102
BEST (2o)	0.1814 × 101/	0.7664 × 10−7/	0.2364 × 10−2/
	0.1628 × 10−5	0.4362 × 10−6	0.1050 × 102
CENT (1o)	0.4786 × 10−2/	0.3480 × 10−7/	0.6335 × 10−3/
	0.4576 × 10−6	0.1274 × 10−6	0.2721 × 101
VORO (1o)	0.4516 × 10−2/	0.2169 × 10−7/	0.5979 × 10−3/
	0.4302 × 10−6	0.9216 × 10−7	0.2667 × 101
BEST (1o)	0.4494 × 10−2/	0.2118 × 10−7/	0.5957 × 10−3/
	0.4278 × 10−6	0.8718 × 10−7	0.2668 × 101
CENT (0.5o)	0.1185 × 10−2/	0.9465 × 10−8/	0.1588 × 10−3/
	0.1164 × 10−6	0.4072 × 10−7	0.6603 × 100
VORO (0.5o)	0.1120 × 10−2/	0.5632 × 10−8/	0.1500 × 10−3/
	0.1051 × 10−6	0.2003 × 10−7	0.6625 × 100
BEST (0.5o)	0.1116 × 10−2/	0.5464 × 10−8/	0.1495 × 10−3/
	0.1042 × 10−6	0.1971 × 10−7	0.6631 × 100

), the $l_2$-norm error is a relative error and is defined as

$$E_2=\sqrt{\frac{{\sum }_i{\left(forecast-truth\right)}^2}{{\sum }_{i}trut{h}^2}},$$

and $E_{\infty }$ norm error is the maximum error over the domain and is given as

$$E_{\infty }=\underset{i}{\max }\left|forecast-truth\right|,$$

where $i$ is an index for all grid cells over a given domain. We reported the results in $E_{\infty }$ to show that error reduction is everywhere over the domain as expected for a second-order scheme, indicating that our implementation of these schemes was correct. Table 2 lists the 14th-day forecast errors of the Rossby–Haurwitz case at two-, one-, and half-degree resolutions against the true solution of the forced shallow water Equations (35)–(37).

The error reduction from coarse to one-level-finer resolutions is close to a factor of four, which matches the second-order accuracy expectation of these models over the rectangular latitudinal and longitudinal grids. Notably, the maximum errors in $E_{\infty }$ are also reduced by a factor of four. This study numerically showed that the schemes are second-order accurate anywhere in the domain and the numerical results reflect the correctness of the implementation. At a resolution of one degree, these models have relative errors less than 5% and 0.6% in terms of vorticity and fluid depth, respectively. Thus, the forecast fields are visually identical to those in Figure 3 and do not need to be plotted (the forecast error of vorticity is plotted in Figure 5).

By applying the averaged boundary conditions introduced in the previous section, we found that the spurious reflection was removed completely. Figure 6 shows a north–south cross-section of the 14th-day vorticity forecast at a longitude of 40° when using the CENT scheme. No noise was observed near the south and north boundary grids, which are labeled near grid indices of −90 (south) and 90 degree (north). The maximum error reductions given in Table 2 indicate that these averaging boundary conditions properly maintain the model accuracy for both the CENT and VORO Z-grid schemes.

For the Galewsky test case, Figure 7 and Figure 8 show the numerical forecasts of the vorticity and absolute vorticity at 96 h, 120 h, and 144 h. Compared to Figure 4 in [12], these figures closely resemble the features at all these hours. Thus, these Z-grid models performed appropriately and reasonably reflected the unstable jet over an isolated mountain. These models produced almost identical plots, and we thus obtained the numerical results from the CENT model’s half-degree resolution outputs.

In summary, we found that the centroidal Z-grid model could utilize time steps that are at least 1.4 times or longer in all of the resolutions for both test cases in this study. Of course, these longer time steps reduce the computing time by 30% or more at half-degree runs.

Note that our numerical experiments are based on the sequential codes of these Z-grid models. The runtime is a single-processor forecast time. As the models are almost identical, this computing-time reduction will certainly hold in a massively parallel computing environment. That is, a more than 30% efficiency improvement of the new model can also be achieved in parallel computing environments.

The new centroidal Z-grid model comprises more grid values used for the normal derivatives, and these may affect efficiency. It is indeed a minor addition as this happens only for calculation of the normal derivatives. Furthermore, the Poisson solver uses the majority of the computing time, even though a multigrid technique is used. For massive parallel computing, the computing involving a large number of grid values is of less concern. Model stability is beneficial in that it is significantly more viable for longer time-integration steps.

6. Conclusions and Remarks

A generalized Z-grid model was proposed and tested in comparison to the Voronoi Z-grid model. The proposed model has improved stability and efficiency. The comparison showed that the new centroidal model reduces the computation time by more than 30% compared to the Voronoi Z-grid model, while retaining the good features of a Voronoi Z-grid model, thus being more efficient. In addition, the small forecast errors at the end of the 14th day show that the centroidal model did not result in much energy loss while conserving all of the scalar variables, such as momentum, divergence, and vorticity, as well as mass fields.

For the weights of the centroidal Z-grid models, we tested different options. Specifically, we found that weights $w_t=w_b=1/8$ and $w_m=3/4$ were more stable than others on rectangular grids, maintaining relatively good dispersion relations. The numerical experiments were performed on a rectangular domain; however, we expect that a similar conclusion would be obtained for global icosahedral triangular Z-grid models. The regional domains of a plane and sphere differ with respect to the continuous spherical derivative operators; this will result in different coefficients of the shallow water equations. In our future study, we will focus on the centroidal Z-grid models with respect to a sphere.

The modification of the Rossby–Haurwitz test case and the methodology of this modification can be applied on a sphere and other test cases of shallow water equations. The advantage of this methodology is that it does not rely on the numerical solutions of a spectral model, and thus enables users to perform a clean comparison to the true solution. However, the realism of adding forcing functions to the governing equations may be of concern. To address this issue, it must be noted that in reality, forecast models or their equations are affected by the radiation of the Sun, for example, along with other physical forces. Hence, having test cases that satisfy the forced shallow water equations is reasonable.

Regarding time-integration schemes, the use of the Runge–Kutta fourth-order scheme in this study was mainly considered for stability comparison. In Z-grid model applications, a third-order Adams–Bashforth or other PBTI methods are more practical for saving computing time.