A 25-member ensemble of relatively high-resolution (75-m horizontal grid spacing) numerical simulations of tornadic supercell storms is used to obtain insight on their intrinsic predictability. The storm environments contain large and directionally varying wind shear, particularly in the boundary layer, large convective available potential energy, and a low lifting condensation level. Thus, the environments are extremely favorable for tornadic supercells. Small random temperature perturbations present in the initial conditions trigger turbulence within the boundary layers. The turbulent boundary layers are given 12 h to evolve to a quasi–steady state before storms are initiated via the introduction of a warm bubble. The spatially averaged environments are identical within the ensemble; only the random number seed and/or warm bubble location is varied. All of the simulated storms are long-lived supercells with intense updrafts and strong mesocyclones extending to the lowest model level. Even the storms with the weakest near-surface rotation probably can be regarded as weakly tornadic. However, despite the statistically identical environments, there is considerable divergence in the finescale details of the simulated storms. The intensities of the tornado-like vortices that develop in the simulations range from EF0 to EF3, with large differences in formation time and duration also being exhibited. The simulation differences only can be explained by differences in how the initial warm bubbles and/or storms interact with turbulent boundary layer structures. The results suggest very limited intrinsic predictability with respect to predicting the formation time, duration, and intensity of tornadoes.
Our ability to distinguish tornadic supercell environments from nontornadic supercell environments is much improved relative to decades ago. Today the vast majority of significant (EF2+) tornadoes occur within tornado watches issued by the Storm Prediction Center, and major outbreaks are frequently forecast days in advance. However, once storms form, we cannot say much about specific behaviors, even in supercell storm environments known to be extremely favorable for tornadoes. Even on tornado outbreak days, typically not all storms are tornadic, and those that are tornadic are not tornadic all of the time. The fact that tornado warning lead time has not improved in the past 15 years may be evidence that we are approaching predictability limits (Brooks and Correia 2018).1
The predictability of convective storms has received considerable attention in recent years, especially from those involved with the development of a “Warn-on-Forecast” system (Stensrud et al. 2009; Lawson et al. 2018). The topic of practical predictability probably has been studied most often at this point, that is, the ability to predict storm behavior using the best-available techniques (Lorenz 1969a; Zhang et al. 2006; Melhauser and Zhang 2012). The practical predictability can be limited by uncertainties in both the numerical model (e.g., its numerics and parameterizations) and the initial conditions. Intrinsic predictability, that is, the ability to predict storm behavior using a nearly perfect model and nearly perfect initial conditions, has received less attention in the convective storms community.
Recent and extensive reviews of both the practical and intrinsic predictability of convective storms, and supercell storms in particular, have been provided by Cintineo and Stensrud (2013), Zhang et al. (2015, 2016), Flora et al. (2018), Lawson (2019), and Snook et al. (2019). Supercell predictability studies have included both idealized studies (e.g., Cintineo and Stensrud 2013; Coffer et al. 2017) and case studies (e.g., Zhang et al. 2015, 2016; Flora et al. 2018; Snook et al. 2019). In the idealized studies, environments tend to be horizontally homogeneous2 and storm initiation is typically accomplished via a warm bubble, which eliminates the complexities of convection initiation as a source of error. Practical predictability is assessed by perturbing a control environment, with the perturbation magnitudes being guided by typical forecast or observation errors (Cintineo and Stensrud 2013; Dahl 2014; Coffer et al. 2017). In case study predictability studies, an ensemble of forecasts, often initially produced in real time, is commonly rerun after modifying the initial conditions, physical parameterizations, or resolution of the original ensemble (Flora et al. 2018; Zhang et al. 2015, 2016; Snook et al. 2019). Both the practical and intrinsic predictability have been assessed in such studies; the latter can be investigated by reducing the initial condition perturbations to magnitudes much smaller than can be detected by an observing system (e.g., Zhang et al. 2016). Generally speaking, smaller-scale aspects of supercells (e.g., mesocyclone location, heavy precipitation regions) are less predictable than the larger-scale aspects (e.g., storm location as delineated by the region of light precipitation).
This article is about the intrinsic predictability of tornadoes within supercell thunderstorms. A 25-member ensemble of relatively high-resolution (75-m horizontal grid spacing) numerical simulations of tornadic supercell storms in neutrally stratified, turbulent boundary layers is generated. The turbulence is initiated via small (0.25-K amplitude) random temperature perturbations in the initial conditions. The turbulent boundary layers are given 12 h to evolve to a quasi–steady state before storms are initiated. Storms are initiated at the 12-h mark via the introduction of a warm bubble. The environments have strong vertical wind shear, especially in the boundary layer, a low lifting condensation level (LCL), and large convective available potential energy (CAPE); thus, the environments are extremely favorable for tornadic supercells (Thompson et al. 2003).
The spatially averaged environments and turbulence statistics are identical within the ensemble; only the random number seed and/or warm bubble location is varied from one simulation to the next. Though it cannot be claimed that the numerical model is perfect (the microphysics, turbulence, and surface physics parameterizations are probably the most significant sources of error), we can treat it as perfect and consider only the effects of different boundary layer realizations or warm-bubble placements on the outcomes in order to assess the intrinsic predictability of tornadic supercell storms.3
The prior predictability studies most applicable to this study are probably those by Dahl (2014), Coffer et al. (2017), Markowski and Richardson (2014a, 2017), Yokota et al. (2018), and Snook et al. (2019), all of which might best be regarded as studies of practical predictability. Dahl (2014) analyzed the formation and characteristics of tornado-like vortices in a 61-member ensemble of simulations with 100-m horizontal grid spacing. The simulations were initialized with a proximity sounding obtained near the 29 May 2004 Geary, Oklahoma, tornadic supercell. The environment was horizontally homogeneous, but small, random errors were added to the vertical profiles of temperature, moisture, and wind following the Cintineo and Stensrud (2013) approach, in that the amplitude of the perturbations was guided by typical model errors. The formation of intense vortices, and even the width and motion of the vortices, was sensitive to the initial condition perturbations. In general, the larger the perturbations, the larger the spread in solutions, though there was some indication of the existence of a threshold “beyond which (the) reduction of error in the initial conditions is unlikely to greatly improve the forecast.”
Coffer et al. (2017) added 2 m s−1 random perturbations to the tornadic and nontornadic supercell wind profiles derived from soundings launched during the Verification of the Origins of Rotation in Tornadoes Experiment 2 (VORTEX2; Wurman et al. 2012; Parker 2014; Coffer and Parker 2017). The thermodynamic profiles were not perturbed. The wind perturbations were intended to represent the effects of turbulent eddies on the wind profile, though the simulations were performed in laminar, horizontally homogeneous environments. The horizontal grid spacing was 125 m. Intense, long-lived, tornado-like vortices developed in each of the 15 tornadic-environment simulations; 40% of the supercells simulated in the perturbed nontornadic environments (15 additional simulations) were described as “weakly tornadic.” Coffer et al. concluded that “chaotic, within-storm details can still play a role and, occasionally, lead to marginally tornadic vortices in suboptimal storms.”
Markowski and Richardson (2014a, 2017) simulated supercell-like, dry, pseudostorms (100-m horizontal grid spacing) and investigated the sensitivity of tornado-like vortex formation to the strength and location of a heat sink. The heat sink emulated the latent chilling that occurs within real storms and was crucial for the development of low-level rotation. Varying the heat sink strength and location is a controlled way of exploring the sensitivity of vortex formation to the downdraft position and strength, which would depend on hydrometeor species and deep-tropospheric wind shear, among other things, in an actual storm, and on the microphysics parameterization in a simulation that includes moist processes. Vortex formation and intensity were found to be extremely sensitive to a horizontal displacement of the heat sink of only a couple of kilometers—distances comparable to a single grid length in today’s operational convection-allowing models. It was concluded that the volatility associated with downdraft position and strength, which affect the baroclinic vorticity generation within the storm, may explain the failure of many supercells to produce tornadoes in seemingly favorable environments.
Yokota et al. (2018) investigated the dynamics of tornadogenesis (50-m horizontal grid spacing) in a 33-member ensemble of supercell simulations for a case that occurred in a landfalling tropical cyclone in Japan. The initial conditions and perturbations came from the Japan Meteorological Agency’s operational mesoscale analysis. The observed storm produced an EF3 tornado, but tornadoes formed in only seven of the 33 ensemble members. The origins of the tornado’s vorticity and mechanisms of vorticity amplification were examined in the seven tornadic cases. Both baroclinic and frictional vorticity generation were found to be important, but curiously the degree of importance varied from simulation to simulation, and tornado intensity was unrelated to the vorticity generation mechanism. What seemed to matter most was the intensity of the dynamically driven low-level updraft, similar to the finding of Coffer and Parker (2017).
Snook et al. (2019) used an ensemble of 10 forecasts of the 20 May 2013 Newcastle–Moore EF5 tornadic supercell (this is the same case studied by Zhang et al. 2015, 2016) to examine the practical predictability of tornadogenesis and tornado characteristics. The horizontal grid spacing was 50 m within the 100 km × 71.5 km inner nest. The ensemble members used different boundary layer parameterizations on the coarse grids (no boundary layer schemes were used on the 50-m grids), as well as perturbed initial and boundary conditions, with the magnitude of the perturbations being governed by typical error magnitudes. Though the supercells were tornadic in all ten members, tornado intensity ranged from EF0 to EF5, and the time of tornadogenesis varied by 80 min, implying limited practical predictability of tornado genesis and characteristics.
The present study differs from prior studies in that the only source of ensemble spread is different realizations of the turbulent boundary layer. (Warm bubble location also is varied, but this is really just a less expensive way to increase the diversity of interactions between storms and boundary layer turbulence without spinning-up additional boundary layers.) This is believed to be the first intrinsic predictability study of tornado genesis and tornado characteristics. Moreover, given the grid spacing, at least in the boundary layer, the simulations can be viewed as so-called large-eddy simulations (LES). This might be the first supercell predictability study using “true” LES (i.e., simulations in which a large fraction of the boundary layer turbulence is explicitly resolved). Though supercell simulations have routinely used LES turbulence schemes since the earliest days (Klemp and Wilhelmson 1978), such schemes are questionable when turbulence is not resolved, either because of insufficient resolution, or because of a lack of turbulence-triggering perturbations or “eddy injection” through open, inflow boundaries. It is likely that resolved turbulence is present in the environment on the finest (50-m) grids in the Yokota et al. (2018) and Snook et al. (2019) studies, but it is unclear whether it would be fully developed by the time inflow air reaches the storms. It would likely have only 20–60 min to develop, given the inflow speeds and proximity of the storms to the boundary of the finest grid.
One aspect of this study that is worth emphasizing is its idealized nature. The storm environments have no mean horizontal gradients, and storm initiation, accomplished via a warm bubble, is extremely idealized. With respect to the latter, the intrinsic predictability of the storms themselves can be examined precisely because the storm initiation is so idealized. In nature, forecast errors are greatly affected, perhaps even dominated in some situations, by what happens during and shortly after convection initiation, such as interactions between horizontal convective rolls and mesoscale boundaries (e.g., Atkins et al. 1995; Xue and Martin 2006), and the formation of multiple updrafts and precipitation cores and their subsequent interactions and mergers (e.g., Hastings et al. 2010; Skinner et al. 2014; Hastings and Richardson 2016; Klees et al. 2016). Given that there is complete control of the timing and location of storm initiation in this study, coupled with the fact that the mature storms simulated in this study are isolated and in environments extremely favorable for tornadoes—as opposed to environments only marginally supportive of tornadoes or even supercell storms—the predictability assessed in this study should probably be regarded as a best-case scenario.
The numerical simulations were performed using Cloud Model version 1 (CM1; Bryan and Fritsch 2002), release 18.3. The relevant model parameters are listed in Table 1; only the most important aspects of the simulations are explained below.
The domain is 127.5 km × 127.5 km × 18.0 km (1700 × 1700 × 121 grid points). The horizontal grid spacing is 75 m throughout the domain. The vertical grid spacing varies from 15 m at the surface to 285 m at the top of the domain. The lateral boundaries are periodic. The top boundary is rigid and free slip. A semislip boundary condition is applied at the bottom of the domain, with the roughness length set to 12 cm, which corresponds to a nondimensional drag coefficient of 0.0094. The vertically implicit time-splitting method of Klemp and Wilhelmson (1978) is used in conjunction with adaptive large and small time steps. Throughout most of the simulations, the large and small time steps are 1.0 and 0.125 s, respectively.
Subgrid-scale turbulence is parameterized using a turbulent kinetic energy (TKE) scheme similar to the one used in Deardorff’s (1980) LES. The National Severe Storms Laboratory (NSSL) double-moment microphysics scheme (Ziegler 1985; Mansell et al. 2010; Mansell and Ziegler 2013) is used (see Table 1 for additional details). The Coriolis acceleration is included (f plane assumed, with f0 = 10−4), but only acts on horizontal velocity perturbations relative to the initial, base state. Surface heat and moisture fluxes are excluded, as is radiative transfer. Though several studies have shown that radiative transfer processes influence supercell storms (Markowski et al. 1998; Markowski and Harrington 2005; Frame and Markowski 2010, 2013; Nowotarski and Markowski 2016), it is unclear whether they might influence the intrinsic predictability of the simulated storms.
The initial vertical profiles of temperature and moisture are similar to those used by the Weisman and Klemp (1982; hereafter WK82), though some modifications were required via trial and error, given the evolution of the sounding during the boundary layer spinup period (Fig. 1a). Without the modifications, both a moist absolutely unstable layer and unwanted/uncontrolled convection initiation occur as a result of the sounding evolution during the spinup period. The initial boundary layer water vapor mixing ratio, the tropopause potential temperature, and the tropopause temperature were set to 0.015, 333, and 203 K, respectively (WK82’s variables qυ0, θtr, and Ttr). Moreover, the exponents in WK82’s analytic functions for potential temperature and relative humidity [see WK82’s Eqs. (1) and (2)] were changed to 1.0 and 0.5, respectively. The initial vertical wind profile was specified using the analytic function given in WK82’s Eq. (4), with Us = 30 m s−1 and zs = 8 km, but with a background velocity of (5, 20) m s−1 superimposed. The initial wind profile is characterized by westerly shear and veering of winds with height, with the shear being strongest at low altitudes and gradually weakening with height (Fig. 1b). Even though there is no base state horizontal pressure gradient, because the Coriolis acceleration only acts on horizontal velocity perturbations relative to the initial wind profile, the wind profile can be considered to represent the base state geostrophic wind profile. Geostrophic wind hodographs resembling the base state hodograph in Fig. 1b are common in the warm sectors of extratropical cyclones (Banacos and Bluestein 2004).
Random temperature perturbations having a maximum amplitude of 0.25 K are imposed in the lowest 1 km at t = 0 in order to excite turbulence. The grid spacing is sufficiently small such that turbulent eddies are explicitly resolved. The boundary layer evolves throughout a spinup period of 12 h, by which time it attains an approximately steady state. The spinup period can be visualized in an animation that is available in the online supplemental material.
Five different boundary layers are developed using five different random number seeds (Fig. 2). The boundary layers have identical horizontal mean fields. The effect of the turbulent motions leads to well-mixed boundary layers in which the horizontally averaged potential temperature and water vapor mixing ratio are constant with height (Figs. 3a,b). The quasi-steady-state boundary layer characteristics are independent of the magnitude of the initial perturbations. Rather, they are determined by the base state vertical profiles of wind and temperature. For example, the boundary layer at t = 12 h obtained using much smaller initial temperature perturbations, having a maximum amplitude of only 0.05 K, is virtually indistinguishable from the boundary layer obtained with the 0.25-K-amplitude initial temperature perturbations. This is the case not only for the vertical profiles of potential temperature and water vapor mixing ratio (Figs. 3a,b), but also for higher-order statistics (e.g., velocity and temperature variances) (Figs. 3c,d).
The combined effects of surface drag, the Coriolis acceleration, and base state geostrophic wind profile lead to a quasi-steady-state hodograph that has significantly more shear than the base state geostrophic wind profile, in addition to substantial curvature in the lowest kilometer. The environment is extremely favorable for tornadoes (Fig. 1), with a surface-based CAPE of 3683 J kg−1, 0–1-km (0–3-km) storm-relative helicity of 227 (300) m2 s−2 given the ensemble-mean storm motion, an LCL of 1.1 km, and a significant tornado parameter4 of 6.1. Because there is no surface heat flux, the quasi-steady boundary layer might best be regarded as a late-day boundary layer near the time of the early evening transition. This is the time of day when tornadoes are most likely in the U.S. Great Plains region anyway (Anderson-Frey et al. 2017; see their Fig. 2a).
At t = 12 h, deep convection is initiated by placing an ellipsoidal, warm, humid bubble within the quasi-steady, turbulent boundary layer, similar to Nowotarski et al. (2015). The bubble has a maximum potential temperature perturbation of 4 K, a relative humidity of 95%, a horizontal radius of 10 km, and a vertical radius of 1.4 km. The bubble is centered at z = 1.4 km. The bubble is inserted at five different horizontal positions within the ensemble. The default is the center of the domain, with four additional positions being 2 km north, south, west, and east of the default position.5 The combination of five different random number seeds (and boundary layers) and five different warm-bubble-insertion locations yields a 25-member ensemble of storm simulations. The warm bubbles quickly trigger supercell storms. Storm-splitting occurs in the first 30 min after initiation, after which an intense, cyclonically rotating supercell dominates.6 The remainder of the paper deals with the behavior and predictability of this storm. The simulations are stopped at t = 14 h.
The numerical simulations were performed on the Institute for Computational and Data Sciences Advanced CyberInfrastructure (ICDS-ACI) system at Penn State. Simulations usually were run using 100 cores across five nodes. The five 12-h boundary layer spinup simulations required approximately 250 000 core hours, and the 25 2-h storm simulations required approximately 400 000 core hours.
a. Overview of storm evolution and variability within the ensemble of simulations
Each simulation contains an intense, cyclonically rotating, supercell storm that persists throughout the entire 2 h of simulation time following the insertion of the warm bubble. Hereafter, references to time (t = 30 min, t = 60 min, etc.) refer to the time elapsed since the insertion of the warm bubble. Select ensemble mean fields are shown in Fig. 4. Prior to averaging, the fields are shifted horizontally in order to account for differences in the location of warm bubble insertion and storm motion. This is accomplished by shifting the fields of the nth ensemble member in order to maximize the correlation between the vertical velocity fields in the 3–10-km layer in the nth ensemble member and first ensemble member. Even the characteristics of the mean fields are strongly suggestive of a supercell storm (e.g., hook echo, rotating updraft; cf. Fig. 7 of Lemon and Doswell 1979). The maximum vertical velocities in each simulation are approximately 80 m s−1 (not shown). The forward-flank precipitation region has considerably less variance than the rear-flank region, both in terms of horizontal velocity (Fig. 5) and potential temperature (Fig. 6), which is not surprising given that the latter region has been long-recognized as being more turbulent (Brandes et al. 1988; also see Fig. 1 in Markowski and Bryan 2016). Animations are included in the online supplement.
The largest velocity variance within the ensemble is found in the low-level mesocyclone region (Fig. 5), and is a consequence of significant differences in the intensity of the tornado-like vortices that form in the simulations. Hereafter, vortices are simply referred to as tornadoes if their vertical vorticity exceeds 0.25 s−1 and the ground-relative wind speed exceeds 29 m s−1 at the lowest grid level (z = 7.5 m), even though “tornado-like vortices” might be more appropriate given the relatively coarse grid spacing, at least relative to grid spacings typically used in tornado-resolving simulations. The differences in the timing and locations of tornado formation are evident in Figs. 7–9 . Although tornadoes develop in all 25 simulations, the peak intensity of the tornadoes ranges from EF0 (29–38 m s−1) to EF3 (60–73 m s−1) on the enhanced Fujita scale, with the peak in the distribution lying within the EF1 range (38–49 m s−1) (Fig. 9). In most of the simulations (24/25), at least a brief, weak tornado develops within 60 min of warm bubble insertion. The early tornado episodes are generally followed by a 20–30-min lull, after which a second, more significant, tornadic phase occurs (Figs. 8, 9). By t = 2 h, the tornado swaths span a region 15-km wide (Fig. 8). The swaths in Fig. 8 are shifted relative to the placement of the initial warm bubbles (i.e., the differences in convection initiation location have been eliminated as a source of tornado-location variability). It is also worth noting that the differences in storm behavior or tornado intensity within the ensemble do not depend on whether an ensemble member was created by perturbing the environment (via a different random number seed) or shifting the location of the warm bubble that initiates each storm.
The ensemble mean storm motion components are (12.2, 12.5) m s−1 from t = 1 h to t = 2 h, with a standard deviation of (0.6, 0.5) m s−1. Instantaneous storm and tornado translational velocities, however, vary by as much as 5 m s−1 within the ensemble (the directional differences in tornado motion are evident in Fig. 8), even though the mean wind profiles are identical in the storm environments. Though the variable storm motions imply variable SRH within the ensemble, the effect on SRH is small (0–1-km SRH varies by less than 10 m2 s−2 as a result of storm-motion variability).
b. Differences between the strongly tornadic and weakly tornadic storms
The nine simulations in which significant tornadoes (EF2+) develop (members 1, 2, 3, 7, 9, 11, 19, 20, 24) are compared to the nine simulations having the weakest tornadoes (members 4, 6, 10, 13, 14, 16, 17, 18, 22). Hereafter, the respective storm types are referred to as “SIGTOR” and “WEAKTOR.” Even though weak tornadoes develop in 16 of the simulated storms, when comparing mean SIGTOR storm characteristics to mean WEAKTOR storm characteristics, it was deemed desirable to average over equal numbers of storms so that the mean fields of both storm types contain similar spatial scales.
Given the design of the numerical experiments, all of the aforementioned differences in storm behavior and tornado intensity are ultimately due to differences in how the initial warm bubbles and/or resulting storms interact with the turbulent boundary layers in their respective storm environments. The turbulent boundary layer structures in the vicinity of the warm bubbles used to initiate what would eventually become SIGTOR and WEAKTOR storms are indeed different (Fig. 10). For example, in the SIGTOR simulations the warm bubbles are introduced in locations where boundary layer vertical velocities tend to be positive on the east side of the bubble and negative on the west side of the bubble (Fig. 10a), and in the WEAKTOR simulations (Fig. 10b) the warm bubbles tend to be centered near local maxima in boundary layer vertical velocity. However, these are not the only differences evident in Fig. 10, and it is impossible to say whether these or any other differences are relevant. There is no obvious reason from comparing Figs. 10a and 10b to expect different storm behaviors as a result of the differences in the kinematic fields where the warm bubbles are introduced.
Storm characteristics are investigated next. In particular, the focus is on the attributes of supercell storms that have been identified in prior studies, both observational and numerical, as being able to discriminate between supercells that spawn significant tornadoes and those that are weakly tornadic or nontornadic. These key storm attributes include the negative buoyancy of the storm outflow (Markowski et al. 2002; Grzych et al. 2007; Hirth et al. 2008; Snook and Xue 2008; Markowski and Richardson 2014a), the strength of the low-level updraft (Markowski and Richardson 2014a; Coffer et al. 2017; Yokota et al. 2018), and the proximity of low-level angular momentum to the dynamically driven low-level updraft (Snook and Xue 2008; Markowski and Richardson 2014a). Negative buoyancy in the outflow implies the presence of a horizontal buoyancy gradient and baroclinic horizontal vorticity generation, but the degree to which this horizontal vorticity can be tilted upward and subsequently intensified via stretching (i.e., the degree to which angular momentum can develop about a vertical axis and subsequently be acted upon by horizontal convergence) also depends on the buoyancy of the air. The horizontal vorticity generated baroclinically tends to increase as the negative buoyancy of the outflow increases. However, too much negative buoyancy can impede both the tilting and stretching of the baroclinically generated horizontal vorticity, and also can cause whatever angular momentum (about a vertical axis) that does develop to be displaced too far from the region of strong updraft forcing, thereby limiting the convergence of angular momentum (Snook and Xue 2008; Markowski and Richardson 2014a). Surface friction can be an additional vorticity source for the tornado (Schenkman et al. 2014, 2016; Markowski 2016; Roberts et al. 2016). There is evidence that its relative importance is greatest when cold pools are especially weak or absent, such as early in the life of a storm (Markowski 2016; Mashiko 2016; Roberts et al. 2016).
Despite the complexities involving internal storm processes summarized above, the larger-scale environment has been found to exert considerable influence on supercell behavior. Specifically, environments with especially strong low-level wind shear (even by supercell storm standards, which have stronger shear than the environments associated with ordinary storms; e.g., WK82) and a low LCL have been found to be conducive to supercells that produce significant tornadoes (e.g., Thompson et al. 2003). Strong low-level wind shear, especially if it is associated with a large streamwise horizontal vorticity component, fosters a stronger dynamically driven low-level updraft (Markowski and Richardson 2014a; Coffer et al. 2017). A low LCL reduces the likelihood of outflow being excessively cold (evaporation of rain is usually the dominant source of negative buoyancy) and has been found to be favorable for significant tornadoes (Markowski et al. 2002; Shabbott and Markowski 2006). It is not yet known whether there are certain environments in which surface friction might be more likely to enhance tornado formation.
In the simulations herein, the environments are statistically identical and extremely favorable for significant tornadoes. But can the differences between the SIGTOR and WEAKTOR storms still be explained in terms of systematic differences in the negative buoyancy of the outflow, the strength of the low-level updrafts, or the amount of angular momentum in proximity to the low-level updrafts? (Again, any differences in the negative buoyancy of the outflow, low-level updraft strength, etc., ultimately would be attributable to the different realizations of a turbulent boundary layer and/or the location of warm-bubble insertion.)
The answer appears to be no. Although there are occasionally differences between the SIGTOR and WEAKTOR storms in their time series of vertical vorticity, cold pool buoyancy, circulation, vertical velocity, and vertical perturbation pressure gradient acceleration (VPPGA), the differences are generally small, not statistically significant,7 and unlikely to be observable (Fig. 11). The cold pool buoyancy is assessed in terms of density potential temperature perturbations (Emanuel 1994, p. 113), where the density potential temperature is θρ ≈ θ(1 + 0.61qυ − qh), qυ is the water vapor mixing ratio, qh is the hydrometeor mixing ratio, and the perturbation density potential temperature, , is related to the buoyancy B via , where g is the gravitational acceleration and is the base state density potential temperature. The circulation is , where is the tangential velocity averaged about a circle of radius r = 1 km surrounding each grid point. The VPPGA is −ρ−1∂p′/∂z, where ρ is the density and p′ is the pressure perturbation relative to the initial hydrostatic, horizontally homogeneous reference state; the VPPGA is the main driver of the low-level updraft and is crucial to tornado formation given the need to accelerate negatively buoyant air upward.
Interpretation of the time series in Fig. 11 is complicated by the fact that they are strongly influenced by tornadoes themselves, particularly the SIGTOR time series. At the rare times when the SIGTOR storms have stronger low-level updrafts, near-surface vertical vorticity, etc. (e.g., 105–120 min), five or more of the SIGTOR storms have significant tornadoes in progress, which are associated with extreme perturbations in the kinematic and pressure fields (cf. Figs. 9, 11). Even circulation, which is proportional to the angular momentum and area-averaged vertical vorticity, is typically greatly enhanced by tornadoes; although C would be conserved (in the inviscid, nonrotating limit) about a material circuit confined to a horizontal plane that is converging upon a vortex (i.e., dC/dt = 0), vertical vorticity stretching increases the C about a fixed-radius circuit (i.e., ∂C/∂t > 0).
Figures 12 and 13 depict cold pools and and near-surface circulation (the low-level updraft is overlaid on both) in each of the SIGTOR storms 5 min prior to their first significant tornadogenesis, and Figs. 14 and 15 show the same fields in each of WEAKTOR storms. It is far from clear when the WEAKTOR storms should most fairly be evaluated and compared to the SIGTOR storms. The WEAKTOR storms are presented at t = 95 min, which is the average time in the simulations when the aforementioned SIGTOR storms have been evaluated. Significant variability characterizes the fields, with no obvious, systematic, differences, at least in the eyes of the author. SIGTOR ensemble member 11 (Figs. 12f, 13f) is particularly unremarkable in its low-level thermodynamic and vertical velocity fields, yet 5 min later, an EF2 tornado develops with peak winds reaching 55.8 m s−1. The appearances of WEAKTOR ensemble members 16 and 22 (Figs. 14f, 14i) might be considered impressive by most readers, given the coiling of the hook echoes and strong low-level updrafts. The near-surface C in WEAKTOR ensemble members 4 and 14 (Figs. 15a,e) is larger than in any of the SIGTOR storms, and in general, high-C air does not appear to be systematically better-placed relative to the updraft maxima in the SIGTOR storms than in the WEAKTOR storms (cf. Figs. 13, 15). The average distances between the locations of wmax and Cmax in the SIGTOR and WEAKTOR storms are 1.9 and 1.8 km, respectively, though the difference is not statistically significant and is sensitive to the exact altitudes at which the maxima are compared.
Comparisons of mean SIGTOR and WEAKTOR fields also fail to yield obvious reasons why the SIGTOR storms would imminently produce significant tornadoes while the WEAKTOR storms would not (Figs. 16–18 ). The maximum updraft at z = 1 km is actually slightly stronger in the WEAKTOR mean (14 m s−1) than in the SIGTOR mean (12 m s−1) (Fig. 16), though the differences between the mean vertical velocity fields are not statistically significant anywhere in the updraft region (or downdraft region, for that matter). The SIGTOR mean cold pool is slightly cooler than the WEAKTOR mean (minimum versus −2.8 K) (Fig. 16). The cold pool differences are not statistically significant either. The gust front is bit more “inflected” in the SIGTOR mean than in the WEAKTOR mean, but it is difficult to know what the implications might be.
Clues as to why the SIGTOR storms produce more intense vortices also are absent in the near-surface mean circulation fields. The circulation in the SIGTOR mean velocity field (maximum C ≈ 10 500 m2 s−1) is actually slightly less than in the WEAKTOR mean velocity field (maximum C ≈ 13 000 m2 s−1), and the high-C air is not clearly better placed relative to the updraft in the SIGTOR mean than in the WEAKTOR mean (Fig. 17). In the SIGTOR mean, a curious “tongue” of high-C air extends from the mesocyclone region into the environment immediately ahead of the gust front (dashed line in Fig. 17a), but it would be difficult to determine whether or not this feature is dynamically consequential. None of the differences in the mean circulation fields are statistically significant.
The region of enhanced streamwise vorticity feeding the low-level updraft from the east or northeast (Klemp and Rotunno 1983; Rotunno and Klemp 1985), on the cool side of the gust front, is another parameter that has received renewed interest lately (e.g., Orf et al. 2017). The low-level streamwise vorticity in the SIGTOR and WEAKTOR mean velocity fields lacks obvious (or statistically significant) differences that might be dynamically important (Fig. 18).
The title of this article asks “What is the intrinsic predictability of supercells?” Based on section 3, the best answer is “it depends.” On one hand, storm mode, intensity, and track are extremely predictable. All of the storms are long-lived supercells with intense updrafts and mesocyclones, have an approximately steady motion toward the northeast, and all probably would warrant the issuance of a tornado warning. All of the storms are at least weakly tornadic, and the envelope of tornado swaths spans just 15 km in width as far out as 2 h from storm initiation. Thus, tornado warnings reliably could be issued at least 2 h in advance for a targeted area smaller than the area of a typical U.S. county.
However, predictability is scale-dependent (e.g., Lorenz 1969b), with the smallest resolved scales of motion—including those associated with the tornadoes—having very short predictability horizons. Tornado genesis (its specific timing, that is), intensity, and longevity within a particular storm have very limited predictability. It is entirely possible, perhaps even likely, that the tornado intensity distribution would be even broader at a smaller grid spacing than used herein, given the likelihood that the strongest tornadoes would be even stronger using finer grid spacings (as explained in section 2, the 75-m horizontal grid spacing is marginal for resolving tornadoes). Although there is uncertainty regarding what might happen to the intensity of the weakest tornadoes at finer grid spacings, there is no reason to expect that the distribution of tornado intensity would narrow as grid spacing decreases.
The differences in outflow buoyancy, low-level updraft strength, and circulation highlighted in section 3b fail to point toward an obvious reason for the more intense vortices in the SIGTOR storms. If anything, based on the ensemble means and what we know from prior studies about the importance of the outflow buoyancy, low-level updraft strength, and circulation in tornadogenesis, we might have been justified in predicting that the WEAKTOR storms would have produced stronger vortices than the SIGTOR storms. Both sets of storms, however, have attributes favorable for tornadoes. For example, the cold pools in both the SIGTOR and WEAKTOR simulated storms are weak relative to the cold pools in observed nontornadic storms (the latter typically have density potential temperature deficits exceeding 4 K; e.g., Markowski et al. 2002; Shabbott and Markowski 2006). Moreover, both the SIGTOR and WEAKTOR storms have plenty of near-surface circulation and sufficiently strong low-level updrafts. The results demonstrate that the attributes of a storm that we commonly regard to be favorable for significant tornadoes—as well as the environmental parameters used to anticipate favorable storm attributes (e.g., strong low-level shear, low LCL)—have a limited influence on the evolution of tornadoes in a particular storm. In other words, even in environments favorable for significant tornadoes, not all supercells produce significant tornadoes. This is not necessarily surprising, and is in basic agreement with the prior findings of Dahl (2014), Markowski and Richardson (2017), Yokota et al. (2018), and Snook et al. (2019). What is perhaps surprising, however, is that the differing behaviors are not so straightforward to explain in terms of our current understanding of tornadogenesis [see section 3b and Markowski and Richardson 2014b], even with the benefit of four-dimensional fields of model output. The “flapping of butterfly wings” might seem to be as good an explanation as any for the observed spectrum of storm behaviors. Perhaps this is a problem well-suited for “deep learning” (LeCun et al. 2015).
Of course, all differences in vortex intensity are due to differences in the degree to which angular momentum can be converged to a small radius, with the angular momentum contraction being accomplished by low-level upward accelerations driven by the VPPGA and retarded by negative buoyancy. In other words, the lack of obvious differences between the SIGTOR and WEAKTOR storms does not mean that there is nothing dynamically different about the SIGTOR storms at substorm scales that favors significant tornadogenesis. The issue here is that differences in upward accelerations, angular momentum, and their superpositioning that lead to more intense vortices in the SIGTOR storms are apparently masked by ensemble averages and/or only arise immediately before or after tornado formation.
It is worth restating that the ensemble spread is solely attributable to how the warm bubbles and/or storms interact with turbulent eddies in the boundary layer, given that the environments are statistically identical and the initial conditions differ only by a random number seed. Moreover, the ensemble spread might be about as small as can be obtained in a suite of convective storm simulations, both because the environment is so favorable for tornadic supercells and convection initiation is idealized. Supercells are arguably the most long-lived and predictable storm type (a possible exception might be a long-lived mesoscale convective system). The supercell storms simulated herein are nearly steady throughout the simulations and would probably last forever if the simulations ran as long. Lilly (1986) hypothesized that the high helicity of supercells (i.e., their rotating updrafts) might enhance their predictability, though Peters et al. (2020) have shown that this helicity effect is not as important as suggested by Lilly. Droegemeier and Levit (1993) found that supercell storms were more predictable than multicell storms, however.
One might wonder why such an extreme part of the convective storm parameter space was explored in this study, as opposed to one marginally favorable for tornadoes, or even supercells. It was decided to first look at an extremely favorable environment because predictability is best there, or at least tornado warning skill is (Anderson-Frey et al. 2016; Anderson-Frey and Brooks 2019). For environments that are more “borderline” between supercells and nonsupercells, or near a “tipping point” between tornadic and nontornadic supercell environments (Coffer and Parker 2018), it is virtually certain that the differences between storm outcomes would be larger. Furthermore, a more realistic convection initiation involving interactions between boundary layer thermals and mesoscale boundaries would be an additional source of variability among the outcomes. In case studies investigating the practical predictability of convective storm forecasts, convection initiation errors are frequently the greatest source of uncertainty (e.g., Crook 1996; Brooks et al. 1993; Martin and Xue 2006; Zhang et al. 2016; Flora et al. 2018).
5. Summary and conclusions
Tornadic supercell storms were simulated in environments that only differed by the random number seed used to initialize the simulations with small-amplitude (0.25 K) temperature perturbations. The initial temperature perturbations were needed to develop realistic turbulent boundary layers, after which the storms were initiated by introducing warm bubbles into the boundary layers. The warm bubble location also was varied from simulation to simulation.
Long-lived, intense supercell storms, and at least weak tornadoes, develop in each of the 25 simulations. However, despite the storm environments being identical in terms of the vertical profiles of mean temperature, humidity, and wind, the timing, duration, and intensity of the tornadoes that develop in any particular simulation have limited predictability. The simulation differences only can be explained by differences in how the initial warm bubbles and/or storms interact with turbulent boundary layer structures.
The results probably represent a “best-case scenario” with respect to storm predictability, because convection initiation was highly idealized and the environment was so favorable for supercells and tornadoes. The range of storm behaviors would likely be larger in environments only marginally supportive of supercells and/or tornadoes, or for less steady modes of moist convection (e.g., multicell storms). More realistic convection initiation, which routinely involves complex interactions between boundary layer thermals and mesoscale airmass boundaries in the real atmosphere, would be an additional source of variability within an ensemble. Though it seems unlikely, this study could not explore whether any physical parameterization in the simulations (e.g., microphysics), or a parameterization omitted from the simulations (e.g., radiation), makes simulated storms less predictable than storms naturally would be.
The results suggest very limited intrinsic predictability with respect to tornado forecasting within supercell storms. Here “tornado forecasting” refers to the prediction of the time of tornadogenesis and tornado duration and intensity. Issuing potentially life-saving tornado warnings for the occurrence of tornadoes within a polygon is a different topic. In the case of the particular long-lived supercell storm simulated herein, warnings could have been successfully issued with lead times for many tens of minutes, if not an hour or longer.
The results also imply that a single, high-resolution simulation is of limited use for investigating environmental controls on tornadogenesis or tornado characteristics (and this article does not even touch upon the sensitivities to the parameterization of microphysics, radiation, turbulence, and surface fluxes). Potvin et al. (2017) have similarly demonstrated the importance of using an ensemble framework. They investigated the short-term forecast impacts of unanalyzed initial-condition scales by performing supercell simulations with differing initial-condition resolutions.
The present study does not exclude the possibility that tornado statistics (e.g., the distributions of tornado intensity and longevity) are much more predictable. More studies like this are warranted to investigate whether such statistics are predictable as a function of the characteristics of the environment. Though such a study is probably unfeasible with present day computing power, it ought to be within reach in the next few decades as tornado-resolving “full-storm” simulations (as opposed to “stripped-down” tornado simulations, either in limited domains without a parent storm, or in larger domains with a highly idealized parent storm) become much more commonplace, owing to increases in both computing power and storage.
One of the anonymous reviewers commented that predictability ought to be higher for a storm that has already developed than for a storm that has not yet developed. For example, the predictability of a storm’s future behavior is much better in short-range, convection-resolving numerical forecasts once the storm has developed and has been assimilated, than is the predictability before the storm has developed. Although the complexities of convection initiation as a source of error were excluded from this study by design, the author could devise no straightforward way to spin up identical mature storms before releasing them into different realizations of turbulent boundary layers. The introduction of small perturbations into “restart” simulations of the storms might be explored in a future study.
It is perhaps fitting to end this article with the following quote by Scorer (1978, p. 248): “This coming to terms with our limitations is actually much more satisfactory than it may seem because it amounts to making our objectives sensible. It is obvious that the deeper we probe into something the more complex it seems to be, and in the case of indescribably complicated motion, the only road to simplification is to decide that there shall be limits to the complexity we are prepared to study. Life is too short to spend it sorting out the infinite details of the flow on a single occasion, and anyone who did analyze them could have no assurance that the next occasion would be the same, nor that anyone else would be interested anyway.”
This article is dedicated to the memory of Prof. Fuqing Zhang, with whom the author discussed the subject of tornado predictability over the years. The author thanks Steve Greybush, Zach Lebo, Dave Stensrud, and Yunji Zhang for their feedback on the direction of this research. Craig Schwartz and Ryan Sobash also are acknowledged for their assistance in interpreting smoothed tornado probabilities. Three anonynous reviewers helped improve the article, and for this the author is grateful. The work was partially supported by National Science Foundation Award AGS-1821885. The numerical simulations were performed on the Institute for Computational and Data Sciences Advanced CyberInfrastructure (ICDS-ACI) system at Penn State. A large fraction of the figures was created using the Grid Analysis and Display System (GrADS), developed by the Center for Ocean–Land–Atmosphere Studies. Figure 1 was created using GrADS software written by Bob Hart.
Here, as in most publications in the severe storms community, “horizontally homogeneous” refers to a laminar storm environment that has no gridpoint-to-gridpoint horizontal variability. In contrast, in the boundary layer community, horizontally homogeneity is defined not by gridpoint-to-gridpoint variability (which is always present when turbulence is present), but by whether averaged quantities vary horizontally (spatial, temporal, or ensemble averages).
It is implicitly assumed that the particular choices of parameterizations do not artificially enhance or limit ensemble spread relative to other parameterization choices.
The displacements are really 2025 m relative to the default location so that the bubble is centered on a scalar grid point within a horizontal plane.
The anticyclonically rotating “left mover” that emerges from the storm splitting process tracks northward and would eventually influence the storm of interest from the south owing to the periodic lateral boundaries. However, the simulations are stopped before this interaction happens.
“Statistical significance” throughout the paper refers to p values less than 0.05.