## Abstract

In this paper, the sensitivity of idealized squall-line simulations to horizontal resolution, subgrid turbulence mixing scheme, and subfilter-scale motion is discussed. Inconsistent results from numerical simulations of convective systems have suggested that there are issues with the behavior of the subgrid turbulent mixing parameterizations with increasing resolution that still need to be understood. WRF is used to perform large-eddy simulation of an idealized squall line with horizontal grid spacings of 4 km, 2 km, 1 km, 500 m, and 250 m. While 4 km grid spacing is able to produce the general structure of the squall line, higher-resolution simulations produce more detailed structures. Individual convective cell size decreases, the maximum cloud top height increases, and the subgrid turbulence kinetic energy (TKE) ratio decreases as resolution increases. As found in past studies, 4 km grid spacing is not recommended as it contains an unreasonable amount of subgrid TKE, is not sufficient to resolve the large energy-containing eddies, and may even suppress propagation of the squall line. While horizontal resolution of 1 km can produce a squall line, there are several discrepancies between the 1 km case and higher resolutions, including trailing banded structures and inhibited three-dimensionalization. These issues at 1 km resolution are investigated by examining the subfilter energy transfer for the highest-resolution simulation filtered to a horizontal scale of 1 km. The subfilter energy transfer rate at a scale of 1 km is dominated by the streamwise and shear components. While dissipation dominates the transfer, a significant amount of backscatter also exists, which is not represented by most subgrid models.

## 1. Introduction

Convective systems have been studied for many years, and great progress has been made in understanding them (e.g., Houze 2004). Numerical simulation is an important tool for testing hypotheses about convective systems. Numerical simulations of idealized convective systems, such as squall lines, are useful for understanding basic physical processes and testing different numerical approaches (e.g., Lilly 1979; Rotunno et al. 1988).

Early successes in storm simulations (e.g., Lilly 1979) have inspired research in different aspects of cloud resolving model (CRM) and large eddy simulation (LES) of mesoscale convective systems, including squall lines, such as the effect of horizontal and vertical resolution (e.g., Bryan et al. 2003; Lebo and Morrison 2015), environmental instability (e.g., Potvin and Flora 2015; Morrison et al. 2015), microphysics parameterization (e.g., Bryan and Morrison 2012; Morrison et al. 2015), vertical wind shear (e.g., Adlerman and Droegemeier 2005), subgrid turbulence mixing scheme (e.g., Verrelle et al. 2015, 2017), cold pool dynamics (e.g., Schlemmer and Hohenegger 2014), environmental humidity (e.g., James and Markowski 2010), and surface conditions (e.g., Peters and Hohenegger 2017) on storm simulations.

The impact of resolution on simulations of mesoscale convective systems is significant (e.g., Weisman et al. 1997; Bélair and Mailhot 2001; Bryan et al. 2003; Fiori et al. 2010; Bryan and Morrison 2012; Lebo and Morrison 2015; Verrelle et al. 2015). Low resolution can lead to a significant delay in storm development (e.g., Weisman et al. 1997; Bryan et al. 2003; Bryan and Morrison 2012). Subkilometer horizontal resolution is necessary for capturing the convective structures of a convective system (e.g., Bryan and Morrison 2012; Verrelle et al. 2015). These overall conclusions have been demonstrated in many studies, which are reviewed below.

Weisman et al. (1997) investigated the effect of horizontal resolution on convective systems using a nonhydrostatic cloud model with horizontal grid sizes from 12 to 1 km. They concluded that 4 km grid spacing was sufficient to resolve the mesoscale structure of the convective system that was produced in the 1 km simulation. Bryan et al. (2003) conducted squall-line simulations with four different grid spacings from 1 km to 125 m. As resolution increased, the results did not converge. Their results suggested that simulations with 1 km grid spacing contained a large amount of subgrid turbulence energy. As a result, simulations with 1 km grid spacing should not be used as a benchmark. Instead, simulations with grid spacings of *O*(100) m should be used.

Bryan and Morrison (2012) conducted idealized squall-line simulations with horizontal resolutions of 4 km, 1 km, and 250 m and different microphysics parameterizations. They showed that the dependence on horizontal resolution was more significant than that of microphysics parameterization. The cloud top height in the 1 km run was higher than those in the 4 km and 250 m runs, suggesting that there is a competition between resolved and subgrid turbulence mixing, which is consistent with other studies (e.g., Waite and Khouider 2010). In a similar study, Morrison et al. (2015) compared idealized deep convective storm simulations with horizontal grid spacing from 2 km to 125 m, three different environmental soundings, and three different microphysics parameterizations. They showed that varying horizontal resolution does not strongly affect precipitation, evaporation and condensation. They also suggested that differences in cloud water evaporation and temporally averaged surface precipitation with Bryan and Morrison (2012) may be related to the sounding profile and wind shear. Lebo and Morrison (2015) conducted idealized simulations of squall lines with seven different horizontal resolutions from 2 km to 33.33 m and vertical resolutions of 500 and 250 m. They found that decreasing the horizontal grid spacing below 250 m did not change convective characteristics such as the mean number and area of convective cores. While the vertical resolution has an effect on simulations of mesoscale convective systems, some studies have found that simulations are qualitatively similar for minimum vertical resolutions of 200 to 50 m (e.g., Adlerman and Droegemeier 2002; Lebo and Morrison 2015; Waite 2016). In particular, Adlerman and Droegemeier (2002) simulated a cyclic mesocyclone and found that decreasing the minimum vertical resolution from 200 to 50 m sped up the cycling process, but that the simulated storms were qualitatively similar. Since our minimum vertical grid spacing is similar, the dependence on vertical resolution will not be discussed in this paper.

In addition to examining the dependence of simulated storm characteristics on grid spacing, some studies have considered the behavior of subgrid turbulence parameterizations at different resolutions. Verrelle et al. (2015) conducted multicellular storm system simulations with grid spacings of 4 km, 2 km, 1 km, and 500 m. They showed that the ratio of subgrid turbulence kinetic energy (TKE) to the total TKE increases as the resolution increases. This finding is consistent with earlier work (e.g., Adlerman and Droegemeier 2002). However, Bryan et al. (2003), who conducted squall-line simulations, found that the ratio of subgrid turbulence kinetic energy to the total turbulence kinetic energy decreases as the resolution increases. Similarly, Potvin and Flora (2015) conducted idealized supercell simulations with horizontal resolutions from 4 km to 333 m. They suggested that horizontal resolution of 1 km could produce a useful forecast, but benefits could be seen when horizontal resolution was further increased to 333 m. Verrelle et al. (2017) conducted LES of a deep convective cloud with 50 m horizontal grid spacing, which they then filtered to 500 m, 1 km, and 2 km. They suggested that thermal effects dominated at resolutions of 2 km and 1 km and dynamical effects dominated at a grid spacing of 500 m or smaller.

The effect of subgrid turbulence mixing schemes on simulations of mesoscale convective systems and general circulation models is significant (e.g., Takemi and Rotunno 2003; Moeng et al. 2010; Kirkil et al. 2012; Schaefer-Rolffs and Becker 2013). Moeng et al. (2010) performed a test with a mixed subgrid turbulence scheme on a tropical deep convection system. The representation of subgrid heat, moisture and momentum flux improved with the mixed subgrid scheme at horizontal resolution of a few kilometers. Kirkil et al. (2012) simulated neutral atmospheric boundary layer flows over flat terrain and a ridge with Smagorinsky, Lagrangian-averaged scale-dependent dynamic, nonlinear backscatter, and anisotropy, and dynamic reconstruction model mixing schemes. The Lagrangian-averaged scale-dependent dynamic and dynamic reconstruction models were found to be more accurate than the Smagorinsky mixing and the nonlinear backscatter models, but also took 25%–28% more CPU time than the Smagorinsky mixing scheme to compute.

The LES studies discussed above primarily used eddy viscosity to represent subgrid turbulence. Eddy viscosity dissipates resolved TKE, consistent with the notion of a downscale energy cascade (e.g., Wyngaard 2010). However, while the energy cascade hypothesis is true on average, it may not always be true locally. Indeed, the transfer of subfilter TKE to resolved scales, known as backscatter, can be important locally in many turbulent flows, including isotropic turbulence, channel flows, and stratified turbulence (e.g., Germano et al. 1991; Piomelli et al. 1991; Domaradzki et al. 1993; Khani and Waite 2013, 2016). The magnitude of forward (dissipation) and backward (backscatter) subfilter energy transfer are often similar (e.g., Germano et al. 1991; Piomelli et al. 1991; Schaefer-Rolffs 2018). Backscatter is not represented in eddy viscosity models, and therefore may account for some of the dependence on grid size discussed above.

In general, there is agreement among these studies that subkilometer resolution is necessary in order to properly resolve mesoscale convective systems, like squall lines (e.g., Bryan et al. 2003; Bryan and Morrison 2012; Verrelle et al. 2015). However, the role of backscatter and the effect of subgrid dissipation on idealized squall-line simulations is still unclear. Previous studies show that the evolution of the squall line can be significantly affected by the representation of subgrid motions. In this work, we further investigate the problems that emerge at low resolution, and the dependence of the results on resolution and subgrid turbulent mixing parameterization, in the context of an idealized squall-line simulation. In addition, and unlike previous studies on this topic, we quantify some of the missing dynamics in the low-resolution simulations by filtering high-resolution output and examining subfilter energy transfers. To better understand the dependence on resolution and mixing scheme, the relationship between the structure of the squall line and subfilter energy transfer rate will be examined.

We address these problems by simulating an idealized squall line with five different horizontal resolutions and three different mixing schemes. In particular, we investigate the sensitivity of idealized squall-line simulations to horizontal grid spacings, turbulent mixing parameterization and subgrid-scale motions. The rest of this paper is organized as follows. In section 2, we present the methodology of this paper. Results and discussion are presented in section 3. In section 3a, the main focus is on the resolution dependence of simulations with TKE mixing. In section 3b, results with the Smagorinsky scheme and no mixing scheme are presented. In section 3c, we discuss the subfilter energy transfer rate. Conclusions are given in section 4.

## 2. Methodology

### a. Experimental setup

The Advanced Research Weather Research and Forecasting (WRF) Model is used in this paper. WRF solves the fully compressible, nonhydrostatic Euler equations using terrain-following hydrostatic pressure as the vertical coordinate (Skamarock et al. 2008). It uses a third-order Runge–Kutta scheme in time, third-order advection scheme in the vertical, fifth-order advection scheme in the horizontal to solve the equations, and a positive-definite advection scheme for scalar variables. The spatial discretization is upwind-biased finite difference on a Arakawa C grid. The setup of the simulation is similar to Bryan and Morrison (2012) except for the initial temperature and moisture profiles. The domain size is 576 km in the streamwise (*x*) direction, 144 km in the spanwise (*y*) direction, and 25 km in the vertical. In all simulations, the vertical grid has 100 levels with nonuniform grid spacing from 220 m at the bottom to 900 m at the top of the domain. The horizontal resolution is Δ*x* = Δ*y* = 4 km, 2 km, 1 km, 500 m, and 250 m. Three subgrid turbulence approaches are used: the prognostic TKE closure (Lilly 1967), the 3D Smagorinsky closure (Smagorinsky 1963) and no subgrid turbulence scheme (Δ*x* = 2 km is only run for TKE and Smagorinsky mixing schemes). All simulations are run for 6 h. A time step of 3 s is used for the Δ*x* = 4 km, 2 km, 1 km, and 500 m simulations, and 1 s is used for the Δ*x* = 250 m simulations.

The Weisman and Klemp (1982) sounding profile is used to generate initial conditions (as in Rotunno et al. 1988; Weisman et al. 1997). The initial cold pool in this paper is a 2 km deep region on the west half of the domain, with a potential temperature perturbation of −8 K set on the surface, which increases linearly to 0 K at 2 km above the surface as in Bryan and Morrison (2012). Random temperature perturbations between ±0.008 K are added to the lowest grid level. Open boundary conditions are used in the *x* direction, and periodic boundary conditions are used in the *y* direction. The initial wind shear, potential temperature and water vapor mixing ratio profiles are shown in Fig. 1.

The Lin et al. (1983) microphysics scheme is used, which is a single-moment scheme with six species of water: vapor, liquid cloud water, ice cloud water, rainwater, snow, and graupel. The drag coefficient is set to be 0.01, and radiative cooling is turned off. Surface heat and moisture fluxes are not included. Rayleigh damping is used in the top 5 km of the domain (from 20 to 25 km, not shown in the following figures) with a damping coefficient set to 0.003.

### b. Diagnostics

#### 1) Intensity

In the following section, several diagnostics are employed to study the squall-line evolution; these diagnostics are briefly introduced here.

The streamwise, spanwise, and vertical velocity are *u*, *υ*, and *w*, respectively. The speed of the squall line is determined by tracking the *x* position of the point with the maximum vertical velocity over the whole domain. The mean propagation speed is calculated for the mature squall line from 200 to 350 min. Mass flux is calculated by the domain sum of *ρ*max(*w*, 0) Δ*x*Δ*y* (e.g., Bryan and Morrison 2012), where *ρ* is the density and *w* is the vertical velocity. While this quantity does not have the units of flux, the term mass flux is used for consistency with previous studies (e.g., Bryan and Morrison 2012). Cloud top height is calculated by finding the maximum height with the horizontal-averaged cloud ice mixing ratio greater than 0.0005 g kg^{−1}. Cold pool intensity *C* is calculated by $C2=\u222b0hB\u2009dz$, where *B* is the buoyancy and *h* is the height of the cold pool, which are defined as in Bryan and Morrison (2012).

To quantify how three-dimensional the simulations are, we introduce ratios of *u*′ and *w*′ to ⟨*u*⟩ and ⟨*w*⟩. Here, ⟨⋅⟩ represents spanwise average and prime denotes fluctuation from the spanwise averages. Specifically, we consider the ratio of the horizontal average of $u\u20322$ and $w\u20322$ to the horizontal average of ⟨*u*⟩^{2} and ⟨*w*⟩^{2}, respectively, at a height of 2 km. If the squall line is primarily two-dimensional, the fluctuation from the spanwise averages will be small compared to the spanwise averages and this ratio will be small.

#### 2) Vertical velocity spectra

The vertical velocity spectra are calculated following Bryan et al. (2003). At height *z* = 5 km, the *y*-average (spanwise-averaged) vertical velocity is computed, and the point with the maximum fluctuation from this average is computed. A one-dimensional *y* spectrum is computed of *w* through this point. Then the spectrum is averaged from *t* = 120 to 270 min.

#### 3) Subgrid TKE ratio

The subgrid TKE ratio *R* is calculated as

where *e* is the subgrid TKE from the LES scheme and *E* is the resolved TKE, defined following Bryan et al. (2003) as

#### 4) Subfilter energy transfer rate

In this work, filtered velocity fields are computed with a low-pass horizontal box filter with scale of Δ = 1 km; for example, for *u*

where the integral is computed with the trapezoidal scheme. Since the filtered fields are computed to investigate the dynamics of different grid spacing, and only horizontal resolution is varied, we do not filter fields in the vertical.

The energy transfer rate *P* from filtered to subfilter length scales in a volume *V* is defined following Piomelli et al. (1991) as

where *ρ* is the density. Below we will compute *P* on a horizontal layer at constant height, in which case *V* is the volume of that layer. Here

is the filtered rate of strain tensor, and the effective subfilter stress tensor is

This term is normally parameterized in LES, but is computed explicitly here from the Δ*x* = 250 m simulation, filtered to a horizontal scale of 1 km. Later in the discussion section, $\u2061(\u22c5)\u02dc$ and $\u2061(\u22c5)\xaf$ will be dropped in the $\tau ij\u02dc$ and $Sij\xaf$ terms.

Furthermore, *P* can be separated into negative (dissipation) and positive (backscatter) contributions, where the dissipation *ε*_{−} and the backscatter *ε*_{+} are defined as in Piomelli et al. (1991) as

The backscatter represents energy transfer from subfilter to filtered scales, and is neglected in eddy viscosity LES models.

To understand the relationship between *P* and the structure of a squall line, we further decompose *P* into six independent contributions

at *t* = 240 min, where *i*, *j* = 1, 2, 3. Note that repeating indices are not summed on the rhs of (3) (i.e., *P*_{13} = *τ*_{13}*S*_{13}).

## 3. Results

### a. Sensitivity to resolution

The squall line has a main updraft at the leading (right) side with clouds trailing behind the main updraft. The cloud top height is approximately 10 km and the propagation speed of the squall line is approximately 10 m s^{−1}. The squall line takes about 90 min to become fully developed. We will focus in particular on the results at *t* = 240 min, at which the flow is mature and not interacting with the boundary on the right. Results at other times in the mature evolution are similar.

#### 1) Structure

The spanwise-averaged cloud water plus ice mixing ratio shows that the TKE mixing scheme needs horizontal resolutions higher than 1 km to correctly resolve the main features of a squall line (Fig. 2). Although a horizontal resolution of 4 km can produce the general structure of the squall line with tilting and trailing stratiform clouds, which can be seen in Fig. 2a, this grid can only produce a large convective plume with a large trailing cloud, without resolving the details of the convective structure in the squall line. When the resolution increases from Δ*x* = 4 to 2 km (Fig. 2b), the effect on the cloud water and ice mixing ratio is significant: the simulation with Δ*x* = 2 km shows more details and discrete structures than the simulation with Δ*x* = 4 km. The Δ*x* = 2 km simulation produces banded structures in the trailing clouds. When the resolution increases from Δ*x* = 2 to 1 km (Fig. 2c), each convective cell in the trailing cloud becomes smaller, and the banded structures are still present. These banded structures in the trailing clouds do not appear at higher resolution. As a result, these structures seem to be an artifact of low resolution and marginally resolved convection. In fact, they disappear when a nonoscillatory advection scheme is used (not shown). The overall convective structures are similar in the simulation with Δ*x* = 2 and 1 km.

As the resolution further increases from Δ*x* = 1 km to Δ*x* = 500 m, the size of the trailing stratiform region decreases, which suggests that only the 1, 2, and 4 km simulations have long trailing clouds. The stratiform region in the simulation with Δ*x* = 500 m is significantly smaller than the lower-resolution simulations, and the peculiar banded structures found with the Δ*x* = 2 and 1 km simulations are absent. At the main updraft core of the squall line, the cloud water and ice mixing ratio is the largest when Δ*x* = 500 m. The cloud height also increases when resolution increases from Δ*x* = 1 km to Δ*x* = 500 m. Such effects may be due to insufficient mixing in the main updraft area. The further increase of the resolution from Δ*x* = 500 to 250 m (Fig. 2d) has a significant effect on the structure of the storm. There is less cloud water in the main updraft core and a larger trailing cloud in the Δ*x* = 250 m case. In general, squall lines simulated at higher resolution show more detailed structures. Lower resolution yields a larger and more attached structure, while a higher resolution gives a smaller and more separated structure. Such dependence of the individual cell size on resolution has also been reported by some studies such as Bryan et al. (2003) and Bryan and Morrison (2012).

To investigate the impact of resolution on the intensity of convection, we consider the spanwise-averaged vertical velocity (Fig. 3) and the domain-maximum vertical velocity (Table 1). The size and number of individual updraft cores are sensitive to the resolution. Lower resolution results in one large updraft core, while high resolution gives more individual and smaller updraft structures. Such dependence of the updraft core area on the resolution has also been shown by Lebo and Morrison (2015), Morrison (2016), and Lebo (2018). We can also see that the maximum vertical velocity approaches 60 m s^{−1} as resolution increases from Δ*x* = 4 km to Δ*x* = 250 m. Similar dependence of the maximum vertical velocity on the resolution has also been reported by Weisman et al. (1997), Morrison (2016), and Lebo (2018), which suggest that the maximum vertical velocity may depend on the updraft width.

To illustrate the dependence of horizontal structure on Δ*x*, we consider horizontal slices of potential temperature at a height of 2 km at different resolutions in Fig. 4. All simulations at different resolutions have a warm leading edge and a cold turbulent region trailing behind. The two-dimensional banded structures behind the leading edge do not exist after the storm has fully developed in cases with Δ*x* ≤ 500 m. As seen in the cloud water slices in Fig. 2, the Δ*x* = 1 and 2 km case show the banded structure behind the leading edge, which is not shown by the simulation with higher or lower resolutions. The banded structures appear to be related to the delayed three-dimensionalization, and they become turbulent at a later time in the simulation. For the higher-resolution simulations, the three-dimensionalization process is faster.^{1} Hence, the banded structures do not persist in simulations with Δ*x* = 500 and 250 m. In addition, the leading edge of the squall line oscillates slightly from north to south, and this meandering is more pronounced at high resolution. The meandering of the leading edge can be observed with grid spacing smaller than Δ*x* = 1 km. As a result, we believe that the meandering is related to the instability and three-dimensionalization of the turbulence in the warm leading edge. Such meandering can also be found in Bryan et al. (2003).

Furthermore, the structure of the warm leading edge depends on resolution. At Δ*x* = 500 m, a solid thick warm edge can be seen, which is not seen with Δ*x* = 250 m. The area of the solid warm edge with Δ*x* = 500 m is similar to the area of the line structure with Δ*x* = 1 km. More details in the structure can be seen in the Δ*x* = 500 m case than the 1 km case at earlier times. The warm leading edge of the Δ*x* = 250 case is slightly colder than the leading edge in the Δ*x* = 500 m case.

To illustrate the dependence of three-dimensionalization on Δ*x*, time series of the three-dimensionalization ratio for *u* and *w* are shown in Fig. 5. Higher-resolution simulations three-dimensionalize more quickly: for both *u* and *w*, as the resolution increases, the time needed for the squall line to three-dimensionalize decreases. Before the squall line matures (*t* ≈ 90 min), the ratios increase as resolution increases. While the ratios for the high-resolution simulations (250 and 500 m) still dominate most of the time after the squall line matures, especially for *w*, the resolution dependence is less obvious at later times.

#### 2) Intensity

According to the magnitude of the mass flux (Fig. 6a), a clear rank on the intensity of the squall line with resolution can be given after the squall line is mature: in ascending order of the magnitude of the mass flux (from least to most): Δ*x* = 4 km, 250 m, 2 km, 1 km, and 500 m. The large jump in mass flux as the grid spacing changes from Δ*x* = 4 to 2 km is likely due to the effect of increased resolution on nonhydrostatic processes as suggested by Weisman et al. (1997). The magnitude of the mass flux for the Δ*x* = 2 and 1 km cases are similar in the first half of the simulation. The magnitude of the mass flux for the Δ*x* = 2 km case drops significantly at *t* = 180–320 min and rebounds at *t* = 320 min. As the resolution increases from Δ*x* = 1 km to 500 m, the mass flux continues to increase significantly. As the grid is further refined to Δ*x* = 250 m, the mass flux drops by about 20%. This nonmonotonic behavior is likely due to the competition between resolved and unresolved processes.

No clear trend in precipitation rate can be identified as the resolution increases (Fig. 6d). Nevertheless, some interesting dependence of the precipitation rate on Δ*x* can be seen. For example, the precipitation peaks at different times for different Δ*x*. The precipitation rate peaks earlier (around *t* = 80 min) in the Δ*x* = 1 km case than in the other cases (around *t* = 100–120 min). The precipitation rate peak for the Δ*x* = 1 km case occurs during the spinup. Moreover, the precipitation rates decrease after 200 min in all cases. The maximum values of the precipitation rates are similar in all cases and they do not depend much on Δ*x*. The Δ*x* = 250 m case has the highest precipitation rate by a small margin. While there is no dependence of precipitation on resolution in our study, Bryan and Morrison (2012) found that increased resolution resulted in decreased precipitation. This discrepancy may be due to the different initial sounding used in Bryan and Morrison (2012).

The maximum cloud top height increases slightly with increasing resolution (see Fig. 2 and Table 1). Past studies have found a more significant impact of resolution on cloud height (e.g., Waite and Khouider 2010; Bryan and Morrison 2012). While our results show that cloud top height increases as resolution increases, Bryan and Morrison (2012) found that Δ*x* = 1 km had the highest cloud top when compared to Δ*x* = 4 km and Δ*x* = 250 m in their idealized squall-line simulations.

The mean propagation speed (Table 1) increases with increasing resolution: for the TKE scheme, the speed increases from 8.9 m s^{−1} at Δ*x* = 4 km to 11.6 m s^{−1} for Δ*x* = 250 m. According to RKW theory (Rotunno et al. 1988), cold pools are important for the development and propagation of convective cells, and a stronger cold pool may create a stronger squall line. As a result, the dependence of the mean propagation speed on resolution may be due to the cold pool intensity, which (at least in the early mature phase) increases as Δ*x* decreases from 4 km to 500 m, after which it approximately converges, like the propagation speed (Fig. 7a).

#### 3) Energy

To investigate how the subgrid TKE scheme performs at each resolution, we consider the ratio of subgrid to total TKE (Fig. 8). Subgrid TKE is mainly localized in the updraft portion of the storm. At each resolution, the highest subgrid TKE ratio occurs at the leading edge. Moreover, the subgrid TKE ratio is significantly higher in the Δ*x* = 2 km case. The subgrid TKE ratios are larger in the Δ*x* = 4, 2, and 1 km cases than the Δ*x* = 500 and 250 m cases. At low resolution (4 and 2 km), the subgrid turbulence scheme is creating more subgrid TKE to compensate for the insufficiently resolved turbulence, as expected. However, the dependence is not monotonic: as the resolution increases from Δ*x* = 4 km to Δ*x* = 2 km, the subgrid TKE ratio increases, while as resolution increases from Δ*x* = 2 km to Δ*x* = 250 m, the subgrid TKE ratio decreases. Note that Bryan et al. (2003) found that as resolution increases from 1 km to 125 m, the subgrid TKE ratio decreases. Overall, the largest subgrid TKE ratio at each resolution occurs at a height below 5 km. More turbulence occurs at these heights due to the interaction between the environmental wind shear and the main updraft.

To investigate the distribution of resolved TKE at different length scales, the vertical velocity spectra in *y*, averaged from *t* = 120 to 270 min, are plotted for different resolutions in Fig. 9. The time interval is chosen after the storm is fully developed. The method for obtaining the spectra is described in section 2b. While the spectra are very noisy, the high-resolution simulations show a broad maximum around *k* = 20–40, corresponding to wavelengths of a few kilometers. These are the length scales of the most energetic convective plumes. The resolution of Δ*x* = 4 km is clearly not sufficient for resolving these scales. Even Δ*x* = 1 km only marginally captures the peak, but the spectrum is very steep beyond *k* = 30 due to eddy dissipation, which is also been noted by Bryan et al. (2003). The spectra for 500 and 250 m are consistent with the Kolmogorov −5/3 spectrum between a narrow range of wavenumbers from *k* = 20 to 60. These spectra, along with the subgrid TKE ratio, suggest that a grid spacing of 500 m or less should be used for squall-line simulations.

### b. Sensitivity to mixing scheme

The above section focuses on the dependence of idealized squall-line simulations on resolution when the TKE scheme is used. In this section, we consider the dependence on subgrid mixing scheme. Additional simulations with Smagorinsky mixing and with no subgrid scheme (i.e., numerical diffusion only) are discussed here.

At a fixed resolution, the squall line is sensitive to subgrid mixing scheme. Simulations with different schemes are able to produce a trailing stratiform squall line within the 6 h simulation time. From the mass flux (Figs. 6a–c) and precipitation rate (Figs. 6d–f), one can see that the squall line in the simulation with Δ*x* = 4 km and Smagorinsky mixing is significantly delayed. Indeed, the Smagorinsky simulation with Δ*x* = 4 km is not able to generate sufficient convection to initiate the squall line during the first half of the simulation time. By contrast, simulations with TKE and no mixing scheme are able to create a propagating squall line at the lowest resolution. The Smagorinsky scheme is known to be too dissipative at low resolution (e.g., Pope 2000), which explains its poor performance at Δ*x* = 4 km. Indeed, simulations with the Smagorinsky scheme can only produce the correct shape of the squall line when Δ*x* ≤ 1 km, which is shown in the spanwise-averaged cloud water plus ice mixing ratio field (Fig. 10). This shows that Smagorinsky mixing should not be used with Δ*x* = 4 km for simulations of convective systems with explicitly resolved convection.

At a fixed resolution, the Smagorinsky mixing simulations have more cloud water and ice in the trailing region behind the main core compared to the other mixing schemes (cf. Figs. 2, 10, 11). On the other hand, the no-mixing simulations have more cloud water and ice in the main convective core, due to the reduced mixing between the core and surroundings. Overall, the squall line is more sensitive to resolution than mixing scheme (except for the coarse Smagorinsky simulation).

Simulations with different mixing schemes have different cloud top heights (Table 1). Except for simulations with Δ*x* = 4 km, the no-mixing simulations have the highest maximum cloud top, the Smagorinsky mixing simulations have the second highest maximum cloud top, and the TKE mixing simulations have the lowest maximum cloud top. The fact that cloud top height is affected by both the resolution and mixing scheme is interesting given the changes of cloud top height with resolutions found in Bryan and Morrison (2012), which implies that cloud height can sometimes depend on subgrid mixing, which is also shown in Waite and Khouider (2010). These results suggest that the dependence of the cloud top height on the resolution and subgrid mixing scheme depends on the environmental condition that the squall line is present in, such as the high CAPE environmental setting in Bryan and Morrison (2012).

The nonmonotonic dependence of mass flux on resolution described in the previous subsection for the TKE scheme also occurs with the Smagorinsky mixing scheme (Fig. 6b). Note that there is almost no mass flux for the first half of the simulation in the Δ*x* = 4 km simulation with the Smagorinsky mixing scheme. With no mixing, the Δ*x* = 4 km simulation has more mass flux than the Δ*x* = 250 m simulations, which is different from the TKE and Smagorinsky cases. On the other hand, the precipitation rate (Figs. 6d–f) is not sensitive to resolution for all schemes, except for the Smagorinsky scheme with Δ*x* = 4 km, for which there is no precipitation in the first half of the simulation.

The propagation speed (Table 1) also depends on the mixing scheme. For high-resolution simulations with Δ*x* = 250 and 500 m, storms with TKE mixing propagate the fastest, the storms with Smagorinsky mixing propagates the second fastest, and those with no mixing propagates the slowest. For low-resolution simulations with Δ*x* = 1 and 4 km, storms with TKE mixing propagate the fastest, the storms with no mixing propagates the second fastest, and those with Smagorinsky mixing propagates the slowest. The propagation speed of the storm may be related to the tilting of the squall line and the cold pool intensity of the squall line. A more tilted squall line propagates faster than a less tilted squall line. Also, the cold pool intensity for TKE simulations is the largest (Fig. 7). According to RKW theory (Rotunno et al. 1988), the tilting will affect the circulation created by the downdraft due to precipitation, which will affect the cold pool intensity. Hence, the tilting of the squall line may modify the propagation speed. With a fixed mixing scheme, the convergence of propagation speed as resolution increases can be seen.

### c. Subfilter energy transfer rate

While a horizontal resolution of 1 km can produce a squall line, the results from the previous sections indicate significant resolution dependence as Δ*x* decreases from 1 km to 500 m, since Δ*x* = 1 km is not small enough to resolve all the convective structures inside a squall line. There are still several discrepancies between the Δ*x* = 1 km case and the higher-resolution cases, such as the fact that the Δ*x* = 1 km case has issues with the trailing banded structures and three-dimensionalization (Fig. 2). Here, we further investigate the effects of these sub-1-km motions using the subfilter energy transfer rate.

Figure 12 shows the time series of the net subfilter energy transfer *P*, calculated from the Δ*x* = 250 m simulation on a horizontal slice at *z* = 2.5 km, filtered in the horizontal to 1 km, along with the subgrid eddy viscosity dissipation from the Δ*x* = 1 km simulation. The evolution of the net subfilter energy transfer rate (black) is similar to the subgrid dissipation from the Δ*x* = 1 km simulation (green). The peak of the subfilter energy transfer curve is approximately 40 min ahead of the subgrid dissipation curve, likely because a horizontal resolution of 1 km is not able to resolve the small-horizontal-scale convective motion at the early stage of the squall line.

The subfilter energy transfer is further decomposed in Fig. 12 into backscatter (red) and dissipation (blue). The subgrid dissipation is decomposed into horizontal (orange) and vertical (yellow) eddy viscosity contributions. After the early developing stage of the squall line (first 100 min), the backscatter is approximately constant, while the dissipation increases to *t* = 120 min and then decreases. The constant backscatter is approximately 15% of the peak dissipation and 50% of the dissipation at the end of the simulation. The dominance of the dissipation implies that the net energy transfer across a horizontal scale of 1 km is downscale, as expected, which is broadly consistent with an eddy viscosity approach; however, the eddy viscosity cannot represent the smaller but significant backscatter that is found. Backscatter may be especially important locally (e.g., Piomelli et al. 1991), and its structure may point to regions where the eddy viscosity model is particularly problematic. While the simulation with Δ*x* = 250 m does not resolve all sub-1-km turbulence, and therefore does not give the full sub-1-km energy transfer rate, it should give most of the dissipation and backscatter if the energy transfers are dominated by turbulence and convective features that are only slightly smaller than 1 km.

The vertical and horizontal eddy dissipation are similar at the early developing stage of the squall line. Both vertical and horizontal eddy dissipation increase significantly when the squall line is mature, which is due to the three-dimensionalization of the mature squall line. The vertical eddy dissipation is double that of the horizontal at its peak. The eddy dissipation is dominated by the vertical eddy viscosity after the squall line has developed.

We further investigate the subfilter energy transfer rate *P* by considering its spanwise average. As we can see from Fig. 13g, the spanwise-averaged *P* is almost entirely dissipation, which is shown in brown. At the low-level leading edge, there is a small patch of backscatter, which is shown in blue. Some backscatter is also visible at the very low levels behind the main updraft (Fig. 13g). Since *P* has the same tilted structure as the squall line, the subfilter energy transfer rate does not appear to be random, despite the turbulent nature of the flow, and the sign of the subfilter energy transfer rate may be related to the structure and tilting of the squall line.

From Fig. 13, we can see that the spanwise-averaged *P*_{11} and *P*_{13} have the largest contributions to *P*; *P*_{11} exhibits extensive dissipation in the lower region and slightly less backscatter in the upper region of the squall line, while *P*_{13} is entirely dominated by strong dissipation except for the strong backscatter at the leading edge. Moreover, *P*_{11} and *P*_{13} show the general structure and tilting of the squall line. Note that while *P*_{33} is also strongly dissipative, the shape of the squall line cannot be seen easily in this field. On the other hand, the contributions from *P*_{12}, *P*_{22}, and *P*_{23} are significantly weaker and less organized than the contribution from *P*_{11}, *P*_{13}, and *P*_{33}. While the squall line is not homogeneous in the spanwise (*y*) direction, the variation in the spanwise direction is smaller than in the streamwise and vertical directions. As a result, the spanwise components of the subfilter energy transfer rate are weaker. From the above discussion, we can conclude that the streamwise (*i*, *j*) = (1, 1) and shear (*i*, *j*) = (1, 3) components are stronger and more organized than the other components. Note that while the overall *P* decreases in time (Fig. 12), the structure of the spanwise-averaged *P*_{ij} and the locations of the positive and negative transfers are similar at earlier and later times in the squall-line evolution.

To better understand the structure of *P*_{11} and *P*_{13}, we consider spanwise-averaged and horizontal slices of these fields, along with the corresponding components of *τ*_{ij} (subgrid-scale momentum flux) and *S*_{ij} (resolved strain rate tensor). Spanwise-averaged and horizontal slices of *P*_{11}, *τ*_{11}, and *S*_{11} are shown in Figs. 14 and 15, respectively. The corresponding (*i*, *j*) = (1, 3) fields (i.e., *τ*_{13} and *S*_{13}) are shown in Figs. 16 and 17 .

Consider first the (*i*, *j*) = (1, 1) components. While the spanwise-averaged *P* (Fig. 13g) is mainly dissipative, *P*_{11} (Fig. 14a) is mainly dominated by backscatter in the upper region. Although there are small amounts of backscatter at the leading edge of the squall line at *z* = 2.5 km, they are small compared to the dissipation in the main updraft (Fig. 15a). Localized weak dissipation is also apparent in the trailing turbulent region. This leading-edge backscatter becomes more significant at higher levels (Fig. 14a), where positive *P*_{11} dominates, as seen in the spanwise average. Note that two distinct regions can be seen in both *τ*_{11} and *S*_{11} at *z* = 2.5 km: there is a stronger region at the front, and a weaker region trailing behind. The stronger region results in stronger dissipation and backscatter mixed at the front of *P*_{11} and the weaker region results in very weak dissipation in the trailing region. While the upper region of *P*_{11} is dominated by backscatter, the lower region is dominated by dissipation. While the spanwise-averaged structure suggests that the lower region of *P*_{11} is dissipative, it is the net result of lots of dissipation and backscatter, which can be seen at *z* = 2.5 km (Fig. 15a).

Both the spanwise-averaged and horizontal slices of *τ*_{11} (Figs. 14b, 15b) indicate that this component of the subfilter stress is mainly positive, implying a subgrid flux of positive zonal momentum to the right and negative momentum to the left. Therefore, the sign of *P*_{11}, which indicates whether these subfilter fluxes contribute to dissipation or backscatter, is determined by *S*_{11} (Figs. 14c, 15c). This component of the filtered strain rate is negative throughout the main updraft, with positive values to the right of the main updraft at mid and upper levels, and at low levels behind the main updraft. The regions of positive *S*_{11} agree with the regions of backscatter in the slices and spanwise average.

The sign of *S*_{11} can be understood using a simple cartoon structure of the squall line (Fig. 18). The main updraft has the maximum vertical velocity at around 3 km (Fig. 3). Therefore, in the neighborhood of the updraft, ∂*w*/∂*z* > 0 below 3 km and ∂*w*/∂*z* < 0 above 3 km. Assuming incompressible or anelastic flow, we must therefore have ∂*u*/∂*x* < 0 below 3 km, where the leftward flow accelerates into the base of the updraft, and ∂*u*/∂*x* > 0 above 3 km (Fig. 14), where the rightward flow accelerates out of the updraft. Therefore, near the updraft, *P*_{11} has dissipation below 3 km and backscatter above 3 km.

Next, we consider the (*i*, *j*) = (1, 3) components. The subfilter energy transfer component *P*_{13} is mainly dissipative (Fig. 16a), which is also shown in horizontal slice (Fig. 17a). The overall shape of the squall line can be seen in the spanwise-averaged *P*_{13} (Fig. 16a). Dissipation is mainly located at the lower and the upper portions of the squall line. There is a small patch of backscatter at low levels to the right of the main updraft, which is also seen in the spanwise-averaged *P* (Fig. 13g).

As for the (*i*, *j*) = (1, 1) component, the sign of *P*_{13}, and therefore the reason that it is dominated by dissipation, can be understood by examining *τ*_{13} and *S*_{13}. The subfilter vertical flux of zonal momentum *τ*_{13} is primarily negative on the right side of the updraft and positive on the left side, as seen in both the horizontal and vertical slices (Figs. 16b, 17b). On the right side of the updraft, subfilter-scale convection transports the leftward momentum upward; on the left side, the rightward momentum is transported upward. The filtered strain rate component *S*_{13} is negative in the turbulent region (Fig. 17b). This turbulent region is the convergence region, where the rear-to-front inflow and the low-level wind shear meet. *S*_{13} is dominated by ∂*u*/∂*z*, which is mainly determined by the flow pattern (Fig. 18). Because of the updraft, a clockwise circulation is created on the right side of the main updraft (Fig. 16c). Hence, ∂*u*/∂*z* and *S*_{13} are positive on the right side of the main updraft. Similarly, the squall line creates a counterclockwise circulation on the left side of the main updraft, leading to negative ∂*u*/∂*z* and *S*_{13}. Putting the signs of *τ*_{13} and *S*_{13} together, *P*_{13} is mainly negative (dissipation), apart from a small region of strong backscatter at the ground, just to the right of the updraft base, where ∂*u*/∂*z* < 0 and *τ*_{13} < 0. Similar to *P*_{11}, while the spanwise-averaged *P*_{13} is dominated by dissipation, it is the net result of lots of backscatter and dissipation in the spanwise direction (Fig. 17a)

## 4. Conclusions

Resolution has great effects on idealized squall-line simulations from Δ*x* = 4 km to 250 m. With TKE mixing, while Δ*x* = 4 km is able to produce the general structures of the squall line, higher-resolution simulations are able to resolve the fine structures of the storms, in agreement with previous studies (e.g., Bryan et al. 2003; Bryan and Morrison 2012). Moreover, higher resolution (Δ*x* = 250 and 500 m) is able to produce squall lines with smaller and more separated convective structures (e.g., Bryan et al. 2003; Bryan and Morrison 2012; Verrelle et al. 2017). Low resolution leads to a longer time to three-dimensionalize and delays the development of the storm, also in agreement with previous studies (e.g., Bryan et al. 2003). The Δ*x* = 1 and 2 km with TKE scheme simulations produce trailing banded clouds, which have not been reported in previous studies. These banded structures do not occur at a lower or higher resolution. The existence of such structures depends on the advection scheme. From the energy spectra, Δ*x* = 4 km is not sufficient to resolve the largest turbulent eddies. Cloud top height increases as resolution increases, which is different from the nonmonotonic trend seen in Bryan and Morrison (2012). This result may due be to the different initial sounding profile since idealized storm simulations are known to be sensitive to environmental conditions (Morrison et al. 2015). The nonmonotonic dependence of the subgrid TKE ratio on resolution is not consistent with past studies such as Bryan et al. (2003), Adlerman and Droegemeier (2002) and Verrelle et al. (2015), where monotonic dependence of subgrid TKE ratio on resolution are found in those studies. From the vertical velocity spectrum, 4 km is also not sufficient to resolve the squall line.

At a fixed resolution, the TKE scheme yields better results than the Smagorinsky scheme in the sense that the lower-resolution TKE results resemble the high-resolution Smagorinsky results. At the coarsest resolution Δ*x* = 4 km, the Smagorinsky scheme simulation does not even produce a propagating squall line, but the TKE scheme does. The Smagorinsky scheme is known to be too dissipative (Pope 2000). In the Δ*x* = 4 km case, excessive mixing from the Smagorinsky scheme is preventing the squall line from properly developing at the earlier stage. Overall, 4 km is not recommended in any mesoscale simulations with explicit convection. Higher resolution yields better simulations. A grid spacing of 250 m is highly recommended for simulations with the TKE and especially the Smagorinsky mixing schemes. Base on the analysis on the subgrid TKE ratio and the vertical velocity spectra, the grid spacing of 250 m produces the best results in both the TKE and Smagorinsky mixing simulations. In general, we recommend using the TKE mixing scheme. Despite these findings, simulations with Δ*x* = 4 km are often used (especially in global simulations) with no convective parameterization (e.g., Done et al. 2004; Weisman et al. 2008).

Significant amounts of backscatter exist across a filter scale of 1 km in the highest resolution simulation. This backscatter cannot be represented by eddy viscosity models such as the Smagorinsky and TKE schemes. Again, these subfilter energy transfers confirm that Δ*x* = 1 km is not ideal for a squall-line simulation because of the backscatter across this scale. We expect less backscatter exists across a filter scale of 500 or 250 m, but this should be confirmed with higher-resolution simulations. The net subfilter energy transfer rate is ahead of the subgrid dissipation with Δ*x* = 1 km, which reflects the delay of three-dimensionalization in the actual 1-km simulation. As expected, dissipation dominates the subfilter energy transfer rate when the squall line is mature. Although the subfilter energy transfer rate is dominated by dissipation, backscatter account for a significant portion of the subfilter energy transfer rate. The subfilter energy transfer rate is not random and is dominated by the streamwise (*i*, *j*) = (1, 1) and shear (*i*, *j*) = (1, 3) components.

Due to the changes in the structure and intensity of the squall line as resolution increases from Δ*x* = 1 km to 250 m, a 1-km filter scale is used in this paper. Further investigations of the subfilter transfer rate across smaller filter scales would be desirable, but require higher resolution. While the subgrid mixing schemes presented in this study do not allow backscatter, models such as the dynamic reconstruction model (DRM) and stochastic subgrid models are alternatives which do allow backscatter in the subgrid turbulence model (Chow et al. 2019). Further work is required to determine whether such models yield improved squall-line simulations.

## Acknowledgments

This research was enabled in part by support provided by the Shared Hierarchical Academic Research Computing Network (SHARCNET: www.sharcnet.ca) and Compute/Calcul Canada (www.computecanada.ca), and funding from the Natural Sciences and Engineering Research Council of Canada (Grant RGPIN-386456-2015) and the Canadian Foundation for Innovation.

## REFERENCES

*Mon. Wea. Rev.*

*Mon. Wea. Rev.*

*Mon. Wea. Rev.*

*Mon. Wea. Rev.*

*Mon. Wea. Rev.*

*Atmosphere*

*Phys. Fluids A Fluid Dyn.*

*Atmos. Sci. Lett.*

*J. Atmos. Sci.*

*Phys. Fluids A Fluid Dyn.*

*Rev. Geophys.*

*Mon. Wea. Rev.*

*J. Turbul.*

*Eur. J. Mech. B/Fluids*

*Mon. Wea. Rev.*

*J. Atmos. Sci.*

*Mon. Wea. Rev.*

*Annu. Rev. Earth Planet. Sci.*

*J. Climate Appl. Meteor.*

*J. Atmos. Sci.*

*J. Atmos. Sci.*

*Mon. Wea. Rev.*

*J. Atmos. Sci.*

*Phys. Fluids A Fluid Dyn.*

*Turbulent Flows*

*Mon. Wea. Rev.*

*J. Atmos. Sci.*

*Meteor. Z.*

*Mon. Wea. Rev.*

*J. Atmos. Sci.*

*Mon. Wea. Rev.*

*Mon. Wea. Rev.*

*Quart. J. Roy. Meteor. Soc.*

*Mon. Wea. Rev.*

*Mon. Wea. Rev.*

*J. Atmos. Sci.*

*Mon. Wea. Rev.*

*Mon. Wea. Rev.*

*Wea. Forecasting*

*Turbulence in the Atmosphere*

## Footnotes

^{1}

Three-dimensionalization is quantified by the ratio of the squared fluctuation of *u* and *w* from their spanwise averages to their squared spanwise averages.

Proc. IBM Scientific Computing Symp. on Environmental Science, Yorktown Heights, NY, Thomas J. Watson Research Center, 195–210