Challenges in Developing Finite-Volume Global Weather and Climate Models with Focus on Numerical Accuracy

Yuanfu XIE and Zilong QIN

1 State Key Laboratory of Severe Weather, Chinese Academy of Meteorological Sciences, China Meteorological Administration, Beijing 100081

2 Guangdong-Hong Kong-Macao Greater Bay Area Weather Research Center for Monitoring Warning and Forecasting (Shenzhen Institute of Meteorological Innovation), Shenzhen 518048

ABSTRACT High-resolution global non-hydrostatic gridded dynamic models have drawn significant attention in recent years in conjunction with the rising demand for improving weather forecasting and climate predictions.By far it is still challenging to build a high-resolution gridded global model, which is required to meet numerical accuracy, dispersion relation, conservation, and computation requirements.Among these requirements, this review focuses on one significant topic—the numerical accuracy over the entire non-uniform spherical grids.The paper discusses all the topic-related challenges by comparing the schemes adopted in well-known finite-volume-based operational or research dynamical cores.It provides an overview of how these challenges are met in a summary table.The analysis and validation in this review are based on the shallow-water equation system.The conclusions can be applied to more complicated models.These challenges should be critical research topics in the future development of finite-volume global models.

Key words: global weather forecast model, finite-volume methods, model accuracy, irregular spherical grids, C-grid and Z-grid schemes

1.Introduction

Modern computing power has enabled the realization of global weather forecasting and climate prediction models at storm- and cloud-resolving scales.Even with success in weather forecasting under present operations by ECMWF (e.g., Simmons et al., 1989) and NOAA(e.g., Kanamitsu et al., 1991), spectral models face difficulties at finer resolutions at which non-hydrostatic phenomena occur.When compared at higher resolution with gridded models, they have significant issues in handling convective storms, mass conservation, nesting, and positive scalar integration (e.g., Williamson and Rasch, 1994),and they are too expensive in computing.Undoubtedly,gridded numerical models are preferred when dealing with the resolutions at the convective scale.

The finite-volume method (FVM) naturally conserves scalar states in gridded numerical models, which is crucial in predicting fine-scale tracers such as aerosols.Instead of discussing finite-difference or finite-element methods (FDM, FEM), this review focuses on the challenges in developing a numerically accurate FVM based global gridded model.The review analyzes issues in meeting accuracy requirements from mathematical and geometrical points of view.We also propose a table that compares different finite-volume schemes, explaining their behaviors observed in numerical experiments (e.g.,Reinecke and Durran, 2009; Yu et al., 2020).Although the discontinuous Galerkin (DG) method (Giraldo et al.,2002; Bao et al., 2014) warrants potential applications in numerical weather forecasting with high-order accuracy and conservation, it is beyond the scope of the traditional FVM scheme discussed in this paper.We also concentrate our discussion on the conventional FVM, while those FVMs with Riemann solvers are also beyond the scope of this study.

Evaluation of the numerical performance of a dynamical core is difficult due to the sheer variety of factors,such as sound wave treatment, thermodynamics, physics,and many more.Even for simple shallow-water equations on a sphere, due to the complexity of spherical grid structures, the challenges in developing a second-order—even a first-order—accurate FVM scheme are still serious as reviewed in this study.Thus, this paper mainly focuses on the accuracy issue, along with a few related challenges of conservation, dispersion, and computation,without involving other more complex issues.

When an FVM based numerical scheme is applied on a non-uniform mesh (e.g., a grid over the global sphere),achieving certain numerical accuracy remains a significant challenge, although it is well-addressed on a plane.This is because no regular (uniform) grid evenly divides a sphere.To avoid any confusion, we define a regular cell as its vertices are concyclic and equally spaced; otherwise, the grid cell is called an irregular cell, and an irregular grid is called if it contains at least one irregular cell.Unlike other reviews (e.g., Saito et al., 2007; Ullrich et al., 2017; Konor and Randall, 2018), this paper mainly discusses the challenges with respect to the numerical accuracy on the spherical irregular grids in the mathematical and geometrical perspectives.Other important properties of dispersion, computational mode,and conservation are also briefly discussed for readers to assess an FVM model comprehensively.

Since the numerical accuracy can be easily confused with model resolution, we briefly clarify the two concepts: the resolution of a gridded model is directly related to the grid cell size, e.g., a 1-km resolution; while accuracy is defined as whether numerical solutions converge to the continuous solution of its corresponding governing equations and how fast they converge as the resolution increases.An increasing resolution could improve the numerical solutions of a convergent scheme.Inaccurate numerical schemes without convergence may yield poorer results with increasing resolutions.We show in the Appendix how an inaccurate scheme introduces the numerical errors in a simple continuous governing equation.Moreover, even if a scheme is nonconvergent over only a small portion of all grid cells, its global solutions will not converge to the true solution of the continuous governing equations either.Moreover, local numerical errors will propagate in a dynamic system.Errors in some regions may contaminate other areas over a certain amount of time.

Proving convergence of numerical solutions is usually difficult.There are two steps to ensure model convergence.The first step is to ensure that a given discretization scheme is consistent with its continuous governing equations, known as consistency.The second step is to stabilize this scheme, which is known as stability.Consistency and stability lead to convergence based on the Lax equivalence theorem (Lax and Richtmyer, 1956).Depending on the order of consistency, a scheme converges at a certain rate, e.g., at a first- or second-order convergence rate.Consistency is a necessary condition of a convergent scheme.For simplicity, we refer to consistency when discussing accuracy, i.e., a numerical model is termed “accurate” when it is consistent everywhere within a domain of interest.Otherwise, it is termed “inaccurate.”

Since no known grid can uniformly mesh a sphere,this limits our ability to extend the application of FVM schemes from a plane to a spherical surface while retaining a certain order of accuracy.Two widely-used quasiuniform spherical grids are the cubed sphere and icosahedral grids (see Fig.1).Their grid cells do not possess the regularity of rectangular, triangular, hexagonal, or pentagonal cells.Therefore, the discussion in this review is based on general irregular grid cells and their related stencils.Also, the conclusions are applicable to any other spherical grid structures as well.

Based on these grids, several well-known FVM schemes are candidates with different state variable staggering.These schemes have long been used in numerical models and are illustrated in Fig.2, which shows the positions of variables in the shallow-water equation system discussed in Section 2.Among these schemes, some have been implemented for global weather forecasting models,including directional splitting (Lin and Rood, 1996; Lin,2004), which estimates the two-dimensional (2D) conservative field with two one-dimensional (1D) directional flux-form operators; A-grid (Satoh et al., 2008), adopted by the Nonhydrostatic Icosahedral Atmospheric Model(NICAM), where a spring dynamic-based technique is used to modify the grid to improve its numerical accuracy; C-grid (Ringler et al., 2010; Wan et al., 2013),where Ringler et al.(2010) used an icosahedral hexagon-pentagon grid to implement their Model for Prediction Across Scales (MPAS), and Wan et al.(2013)developed the Icosahedral Nonhydrostatic model (ICON)with an icosahedral triangle grid; C-D grid (Adcroft et al., 1999), which places both wind components at an edge center; and Z-grid (Heikes and Randall, 1995a),which uses the divergence and vorticity as prognostic variables, and poses them all at a cell center along with other conservative scalars.We discuss how the accuracy challenges are met for these schemes, which may help readers assess and improve the existing models or develop future models.

In Section 2, we review challenges in grid geometry that affect numerical accuracy.Section 3 summarizes challenges associated with dispersion, computational modes, and conservation.We build in Section 4 a matrix that evaluates some of the widely used FVM schemes.Section 5 provides a summary of this review.

Fig.1.Two quasi-uniform grids on a sphere, with the icosahedral grid on the left and the cubed sphere on the right.

2.Challenges with spherical geometry

In this study, we present the challenges for constructing an accurate FVM for a shallow-water system on a spherical grid.First, we review a shallow-water equation system in two different forms.

Fig.2.Finite-volume scheme definitions for the A-, B-, C-, C-D, D-, and Z-grid schemes, where h is the fluid thickness, V is the velocity vector, un is the normal velocity, uτ is the tangential velocity component, and ζ and δ are the vorticity and divergence, respectively.

A velocity vector-invariant form of the shallow-water equations (Ringler et al., 2010) is given as:

2.1 Irregularity of cells reduces accuracy(“accuracy drop”)

Equations (10) and (11) transform the operators into the line integrals over a curved cell polygon (spherical patch on the global grid).For each cell on the grid, a set of piecewise line integrals along the edges can be used to approximate the right-hand side (RHS) of Eq.(10) or (11):

Table1.Line integral errors in the approximation of a Laplacian operator, where the top row gives the edge lengths and the remaining rows give the errors in each scheme

2.2 Accuracy challenges

Since there is no mesh that can evenly subdivide a sphere, every known grid is irregular, and an FVM designed on these irregular grids suffers the “accuracy drop.” Designing an accurate FVM scheme on a sphere is challenging.We review the challenges due to the irregular spherical geometry.

Challenge 2.1: Which cell center is the best choice to pose the state variable?

For an FVM based scheme on a grid of convex polygonal cells, one has to select a location of its discrete cells to pose its function value as an approximation of the mean value, e.g.,defined above.Over a polygon, several centers exist, i.e., centroid, icosahedral, Voronoi, circumcenter, and mid-centers, as shown in Fig.3.An icosahedral center is the icosahedral grid point location,which derives from the recursive subdivision of an icosahedron with triangular grids on a sphere.The mid-center is the intersection point between two line-intervals that connect the middle points of two opposite edges of a quadrangle-based grid.For modern global models,MPAS uses centroid center through a centroidal Voronoi tessellation (Skamarock et al., 2012), NICAM uses icosahedral center (Satoh et al., 2008), a Z-grid model uses Voronoi center (Heikes and Randall, 1995a), ICON(Wan et al., 2013) uses circumcenter, and Finite-Volume Cubed-Sphere Dynamical Core (FV3; Lin and Rood,1996) uses mid-center.These centers may be at different locations if a polygon is irregular, which leads toChallenge 2.1for an FVM scheme on a sphere (note that all of these centers coincide at the same point if the underlying polygon is regular).

Fig.3.Grids of centroidal center (top left), Voronoi center (top right), icosahedral center (bottom left), circumcenter (bottom middle), and midcenter (bottom right).

Because only a function value at a centroidal center is a second-order approximation of the mean value over the cell,the mean values are merely first-order accurate at any other centers.The distance,d, between these centers may be bounded by a positive constant, ?, multiplied by the edge length,l, for an irregular cell (e.g., some cells inside either an icosahedral or cubed sphere grid):

So, the ratio ofdandlmeasures how uniform a polygon is.Based on the mean value theorem for integration,a centroid center is the only one that guarantees a second-order accuracy of the FVM mean valueover a grid cell if it is selected to place the function values.Other such centerings for integrand values can yield first-order accuracy at most.This issue rises asChallenge 2.1.

Challenge 2.2: The center-arc may not bisect the edge.How to construct approximations of the integrands of Eqs.(10) and (11) to be accurate?

For two convex cells that share an edge, the line interval that connects their centers is referred to as a centerarc, which normally intersects the edge.In regular cells,the edge bisects and is perpendicular to the center-arc.For irregular cells, these properties no longer hold, which causes other accuracy issues for FVM schemes.

First of all, only a center-arc of two circumcenters bisects the shared edge.We use the definition of λ: the distance between the intersection and middle point of this edge as shown in Fig.7 of Heikes et al.(2013) copied here as Fig.4.For standard icosahedral or centroidal Voronoi tessellation-based grids (CVT), the ratio between λ and the edge length,l, has a positive lower bound approximately equal to 0.088 (Heikes et al.,2013):

for non-Voronoi center-arcs, which may reduce the accuracy of non-Voronoi center schemes to zeroth-order as well.This also poses a challenge for maintaining accuracy of finite-difference derivatives if function values at only two centers are used.

Challenge 2.4: A center-arc may not be perpendicular to the edge.How to construct an accurate FVM?

The center-arcs of the Voronoi centers and circumcenters are perpendicular to the edges, whereas the centerarcs for the centroidal, icosahedral, or mid-centers are not perpendicular.This is another challenge in designing an accurate FVM.Directly using function values between these centers yields an accurate first-order scheme.Following the accuracy drop discussed earlier, it may also result in a zeroth-order FVM scheme.

Fig.4.Illustration of center-arc and edge.The center-arc may not bisect the edge.

Fig.5.Definition of λ⊥ of the distance between an arc-center and its edge intersection.

Challenge 2.5: For irregular polygons,how to construct accurate information along with cell edge’s tangential directions,such as tangential velocity for Eq.(1)or Jacobian in Eqs.(3)-(5)?

In general, no scheme defines the vector state variables at vertices on a sphere due to tangential velocity along one edge that does not match either the tangential or normal velocity along its adjacent edges if the underlying grid is irregular.Similarly, tangential derivatives along edges do not match either tangential or normal derivatives from other edges (see in Fig.6 the black normal and tangential vectors and the thin red normal and tangential ones).Thus, the ability to construct tangential information accurately, e.g.,k×Vfor Eq.(1) orJfor Eqs.(3)-(5) at the edge centers is another challenge.Both constructing tangential information through a linear combination of adjacent normal and tangential vectors at other edges like the Thuburn-Ringler-Skamarock-Klemp (TRiSK) scheme for MPAS (Skamarock et al.,2012), and interpolating function values at neighbor cells to vertices for tangential information (a Z-grid model;Heikes and Randall, 1995a) result in significant loss of numerical accuracy.This challenge seems the most difficult one for all FVM models.

An accurate global numerical forecasting model needs to address all these challenges.Several state-of-the-art FVM schemes are designed to overcome a few of these challenges, but not all of them.Thus, achieving numerical accuracy of an FVM scheme on a sphere remains a significant challenge to the numerical modeling community.

Fig.6.Challenge of defining tangential derivatives on irregular grid because the normal and tangential vectors (black and red vectors) do not match and interpolation to vertex location (blue dots) causes accuracy loss.

3.Challenges with physical requirements

Discretization schemes for the equations of atmospheric dynamics require that physical properties encapsulated by the equations be reflected in whole or in part by the numerics.Properties such as the proper dispersion relations and conservation at least are often crucial for numerical accuracy requirements.Minimization of computational modes is usually required of discretization in the physical space to prevent spurious, non-physical artifacts.We focus on these specific properties in this section, while cautioning that additional mimetic properties like conservation of angular momentum or potential enstrophy can also be important in certain applications.We summarize these challenges based on previous other studies on regular grids for convenience in assessing an FVM model’s quality.These challenges are even more serious for spherical grids, and more theoretical research is required on this topic.We begin with the dispersion relation, which is one of the most crucial physical properties to be considered while maintaining model numerical discretization accuracy.

Challenge 3.1: How to design a numerical scheme with a good dispersion relation independent of horizontal grid sizes?

A dispersion relation describes how close a temporal and spatial spectral relationship of an FVM scheme approximates the wave solution of the continuous governing equations.A poor dispersion relation causes the initial condition to rapidly deviate from accurate predictions (Adcroft et al., 1999; Thuburn et al., 2009; Konor and Randall, 2018).

Dispersion relations are complicated to analyze for a full atmospheric dynamic model.Linearized rotating shallow-water equations are usually adopted to analyze a dispersion relation (Randall, 1994; Rajpoot et al., 2012).Recently, Konor and Randall (2018) performed dispersion analysis of a more sophisticated linearized non-hydrostatic model for inertia-gravity modes, which resulted in similar but more precise conclusions, which is not discussed in this review.Dispersion analysis with linearized rotating shallow-water-equation system starts from Eqs.(1)-(2):

is the deformation radius (Gill, 1982), andfis the Coriolis parameter.For a discretization scheme applied to Eqs.(18)-(20), how close its discrete form approximates Eq.(23) may or may not depend on the γ/dratio, wheredis the grid size or distance between centers.Dispersion relations become unsatisfactory if γ/dis too small for certain schemes (Randall, 1994).

For modern numerical prediction models, grid sizes ofdare small enough to resolve the deformation radius, i.e.,the γ/dratio is sufficiently large if γ is a constant as defined in Eq.(24).However, a shallow-water equation system can be regarded as an arbitrary vertical mode of the full primitive equation (Salmon, 2007).The general deformation radius defined by Eq.(24) is also proportional tokz?1, wherekzis the vertical wavenumber (Gill,1982; Salmon, 1998).This indicates that the γ/dratio can be small, even for small grid sizes.Konor and Randall (2018) show how the vertical modes impact dispersion relations for a full three-dimensional inertia-gravity model.So far, these dispersion relations of various FVM are analyzed on regular grids.It would be an interesting topic to examine the FVM on irregular grids.

Challenge 3.2: Different numbers of state equations may cause additional computation modes.How to design an FVM with an equal number of equations for all state variables?

The computational mode is distinguished from the physical mode, which is the misrepresentation of linear wave propagation by numerical discretization schemes.Several sources cause computational modes, but certain sources may not be directly related to a spherical grid.This review is mainly concerned with the computational modes caused by spherical-grid structures.

To reduce computational modes, we should have an equal number of equations for all state variables; otherwise, we may introduce computational modes when different numbers of equations are discretized for different model state variables.Heikes et al.(2013) have already discussed the computational modes for some FVM schemes.For example, over an icosahedral grid on a sphere, there are three times as many edges as cells.If one model discretizes some of its model states on edges and some others on cell centers, they transfer more or less information from one state to another.This allows computational modes and should be avoided.This challenge essentially prohibits staggering grids with state variables at different geometric locations, say edges and centers.

Challenge 3.3: How to design an FVM conserving total energy and enstrophy for shallow-water equations?

With increasing societal demand for the environmental predictions of aerosols or other chemical tracers, scalar conservation is necessary for numerical models.FVM schemes are derived from a conservation form of the continuous equations, which naturally conserve scalar quantities, such as moisture, mass field, and other chemical tracers.However, the most challenging task is to conserve energy and enstrophy for shallow-water equations.

There is no full three-dimensional global gridded model that currently meets this challenge.In fact, the real atmosphere is not a perfect energy-conserving system either because of friction, radiation, and microscale physics.For shallow-water equations, a numerical FVM meeting this challenge would be of significant interest(Sadourny et al., 1968).

4.State-of-the-art finite-volume schemes

We examine several state-of-the-art FVM schemes,such as A-grid, C-grid, and Z-grid, and evaluate how these schemes resolve specific challenges and the areas that require improvement.Finally, we summarize them in a matrix listed in Table 2 and explain how they meet and fail the challenges in Table 3.

4.1 A-grid

An A-grid scheme places all state variables at the same location, e.g., the icosahedral centers, such as in the NICAM (Satoh et al., 2008).Their scheme successfully runs at a global 3-km resolution, with a grid modified by an algorithm called spring dynamics (Tomita et al., 2001,2002).Spring dynamics helps the scheme overcomeChallenge 2.1by moving the icosahedral centers toward the centroidal centers.When calculating the flux across the edges, Tomita et al.(2001) projected the vectors for the velocity and pressure gradient onto the normal direction to prevent accuracy loss at the edges; this scheme passesChallenge 2.4.The center-arc and edge are not perpendicular to each other.Challenges 2.2,2.3,and 2.5are problematic for an A-grid scheme.The center-arc and edge do not bisect each other and cannot overcomeChallenges 2.2 and 2.3, as listed in Table 2 of Heikes et al.(2013); i.e., there are positive constants of , so that Eqs.(16) and (17) hold.This results in the first-order accuracy of these values or vectors.Due to the accuracy drop,the scheme has zeroth-order accuracy at those distorted cells.For the tangential velocity components, Tomita et al.(2001) first used a weighted interpolation to approximate the velocities at the vertices [see Eq.(10) in their study].Next, they took an average of these velocities to approximate the tangential velocity at the edge center.Overall, this tangential velocity calculation is also firstorder accurate.Thus, an A-grid scheme failsChallenges 2.2,2.3,and 2.5.

The A-grid schemes do not possess a good dispersion relation, as shown by Randall (1994).Since A-grid schemes place all model states at the same location, e.g.,the icosahedral centers, A-grid schemes can passChal-lenge 3.2.There are no previous studies that demonstrate A-grid’s energy and enstrophy conservation for shallowwater equations.The A-grid row of Table 2 lists how well an A-grid scheme overcomes these challenges.

?

Table2.Matrix of all the discussed schemes

4.2 C-grid

A C-grid scheme places its mass field at the cell centers and the normal velocity field on the edges.Two Cgrid models have been implemented in previous studies.The ICON (Wan et al., 2013) uses circumcenters and the MPAS (Skamarock et al., 2012) uses Voronoi centers.MPAS based on the CVT C-grid scheme can meet some more challenges than the circumcenter one.The circumcenter C-grid failsChallenge 2.1, but the CVT C-grid passes the challenge.Although a CVT is used in MPAS,it may not always be obtained using the algorithm reported in Du et al.(2003).

CVT schemes failChallenge 2.2, while circumcenter schemes meet it.Conversely, CVT schemes meetChallenge 2.3, but circumcenter schemes fail.Both of them do well forChallenge 2.4.Although Thuburn et al.(2009) and Ringler et al.(2010) used special combinations of normal velocities or fluxes at the adjacent edges to construct the tangential derivatives to satisfy certain physical requirements (called the TRiSK method), the tangential information obtained was not accurate.Challenge 2.5remains problematic for both of the C-grid schemes.

The main issue with the C-grid schemes regarding the dispersion relation is that this relation depends on the deformation radius.Thus, both C-grid schemes failChallenge 3.1.

As Heikes et al.(2013) pointed out, a C-grid scheme(except the quadrilateral grid) has three times the number of discrete equations for velocity than that of the thickness field because of the cell and edge ratio if the underlying grid is an icosahedral grid.Obviously, it failsChallenge 3.2.For a quadrilateral grid, the velocity state and mass state possess the same numbers of equations because the numbers of edges and cells are basically the same.

To conserve energy and enstrophy, Ringler et al.(2010) constructed a tangential flux by selecting a special combination of normal fluxes from neighboring cells.This particular C-grid conserves both energy and enstrophy.The matrices of the C-grid models are given in the C-CVT and C-circum rows in Table 2.

4.3 C-D grid

Adcroft et al.(1999) introduced a C-D grid to avoid grid-scale noise associated with an early version of the C-grid used for coarse resolution oceanic models.The C-D grid is also adopted by the well-known Geophysical Fluid Dynamics Laboratory (GFDL) FV3 (Lin and Rood,1996; Lin, 2004; Putman and Lin, 2007).With the tangential velocity reconstruction scheme reported in Thuburn et al.(2009), the C-grid scheme has overcome the noise issue, rendering it more preferable to the C-D grid.Thus, we do not discuss the C-D grid here.

4.4 Z-grid

Unlike other FVM schemes, a Z-grid scheme solves Eqs.(3)-(5) and places the vorticity and divergence, as well as other variables, at the Voronoi centers of the icosahedral hexagon/pentagon grid.Hereinafter, this implementation is referred to as the Voronoi Z-grid scheme.This scheme solves Poisson’s equations [i.e., Eqs.(6)and (7)] to convert vorticity and divergence into streamfunction and velocity potential at every time step, respectively.

A Voronoi Z-grid does not meetChallenge 2.1, which introduces the first-order error.Heikes and Randall(1995b) developed a tweaking algorithm that modifies the underlying grid to overcomeChallenge 2.2.They showed that the tweaking algorithm improved the accuracy of certain operators, e.g., the Laplace and flux-divergence operators but does not help the Jacobian operator.A Voronoi Z-grid scheme certainly meetsChallenges 2.3 and 2.4.However, the major issue with a Voronoi Z-grid scheme is its ability to overcomeChallenge 2.5, as explained in the following paragraph.

Old Tnate Sanna would open the door to the rather frightened little messenger and would usher5 him-or her - into the dark voor-kamer, where the shutters6 were always closed to keep out the heat and the flies

A Z-grid scheme calculates a Jacobian operator using the following relation:

Table3.Challenge list

To approximate the tangential derivative, a Voronoi Zgrid scheme uses either an equal or area-weighted interpolation to approximate the function values at the vertices of an edge, after which it uses a centered finite difference to approximate the tangential derivative at the edge center.The function interpolation of these equal or area-weighted schemes results in no higher than a firstorder accuracy.Due to the accuracy drop, the Jacobian operator in Eq.(25) will only result in zeroth-order accuracy.The tangential derivative remains a significant challenge for the Voronoi Z-grid scheme, as discussed in Eldred and Randall (2017).Thus, it failsChallenge 2.5.

Previous studies have demonstrated that a Z-grid scheme possesses the best dispersion relation among all finite-volume schemes (Randall, 1994; Konor and Randall, 2018) and does not depend on the deformation radius,which meetsChallenge 3.1.Since it has the same number of discretized equations for all the state variables, Zgrid models overcomeChallenge 3.2.

Canonical implementations of Z-grid schemes conserve both the energy and enstrophy for shallow-water equations (Salmon, 2007; Eldred and Randall, 2017), although they use slightly different forms of Eqs.(3)-(7).The last row of Table 2 lists the Voronoi Z-grid schemes that meet the majority of the challenges.

Solving Poisson’s equations at every time step is computationally more expensive when using a Z-grid scheme.Konor and Randall (2018) found that this would require 20% more total integration time with the multigrid technique.For future exascale computing, this may become less of a concern if improvements to the Z-grid scheme conformed to the numerical accuracy-related challenges associated withChallenges 2.1 and 2.5.

Xie (2019) proposed a generalization of the Voronoi Z-grid model by defining the state variables at the centroid centers.This generalization improved the Z-grid model and allowed it to overcomeChallenge 2.1on a sphere if they can demonstrate numerical stability.

Among all these schemes, the C- and Z-grid are the two schemes that may be preferred in the future development of a global model based on Table 2.Comparing the two schemes, Z-grid automatically maintains stationary geostrophic mode (Randall, 1994), while the C-grid Coriolis terms may cause unphysical behavior in its numerical solutions.Even though the techniques of adding linear combination weights had been introduced to ensure steady geostrophic modes and energy neutrality(Thuburn et al., 2009) to patch the C-grid, the drawback is to sacrifice the local precision on the cell edges, which may lead to the non-convergence of the Coriolis terms.The most recent comparison between A- and C-grid exposes a serious issue with the C-grid scheme (Yu et al.,2020), which seems difficult to overcome.In summary, a Z-grid scheme may have more flexibility to improve when compared to a C-grid scheme.Overcoming numerical accuracy-relatedChallenges 2.1 and 2.5should be the major research focus of future Z-grid research.

5.Conclusions

Previous studies have made significant progress toward improving FVM schemes to meet the accuracy requirements as we discussed in Section 2 by modifying grid structures.However, the modified grids cannot meet all these challenges (i.e.,2.1-2.5).Both CVT and spring dynamics addressChallenge 2.1 but not 2.2 and 2.5,whereas the tweaking algorithm focuses onChallenge 2.2 but not 2.1 and 2.5.Among all these schemes, the Zgrid scheme is the only one that meets all three physicalChallenges 3.1-3.3( Randall, 1994; Salmon, 2007;Eldred and Randall, 2017; Konor and Randall, 2018) if a Hamiltonian form of Z-grid scheme is used forChallenge 3.3based on the existing analysis.The A-grid scheme meetsChallenge 3.2only and C-grid meetsChallenge 3.3only.

Besides modifying the grid structure on a sphere, it is still possible to construct new finite-volume schemes that satisfy all the conditions and requirements discussed in this review.By reviewing the state-of-the-art global prediction models, we conclude that the Z-grid schemes may be the best candidates for improvements in numerical accuracy in future global weather forecasting models.

Acknowledgments.The authors would like to thank Dr.Ning Wang from the Earth System Research Laboratory, NOAA, for several in-depth discussions on general numerical modeling.The authors also appreciated the great effort of the anonymous reviewers for their comments and suggestions to make this paper in better quality.

Appendix: Impact of Inaccurate Discretization

An inaccurate finite-volume scheme does not provide numerical solutions that converge to the solution of the continuous governing equations.An experiment is provided here for readers to understand the importance of accuracy and the fact that we cannot alleviate numerical error solely by increasing the resolution.

Assessing a nonlinear numerical prediction model is difficult and involves various computational issues.Here,we use a simple governing equation to illustrate the impact of an inaccurate finite-volume scheme.A more complex system will undoubtedly face a worse impact.

The most straightforward conservative system is linear, and has the following form:

with a periodic boundary condition to avoid any potential problems associated with the boundary and a function of sin(x?t) as an initial condition over a domain(0,2π)×(0,8π) of (x,t).We investigate the forecast errors at timet= 8π using different discretization schemes.First, we test three time-integration schemes to avoid any issues with respect to temporal integration: leapfrog, and second- ( Süli and Mayers, 2003) and third-order Runge-Kutta schemes (Wicker and Skamarock, 2002).The solutions are shown in Fig.A1b, and they get the almost identical solution.

Equation (A2) assumes that uniform grid cells are used with a grid spacing of Δx.

Two zeroth-order-accurate schemes are introduced to simulate the loss of accuracy due to grid-cell irregularity on a sphere.The first scheme is conservative, but has the zeroth-order accuracy at each grid point and is as follows:

where ? is a small constant.We refer to this scheme as scheme A.The second scheme that can emulate a portion of the grid cells with distortion involves the introduction of a small constant, ?, at grid pointi:

is shown in Fig.A1a with different σ2values, i.e., 0.001 and 0.0001.Using Taylor expansion, we obtain:

Equation (A6) is also a zeroth-order accuracy scheme,which we refer to scheme B, in which we emulate distortion in a large portion of cells and to scheme C in which we emulate distortion in a small portion of cells.

Fig.A1.(a) The distribution of cell distortion ? (x) in the x-direction with two different values of σ2 for scheme B.(b) Time integrations, Heun second-order, Runge-Kutta, Wick-Skamarock (WS) third-order Runge-Kutta, and leapfrog all yield the same solution at a resolution of 161 grid points.A second-order scheme with different time integration schemes yields identical results with the maximum value of 1.0 while the zerothorder ones have much reduced values.(c) A comparison of the second-order scheme with the zeroth-order scheme A.The conserving schemes with zeroth-order schemes have significant phase errors.(d) A different number of distorted cells in schemes B (blue) and C (yellow and red) result in deviate solutions.Smaller distortion portions yield solutions that are closer to the true solution.(e) Comparison of scheme C with different signs of ?, and the true solution.(f) The numerical solutions for the second-order accurate and zeroth-order schemes under different resolutions.Inaccurate schemes do not converge to the true solution with increases in resolution.

To avoid any impact from temporal integration, we use the time integration schemes mentioned above to integrate the accurate Eq.(A2) and inaccurate A and B schemes to solve Eq.(A1).Figure A1b shows the second- and zeroth-order schemes (B) with three different temporal integration schemes, where ? (x) ≡awith a value of 0.1.The second-order scheme solutions converge to the true solution, whereas the zeroth-order scheme solutions all converge to a state different from the true state.Large errors associated with the schemes derive solely from the inaccuracy of the zeroth-order schemes and not from the temporal integration schemes.

Figure A1c shows the results of the second-order and zeroth-order schemes.Scheme A conserves energy but transports the solutions faster than the true solution.This indicates that energy conservation is meaningless if a scheme is not accurate.Figure A1d shows the second-order scheme and zeroth-order scheme B with differentσ values used in Eq.(A5).The results show that an increase in the number of distorted grid cells yields poorer numerical solutions.

We also plotted the numerical solutions for scheme C with a different sign for ? in Fig.A1e.The solution where ? > 0 is known as an expanding cell and the solution where ? < 0 is a shrinking cell.None of these solutions tend to be close to the true solution.

Finally, to show that high resolutions do not improve the solution, Fig.A1f illustrates the second-order and zeroth-order accurate scheme B at three resolutions.Higher resolutions cannot help an inaccurate scheme.

In summary, an inaccurate finite-volume scheme results in significant errors.Higher resolutions do not yield numerical solutions that are closer to the true solution.Model accuracy is critical.