## Abstract

The sensitivity of an idealized squall line to horizontal and vertical grid spacing is investigated using a new approach. Simulations are first performed at a horizontal grid spacing of 1 km until the storm reaches its mature stage. The model output is then interpolated to smaller (and larger) grid spacings, and the model is restarted using the interpolated state plus small thermodynamic perturbations to spin up small-scale motions. This framework allows an investigation of the sensitivity of the storm to changes in without complications from differences in storm initiation and early evolution. The restarted simulations reach a quasi steady state within approximately 1 h. Results demonstrate that there are two -dependent regimes with the transition between regimes occurring for between 250 and 500 m. Some storm characteristics, such as the mean convective core area, change substantially for 250 m but show limited sensitivity as is decreased below 250 m, despite better resolving smaller-scale turbulent motions. This transition is found to be independent of the chosen . Mixing in the context of varying and is also investigated via passive tracers that are initialized 1 h after restarting the simulations (i.e., after the spin up of small-scale motions). The tracer field at the end of the simulations reveals that entrainment and detrainment are suppressed in the simulations with 500 m. For decreasing , entrainment and detrainment are substantially more important, limiting the flux of low-level tracer to the upper troposphere, which has important implications for modeling studies of convective transport from the boundary layer through the troposphere.

## 1. Introduction

Weather in the tropics and midlatitudes is strongly affected by organized deep convective systems, resulting in a substantial fraction of precipitation in these regions (e.g., Houze 1993; Markowski and Richardson 2010). Such systems are typically studied using so-called cloud-resolving or cloud system–resolving models (CRMs and CSRMs, respectively). These models have horizontal grid spacings of less than approximately 5–8 km, which is needed to begin to resolve nonhydrostatic motion (Weisman et al. 1997). Recently, these models have been used with increasing frequency for regional-scale numerical weather prediction (NWP), especially in the context of summertime continental convection (e.g., Done et al. 2004; Kain et al. 2008; Weisman et al. 2008; Schwartz et al. 2009; VandenBerg et al. 2014). Moreover, high-resolution CRMs are being used within coarse-resolution models (e.g., general circulation models, which are commonly referred to as the multiscale modeling framework or “superparameterization”) (e.g., Grabowski 2001; Randall et al. 2003; Tao et al. 2009).

Several studies have shown consistent differences in simulated deep convective systems with changes in horizontal grid spacing from *O*(100) m to *O*(10) km (e.g., Weisman et al. 1997; Grabowski et al. 1998; Belair and Mailhot 2001; Petch and Gray 2001; Adlerman and Droegemeier 2002; Bryan et al. 2003; Gentry and Lackmann 2010; Fiori et al. 2010; Bryan and Morrison 2012; Verelle et al. 2014). For example, horizontal grid spacings of 4 and 1 km were shown by Weisman et al. (1997) to capture similar mesoscale structures and nonhydrostatic motions, while larger grid spacings resulted in systems that evolved more slowly. However, several studies have shown that a horizontal grid spacing of *O*(1) km is insufficient to resolve turbulent eddies and cloud-scale motions (e.g., Grabowski 1998; Petch et al. 2002; Bryan et al. 2003; Rotunno et al. 2009; Dawson et al. 2010; Bryan and Morrison 2012; Morrison et al. 2015). Bryan et al. (2003) showed that grid spacings of 250 m or less were required to resolve an inertial subrange and for subgrid-scale closures designed for large-eddy simulation (LES) to perform according to their design. However, the statistical convergence of deep convection simulations with decreased grid spacing has not been well characterized. Craig and Dornbrack (2008) defined convergence to occur when differences in statistical quantities with different grid spacings were comparable to differences induced by small perturbations in the initial conditions. Based on this definition, they showed that convergence occurs for 0.2, where was defined as a length scale related to the ratio of thermal buoyancy and stratification (typically a few 100 m). Several other studies have found convergence in various bulk quantities, for example, horizontally averaged vertical fluxes, as was reduced to 200–500 m, using various definitions of statistical convergence (e.g., Fiori et al. 2010, 2011; Langhans et al. 2012; Verelle et al. 2014; Morrison et al. 2015). However, Bryan et al. (2003) found a lack of statistical convergence of the turbulent flow associated with deep convection with decreased to 125 m and argued that convergence may not be expected until the grid spacing is reduced to approximately 100 times smaller than the large-eddy scale (peak kinetic energy scale associated with vertical motion, approximately a few kilometers in deep convection).

Changes in grid spacing can also have large effects on transport and mixing in deep convection simulations. Entrainment is often represented in convective parameterizations by assuming that the fractional rates exhibit an inverse dependence on the convective core radius *r* (e.g., Simpson 1971).^{1} Because the average radius of convective towers decreases substantially as the horizontal grid spacing of the model decreases to a few hundred meters (Bryan and Morrison 2012), this suggests an increase in entrainment, all else being equal. Bryan and Morrison (2012) found a pronounced difference in the transport of passive tracers with a decrease in horizontal grid spacing from 1 km to 250 m. This difference was largely attributable to differences in entrainment and detrainment of convective updrafts, at least partly because of a reduction in updraft core size. However, the extent to which this increase in entrainment is important for the simulation of continental convection, especially squall lines, remains uncertain and is a focal point of this work.

Based on the aforementioned studies, simulations with grid spacings of *O*(10) m may be necessary to investigate the convergence of simulated convective characteristics and mixing (e.g., Craig and Dornbrack 2008; Bryan et al. 2003), that is, two orders of magnitude smaller in the horizontal plane compared to typical CRM grid spacings (i.e., 1 km). This grid spacing increases the computational expense by at least six orders of magnitude (assuming a decrease in time step commensurate with the decrease in horizontal grid spacing), excluding any commensurate decrease in . Even with recent advances in computational power, such simulations are on the edge of what is now possible, especially for simulations exceeding 2 h, which is typically needed to simulate the life cycle of organized convective systems. Moreover, even if the computational capabilities were in place for such simulations, the resulting sensitivity may be a manifestation of differences in the initialization and early storm evolution and not necessarily due to differences in representing the mature stage of convective systems. In particular, slower storm initiation and evolution as is decreased (Weisman et al. 1997) can lead to challenges in interpreting this sensitivity. To address these issues, we propose a new approach for simulating organized deep convection with very small grid spacings. Briefly, this approach requires a simulation run at a relatively modest for a sufficient amount of time to allow the modeled system to spin up, which is determined based on the time required for the system to reach a quasi steady state (e.g., nearly constant precipitation rate and/or mean updraft velocity). Model output fields are then interpolated to smaller (or larger) grid spacings, small random perturbations to the potential temperature field are added to facilitate the generation of small-scale motion, and the model is restarted. Such an approach is found to be far more efficient than the traditional approach of using a very small for an entire simulation and alleviates the aforementioned issue regarding differences in model initialization and early system evolution. This approach allows for a direct comparison of mature, quasi-steady convective systems simulated with various grid spacings.

In this study, we expand upon Bryan et al. (2003) and Bryan and Morrison (2012) by exploring the sensitivity of deep convection simulations to a wider range of horizontal grid spacings using the new approach and by examining the sensitivity to changes in vertical grid spacing . We also provide a detailed tracer analysis to examine how entrainment, detrainment, and mixing vary as a function of grid spacing (both horizontal and vertical). The overall aim is to understand how the ability to resolve the convective core size and smaller-scale turbulent motions affects entrainment and detrainment pertinent to convective transport in squall lines.

The remainder of this paper is organized as follows: The new technique for very small simulations of deep convective systems and the modeling framework are presented in section 2. The results of the simulations are presented in section 3. A discussion and the main conclusions are provided in sections 4 and 5, respectively.

## 2. Methods

### a. Experimental setup

We employ the Weather Research and Forecasting (WRF) Model, version 3.3.1 (Skamarock et al. 2008), to test the sensitivity of simulated squall lines to and . WRF is a compressible, nonhydrostatic model that can be run at drastically different grid spacings, thus allowing both CRM and LES to be performed. The model domain is set to 300 km by 60 km in the line-normal and line-parallel directions, respectively, and 24 km in the vertical. A 5-km Rayleigh damping layer is applied at the top of the model domain. The boundaries are assumed to be open in the line-normal direction and periodic in the line-parallel direction. The 1.5-order turbulent kinetic energy (TKE) closure scheme is used (e.g., Skamarock et al. 2008). Stein et al. (2015) showed sensitivity to the assumed subgrid-scale turbulence scheme. However, for consistency among the simulations discussed herein, we examine the results using a single scheme. For simplicity in these idealized simulations, we neglect the effects of radiation, Coriolis acceleration, and surface fluxes. All scalars are advected following standard fifth- and third-order monotonic advection schemes in the horizontal and vertical, respectively (Wang et al. 2009). Surface drag is neglected in the simulations (i.e., a free-slip condition is applied at the lowest model level).

To initiate convection, the vertical velocity *w* is forced over the first hour of the simulations (Ziegler et al. 2010). The forcing is applied within a half cylinder with a radius of 10 km in the line-normal direction and 2.5 km in the vertical (the flat side of the half cylinder is located at *z * 0); the maximum acceleration of 0.5 m s^{−2} is located at the center of the half cylinder. (Note that the *w* forcing is uniform in the line-parallel dimension.) The *w* forcing decays radially from the center assuming a cosine function of the radius. Three-dimensional (3D) motion is initiated by applying random potential temperature perturbations (amplitude 0.1 K) within a 40-km-wide region positioned at the center of the domain in the line-parallel direction and within the lowest 4 km.

The model is initialized with the thermodynamic sounding shown in Fig. 1. This sounding is similar to that used in the continental (i.e., Oklahoma) squall-line case study of the Eighth International Cloud Modeling Workshop held in July 2012 in Warsaw, Poland (Muhlbauer et al. 2013). Additional details of the sounding can be found in Morrison et al. (2012) and Lebo and Morrison (2014). One important difference between the sounding implemented here and that used in prior studies is that we limit the convective available potential energy (CAPE) to 4200 J kg^{−1} based on the most energetic parcel (as opposed to approximately 6800 J kg^{−1} for the original sounding) by reducing the initial water vapor mixing ratio below a height of 2.5 km. We prescribe a base state low-level wind shear of 0.0024 s^{−1} over the lowest 5 km; the line-normal wind component (i.e., toward the line) decreases linearly with height from the surface to the top of the shear layer. Note that the mean wind is subtracted from the profile and thus the line stays nearly stationary throughout the duration of the simulations (equivalent to a constant horizontal translation of the domain). As a result, a smaller domain relative to that used in Lebo and Morrison (2014) is permitted (i.e., 300 km × 60 km vs 714 km × 124 km), and features of interest stay well away from the line-normal lateral domain boundaries.

We employ a two-moment bulk microphysics scheme (Morrison et al. 2009) for the grid spacing sensitivity simulations. Morrison et al. (2012) presented a detailed comparison of WRF simulations using the Morrison et al. (2009) microphysics scheme with observations for this case using a similar model setup. In particular, the model was able to reproduce the overall observed reflectivity field with a well-defined convective core and trailing stratiform region.

### b. Sensitivity to grid spacing

Allowing the simulations at various to spin up independently presents an interesting problem in and of itself. In doing so, the systems may evolve in different ways. To ensure that we focus solely on the effects of on a mature squall line, a spinup simulation is performed using a modest (i.e., 1 km; a grid spacing of 500 m was also tested, although the suite of simulations resulted in qualitatively similar findings). This spinup simulation is performed for 4 h as suggested by a previous study using a similar domain setup and sounding (Lebo and Morrison 2014). After the initial 4 h, the model fields, including all microphysical variables, are interpolated (bilinear interpolation) to smaller horizontal grid spacings, namely, 500, 250, 100, 66.67, and 33.33 m; an additional test is performed by increasing the grid spacing to 2 km. To facilitate the generation of smaller-scale 3D motions, a random ±0.05-K potential temperature perturbation is applied to the interpolated fields throughout the domain. A WRF restart file is generated from the interpolated model domain and used to restart the model with the new . This process is performed for two (i.e., 250 and 125 m), although the 33.33 m simulation is not included in the 125 m simulations because of computational constraints. The vertical grid spacings are chosen for two reasons: 1) Using a constant for each set of simulations allows us to explore the effects of changes in the vertical and horizontal grid spacings independently. 2) Using becomes computationally unfeasible for small grid spacings. The chosen values are nearly halfway between the maximum and minimum . Using this new approach, we are assured that all of the simulated convective systems are starting at nearly the same mature stage of the storm; therefore, we can more easily determine the effects of changes in . After the restart, the model is run for 2 h. We restrict most of the analysis to the final hour of the simulations, after smaller-scale motions are spun up in the restarted simulations. However, we do analyze the first 1 h of spinup after the restart to investigate the dynamical effects of differences in convective core sizes before small-scale motions become important.

To investigate the effects of changes in on mixing, passive tracers are added to the domain after the 1-h spinup period within two different regions: 1) below 3 km and 2) between 3 and 10.5 km. These regions are chosen to address the issue of mixing of low-level, midtropospheric, and upper-tropospheric air in squall lines, motivated, for example, by the potential effects of boundary layer/low-level aerosols on mid- to upper-troposphere convective processes. The tracers are only added for the final hour of the simulations to ensure that the transport occurs after small-scale motion is spun up. Additional tests suggest that longer simulations provide little additional information about mixing in the squall lines because the systems are in a quasi steady state during the final hour. The tracers are treated like all other scalar (microphysical) variables in the model, that is, being affected by advection, numerical diffusion, and subgrid-scale mixing, except that they do not affect the buoyancy.

## 3. Results

Previous studies have shown that the width of convective cores decreases as decreases (i.e., increased resolution; see Bryan and Morrison 2012; Morrison et al. 2015). Moreover, as decreases, the model is more capable of capturing smaller-scale turbulent motions that act to detrain (entrain) cloudy (dry) air from (into) convective cores, effectively mixing the perimeter of the core more efficiently. Both of these aspects are expected to affect mixing, entrainment, and detrainment in convective updrafts and are investigated below.

### a. General depiction of the simulated system

Before delving into results from the detailed analysis, a depiction of the structure of the simulated squall line at the various horizontal grid spacings is useful for both suites of simulations (i.e., both ). Figure 2 shows excerpts (20 km × 20 km) of the plan view of the squall-line state at the end of the restarted simulations. Both the total mass mixing ratio (i.e., sum of the cloud, rain, ice, snow, and graupel mixing ratios) and the vertical velocity are shown at a height of 5 km for the suite of simulations performed with 250 m. The updraft area per convective core clearly decreases for smaller . The updrafts are less frequent and substantially larger at the largest . Interestingly, we find little difference in the plan views for ranging from 33.33 to 100 m. The plan views for of 125 m are qualitatively similar to those shown in Fig. 2 (not shown).

Figures 3 and 4 show transects of the vertical velocity (colored) and line-normal wind (contours) for of 250 and 125 m, respectively. The presented transects are at the location of the maximum vertical velocity in the domain, providing a representation of the vertical structure of the squall line at its strongest point. Warm (cool) colors correspond to line-parallel (i.e., average in the *y* direction) upward (downward)-moving air. Analogous to the results presented in Fig. 2, regardless of , there is a trend toward narrower updrafts for decreasing throughout the vertical, although this trend is most dramatic for a change in from 500 to 250 m.

Figures 2, 3, and 4 provide our first evidence of two important results that will be further analyzed and discussed in the remainder of this section: 1) The simulated structures for a given appear to be mostly independent of the change in , that is, a change in by a factor of 2 has little effect on the results presented herein. However, the vertical grid spacing may be more important in less idealized simulations that do not neglect surface fluxes, which can also affect the convective characteristics. 2) There are two distinct regimes, one for large (i.e., 500 m and greater) and for small (i.e., 250 m and smaller), with a sharp transition between these regimes.

### b. Time series of convective parameters

To illustrate the performance of the proposed downscaling method for very small , large-eddy simulations of organized deep convection (specifically squall lines), time series of several convective parameters are useful. Figure 5 ( 250 m) and Fig. 6 ( 125 m) depict the domain maximum vertical velocity and the mean vertical velocity conditionally averaged over only locations in which *w * 2 m s^{−1} (to constrain the results to only regions of rapid ascent). Similar trends appear for both as a function of changes in . For example, immediately following the restart, both the maximum and average updraft velocities increase rapidly, reaching a maximum after approximately 5 min (except in the 2 km simulation, for which a minimum occurs at this time). This rapid increase in updraft strength likely reflects the change in the perturbation pressure forcing that is consistent with a decrease in convective core size caused by the decreased (e.g., Morrison 2015a,b). Moreover, the maximum vertical velocities that are attained approximately 5 min into the restarted simulations systematically increase as the horizontal grid spacing decreases, reaching 65 (for 33.33 m) and 69.5 m s^{−1} (for 66.67 m) for of 250 and 125 m, respectively. For comparison, the thermodynamic maximum *w*, that is, , for these simulations is 92.3 m s^{−1} (here, CAPE 4259 J kg^{−1}). Therefore, even with a grid spacing of 100 m (or less), the simulated vertical velocity only attains approximately 75% of the thermodynamic maximum, which is likely related to condensate loading in the updrafts (which decreases buoyancy).

After the initial peak in the maximum and average vertical velocities in the 1 km simulations, entrainment acts to reduce the vertical velocity in the convective cores; therefore, both the maximum and mean vertical velocities decrease rapidly after 5 min into the restarted simulations. Thereafter, the maximum vertical velocities are very noisy (Figs. 5a, 6a); after approximately 1 h, systematic changes in the maximum vertical velocity are no longer apparent. However, this is not the case for the mean updraft velocity. Figures 5b and 6b demonstrate that there is a rapid reduction in the mean vertical velocity after the peak until approximately 20 min into the restarted simulations for 1 km and an increase during this time period in the 2-km simulation, except for the simulation with 33.33 m, which requires additional time to spin up after the model is restarted. Overall, this finding suggests that convection rapidly adjusts to the new . Furthermore, while the initial effect for a decrease in is to increase the mean vertical velocity, after the system relaxes to its new quasi steady state, the mean updraft velocity decreases for a decrease in . With that said, Fig. 5b suggests convergence in the mean updraft velocity because it is similar as is decreased from 100 to 33.33 m. Interestingly, the mean vertical velocity approaches approximately 4–4.5 m s^{−1} as is increased for both (Figs. 5b, 6b). Similarly, regardless of , the maximum vertical velocities cluster primarily between 35 and 45 m s^{−1} over the final 30 min of the simulations (Figs. 5a, 6a). Thus, changes in do not appear to play a critical role in the sensitivity of vertical velocity to over the range of tested here.

### c. Energy spectra

An important component of atmospheric modeling is the ability to accurately represent 3D air motions. Typically, the ability of a model to represent these motions is assessed in terms of the kinetic energy power spectra. A fast Fourier transformation (FFT) is utilized to convert the vertical velocity field *w* from spatial coordinates to frequencies (or wavelengths). Previous studies have computed the energy spectra in a single direction within the domain (e.g., Bryan et al. 2003). However, the variability in a single dimension may not necessarily characterize the 3D characteristics of the vertical velocity field for sheared systems (a squall line), especially at the mesoscale. Therefore, we use a 1D FFT to determine the vertical velocity energy spectra in both the *x* [] and *y* [] directions, which are functions of the wavenumber *κ*. Then, we define the 2D energy spectra, that is, , as

The spectra are vertically averaged between heights of 4 and 5 km. The results of these calculations are shown in Fig. 7a after 1 h into the restarted simulations and in Fig. 7b at the end of the 2-h restarted simulations; the spectra are truncated at 6, which is the approximate effective model resolution. The choice of truncating the spectra at follows the analysis of Bryan et al. (2003) and is close to the effective grid spacing of 7 reported for a similar WRF configuration by Skamarock (2004). The FFT implicitly assumes that the boundaries are periodic, which is true in the *y* direction; however, this is not the case in the *x* direction. For the system simulated herein, this is found to have little effect on the power spectra generated using the FFT (i.e., the results in the *x* direction are qualitatively consistent with those obtained in the *y* direction; not shown). The spectra are shown for two different times in Fig. 7 to demonstrate that the model has reached a quasi-steady kinetic energy state after only 1 h, which is seen by the similarity between the spectra after 1 and 2 h.

The energy spectra depicted in Fig. 7 demonstrate several important concepts. Regardless of , the energy spectra do not approach a slope of − for 500 m. For 250 m, a − slope in the vertical velocity energy spectra is obtained, at least for a portion of the smallest resolved wavelengths. The fact that the peak in the energy spectra for the large simulations occurs at greater wavelengths than in the simulations with of 250 m or less suggests that the simulations performed at these cannot properly resolve the most important scales of vertical motion, that is, wavelengths less than approximately 1.5 km. Given the expectation that the spectral peak is associated with vertical motions having a horizontal scale similar to that of the updraft cores, these results suggest an inability to properly resolve the updraft core size. This is confirmed by an analysis of updraft size discussed later.

Interestingly, Bryan et al. (2003) noted the existence of an inertial subrange for horizontal grid spacings of 250 and 125 m. However, they found that even at of 125 m, the slope of the power spectra was slightly smaller than the expected −. They speculated that this reduction in the slope may have been due to a grid spacing that was still not small enough to produce an adequately high Reynolds number. Bryan et al. (2003) suggested that a horizontal grid spacing of approximately 30 m may be necessary to achieve a − relationship. Figure 7a suggests that a nearly − slope is attained for a rather broad range of wavenumbers at a horizontal grid spacing of 33 m, while the same slope occurs for 66 and 100 m, except over a shorter range of wavenumbers. However, as shown in section 3d, the convective characteristics exhibit only small changes as a result of resolving additional details of turbulent motion with decreasing at very small . These results also suggest that grid anisotropy (or the grid cell aspect ratio) is not a critical factor controlling vertical motion, at least for scales larger than 6 and for the range of and tested here, because the spectra have a nearly − slope down to the truncation scale of 6 for of 100 m and less despite the increasing grid anisotropy. There is also little sensitivity in the spectra to a decrease in from 250 to 125 m.

### d. Squall-line statistics

To obtain a comprehensive perspective of the sensitivity of the simulated squall-line structure to and , several fundamental characteristics of deep convection are averaged over the final hour of the restarted simulations. The first hour is neglected in these calculations because it was shown in the previous section that this period of time is largely characterized by the system adjusting to a new quasi steady state. Figure 8 ( 250 m) and Fig. 9 ( 125 m) show the following characteristics (the vertical lines represent one standard deviation from the temporal mean based on 1-min model output):

Mean number of cores: The average number of convective cores. A convective core is defined using the vertically averaged vertical velocity between 5 and 10 km in height. This variable was chosen to define convective cores to ensure that only the deepest convective cores are considered, that is, cumulus congestus along the gust front is not included. Then, a convective core is defined by a region in which the local maximum exceeds a predefined threshold , covering an area in which exceeds a different threshold that defines the periphery of the cores . Various combinations of thresholds (i.e., 8 or 10 m s

^{−1}and 0 or 3 m s^{−1}) were tested to ensure the qualitative robustness of the results. The results presented herein are for 8 m s^{−1}and 3 m s^{−1}(Figs. 8a, 9a).Average core area: The total convective area divided by the number of convective cores (Figs. 8b, 9b).

Total core area: The total horizontal area that was defined according to the convective core thresholds (Figs. 8c, 9c).

Mean convective updraft mass flux: The product of the vertical velocity and air density for points with vertical velocities exceeding a predefined threshold (here, 2 m s

^{−1}is used, although little sensitivity to this threshold was found). The sum of these points is computed and divided by the total number of grid points in the domain (Figs. 8d, 9d).

The results suggest that the system exhibits two -dependent regimes. For of 500 m to 2 km, some of the convective parameters tend to exhibit a monotonic change with , for example, the total core area, while other characteristics, for example, the mean convective updraft mass flux, do not. Moreover, the average number of cores increases and the average core area decreases as is decreased in this regime, leading to small decreases in the total core area (Fig. 8c). However, the system enters a qualitatively different regime at approximately 250 m. At this scale, the convective parameters become fairly constant (or do not show a consistent trend) with further decreases in . This is especially true for 100 m. For all , there is a close correspondence between the wavelength of the peak in the *w* kinetic energy spectra (Fig. 7) and the mean updraft core radius, derived from the core area by dividing by *π* and taking the square root, which implies that the most energetic eddies are those on the scale of the updraft cores. Thus, the wavelength of the peak in the *w* energy spectrum provides a useful and precise definition of the mean convective core width scale, which does not depend on vertical velocity or buoyancy thresholds traditionally used to define updraft cores. While it was shown in section 3b that for decreasing the inertial subrange with a − spectral slope encompassed a larger range of resolved wavenumbers, this improvement in the representation of finer-scale turbulent motions seems to have little effect on the mean convective characteristics of the storm, suggesting that the subgrid-scale mixing scheme is functioning appropriately at of less than or equal to 100–250 m. These results are found to be insensitive to changes in (cf. Figs. 8 and 9).

### e. Tracer analysis

As discussed in section 3b, tracers are initialized at two levels 1 h after restarting the simulations. The 1-h delay allows the simulations to adjust to the new grid spacings and ensures that the resulting tracer fields are analyzed for quasi-steady conditions after small-scale motion is spun up. Here, we focus on the line-parallel state (i.e., average in the *y* direction) of the tracer fields (Fig. 10), joint PDFs of the tracer quantities and *w* (Figs. 11, 12), and PDFs of for *z * 10 km (Fig. 13). For reference, is the mixing ratio of the tracer initialized below a height of 3 km, and is the mixing ratio of the tracer initialized at heights between 3 and 10.5 km. In Fig. 10, regions in which the line-averaged *w* exceeds 0.5 m s^{−1} are stippled for reference. The stippled region clearly portrays both the region of convective updrafts at the leading edge of the squall line and the mesoscale updraft of the trailing stratiform region. At first glance, the panels shown in Fig. 10 appear to be quite similar. However, there are important differences in the simulations that correspond to the -dependent regimes discussed in the preceding subsection. For example, for 500 m, there is a pronounced, fairly narrow region near the leading edge of convection in which the ratio of to exceeds 0.5 in the mid- to upper levels (very light blue regions above approximately 9 km). This region is absent in the simulations performed with smaller , suggesting that the low-level tracer is more readily transported upward in the large simulations due to a lack of entrainment between the convective cores and the surrounding environment at midlevels and is also less readily detrained above 9 km (cf. Figs. 10a,b with, e.g., Fig. 10f).

To examine further the differences in entrainment and detrainment as a function of , joint PDFs are presented for and *w* (Fig. 11) and for and *w* (Fig. 12) for heights between 4 and 9 km. These levels are chosen because they contain no low-level tracer when the tracers are initially added to the domain. Therefore, if is nonzero, low-level tracer was transported upward by the system. The colors indicate frequency and are logarithmically spaced. (Note that the lack of blue points for the large simulations is a direct result of simply having fewer points.) Figure 11 corroborates the results presented above regarding the changes in entrainment/detrainment, especially for the change in from 500 to 250 m. Specifically, the joint PDFs can be separated into two groups, that is, Figs. 11a–c and Figs. 11d–g. For large , there are higher frequencies for at high *w* (i.e., exceeding 20 m s^{−1}), especially for of 500 m and 1 km, indicating that more low-level tracer remains in the convective cores at these larger . Interestingly, for the largest , there is an absence of points for which *w* exceeds 25 m s^{−1}, and these points increase in frequency as the grid spacing is reduced. The lack of high *w* for the largest is likely because of perturbation pressure effects (e.g., Morrison 2015a,b). Because the updraft size (specifically, the ratio of width to height) is large at a grid spacing of 2 km, it follows that the downward perturbation pressure force is large relative to the smaller simulations, which hinders the formation of strong updrafts. For the smaller simulations, that is, Figs. 11d–g, the frequency of high values (greater than 0.7) decreases relative to that in the larger simulations for *w* greater than 10 m s^{−1}, suggesting that convective cores are more diluted for smaller than approximately 250 m. Again, these results are found to exhibit only small quantitative changes and little qualitative differences for changes in (not shown).

Figure 12 provides joint PDFs for and *w* between 4 and 9 km. Here, we see a systematic decrease in for increasing *w* for all . Again, however, there is a change in the tails of the PDFs when is changed from 500 to 250 m, that is, the high frequency of low and high *w* decreases substantially (there is a lack of green boxes for *w* exceeding 25 m s^{−1}). This supports the results from the analysis of and indicates increased entrainment and dilution of convective cores as is decreased, especially from 500 to 250 m.

To reinforce these findings, Fig. 13 presents PDFs of at *z * 10 km for both . This figure confirms that decreasing acts to substantially decrease the amount of low-level tracer that is transported to the upper troposphere. The reduction of low-level tracer at upper levels, especially between of 500 and 250 m, is much greater than the corresponding reduction in the mean convective updraft mass flux (Figs. 8d, 9d). For example, the decrease in the frequency of large is more than an order of magnitude (Fig. 13), compared to a less than about 30% decrease in the mean convective updraft mass flux (Figs. 8d, 9d). This result implies that large increases in entrainment and detrainment at low and midlevels are the primary cause of the reduced occurrence of large above 10 km as is decreased. Moreover, the decrease in convective mass flux as is decreased is itself likely caused by the increased entrainment and mixing of updraft cores.

## 4. Discussion

According to traditional parameterizations of entrainment, it is commonly assumed that the entrainment/detrainment rate is proportional to (e.g., Simpson 1971), that is, decreased core size is related to increased entrainment. This is generally not consistent with the results presented herein because there is a sharp increase in entrainment and mixing between of 500 to 250 m, which is shown by the passive tracer analysis (more than an order of magnitude change in the occurrence of large low-level tracer values in the upper troposphere), yet only a small change in updraft core radius, which is shown by the mean updraft core area in Figs. 8b and 9b. There is some evidence for an increase in entrainment and mixing in convective cores between of 2 and 1 km (Fig. 13), which may be because of the decrease in the core radius for these grid spacings. However, we cannot definitively say that enhanced entrainment is caused by the change in core radius for between 1 and 2 km because there are other -dependent factors that also affect entrainment, notably computational and subgrid-scale mixing.

The abrupt changes in several of the convective characteristics between of 500 and 250 m, that is, mean convective updraft mass flux and tracer transport, suggest two distinct regimes. The sharp transition between these regimes does not appear to be directly associated with changes in the mean updraft core area (radius) or the kinetic energy spectra, which show generally smooth changes as is modified. A possible explanation for this behavior is a shift from laminar to turbulent flow between these two -dependent regimes, with laminar flow primarily occurring for larger than 250 m and turbulent flow for smaller , which we will refer to as CRM and LES regimes. The sharp increase in entrainment and mixing between these regimes based on the tracer analysis in particular is consistent with a shift to turbulent flow.

To illustrate this transition between regimes, line-parallel cross sections of the equivalent potential temperature are shown in Fig. 14. By comparing Figs. 14a and 14b ( 1 km) with Figs. 14c and 14d ( 100 m), a distinct change in the nature of the convective updrafts can be seen. For large , the flow consists of wide, straight, and continuous updrafts, whereas for small , there is substantial small-scale variability superimposed on the updraft structures and less coherency, especially in the vertical.

To explore the idea of a laminar-to-turbulent transition further, we refer to the scaling of the model turbulent Reynolds number in Bryan et al. (2003):

where *r* is the energy-containing large-eddy scale, and 1 is required for turbulent flow (i.e., *r *). We take *r* as the wavelength of the peak in the *w* kinetic energy spectrum, which is analogous to the mean convective core size as discussed above. As shown in Fig. 7, *r* decreases as decreases to 250 m but changes little as is decreased further. This implies that is approximately constant in the CRM regime ( greater than 250 m) but increases as is decreased in the LES regime. Qualitatively, this increase in is consistent with a transition to turbulent flow occurring between these regimes. Thus, resolving *r*, that is, resolving the peak in the *w* kinetic energy spectra or the mean updraft core size, seems to be the key for simulating a sufficiently large and transition to turbulent flow. In our simulations, *r* is resolved when it is approximately a factor of 10 larger than (or of approximately 250 m for the case simulated here), which is in agreement with the turbulence study of Wyngaard (2004). There is also little sensitivity of the mean convective characteristics with a further decrease in in the LES regime, consistent with the idea of Reynolds number similarity in turbulent flow. However, Bryan et al. (2003) suggested that a ratio of approximately 100 may be necessary for statistical convergence of mean deep convective characteristics based on LES studies of the planetary boundary layer. Differences between the findings shown here and those of, for example, Bryan et al. (2003) and Stein et al. (2015), could be related to the use of different models, cloud systems, and so on. Notably, Stein et al. (2015) found that the convective core width continued to change for decreasing from 200 to 100 m, although the 200-m simulation was found to best represent the observations. These simulations were performed with a different model and for a different case study, focusing on numerical weather prediction modeling. Thus, it is difficult to directly compare the results presented herein with those of Stein et al. (2015).

The grid spacing for which statistical convergence occurs depends somewhat on how one defines “convergence.” Here, we use a qualitative definition of convergence, that is, the point at which changes in the studied parameter, which is in this case, result in relatively small changes in the analyzed macroscale convective characteristics. For example, in Fig. 8d, the range defined by the mean plus and minus one standard deviation exhibits overlapping regions for from 2 km to 500 m and from 100 to 33.33 m. The range for of 250 m does not overlap with the range for any other . We treat the overlapping ranges as being indicative of convergence, while the lack of overlap for of 250 m suggests a shift in the storm structure. However, whether or not convergence is achieved depends on the characteristic of interest. For example, if one were to use the maximum updraft velocity, it is clear there are no systematic differences between the simulations performed with different . However, systematic differences between the simulations performed with different can be observed when looking at, for example, the average core area or the mean convective mass flux, which are aggregated quantities that are intended to represent the general structure and characteristics of convection.

## 5. Conclusions

A framework for performing numerical simulations of mesoscale systems with very small is presented and applied in the context of an idealized midlatitude continental squall line. This approach is based on performing a simulation with a relatively modest to spin up 3D fields, interpolating the model output to smaller (or larger) , adding small perturbations to the potential temperature field, and subsequently restarting the model at the chosen . It is shown that the restarted simulations require less than 1 h of simulation time to reach a quasi steady state (for convective characteristics and kinetic energy spectra), which is less than 25% of the time needed to spin up the same simulation using the chosen at the outset. Moreover, this framework allows for a direct comparison between simulations at different for mature systems while ensuring consistent initialization and early storm evolution.

Pronounced changes in convective characteristics due to decreasing were identified, especially when was changed from 500 to 250 m. This shift in the convective characteristics broadly agrees with the results of Bryan et al. (2003). Here, we have extended the findings of Bryan et al. (2003) with additional analysis and by exploring simulations with smaller . Interestingly, further decreasing below 250 m did not produce much change in the convective characteristics, such as the mean number and area of convective cores and mean convective mass flux, despite resolving a greater portion of the inertial subrange (with a nearly − slope in the kinetic energy spectra). Moreover, decreasing from 250 to 125 m had negligible effects on the structure of the simulated squall line and the horizontal grid spacing at which the structure of the system appeared to converge. This change in also had a limited impact on the kinetic energy spectra for scales greater than 6. These results suggest grid cell anisotropy may not be a critical factor, at least for this idealized squall line and for the range of and tested here.

A passive tracer analysis demonstrated that the shift in convective characteristics between of 500 and 250 m corresponded with a shift in vertical transport and entrainment/detrainment. For 500 m, the low-level tracer was more readily lofted to higher levels within convective cores even though the updrafts tended to be somewhat weaker. This lack of dilution was because of the suppression of entrainment/detrainment processes. However, for 500 m, there was a substantial decrease in the vertical transport of the low-level tracer because of the detrainment at low levels and entrainment at midlevels, suggesting that the cores were highly diluted by the time air parcels reached the upper troposphere. For 250 m, tracer transport was fairly insensitive to , even though the energy spectra clearly confirm that smaller-scale motions were better resolved with smaller . An increase in entrainment, detrainment, and mixing with a decrease in from 2 km to 250 m is not surprising given the decrease in mean updraft core size, meaning a decrease in the core volume to lateral surface area. However, this cannot explain the abrupt shift in the mean convective mass flux, entrainment, and mixing between of 250 and 500 m because the mean updraft core size decreased smoothly from of 2 km to 250 m. Instead, it was proposed that a transition from laminar to turbulent flow led to the shift in properties, consistent with an increase in the turbulent Reynolds number with a decrease in when the peak in the vertical velocity kinetic energy spectrum (or analogously the mean updraft core size) is resolved.

Overall, these results suggest that numerical simulations of moist deep convection have two -dependent regimes, one for large (CRM regime) and one for small (LES regime), with a sharp transition between the two. These two regimes can be summarized as follows:

CRM regime ( 500 m): The peak in the vertical velocity energy spectrum is not resolved. The mean convective core size decreases, while the number of cores increases for smaller . The mean convective updraft mass flux is relatively large and updrafts remain fairly undiluted, resulting in substantial transport of low-level tracer to the upper troposphere.

LES regime ( 250 m): The peak in the vertical velocity energy spectrum is resolved. The mean convective core size and number remain nearly constant for changes in . The mean convective updraft mass flux is relatively small and updraft cores are more diluted, resulting in less transport of low-level tracer to the upper troposphere than in the CRM regime.

These -dependent regimes are found to be insensitive to a change in by a factor of 2, although the inclusion of interactive surface fluxes may alter this finding.

The framework presented herein could be applied in the context of other convective systems, for example, supercells, and is amenable to real case studies, although it is possible that the convergence of the convective statistics may differ between storm types and storm environments. Moreover, this framework can facilitate future investigations into more detailed entrainment/detrainment-related effects and ultimately improved representations of entrainment/detrainment that can be implemented in coarser-scale models. This could include simulations performed at even smaller , employing and of *O*(10) m.

## Acknowledgments

The case presented here stems from a squall-line case study that was part of the Eighth International Cloud Modeling Workshop held in Warsaw, Poland. ZL was supported by the Advanced Study Program at the National Center for Atmospheric Research (NCAR) and the National Oceanic and Atmospheric Administration’s Climate Goal funding. HM was partially supported by U. S. DOE ASR DE-SC0008648, NASA NNX14AO85G, and the NSF Science and Technology Center for Multiscale Modeling of Atmospheric Processes (CMMAP), managed by Colorado State University under Cooperative Agreement ATM-0425247. We would also like to acknowledge high performance computing support from Yellowstone provided by NCAR’s Computational and Information Systems Laboratory, sponsored by the National Science Foundation.

## REFERENCES

*Cloud Dynamics.*International Geophysics Series, Vol. 53, Academic Press, 573 pp.

*Mesoscale Meteorology in Midlatitudes.*Wiley-Blackwell, 407 pp.

## Footnotes

Current affiliation: Department of Atmospheric Science, University of Wyoming, Laramie, Wyoming.

^{+}

The National Center for Atmospheric Research is sponsored by the National Science Foundation.

^{1}

More recent investigations have suggested that the scaling of entrainment is inversely proportional to height within a cloud (e.g., de Rooy and Siebesma 2008).