Configuration and evaluation of a global unstructured mesh atmospheric model (GRIST-A20.9) based on the variable-resolution approach
- 1State Key Laboratory of Numerical Modeling for Atmospheric Sciences and Geophysical Fluid Dynamics (LASG), Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing, China
- 2University of Chinese Academy of Sciences, Beijing, China
- 3State Key Laboratory of Severe Weather (LaSW), Chinese Academy of Meteorological Sciences, China Meteorological Administration, Beijing, China
- 4National Supercomputing Center, Wuxi, Jiangsu, China
Correspondence: Yi Zhang (email@example.com)
Targeting a long-term effort towards a variable-resolution (VR) global weather and climate model, this study systematically configures and evaluates an unstructured mesh atmospheric model based on the multiresolution approach. The model performance is examined from dry dynamics to simple physics and full physics scenarios. In the dry baroclinic wave test, the VR model reproduces comparable fine-scale structures in the refined regions as a fine-resolution quasi-uniform (QU) mesh model. The mesh transition zone does not adversely affect the wave pattern. Regional kinetic energy spectra show that the fine-scale resolving ability improves as the fine resolution increases. Compared to a QU counterpart that has equivalent degrees of freedom, the VR model tends to increase the global errors, but the errors can be reduced when the resolution of the coarse region is increased. The performance over the coarse region is generally close to that of a low-resolution QU counterpart. Two multi-region refinement approaches, the hierarchical and polycentric refinement modes, further validate the model performance under the multiresolution refinement. Activating hyperdiffusion for horizontal velocity is helpful with respect to VR modeling. An idealized tropical cyclone test is further used to examine its ability to resolve fine-scale structures. In the simple physics environment, the VR model can have the tropical cyclone stably pass the transition zone in various configurations. A series of sensitivity tests examines the model performance in a hierarchical refinement mode. The simulations exhibit consistency even when the VR mesh is slightly perturbed by one of the three parameters that control the density function. The tropical cyclone, starting from the second refinement region and passing through the inner transition zone, gets intensified and covers a smaller area in the refined regions. Such variations are consistent with the behavior that one may observe when uniformly refining the QU mesh. In the full physics environment with a highly variable mesh that reaches sub-10 km resolution, the VR model also produces a reasonable evolution for the tropical cyclone. The explicit diffusion shows its usefulness in terms of suppressing some unrealistic isolated-scale structures that are far away from the initial vortex and does not adversely affect the physically important object. The fine-scale structure is determined mainly by the fine-resolution area, although the systems may have larger differences before they move into the fine-resolution area. Altogether, this work demonstrates that the multiresolution configuration is a reliable and economic alternative to high-resolution global modeling. The adverse impact due to mesh transition and the coarse region can be controlled well.
Increasing resolution is generally regarded as an effective way to improve global weather and climate modeling (Jung et al., 2012; Wehner et al., 2014; Zhang et al., 2014; Yu et al., 2019). It is apparent that more computational and storage resources are required for higher-resolution models. This leads to a major challenge for efficient model development and application. The emergence of the locally refined, variable-resolution (VR) modeling approach offers a complementary route. The term VR is a broad concept. It may be realized with different styles, such as nested regional modeling with multiple grids, abrupt nonconforming mesh division, stretched grids, and the multiresolution approach. The stretched grid (e.g., Hourdin et al., 2006; Harris et al., 2016) and multiresolution approaches (e.g., Ringler et al., 2011; Guba et al., 2014) are close in terms of their conforming style. They maintain the global modeling configuration while permitting increased resolution for certain regions. The multiresolution approach is usually realized by an unstructured mesh model, such that a more flexible resolution choice can be achieved by considering multiple regions. Such a global-to-regional approach is the VR style that will be investigated in this study.
Numerical weather and climate modeling has shown that the VR approach can preserve the benefits of high-resolution applications for certain regions at a lower computational cost, as the total number of grid points can be greatly reduced (e.g., Sakaguchi et al., 2015; Skamarock et al., 2018; Gettelman et al., 2018). This advantage is especially valuable for high-resolution modeling that may reach the convection-permitting regime. While the global cloud- and storm-resolving (i.e., convection-permitting or allowing) modeling approach has been widely adopted (e.g., Stevens et al., 2019), it is still expensive and inefficient to frequently run such models for routine model development, research, and application. The VR model provides an efficient test bed for evaluating global model configurations and testing scale-aware physics. It offers flexible resolution configurations that may depend on physical interests. When properly formulated, it can be an intermediate and transitional step before establishing global convection-permitting modeling.
While the VR approach has shown some benefits, it may potentially suffer from some problems1. The nonuniform mesh, though it can be gradually refined, hardly decreases the global errors (as compared to its uniform counterpart that has the same degrees of freedom) because the truncation error is controlled mainly by the coarse-resolution region (Weller et al., 2009; Ringler et al., 2011); the numerical convergence rate may also be affected (e.g., Düben and Korn, 2014). Mesh refinement also tends to create artificial wave distortion and reflection. This issue is more challenging to the staggered finite-volume methods (Ullrich and Jablonowski, 2011), which are widely employed in today's weather and climate models due to their cost-effectiveness. At first glance, this seems to pose some disadvantages. Fortunately, the primary motivation of increasing resolution is to accurately resolve meteorologically important fine-scale structures. This implies that the solutions, in particular at high wavenumbers, change as the resolution increases. As long as one can maximize the model performance over the refined region, and have good control over the adverse impact due to the non-refined regions, the VR approach based on the staggered finite-volume method is extremely promising. This statement is not intended to diminish the importance of pursuing numerical precision, which is one of several important properties when developing model dynamics. It is hoped that a balanced compromise can be achieved, and we can thus take advantage of this promising approach.
Previous numerical studies have investigated the impact of grid refinement on the solution error in shallow-water models (Ringler et al., 2011; Guba et al., 2014), mainly based on single-region refinement. The impact of the width of the grid transition zone and the densification ratio has been emphasized. On the basis of spherical centroid Voronoi tessellation (SCVT), Ringler et al. (2011) demonstrated that the solution error is controlled primarily by the coarse-resolution region, and suggested that this can help to specify the coarse-mesh resolutions by determining what is an acceptable level of accuracy. They also suspected that the width of the transition zone may lead to increased errors. Liu and Yang (2017) suggested that the width of the transition zone may cause smaller additional errors compared to the increase in the densification ratio (for the tests they examined). Within a unified global model, these problems may potentially reduce the performance in the refined region.
In terms of resolving fine-scale fluid structures, the tropical cyclone is a useful test bed, and has been frequently used to examine resolution sensitivity. The global VR approach has been employed to simulate tropical cyclones for cost-effective climate simulations (Zarzycki and Jablonowski 2014). The VR model can capture smoother cloud patterns and smoother mid-level jet structures across the grid-refined region, leading to enhanced tropical cyclone activities compared to a nested model in which boundary forcing may have an adverse influence (Hashimoto et al., 2015). Based on a VR configuration, the Community Atmosphere Model (CAM) with a spectral element core (Taylor, 2011) maintains tropical cyclones crossing the transition zone well without discernable wave reflection (Zarzycki et al., 2014).
While these earlier studies have reported the benefits of VR modeling, a proper utilization of this technique is still challenging and deserves ongoing exploration. In this study, we systematically configure and evaluate the Global-to-Regional Integrated forecast SysTem (GRIST) atmosphere model based on the VR approach. GRIST is a new modeling system developed on an unstructured mesh, independent of existing atmospheric models available in the community. Previous studies have described the model formulation and evaluated its performance in shallow-water model tests (Zhang 2018; Wang et al., 2019), 3-D dry dynamical core (dycore2 hereafter; Zhang et al., 2019; Z19 hereafter) tests, and multiscale moist-atmosphere tests forced by simple physics (Zhang et al., 2020; Z20 hereafter). These studies mainly considered the quasi-uniform (QU) mesh. The model configuration and performance remain underexplored when local mesh refinement is considered.
In this study, we describe the model configuration for VR modeling. In particular, we will detail the explicit diffusion option and demonstrate its impact. We then examine the model behavior to understand its strengths and weaknesses under various mesh-refinement styles. This work is intended to provide a basis for utilizing GRIST-VR for more realistic modeling in future. To achieve these goals, we adopt two idealized initial atmospheric conditions, endorsed by the Dynamical Core Model Intercomparison Project (DCMIP; Ullrich et al., 2017). They drive the model towards some well-expected behaviors, facilitating a basic understanding of VR modeling. The model is forced from zero physics (i.e., pure dynamics) to simple and full physics, so as to represent an increasing degree of complexity.
The remainder of this paper is organized as follows. Section 2 presents the model description and its configuration for VR modeling. Section 3 describes the mesh configuration. Section 4 examines the VR performance in the dry baroclinic wave test. Section 5 investigates the model sensitivity in the tropical cyclone test. Section 6 presents a summary.
2.1 Model framework and dynamics
The model evaluated here is a frozen version of the GRIST-Atmosphere model. A20.9 denotes the version frozen at September 2020. The major descriptions (dynamical framework and component coupling) are in Z20 and Z19. GRIST is formulated on an unstructured mesh, which permits the use of SCVT (Ringler et al., 2008; Jacobsen et al., 2013) that enables VR modeling. A dry-mass-based generalized vertical coordinate is used. It allows flexible switching between the hydrostatic (HDC) and nonhydrostatic (NDC) cores. The moist-atmospheric model exactly conserves the dry air mass to within machine roundoff. The sink of the moist total energy is limited to a quite small value. The flux-form scalar variables are formulated in a layer-averaged manner. The momentum variables are formulated in their primitive forms. A vertically semi-implicit approach is used for solving the acoustic equations in the NDC, with explicit Eulerian vertical advection. The HDC is fully explicit.
The horizontal discretization is formulated on a hexagonal-C grid, i.e., using a staggered finite-volume method. Thuburn et al. (2009) proposed the key construction of the Coriolis term to achieve desirable mimetic properties. Ringler et al. (2010) formulated a set of spatial operators for the nonlinear shallow-water equations, under the constraints of integral invariant conservation and compatible vorticity dynamics. This approach has been used and examined by the Model for Prediction Across Scales (MPAS; Skamarock et al., 2012; Ringler et al., 2013), ICON-IAP (Icosahedral Nonhydrostatic model at the Institute of Atmospheric Physics, Gassmann, 2013), and DYNAMICO (Dubos et al., 2015). GRIST adopts two variations that differ from that in Ringler et al. (2010). Zhang (2018) extended a set of high-order upwind and center flux operators (Skamarock and Gassmann, 2011) for approximating the edge-based potential vorticity flux and demonstrated that both the higher nominal order and implicit upwind damping improve the simulated vorticity field. A pure third-order upwind formulation is used in GRIST. The other variation is a redefinition of the kinetic energy term by blending the original primal-cell value with a reconstructed value from the dual cell (Gassmann, 2013). This helps to alleviate the noise associated with the Hollingsworth instability (Hollingsworth et al., 1983) according to earlier QU model tests. A default coefficient of 0.9 is used following Eq. (20) in Z19.
The C-grid is cost-effective in dealing with the flux-divergence and gradient operators, which constitute the major horizontal computation involved in a fully fledged atmospheric dynamical core. The potentially adverse dispersion issue due to increasing mesh discontinuity in the VR mode (Ullrich and Jablonowski, 2011) can be well controlled by (i) using a smooth and gradual mesh transition (e.g., SCVT) and (ii) using a slight amount of explicit diffusion (as will be discussed). The basic horizontal operators are nominally accurate to the second-order, while the flux operator can be approximated using higher-order extensions. GRIST has several options for the flux operator. Among these, a nominal fifth-order upwind formulation can generate the smallest numerical error (see results in Zhang, 2018; Wang et al., 2019) but is not used for our default configuration as it requires three halo layers (the default number is two)3. This formulation is still instrumental in model development, helping to validate that the parallel computing infrastructure (Liu et al., 2020) is working correctly when using different minimum halo layers. The pure upwind formulation of Skamarock and Gassmann (2011) is used for the dycore. When combined with a flux limiter, this scheme can be used for tracer transport (as in Sect. 5.2). A Two-step Shape-Preserving Advection Scheme (TSPAS; Zhang et al., 2017; Yu, 1994) is also a major option for tracer transport (as in Sect. 5.1). Two-time-level single or multistage forward-in-time integration is used for dycore and tracer transport such that dry air mass and tracer mass are coupled in a consistent manner.
Two initial conditions are used in this study. The baroclinic wave (Jablonowski and Williamson, 2006; JW06) examines the adiabatic behaviors in a dry environment. The solution from a high-resolution run can be used as a reference solution. The idealized tropical cyclone is initialized following Reed and Jablonowski (2011). It is available from the DCMIP testing scripts. This test does not support a reference solution. As GRIST uses a dry-mass vertical coordinate, obtaining the moist–atmosphere state requires some special treatment (see Z20 for details). All the simulations use 30 full vertical levels with a top at ∼ 2.25 hPa, basically identical to the default CAM5 setup (e.g., Reed and Jablonowski, 2012). Both two tests are short-term deterministic tests (i.e., weather forecasting style), helping to validate the model configuration. Long-term climate modeling based on the VR configuration will be reported elsewhere. Also note that while some studies have pointed out that the vertical resolution should increase with horizontal resolution, we keep it unchanged in this study.
2.2 Model physics
GRIST provides a general physics–dynamics coupling interface to incorporate various physics packages. A tailored package can be used as a plugin, and its development can benefit from the broad community resources. One may add a specific physics scheme to an existing physics package or create an entirely new physics package, as long as it is compatible with the current workflow and is scientifically reliable. The surface model (e.g., land or a mixed-layer ocean model), though not used in this study, is coupled in a point-to-point style and can be shared by different physics configurations. Three physics packages are currently available as a basis for continuous research and development. These packages are separate in the sense that they have different physics drivers and data structures. For completeness, we describe them in this section.
The DCMIP simple physics package. In this study, the suite of Reed and Jablonowski (2012) is used. It contains a large-scale condensation process, a surface flux scheme, and a boundary layer process. The sea surface temperature is 29 ∘C globally. These processes are coupled in a time-splitting manner within the package, and the package is coupled to GRIST in a pure operator-splitting approach (ptend_f2_sudden; see Z20).
A climate physics package adopted from CAM5. This package is not used in this study, but the details can be found in Li et al. (2020). It is currently being tuned for the HDC that targets long-term climate modeling and can also be used for short-term integration in a weather forecast mode (e.g., Zhang et al., 2015).
A set of parameterization schemes from the Weather Research and Forecast (WRF) model (Powers et al., 2017). This package is being tuned for GRIST-NDC, which is targeted at simulating nonhydrostatic dynamics in a nonhydrostatic regime. The detailed schemes used in this paper include a six-species cloud microphysics scheme (Lin et al., 1983) from WRF version 2.0; the Tiedtke cumulus scheme (Tiedtke, 1989) from WRF version 3.7.1; the YSU (Yonsei University) planetary boundary layer scheme (Hong et al., 2006) and a surface scheme from WRF version 2.0; the RRTM (Rapid Radiative Transfer Model) longwave radiation scheme from WRF version 2.0 (Mlawer et al., 1997); and the CAM radiation module shortwave scheme from WRF version 3.4.1.
The internal coupling of this package is process splitting (see, e.g., a review in Gross et al., 2018, for details). All processes start from the dynamics-updated state and send back their respective tendencies to the dynamical model without modifications of the physics state variables. The exception is that microphysics will update the local physics state variables, so the calling sequence still matters. The physics–dynamics coupling uses a hybrid approach that combines the tendency method (ptend_rk) and the operator-splitting approach, as described in Z20. In particular, radiation heating is carried over the internal integration of the dycore as in a tendency method, and other tendencies (microphysics, boundary layer, cumulus) are updated as in a pure operator-splitting approach (ptend_f2_sudden in Z20). We emphasize that this package (including the choice of different schemes and its internal coupling) is still experimental and preliminary. It has not been comprehensively tuned. Ongoing tests, modifications, and additions are required to refine the performance. In this study, it is only intended to evaluate the particular performance of VR modeling and the role of explicit diffusion in a full physics scenario.
2.3 Time-step choices and explicit diffusion
For a VR model, the time step is theoretically restricted by its fine-resolution region. As emphasized by one reviewer and based on our own experience, some improper VR configurations may lead to higher stability restriction than their equivalent QU counterparts. For example, with regard to the acoustic-mode filter, Klemp et al. (2018) suggested that their original formulation is more problematic in VR applications, and cannot effectively remove the acoustic noise. Uncontrolled acoustic modes may artificially accumulate energy and thus impose a higher stability restriction. For the configurations and tests that we have examined, a suitable configuration of explicit diffusion is mostly relevant to this issue. The details will be given in the following section.
The explicit diffusion tendencies are generated at the largest time interval of each model step (i.e., physics step). They are coupled to model dynamics as in a tendency method (ptend_rk). Z20 showed that the tendency method has a slightly higher stability restriction than the pure operator-splitting approach (ptend_f2_sudden), but it can benefit from the higher accuracy of the time integrator. The tendency method does not require additional data communication. No explicit diffusion is activated for tracer transport in this paper.
2.3.1 Smagorinsky diffusion
Flow-dependent Smagorinsky diffusion (Smagorinsky, 1963) is used for velocities and potential temperature, following previous QU model tests. It is activated in all experiments, except the baroclinic wave test at the quasi-uniform G6 resolution. This scheme uses a second-order Laplacian operator multiplied by a flow-dependent eddy viscosity. The eddy viscosity is defined at the edge point. The Smagorinsky diffusion is not very scale selective, but its flow-dependent feature makes it selective in terms of where and how much to diffuse (see e.g., Fig. 9 in Gassmann, 2013). The diffusion strength is acceptable overall, as evidenced by the sharp gradient in the QU baroclinic wave test.
In the QU mode, a mean length scale is used for calculating the eddy viscosity. For the VR mode, doing so implies a stronger diffusion (than typically required) for the refinement region. This also leads to a higher stability restriction, especially when the refinement ratio becomes large. In this version, the square of this mean length scale is replaced by the local length product of a pair of crossing edges. In the tropical cyclone test (simple physics), this local approach increases the maximum wind magnitude in the eyewall (as compared to the unscaled version in the preprint), because a smaller amount of diffusion is imposed on the refinement area (see Sect. 5.1). It also reduces the parametric sensitivity to the Smagorinsky coefficient as found in the preprint. When varying the coefficient, the present version generates more consistent solutions, similar to a QU model. This local scaling approach was used for the VR mode but not for the QU mode.4
The hyperdiffusion option was not used in the initial preprint. Based on some exploration and experiments, we have found that activating scale-selective fourth-order hyperdiffusion for the horizontal velocity shows demonstrable added value for VR modeling. First, in the baroclinic wave test, in the absence of hyperdiffusion, a higher Smagorinsky coefficient (compared to the equivalent QU test) is required to suppress grid scale noise due to the mesh refinement, which in turn restricts the numerical stability. When the hyperdiffusion is activated, we can use a moderate Smagorinsky coefficient that does not challenge the numerical stability and maintains a quality solution. The DTP (dycore–tracer–physics) splitting mode benefits from this most because diffusion is called at the step of the model physics. Second, even with a higher coefficient, the Smagorinsky diffusion is inactive over certain regions where weak flow deformation dominates. In the baroclinic wave test, some grid scale oscillations are more conspicuous over these regions (see the preprint). Such noise is due to the mesh transition, and some of it is akin to the Hollingsworth instability. As will be shown, activating a background hyperdiffusion for the horizontal wind successfully removes such noise. The solutions are overall less oscillatory than those in the preprint. Third, in the tropical cyclone tests with full physics, which presents more nonlinear feedback, the hyperdiffusion option is also effective in suppressing some isolated-scale structures. Unlike the minor disturbances that may be generated near the major tropical cyclone in the simple physics test, these systems are far away from the initial vortex, and they are thus highly unrealistic. For these reasons, the hyperdiffusion option is used for the VR mode, but not for the QU test. Thus, QU and VR models are applied by different forms of explicit diffusion in this paper. Future work may also need to examine the possible impact of hyperdiffusion in a QU model to better isolate its effect.
The hyperdiffusion operator is formulated by recursively using the Laplacian operator (see Z19 for details). The diffusion coefficient can be determined in a relatively empirical way. For VR modeling, we adopt the approach documented in Zarzycki et al. (2014) for scaling the coefficient:
The reference length Δxref and reference viscosity coefficient K4(Δxref) are empirically determined. This formulation reduces K4(Δx) by a factor of 10 for every halving of resolution. A similar scaling approach is also used in MPAS (Skamarock, 2016). We typically use the configuration for a G-level resolution that is close to the finest resolution on the mesh. Some typical values for GRIST are documented in Table S1 in the Supplement. These values5 are smaller than those in Zarzycki et al. (2014) for the corresponding length scale by a factor of 5. The local grid distance (Δx) is an average distance between the grid point and all its nearest neighbors, i.e., a cell-based value. The edge-based value used for hyperdiffusion is an average of two neighboring cell values.
The properties and generation of the SCVT are detailed in Ringler et al. (2011) and Ju et al. (2011). We focus on two key elements of the SCVT: generators and density function. A spherical Voronoi tessellation is a spatial subdivision of a sphere Ω based on a set of distinct points on Ω. For each point xi, , the corresponding Voronoi region Vi, , is defined by
where ∥⋅∥ denotes the geodesic distance. Each point xi is called a generator and its corresponding Voronoi region Vi is called the Voronoi cell. A spherical Voronoi tessellation becomes an SCVT when the generators are also the centroids of the Voronoi cells. In this study, the SCVT mesh is constructed by an iterative process based on Lloyd's algorithm (Du et al., 1999). In particular, a parallel algorithm is used (Jacobsen et al., 2013) to avoid time-consuming serial construction. That said, generating a quality VR SCVT is still nontrivial (but only done once). In our implementation, the iteration stops when two criteria are satisfied: (i) it reaches an empirically determined minimum step and (ii) the circumcenter of each triangle falls within its shape.
Three original generators are used in this study.
Icosahedron bisection. This approach benefits from the excellent uniform properties due to bisections of a regular icosahedron. The mesh resolution is referred to as G-level/Gn, where n denotes the number of bisections. After each bisection, the total grid number is approximately four times greater than the previous one.
Icosahedron bisection with a final-step trisection. Instead of using recursive bisection, a trisection is used at the final step, to achieve an intermediate resolution between two neighboring G-level resolutions. This mesh is referred to as GnB3 (i.e., n bisections plus one trisection). For instance, the resolution of G5B3 (∼ 80 km) is between that of G6 (∼ 120 km) and that of G7 (∼ 60 km). The number of added primal cells (mainly hexagons) from G5 to G5B3 is equal to four times the number of dual cells (triangles) in G5.
Spherical uniform random (SUR) set of points. Initial points are created uniformly on the sphere by the Monte Carlo method. In this way, the original generators can be obtained by using an arbitrary number of points, not restricted to a sub-divided icosahedron.
3.2 Density function
By specifying the density function, the SCVT is able to precisely control the distribution of the local resolution. For any two Voronoi regions indexed by i and j, the conjecture is
where ρ(xi) is the density function evaluated at xi, and dxi measures the local mesh resolution. This relation is valid in a theoretical sense.
A QU mesh can be constructed when the density is one on the sphere, and the Voronoi regions are approximately equivalent to each other. For the VR mesh, the basic density function used in this study is
where ∥xrc−xi∥ denotes the geodesic distance between the location of the refinement center and each generator. α indicates the width of the transition zone between the fine-resolution and coarse-resolution regions, β defines the coverage radius of the fine-resolution region, and γ measures the densification ratio between the finest and coarsest resolutions. A sample X4 mesh based on this density function with is shown in Fig. 1a.
Because the basic density function is fixed to a single-region refinement, we adjust it for multi-region refinement. The multi-region refinement is divided into two styles based on the refinement centers. In the hierarchical refinement way (Fig. 1b), we add a uniform intermediate-resolution region between the inner fine-resolution and the outer coarse-resolution regions (see Fig. 2). This helps to avoid an overly high densification ratio between two neighboring regions. Equation (4) can be generalized to a form that allows us to control the resolution of the intermediate region:
λ is designed to control the resolution of the intermediate-resolution region, also referred to as the second-refinement region () that is located between the first-refinement () and the coarse-resolution regions (dxc):
Generally, . The function of γ is similar to that in the previous single-region refinement, except that the fine-resolution region is referred to as the first-refinement region here:
Corresponding to γ, we refer to the meshes that are generated based on λ values of (1)4, (1∕2)4, and (1∕3)4 as XL1, XL2, and XL3 meshes, since the resolutions of the first-refinement and the second-refinement regions vary according to the inner densification ratios of 1, 2, and 3, respectively. For example, when γ is fixed at X4, the hierarchical meshes based on G6 are called G6X4L1, G6X4L2, and G6X4L3 meshes.
In a polycentric refinement mode (Fig. 1c), by adding a different refinement center xrc2, the density function of the polycentric refinement mode is defined as follows:
The geodesic distance between the two refinement centers should satisfy .
4.1 Single-region refinement
The dry-atmosphere test examines the pure numerical solution of the model. It does not include the nonlinear interaction between dynamics, moisture transport, and parameterization. Based on the JW06 baroclinic wave test, we first compare the VR and QU simulations. Previous studies that employed this test for a VR model include Gettelman et al. (2018; using the multiresolution approach) and Harris and Lin (2012; using the multi-grid approach), based on different evaluation metrics. The QU grids include G6 (∼ 120 km; 40 962 cells), G7 (∼ 60 km; 163 842 cells), and G8 (∼ 30 km; 655 362 cells). The VR grids examine all three generators: (i) icosahedral bisection (e.g., G6X4), (ii) final-step trisection (G5B3X4), and (iii) spherical uniform random points (SURX4). The minimum iteration number is 300 000 for G6X4 and 1 000 000 for G5B3X4 and SURX4. The refinement center is placed at 35∘ N, 180∘ E, with and . The detailed model configuration is given in Table 1. The time step of the VR model is limited by its fine-resolution region and is set accordingly, although the currently used time steps do not represent the maximum allowable step.
We first examine the ability of resolving the fine-scale structures. Figure 3 shows the relative vorticity field at the model level nearest to 850 hPa (level 24) after 10 d. These values have been remapped to the Voronoi cell from the raw triangular grid, so as to avoid the aliasing of certain oscillatory patterns. Such patterns (see preprint) actually reflect the mesh shape and are more conspicuous for coarse resolution. Figure S1 in the Supplement further shows the QU model results interpolated to the regular longitude–latitude grid as a reference. As the resolution increases, the QU models simulate stronger vortices with a clear filament structure (Fig. 3a–c). The VR model can simulate the smooth structure of the waves in the refined region, as in the high-resolution QU model. In the VR mode, two vortices in the west fall within the fine-resolution region. The structure of the westernmost vortex in G6X4 (Fig. 3d) is close to G7. G5B3X4 and SURX4 (Fig. 3e and f) further produce finer-scale structures than G6X4, closer to G8. The easternmost wave falls within the transition zone in three VR runs. The variation in the mesh sizes there does not distort the wave pattern, and the fine-scale structure can be improved as the resolution of the transition zone increases (e.g., from G6X4 to G5B3X4). In the SURX4 test, a minor roughness is found on the tails of the vortices. This deficiency is largely due to local mesh irregularities. Compared to the icosahedron-based SCVT generators, the random generators tend to slightly degrade the mesh quality and the simulation performance.
To present a more quantitative evaluation of the momentum field, we examine regional kinetic energy spectra6 over the refinement area (red box in Fig. 3). We use the discrete cosine transform to perform this analysis (Denis et al., 2002). The computational procedure is briefly documented in the Appendix. The decomposition is made for relative vorticity and divergence. Figure 4a–c show the results from different tests. It is clear that the rotational mode dominates the kinetic energy. The hydrostatic model (Fig. S2 in the Supplement) produces similar results to the nonhydrostatic model; thus only the nonhydrostatic core is discussed in the main text.
At the first wavenumber, runs with different resolutions show larger discrepancies. This is because this lowest mode absorbs most of the large-scale trend and corresponds to only a half-cosine wave. Denis et al. (2002) suggested that one may remove this mode if desired, but we keep it here. From the 2nd to the 10th wavenumbers, all the curves are consistent overall, suggesting that the well-resolved structures are robust to various tests. The major difference lies in the high wavenumbers. For the VR results, G6X4 is close to G7. This is expected because they have similar resolution over the selected domain (∼ 60 km). For the same reason, G5B3X4 and SURX4 produce better spectral tails than G7 and G6X4, closer to G8. This confirms that increasing the fine resolution of the VR model is able to improve its fine-scale resolving ability. At such high wavenumbers, SURX4 is even closer to G8, as it has slightly higher energy in the tail. However, this actually reflects that SURX4 has slightly more grid scale oscillations than G5B3X4. Therefore, an examination of kinetic energy spectra should be accompanied by a close look at the real field.
In the context of kinetic energy spectra, it is useful to demonstrate the impact of the hyperviscosity coefficient. For G5B3X4, we vary the reference hyperviscosity coefficient under a fixed reference length (30 km). Note that 2 × 1012 is the default configuration for this test. As shown in Fig. 4d–f, when using 2 × 1013, spectra are seriously damped at the high wavenumbers in both the rotational and divergent components. A flat tail is generated. This indicates that the diffusion strength is overly strong for the fine-resolution area. The other four tests generally produce consistent results: reducing the coefficient tends to slightly uplift the tail, i.e., increase the kinetic energy. While the lowest coefficient (2 × 1011) seems to produce a nicer tail than the default run, this again reflects the fact that slightly more small-scale oscillations are generated and that the solutions are less clean (but still acceptable). Therefore, when tuning a VR model, one should achieve a minimally required hyperviscosity that neither significantly damps the field nor becomes unable to suppress certain grid scale oscillations. For this test, 2 × 1012 is fairly close to the optimal choice.
Next, we assess the solution errors by first examining the surface pressure field. Its global distribution on day 9 simulated by GRIST-NDC is shown in Fig. S3 in the Supplement. The locations and magnitudes of the high- and low-pressure centers are consistent overall in the QU (Fig. S3a–c) and VR (Fig. S3d–f) runs. There are some nonzero wave patterns in all VR simulations over the Southern Hemisphere, reflecting that the nonuniform grid structure leads to higher truncation errors. G5B3X4 has smaller wave patterns than G6X4. This is expected as G5B3X4 has a higher resolution than G6X4 there. Moreover, G5B3X4 is also better than SURX4 with regard to these wave patterns. Considering that SURX4 has more mesh cells, this suggests that generators based on a regular icosahedron produce higher mesh quality than random generators, given roughly the same number of iterations.
Figure 5 shows a quantitative comparison. The global l2 error norm of each test is computed against the highest resolution (G8). The errors of three VR meshes (G6X4, G5B3X4, and SURX4) are higher overall than those of the QU meshes. Note that these errors are also slightly higher than those in the preprint because a larger time step is used in this version (we performed additional tests to check this sensitivity; figure not shown). The nonhydrostatic solver produces slightly smaller errors using three VR meshes, but the overall accuracy of the two solvers is comparable. For three VR meshes, G5B3X4 shows the smallest error and is generally close to the results of G6 during the first 10 d. G6X4 shows the largest error and SURX4 lies in between them. Again, the fact that SURX4 is less accurate than G5B3X4 implies that random generators are more likely to be trapped into the local area during the iterative procedure of mesh generation, leading to a higher degree of local irregularity.
For GRIST-NDC, we further tested G7X4 and G8X4. As shown (Fig. 5b), G7X4 produces smaller errors overall than G5B3X4, although the difference is small. G8X4 produces higher errors than G7X4 during the first 4 d. We suspect this is because certain initial imbalance becomes more active in a high-resolution VR configuration. Such imbalance can be caused by the discrete initialization, for example, the continuous properties imposed by the analytic wind field will not be exactly satisfied by the discrete normal velocity, unless the velocity is obtained based on some constraints. After day 5, G8X4 produces smaller errors than G7X4, and the increment in error reduction is close to the difference between G6X4 and G7X4. G8X4 also produces clearly smaller errors than G6 from day 5 to day 12. After day 13, G8X4 has higher errors than G6, although its coarsest resolution is finer than that of G6. These results suggest that the VR models indeed increase the global errors, but the errors can be reduced by increasing the coarse resolution of the mesh. This is not surprising but a reconfirmation of the conclusion in Ringler et al. (2011), using a 3-D atmospheric dynamical core.
The JW06 test was originally proposed in the era that models based on a regular latitude–longitude or Gaussian grid were predominant. For a quasi-uniform grid, it is well known that the steady state in the Southern Hemisphere (or in an unperturbed condition) cannot be perfectly maintained due to mesh irregularity (Lauritzen et al., 2010). In a baroclinic environment, such errors will grow exponentially to break the steady state. Because mesh irregularity increases in a VR mode, the inability to maintain the steady state can be further amplified, ultimately contributing to the increased errors.
To reduce the impact of this issue, we add another perturbation over the Southern Hemisphere to excite two wave trains, as was done in some earlier studies (e.g., Gassmann, 2013). Figure 6 shows the relative vorticity field on day 10 from three selected tests: G6, G6X4, and G5B3X4. The sign of the relative vorticity in the Southern Hemisphere is flipped to facilitate a visual comparison. The QU model (Fig. 6a) produces a comparable wave train in each hemisphere. The wave trains are not exactly symmetrical because the mesh cells are not exactly symmetrical across the Equator. In the Northern Hemisphere, G6X4 and G5B3X4 (Fig. 6b and c) produce similar solutions to Fig. 3d and e. In the Southern Hemisphere, the solution of G5B3X4 is closer to that of G6 than G6X4 in terms of wave pattern and magnitude. This is to be expected as G5B3X4 has a higher resolution in the Southern Hemisphere. Although G6X4 cannot simulate the structure in the Southern Hemisphere as G5B3X4, no serious problem is found for that region.
Figure 7 presents a quantitative estimation of the solution error. Generally, the magnitude and growth of the errors are close to those in Fig. 5b. A notable difference is that G8X4 produces smaller relative errors, closer to those of G7 from day 8 to day 12. On day 15, G8X4 still shows comparable errors to G6. The error reduction from G6X4 to G7X4 also becomes larger than that in Fig. 5b, when the waves become more developed (e.g., day 8 to 11). This implies that the inability to maintain a steady state in the Southern Hemisphere indeed worsens the error estimation for a VR model.
4.2 Multi-region refinement
To increase the application range of the single-region refinement, two multi-region refinement modes are examined to obtain desired resolutions in multiple regions. This represents a unique aspect of the multiresolution approach. The first way is a hierarchical style with one refinement center. It contains three consecutive uniform sub-regions outside the refinement center: the first-refinement region, the second-refinement region, and the coarse-resolution region. The second-refinement region provides an intermediate resolution between the first-refinement and coarse-resolution regions. Compared to the single-region refinement, this second-refinement region provides more uniform resolution between the refinement center and the coarse-resolution area.
The symmetrical perturbation test was performed using a G8X4L2 mesh. On day 10, the first wave (the strongest) over the Northern Hemisphere is experiencing a higher gradient of resolution than that over the Southern Hemisphere. The fine-scale structures are well captured by the VR model (Fig. 8a). The second and third waves in the Northern Hemisphere have stronger magnitudes than their equivalents in the Southern Hemisphere. The difference in each wave train generally reflects the differences in the local resolution.
The second way uses a polycentric style with multiple refinement centers. The same double-baroclinic wave trains are generated using a G7X4 mesh, with two refinement centers at 35∘ N, 180∘ E and 35∘ S, 180∘ E, respectively. On day 10, the model well simulates the fine-scale structures over two selected regions (Fig. 8b). The transitions between the fine and coarse regions are smooth and stable. This refinement mode may provide an effective way to simultaneously improve the simulations over different regions. The wave train in the Southern Hemisphere coincides well with that in the Northern Hemisphere. Again, they are not exactly identical because the mesh cells are not exactly symmetrical across the Equator.
In this test case, it is useful to demonstrate the impact of activating the hyperdiffusion option. As shown in Fig. 8c and d, when hyperdiffusion is turned off, the simulation results become rather oscillatory. The noise pattern at the tail of the wave is akin to the Hollingsworth instability and is amplified by a VR model due to mesh transition. The noise is more conspicuous in the first-refinement region of G8X4L2, due to higher resolution and rapid mesh transition. Also note that these results are more oscillatory than those in the preprint because the Smagorinsky diffusion in that version is stronger. As mentioned in Sect. 2.3, the Smagorinsky option is still inactive over certain regions. Using a background hyperdiffusion is a good complement to this deficiency.
5.1 Simple physics
Moist-atmosphere modeling includes the nonlinear interaction of dynamics, moisture transport, and model physics. We use the idealized tropical cyclone test (Reed and Jablonowski, 2011) because of its clear resolution sensitivity. This test is useful to examine the VR performance because the solution does not fully converge even at 10 km (Z20). The simulated tropical cyclone is very sensitive to the mesh size. A higher-resolution model will produce more intense storms. An initial vortex is placed 10∘ N, 180∘ E with no background flow. The tropical cyclone moves northwestward due to beta drift. The initial virtual temperature profiles are designed to be conditionally unstable in the troposphere. A small perturbation (either physical or computational) is more likely to excite new storms. Whether these additional signals are realistic depend on the situations.
Previous studies have shown that model dynamics have a clear impact on the simulations of tropical cyclones. For example, Zhao et al. (2012) showed that increasing the strength of 2-D divergence damping in the finite-volume dynamical core (Lin, 2004) leads to more occurrences of tropical cyclones. Reed et al. (2015) showed that in CAM5 the spectral element core produces stronger tropical cyclones than the finite-volume core, when the parameterization suite remains almost unchanged. For GRIST, a known sensitivity is that different tracer transport options may lead to different wind magnitudes in the eyewall, as the full pressure gradient term is, by design, related to tracer mixing ratios (Z20).
We first examine the model behaviors under a simple-physics environment (Reed and Jablonowski, 2012). In the preprint, we compared two representative groups of numerical tests: one group based on the NDC with DTP splitting enabled and one based on the HDC with no DTP splitting (i.e., dycore, tracer transport, and physics use the same time step). Results from these two groups are consistent overall. The impacts of DTP splitting and hydrostatic and nonhydrostatic options do not generate discernable differences across various VR meshes. This is consistent with the previous QU model tests (Z19; Z20) in that (i) the nonhydrostatic solver behaves similarly to its hydrostatic counterpart under the hydrostatic regime and (ii) the DTP splitting does not cause a degeneration in performance when it is properly configured. In this version, which uses an updated configuration, only the NDC is tested. The configuration for the simple physics test is given in Table 2.
We first examine the evolution of the tropical cyclone when it moves across the transition zone in the hierarchical refinement mode. The mesh is fixed at X4, with and . In the control run (G6X4L2; Fig. 9a), , , and . The first-refinement region of this VR mesh has ∼ 40 km resolution. Two cases are used to examine the impact of the mesh distribution. In the first case, we choose to use more rapid resolution changes in the inner transition zone based on the two parameters α1 and λ. α1 controls the width of the inner transition zone, and λ represents the densification ratio between the first- and second-refinement regions. Either narrowing the width of the transition zone (Fig. 9b) or enlarging the inner densification ratio (Fig. 9c) generates a more abrupt transition zone. The other way is to have the transition zone affect the tropical cyclone at an earlier stage. The initial cyclone is placed closer to or even partly within the transition zone by increasing β1 (Fig. 9d), a parameter that controls the radius of the first-refinement region. The tropical cyclone is initialized at 10∘ N, 180∘ E over the second-refinement region, near the transition zone between the first-refinement and the second-refinement regions.
On day 10, all four tests simulate the shape of the tropical cyclone well. During its movement from the second-refinement into the first-refinement region, the change in the grid size leads to little distortion on the tropical cyclone. Compared to the preprint, two major differences can be found. First, the minor disturbances near the major tropical cyclones almost disappear in this version because activating hyperdiffusion damps the wind field more effectively. Though it can be damped, this minor disturbance is not that unrealistic because it is near the major cyclone. During the movement, it is possible that the nonlinear feedback can generate new small-scale systems. Similar minor disturbances have also been observed in some models participating in DCMIP2016 (see e.g., https://www.earthsystemcog.org/projects/dcmip-2016/, last access: 6 October 2020). The other difference is that the maximum wind magnitude in the eyewall is stronger overall in this version. This is due to the locally scaled Smagorinsky formulation, as mentioned in Sect. 2.3.1.
Figure 10 presents the evolution of the tropical cyclone in each test on days 2, 4, 6, and 8. The tropical cyclones in all tests are consistent. No discernable difference is found when they move across the transition. A slightly larger difference is more evident in the early stage (day 2), but such differences diminish as the tropical cyclones move into the refinement center. The relative vorticity field also looks smooth and does not show any artifacts (Fig. S4 in the Supplement). This indicates that the mesh transition does not create notable problems.
To further examine possible sensitivity, we performed a group of experiments by altering one of the three parameters: α1, β1, and λ. Figure 11 shows the minimum surface pressure (Fig. 11a–c) and the maximum wind speed at 850 hPa (Fig. 11e and f). All these tests use a DTP splitting number , with a dycore step of 60 s (Table 2). Only one test with the smallest α1 () is unstable at this DTP splitting number, so we adjusted it to 1 : 2 : 4 with the same dycore step. This suggests that an overly narrow inner transition zone can impose a higher stability restriction. The tropical cyclone rapidly strengthens during the first 2 d. After day 2, it enters into a moderately developing stage in each experiment. The evolution of intensity is diverse. All the VR runs tend to produce stronger tropical cyclones than G6 or G7 because the fine-resolution area determines the ultimate strength. The tests have run-to-run differences, but they are small overall, showing consistent results and robustness.
Figure 12 shows a further comparison between two VR runs (G6X4L2 and G7X4L2) and two QU runs (G7 and G8). The highest resolution of G6X4L2 and G7X4L2 is slightly higher than that of G7 and G8, respectively (see Table 2). As shown, in each comparison (G6X4L2 vs. G7, G7X4L2 vs. G8) the VR model produces stronger wind magnitudes in the eyewall and a more compact size featuring a smaller area coverage. The eyewall of the cyclone converges towards its center, almost within 1∘ from the center. The maximum winds develop to higher vertical levels as the local resolution increases. Overall, the difference between the VR and QU tests is attributed to different local resolutions.
In the tropical cyclone test, if an unscaled Smagorinsky eddy viscosity is used (Fig. S5 in the Supplement; preprint), the VR model will show a higher parametric sensitivity for the Smagorinsky coefficient than the QU model. The current version has largely reduced such sensitivity, closer to the behavior of the QU model (Fig. S5). This confirms the reasonable behavior of the local length formulation.
5.2 Full physics
The simple physics test, while insightful, is limited in the sense that the physics processes are simplified and do not support enough of the nonlinear feedback typical of a real-world model. A VR model may introduce a higher degree of nonlinear feedback due to mesh refinement. Thus, it is useful to check its behavior given a full parameterization suite. We use a more highly variable mesh: G6B3X16L4. The finest part on this mesh reaches ∼ 7–10 km. The refinement center is placed at 35∘ N, 165∘ E. The sea surface temperature is identical to that in the simple physics test (29 ∘C uniformly). The starting date is the first day of June. The solar constant is 1370 W m−2. The DTP splitting number is 1 : 5 : 10 with a dycore step of 10 s. If activated, the square of the Smagorinsky coefficient is 0.005, and the reference hyperviscosity is 2 × 1010 m4 s−1 with a reference length scale of 7000 m.
We first focus on the impact of the explicit diffusion process. Figure 13 shows the tropical cyclones on day 10 from three tests: no explicit diffusion (noDiff), using hyperdiffusion only for the horizontal velocity (hyper), and hyperdiffusion plus Smagorinsky (hyper + smg). The wind speed is shown at the model level nearest to 850 hPa (level 24). On day 10, the tropical cyclone center just moves into the refinement center (35∘ N, 165∘ E). In general, the results in the three tests are consistent. The major difference lies in the eyewall. The noDiff and hyper tests produce slightly higher wind maxima than hyper+smg. This suggests that Smagorinsky diffusion becomes active in the eyewall, diffusing the solutions a little bit more strongly. Except for this minor difference (see the zoomed results), hyper and hyper+smg produce very close solutions.
Explicit diffusion, while it seems to be irrelevant to the performance over the fine-resolution region, plays a more important role in the coarse-resolution region. Due to the globally uniform sea surface temperature, near the South Pole, some isolated systems can be generated in the noDiff test (Fig. 14a). Unlike the possible minor disturbances near the major tropical cyclone in the simple physics test, these systems are not expected there because they are far away from the initial vortex. They are likely to be caused by the nonlinear feedback between dynamics and physics. Only using the Smagorinsky diffusion can suppress them to a certain extent, but it is not enough due to its flow-dependent nature (Fig. 14b). By further activating the hyperdiffusion option, these systems can be effectively removed (Fig. 14c). This suggests that the explicit diffusion configuration is able to achieve a balance between underdiffusion and overdiffusion.
With this X16 mesh, it is also useful to examine whether a single refinement with a high densification ratio (G6B3X16) tends to create serious problems during the simulation, especially when the tropical cyclone moves across the transition. The refinement center is rotated to 60∘ N, 165∘ E. The tropical cyclone will encounter more mesh transitions in G6B3X16. The results are shown in Fig. 15. On day 2, the tropical cyclone in G6B3X16L4 is experiencing mesh cells between ∼ 30 and ∼ 50 km. In G6B3X16, the encountered mesh sizes are above 70 km. Hence, G6B3X16L4 has a stronger and more developed eyewall at this stage. On day 6, the eyewall enters into a region with cell sizes ∼ 25–30 km in G6B3X16L4 and ∼ 50 km in G6B3X16. G6B3X16L4 still has a more compact eyewall, but the difference between the two runs becomes smaller. On day 10, the eyewalls in both tests enter into an area with cell sizes ∼ 20–25 km. At this stage, both cyclone systems show similar maximum wind magnitudes in the eyewall, and exhibit a similar distribution. Clearly, the fine-scale structure is determined mainly by the fine resolution, although two systems can have larger differences before they enter into the fine-resolution area. These results also demonstrate that a single refinement with a high densification ratio, though more challenging, can perform competitively with a more gradually refined mesh when properly configured (e.g., an overly narrow transition is undesirable). The choice of different mesh styles thus demands more numerical modeling experience and depends on the target issue. This needs to be further examined with more realistic weather and climate system.
In this study, an atmospheric model formulated on an unstructured mesh has been systematically configured and evaluated in its VR mode. Different mesh-refinement styles have been utilized to evaluate the model performance under increasing degrees of complexity, from dry dynamics to simple physics and finally full physics. Based on some additional NDC tests with full physics, it is found that the QU and VR models (with the same degrees of freedom) have comparable parallel efficiency; i.e., the VR configuration does not cause a degeneration in strong scaling performance. This is to be expected as the domain decomposition does not distinguish between QU and VR grids. Other scaling performances (e.g., whether the speedup ratio can achieve the ideal value by simply reducing the total grid number, with all other things being equal) can be found in Liu et al. (2020). Overall, these results demonstrate that the VR configuration of GRIST is a reliable and economical alternative to high-resolution quasi-uniform modeling. The adverse impact due to the mesh transition and the coarse-resolution area can be controlled well. The major conclusions may be summarized into three aspects.
6.1 On the overall performance
Based on the dry baroclinic wave test, all VR styles can capture the fine-scale wave structures well. Such fine-scale resolving capability is supported by analysis of regional kinetic energy spectra and is further verified in the multi-region refinement mode. In the transition zone, the waves are not adversely affected by the mesh refinement. In the coarse-resolution region, the VR model can also simulate an equivalent distribution of waves to its low-resolution counterpart. A VR model indeed produces greater solution errors compared to its QU counterpart that has the same degrees of freedom. However, the solution error can be reduced when the coarse region increases its resolution. Thus, its impact can be controlled. In the tropical cyclone test, the VR model in simple physics and full physics can simulate the gradual evolution of the tropical cyclone, showing reasonable resolution sensitivity. Importantly, the simulation of the fine-scale structure is mainly controlled by the fine resolution, although larger differences may exist before the systems move into the refinement area.
6.2 On the impact of explicit diffusion
Comparing this version with its earlier version (preprint), the impact of explicit diffusion can be clearly demonstrated. It has been shown that activating scaled hyperdiffusion for the horizontal velocity as background diffusion is helpful in alleviating grid scale oscillations due to mesh transitions. It also suppresses some highly unrealistic disturbances found in the full physics scenario. When activating hyperdiffusion, it is suggested that a minimally required hyperviscosity should be used, such that it neither significantly damps the field nor becomes unable to suppress undesirable grid scale disturbances. Using only Smagorinsky diffusion would require higher coefficients to remove oscillations, which, in turn, would restrict the numerical stability. Even so, the Smagorinsky-only configuration is not enough for the VR mode because it becomes inactive over certain regions that are dominated by weak flow deformation. It is also shown that, in the tropical cyclone test, the scaled Smagorinsky diffusion has much lower parametric sensitivity than its unscaled counterpart. In general, the use of (properly scaled) explicit diffusion plays a positive role in the VR configuration, although there were concerns that their own discretization may introduce additional problems. The increased demands of explicit diffusion from the QU to VR configuration is in accordance with increasing mesh discontinuity due to transitions, and thus this should not be viewed as a disadvantage.
6.3 On the impact of the mesh styles
Different mesh styles are all able to produce fine-scale structures. They maintain solution quality over the mesh transition and the coarse-resolution area. In the dry test, it has been shown that G5B3X4, while it has fewer cells than SURX4, leads to smaller solution errors and less oscillatory solutions. As these two meshes are generated given a similar number of iterations, this might suggest that the random generators need more computation to achieve a quality mesh. In general, when generating a VR mesh, we empirically suggest that one would be better to start from a subdivided icosahedron. In the tropical cyclone case, a series of sensitivity tests examined the model performance in a hierarchical refinement style. The solutions exhibit consistency even when the VR mesh is slightly perturbed by one of the three parameters that control the density function. It is suggested that an overly narrow transition zone should be avoided. In the full physics test, G6B3X16L4 and G6B3X16 produce consistent results on day 10, although the simulations at the initial stage (e.g., day 2) exhibit larger differences.
The model data are first interpolated to the Gaussian grid. For the QU model, we use the Gaussian grid that is close to the nominal resolution of the icosahedral mesh (T106 for G6, T213 for G7, and T426 for G8). For the VR model, the data are interpolated based on the fine resolution of the mesh (T221 for G6X4, T328 for G5B3X4, and T328 for SURX4). The selected regional domain (red box in Fig. 3) ranges from 30∘ N to 65∘ N and from 150∘ E to 150∘ W with two vortices.
We follow the approach documented by Denis et al. (2002) for computing kinetic energy spectra on a selected regional domain. Let Fζ(m, n) be the discrete cosine transform (DCT) of a 2-D relative vorticity field ζ(i, j) of Ni-by-Nj grid points. The variance array can be computed from the DCT field as follows:
where , , and . Each 2-D wavenumber pair (m, n) is associated with a wavelength:
where Δ is the grid spacing. μ is a normalized 2-D wavenumber defined as follows:
To construct a spectrum, the variance contributions of Fζ(m, n) needs to be binned by bands of μ. The ranges of μ for a given wavelength band are determined as follows:
where denotes the wavenumber.
The rotational part of a spectrum as a function of wavenumber k can be obtained by binning each element in the variance array. For a given k, Eqs. (A4) and (A5) are used to calculate the lower and upper bounds. For each variance element, if the corresponding μ satisfies , the variance will be summed to give the spectrum for this k, i.e.
where is the circular wavenumber:
Similarly, the divergent part can be evaluated based on the divergence field:
where is the variance array of divergence. The total kinetic energy as a function of wavenumber k is given as follows:
In this study, kinetic energy spectra are displayed as a function of wavelength A and wavenumber k. It should be noted that the lowest mode (k=1) absorbs most of the large-scale trend and can be excluded if needed (Denis et al., 2002).
GRIST is available at https://github.com/grist-dev (last access: 7 December 2020). A version of the model code and running and postprocessing scripts for supporting this paper are available at https://doi.org/10.5281/zenodo.3930643 (GRIST-Dev., 2020a). Version 2 is for this paper, and version 1 is for the preprint. The preprint is available from the online link of this paper. The grid data used to enable the tests are located at https://doi.org/10.5281/zenodo.3817060 (GRIST-Dev., 2020b). Public code access is available after authorization, following https://github.com/GRIST-Dev/TermsAndConditions (last access: 7 December 2020).
The supplement contains Table S1 and Figures S1–S5. The supplement related to this article is available online at: https://doi.org/10.5194/gmd-13-6325-2020-supplement.
YihZ performed the mesh generation, most of the data analysis, and an initial exploration of VR modeling, with input and guidance from YiZ. RY supervised YihZ's work, provided impetus and resources, and guided the experimental design. JL supervised the model development project. ZL is responsible for parallel computing and customized the parallel mesh generation software. YiZ designed this study, developed and maintained the model, performed production runs, and analyzed the simulations. YiZ and YihZ wrote the paper, with contributions from all authors. All the authors continuously discussed the model development.
The authors declare that they have no conflict of interest.
The constructive comments from two anonymous reviewers are highly appreciated.
This study was supported by the National Natural Science Foundation of China (grant nos. 41875135 and 41675075), the National Key R&D Program of China (grant nos. 2017YFC1502202 and 2016YFA0602101), and the S&T Development Fund of CAMS (grant no. 2019KJ011).
This paper was edited by Richard Neale and reviewed by two anonymous referees.
Denis, B., Côté, J., and Laprise, R.: Spectral Decomposition of Two-Dimensional Atmospheric Fields on Limited-Area Domains Using the Discrete Cosine Transform (DCT), Mon. Weather Rev., 130, 1812–1829, https://doi.org/10.1175/1520-0493(2002)130<1812:SDOTDA>2.0.CO;2, 2002.
Du, Q., Faber, V., and Gunzburger, M.: Centroidal Voronoi tessellations: Applications and algorithms, Siam Rev., 41, 637–676, 1999.
Düben, P. D. and Korn, P.: Atmosphere and Ocean Modeling on Grids of Variable Resolution – A 2D Case Study, Mon. Weather Rev., 142, 1997–2017, https://doi.org/10.1175/mwr-d-13-00217.1, 2014.
Dubos, T., Dubey, S., Tort, M., Mittal, R., Meurdesoif, Y., and Hourdin, F.: DYNAMICO-1.0, an icosahedral hydrostatic dynamical core designed for consistency and versatility, Geosci. Model Dev., 8, 3131–3150, https://doi.org/10.5194/gmd-8-3131-2015, 2015.
Gassmann, A.: A global hexagonal C-grid non-hydrostatic dynamical core (ICON-IAP) designed for energetic consistency, Q. J. Roy. Meteor. Soc., 139, 152–175, https://doi.org/10.1002/qj.1960, 2013.
Gettelman, A., Callaghan, P., Larson, V. E., Zarzycki, C. M., Bacmeister, J. T., Lauritzen, P. H., Bogenschutz, P. A., and Neale, R. B.: Regional Climate Simulations With the Community Earth System Model, J. Adv. Model. Earth Sy., 10, 1245–1265, https://doi.org/10.1002/2017ms001227, 2018.
GRIST-Dev.: ParGRIST-A20-0705, Zenodo, https://doi.org/10.5281/zenodo.3930643, 2020a.
GRIST-Dev.: GRIST-QUVR-Mesh (Version 1.0), Zenodo, https://doi.org/10.5281/zenodo.3817060, 2020b.
Gross, M., Wan, H., Rasch, P. J., Caldwell, P. M., Williamson, D. L., Klocke, D., Jablonowski, C., Thatcher, D. R., Wood, N., Cullen, M., Beare, B., Willett, M., Lemarié, F., Blayo, E., Malardel, S., Termonia, P., Gassmann, A., Lauritzen, P. H., Johansen, H., Zarzycki, C. M., Sakaguchi, K., and Leung, R.: Physics–Dynamics Coupling in Weather, Climate, and Earth System Models: Challenges and Recent Progress, Mon. Weather Rev., 146, 3505–3544, https://doi.org/10.1175/MWR-D-17-0345.1, 2018.
Guba, O., Taylor, M. A., Ullrich, P. A., Overfelt, J. R., and Levy, M. N.: The spectral element method (SEM) on variable-resolution grids: evaluating grid sensitivity and resolution-aware numerical viscosity, Geosci. Model Dev., 7, 2803–2816, https://doi.org/10.5194/gmd-7-2803-2014, 2014.
Harris, L. M. and Lin, S.-J.: A Two-Way Nested Global-Regional Dynamical Core on the Cubed-Sphere Grid, Mon. Weather Rev., 141, 283–306, https://doi.org/10.1175/MWR-D-11-00201.1, 2012.
Harris, L. M., Lin, S.-J., and Tu, C.: High-Resolution Climate Simulations Using GFDL HiRAM with a Stretched Global Grid, J. Climate, 29, 4293–4314, https://doi.org/10.1175/JCLI-D-15-0389.1, 2016.
Hashimoto, A., Done, J. M., Fowler, L. D., and Bruyère, C. L.: Tropical cyclone activity in nested regional and global grid-refined simulations, Clim. Dynam., 47, 497–508, https://doi.org/10.1007/s00382-015-2852-2, 2015.
Hollingsworth, A., Kållberg, P., Renner, V., and Burridge, D. M.: An internal symmetric computational instability, Q. J. Roy. Meteor. Soc., 109, 417–428, https://doi.org/10.1002/qj.49710946012, 1983.
Hong, S.-Y., Noh, Y., and Dudhia, J.: A New Vertical Diffusion Package with an Explicit Treatment of Entrainment Processes, Mon. Weather Rev., 134, 2318–2341, https://doi.org/10.1175/MWR3199.1, 2006.
Hourdin, F., Musat, I., Bony, S., Braconnot, P., Codron, F., Dufresne, J.-L., Fairhead, L., Filiberti, M.-A., Friedlingstein, P., Grandpeix, J.-Y., Krinner, G., LeVan, P., Li, Z.-X., and Lott, F.: The LMDZ4 general circulation model: climate performance and sensitivity to parametrized physics with emphasis on tropical convection, Clim. Dynam., 27, 787–813, https://doi.org/10.1007/s00382-006-0158-0, 2006.
Jablonowski, C. and Williamson, D. L.: A baroclinic instability test case for atmospheric model dynamical cores, Q. J. Roy. Meteor. Soc., 132, 2943–2975, https://doi.org/10.1256/qj.06.12, 2006.
Jacobsen, D. W., Gunzburger, M., Ringler, T., Burkardt, J., and Peterson, J.: Parallel algorithms for planar and spherical Delaunay construction with an application to centroidal Voronoi tessellations, Geosci. Model Dev., 6, 1353–1365, https://doi.org/10.5194/gmd-6-1353-2013, 2013.
Ju, L., Ringler, T., and Gunzburger, M.: Voronoi tessellations and their application to climate and global modeling, in: Numerical Techniques for Global Atmospheric Models, Lecture Notes in Computational Science and Engineering, 80, edited by: Lauritzen, P. H., Jablonowski, C., Taylor, M. A., and Nair, R. D., Springer-Verlag, 313–342, 2011.
Jung, T., Miller, M. J., Palmer, T. N., Towers, P., Wedi, N., Achuthavarier, D., Adams, J. M., Altshuler, E. L., Cash, B. A., Kinter, J. L., Marx, L., Stan, C., and Hodges, K. I.: High-Resolution Global Climate Simulations with the ECMWF Model in Project Athena: Experimental Design, Model Climate, and Seasonal Forecast Skill, J. Climate, 25, 3155–3172, https://doi.org/10.1175/jcli-d-11-00265.1, 2012.
Klemp, J. B., Skamarock, W. C., and Ha, S.: Damping Acoustic Modes in Compressible Horizontally Explicit Vertically Implicit (HEVI) and Split-Explicit Time Integration Schemes, Mon. Weather Rev., 146, 1911–1923, https://doi.org/10.1175/MWR-D-17-0384.1, 2018.
Lauritzen, P. H., Jablonowski, C., Taylor, M. A., and Nair, R. D.: Rotated Versions of the Jablonowski Steady-State and Baroclinic Wave Test Cases: A Dynamical Core Intercomparison, J. Adv. Model. Earth Sy., 2, https://doi.org/10.3894/JAMES.2010.2.15, 2010.
Li, X., Zhang, Y., Peng, X., and Li, J.: Using a single column model (SGRIST1.0) for connecting model physics and dynamics in the Global-to-Regional Integrated forecast SysTem (GRIST-A20.8), Geosci. Model Dev. Discuss., https://doi.org/10.5194/gmd-2020-254, in review, 2020.
Lin, S.-J.: A “Vertically Lagrangian” Finite-Volume Dynamical Core for Global Models, Mon. Weather Rev., 132, 2293–2307, https://doi.org/10.1175/1520-0493(2004)132<2293:AVLFDC>2.0.CO;2, 2004.
Lin, Y.-L., Farley, R. D., and Orville, H. D.: Bulk Parameterization of the Snow Field in a Cloud Model, J. Clim. Appl. Meteorol., 22, 1065–1092, https://doi.org/10.1175/1520-0450(1983)022<1065:BPOTSF>2.0.CO;2, 1983.
Liu, Y. and Yang, T.: Impact of Local Grid Refinements of Spherical Centroidal Voronoi Tessellations for Global Atmospheric Models, Commun. Comput. Phys., 21, 1310–1324, https://doi.org/10.4208/cicp.050815.020916a, 2017.
Liu, Z., Zhang, Y., Huang, X., Li, J., Wang, D., Wang, M., and Huang, X.: Development and performance optimization of a parallel computing infrastructure for an unstructured-mesh modelling framework, Geosci. Model Dev. Discuss., https://doi.org/10.5194/gmd-2020-158, in review, 2020.
Mlawer, E. J., Taubman, S. J., Brown, P. D., Iacono, M. J., and Clough, S. A.: Radiative transfer for inhomogeneous atmospheres: RRTM, a validated correlated-k model for the longwave, J. Geophys. Res., 102, 16663–16682, https://doi.org/10.1029/97jd00237, 1997.
Powers, J. G., Klemp, J. B., Skamarock, W. C., Davis, C. A., Dudhia, J., Gill, D. O., Coen, J. L., Gochis, D. J., Ahmadov, R., Peckham, S. E., Grell, G. A., Michalakes, J., Trahan, S., Benjamin, S. G., Alexander, C. R., Dimego, G. J., Wang, W., Schwartz, C. S., Romine, G. S., Liu, Z., Snyder, C., Chen, F., Barlage, M. J., Yu, W., and Duda, M. G.: The Weather Research and Forecasting Model: Overview, System Efforts, and Future Directions, B. Am. Meteorol. Soc., 98, 1717–1737, https://doi.org/10.1175/BAMS-D-15-00308.1, 2017.
Reed, K. A. and Jablonowski, C.: An Analytic Vortex Initialization Technique for Idealized Tropical Cyclone Studies in AGCMs, Mon. Weather Rev., 139, 689–710, https://doi.org/10.1175/2010mwr3488.1, 2011.
Reed, K. A. and Jablonowski, C.: Idealized tropical cyclone simulations of intermediate complexity: A test case for AGCMs, J. Adv. Model. Earth Sy., 4, https://doi.org/10.1029/2011MS000099, 2012.
Reed, K. A., Bacmeister, J. T., Rosenbloom, N. A., Wehner, M. F., Bates, S. C., Lauritzen, P. H., Truesdale, J. E., and Hannay, C.: Impact of the dynamical core on the direct simulation of tropical cyclones in a high-resolution global model, Geophys. Res. Lett., 42, 3603–3608, https://doi.org/10.1002/2015GL063974, 2015.
Ringler, T., Ju, L., and Gunzburger, M.: A multiresolution method for climate system modeling: application of spherical centroidal Voronoi tessellations, Ocean Dynam., 58, 475–498, https://doi.org/10.1007/s10236-008-0157-2, 2008.
Ringler, T., Petersen, M., Higdon, R. L., Jacobsen, D., Jones, P. W., and Maltrud, M.: A multi-resolution approach to global ocean modeling, Ocean Model., 69, 211–232, https://doi.org/10.1016/j.ocemod.2013.04.010, 2013.
Ringler, T. D., Thuburn, J., Klemp, J. B., and Skamarock, W. C.: A unified approach to energy conservation and potential vorticity dynamics for arbitrarily-structured C-grids, J. Comput. Phys., 229, 3065–3090, https://doi.org/10.1016/j.jcp.2009.12.007, 2010.
Ringler, T. D., Jacobsen, D., Gunzburger, M., Ju, L., Duda, M., and Skamarock, W.: Exploring a Multiresolution Modeling Approach within the Shallow-Water Equations, Mon. Weather Rev., 139, 3348–3368, https://doi.org/10.1175/mwr-d-10-05049.1, 2011.
Sakaguchi, K., Leung, L. R., Zhao, C., Yang, Q., Lu, J., Hagos, S., Rauscher, S. A., Dong, L., Ringler, T. D., and Lauritzen, P. H.: Exploring a Multiresolution Approach Using AMIP Simulations, J. Climate, 28, 5549–5574, https://doi.org/10.1175/JCLI-D-14-00729.1, 2015.
Skamarock, W.: NGGPS/DTG briefing on MPAS configuration options, available at: https://www.weather.gov/media/sti/nggps/20160122_MPAS_configuration_overview.pdf (last access: 7 December 2020), 2016.
Skamarock, W. C. and Gassmann, A.: Conservative Transport Schemes for Spherical Geodesic Grids: High-Order Flux Operators for ODE-Based Time Integration, Mon. Weather Rev., 139, 2962–2975, https://doi.org/10.1175/MWR-D-10-05056.1, 2011.
Skamarock, W. C., Klemp, J. B., Duda, M. G., Fowler, L. D., Park, S.-H., and Ringler, T. D.: A Multiscale Nonhydrostatic Atmospheric Model Using Centroidal Voronoi Tesselations and C-Grid Staggering, Mon. Weather Rev., 140, 3090–3105, https://doi.org/10.1175/MWR-D-11-00215.1, 2012.
Skamarock, W. C., Duda, M. G., Ha, S., and Park, S.-H.: Limited-Area Atmospheric Modeling Using an Unstructured Mesh, Mon. Weather Rev., 146, 3445–3460, https://doi.org/10.1175/MWR-D-18-0155.1, 2018.
Smagorinsky, J.: General Circulation Experiments with the Primitive Equations, Mon. Weather Rev., 91, 99–164, https://doi.org/10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2, 1963.
Stevens, B., Satoh, M., Auger, L., Biercamp, J., Bretherton, C. S., Chen, X., Düben, P., Judt, F., Khairoutdinov, M., Klocke, D., Kodama, C., Kornblueh, L., Lin, S.-J., Neumann, P., Putman, W. M., Röber, N., Shibuya, R., Vanniere, B., Vidale, P. L., Wedi, N., and Zhou, L.: DYAMOND: the DYnamics of the Atmospheric general circulation Modeled On Non-hydrostatic Domains, Progr. Earth Planet. Sci., 6, 61, https://doi.org/10.1186/s40645-019-0304-z, 2019.
Taylor, M. A.: Conservation of Mass and Energy for the Moist Atmospheric Primitive Equations on Unstructured Grids, in: Numerical Techniques for Global Atmospheric Models, edited by: Lauritzen, P., Jablonowski, C., Taylor, M., and Nair, R., Springer, Berlin, Heidelberg, 357–380, 2011.
Thuburn, J., Ringler, T. D., Skamarock, W. C., and Klemp, J. B.: Numerical representation of geostrophic modes on arbitrarily structured C-grids, J. Comput. Phys., 228, 8321–8335, https://doi.org/10.1016/j.jcp.2009.08.006, 2009.
Tiedtke, M.: A Comprehensive Mass Flux Scheme for Cumulus Parameterization in Large-Scale Models, Mon. Weather Rev., 117, 1779–1800, https://doi.org/10.1175/1520-0493(1989)117<1779:acmfsf>2.0.co;2, 1989.
Ullrich, P. A. and Jablonowski, C.: An analysis of 1D finite-volume methods for geophysical problems on refined grids, J. Comput. Phys., 230, 706–725, https://doi.org/10.1016/j.jcp.2010.10.014, 2011.
Ullrich, P. A., Jablonowski, C., Kent, J., Lauritzen, P. H., Nair, R., Reed, K. A., Zarzycki, C. M., Hall, D. M., Dazlich, D., Heikes, R., Konor, C., Randall, D., Dubos, T., Meurdesoif, Y., Chen, X., Harris, L., Kühnlein, C., Lee, V., Qaddouri, A., Girard, C., Giorgetta, M., Reinert, D., Klemp, J., Park, S.-H., Skamarock, W., Miura, H., Ohno, T., Yoshida, R., Walko, R., Reinecke, A., and Viner, K.: DCMIP2016: a review of non-hydrostatic dynamical core design and intercomparison of participating models, Geosci. Model Dev., 10, 4477–4509, https://doi.org/10.5194/gmd-10-4477-2017, 2017.
Wang, L., Zhang, Y., Li, J., Liu, Z., and Zhou, Y.: Understanding the Performance of an Unstructured-Mesh Global Shallow Water Model on Kinetic Energy Spectra and Nonlinear Vorticity Dynamics, J. Meteorol. Res.-PRC, 33, 1075–1097, https://doi.org/10.1007/s13351-019-9004-2, 2019.
Wehner, M. F., Reed, K. A., Li, F., Prabhat, Bacmeister, J., Chen, C.-T., Paciorek, C., Gleckler, P. J., Sperber, K. R., Collins, W. D., Gettelman, A., and Jablonowski, C.: The effect of horizontal resolution on simulation quality in the Community Atmospheric Model, CAM5.1, J. Adv. Model. Earth Sy., 6, 980–997, https://doi.org/10.1002/2013ms000276, 2014.
Weller, H., Weller, H. G., and Fournier, A.: Voronoi, Delaunay, and Block-Structured Mesh Refinement for Solution of the Shallow-Water Equations on the Sphere, Mon. Weather Rev., 137, 4208–4224, https://doi.org/10.1175/2009mwr2917.1, 2009.
Yu, R.: A Two-Step Shape-Preserving Advection Scheme, Adv. Atmos. Sci., 11, 479–490, 1994.
Yu, R., Zhang, Y., Wang, J., Li, J., Chen, H., Gong, J., and Chen, J.: Recent Progress in Numerical Atmospheric Modeling in China, Adv. Atmos. Sci., 36, 938–960, https://doi.org/10.1007/s00376-019-8203-1, 2019.
Zarzycki, C. M. and Jablonowski, C.: A multidecadal simulation of Atlantic tropical cyclones using a variable-resolution global atmospheric general circulation model, J. Adv. Model. Earth Sy., 6, 805–828, https://doi.org/10.1002/2014ms000352, 2014.
Zarzycki, C. M., Jablonowski, C., and Taylor, M. A.: Using Variable-Resolution Meshes to Model Tropical Cyclones in the Community Atmosphere Model, Mon. Weather Rev., 142, 1221–1239, https://doi.org/10.1175/mwr-d-13-00179.1, 2014.
Zhang, Y.: Extending High-Order Flux Operators on Spherical Icosahedral Grids and Their Applications in the Framework of a Shallow Water Model, J. Adv. Model. Earth Sy., 10, 145–164, https://doi.org/10.1002/2017MS001088, 2018.
Zhang, Y., Chen, H., and Yu, R.: Simulations of Stratus Clouds over Eastern China in CAM5: Sensitivity to Horizontal Resolution, J. Climate, 27, 7033–7052, https://doi.org/10.1175/jcli-d-13-00732.1, 2014.
Zhang, Y., Chen, H., and Yu, R.: Simulations of Stratus Clouds over Eastern China in CAM5: Sources of Errors, J. Climate, 28, 36–55, https://doi.org/10.1175/JCLI-D-14-00350.1, 2015.
Zhang, Y., Yu, R., and Li, J.: Implementation of a conservative two-step shape-preserving advection scheme on a spherical icosahedral hexagonal geodesic grid, Adv. Atmos. Sci., 34, 411–427, https://doi.org/10.1007/s00376-016-6097-8, 2017.
Zhang, Y., Li, J., Yu, R., Zhang, S., Liu, Z., Huang, J., and Zhou, Y.: A Layer-Averaged Nonhydrostatic Dynamical Framework on an Unstructured Mesh for Global and Regional Atmospheric Modeling: Model Description, Baseline Evaluation, and Sensitivity Exploration, J. Adv. Model. Earth Sy., 11, 1685–1714, https://doi.org/10.1029/2018ms001539, 2019.
Zhang, Y., Li, J., Yu, R., Liu, Z., Zhou, Y., Li, X., and Huang, X.: A Multiscale Dynamical Model in a Dry-mass Coordinate for Weather and Climate Modeling: Moist Dynamics and its Coupling to Physics, Mon. Weather Rev., 148, 2671–2699, https://doi.org/10.1175/mwr-d-19-0305.1, 2020.
Zhao, M., Held, I. M., and Lin, S.-J.: Some Counterintuitive Dependencies of Tropical Cyclone Frequency on Parameters in a GCM, J. Atmos. Sci., 69, 2272–2283, https://doi.org/10.1175/JAS-D-11-0238.1, 2012.
We only consider the issues related to model dynamics in this study.
In the context of GRIST, dycore specifically refers to the dry part of the governing equations excluding tracer transport; this should be distinguished from the typical usage of dycore in the literature.
Using two halo layers may be possible but would require a more complicated communication rule than the current one, which is undesirable.
To achieve consistency with earlier QU tests.
Based on some tests, further halving the coefficient for a given resolution is also acceptable; see Fig. 4.
Note that in Z19, due to an incorrect display setting, there is a plotting mistake in the kinetic energy spectra (their Fig. 10): the tick marks of the entire top x axis (representing wavelength) do not correspond to the actual wavelength of the data. The bottom x axis is correct.