A Statistical Investigation of Mesoscale Precursors of Significant Tornadoes: The Italian Case Study

Abstract

In this study, mesoscale environments associated with 57 significant tornadoes occurring over Italy in the period 2000–2018 are analyzed. The role of the vertical Wind Shear in the lower and middle troposphere, in terms of low-level shear (LLS) and deep-level shear (DLS), and of the convective available potential energy (CAPE) as possible precursors of significant tornadoes is statistically investigated. Wind shear and CAPE data are extracted from the ERA-5 and ERA-Interim reanalyses. Overall, the study indicates that: (a) values of these variables in the two uppermost quartiles of their statistical distribution significantly increases the probability of tornado occurrences; (b) the probability increases for increasing values of LLS and DLS, and (c) is maximum when either wind shear or CAPE are large. These conclusions hold for both the reanalysis datasets and do not depend upon the season and/or the considered area. With the possible exception of weak tornadoes, which are not included in our study, our results show that large wind shear, in the presence of medium-to-high values of CAPE, are reliable precursors of tornadoes.

1. Introduction

Tornadoes are among the severe weather phenomena producing the greatest damage and casualties [1,2]. Although the highest frequency of tornado occurrence is recorded in the United States, they are also common in southern Argentina, southern China, South Africa and Europe [3,4,5]. Tornadoes are also relatively common in Italy, occurring mostly in flat terrains in northern and southern plain areas, and along the coast of the Tyrrhenian Sea [6,7,8,9]. The intensity of tornadoes can be assessed using different scales. In 1971, the Fujita (F) scale was introduced [10], based on the damage investigation. Unfortunately, it is not suitable for European buildings [11,12,13]. In 1975, the Torro (T) scale was proposed for constructions in the United Kingdom [14]. Although it is more detailed than the F-scale, the precise assignments of damage to this scale have been considered doubtful [15]. The Enhanced Fujita scale (EF; [16]) is an update of the F-scale, providing a better estimate of the damage, but focuses mainly on buildings in the USA and Canada; moreover, recent studies have shown that wind values, corresponding to a prescribed damage, are sometimes underestimated. The EF-scale classifies tornadoes into six categories, based on their impact on vegetation, buildings and vehicles, ranging from EF0 (weak damage) to EF5 (almost complete destruction). Other solutions, e.g., the definition of specific damage indicators for Europe in the EF-scale [11] or extensions of the F-Scale outside the USA and Canada [12,13] are not yet in widespread use. Hence, in the absence of a better solution, EF-scale is still adopted here.

The percentage of “significant” (meaning EF2 or stronger—hereinafter EF2+) tornadoes in Italy is smaller than in the USA and in continental Europe [4,8]. However, some cases have caused fatalities and relevant damage in the recent past. For instance, an EF3 Tornado in the Taranto area (southern Italy), in November 2012, was responsible for one casualty and estimated damage amounting to about 60 million euros [17,18,19]. An EF4 Tornado affecting Mira and Dolo (near Venice, northern Italy), in July 2015, caused one death and 20 million euros of damage [20,21,22].

In the past, several studies have investigated the environmental conditions associated with the occurrence of tornadoes. Some variables such as the vertical wind shear (WS) in the lower and middle troposphere, the convective available potential energy (CAPE), the tropospheric relative helicity and the low-level moisture have been considered to be relevant indicators of the presence of environmental conditions favorable to the development of severe convective weather [23,24]. Rotunno and Klemp [25] and Klemp [26] stressed the important role of the pre-convective vertical WS in triggering mesocyclones and in their evolution. The tilting of the horizontal vorticity associated with the WS is an important factor to create the mid-level rotation typical of mesocyclones [27,28,29,30]. Weisman and Rotunno [31] suggest that the life of the storm, its rotation and propagation are connected to nonlinear interactions between the updraft and the environmental shear, while the helicity-streamwise vorticity is less relevant. Markowski and Richardson [32] show that tornadic supercells are favored by large low-level wind shear and high boundary layer relative humidity.

Other works have investigated whether the actual occurrence of tornadoes is associated with large values of relevant environmental variables, using a variety of datasets. Brooks et al. [3] found a reasonable consistency between the spatial distribution of large CAPE, mixing ratio, lapse rate values, obtained by combining data from the NCAR/NCEP reanalysis, and the observed distribution of tornadoes.

A statistically significant relationship between the deep-layer WS and tornado intensity is found in Colquhoun and Shepherd [33], while the tropospheric relative helicity and midtropospheric relative humidity appear to be less important. Thompson et al. [34] find larger values of WS in tornadic supercells as compared to non-tornadic supercells and non-supercell storms. Large values of vertical WS during significant tornadic storms in Canada are also found in [35].

Although Dupilka and Reuter [35] do not identify any correlation for tornadoes in Canada, most studies support the role of CAPE, though with different levels of relevance. Taszarek et al. [36] find that the probability of tornado occurrence is higher when values of CAPE (in terms of WMAX) and deep-layer shear are larger. Grünwald and Brooks [37] highlight that the combination of WS and CAPE could be a valuable discriminator between weak and significant tornadoes in the European region, since intense tornadoes occur when both WS and CAPE take on large values. Higher values of vertical wind shear, as well as of storm relative helicity, are detected in Púčik et al. [38] for significant tornadoes in Europe, whereas CAPE seems to discriminate only between no tornado and weak tornado occurrences. Rodriguez and Bech [39] focus on tornadoes in Catalonia, considering different categories of upstream vertical profiles, associated with dry conditions, non-tornadic thunderstorms and tornadic thunderstorms. They find that large values of CAPE cannot explain the development of tornadoes and waterspouts, but are useful to discriminate with respect to conditions favorable to deep convection. Moreover, they show a direct proportionality between the values of WS and tornado intensity. Similar results are shown by Groenemeijer and van Delden [40] for Dutch tornadoes. On the other hand, Thompson et al. [41] stress that CAPE may not represent an effective discriminator between storm modes, associated with tornadoes, and disorganized cells. However, Thompson et al. [34] and Craven and Brooks [42] observe a monotonic increase of CAPE from non-tornadic to tornadic supercells, and note that large values of CAPE are associated with supercells characterized by strong tornadoes.

The heterogeneous evaluations of the role of CAPE may be related to the different environmental conditions associated with severe weather in the USA, where typically values of CAPE are larger than in Europe. Indeed, Brooks [24] shows that larger values of WS and CAPE generally characterize severe weather conditions in the USA as compared to Europe. At the same time, he also claims that in the USA, WS appears to be more discriminant for the development of tornadoes than CAPE. The lower values of CAPE during tornadoes in Europe is also stressed in several other papers [43,44]. Considering the Italian case study investigated here, the importance of large values of low-level WS and CAPE has been discussed for the event that affected Taranto (southeastern Italy) on 28 November 2012 [17,18].

This work intends to continue the above-mentioned investigations and aims at checking and clarifying the relative importance of CAPE and WS on the example of significant tornado reports in Italy. A thorough statistical investigation is carried out to assess whether increasing values of WS and CAPE may correspond to increasing probabilities of tornadoes. Furthermore, an original copula analysis of the multivariate random behavior of the possible precursors of tornado occurrence is carried out (see Section 3.4). As a benefit, for the first time, it is possible to both (i) qualify and (ii) quantify the dependence (or independence) of these variables. In fact, on the one hand, an appropriate copula family (passing suitable goodness-of-fit tests) is identified as a dependence structure modeling the joint random behavior of the variables of interest. On the other hand, the value of the parameter of the fitted copula objectively measures the amount of dependence between the variables, and provides valuable information about its strength. We anticipate here that as a result and a benefit of the copula analysis, the particular structure of the fitted copula indicates that the WS precursors may be equally important, being the corresponding copula symmetric. In addition, since WS and CAPE are statistically independent (via the so-called Independence or Product copula), it turns out that the tornado-genesis is independently controlled by WS and CAPE. The analyses are based on a set of tornadoes recorded in Italy in the period 2000–2018. These events are associated with the values of WS and CAPE observed in the hours immediately before their occurrence, by using two different reanalysis datasets, both developed at the European Center for Medium-Range Weather Forecasts (ECMWF): ERA-Interim [45] and ERA-5 [46]—see Section 2.2.

The data used, and the methodology adopted here, are introduced and sketched in Section 2. The analyses and the discussion of the results are presented in Section 3. Conclusions are drawn in Section 4.

2. Materials and Methods

2.1. Records of Tornadoes

The data used in this study have already been validated in [8] for the analysis of tornadoes in Italy (limited to the period 2007–2016), and they represent an extension of this database for the period 2000–2018. These data integrate the European Severe Weather Database (ESWD)—managed by the European Severe Storm Laboratory (ESSL) (https://www.essl.org/cms) with other reports by weather amateurs, national and local newspapers, and websites including videos and pictures, as discussed in Section 3 of Miglietta and Matsangouras [8].

This analysis is restricted to EF2+ tornadoes, since we expect sufficiently large deviations of the environmental conditions from the mean climate if only the most intense events are selected. Furthermore, we only consider tornadoes for which the time of occurrence is available.

Figure 1 shows the geographical distribution, category and season of the tornadoes considered here. The dataset consists of 57 EF2+ events: 15 occurred in spring, 19 in summer, and 23 in autumn. The winter season is discarded because too few significant tornadoes are present, which makes the corresponding seasonal characterization statistically very weak.

2.2. Reanalysis Datasets

The data used in this work are extracted from two global reanalysis datasets: ERA-Interim [45] and ERA-5 [46], which, respectively, have an approximate resolution of 80 km (0.75°) and 31 km (0.25°), with 60 and 137 vertical levels from the surface up to 0.1 hPa. The data are extracted every 6 h for the u (West-East) and v (South-North) components of WS and for CAPE. Although in principle shorter time steps are possible (since hourly ERA-5 data are available), for the sake of comparison between the two datasets, a 6 h time resolution of CAPE in ERA-Interim [47] dictated such a choice. A known error in the way ERA-Interim calculated CAPE produces zero values in the 3 h forecast, and prevents to download CAPE at regular 3 h time steps https://confluence.ecmwf.int/display/CKB/ERA-Interim+known+issues).

2.3. Meteorological Parameters

As shown in Section 1, the scientific literature indicates the vertical wind shear (WS) and the convective available potential energy (CAPE) among the most relevant variables for the formation of tornadoes; for this reason, in the following, we will consider them as the variables of interest. Indeed, in mesocyclonic tornadoes, the vertical vorticity is partly generated by the tilting of the horizontal component of vorticity that is associated with vertical WS, while CAPE is a good predictor of the maximum updraft speed [23]. This represents a first step of the analysis, with the aim of extending the investigation to other variables that may be equally informative about the possible occurrence of tornadoes (e.g., storm relative helicity).

For each event, WS and CAPE are calculated at the closest grid point on the upstream side of the location where the tornado has occurred. Four Time-Steps (TS) before the event are considered for each tornado. Taking into account the 6 h resolution at which data are extracted, this allows checking the temporal evolution of the parameters of interest in the 24 h before the tornado. The Time Step closest to the tornado occurrence (i.e., from 0 to 6 hours before) is denoted as “TS1”, while “TS2”, “TS3”, and “TS4” correspond to, respectively, time ranges of 6–12, 12–18, and 18–24 h before the event. While 24 h may appear too long for fast-changing variables, such as CAPE, one should consider that we plan to extend our analysis to climate simulations; thus, it appears natural to work with the same time resolution available in such runs, and see if any signal is present within this time frame.

• Wind shear (LLS, DLS). WS is computed as the magnitude of the vector representing the difference between the wind near the surface and at a fixed level in the air column. Two parameters are considered, the low-level shear (LLS) and the deep-level shear (DLS), representing the differences between the wind speed near the surface (1000 hPa) and at about 1 km and 6 km (900 hPa and 500 hPa, respectively). Although LLS and DLS are not calculated within the generally used 0–1km and 0–6 km layers, our different choice does not affect the results, also considering that all the analyzed tornadoes occur in flat regions (only one case is recorded above 300 m). It is important to highilight that previous studies [48,49] reveal that atmospheric reanalysis tends to underestimate the low-level shear compared to sounding observations.
• Updraft Maximum Vertical Velocity (WM). CAPE [50,51,52] represents the available potential energy. It is defined as $\mathrm{CAPE}=g{\int }_{LFC}^{EL}Bdz$, where z is the height, $B=g\frac{T_{Vp}-T_{Ve}}{T_{Ve}}$ is the buoyancy force, $LFC$ is the Level of Free Convection, $EL$ is the Equilibrium Level, and $T_{Vp}$ and $T_{Ve}$ are the virtual temperature of the air parcel and of the surrounding environment, respectively [23,53].
  The buoyancy produces the vertical acceleration of the air parcels: $B=dw/dt$, where w is the vertical velocity component and t is time. In the reanalyses used in this work, CAPE is extracted from ECMWF, calculated as the maximum value among the air parcels lifted from different model levels below 350 hPa (MUCAPE). In this study, CAPE is expressed via the Updraft Maximum Vertical Velocity (WM), supposing that the buoyancy is the only active force according to the Parcel Theory: i.e., $WM=\sqrt{2·\mathrm{CAPE}}$ [23,53]. WM is computed at the same grid points and times as WS01 and WS06.
• Standardized values. To obtain comparable dimensionless indices, a standardization is applied. The generic index I is defined as $I=\frac{X-\overline{x}}{\sigma }$, where X is the value of the parameter of interest (WM, LLS, or DLS). The mean value $\overline{x}$ (climatological mean), as well as the standard deviation $\sigma$, are computed separately for each tornado, considering the 19 values of the parameter I at the same hour and calendar day when the tornado occurred for the years 2000–2018. Therefore, $\overline{x}$ and $\sigma$ have different values for each individual tornado, both for the ERA-5 and the ERA-Interim datasets.

The time series of WM, LLS, and DLS at TS1, are plotted in Figure 2 for the ERA-5 dataset (the ERA-Interim data show similar features—see Figure A1).

2.4. Statistical Tools and Procedures

Below, we briefly describe the adopted statistical tools and procedures.

2.4.1. Homogeneity Tests

The Kolmogorov—Smirnov (KS) and the Anderson–Darling (AD) are two-sample homogeneity tests commonly used to check whether two variables have the same probability distribution [54,55]. They allow rejection of the Null Hypothesis “${\mathcal{H}}_0$: the probability distribution of the variables extracted from two different samples is the same”. The KS test is indicated to spot possible differences in the body of the distributions, while AD is suggested for detecting possible dissimilarities in the tails. Instead, the Mann–Whitney (MW) test [56] is used to check whether the difference in central tendency (such as the mean, the median and the mode) of two populations differs from zero. These tests are used in Section 3.1 and Section 3.2.

2.4.2. Conditional Probabilities

An attempt to estimate the probability of tornado occurrence, conditional to the fact that WM (or LLS, or DLS) takes on values in a given range just before the event (i.e., at time step TS1), is carried out.

For each parameter (WM, or LLS, or DLS) the range spanned by the variable is first partitioned into K sub-intervals ${\Delta }_i$’s, with $i=1,\dots ,K$. Then, the number of times $N_i^P$ that a specific parameter takes on a value in the interval ${\Delta }_i$ is computed, for $i=1,\dots ,K$. Finally, the number of tornadoes $N_i^T$ when the parameter has taken on a value in ${\Delta }_i$ is calculated. In turn, the conditional probability of interest is given by

$$p_i^X=\frac{N_i^T}{N_i^P},i=1,\dots ,K,$$

(1)

where X denotes either WM, or LLS, or DLS. The uncertainty affecting the estimates of the $p_i$’s is assessed via suitable Bootstrap (Monte Carlo) procedures [57].

Since no physical indications exist that could help choosing the values for ${\Delta }_i$’s, here a reasonable statistical strategy is adopted. Specifically, the K sub-intervals are picked out in such a way that they contain approximately the same number of tornadoes. The choice $K=4$ represents a compromise between a significant number of bins and a sufficient number of tornadoes in each sub-interval: practically, they correspond to the first, second, third, and fourth inter-quartile range of the parameter of interest. For convenience, the bins are labelled as “Low”, “Medium”, “High”, and “Extreme”. The ${\Delta }_i$’s used in this work are shown in

Table 1. The sub-intervals ${\Delta }_i$’s used to partition the ranges spanned by WM, LLS, and DLS, for both the ERA-5 and the ERA-Interim datasets: the units are in m/s. The values in brackets indicate the number of tornadoes recorded in each bin.

ERA-5
	${\Delta }_1$Low	${\Delta }_2$Medium	${\Delta }_3$High	${\Delta }_4$Extreme
	[15]	[14]	[14]	[14]
WM	$(-\infty ,28.9]$	$(28.9,39.9]$	$(39.9,51.8]$	$(51.8,+\infty )$
LLS	$(-\infty ,5.1]$	$(5.1,7.9]$	$(7.9,11.4]$	$(11.4,+\infty )$
DLS	$(-\infty ,15.7]$	$(15.7,19.5]$	$(19.5,23.9]$	$(23.9,+\infty )$
ERA-Interim
	${\Delta }_1$Low	${\Delta }_2$Medium	${\Delta }_3$High	${\Delta }_4$Extreme
	[14]	[14]	[13]	[14]
WM	$(-\infty ,22.1]$	$(22.1,36.3]$	$(36.3,45.8]$	$(45.8,+\infty )$
LLS	$(-\infty ,3.7]$	$(3.7,6.9]$	$(6.9,9.1]$	$(9.1,+\infty )$
DLS	$(-\infty ,12.6]$	$(12.6,16.6]$	$(16.6,20.7]$	$(20.7,+\infty )$

. (A different number of events characterizes the datasets ERA-5 and ERA-INTERIM: respectively, 57 and 55. The reason is that the lower resolution in ERA-INTERIM does not allow differentiation of some events that occurred at close times in similar locations.)

2.4.3. Multivariate Analysis via Copulas

Copulas represent a convenient and efficient tool to model the joint random behavior of a set $X_1,\dots ,X_d$ of variables (here, the trivariate vector (WM, LLS, DLS)) ruling the statistics of the phenomenon of interest (here, the occurrence of a tornado). For a theoretical introduction to copulas, see Nelsen [58], Joe [59], Durante and Sempi [60]; for a practical approach see Salvadori et al. [61].

In particular, copulas are multidimensional functions mathematically formalizing the statistical dependence among the variables of interest. According to the Sklar’s Theorem representation [62]—i.e., Equation (2) below, copulas allow construction of the joint distribution of the $X_i$’s, with marginal distributions $F_i$’s, as

$$\mathbf{F}(x_1,\dots ,x_d)=\mathbf{P}\left(X_1\le x_1,\dots ,X_d\le x_d\right)=\mathbf{C}\left(F_1\left(x_1\right),\dots ,F_d\left(x_d\right)\right),$$

(2)

for all $(x_1,\dots ,x_d)\in {\mathbf{R}}^d$, where $\mathbf{C}$ is the copula associated with the random vector $(X_1,\dots ,X_d)$; here, $d=3$.

In the present case, according to the Kendall and Spearman independence tests [58,61], WM is statistically independent of both LLS and DLS, whereas LLS and DLS are dependent on each other. As a consequence, a trivariate copula model $\mathbf{C}$ for the vector (WM, LLS, DLS) can be set up as [61]

$$\mathbf{C}(a,b,c)={\mathbf{C}}_{\mathrm{WS}}(a,b)·c,$$

(3)

for all $(a,b,c)\in {[0,1]}^3$ (the natural domain of a trivariate copula, since its arguments are probability distribution functions ranging in $[0,1]$), where ${\mathbf{C}}_{\mathrm{WS}}$ denotes the (marginal) bivariate copula of (LLS, DLS).

Here, ${\mathbf{C}}_{\mathrm{WS}}$ can be taken as a member of the Frank family of copulas [58,61], since it adequately fits the (LLS, DLS)’s data: in fact, the Goodness-of-Fit p-Value is larger than 10%, for both the ERA-5 and the ERA-Interim reanalyses. The analytical expression of ${\mathbf{C}}_{\mathrm{WS}}$ is presented in Appendix A.

3. Results

3.1. Relevance of Geographical and Seasonal Features

Figure 1 shows that in northern Italy, tornadoes typically occur in spring and summer, whereas those affecting the Tyrrhenian coast and southern Italy are more common in autumn. However, according to the KS test, WM (as well as LLS and DLS) shows no statistically different behavior neither among the three sub-regions (as in Figure 3-left), nor among seasons (as in Figure 3-right), for both the ERA-5 and the ERA-INTERIM dataset (the latter plotted in Figure A2). The p-values reported in

Table 2. Statistical comparison of the samples of WM, LLS, and DLS collected in different regions and seasons: both ERA-5 and ERA-Interim reanalyses are considered. Shown are the p-values (in %) of (i) the KS and AD two-sample homogeneity tests (the Null assumption is that the distributions of the samples are statistically compatible), and (ii) the MW central tendency test (the Null assumption is that it is equally likely that a randomly selected value from one population is less, or greater, than a randomly selected value from a second population). Values in bold are smaller than 1%, while those in italic are smaller than 5%. The Null assumptions could be rejected at $\alpha$-levels larger than the p-values indicated.

Regional Analysis
		KS		AD		MW
	Region	ERA5	ERAint.	ERA5	ERAint.	ERA5	ERAint.
	Tyr.-South	14.30	24.59	2.10	7.03	1.70	7.74
WM	Tyr.-North	5.70	80.88	4.80	67.24	4.20	69.62
	South-North	27.70	0.99	25.50	1.73	92	4.12
	Tyr.-South	99.30	43.23	98.70	56.19	89.70	39.14
LLS	Tyr.-North	45.70	0.39	45.30	0.32	27.30	0.37
	South-North	27.70	0.06	11.80	0.01	8.60	0.01
	Tyr.-South	74.90	92.31	83.90	92.58	73.70	93.82
DLS	Tyr.-North	75.20	79.59	62.10	70.94	44.50	51.39
	South-North	26	46.31	24.50	33.06	20.70	26.64
Seasonal Analysis
		KS		AD		MW
	Season	ERA5	ERAint.	ERA5	ERAint.	ERA5	ERAint.
	Spring-Summer	86.35	87.83	66.18	91.50	51.49	100
WM	Spring-Autumn	20.68	25.62	23.03	13.24	27.29	74.25
	Summer-Autumn	31.18	13.17	11.82	5.05	10.09	27.72
	Spring-Summer	83.08	67.56	69.92	56.88	60.74	94.47
LLS	Spring-Autumn	5.69	0.06	2.20	0.02	2.38	0.02
	Summer-Autumn	6.18	0.18	1.28	0.03	0.96	0.02
	Spring-Summer	93.18	67.56	91.55	70.01	75.81	91.71
DLS	Spring-Autumn	63.24	96.53	61.28	85.21	100	90.49
	Summer-Autumn	99.31	74.32	98.42	70.73	86.11	70.76

are generally larger than 5%. Note that the AD and the MW tests lead to similar conclusions.

There is one main exception: LLS takes on values generally smaller in the northern region than in southern and Tyrrhenian ones (Figure 3-left and Figure A2-left), and it takes on values in Autumn generally larger than in Spring and Summer (Figure 3-right and Figure A2-right). These results seem consistent, since Autumn cases occur mainly in Tyrrhenian and southern regions. The visual analysis is supported by the combined p-values of the three statistical tests (KS, AD, and MW) reported in Table 2. Regarding the regions, ERA-Interim shows significant differences in LLS between the north and the other regions, which is not the case for ERA-5 (p-values are high). Therefore, in the following, the data will be pulled together and investigated, irrespectively of the region and the season of tornado occurrences analyzed. This strategy will also increase the sample size, leading to more robust statistical results.

The analysis of the Standardized indices shows to what extent the values of WM, LLS, and DLS corresponding to the occurrence of tornadoes differ from regional and seasonal mean climate conditions. In other words, the standardized indices provide somewhat different information as compared to the actual values of these parameters, since they may spotlight possible differences with respect to the mean climate when a tornado occurs. Here, the main point is that values of a standardized index larger than zero indicate that the corresponding variable takes on values larger than the average climate, and thus the crucial question is simply whether or not (and to what extent) an index is significantly larger than zero.

The boxplots presented in Figure 4 and Figure A3 show that for both the ERA-5 and the ERA-Interim datasets, the parameters WM, LLS, and DLS generally take on values (much) larger than the mean climate conditions just before the occurrence of a tornado. In some cases, all the index values are strictly positive (e.g., WM in Tyrrhenian and southern regions, and in Autumn). In general, (at least) 75% of the index values are larger than zero.

3.2. Evolution of WS and CAPE before a Tornado Occurrence

The distribution of the parameters LLS, DLS, and WM as observed at time steps TS1 and TS4—i.e., respectively, in the intervals 0–6 h and 18–24 h before the occurrence of a tornado—are statistically compared in

Table 3. Estimated p-values of the KS and AD two-sample homogeneity tests considering the distributions of the variables WM, LLS, and DLS (both actual values and standardized indices) at time-steps TS1 and TS4, for both the ERA-Interim and the ERA-5 datasets.

Actual Values
	ERA-Interim		ERA-5
	KS	AD	KS	AD
WM	0.0006%	6.1 × 10−6%	0.019%	0.00057%
LLS	0.18%	0.82%	3.8%	2.6%
DLS	3.9%	4.6%	6.4%	1.1%
Standardised indices
	ERA-Interim		ERA-5
	KS	AD	KS	AD
WM	0.0082%	0.00057%	0.16%	0.053%
LLS	1.3%	0.12%	0.016%	0.019%
DLS	3.9%	2.4%	0.65%	0.87%

, where the p-values of the KS and AD Homogeneity tests are shown, for both ERA-Interim and ERA-5 reanalyses. The results show that all three parameters take on values that immediately before the occurrence of a tornado (i.e., at TS1) are significantly different from the ones at the earlier time step TS4. In fact, all the p-values are smaller than 5%, with the only exception the KS test for DLS in ERA-5 (but the corresponding AD test yields 1.1%).

Consequently, the distributions of WM, LLS, and DLS at TS1 and TS4 turn out to be statistically different. These results, put together with those presented in Section 3.1, support the claim that large values of WM, LLS, and DLS may be robust precursors of the genesis of a tornado, especially if their magnitude is assessed with respect to the background local climatology.

3.3. Conditional Probability of Tornado Occurrence

The estimate of the probability of tornado occurrence p, conditional to the fact that WM (or LLS, or DLS) takes on a value in a given sub-interval ${\Delta }_i$ at TS1 (see Section 2.4), is plotted in Figure 5 for the ERA-5 dataset—the ERA-Interim case is illustrated in Figure A4.

A main goal of our study is to provide a statistical relation between the probability $p_i$’s and the magnitude of the considered parameters (WM, LLS, DLS). For this purpose, a linear regression is carried out, using the $p_i$’s as a response, and the sub-intervals index $i=1,\dots ,K$ as a predictor. Then, the (statistical) admissibility of a non-constant linear trend is checked via the p-Value associated with the slope of the regression line. Figure 5 shows the fitted lines, as well as the corresponding p-values. As a result, at a standard 5% level, the (conditional) probability of tornado occurrence is likely to increase as LLS (or DLS) increases (regression p-values, respectively, ≈ 4.4% and ≈ 5%), whereas the statistical evidence of a similar growth for WM (p-Value ≈ 11%) might be questionable. The analysis of the ERA-Interim datasets leads to similar conclusions.

The probability is estimated only in locations where tornadoes occurred, i.e., it is valid where morphological conditions allow their formation. In idealized conditions, Markowski and Dotzek [63] showed that the orography has the effect of enhancing convective inhibition and reducing relative humidity; accordingly, there is a tendency for the supercells to weaken in mountainous areas. Hence, although the development of tornadoes may be occasionally favored by the presence of mountains [18], they preferentially develop over flat terrain.

3.4. Multivariate Analysis of Tornado Occurrence

As anticipated in Section 2.4, the construction of a trivariate probabilistic model of the joint random behavior of the vector (WM, LLS, DLS) is carried out by using Copulas. As a novelty, the (physically) expected fact that low-level WS is more correlated with deep-layer wind shear than CAPE is given, for the first time, a (statistical) analytical multivariate model. On the one hand, the copula framework proposed here is able to “qualify” the type of dependence, via the selection of a suitable admissible copula family (i.e., different copulas identify different dependence models). On the other hand, it possible to objectively “quantify” the strength of the statistical dependence: e.g., via the estimate of the copula parameter, which is in a one-to-one relation with measures of concordance such as the Kendall and Spearman rank correlation coefficients. In summary, the copula framework outlined here provides a valuable model of the (analytical/statistical) relationships between the variables at play, and furnishes easy formulas for, e.g., setting up models of their joint behavior, including conditional relations (see Equation (3.26) in Salvadori et al. [61]).

Figure 6 shows the bivariate copulas fitted on the pairs (WM, LLS), (WM, DLS), and (LLS, DLS) for the ERA5 dataset—the ERA-Interim case is presented in Figure A5. Obviously, the full trivariate copula ${\mathbf{C}}_{\mathrm{L}\mathrm{L}\mathrm{S},\mathrm{D}\mathrm{L}\mathrm{S},\mathrm{W}\mathrm{M}}$ cannot be plotted, since its visualization would require a 4-dimensional space. Also shown are the corresponding empirical copulas, i.e., the non-parametric approximations of the true (but unknown) copulas actually ruling the joint random behavior of the pairs of interest. A comparison of the three bivariate fitted copulas is shown in the bottom-right panel (note that the two copulas fitted on the pairs (WM, LLS) and (WM, DLS) are identical, and overlap).

Overall, it is clear that there is no association between WM and LLS, as well as between WM and DLS: actually, the observations are uniformly arranged in the unit square, entailing that the corresponding copula is the one modeling independent variables—this is also supported by the large p-values of the Kendall and Spearman independence tests [61]. Instead, LLS and DLS turn out to be fairly concordant (i.e., positively correlated): small/large values of LLS are likely to be associated with small/large values of DLS (and vice-versa). The bottom-right panel of Figure 6 shows the remarkable differences between the Independence copula modeling the pairs (WM, LLS) and (WM, DLS), and the Frank copula modeling (LLS, DLS).

Figure 6 suggests that the probability of occurrence of tornadoes is independently controlled by CAPE (represented by WM) and WS (represented either by LLS or DLS). In addition, the roles played by LLS and DLS are to some extent “exchangeable”, since the corresponding copula is symmetric. In other words, there is no indication that the role of the WS as a precursor of tornadoes is better described by LLS than by DLS. This suggests that when tornadoes are formed, generally, a strong horizontal shear is present in the lower half of the troposphere, as well as in the lowest layers. The analysis of the ERA-Interim case leads to similar conclusions.

4. Conclusions

In this study, we have estimated the probability of occurrence of significant (EF2+) tornadoes as a function of three parameters: LLS, DLS, and WM. What we call here LLS and DLS are the difference between the wind vector near the surface and at the 1 and 6 km level, respectively, and are representative of the WS in the lower and middle troposphere. WM, the updraft maximum vertical velocity, is derived from CAPE. The values of these three parameters have been extracted from the ERA-5 and ERA-Interim reanalyses. The information on location and intensity of tornadoes has been provided by a dataset containing Tornadoes that occurred in Italy in the period 2000–2018.

Regional and seasonal differences of these three parameters at the time closest to the tornado occurrence appears not to be significant, and allow us to analyze all events independently of the geographical location and season.

All three parameters (WM, LLS and DLS) are important precursors of significant tornadoes. In fact, at the time closest to the tornado occurrence, their values are significantly different from those observed one day before and from the climatological means.

The values of WM, LLS, DLS have been partitioned into four bins, qualitatively representing “Low”, “Medium”, “High”, and “Extreme” values. In particular, the (conditional) probability of occurrence of tornadoes is significantly different from zero when the parameters take on “Medium” (or larger) values, and especially when WM, LLS and DLS take on values in the “Extreme” bin. Similar results are obtained with a different methodology by Taszarek et al. [36], where the highest probability of tornado occurrence corresponds to environments with high deep-layer shear and moderate CAPE. Brooks [24] found that the probability of tornado-spawning thunderstorms increases with increasing values of DLS, in agreement with our results. Furthermore, the (conditional) probability of occurrence significantly increases with the magnitude of LLS and DLS, showing that the probability increases for larger WS. The same evidence is uncertain for WM. In this case, the analysis suggests that the probability of tornado occurrence becomes significant only for “Medium” to “Extreme” values of WM, without showing an increasing trend. This result suggests that tornadoes require WM to be larger than a minimum threshold, but their probability does not increase with increasing WM.

The joint probability distribution of the parameters (WM, LLS, DLS) suggests that LLS and DLS are dependent on one other, but independent of WM. Since LLS and DLS are (statistically) exchangeable, either of the two parameters can be used a predictor of tornadoes. In addition, given the results of the multivariate analysis, the role of CAPE “alone” regarding the generation of tornadoes is somewhat peculiar, and may deserve further investigation.

This study has implicitly allowed a comparison between the new ERA-5 and the former ERA-Interim datasets, suggesting that the latter underestimates the actual values of the considered precursors [49]. Although, in general, ERA-5 tends to produce larger values of all parameters than ERA-Interim before the formation of tornadoes, the results of the statistical analysis are substantially equivalent for both databases in terms of relevance and timing of the characteristics of the considered precursors. This suggests that our results are robust and do not depend on the reanalysis dataset that is used. Many climate change simulations cannot afford the space resolution of ERA5 and are comparable to ERA-Interim. The validity of our conclusions for both datasets is a useful support for their future use in climate change simulations.

There are open issues deserving further investigation, such as refining the assessment of the dependence of tornado probability on CAPE, eventually identifying its minimum value for their occurrence. Furthermore, the validation of this probability model with a larger tornado dataset, representative of other regions in the world, and a higher temporal resolution, as allowed by ERA-5 reanalysis, will be important.

The study has produced a quantitative estimate of the increasing probability of significant tornadoes with increasing WS. It has, moreover, confirmed the relevance of CAPE and assessed quantitatively the probability of tornadoes when CAPE takes on “Medium” to “Extreme” values. Finally, it has shown the relevance of the combined values of WS and CAPE on the probability of tornadoes, showing that large values can independently support their formation.