Brief communication: Towards a universal formula for the probability of tornadoes

A methodological approach is proposed to provide an analytical (exponential-like) expression for the probability of occurrence of tornadoes as a function of the convective available potential energy and the wind shear (or, alternatively, the storm relative helicity). The resulting expression allows the probability of tornado occurrence to be calculated using variables that are computed by weather prediction and climate models, thus compensating for the lack of resolution needed to resolve these phenomena in numerical simulations.

1 Introduction

Tornadoes are rapidly rotating columns of air (American Meteorological Society, 2020), extending vertically from the surface to the base of a cumuliform cloud, and represent one of the most severe weather phenomena in terms of victims and damage. Considering only the USA, every year about 500 tornadoes (Kunkel et al., 2013) of intensity EF1 (enhanced Fujita scale; Fujita, 1971; Potter, 2007) or stronger occur, producing an average of 125 victims and huge devastation (Ashley, 2007). Numerical simulations of the very fine spatial and temporal scale of tornadoes (typically with a diameter of less than 2 km and a duration of less than 1000 s) require resolutions that are orders of magnitude smaller than those currently available in operational weather prediction and climate models (Yokota et al., 2018). Further, the chaotic dynamics of these vortices limit their deterministic prediction (Markowski, 2020). Consequently, climatological studies focused on the identification of the environmental conditions favourable to tornado-spawning severe convective storms. Several thermodynamic and kinematic meteorological parameters have been analysed, either individually or considering combined instability indices, to identify the conditions most favourable to the genesis of tornadoes (Brooks et al., 2003; Romero et al., 2007; Taszarek et al., 2018, 2020; Ingrosso et al., 2020; Bagaglini et al., 2021). This approach is consistent with the basic idea that tornadoes result from a multi-stage process, which takes into account that the tilting of the horizontal vorticity near the ground by a violent updraught plays a basic role (Rotunno, 2013; Davies-Jones, 2015). Such a conceptual model is used here as a background framework for introducing an analytical formula for the probability of tornado occurrence. A previous study defined a tornado index limited to the USA based on a Poisson regression between the observed U.S. climatology of tornadoes and monthly averaged environmental parameters from reanalysis (Tippett et al., 2012). Other studies limited their conclusions to the identification of the conditions that are associated with mesoscale convective hazards (Brooks, 2013; Diffenbaugh et al., 2013). The expression that we propose in this study is meant to provide a tool for supporting tornado warning in operational weather predictions and estimating changes in the frequency of tornado occurrence in climate projections.

2 Data and methods

Our analysis is based on tornadoes that occurred in the USA (dataset provided by the Storm Prediction Center (SPC), https://www.spc.noaa.gov/wcm/#dat, last access: 4 June 2023) and in Europe (dataset provided by the European Severe Weather Database (ESWD), https://www.essl.org, last access: 4 June 2023, managed by the European Severe Storm Laboratory (ESSL); Dotzek et al., 2009). We considered only tornadoes of category 2 or higher (F2+), following the idea that weak events might have an uncertain signature in the environmental conditions and their reporting in official databases is less accurate. A total number of 3073 tornadoes have been considered in this study (2632 for the USA and 441 for Europe; see the Supplement for density plots) during the period 2000–2018. Unfortunately, our dataset does not allow us to differentiate supercellular tornadoes from landspouts in most cases. The hourly fields of ERA5 (ECMWF Reanalysis 5; Hersbach et al., 2020) are used to establish a statistical link between the occurrence of tornadoes and a set of meteorological variables, allowing a straightforward physical interpretation of the results: the updraught maximum parcel vertical velocity (WMAX), which depends on the convective available potential energy (CAPE), the mid-level wind shear (WS₇₀₀), the low-level storm relative helicity (SRH₉₀₀), and the lifting condensation level (LCL; Kaltenböck et al., 2009). The Supplement reports the expressions defining the variables used in this study. The values of these variables have been extracted in the period 2000–2018 in all cells where at least one tornado occurred, considering the hourly reanalysis fields at 25 km resolution. The values corresponding to the occurrence of tornadoes have been selected considering the time step closest to the recorded time of the tornado onset in the database.

The univariate analysis of the (conditional) probability P of tornado occurrence is carried out by partitioning the observed range spanned by each variable into 17 equiprobable sub-intervals (bins). Such a number has been chosen as a compromise between the need of a number of bins sufficient for robust regressions and of a number of observations in each bin sufficient for a robust statistical analysis. An empirical estimate of the probability of tornado occurrence, conditional to the fact that the value of the variable lies in a given bin, is computed as the relative frequency of tornadoes in the bin. Its uncertainty is estimated via a suitable bootstrap (Monte Carlo) procedure. An analytical expression of y=log ₁₀P is found by a simple linear regression for WS₇₀₀, SRH₉₀₀, and LCL, as well as by a non-linear regression for WMAX (see the Supplement). Notice that first the climatology of the variable of interest is calculated via the partition mentioned above, and then it is compared with the tornadic cases (an approach similar to the one adopted in Romero et al., 2007).

3 Results

The univariate analysis shows that all the four variables considered in our study (i.e. WMAX, WS₇₀₀, SRH₉₀₀, LCL) are significantly linked to the formation of tornadoes. However, the formulas involving WS₇₀₀ and WMAX, i.e.

$$\begin{array}{}\text{(1)}& {\displaystyle }{\log }_{\mathrm{10}}P=-\mathrm{6.8}+\mathrm{0.11}{\text{WS}}_{\mathrm{700}},\text{(2)}& {\displaystyle }{\log }_{\mathrm{10}}P=-\mathrm{6.9}+{\displaystyle \frac{\text{WMAX}}{\mathrm{3}+\mathrm{0.32}\text{WMAX}}},\end{array}$$

describe a range of probabilities (from 10⁻⁷ to 10⁻⁴) wider than that spanned by SRH₉₀₀ and LCL. In the case of WS₇₀₀, the probability increases exponentially over the whole range. Instead, the behaviour of log ₁₀P as a function of WMAX is non-linear and shows a hyper-exponential increase in P for low values (WMAX <10 m s⁻¹), when the probability is small (about 10⁻⁷); in the intermediate range the growth gradually slows down, and P becomes quasi-constant for large values (WMAX >30 m s⁻¹), where the probability tends to $\approx {\mathrm{10}}^{-\mathrm{4}}$. For LCL and SRH₉₀₀, the exponential decrease and increase, respectively, only describe a narrow range of probability (approximately from 10⁻⁶ to 10⁻⁵). In other words, variations of these two variables do not allow us to discriminate between the low and high probability of occurrence of tornadoes as effectively as in the case of WS₇₀₀ and WMAX (see Fig. 1).

Figure 1 Univariate probability distribution for WMAX, WS₇₀₀, SRH₉₀₀, and LCL. Markers and whiskers denote the empirical probabilities with uncertainty range. Lines denote the empirical estimates (continuous) with uncertainty ranges (dashed). Different colours represent values based on the full dataset (USA and EU, black), the USA data only (red), and the European data only (EU, blue). Uncertainty ranges correspond to a 95 % confidence level.

Concerning the bivariate analysis (i.e. considering the joint behaviour of pairs of predictors), in analogy with the univariate case, a 17×17 grid matrix is constructed to partition the whole two-dimensional domain in cells. The empirical estimate of the (conditional) probability P of tornado occurrence, provided that the pair of variables lie in a given cell, is empirically computed as above via the estimate of the relative frequency of occurrence. Six different bivariate analyses are carried out considering all possible pair combinations of WMAX, WS₇₀₀, SRH₉₀₀, and LCL. For the bivariate probability, non-linear expressions have been adopted for all the pairs of variables involving WMAX and a multiple linear expression for the remaining pairs (see the Supplement). The values of the parameters of the bivariate probability functions have been estimated by a regression of the proposed expressions over the empirical probabilities.

Considering the bivariate expression of P as a function of the pairs (WMAX, LCL) and (WS₇₀₀, SRH₉₀₀), the second variable lacks significance, meaning that it provides information analogous to the first one of the pair (in fact, they are fairly correlated), but the first variable provides more (univariate) informative details than the second one in terms of the range of P. Considering the pairs (WMAX, SRH₉₀₀), (WS₇₀₀, LCL), and (SRH₉₀₀, LCL), the probability of tornadoes significantly depends on both variables, but they describe variation in P only over 2–3 orders of magnitude, whereas using the pair (WMAX, WS₇₀₀) shown in Fig. 2 it is possible to discriminate between conditions where the probability ranges from 10⁻⁷ to 10⁻³ (see the Supplement for the figures regarding all the other pairs). In conclusion, a valuable fit of the probability of occurrence of tornadoes over the range 10⁻⁷–10⁻³ is

$$\begin{array}{}\text{(3)}& {\log }_{\mathrm{10}}P=-\mathrm{6.6}+{\displaystyle \frac{\text{WMAX}}{\mathrm{3.1}+\mathrm{5.2}\text{WMAX}/{\text{WS}}_{\mathrm{700}}}}.\end{array}$$

Figure 2 Bivariate probability distribution for $(X_{\mathrm{1}}=\text{WMAX},X_{\mathrm{2}}={\text{WS}}_{\mathrm{700}})$. The coloured surface shows the empirical fit of y=log ₁₀P. Upward and downward triangles represent empirical estimates located above and below the fitted surface. All values are reported according to the colour bar.

All parameters of the univariate fits in Fig. 1 and bivariate ones in Fig. 2 are statistically significant and significantly different from zero, since the p values of the corresponding tests are (much) smaller than 1 %. For all univariate linear regressions, the adjusted R² is larger than 90 %, and, in general, the goodness of the fits is visually confirmed by the overwhelming fraction (from 90 % to 100 %) of probability values within the 95 % confidence bands. In the bivariate case, considering the multiple linear regressions of the pairs (WS₇₀₀, SRH₉₀₀), (WS₇₀₀, LCL), and (SRH₉₀₀, LCL), R² is, respectively, 70 %, 72 %, and 54 %: in general, these are smaller than in the single-variable case, but this is justified by the fact that the residual variances are about 3 times larger than those estimated in the univariate case. For the three pairs involving WMAX, R² cannot be used to assess the goodness-of-fit because the regression is non-linear. However, a slice analysis of the fits (see the Supplement for details) shows that the proposed models provide valuable fits over the whole domain of interest.

4 Discussion

Further investigations are required to ensure the validity of the expressions in Eqs. (1), (2), and (3) in different environmental and geomorphological conditions. Hypothesis-testing the similarity of the populations of tornado probabilities P_EU and P_USA, obtained using only EU and only USA data, respectively, has been carried out by using a Kolmogorov–Smirnov-like (KS) approach (Lopes, 2011) adopting the metric $d_{\mathrm{0}}=max|P_{\mathrm{EU}}-P_{\mathrm{USA}}|$. The significance level of the difference is assessed by computing the fraction of statistics exceeding d₀ using a Monte Carlo permutation procedure. Considering the univariate models, the null hypothesis that P_EU and P_USA, as a function of WMAX and WS₇₀₀, are statistically compatible cannot be rejected at 95 % and 99 % levels (suggesting that Eqs. 1 and 2 are acceptable in different geographical domains), whereas it is rejected at a level larger than 99 % for P_EU and P_USA as a function of SRH₉₀₀ and LCL. Considering the bivariate conditional probabilities, the null hypothesis – that P_EU and P_USA are statistically compatible – could not be rejected (at a 90 % level) only for the pair (WMAX, SRH₉₀₀). In this case, the overall conditional probability (combining USA and EU data) is

$$\begin{array}{}\text{(4)}& {\log }_{\mathrm{10}}P=-\mathrm{6.6}+\mathrm{0.34}{\text{WMAX}}^{\mathrm{0.37}}|{\text{SRH}}_{\mathrm{900}}{|}^{\mathrm{0.12}}.\end{array}$$

For all other pairs the null hypothesis could be rejected at the 99 % level.

Possible explanations of the lack of compatibility between conditional probabilities obtained using the EU and USA datasets alone could be different tornado damage-reporting practices (leading to different counting and attributions of tornadoes to the EF/F scale) and different meteorological and/or morphological conditions in the two domains. In spite of these limitations, as well as the need for further investigations, the proposed statistical models suitably fit the conditional probabilities of tornado occurrence. In particular, Eq. (3) has the merit of fitting the bulk of all available data and Eqs. (1), (2), and (4) of being robust with respect to the considered geographical domains.

The formulas of Eqs. (1)–(4), and particularly the bivariate expressions of Eqs. (3) and (4), outline a new statistical tool that can be used for diagnosing the likelihood of tornadoes with potential applications to short–medium range weather predictions and future changes in their frequency in climate projections. Former results considered monthly average probability (Tippett et al., 2012) or provided a modest fit to the data and were based on a smaller dataset (Cohen et al., 2018). The closest analogue to our approach is the formula of tornado probability in Grieser and Haines (2020), who considered two parameters: one describing vertical changes in temperature and a composite parameter merging CAPE and wind shear. Our results differ from Grieser and Haines (2020) in the adopted methodology for estimating the probability of occurrence of tornadoes. Grieser and Haines (2020) propose a linear regression of the logistic function, whereas we propose a non-linear bivariate fit of the logarithm of the probability. In addition, our study shows that the relationship of CAPE to the probability of tornado occurrence departs significantly from a linear dependence and that the interaction between the action of CAPE and wind shear in the lower troposphere cannot be adequately represented by their additive combination, further expanding the outcomes of Grieser and Haines (2020). Finally, Grieser and Haines (2020) used their formula for estimating past occurrence rates of tornado occurrences, while, to our best knowledge, this is the first time that analytical expressions in the form of Eqs. (3) and (4) are proposed in the scientific literature with the general aim of describing probability of tornadoes at high time and space resolution with applications in weather forecasting and climate projections.