DCMIP2016: The Splitting Supercell Test Case

This paper describes the splitting supercell idealized test case used in the 2016 Dynamical Core Model Intercomparison Project (DCMIP2016). These storms are useful testbeds for global atmospheric models because the horizontal scale of convective plumes is O(1km), emphasizing non-hydrostatic dynamics. The test case simulates a supercell on a reduced radius sphere with nominal resolutions ranging from 4km to 0.5km and is based on the work of Klemp et al. (2015). Models are initialized with an atmospheric environment conducive to supercell formation and forced with a small thermal perturbation. 5 A simplified Kessler microphysics scheme is coupled to the dynamical core to represent moist processes. Reference solutions for DCMIP2016 models are presented. Storm evolution is broadly similar between models, although differences in final solution exist. These differences are hypothesized to result from different numerical discretizations, physics-dynamics coupling, and numerical diffusion. Intramodel solutions generally converge as models approach 0.5km resolution, although exploratory simulations at 0.25km imply some dynamical cores require more refinement to fully converge. These results can be used as a 10

Supercells are strong, long-lived convective cells containing deep, persistent rotating updrafts that operate on spatial scales O(10km). They can persist for many hours and frequently produce large hail, tornados, damaging straight line winds, cloudto-ground lightning, and heavy rain (Browning, 1964;Lemon and Doswell, 1979;Doswell and Burgess, 1993). Therefore, accurate simulation of these features is of great societal interest and critical for atmospheric models. 5 The supercell test applied in DCMIP2016 (Ullrich et al., 2017) permits the study of a non-hydrostatic moist flow field with strong vertical velocities and the associated precipitation. This test is based on the work of Klemp and Wilhelmson (1978) and Klemp et al. (2015) and assesses the performance of global models at extremely high spatial resolution. It has recently been used in the development of next-generation numerical weather prediction systems (Ji and Toepfer, 2016).
While some recent work has targeted the impact of dynamical cores on extreme weather -primarily tropical cyclones 10 (e.g., Zhao et al. (2012); Reed et al. (2015)) -these studies have almost exclusively employed hydrostatic dynamical cores at grid spacings approximately 0.25 • and coarser. Conversely, the supercell test here emphasizes resolved, non-hydrostatic dynamics. In this regime the effective grid spacing is very similar to the horizontal scale of convective plumes. Further, the addition of simplified moist physics injects energy near the grid-scale in a conditionally-unstable atmosphere, which imposes significant stress on model numerics. The supercell test case therefore sheds light on the interplay of the dynamical core and 15 subgrid parameterizations and highlights the impact of both implicit and explicit numerical diffusion on model solutions. It also demonstrates credibility of a global modeling framework to simulate extreme phenomena, essential for future weather and climate simulations.

Description of test
The test case is defined as follows. All test cases are run on a non-rotating, reduced radius sphere with scaling factor X = 120. 20 Reducing the model's planetary radius allows for fine grid spacing to be achieved without the increased computational expense associated with adding grid cells to a standard global mesh in order to achieve non-hydrostatic resolutions (Kuang et al., 2005). Wedi and Smolarkiewicz (2009) provide a detailed overview of the reduced-radius framework for testing global models. For a 1 • mesh, the grid spacing of the reduced radius sphere is approximately 1 • /X ∼ 111km/X ∼ 111km/120 ∼ 1km near the equator. Klemp et al. (2015) demonstrated excellent agreement between simulations using this value of X and those completed 25 on a flat plane with equivalent resolution. The model top is placed at 20km with uniform vertical grid spacing (∆z) equal to 500m, resulting in 40 full vertical levels. No surface drag is imposed at the lower boundary (free slip condition). Water vapor (q v ), cloud water (q c ) and rain water (q r ) are handled by a simple Kessler microphysics routine (Kessler, 1969). In particular, the Kessler microphysics used here is outlined in detail in Appendix C of Klemp et al. (2015) and code for reproducing this configuration is available via the DCMIP2016 repository (http://dx.doi.org/10.5281/zenodo.1298671).

30
All simulations are integrated for 120min. Outputs of the full three-dimensional prognostic fields as well as all variables pertaining to the microphysical routines were stored for post-processing with a frequency of every 15min or finer. Four different horizontal resolutions were specified; 4 • , 2 • , 1 • , and 0.5 • . When considering the radius reduction, this results in approximate grid spacings of 4km, 2km, 1km, and 0.5km, respectively. Note that here we use '(nominal) resolution ' and 'grid spacing' interchangeably to refer to the length of a single grid cell or distance between gridpoints. All relevant constants mentioned here and in the following section are defined in Table 1.

Mean atmospheric background
The mean atmospheric state is designed such that it consists of large instability (convective available potential energy (CAPE) 5 of approximately 2200 m 2 s −1 ) and strong low-level wind shear, both of which are strong precursors of supercell formation (Weisman and Klemp, 1982).
The definition of this test case relies on hydrostatic and cyclostrophic wind balance, written in terms of Exner pressure π and virtual potential temperature θ v as (1)

10
Defining u = u eq cos ϕ to maintain solid body rotation, where u eq is the equatorial wind velocity, these equations can be combined to eliminate π, leading to The wind velocity is analytically defined throughout the domain. Meridional and vertical wind is initially set to zero. The zonal wind is obtained from Pressure and temperature at the equator are obtained by iterating on hydrostatic balance with initial state and iteration procedure This iteration procedure generally converge to machine epsilon after approximately 10 iterations. The equatorial moisture profile is then extended through the entire domain, 10 q(z, ϕ) = q eq (z).
Once the equatorial profile has been constructed, the virtual potential temperature through the remainder of the domain can be computed by iterating on (2), Again, approximately 10 iterations are needed for convergence to machine epsilon. Once virtual potential temperature has 15 been computed throughout the domain, Exner pressure throughout the domain can be obtained from (1), and so 20 Note that, for (13-14), Smolarkiewicz et al. (2017) also derive an analytic solution for the meridional variation of the initial background state for shallow atmospheres.

Potential temperature perturbation
To initiate convection, a thermal perturbation is introduced into the initial potential temperature field: An additional iterative step is then required to bring the potential temperature perturbation into hydrostatic balance. Without this additional iteration, large vertical velocities will be generated as the flow rapidly adjusts to hydrostatic balance since the test does not possess strong non-hydrostatic characteristics at initialization. The test case is designed such that the thermal perturbation will induce a convective updraft immediately after initialization.
As rain water is generated by the microphysics, reduced buoyancy and a subsequent downdraft at the equator in combination with favorable vertical pressure gradients near the peripheral flanks of the storm will cause it to split into two counterrotating 10 cells that propagate transversely away from the equator until the end of the test (Rotunno andKlemp, 1982, 1985;Rotunno, 1993;Klemp et al., 2015).

Physical and Numerical Diffusion
As noted in Klemp et al. (2015), dissipation is an important process near the grid-scale, particularly in simulations investigating convection in unstable environments such as this. To represent this process and facilitate solution convergence as resolution is 15 increased for a given model, a second-order diffusion operator with a constant viscosity (value) is applied to all momentum equations (ν = 500 m 2 s −1 ) and scalar equations (ν = 1500 m 2 s −1 ). In the vertical, this diffusion is applied to the perturbation from the background state only in order to prevent the initial perturbation from mixing out.
Models that contributed supercell test results at DCMIP2016 are listed in Table 2. They are formally described in Ullrich et al. (2017) and the references therein. Further, specific versions of the code used in DCMIP2016 and access instructions are 20 also listed in Ullrich et al. (2017). Note that not all DCMIP2016 participating groups submitted results for this particular test.
Due to the multitude of differing implicit and explicit diffusion in the participating models, some groups chose to apply variations in how either horizontal or vertical diffusion were treated in this test case. Deviations from the above specified diffusion are as follows. CSU applied uniform three-dimensional second order diffusion with coefficients of ν = 1500 m 2 s −1 for q v and θ v , ν = 1000 m 2 s −1 for q c and q r , and ν = 500 m 2 s −1 for divergence and relative vorticity. FV 3 applied divergence 25 and vorticity damping separately to the velocity fields along the floating Lagrangian surface. A Smagorinsky diffusion is also applied to the horizontal wind. ICON applied constant horizontal second-order diffusion to the horizontal and vertical velocity components (ν = 500 m 2 s −1 ) as well as the scalar variables θ v and q v,c,r (ν = 1500 m 2 s −1 ). No explicit diffusion was applied in the vertical. NICAM applied a dynamically-defined fourth-order diffusion to all variables in the horizontal with vertical dissipation being implicitly handled by the model's vertical discretization.

30
The following section describes the results of the supercell test case at DCMIP2016, both from a intermodel time evolution perspective and intramodel sensitivity to model resolution and ensuing convergence. Note that there is no analytic solution for the test case, but features specific to supercells should be observed and are subsequently discussed. It is not the intent of this manuscript to formally explore the precise mechanisms for model spread or define particular solutions as superior, but rather, 5 to publish an overview set of results from a diverse group of global, non-hydrostatic models to be used for future development endeavours. Future work employing this test case in a more narrow sense can isolate some of the model design choices that impact supercell simulations. Structural differences also begin to emerge at 60min. For example, FVM, GEM, ACME-A, and TEMPEST all exhibit three local maxima in vertical velocity; two large updrafts mirrored about the equator with one small maximum still located over equator centered near the initial perturbation. Similar behavior is noted in the q r fields. This is in contrast with other models 25 which lack a third updraft on the equatorial plane. Generally speaking, q r maxima are collocated with the locations of maximum updraft velocities, and thereby conversion from q v and q c to q r in the Kessler microphysics.

Time evolution of supercell at control resolution
While the aggregate response of a single updraft eventually splitting into poleward-propagating symmetric storms about the equator is well-matched between the configurations, notable differences exist, particularly towards the end of the runs. At and TEMPEST both produce longitudinally-transverse storms that stretch towards the equator in addition to the two main cells.
Each of the splitting supercells split a second time in ICON, forming, in conjunction with a local maximum at the equator, five maxima of vertical velocity (and correspondingly rainwater). NICAM produces two core supercells (as more clearly evident in the q r field at 120min), but has noticeable alternating weak updrafts and downdrafts in the north-south space between the two storm cores.
The relative smoothness of the storms as measured by the vertical velocity and rain water fields also varies between models, particularly at later times. ACME-A, FVM, GEM, OLAM, and MPAS produce updrafts that are relatively free of additional, small-scale local extrema in the vicinity of the core of the splitting supercell. Conversely, CSU, FV 3 , ICON, NICAM, and 5 TEMPEST all exhibit solutions with additional convective structures, with multiple updraft maxima versus two coherent cells.
This spread is somewhat minimized when looking at rain water, implying that the overall dynamical character of the cells as noted by precipitation generation is more similar, with all models other than ICON showing cohesive rain water maxima O(10 g/kg).

Resolution sensitivity of supercell
10 Fig. 4 shows the same cross-section variables as Fig. 3 except across the four specified test resolutions (nominally 4km, 2km, 1km, 0.5km, from left to right) at test termination of 120min. Therefore, the third panel from the left for each model (1km) should match the fourth panel from the left for each model in Fig. 3.
As resolution increases (left to right) models show increasing horizontal structure in both the vertical velocity and rain water fields. Updraft velocity generally increases with resolution, particularly going from 4km to 2km, implying that the 15 supercell is underresolved at 4km resolution. This is supported by previous mesoscale simulations investigating supercells in other frameworks (Potvin and Flora, 2015;Schwartz et al., 2017), although it should be emphasized that this response is also subject to each numerical scheme's effective resolution (Skamarock, 2004) and that the resolvability of real-world supercells can depend on the size of individual storms.
At the highest resolutions, there is a distinct group of models that exhibit more small-scale structure, particularly in vertical 20 velocity, at +120min at higher resolutions. CSU, GEM, and NICAM appear to have the largest vertical velocity variability at 0.5km, while ACME-A, FVM, MPAS, and TEMPEST appear to produce the smoothest solutions. This result is likely due to the differences in explicit diffusion treatment as noted before, as well as differences in the numerical schemes' implicit diffusion, particularly given the large impact of dissipation on kinetic energy near the grid scale (Skamarock, 2004;Jablonowski and Williamson, 2011). Additional focused sensitivity runs varying explicit diffusion operators and magnitude may be insightful 25 for developers to explore. It is also hypothesized that differences in the coupling between the dynamical core and subgrid parameterizations may lead to some of these behaviors (e.g., Staniforth et al. (2002); Malardel (2010); Thatcher and Jablonowski (2016); Gross et al. (2018)) although more constrained simulations isolating physics-dynamics coupling in particular modeling frameworks is a target for future work. As before, rain water cross-sections tend to be less spatially variable at 0.5km than vertical velocity, although CSU and NICAM both show some additional local maxima in the field associated with some of the 30 aforementioned w maxima.

Convergence of global supercell quantities with resolution
While Fig. 4 highlights the structural convergence with resolution more storm-wide measures of supercell intensity are also of interest. Fig. 5 shows the maximum resolved updraft velocity over the global domain as a function of time for each dynamical core and each resolution (finer model resolution is denoted by progressively darker lines). Maximum updraft velocity is chosen as a metric of interest due to its common use in both observational and modeling studies of supercells. All models show 5 increasing updraft velocity as a function of resolution, further confirming that, at 4km, the supercell is underresolved dynamically. For the majority of models and integration times, the gap between 4km and 2km grid spacing is the largest in magnitude, with subsequent increases in updraft velocity being smaller as models further decrease horizontal grid spacing. At 0.5km, the majority of models are relatively converged, with FV 3 , ICON, and MPAS showing curves nearly on top of one another at these resolutions. Other models show larger differences between 0.5km and 1km curves, implying that these configuration may not 10 yet be converged in this bulk sense. Further grid refinement or modifications to the dissipation schemes are necessary to achieve convergence; this is left to the invidual modeling groups to verify.
The maximum updraft velocity as a function of resolution for particular model configurations varies quite widely. NICAM produces the weakest supercell, with velocities around 30 m s −1 at 0.5km, while ACME-A, TEMPEST, GEM, and CSU all produce supercells that surpass 55 m s −1 at some point during the supercell evolution. Models that have weaker supercells 15 at 0.5km tend to also have weaker supercells at 4km (e.g., NICAM) while the same is true for stronger supercells (e.g., TEMPEST), likely due to configuration sensitivity. This agrees with the already discussed structural plots (Fig. 4) which demonstrated model solutions were generally converging with resolution on an intramodel basis but not necessarily across models. Fig. 6 shows the same analysis except for area-integrated precipitation rate for each model and each resolution. Similar 20 results are noted as above -with most models showing large spread at the coarsest resolutions, but general convergence in precipitation by 0.5km. All models produce the most precipitation at 120min with the 4km simulation. This is consistent with Klemp et al. (2015), who postulated this behavior is due to increased spatial extent of available q r to fall out of the column at these grid spacings, even though updraft velocities are weaker at coarser resolutions. Unlike maximum vertical velocity, the integrated precipitation rate does not monotonically increase with resolution for most models. At 120min, integrated rates 25 at 0.5km range by approximately a factor of two, from a low of 70-75 x10 5 kg s −1 (ACME-A, FVM, OLAM) to a high of 150-160 x10 5 kg s −1 (CSU, FV 3 ), highlighting the sensitivity of final results that have been already been discussed.

Conclusions
Non-hydrostatic dynamics are required for accurate representation of supercells. The results from this test case show that clear differences and uncertainties exist in storm evolution when comparing identically initialized dynamical cores at similar 30 nominal grid resolutions. Intramodel convergence in bulk, integrated quantities appears to generally occur at approximately 0.5km grid spacing. However, intermodel differences are quite large even at these resolutions. For example, maximum updraft velocity within a storm between two models may vary by almost a factor of two even at the highest resolutions assessed at DCMIP2016.
Structural convergence is weaker than bulk integrated metrics. Two-dimensional horizontal cross-sections through the supercells at various times show that some models are well-converged between 1km and 0.5km, while results from other models imply that finer resolutions are needed to assess whether convergence will occur with a particular test case formulation and 5 model configuration. Interestingly, in some cases maximum, bulk quantities converge faster than snapshots of cross-sections.
We postulate that these differences and uncertainties likely stem from not only the numerical discretization and grid differences outlined in Ullrich et al. (2017), but also from the form and implementation of filtering mechanisms (either implicit or explicit) specific to each modeling center. The simulation of supercells at these resolutions are particularly sensitive to numerical diffusion since damping of prognostic variables in global models is occurring at or near the scales required for resolvability 10 of the storm. This is different from other DCMIP2016 tests (baroclinic wave and tropical cyclone), which produced dynamics that were less non-hydrostatic in nature and required resolvable scales well coarser than the grid cell level. Further, since DCMIP2016 did not formally specify a particular physics-dynamics coupling strategy, it would not be surprising for particular design choices regarding how the dynamical core is coupled to subgrid parameterizations to also impact results.
Given the lack of an analytic solution, we emphasize that the goal of this paper is not to define particular supercells as 15 optimal answers. Rather, the main intention of this test at DCMIP2016 was to produce a verifiable database for models to use as an initial comparison point when evaluating non-hydrostatic numerics in dynamical cores. Pushing grid spacings to 0.25km and beyond to formalize convergence would be a useful endeavour in future application of this test, either at the modeling center level or as part of future iterations of DCMIP. Variable-resolution or regionally-refined dynamical cores may reduce the burden of such simulations, making them more palatable for researchers with limited computing resources. 20 We acknowledge that as groups continue to develop non-hydrostatic modeling techniques that small changes in the treatment of diffusion in the dynamical core will likely lead to changes in their results from DCMIP2016. We recommend modeling centers developing or optimizing non-hydrostatic dynamical cores perform this test and compare their solutions to the baselines contained in this manuscript as a check of sanity relative to a large and diverse group of next-generation dynamical cores actively being developed within the atmospheric modeling community.

Code availability
Information on the availability of source code for the models featured in this paper can be found in Ullrich et al. (2017)