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

构建具有数值计算精度的全球有限体积天气气候模式面临的挑战

+ Author Affiliations + Find other works by these authors
  • Corresponding author: Yuanfu XIE, xieyf@cma.gov.cn
  • Funds:

    Supported by the National Key Research and Development Program of China (2017YFC1502201) and Basic Scientific Research and Operation Fund of Chinese Academy of Meteorological Sciences (2017Z017)

  • doi: 10.1007/s13351-021-0202-3

PDF

  • 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.
    近年来,随着对更精细化天气预报和更精准气候预测的需求的不断增加,高分辨率全球网格非静力平衡预报模式的构建一直是大气科学领域研究与应用的热点。目前构建满足各项数学、物理及计算要求的网格模式仍然遇到严峻的挑战。本综述深入讨论了在构造基于有限体积法的全球数值模式动力框架中面临的主要挑战,即如何保证球面上各种算子的离散数值精度、模式的频散关系、预测变量的守恒特性以及计算效率问题。其中,如何在非均匀球面网格上评估以及改进数值精度是本文讨论的重点。文章回顾了目前已经广泛应用的基于有限体积法的全球数值模式动力框架,逐个分析了其数值精度、频散以及守恒问题,并对各个框架在各项指标上的优缺点进行评估,详尽展示基于有限体积法的全球动力模式的发展状况。尽管文中的分析和讨论主要基于浅水方程系统,但所涉及的大多数概念仍可广泛地用于评估其他更复杂的动力框架模式。优化并解决这些挑战与难点,是未来高精度全球动力模式研究的主要方向。
  • 加载中
  • Fig. 1.  Two quasi-uniform grids on a sphere, with the icosahedral grid on the left and the cubed sphere on the right.

    Fig. 2.  Finite-volume scheme definitions for the A-, B-, C-, C–D, D-, and Z-grid schemes, where ${h}$ is the fluid thickness, ${\boldsymbol{V}}$ is the velocity vector, ${{u}}_{\rm{n}}$ is the normal velocity, ${{u}}_{{\tau }}$ is the tangential velocity component, and $ \mathrm{\zeta } $ and $ \mathrm{\delta } $ are the vorticity and divergence, respectively.

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

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

    Fig. 5.  Definition of $ {\lambda }^{\perp } $ of the distance between an arc-center and its edge intersection.

    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.

    Fig. A1.  (a) The distribution of cell distortion $ \epsilon \left(x\right) $ in the x-direction with two different values of $ {\sigma }^{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 zeroth-order 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 $ \epsilon $, 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.

    Table 1.  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

    $ l $0.50.250.1250.06250.03125
    M10.021657070.005432330.001359240.000339910.00008502
    M20.017484290.034351910.038586720.039646540.03991157
    M30.061043280.045278270.041320660.040330240.04008257
    Download: Download as CSV

    Table 2.  Matrix of all the discussed schemes

    Challenge2.12.22.32.42.53.13.23.3
    A-gridYesNoNoYesNoNoYesNo
    C-CVTYesNoYesYesNoNoNoYes
    C-circumNoYesNoYesNoNoNoNo
    Z-gridNoYesYesYesNoYesYesYes
    Note that no staggered schemes overcome all challenges now, but the Z-grid scheme meets more challenges than the others. We make a list in Table 3 for readers to recall the challenges.
    Download: Download as CSV

    Table 3.  Challenge list

    ChallengeContent
    2.1Choice of cell center to pose the state variable
    2.2Center-arc may not bisect the edge for irregular grid
    2.3Edge may not bisect the center-arc for irregular grid
    2.4Center-arc may not be perpendicular to the edge for irregular grid
    2.5Construct the tangential derivatives with high-order accuracy
    3.1Preserving the dispersion relation independent of horizontal grid size
    3.2Avoid computation mode
    3.3Energy and enstrophy conservation for shallow-water equations
    Download: Download as CSV
  • [1]

    Adcroft, A. J., C. N. Hill, and J. C. Marshall, 1999: A new treatment of the Coriolis terms in C-grid models at both high and low resolutions. Mon. Wea. Rev., 127, 1928–1936. doi: 10.1175/1520-0493(1999)127<1928:ANTOTC>2.0.CO;2.
    [2]

    Bao, L., R. D. Nair, and H. M. Tufo, 2014: A mass and momentum flux-form high-order discontinuous Galerkin shallow water model on the cubed-sphere. J. Comput. Phys., 271, 224–243. doi: 10.1016/j.jcp.2013.11.033.
    [3]

    Du, Q., M. D. Gunzburger, and L. L. Ju, 2003: Constrained centroidal Voronoi tessellations for surfaces. SIAM J. Sci. Comput., 24, 1488–1506. doi: 10.1137/S1064827501391576.
    [4]

    Eldred, C., and D. Randall, 2017: Total energy and potential enstrophy conserving schemes for the shallow water equations using Hamiltonian methods—Part 1: Derivation and properties. Geosci. Model Dev., 10, 791–810. doi: 10.5194/gmd-10-791-2017.
    [5]

    Gill, A. E., 1982: Atmosphere–Ocean Dynamics. Academic Press, San Diego, 662 pp.
    [6]

    Giraldo, F. X., J. S. Hesthaven, and T. Warburton, 2002: Nodal high-order discontinuous Galerkin methods for the spherical shallow water equations. J. Comput. Phys., 181, 499–525. doi: 10.1006/jcph.2002.7139.
    [7]

    Heikes, R., and D. A. Randall, 1995a: Numerical integration of the shallow-water equations on a twisted icosahedral grid. Part I: Basic design and results of tests. Mon. Wea. Rev., 123, 1862–1880. doi: 10.1175/1520-0493(1995)123<1862:NIOTSW>2.0.CO;2.
    [8]

    Heikes, R., and D. A. Randall, 1995b: Numerical integration of the shallow-water equations on a twisted icosahedral grid. Part II: A detailed description of the grid and an analysis of numerical accuracy. Mon. Wea. Rev., 123, 1881–1887. doi: 10.1175/1520-0493(1995)123<1881:NIOTSW>2.0.CO;2.
    [9]

    Heikes, R. P., D. A. Randall, and C. S. Konor, 2013: Optimized icosahedral grids: Performance of finite-difference operators and multigrid solver. Mon. Wea. Rev., 141, 4450–4469. doi: 10.1175/MWR-D-12-00236.1.
    [10]

    Kanamitsu, M., J. C. Alpert, K. A. Campana, et al., 1991: Recent changes implemented into the global forecast system at NMC. Wea. Forecasting, 6, 425–435. doi: 10.1175/1520-0434(1991)006<0425:RCIITG>2.0.CO;2.
    [11]

    Konor, C. S., and D. A. Randall, 2018: Impacts of the horizontal and vertical grids on the numerical solutions of the dynamical equations—Part 1: Nonhydrostatic inertia–gravity modes. Geosci. Model Dev., 11, 1753–1784. doi: 10.5194/gmd-11-1753-2018.
    [12]

    Lax, P. D., and R. D. Richtmyer, 1956: Survey of the stability of linear finite difference equations. Commun. Pure Appl. Math., 9, 267–293. doi: 10.1002/cpa.3160090206.
    [13]

    Lin, S.-J., 2004: A “vertically Lagrangian” finite-volume dynamical core for global models. Mon. Wea. Rev., 132, 2293–2307. doi: 10.1175/1520-0493(2004)132<2293:AVLFDC>2.0.CO;2.
    [14]

    Lin, S.-J., and R. B. Rood, 1996: Multidimensional flux-form semi-Lagrangian transport schemes. Mon. Wea. Rev., 124, 2046–2070. doi: 10.1175/1520-0493(1996)124<2046:MFFSLT>2.0.CO;2.
    [15]

    Putman, W. M., and S.-J. Lin, 2007: Finite-volume transport on various cubed-sphere grids. J. Comput. Phys., 227, 55–78. doi: 10.1016/j.jcp.2007.07.022.
    [16]

    Rajpoot, M. K., S. Bhaumik, and T. K. Sengupta, 2012: Solution of linearized rotating shallow water equations by compact schemes with different grid-staggering strategies. J. Comput. Phys., 231, 2300–2327. doi: 10.1016/j.jcp.2011.11.025.
    [17]

    Randall, D. A., 1994: Geostrophic adjustment and the finite-difference shallow-water equations. Mon. Wea. Rev., 122, 1371–1377. doi: 10.1175/1520-0493(1994)122<1371:GAATFD>2.0.CO;2.
    [18]

    Reinecke, P. A., and D. Durran, 2009: The overamplification of gravity waves in numerical solutions to flow over topography. Mon. Wea. Rev., 137, 1533–1549. doi: 10.1175/2008MWR2630.1.
    [19]

    Ringler, T. D., J. Thuburn, J. B. Klemp, et al., 2010: A unified approach to energy conservation and potential vorticity dynamics for arbitrarily-structured C-grids. J. Comput. Phys., 229, 3065–3090. doi: 10.1016/j.jcp.2009.12.007.
    [20]

    Sadourny, R., A. Arakawa, and Y. Mintz, 1968: Integration of the nondivergent barotropic vorticity equation with an icosahedral-hexagonal grid for the sphere. Mon. Wea. Rev., 96, 351–356. doi: 10.1175/1520-0493(1968)096<0351:IOTNBV>2.0.CO;2.
    [21]

    Saito, K., J.-I. Ishida, K. Aranami, et al., 2007: Nonhydrostatic atmospheric models and operational development at JMA. J. Meteor. Soc. Japan, 85B, 271–304. doi: 10.2151/jmsj.85B.271.
    [22]

    Salmon, R., 1998: Lectures on Geophysical Fluid Dynamics. Oxford University Press, New York, 378 pp.
    [23]

    Salmon, R., 2007: A general method for conserving energy and potential enstrophy in shallow-water models. J. Atmos. Sci., 64, 515–531. doi: 10.1175/JAS3837.1.
    [24]

    Satoh, M., T. Matsuno, H. Tomita, et al., 2008: Nonhydrostatic icosahedral atmospheric model (NICAM) for global cloud resolving simulations. J. Comput. Phys., 227, 3486–3514. doi: 10.1016/j.jcp.2007.02.006.
    [25]

    Simmons, A. J., D. M. Burridge, M. Jarraud, et al., 1989: The ECMWF medium-range prediction models development of the numerical formulations and the impact of increased resolution. Meteor. Atmos. Phys., 40, 28–60. doi: 10.1007/BF01027467.
    [26]

    Skamarock, W. C., J. B. Klemp, M. G. Duda, et al., 2012: A multiscale nonhydrostatic atmospheric model using centroidal Voronoi tesselations and C-grid staggering. Mon. Wea. Rev., 140, 3090–3105. doi: 10.1175/MWR-D-11-00215.1.
    [27]

    Süli, E., and D. F. Mayers, 2003: An Introduction to Numerical Analysis. Cambridge University Press, Cambridge, 444 pp.
    [28]

    Thuburn, J., T. D. Ringler, W. C. Skamarock, et al., 2009: Numerical representation of geostrophic modes on arbitrarily structured C-grids. J. Comput. Phys., 228, 8321–8335. doi: 10.1016/j.jcp.2009.08.006.
    [29]

    Tomita, H., M. Tsugawa, M. Satoh, et al., 2001: Shallow water model on a modified icosahedral geodesic grid by using spring dynamics. J. Comput. Phys., 174, 579–613. doi: 10.1006/jcph.2001.6897.
    [30]

    Tomita, H., M. Satoh, and K. Goto, 2002: An optimization of the icosahedral grid modified by spring dynamics. J. Comput. Phys., 183, 307–331. doi: 10.1006/jcph.2002.7193.
    [31]

    Ullrich, P. A., C. Jablonowski, J. Kent, et al., 2017: DCMIP2016: a review of non-hydrostatic dynamical core design and intercomparison of participating models. Geosci. Model Dev., 10, 4477–4509. doi: 10.5194/gmd-10-4477-2017.
    [32]

    Wan, H., M. A. Giorgetta, G. Zängl, et al., 2013: The ICON-1.2 hydrostatic atmospheric dynamical core on triangular grids—Part 1: Formulation and performance of the baseline version. Geosci. Model Dev., 6, 735–763. doi: 10.5194/gmd-6-735-2013.
    [33]

    Wicker, L. J., and W. C. Skamarock, 2002: Time-splitting methods for elastic models using forward time schemes. Mon. Wea. Rev., 130, 2088–2097. doi: 10.1175/1520-0493(2002)130<2088:TSMFEM>2.0.CO;2.
    [34]

    Williamson, D. L., and P. J. Rasch, 1994: Water vapor transport in the NCAR CCM2. Tellus A Dyn. Meteor. Oceanogr., 46, 34–51. doi: 10.3402/tellusa.v46i1.15426.
    [35]

    Xie, Y. F., 2019: Generalized Z-grid model for numerical weather prediction. Atmosphere, 10, 179. doi: 10.3390/atmos10040179.
    [36]

    Yu, Y. G., N. Wang, J. Middlecoff, et al., 2020: Comparing numerical accuracy of icosahedral A-grid and C-grid schemes in solving the shallow-water model. Mon. Wea. Rev., 148, 4009–4033. doi: 10.1175/MWR-D-20-0024.1.
  • 加载中
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

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

    Corresponding author: Yuanfu XIE, xieyf@cma.gov.cn
  • 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
Funds: Supported by the National Key Research and Development Program of China (2017YFC1502201) and Basic Scientific Research and Operation Fund of Chinese Academy of Meteorological Sciences (2017Z017)

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.

构建具有数值计算精度的全球有限体积天气气候模式面临的挑战

近年来,随着对更精细化天气预报和更精准气候预测的需求的不断增加,高分辨率全球网格非静力平衡预报模式的构建一直是大气科学领域研究与应用的热点。目前构建满足各项数学、物理及计算要求的网格模式仍然遇到严峻的挑战。本综述深入讨论了在构造基于有限体积法的全球数值模式动力框架中面临的主要挑战,即如何保证球面上各种算子的离散数值精度、模式的频散关系、预测变量的守恒特性以及计算效率问题。其中,如何在非均匀球面网格上评估以及改进数值精度是本文讨论的重点。文章回顾了目前已经广泛应用的基于有限体积法的全球数值模式动力框架,逐个分析了其数值精度、频散以及守恒问题,并对各个框架在各项指标上的优缺点进行评估,详尽展示基于有限体积法的全球动力模式的发展状况。尽管文中的分析和讨论主要基于浅水方程系统,但所涉及的大多数概念仍可广泛地用于评估其他更复杂的动力框架模式。优化并解决这些挑战与难点,是未来高精度全球动力模式研究的主要方向。
    • 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 traditio-nal 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 dynami-cal 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 quasi-uniform 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.

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

      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.

      Figure 2.  Finite-volume scheme definitions for the A-, B-, C-, C–D, D-, and Z-grid schemes, where ${h}$ is the fluid thickness, ${\boldsymbol{V}}$ is the velocity vector, ${{u}}_{\rm{n}}$ is the normal velocity, ${{u}}_{{\tau }}$ is the tangential velocity component, and $ \mathrm{\zeta } $ and $ \mathrm{\delta } $ are the vorticity and divergence, respectively.

      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.

    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.

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

      $$ \frac{\partial \boldsymbol{V}}{\partial t}+\eta \boldsymbol{k}\times \boldsymbol{V}+\nabla \left[K+{\rm{g}}\left(h+{h}_{{\rm{s}}}\right)\right]=0, $$ (1)
      $$ \frac{\partial h}{\partial t}+\nabla \cdot \left(h\boldsymbol{V}\right)=0, $$ (2)

      where $ h $ is the fluid thickness, $\boldsymbol{V}$ is the velocity vector, $\eta =\boldsymbol{k}\cdot \nabla \times \boldsymbol{V}+f$ is the absolute vorticity, $ f $ is the Coriolis parameter, $\boldsymbol{k}$ is the vertical unit vector, $K=\dfrac{1}{2}\boldsymbol{V}\cdot \boldsymbol{V}$ is the kinetic energy, ${\rm{g}}$ is gravitational acceleration, and ${h}_{{\rm{s}}}$ is the surface topography. These equations can be rewritten in their vorticity and divergence forms (Heikes and Randall, 1995a) as:

      $$ \frac{\partial \eta }{\partial t}+\nabla \cdot \left(\eta \nabla \chi \right)-J\left(\eta,\psi \right)=0, $$ (3)
      $$ \frac{\partial \delta }{\partial t}-\nabla \cdot \left(\eta \nabla \psi \right)-J\left(\eta,\chi \right)+{\nabla }^{2}\left[K+{\rm{g}}\left(h+{h}_{{\rm{s}}}\right)\right]=0, $$ (4)
      $$ \frac{\partial h}{\partial t}+\nabla \cdot \left(h\nabla \chi \right)-J\left(h,\psi \right)=0, $$ (5)

      where $ \eta $ is the absolute vorticity, $ \delta $ is the horizontal divergence, and J is the Jacobi operator [$J\left(\eta,\psi \right)$ = ${\boldsymbol{ k}} \cdot ( \nabla \eta \times \nabla \psi) $], with their corresponding streamfunction $ \psi $ and velocity potential $\; \chi $ satisfying the following relations:

      $$ {\nabla }^{2}\psi =\eta -f, $$ (6)
      $$\hspace{-15pt} {\nabla }^{2}\chi =\delta . $$ (7)

      The kinetic energy has the form of a streamfunction and velocity potential as:

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

      To solve Eqs. (1)–(2) or (3)–(7), finite-volume schemes integrate Eqs. (1)–(2) or (3)–(5) over each cell and yield equations for the mean values,

      $$ \bar{\alpha }=\frac{1}{A}{\int }_{A}^{}\alpha {\rm{d}}s, $$ (9)

      where $ \bar{\alpha } $ is an average value of $ \alpha $ over the cell A. Based on the Gauss theorem, these integrals can be converted into line integrals along edges surrounding the cell and further written into the linear combination of function values at the stencil. For example, the divergence operator

      $$ \frac{1}{A}{\int }_{A}^{}\nabla \cdot \left(h\boldsymbol{V}\right){\rm{d}}s=\frac{1}{A}{\oint}_{\partial A}^{}\boldsymbol{n}\cdot \left(h\boldsymbol{V}\right){\rm{d}}l, $$ (10)

      or the Laplacian operator

      $$ \frac{1}{A}{\int }_{A}^{}{\nabla }^{2}\alpha {\rm{d}}s=\frac{1}{A}{\oint }_{\partial A}^{}\frac{\partial \alpha }{\partial n}{\rm{d}}l, $$ (11)

      where $ A $ is the area of an underlying cell, $ \partial A $ is the edge that surrounds $ A $, and the normal unit vector, $\boldsymbol{n}$, points outwards from the cell from a point on the edge center.

    • 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):

      $$ \frac{1}{A}{\oint }_{\partial A}^{}\boldsymbol{n}\cdot \left(h\boldsymbol{V}\right){\rm{d}}l \approx \frac{1}{A}\sum _{i}{l}_{i}{\left.\left[\boldsymbol{n}\cdot \left(h\boldsymbol{V}\right)\right]\right|}_{c}, $$ (12)

      or

      $$ \frac{1}{A}{\oint }_{\partial A}^{}\frac{\partial \alpha }{\partial n}{\rm{d}}l \approx \frac{1}{A}\sum _{i}{l}_{i}{\left.\frac{\partial \alpha }{\partial n}\right|}_{c}, $$ (13)

      where $ {l}_{i} $ is the length of the $ i $th edge of the cell, and the subscript $ c $ denotes the edge center. For any regular polygon with more than three vertices (e.g., a square cell), if a discretization scheme achieves the approximation of the integrand for either ${\left.\boldsymbol{n}\cdot \left(h\boldsymbol{V}\right)\right|}_{c}$ or $ {\left.\dfrac{\partial \alpha }{\partial n}\right|}_{c} $ to second-order accuracy, it also results in a second-order accurate approximation of $\dfrac{1}{A}\sum _{i}{l}_{i}{\left.\left[\boldsymbol{n}\cdot \left(h\boldsymbol{V}\right)\right]\right|}_{c}$ or $\dfrac{1}{A} \sum _{i}{l}_{i} $${\left.\dfrac{\partial \alpha }{\partial n}\right|}_{c} $. That is, the order of accuracy is preserved.

      It is important to note, however, that this preservation of accuracy does not hold for irregular cells. A second-order accurate approximation of either ${\left.\boldsymbol{n}\cdot \left(h\boldsymbol{V}\right)\right|}_{c}$ or $ {\left.\dfrac{\partial \alpha }{\partial n}\right|}_{c} $ can only result in a first-order accurate approximation of either $\dfrac{1}{A}\sum _{i}{l}_{i}{\left.\left[\boldsymbol{n}\cdot \left(h\boldsymbol{V}\right)\right]\right|}_{c}$ or $ \dfrac{1}{A}\sum _{i}{l}_{i} $${\left.\dfrac{\partial \alpha }{\partial n}\right|}_{c} $, at most in general [see the analysis in Heikes and Randall (1995b)]. We call this phenomenon “accuracy drop” due to the irregularity of cells. Moreover, for any irregular polygon, a first-order approximation of these integrands will result in a zeroth-order accuracy, which leads to a non-convergent scheme.

      To illustrate this accuracy drop, we performed a simple numerical experiment on Eq. (11) to demonstrate the accuracy losses due to the approximations of the above integrands. We prepared three approximation schemes: one second-order scheme and two first-order schemes. The second-order scheme is based on a regular hexagon cell, while the two first-order schemes are formed by moving the finite-difference point $ \dfrac{\partial \alpha }{\partial n} $ from the edge center to a distance of $ 0.1\times l $ in either the normal or tangential direction, to simulate irregular cells. These schemes are labeled as M1 (regular), M2 (irregular, variation in the normal direction), and M3 (irregular, variation in the tangential direction). The function is assumed to be $ \alpha =\mathrm{cos}\left(x+y\right) $. In all three schemes, $ \dfrac{\partial \alpha }{\partial n} $ adopts the value of the continuous function values at its location. Table 1 lists the numerical errors for $ \dfrac{1}{A}\sum _{i}{l}_{i} $${\left.\dfrac{\partial \alpha }{\partial n}\right|}_{c} $ by approximating ${\nabla }^{2}\alpha =\dfrac{1}{A}{\oint }_{\partial A}^{}\dfrac{\partial \alpha }{\partial n}{\rm{d}}l$ on these hexagonal cells.

      $ l $0.50.250.1250.06250.03125
      M10.021657070.005432330.001359240.000339910.00008502
      M20.017484290.034351910.038586720.039646540.03991157
      M30.061043280.045278270.041320660.040330240.04008257

      Table 1.  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

      The errors in Table 1 clearly show that the M1 scheme (second-order) reduces the error by a factor of four whenever the edge length is reduced by a half, indicating that the approximation is also second-order accurate on this regular hexagon cell. However, both M2 and M3 (first-order) do not yield any improvement in accuracy. They are zeroth-order accurate approximations of $ {\nabla }^{2}\alpha $. This fact also holds for other operators, such as Eq. (10).

      This is a crucial point regarding FVM schemes constructed on a sphere: one must construct a higher accuracy approximation of the line integrand in order to achieve accuracy of one-order lower in the approximation of the integrals. For example, if even a first-order scheme to approximate $ {\nabla }^{2}\alpha =\dfrac{1}{A}{\oint }_{\partial A}^{} $$\dfrac{\partial \alpha }{\partial n}{\rm{d}}l$ is required, one has to construct at least a second-order accurate approximation to $ \dfrac{\partial \alpha }{\partial n} $. With our simple cases, we can reveal several challenges that must be addressed in order to achieve accuracy for FVM on the irregular grids.

    • 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., $ \bar{\alpha } $ 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 to Challenge 2.1 for an FVM scheme on a sphere (note that all of these centers coincide at the same point if the underlying polygon is regular).

      Figure 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,

      $$ \bar{\alpha }=\frac{1}{A}{\iint }_{}^{}\alpha {\rm{d}}s+O\left({l}^{2}\right), $$ (14)

      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, $ \epsilon $, multiplied by the edge length, $ l $, for an irregular cell (e.g., some cells inside either an icosahedral or cubed sphere grid):

      $$ \frac{d}{l}\geqslant \epsilon >0. $$ (16)

      So, the ratio of d and l measures 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 value $ \bar{\alpha } $ over 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 as Challenge 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 center-arc, 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 $ \lambda $: 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 $ \lambda $ and the edge length, $ l $, has a positive lower bound approximately equal to 0.088 (Heikes et al., 2013):

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

      $$ \frac{\lambda }{l}\geqslant \varepsilon >0. $$ (17)

      The constant, $ \varepsilon $, is greater than 3.6 for a grid modified by spring dynamics (Tomita et al., 2001). Although the tweaking optimization algorithm (Heikes and Randall, 1995b; Heikes et al., 2013) minimizes this ratio, $ \varepsilon $ still shows no sign of convergence to zero after the grid level G13, as listed in their Table 1 of Heikes et al. (2013).

      If one simply uses a finite-difference of the function values at their two centers to approximate the $ \dfrac{\partial \alpha }{\partial n} $, the $ {\left.\dfrac{\partial \alpha }{\partial n}\right|}_{c} $ will be only first-order accurate, because of Eq. (16). Moreover, the accuracy drop makes this scheme a zeroth-order accuracy approximation.

      Challenge 2.3: An edge may not bisect the center-arc. How to construct an accurate FVM?

      Edges only bisect their Voronoi center-arcs. We similarly denote the distance of intersecting point and middle of a center-arc by $ {\lambda }^{\perp } $ in Fig. 5. Obviously, on a sphere, a positive constant may exist, such that

      Figure 5.  Definition of $ {\lambda }^{\perp } $ of the distance between an arc-center and its edge intersection.

      $$ \frac{{\lambda }^{\perp }}{l}\geqslant \epsilon >0, $$ (18)

      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 center-arcs 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.

      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., $\boldsymbol{k}\times \boldsymbol{V}$ for Eq. (1) or $ J $ for 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.

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

      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.

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

      $$\hspace{-35pt} \frac{\partial u}{\partial t}-fv+{\rm{g}}\frac{\partial h}{\partial x}=0, $$ (19)
      $$\hspace{-35pt} \frac{\partial v}{\partial t}+fu+{\rm{g}}\frac{\partial h}{\partial y}=0, $$ (20)
      $$\hspace{-25pt} \frac{\partial h}{\partial t}+{\rm{H}}\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)=0. $$ (21)

      This results in an equation for the divergence:

      $$\hspace{-35pt} \frac{\partial \delta }{\partial t}-f\zeta +{\rm{g}}{\nabla }^{2}h=0, $$ (22)

      or its equivalent in terms of thickness $ h $:

      $$ -\frac{{\partial }^{3}h}{\partial {t}^{3}}-{f}^{2}\frac{\partial h}{\partial t}+{\rm{g}}{\rm{H}}\frac{\partial {\nabla }^{2}h}{\partial t}=0. $$ (23)

      Assuming a plane wave solution $h={h}_{0}{{\rm{e}}}^{{\rm{i}}\left({k}_{x}x+{k}_{y}y-\omega t\right)}$, where $ {k}_{x} $ and $ {k}_{y} $ are the horizontal wavenumber components and H is the constant depth of the “water,” by substituting the solution into Eq. (22), we arrive at the dispersion relation for the continuous system:

      $$ \omega \left[{\left(\frac{\omega }{f}\right)}^{2}-1-{\gamma }^{2}\left({{k}_{x}^{2}}+{{k}_{y}^{2}}\right)\right]=0, $$ (24)

      where

      $$ \gamma =\frac{\sqrt{{\rm{g}}{\rm{H}}}}{f}, $$ (25)

      is the deformation radius (Gill, 1982), and $ f $ is 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 $ \gamma /d $ ratio, where $ d $ is the grid size or distance between centers. Dispersion relations become unsatisfactory if $ \gamma /d $ is too small for certain schemes (Randall, 1994).

      For modern numerical prediction models, grid sizes of $ d $ are small enough to resolve the deformation radius, i.e., the $ \gamma /d $ ratio is sufficiently large if $ \gamma $ 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 to $ {{k}_{z}^{-1}} $, where $ {k}_{z} $ is the vertical wavenumber (Gill, 1982; Salmon, 1998). This indicates that the $ \gamma /d $ ratio 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.

      Challenge2.12.22.32.42.53.13.23.3
      A-gridYesNoNoYesNoNoYesNo
      C-CVTYesNoYesYesNoNoNoYes
      C-circumNoYesNoYesNoNoNoNo
      Z-gridNoYesYesYesNoYesYesYes
      Note that no staggered schemes overcome all challenges now, but the Z-grid scheme meets more challenges than the others. We make a list in Table 3 for readers to recall the challenges.

      Table 2.  Matrix of all the discussed schemes

      ChallengeContent
      2.1Choice of cell center to pose the state variable
      2.2Center-arc may not bisect the edge for irregular grid
      2.3Edge may not bisect the center-arc for irregular grid
      2.4Center-arc may not be perpendicular to the edge for irregular grid
      2.5Construct the tangential derivatives with high-order accuracy
      3.1Preserving the dispersion relation independent of horizontal grid size
      3.2Avoid computation mode
      3.3Energy and enstrophy conservation for shallow-water equations

      Table 3.  Challenge list

    • 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 overcome Challenge 2.1 by 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 passes Challenge 2.4. The center-arc and edge are not perpendicular to each other. Challenges 2.2, 2.3, and 2.5 are problematic for an A-grid scheme. The center-arc and edge do not bisect each other and cannot overcome Challenges 2.2 and 2.3, as listed in Table 2 of Heikes et al. (2013); i.e., there are positive constants of $ \epsilon $, 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 first-order accurate. Thus, an A-grid scheme fails Challenges 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 pass Challenge 3.2. There are no previous studies that demonstrate A-grid’s energy and enstrophy conservation for shallow-water equations. The A-grid row of Table 2 lists how well an A-grid scheme overcomes these challenges.

    • A C-grid scheme places its mass field at the cell centers and the normal velocity field on the edges. Two C-grid 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 fails Challenge 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 fail Challenge 2.2, while circumcenter schemes meet it. Conversely, CVT schemes meet Challenge 2.3, but circumcenter schemes fail. Both of them do well for Challenge 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.5 remains 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 fail Challenge 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 fails Challenge 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.

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

    • 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 meet Challenge 2.1, which introduces the first-order error. Heikes and Randall (1995b) developed a tweaking algorithm that modifies the underlying grid to overcome Challenge 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 meets Challenges 2.3 and 2.4. However, the major issue with a Voronoi Z-grid scheme is its ability to overcome Challenge 2.5, as explained in the following paragraph.

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

      $$ {\int }_{c}^{}J\left(\sigma,\beta \right){\rm{d}}a={\oint }_{\partial c}^{}\alpha \frac{\partial \beta }{\partial \tau }{\rm{d}}l. $$ (26)

      To approximate the tangential derivative, a Voronoi Z-grid 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 first-order 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 fails Challenge 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 meets Challenge 3.1. Since it has the same number of discretized equations for all the state variables, Z-grid models overcome Challenge 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 with Challenges 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 overcome Challenge 2.1 on 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-related Challenges 2.1 and 2.5 should be the major research focus of future Z-grid research.

    5.   Conclusions
    • By reviewing several of the primary challenges that plague the development of a global numerical prediction model using a finite-volume method, we are able to clearly understand the major issues that we face when attempting to construct a numerically and physically accurate prediction model on the global sphere. Such models have to be numerically accurate, physically meaningful, and computationally efficient. However, there are no perfect FVM schemes yet that can overcome all the challenges discussed in this review.

      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.12.5). Both CVT and spring dynamics address Challenge 2.1 but not 2.2 and 2.5, whereas the tweaking algorithm focuses on Challenge 2.2 but not 2.1 and 2.5. Among all these schemes, the Z-grid scheme is the only one that meets all three physical Challenges 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 for Challenge 3.3 based on the existing analysis. The A-grid scheme meets Challenge 3.2 only and C-grid meets Challenge 3.3 only.

      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:

      $$ \frac{\partial u}{\partial t}+\frac{\partial u}{\partial x}=0, \tag{A1} $$

      with a periodic boundary condition to avoid any potential problems associated with the boundary and a function of $ \mathrm{sin}\left(x-t\right) $ as an initial condition over a domain $ \left(0,2\pi \right)\times \left(0,8\pi \right) $ of $ \left(x,t\right) $. We investigate the forecast errors at time $ t=8\pi $ 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.

      Figure A1.  (a) The distribution of cell distortion $ \epsilon \left(x\right) $ in the x-direction with two different values of $ {\sigma }^{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 zeroth-order 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 $ \epsilon $, 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.

      For the given spatial derivative, $ \dfrac{\partial u}{\partial x} $, a second-order accurate finite-volume scheme is as follows:

      $$ {\left.\frac{\partial u}{\partial x}\right|}_{i}^{n}=\frac{{u}_{i+1}^{n}-{u}_{i-1}^{n}}{2\Delta x}+O\left({\Delta x}^{2}\right). \tag{A2} $$

      Equation (A2) assumes that uniform grid cells are used with a grid spacing of $ \Delta 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:

      $$ \frac{{u}_{i+1}^{n}-{u}_{i-1}^{n}}{2\Delta x}\left(1+\epsilon \right)\approx {\left.\frac{\partial u}{\partial x}\right|}_{i}^{n}+O\left(\epsilon \right), \tag{A3} $$

      where $ \epsilon $ 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, $ \epsilon $, at grid point $ i $:

      $$ \frac{{\left[1+\epsilon \left({x}_{i}\right)\Delta x\right]u}_{i+1}^{n}-{u}_{i-1}^{n}}{2\Delta x}, \tag{A4} $$

      where

      $$ \epsilon \left(x\right)=a{{\rm{e}}}^{\frac{-{\left(x-0.5\right)}^{4}}{{\sigma }^{2}}}, \tag{A5} $$

      is shown in Fig. A1a with different $ {\sigma }^{2} $ values, i.e., 0.001 and 0.0001. Using Taylor expansion, we obtain:

      $$ \frac{{\left(1+\epsilon \Delta x\right)u}_{i+1}^{n}-{u}_{i-1}^{n}}{2\Delta x}={\left.\frac{\partial u}{\partial x}\right|}_{i}^{n}+O\left({\Delta x}^{2}\right)+\frac{\epsilon }{2}{u}_{i+1}^{n}={\left.\frac{\partial u}{\partial x}\right|}_{i}^{n}+O\left(\epsilon \right). \tag{A6} $$

      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.

      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 $ \epsilon \left(x\right)\equiv a $ with 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 $ \sigma $ 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 $ \epsilon $ in Fig. A1e. The solution where $ \epsilon >0 $ is known as an expanding cell and the solution where $ \epsilon < 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.

Reference (36)

Catalog

    /

    DownLoad:  Full-Size Img  PowerPoint
    Return
    Return