Configuration and evaluation of a global unstructured mesh atmospheric model (GRIST-A20.9) based on the variable-resolution approach

Targeting a long-term effort towards a variableresolution (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 fineresolution 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. Published by Copernicus Publications on behalf of the European Geosciences Union. 6326 Y. Zhou et al.: GRIST-A20.9

environment, the VR model can have the tropical cyclone stably pass the transition zone in various configurations. A 23 series of sensitivity tests examines the model performance in a hierarchical refinement mode. The simulations exhibit 24 consistency even when the VR mesh is slightly perturbed by one of the three parameters that control the density function. 25 The tropical cyclone, starting from the 2 nd -refinement region and passing through the inner transition zone, gets 26 intensified and covers a smaller area in the refined regions. Such variations are consistent with the behavior that one 27 may observe when uniformly refining the QU mesh. In the full physics environment with a highly variable mesh that 28 reaches sub-10-kilometer resolution, the VR model also produces a reasonable evolution for the tropical cyclone. The 29 explicit diffusion shows its usefulness in terms of suppressing some unrealistic isolated-scale structures that are far 30 away from the initial vortex, and does not adversely affect the physically important object. The fine-scale structure is 31 determined mainly by the fine-resolution area, although the systems may have larger differences before they move into 32 the fine-resolution area. Altogether, this work demonstrates that the multiresolution configuration is a reliable and 33 economic alternative to high-resolution global modeling. The adverse impact due to mesh transition and the coarse 34 region can be well controlled. 35 36

Introduction 37
Increasing resolution is generally regarded as an effective way to improve global weather and climate modeling 38 (Jung et al. 2012;Wehner et al. 2014;Zhang et al. 2014;Yu et al. 2019). It is apparent that more computational and 39 storage resources are required for higher resolution models. This leads to a major challenge for efficient model 40 development and application. The emergence of the locally refined, variable-resolution (VR) modeling approach offers 41 a complementary route. The term VR is a broad concept. It may be realized with different styles, such as nested regional 42 modeling with multiple grids, abrupt nonconforming mesh division, stretched grids, and the multiresolution approach. 43 The stretched grid (e.g., Hourdin et al. 2006; Harris et al. 2016) and the multiresolution approaches (e.g., Ringler et al. 44 2011; Guba et al. 2014) are close in terms of their conforming style. They maintain the global modeling configuration 45 while permitting increased resolution for certain regions. The multiresolution approach is usually realized by an 46 unstructured mesh model, such that a more flexible resolution choice can be achieved by considering multiple regions. 47 Such a global-to-regional approach is the VR style that will be investigated in this study. 48 Numerical weather and climate modeling has shown that the VR approach can preserve the benefits of high-49 resolution applications for certain regions at a lower computational cost, as the total number of grid points can be greatly 50 reduced (e.g., Sakaguchi et al. 2015;Skamarock et al. 2018;Gettelman et al. 2018). This advantage is especially 51 valuable for high-resolution modeling that may reach the convection-permitting regime. While the global cloud/storm-52 resolving (a.k.a., convection-permitting/allowing) modeling approach has been widely adopted (e.g., Stevens et al. 53 2019), it is still expensive and inefficient to frequently run such models for routine model development, research and 54 application. The VR model provides an efficient testbed for evaluating global model configurations and testing scale-55 aware physics. It offers flexible resolution configurations that may depend on physical interests. When properly 56 formulated, it can be an intermediate and transitional step before establishing global convection-permitting modeling. 57 While the VR approach has shown some benefits, it may potentially suffer from some problems 1 . The nonuniform 58 mesh, though it can be gradually refined, hardly decreases the global errors (as compared to its uniform counterpart 59 that has the same degrees of freedom) because the truncation error is controlled mainly by the coarse-resolution region 60 temperature is 29℃ globally. These processes are coupled in a time-splitting manner within the package, and the 162 package is coupled to GRIST in a pure operator-splitting approach (ptend_f2_sudden; see Z20). 163 (ii) A climate physics package adopted from CAM5. This package is not used in this study, but the details can be 164 found in Li et al. (2020). It is currently being tuned for the HDC that targets long-term climate modeling, and can also 165 be used for short-term integration in a weather forecast mode (e.g., Zhang et al. 2015). 166 (iii) For GRIST-NDC, which is targeted at simulating nonhydrostatic dynamics in a nonhydrostatic regime, a set 167 of parameterization schemes from the Weather Research and Forecast (WRF) model (Powers et al. 2017) has been 168 implemented, and is used in this study. The detailed schemes include: a six-species cloud microphysics scheme (Lin et 169 al. 1983) from WRF version 2.0; the Tiedtke cumulus scheme (Tiedtke 1989) from WRF version 3.7.1; the YSU (Yonsei 170 University) planetary boundary layer scheme (Hong et al. 2006) and a surface scheme from WRF version 2.0; the 171 longwave radiation scheme is the RRTM (Rapid Radiative Transfer Model) from WRF version 2.0 (Mlawer et al. 1997); 172 and the shortwave scheme is a CAM radiation module from WRF version 3.4.1. 173 The internal coupling of this package is process splitting. All processes start from the dynamics-updated state and 174 send back their respective tendencies to the dynamical model without modifications of the physics state variables. The 175 exception is that microphysics will update the local physics state variables, so the calling sequence still matters. The 176 physics-dynamics coupling uses a hybrid approach that combines the tendency method (ptend_rk) and the operator-177 splitting approach, as described in Z20. In particular, radiation heating is carried over the internal integration of the 178 dycore as in a tendency method, and other tendencies (microphysics, boundary layer, cumulus) are updated as in a pure 179 operator-splitting approach (ptend_f2_sudden in Z20). We emphasize that this package (including the choice of 180 different schemes and its internal coupling) is still experimental and preliminary. It has not been comprehensively tuned. 181 Ongoing tests and modifications are required to refine the performance. In this study, it is only intended to evaluate the 182 particular performance of VR modeling and the role of explicit diffusion in a full physics scenario. 183

Time-step choices and explicit diffusion 184
For a VR model, the time step is theoretically restricted by its fine-resolution region. As emphasized by one reviewer 185 and based on our own experience, some improper VR configurations may lead to higher stability restriction than their 186 equivalent QU counterparts. For example, with regard to the acoustic-mode filter, Klemp et al. (2018) suggested that 187 their original formulation is more problematic in VR applications, and cannot effectively remove the acoustic noise. 188 Uncontrolled acoustic modes may artificially accumulate energy and thus impose a higher stability restriction. For the 189 configurations and tests that we have examined, a suitable configuration of explicit diffusion is mostly relevant to this 190 issue. The details will be given in the following section. 191 The explicit diffusion tendencies are generated at the largest time interval of each model step (i.e., physics step). 192 They are coupled to model dynamics as in a tendency method (ptend_rk). Z20 showed that the tendency method has a 193 slightly higher stability restriction than the pure operator-splitting approach (ptend_f2_sudden), but it can benefit from 194 the higher accuracy of the time integrator. The tendency method does not require additional data communication. No 195 explicit diffusion is activated for tracer transport in this paper. 196

Smagorinsky diffusion 197
Flow-dependent Smagorinsky diffusion (Smagorinsky 1963) is used for velocities and potential temperature, 198 following previous QU model tests. It is activated in all experiments, except the baroclinic wave test at the quasi-199 uniform G6 resolution. This scheme uses a second-order Laplacian operator multiplied by a flow-dependent eddy 200 viscosity. The eddy viscosity is defined at the edge point. The Smagorinsky diffusion is not very scale selective, but its 201 flow-dependent feature makes it selective in terms of where and how much to diffuse (see e.g., Fig. 9 in Gassmann 202 2013). The diffusion strength is acceptable overall, as evidenced by the sharp gradient in the QU baroclinic wave test. 203 a stronger diffusion (than typically required) for the refinement region. This also leads to a higher stability restriction, 205 especially when the refinement ratio becomes large. In this version, the square of this mean length scale is replaced by 206 the local length product of a pair of crossing edges. In the tropical cyclone test (simple physics), this local approach 207 increases the maximum wind magnitude in the eyewall (as compared to the unscaled version in the preprint), because 208 a smaller amount of diffusion is imposed on the refinement area (see Section 5.1). It also reduces the parametric 209 sensitivity to the Smagorinsky coefficient as found in the preprint. When varying the coefficient, the present version 210 generates more consistent solutions, similar to a QU model. This local scaling approach was used for the VR mode but 211 not for the QU mode. 4 212

Hyperdiffusion 213
The hyperdiffusion option was not used in the initial preprint. Based on some exploration and experiments, we have 214 found that activating scale-selective fourth-order hyperdiffusion for the horizontal velocity shows demonstrable added 215 value for VR modeling. First, in the baroclinic wave test, in the absence of hyperdiffusion, a higher Smagorinsky 216 coefficient (compared to the equivalent QU test) is required to suppress grid-scale noise due to the mesh refinement, 217 which, in turn, restricts the numerical stability. When the hyperdiffusion is activated, we can use a moderate 218 Smagorinsky coefficient that does not challenge the numerical stability, and maintains a quality solution. The DTP 219 (dycore-tracer-physics) splitting mode benefits from this most, because diffusion is called at the step of the model 220 physics. Second, even with a higher coefficient, the Smagorinsky diffusion is inactive over certain regions where weak 221 flow deformation dominates. In the baroclinic wave test, some grid-scale oscillations are more conspicuous over these 222 regions (see preprint). Such noise is due to the mesh transition, and some of it is akin to the Hollingsworth instability. 223 As will be shown, activating a background hyperdiffusion for the horizontal wind successfully removes such noise. The 224 solutions are overall less oscillatory than those in the preprint. Third, in the tropical cyclone tests with full physics, 225 which presents more nonlinear feedback, the hyperdiffusion option is also effective in suppressing some isolated-scale 226 structures. Unlike the minor disturbances that may be generated near the major tropical cyclone in the simple physics 227 test, these systems are far away from the initial vortex, and thus highly unrealistic. For these reasons, the hyperdiffusion 228 option is used for the VR mode, but not for the QU test. 229 The hyperdiffusion operator is formulated by recursively using the Laplacian operator (see Z19 for details). The 230 diffusion coefficient can be determined in a relatively empirical way. For VR modeling, we adopt the approach 231 documented in Zarzycki et al. (2014) for scaling the coefficient: 232 (1) The reference length ∆ "#$ and reference viscosity coefficient ! (∆ "#$ ) are empirically determined. This 233 formulation reduces ! (∆ ) by a factor of 10 for every halving of resolution. A similar scaling approach is also used 234 in MPAS (Skamarock 2016). We typically use the configuration for a G-level resolution that is close to the finest 235 resolution on the mesh. Some typical values for GRIST are documented in Table 1. These values 5 are smaller than 236 those in Zarzycki et al. (2014) for the corresponding length scale by a factor of 5. The local grid distance (∆ ) is an 237 average distance between the grid point and all its nearest neighbors, that is, a cell-based value. The edge-based value 238 used for hyperdiffusion is an average of two neighboring cell values. two key elements of the SCVT: generators and density function. A spherical Voronoi tessellation is a spatial subdivision 242 4 To achieve consistency with earlier QU tests. 5 Based on some tests, further halving the coefficient for a given resolution is also acceptable; see where ‖•‖ denotes the geodesic distance. Each point , is called a generator and its corresponding Voronoi region , 245 is called the Voronoi cell. A spherical Voronoi tessellation becomes an SCVT when the generators are also the centroids 246 of the Voronoi cells. In this study, the SCVT mesh is constructed by an iterative process based on Lloyd's algorithm 247 (Du et al. 1999). In particular, a parallel algorithm is used  to avoid time-consuming serial 248 construction. That said, generating a quality VR SCVT is still nontrivial (but only done once). In our implementation, 249 the iteration stops when two criteria are satisfied: (i) it reaches an empirically determined minimum step; and (ii) the 250 circumcenter of each triangle falls within its shape. 251

Generators 252
Three original generators are used in this study: 253 Icosahedron bisection. This approach benefits from the excellent uniform properties due to bisections of a 254 regular icosahedron. The mesh resolution is referred to as G-level/G , where denotes the number of 255 bisections. After each bisection, the total grid number is approximately four times greater than the previous 256 one. 257 Icosahedron bisection with a final-step trisection. Instead of using recursive bisection, a trisection is used 258 at the final step, to achieve an intermediate resolution between two neighboring G-level resolutions. This mesh 259 is referred to as G B3 (i.e., bisections plus one trisection). For instance, the resolution of G5B3 (~80 km) 260 is between that of G6 (~120 km) and that of G7 (~60 km). The number of added primal cells (mainly hexagons) 261 from G5 to G5B3 is equal to four times the number of dual cells (triangles) in G5. 262 Spherical uniform random (SUR) set of points. Initial points are created uniformly on the sphere by the 263 Monte Carlo method. In this way, the original generators can be obtained by using an arbitrary number of 264 points, not restricted to a sub-divided icosahedron. 265

Density function 266
By specifying the density function, the SCVT is able to precisely control the distribution of the local resolution. 267 For any two Voronoi regions indexed by and , the conjecture is: 268 where ( , ) is the density function evaluated at , and , measures the local mesh resolution. This relation is valid 269 in a theoretical sense. 270 A QU mesh can be constructed when the density is one on the sphere, and the Voronoi regions are approximately 271 equivalent to each other. For the VR mesh, the basic density function used in this study is: 272 where ‖ "9 − , ‖ denotes the geodesic distance between the location of the refinement center and each generator. 273 indicates the width of the transition zone between the fine-resolution and coarse-resolution regions; defines the 274 coverage radius of the fine-resolution region; measures the densification ratio between the finest and coarsest 275 resolutions. A sample X4 mesh based on this density function with = (1/4) ! is shown in Fig. 1a. 276 Because the basic density function is fixed to a single-region refinement, we adjust it for multi-region refinement. 277 The multi-region refinement is divided into two styles based on the refinement centers. In the hierarchical refinement 278 way (Fig. 1b), we add a uniform intermediate-resolution region between the inner fine-resolution and the outer coarse-279 resolution regions (see Fig. 2). This helps to avoid an overly high densification ratio between two neighboring regions. 280 Eq. (4) can be generalized to a form that allows us to control the resolution of the intermediate region: 281 is designed to control the resolution of the intermediate-resolution region, also referred to as the 2 nd -refinement region 282 ( " ( ) that is located between the 1 st -refinement ( " ' ) and the coarse-resolution regions ( 9 ): 283 Generally, ∈ [ , 1]. The function of is similar to that in the previous single-region refinement, except that the fine-284 resolution region is referred to as the 1 st -refinement region here: 285 Corresponding to , we refer to the meshes that are generated based on values of (1) ! , (1/2) ! , and (1/3) ! as 286 XL1, XL2, and XL3 meshes, since the resolutions of the 1 st -refinement and the 2 nd -refinement regions vary according 287 to the inner densification ratios of 1, 2, and 3, respectively. For example, when is fixed at X4, the hierarchical meshes 288 based on G6 are called G6X4L1, G6X4L2, and G6X4L3 meshes. 289 In a polycentric refinement mode (Fig. 1c), by adding a different refinement center "9) , the density function of 290 the polycentric refinement mode is defined as: 291 The geodesic distance between the two refinement centers should satisfy ‖ "9* − "9) ‖ > 2 . 292

Single-region refinement 294
The dry-atmosphere test examines the pure numerical solution of the model. It does not include the nonlinear 295 interaction between dynamics, moisture transport, and parameterization. Based on the JW06 baroclinic wave test, we 296 first compare the VR and QU simulations. Previous studies that employed this test for a VR model include Gettelman 300,000 for G6X4, and 1,000,000 for G5B3X4 and SURX4. The refinement center is placed at 35°N, 180°E, with = 302 /20 and = /6. The detailed model configuration is given in Table S1. The timestep of the VR model is limited 303 by its fine-resolution region and is set accordingly, although the currently used timesteps do not represent the maximum 304 allowable step. 305 We first examine the resolving ability of the fine-scale structures. Figure 3 shows the relative vorticity field at the 306 model level nearest to 850 hPa (level 24) after 10 days. These values have been remapped to the Voronoi cell from the 307 raw triangular grid, so as to avoid the aliasing of certain oscillatory patterns. Such patterns (see preprint) actually reflect 308 the mesh shape, and are more conspicuous for coarse resolution. Figure S1 further shows the QU model results 309 interpolated to the regular longitude-latitude grid as a reference. As the resolution increases, the QU models simulate 310 stronger vortices with a clear filament structure (Fig. 3a-3c). The VR model can simulate the smooth structure of the 311 waves in the refined region, as in the high-resolution QU model. In the VR mode, two vortices in the west fall within 312 the fine-resolution region. The structure of the westernmost vortex in G6X4 (Fig. 3d) is close to G7. G5B3X4 and 313 SURX4 (Fig. 3e, 3f) further produce finer scale structures than G6X4, closer to G8. The easternmost wave falls within 314 the transition zone in three VR runs. The variation in the mesh sizes there does not distort the wave pattern, and the 315 fine-scale structure can be improved as the resolution of the transition zone increases (e.g., from G6X4 to G5B3X4). 316 In the SURX4 test, a minor roughness is found on the tails of the vortices. This deficiency is largely due to local mesh 317 irregularities. Compared to the icosahedron-based SCVT generators, the random generators tend to slightly degrade the 318 To present a more quantitative evaluation of the momentum field, we examine regional kinetic energy spectra 6 320 over the refinement area (red box in Fig. 3). We use the discrete cosine transform to perform this analysis (cf. Denis et 321 al. 2002). The computational procedure is briefly documented in the Appendix. The decomposition is made for relative 322 vorticity and divergence. Fig. 4a-4c show the results from different tests. It is clear that the rotational mode dominates 323 the kinetic energy. The hydrostatic model (Fig. S2) produces similar results to the nonhydrostatic model; thus only the 324 nonhydrostatic core is discussed in the main text. 325 At the 1 st wavenumber, runs with different resolutions show larger discrepancies. This is because this lowest mode 326 absorbs most of the large-scale trend, and corresponds to only a half-cosine wave. Denis et al. (2002) suggested that 327 one may remove this mode if desired, but we keep it here. From the 2 nd to the 10 th wavenumber, all the curves are 328 consistent overall, suggesting that the well-resolved structures are robust to various tests. The major difference lies in 329 the high wavenumbers. For the VR results, G6X4 is close to G7. This is expected because they have similar resolution 330 over the selected domain (~60 km). For the same reason, G5B3X4 and SURX4 produce better spectral tails than G7 331 and G6X4, closer to G8. This confirms that increasing the fine resolution of the VR model is able to improve its fine-332 scale resolving ability. At such high wavenumbers, SURX4 is even closer to G8 as it has slightly higher energy in the 333 tail. However, this actually reflects that SURX4 has slightly more grid-scale oscillations than G5B3X4. Therefore, an 334 examination of kinetic energy spectra should be accompanied by a close look at the real field. 335 In the context of kinetic energy spectra, it is useful to demonstrate the impact of the hyperviscosity coefficient. On 336 the basis of G5B3X4, we vary the reference hyperviscosity coefficient under a fixed reference length (30 km). Note 337 that 2 × 10 *) is the default configuration for this test. As shown in Fig. 4d-4f, when using 2 × 10 *' , spectra are 338 seriously damped at the high wavenumbers in both the rotational and divergent components. A flat tail is generated. 339 This indicates that the diffusion strength is overly strong for the fine-resolution area. The other four tests generally 340 produce consistent results: reducing the coefficient tends to slightly uplift the tail, that is, increase the kinetic energy. 341 While the lowest coefficient (2 × 10 ** ) seems to produce a nicer tail than the default run, this, again, reflects the fact 342 that slightly more small-scale oscillations are generated and the solutions are less clean (but still acceptable). Therefore, 343 when tuning a VR model, one should achieve a minimally required hyperviscosity that neither significantly damps the 344 field nor becomes unable to suppress certain grid-scale oscillations. For this test, 2 × 10 *) is fairly close to the optimal 345 choice. 346 Next, we assess the solution errors by first examining the surface pressure field. Its global distribution on day 9 347 simulated by GRIST-NDC is shown in Fig. S3. The locations and magnitudes of the high-and low-pressure centers are 348 consistent overall in the QU (Fig. S3a-3c) and VR (Fig. S3d-3f) runs. There are some nonzero wave patterns in all VR 349 simulations over the Southern Hemisphere, reflecting that the nonuniform grid structure leads to higher truncation 350 errors. G5B3X4 has smaller wave patterns than G6X4. This is expected as G5B3X4 has a higher resolution than G6X4 351 there. Moreover, G5B3X4 is also better than SURX4 with regard to these wave patterns. Considering that SURX4 has 352 more mesh cells, this suggests that generators based on a regular icosahedron produce higher mesh quality than random 353 generators, given roughly the same number of iterations. 354 Figure 5 shows a quantitative comparison. The global ) error norm of each test is computed against the highest 355 resolution (G8). The errors of three VR meshes (G6X4, G5B3X4, and SURX4) are higher overall than those of the QU 356 meshes. Note that these errors are also slightly higher than those in the preprint because a larger time step is used in 357 this version (we performed additional tests to check this sensitivity; figure not shown). The nonhydrostatic solver 358 produces slightly smaller errors using three VR meshes, but the overall accuracy of the two solvers is comparable. For 359 three VR meshes, G5B3X4 shows the smallest error, and is generally close to the results of G6 during the first 10 days. 360 G6X4 shows the largest error and SURX4 lies in between them. Again, the fact that SURX4 is less accurate than 361 G5B3X4 implies that random generators are more likely to be trapped into the local area during the iterative procedure 362 of mesh generation, leading to a higher degree of local irregularity. 363 For GRIST-NDC, we further tested G7X4 and G8X4. As shown (Fig. 5b), G7X4 produces smaller errors overall 364 than G5B3X4, although the difference is small. G8X4 produces higher errors than G7X4 during the first four days. We 365 suspect this is because certain initial imbalance becomes more active in a high-resolution VR configuration. Such 366 imbalance can be caused by the discrete initialization, for example, the continuous properties imposed by the analytic 367 wind field will not be exactly satisfied by the discrete normal velocity, unless the velocity is obtained based on some 368 constraints. After day 5, G8X4 produces smaller errors than G7X4, and the increment in error reduction is close to the 369 difference between G6X4 and G7X4. G8X4 also produces clearly smaller errors than G6 from day 5 to day 12. After 370 day 13, G8X4 has higher errors than G6, although its coarsest resolution is finer than that of G6. These results suggest 371 that the VR models indeed increase the global errors, but the errors can be reduced by increasing the coarse resolution 372 of the mesh. This is not surprising but a reconfirmation of the conclusion in Ringler et al. (2011), using a 3D atmospheric 373 dynamical core. 374 The JW06 test was originally proposed in the era that models based on a regular latitude-longitude or Gaussian 375 grid are still popular. For a quasi-uniform grid, it is well known that the steady state in the Southern Hemisphere (or in 376 an unperturbed condition) cannot be perfectly maintained due to mesh irregularity (Lauritzen et al. 2010). In a baroclinic 377 environment, such errors will grow exponentially to break the steady state. Because mesh irregularity increases in a VR 378 mode, the inability to maintain the steady state can be further amplified, ultimately contributing to the increased errors. 379 To reduce the impact of this issue, we add another perturbation over the Southern Hemisphere to excite two wave 380 trains, as was done in some earlier studies (e.g., Gassmann 2013). Figure 6 shows the relative vorticity field on day 10 381 from three selected tests: G6, G6X4, and G5B3X4. The sign of the relative vorticity in the Southern Hemisphere is 382 flipped to facilitate a visual comparison. The QU model (Fig. 6a) produces a comparable wave train in each hemisphere. 383 The wave trains are not exactly symmetric because the mesh cells are not exactly symmetric across the equator. In the 384 Northern Hemisphere, G6X4 and G5B3X4 ( Fig. 6b and 6c) produce similar solutions to Fig. 3d and 3e. In the Southern 385 Hemisphere, the solution of G5B3X4 is closer to that of G6 than G6X4, in terms of wave pattern and magnitude. This 386 is to be expected as G5B3X4 has higher resolution in the Southern Hemisphere. Although G6X4 cannot simulate the 387 structure in the Southern Hemisphere as G5B3X4, no serious problem is found for that region. 388 Figure 7 presents a quantitative estimation of the solution error. Generally, the magnitude and the growth of the 389 errors are close to those in Fig. 5b. A notable difference is that G8X4 produces smaller relative errors, closer to those 390 of G7 from day 8 to day 12. On day 15, G8X4 still shows comparable errors to G6. The error reduction from G6X4 to 391 G7X4 also becomes larger than that in Fig. 5b, when the waves become more developed (e.g., day 8 to 11). This implies 392 that the inability to maintain a steady state in the Southern Hemisphere indeed worsens the error estimation for a VR 393 model. 394

Multi-region refinement 395
To increase the application range of the single-region refinement, two multi-region refinement modes are examined 396 to obtain desired resolutions in multiple regions. This represents a unique aspect of the multiresolution approach. The 397 first way is a hierarchical style with one refinement center. It contains three consecutive uniform sub-regions outside 398 the refinement center: the 1 st -refinement region, the 2 nd -refinement region, and the coarse-resolution region. The 2 nd -399 refinement region provides an intermediate resolution between the 1 st -refinement and coarse-resolution regions. 400 Compared to the single-region refinement, this 2 nd -refinment region provides more uniform resolution between the 401 The symmetrical perturbation test was performed using a G8X4L2 mesh. On day 10, the first wave (the strongest) 403 over the Northern Hemisphere is experiencing a higher gradient of resolution than that over the Southern Hemisphere. 404 The fine-scale structures are well captured by the VR model (Fig. 8a). The second and third waves in the Northern 405 Hemisphere have stronger magnitudes than their equivalents in the Southern Hemisphere. The difference in each wave 406 train generally reflects the differences in the local resolution. 407 The second way uses a polycentric style with multiple refinement centers. The same double-baroclinic wave trains 408 are generated using a G7X4 mesh, with two refinement centers at 35°N, 180°E and 35°S, 180°E, respectively. On day 409 10, the model well simulates the fine-scale structures over two selected regions (Fig. 8b). The transitions between the 410 fine and coarse regions are smooth and stable. This refinement mode may provide an effective way to simultaneously 411 improve the simulations over different regions. The wave train in the Southern Hemisphere coincides well with that in 412 the North. Again, they are not exactly identical because the mesh cells are not exactly symmetric across the equator. 413 In this test case, it is useful to demonstrate the impact of activating the hyperdiffusion option. As shown in Fig. 8c  414 and 8d, when hyperdiffusion is turned off, the simulation results become rather oscillatory. The noise pattern at the tail 415 of the wave is akin to the Hollingsworth instability, and is amplified by a VR model due to mesh transition. The noise 416 is more conspicuous in the 1 st -refinement region of G8X4L2, due to higher resolution and rapid mesh transition. Also 417 note that these results are more oscillatory than those in the preprint, because the Smagorinsky diffusion in that version 418 is stronger. As mentioned in Section 2.3, even so, the Smagorinsky option is still inactive over certain regions. Using a 419 background hyperdiffusion is a good complement to this deficiency. 420

Simple physics 422
Moist-atmosphere modeling includes the nonlinear interaction of dynamics, moisture transport, and model physics. 423 We use the idealized tropical cyclone test (Reed and Jablonowski 2011) because of its clear resolution sensitivity. This 424 test is useful to examine the VR performance because the solution does not fully converge even at 10 km (Z20). The 425 simulated tropical cyclone is very sensitive to the mesh size. A higher resolution model will produce more intense 426 storms. An initial vortex is placed 10°N, 180°E with no background flow. The tropical cyclone moves northwestward 427 due to beta drift. 428 Previous studies have shown that model dynamics has a clear impact on the simulations of tropical cyclones. For 429 example, Zhao et al. (2012) showed that increasing the strength of 2D divergence damping in the finite-volume 430 dynamical core (Lin 2004) leads to more occurrences of tropical cyclones. Reed et al. (2015) showed that, in CAM5, 431 the spectral element core produces stronger tropical cyclones than the finite-volume core, when the parameterization 432 suite remains almost unchanged. For GRIST, a known sensitivity is that different tracer transport options may lead to 433 different wind magnitudes in the eyewall, as the full pressure gradient term is, by design, related to tracer mixing ratios 434 (Z20). 435 We first examine the model behaviors under a simple-physics environment (Reed and Jablonowski 2012). In the 436 preprint, we compared two representative groups of numerical tests: one group based on the NDC with DTP splitting 437 enabled; and one based on the HDC with no DTP splitting (i.e., dycore, tracer transport, and physics use the same time 438 step). Results from these two groups are consistent overall. The impacts of DTP splitting and hydrostatic/nonhydrostatic 439 options do not generate discernable differences across various VR meshes. This is consistent with the previous QU 440 model tests (Z19; Z20) in that: (i) the nonhydrostatic solver behaves similarly to its hydrostatic counterpart under the 441 hydrostatic regime; and (ii) the DTP splitting does not cause a degeneration in performance when it is properly 442 configured. In this version, which uses an updated configuration, only the NDC is tested. The configuration for the 443 simple physics test is given in Table S2. 444 refinement mode. The mesh is fixed at X4, with ) = /36 and ) = /4. In the control run (G6X4L2; Fig. 9a), 446 * = /36, = (1/2) ! , and * = /12. The 1 st -refinement region of this VR mesh has ~40-km resolution. Two 447 cases are used to examine the impact of the mesh distribution. In the first case, we choose to use more rapid resolution 448 changes in the inner transition zone based on the two parameters * and . * controls the width of the inner 449 transition zone, and represents the densification ratio between the 1 st -and 2 nd -refinement regions. Either narrowing 450 the width of the transition zone (Fig. 9b) or enlarging the inner densification ratio (Fig. 9c) generates a more abrupt 451 transition zone. The other way is to have the transition zone affect the tropical cyclone at an earlier stage. The initial 452 cyclone is placed closer to or even partly within the transition zone by increasing * (Fig. 9d), a parameter that controls 453 the radius of the 1 st -refinement region. The tropical cyclone is initialized at 10°N, 180°E over the 2 nd -refinement region, 454 near the transition zone between the 1 st -refinement and the 2 nd -refinement regions. 455 On day 10, all four tests well simulate the shape of the tropical cyclone. During its movement from the 2 nd -456 refinement into the 1 st -refinement region, the change in the grid size leads to little distortion on the tropical cyclone. 457 Compared to the preprint, two major differences can be found. First, the minor disturbances near the major tropical 458 cyclones almost disappear in this version, because activating hyperdiffuion damps the wind field more effectively. 459 Though it can be damped, this minor disturbance is not that unrealistic because it is near the major cyclone. https://www.earthsystemcog.org/projects/dcmip-2016/). The other difference is that the maximum wind magnitude in 463 the eyewall is stronger overall in this version. This is due to the locally scaled Smagorinsky formulation, as mentioned 464 in Section 2.3.1. 465 Figure 10 presents the evolution of the tropical cyclone in each test on days 2, 4, 6, and 8. The tropical cyclones in 466 all tests are consistent. No discernable difference is found when they move across the transition. A slightly larger 467 difference is more evident in the early stage (day 2), but such differences diminish as the tropical cyclones move into 468 the refinement center. The relative vorticity field also looks good and does not show any artifacts (Fig. S4). This 469 indicates that the mesh transition does not create notable problems. 470 To further examine possible sensitivity, we performed a group of experiments by altering one of the three parameters: 471 * , * , and . Figure 11 shows the minimum surface pressure (Fig. 11a-11c) and the maximum wind speed at 850 472 hPa ( Fig. 11e-11f). All these tests use a DTP splitting number 1 : 4 : 8, with a dycore step of 60 s (Table S2). Only one 473 test with the smallest * ( ! "# ) is unstable at this DTP splitting number, so we adjusted it to 1 : 2 : 4 with the same dycore 474 step. This suggests that an overly narrow inner transition zone can impose a higher stability restriction. The tropical 475 cyclone rapidly strengthens during the first two days. After day 2, it enters into a moderately developing stage in each 476 experiment. The evolution of intensity is diverse. All the VR runs tend to produce stronger tropical cyclones than G6 477 or G7, because the fine-resolution area determines the ultimate strength. The tests have run-to-run differences, but they 478 are small overall, showing consistent results and robustness. 479 Figure 12 shows a further comparison between two VR runs (G6X4L2 and G7X4L2) and two QU runs (G7 and 480 G8). The highest resolution of G6X4L2 and G7X4L2 is slightly higher than that of G7 and G8, respectively (see Table  481 S2). As shown, in each comparison (G6X4L2 vs G7, G7X4L2 vs G8), the VR model produces stronger wind 482 magnitudes in the eyewall and a more compact size featuring a smaller area coverage. The eyewall of the cyclone 483 converges towards its center, almost within 1 degree from the center. The maximum winds develop to higher vertical 484 levels as the local resolution increases. Overall, the difference between the VR and QU tests is attributed to different 485 local resolutions. 486 In the tropical cyclone test, if an unscaled Smagorinsky eddy viscosity is used (Fig. S5; preprint), the VR model 487 largely reduced such sensitivity, closer to the behavior of the QU model (Fig. S5). This confirms the reasonable behavior 489 of the local length formulation. 490

Full physics 491
The simple physics test, while insightful, is limited in the sense that the physics processes are simplified and do not 492 support enough of the nonlinear feedback typical of a real-world model. A VR model may introduce a higher degree of 493 nonlinear feedback due to mesh refinement. Thus, it is useful to check its behavior given a full parameterization suite. 494 We use a more highly variable mesh: G6B3X16L4. The finest part on this mesh reaches ~7-10 km. The refinement 495 center is placed at 35°N, 165°E. The sea surface temperature is identical to that in the simple physics test (29℃ 496 uniformly). The starting date is the first day of June. The solar constant is 1370 W.m −2 . The DTP splitting number is 1 : 497 5 : 10 with a dycore step of 10 s. If activated, the square of the Smagorinsky coefficient is 0.005, and the reference 498 hyperviscosity is 2 × 10 *; m 4 .s −1 with a reference length scale of 7000 m. 499 We first focus on the impact of the explicit diffusion process. Figure 13 shows the tropical cyclones on day 10 from 500 three tests: no explicit diffusion (noDiff), using hyperdiffusion only for the horizontal velocity (hyper), and 501 hyperdiffusion plus Smagorinsky (hyper+smg). The wind speed is shown at the model level nearest to 850 hPa (level 502 24). On day 10, the tropical cyclone center just moves into the refinement center (35°N, 165°E). In general, the results 503 in the three tests are consistent. The major difference lies in the eyewall. The noDiff and hyper tests produce slightly 504 higher wind maxima than hyper+smg. This suggests that Smagorinsky diffusion becomes active in the eyewall, 505 diffusing the solutions a little bit more strongly. Except for this minor difference (one may need to zoom in to see it), 506 hyper and hyper+smg produce very close solutions. 507 Explicit diffusion, while it seems to be irrelevant to the performance over the fine-resolution region, plays a more 508 important role in the coarse-resolution region. Due to the globally uniform sea surface temperature, near the South Pole, 509 some isolated systems can be generated in the noDiff test (Fig. 14a). Unlike the possible minor disturbances near the 510 major tropical cyclone in the simple physics test, these systems are not expected there because they are far away from 511 the initial vortex. They are likely to be caused by the nonlinear feedback between dynamics and physics. Only using 512 the Smagorinsky diffusion can suppress them to a certain extent, but not enough due to its flow-dependent nature (Fig.  513 14b). By further activating the hyperdiffusion option, these systems can be effectively removed (Fig. 14c). This suggests 514 that the explicit diffusion configuration is able to achieve a balance between underdiffusion and overdiffusion. 515 With this X16 mesh, it is also useful to examine whether a single refinement with a high densification ratio 516 (G6B3X16) tends to create serious problems during the simulation, especially when the tropical cyclone moves across 517 the transition. The refinement center is rotated to 60°N, 165°E. The tropical cyclone will encounter more mesh 518 transitions in G6B3X16. The results are shown in Fig. 15. On day 2, the tropical cyclone in G6B3X16L4 is experiencing 519 mesh cells between ~30 km and ~50 km. In G6B3X16, the encountered mesh sizes are above 70 km. Hence, 520 G6B3X16L4 has a stronger and more developed eyewall at this stage. On day 6, the eyewall enters into a region with 521 cell sizes ~25-30 km in G6B3X16L4, and ~50 km in G6B3X16. G6B3X16L4 still has a more compact eyewall, but 522 the difference between the two runs becomes smaller. On day 10, the eyewalls in both tests enter into an area with cell 523 sizes ~20-25 km. At this stage, both cyclone systems show similar maximum wind magnitudes in the eyewall, and 524 exhibit a similar distribution. Clearly, the fine-scale structure is determined mainly by the fine resolution, although two 525 systems can have larger differences before they enter into the fine-resolution area. These results also demonstrate that 526 a single refinement with a high densification ratio, though more challenging, can perform competitively with a more 527 gradually refined mesh when properly configured (e.g., an overly narrow transition is undesirable). The choice of 528 different mesh styles thus demands more numerical modeling experience and depends on the target issue. This needs 529 to be further examined with the aid of more realistic weather and climate system. 530 In this study, an atmospheric model formulated on an unstructured mesh has been systematically configured and 532 evaluated in its VR mode. Different mesh-refinement styles have been utilized to evaluate the model performance under 533 increasing degrees of complexity, from dry dynamics to simple and then to full physics. Based on some additional NDC 534 tests with full physics, it is found that the QU and VR models (with the same degrees of freedom) have comparable 535 parallel efficiency, that is, the VR configuration does not cause a degeneration in strong scaling performance. This is to 536 be expected as the domain decomposition does not distinguish between QU and VR grids. Other scaling performance 537 (e.g., whether the speedup ratio can achieve the ideal value by simply reducing the total grid number, with all other 538 things being equal) can be found in Liu et al. (2020). Overall, these results demonstrate that the VR configuration of 539 GRIST is a reliable and economic alternative to high-resolution quasi-uniform modeling. The adverse impact due to 540 the mesh transition and the coarse-resolution area can be well controlled. The major conclusions may be summarized 541 in three aspects. 542 On the overall performance. Based on the dry baroclinic wave test, all VR styles can well capture the fine-scale 543 wave structures. Such fine-scale resolving capability is supported by analysis of regional kinetic energy spectra, and is 544 further verified in the multi-region refinement mode. In the transition zone, the waves are not adversely affected by the 545 mesh refinement. In the coarse-resolution region, the VR model can also simulate an equivalent distribution of waves 546 to its low-resolution counterpart. A VR model indeed produces greater solution errors compared to its QU counterpart 547 that has the same degrees of freedom. However, the solution error can be reduced when the coarse region increases its 548 resolution. Thus, its impact can be controlled. In the tropical cyclone test, the VR model in simple and full physics can 549 simulate the gradual evolution of the tropical cyclone, showing reasonable resolution sensitivity. Importantly, the 550 simulation of the fine-scale structure is controlled mainly by the fine resolution, although larger differences may exist 551 before the systems move into the refinement area. Overall, the adverse impact due to mesh transition and the coarse-552 mesh region can be well controlled. 553 On the impact of explicit diffusion. Comparing this version with its earlier version (preprint), the impact of explicit 554 diffusion can be clearly demonstrated. It has been shown that activating scaled hyperdiffuion for the horizontal velocity 555 as background diffusion is helpful in alleviating grid-scale oscillations due to mesh transitions. It also suppresses some 556 highly unrealistic disturbances found in the full physics scenario. When activating hyperdiffuion, it is suggested that a 557 minimally required hyperviscosity should be used, such that it neither significantly damps the field nor becomes unable 558 to suppress undesirable grid-scale disturbances. Using only Smagorinsky diffusion would require higher coefficients to 559 remove oscillations, which, in turn, would restrict the numerical stability. Even so, the Smagorinsky-only configuration 560 is not enough for the VR mode, because it becomes inactive over certain regions that are dominated by weak flow 561 deformation. It is also shown that, in the tropical cyclone test, the scaled Smagorinsky diffusion has much lower 562 parametric sensitivity than its unscaled counterpart. In general, the use of explicit diffusion (properly scaled) plays a 563 positive role in the VR configuration, although there were concerns that their own discretization may introduce 564 additional problems. The increased demands of explicit diffusion from the QU to VR configuration is in accordance 565 with increasing mesh discontinuity due to transitions, and thus should not be viewed as a disadvantage. 566 On the impact of the mesh styles. Different mesh styles are all able to produce fine-scale structures. They maintain 567 solution quality over the mesh transition and the coarse-resolution area. In the dry test, it has been shown that G5B3X4, 568 while it has fewer cells than SURX4, leads to smaller solution errors and less oscillatory solutions. As these two meshes 569 are generated given a similar number of iterations, this might suggest that the random generators need more computation 570 to achieve a quality mesh. In general, when generating a VR mesh, we empirically suggest that one would be better to 571 start from a subdivided icosahedron. In the tropical cyclone case, a series of sensitivity tests examined the model 572 performance in a hierarchical refinement style. The solutions exhibit consistency even when the VR mesh is slightly 573 zone should be avoided. In the full physics test, G6B3X16L4 and G6B3X16 produce consistent results on day 10, 575 although the simulations at the initial stage (e.g., day 2) exhibit larger differences. 576 Code and data availability: GRIST is available at https://github.com/grist-dev. A version of the model code, and 577 running and postprocessing scripts for supporting this paper are available at: https://zenodo.org/record/3930643. 578 Version 2 is for this manuscript, and version 1 is for the preprint. The grid data used to enable the tests are located at 579 https://zenodo.org/record/3817060. Public code access is available after authorization, following 580 https://github.com/GRIST-Dev/TermsAndConditions. 581 Appendix: Regional kinetic energy spectra analysis 582 The model data are first interpolated to the Gaussian grid. For the QU model, we use the Gaussian grid that is 583 close to the nominal resolution of the icosahedral mesh (T106 for G6, T213 for G7, and T426 for G8). For the VR 584 model, the data are interpolated based on the fine resolution of the mesh (T221 for G6X4, T328 for G5B3X4, and T328 585 for SURX4). The selected regional domain (red box in Fig. 3) ranges from 150°E to 150°W and from 30°N to 65°N, 586 with two vortices. 587 We follow the approach documented by Denis et al. (2002) for computing kinetic energy spectra on a selected 588 regional domain. Let < ( , ) be the discrete cosine transform (DCT) of a 2D relative vorticity field ( , ) of , -589 by--grid points. The variance array can be computed from the DCT field as: 590 where L ) ( , ) is the variance array of divergence. The total kinetic energy as a function of wavenumber is given 603  shows the results on days 2, 6, and 10 on the G6B3X16L4 mesh, the right-hand column shows the corresponding 830 results on the G6B3X16 mesh. The results are shown for the wind speed (m.s −1 ) at the model level nearest to 850 hPa 831 (level 24). 832