An Accurate Multimoment Constrained Finite Volume Transport Model on Yin-Yang Grids

With the rapid development of computer hardware involving hundreds of thousands of processors, it is now possible to use high-resolution general circulation models (GCM) to simulate weather and climate changes. In developing high-performance global models for such simulations, global computational meshes with quasi-uniform grid spacing are preferred. Representative meshes of this kind include the cubed-sphere grid (Sadourny, 1972), the icosahedral grid (Sadourny et al., 1968; Williamson, 1968), and the Yin-Yang grid (Kageyama and Sato, 2004).

Figure 1.  Schematic diagram of the Yin-Yang grid, which consists of two notched LAT-LON grids, i.e., the Yang grid and Yin grid, oriented perpendicularly to each other. Each component grid covers a domain of 270° and 90° in longitude and latitude directions.

By projecting an icosahedron onto a sphere, spherical icosahedral grids possess considerable uniformity. However, the icosahedral grid is unstructured, and numerical models on unstructured meshes require complicated numerical techniques that are very different from those used for structured grids. As a result, numerical models on unstructured grids are often more computationally expensive. Further, on unstructured grids, it is usually nontrivial to develop high-order schemes (more than second-order accuracy). The gnomonic cubic grid can be constructed by projecting an inscribed cube onto a sphere. It covers the whole globe with six identical patches. The major disadvantage of this grid is that the coordinate system is not orthogonal and discontinuous local coordinates along patch boundaries often introduce extra numerical errors. The conformal cubic grid was proposed by (Rancic et al., 1996) to obtain continuous local coordinates at patch boundaries, but at the cost of losing uniformity in grid spacing and analytical relations for transformations.

The Yin-Yang grid (Fig. 1) is a composite system, originally introduced by (Kageyama and Sato, 2004). A Yin-Yang grid is constructed with two perpendicularly oriented LAT-LON grid components of low-latitude regions; it is quasi-uniform and free of polar singularities. Both zones have grids in the local LAT-LON system, thereby providing a convenient platform for directly transplanting existing numerical models onto the LAT-LON grid. A Yin-Yang grid needs to communicate data in the overlapping area between Yin and Yang components, which means that conservation of total mass over the sphere is not automatically guaranteed. However, conservative schemes can be adopted within two-component grids. A transport model with rigorous conservation on Yin-Yang grids was proposed by (Peng et al., 2006); in their approach, flux computations for cells in the overlapping area are adjusted based on a piecewise constant reconstruction. Alternatively, the scheme adopted by (Lauritzen et al., 2010) for high-order mass integration over spherical polygons might be a potential choice for high-order implementations of the flux adjustment for global conservation of mass.

In this paper, a fourth-order transport model is built on the Yin-Yang grid using the recently proposed multimoment constrained finite volume (MCV) method (Ii and Xiao, 2009). This approach introduces a novel numerical framework for developing schemes of arbitrarily high order. Instead of using multimoments as model variables directly, degrees of freedom (DOFs) for high-order reconstruction are locally defined as pointwise values at specified points within each cell. Multimoments (constraints), including line-integrated averages (LIA), point values (PVs), and derivatives of different orders (DVs), are adopted to provide the spatial discretization formulations for local DOFs. The MCV model guarantees mass conservation within a Yin and Yang grid. Although there is no extra manipulation for conservation correction in the overlapping areas of a Yin-Yang grid, we have found (Li et al., 2008) that numerical errors, including those in mass conservation, can be markedly reduced by implementing high-order local reconstructions using multimoments. The numerical model is designed here with adequate accuracy for medium-range weather predictions or short-range climate simulations. A quasi-monotone model is constructed by using a TVB limiter in the present model. The TVB limiter designs a smooth indicator, which is useful for preserving high-order accuracy of the present model in smooth areas, while nonphysical oscillations in areas with a discontinuity or a large gradient are effectively removed.

The remainder of this paper is organized as follows. In section 2 we describe the fourth-order MCV formulations for solving transport equations on a Yin-Yang grid. Section 3 contains the numerical results from applying the proposed model to widely used benchmark tests; the proposed model is also compared with other representative global transport models for verification. A short conclusion is given in section 4.

s The Yin-Yang overset grid consists of two components: each is merely the low-latitude part of a conventional LAT-LON grid. As a result, the grid spacing of the Yin-Yang grid is quite uniform. To cover the whole globe, both Yin and Yang grids should cover at least the region from 45°S to 45°N in latitude and from 45°E to 45°W in longitude (that is, a domain of 270° in longitude direction). The two components of a Yin-Yang grid are normal to each other, so the transformation equations between Yin coordinates (λ’, φ’) and Yang coordinates (λ,φ) are easily written as

Bearing in mind that the Yin and Yang components are identical, the same numerical formulations can be independently used on either the Yin or Yang component only if the required data in the ghost cells are interpolated from the adjacent component. Locations of ghost cells are calculated and stored at the beginning of the computation. The MCV scheme is constructed based on a single-cell stencil. As a result, only one-layer of ghost cells is needed. Such a compact stencil not only improves computational efficiency, but it also suppresses extra numerical errors generated by nonconforming connections between Yin and Yang components. For the sake of brevity, we restrict discussions in the following sections to only the Yang grid.

In this subsection, we describe numerical formulations of the fourth-order MCV model for global transport computations on a Yin-Yang grid.

2.2.1 Governing equation

The two-dimensional spherical transport equation is written on the LAT-LON grid in conservative form as

where ,φ is the transported field, is the Jacobian of the transformation, α is the radius of the Earth, λ and φ denote the longitude and latitude directions on a sphere, e and f are the flux components inλand φ directions, is the flux vector, and is the angular velocity on the LAT-LON grid.

2.2.2 Configuration of DOFs

A cubic spatial reconstruction over each control volume is required to construct a fourth-order MCV model. As a result, for single-cell reconstruction, 16 local DOFs should be defined within a 2D rectangular control volume. As shown in the computational space in Fig. 2, within control volume , 16 local DOFs for model variable ψ(λ,φ,t) are defined at solution points P_i,j,m,n (denoted by solid circles) as

ψ_i,j,m,n (t)= ψ(λ_i,m,φ_j,n, t) (m=1…4, n=1…4), (3)

where i and j are indices of control volume, m and n are local indices of solution points P_i,j,m,n inλ andφ directions within each control volume, and are the locations of solution points.

The discretization formulations for these local DOFs, based on constraint conditions on multimoments, were developed in (Ii and Xiao, 2009). Hereafter we describe the implementation of a fourth-order MCV scheme for global transport computations on a Yin-Yang overset grid.

To construct a high-order model, the derivatives of the flux vectors [the second and third terms on the left-hand side of Eq, (2)] are first computed pointwise at the solution points, using the MCV scheme, to obtain semi-discrete equations. Then time integration is performed using a fourth-order Runge-Kutta method.

2.2.3 Spatial discretization of DOFs

On structured grids, a fully multidimensional MCV scheme can be constructed by efficiently implementing one-dimensional formulations in different directions. Without losing generality, we consider the computation of the derivative of flux component e with respect toλto discretize the governing equation in theλ-direction as

A similar procedure can be used along grid lines in the φ-direction by exchanging e with f and λ withφ.

Figure 2.  For the fourth-order MCV model, equidistant local DOFs are defined within each control volume . where 16 local DOFs are denoted by the solid circles. Within each control volume there are four rows (grid lines) in the φ-direction and four columns (grid lines) in the λ-direction.

Numerical formulations are described by considering the (nth=1…4) row of control volume C_i,j (Fig. 2). In the following discussion, we omit the subscripts j and n for sake of the conciseness. A cubic reconstruction polynomial can be built for any physical field using four local DOFs. For example, for model variableψ(λ,φ_j,n, t) we get

, where l_m(λ) is a Lagrange basis function,

To build a MCV model, we first define the multimoment constraints, which are used to derive the spatial discretization equations for the DOFs. Three kinds of moments are adopted to construct a fourth-order scheme:Line-integrated averages (LIA) over the grid line are defined by

Point values (PVs) at the two endpoints of the grid line are defined by

First-order derivative values (DV) with respect toλat the center of the grid line are defined by

where subscript C denotes the center of the grid line.

The four DOFs and multimoment constraints are connected through the reconstruction polynomial Eq. (5)

Once we obtain formulations for updating the multimoment constraints in terms of the DOFs at the present time step, semidiscrete equations to predict the unknown DOFs can be derived by differentiating Eq. (9) with respect to time. We describe the discretization of the governing equations for the constraints as follows.

where e_i1 and e_i4 are pointwise values of flux component e at the two endpoints of the grid line.

Using the reconstruction profile Eq. (5), the transported field is continuous at cell interfaces. As a result, e_i1 and e_i4 can be directly computed by DOFs defined there.

The PV moment is updated by the differential-form Eq. (4) written at the two endpoints of the grid line. For example, at the left endpoint it is written as

where E is the spatial reconstruction for flux component e having the form in Eq. (5).

and e_ic=E_i(λ_i).

Finally, the time evolution of DOFs is computed by

To solve the problem with discontinuities or large gradients, a limiting projection using the TVB concept (Shu, 1987) is adopted in the MCV scheme to modify the DOFs to

where L(λ) is a linear reconstruction function written as

L(λ)=b₀+b(λ-λ_i) (17)

Here, b₀ equals the value of LIA to assure conservation of the MCV model, and the slope b is determined using the Superbee limiter,

with functions and defined as

and , , .

In the control volume where DOFs are modified by the above limiting projection, the continuity of transported field at cell interfaces cannot be preserved. Under this circumstance, we solve the Riemann problem to determine the pointwise value of flux component e. For example, e_i1, mentioned above, is computed by

2.2.4 Time stepping

After accomplishing spatial discretization, the ordinary differential equation are obtained from

So far, we have developed a fourth-order transport model on a Yin-Yang grid using the MCV scheme, which is proposed under the Eulerian framework. In earlier work, we implemented numerical models on a Yin-Yang grid for a transport model (Li et al., 2006) and for a shallow-water model (Li et al., 2008) using semi-Lagrangian schemes. Although semi-Lagrangian schemes are widely used in GCMs, it is not easy to treat the integration of source terms along the trajectory with uniformly high-order accuracy. Therefore, for general systems of conservation laws, we prefer to adopt high-order schemes developed under the Eulerian framework.

In this section we present results from convergence tests and from some standard and widely used benchmark tests to evaluate the MCV transport model on the Yin-Yang grid. Using the error measurements defined in (Williamson et al., 1992), the global integral on the unit earth is defined by

where φ_Ya and φ_Yi, are numerical solutions of the transported field on Yang and Yin components, and the normalized global errors are,

where the subscript T denotes the true value, and φ_Ya,T and φ_Yi,T are the true solutions on Yin and Yang components.

A smooth initial distribution was specified for accuracy tests by

where the parameter α is the angle between the axis of the solid body rotation and the Earth's axis, and u₀=2πα/12.

Normalized errors (defined following Williamson1992) in the present model on a series of refined grids are given in Tables 1-3 for solid rotations in different directions. According to these results, fourth-order accuracy of the present model is achieved, as expected.

Figure 3.  Numerical and exact solutions for the cosine bell advection test at day 12 on a grid with a DOF resolution of (4/3)°. Solid lines show numerical results and dashed lines show exact solutions. The contour intervals are 100 from 100 to 900.

y Advection of a cosine bell, the first case in Williamson's set of benchmark tests (Williamson et al., 1992), is widely checked to evaluate global models. The cosine bell is initially specified by

where h₀=1000m and r=acos^-1[sinφ₀sinφ+cosφ₀ cosφcos(λ-λ₀)] is the great circle distance between initial center (λ₀, φ₀)=( π/2,0) and point (λ, φ). The driving wind field was the same as in the previous test. Grid spacing was 4° in this case, which corresponds to a grid resolution of (4/3)° in terms of DOFs.

Figure 3 shows that the cosine bell returns to its initial position after one revolution (12 days). The numerical contour lines (dashed lines) are visually identical to the analytic solution (solid lines). The history of normalized errors is shown in Fig. 4 where fluctuations in the normalized errors are observed when the cosine bell moves across the overlapping region of the Yin-Yang grid. Table 4 gives the normalized errors after one revolution. Normalized errors along the polar direction are larger than the others because, in this case, the advection path has more overlapping regions.

Compared with other models developed on a LAT-LON grid [semi-Lagrangian inherently conserving and efficient (SLICE) scheme proposed by Zerroukat et al., 2004), our result is competitive along the polar direction. Although a smaller Courant-Friedrichs-Lewy (CFL) number is used compared to that in the SLICE scheme, the present model easily treats source terms with uniformly high-order accuracy when further developing the shallow-water model and a three-dimensional dynamic core on a sphere. In addition, polar problems are avoided by utilizing the quasi-uniform Yin-Yang grid.

Numerical results from the fourth-order MCV model are better than those from the third-order conservative semi-Lagrangian model (Li et al., 2008) with an equivalent spatial resolution; for example, compared with the numerical tests under a 1.40625° resolution in (Li et al., 2008). In addition, numerical results from the present model are obviously better than those in (Baba et al., 2010) with a resolution of 1° on a Yin-Yang grid and better than those published in (Li et al., 2006) (see the bottom panel in their Fig. 5) with a resolution of 1.25° on a Yin-Yang grid using the cubic Lagrange interpolation polynomial. Our simulation re sults on the Yin-Yang grid are also competitive with other multimoment fourth-order global models, such as on the icosahedral grid (Ii and Xiao, 2010) and on the cubed sphere (Chen and Xiao, 2008) using similar degrees of freedom.

A more practical and challenging advection test was introduced by (Nair and Jablonowski, 2008). In this test, a complex velocity field is generated by superposing two vortical deformational flow fields whose vortex centers move along the trajectory of the solid rotational flow specified by Eq. (23). A detailed description of this test can be found in (Nair and Jablonowski, 2008).

We ran this test on a grid with a resolution equivalent to 2.5° in terms of DOFs. Rotations of the vortex centers in three different directions, corresponding toα=0，α=π/4, and α=π/2, were checked. The time evolutions of normalized errors along the three different directions are depicted in Fig. 5, the errors are hardly influenced by the direction of rotation.

Figure 4.  Normalized errors from the solid-body test on a grid with a DOF resolution of (4/3)°. Normalized global errors are defined in Eq.(23).

As shown in (Nair and Jablonowski, 2008), a discontinuous Galerkin (DG) method, a semi-Lagrangian method, and a finite volume (FV) method have been evaluated using this moving vortex test. Compared with the semi-Lagrangian method (see their Fig. 5 for normalized errors), the normalized errors of the MCV scheme are obviously better when equivalent resolutions are used with respect to the conventional LAT-LON sphere. The present model also gives smaller errors than the FV solution on a grid with the same resolution (see their Fig. 7 for normalized errors). Although the DG method shows a somewhat smaller normalized error under an equivalent resolution, the MCV scheme allows a larger CFL number for computational stability (see Table 1 in Ii,et al., 2009). Our MCV scheme also has a competitive computational efficiency since it does not involve any time-consuming numerical integrations. Contour plots of the transported field φ and absolute errors are shown in Fig. 6 at different days. The location of moving vortex centers and the phase of the front of deforming vortices were accurately captured.

Figure 5.  Normalized errors for the moving vortex test on a grid with a DOF resolution of 2.5°. Normalized global errors are defined in Eq.(23).

Finally, we present numerical results of a nonsmooth profile transported by a deformational flow (Nair et al., 1999). As expected, a high-order scheme always needs a limiting projection to suppress numerical oscillations. Hereafter, we refer to the model with limiting projection Eq. (16) as the limited model and that without Eq. (16) as the unlimited model.

The wind field was as specified in previous tests, except the two vortex centers were fixed. This test includes a sharp step at the initial time t=0. Details of the setup of this test can be found in (Nair et al., 1999). The computational parameters adopted in this study wereγ=3, δ=0.01, with a grid resolution of 3°.

Figure 7 shows the numerical solutions at t=3 from the unlimited (Panel a) and limited (Panel b) models, as well as the exact solution (Panel c). Both unlimited and limited models give the phase of the moving front in good agreement with the exact solution. However, the numerical solution from the unlimited model contains nonphysical oscillations near the sharp step. Those oscillations are effectively removed in the limiting projection Eq. (16), as shown in Panel b. The maximum and minimum values are φ_max=2.12 andφ_min=-0.066 in the numerical solution from unlimited model and φ_max=2.018 andφ_min=-0.0019 in that from the limited one.

As a clear illustration of the performance of the limiting projection adopted in present study, Figure 8 shows numerical solutions on a cross section from both the limited and unlimited models. Nonphysical oscillations are invisible on the cross section from the limiting projection, while the steepness is well-preserved, compared to the exact solution.

A fourth-order transport model on the spherical Yin-Yang grid has been developed using the multimoment constrained finite volume (MCV) scheme. In the present model, the polar singularities and convergence of meridians in polar regions of conventional LAT-LON girds are completely avoided by using the quasi-uniform Yin-Yang grid. As model variables, the MCV method uses local DOFs defined within single elements at equally spaced points. In our transport model, the evolution formulations for local DOFs are derived through constraint conditions on moments of line-integrated averages (LIA), point values (PV), and first-order derivative values (DV). The compact local stencils used by the MCV scheme not only contribute to computational efficiency when dealing with information exchange between Yin and Yang components, but they also help reduce extra numerical errors due to overlapping connections between two grid components. Compared with other local high-order schemes, such as the discontinuous Galerkin scheme and the spectral volume scheme, an MCV scheme offers improved computational efficiency, since it allows less restrictive CFL stability conditions, and it does not require explicit volume quadrature. In addition, a limiting projection can be adopted to control spurious numerical oscillations when a discontinuity or large gradient exists during tracer transport. Widely used benchmark tests show that, compared with most existing models, the present model is competitive both in high-order numerical accuracy and computational efficiency.