Trend Analysis of U. S. Tornado Activity Frequency

Abstract

As one of the most severe and high-impact weather phenomena, tornadoes and their long-term frequency trends have received lots of attention from both the scientific community and the public. Here, we show that over the last six decades (1954–2018), U.S. (E)F1 tornadoes have a statistically significant upward trend but (E)F2–E)F4 tornadoes have statistically significant downward trends based on both solid trend analyses of three independent methods and robust verifications of reported tornado data using a recently developed approach called sample generation by replacement (SGR). Furthermore, we develop a statistical framework to quantitatively explain two-way interconnections between long-term climate trends and internal variabilities. With the support of quantile–quantile plots, we find that a large positive trend of U.S. (E)F1 tornadoes over the last six decades statistically results from a small internal variability of (E)F1 tornado activities. The long-term trends of U.S. (E)F2–(E)F4 tornadoes are also inversely proportional to their internal variabilities, as anticipated by the two-way interconnection theory developed in this study.

1. Introduction

Tornado activities in the United States have significant impacts on society and economics in terms of loss of life and property damages [1,2,3]. Whenever severe tornado events occur, such as the recent tornado outbreaks that happened on 10–11 December 2021 over the U.S. (Tornado outbreak of 10–11 December 2021—Wikipedia), the public and media often raise numerous questions related to tornadoes and their relationship with climate change. The Intergovernmental Panel on Climate Change (IPCC) and the U.S. Climate Change Science Program have also drawn their attention to these questions, especially on how localized severe weather activities, such as tornadoes [4,5], are changed under a changing climate, even though it is challenging to identify the influences of climate warming on tornado activity [6,7]. This research helps answer some unknown or unresolved scientific questions of concern to the IPCC, U.S. climate change program, and public. In particular, we will investigate long-term temporal changes of tornado activity at various Enhanced Fujita [(E)F]-scales (

Table 1. Enhanced F scale of tornado and wind speed (m s⁻¹).

EF Scale	(E)F 0	(E)F1	(E)F2	(E)F3	(E)F4	(E)F5
Vmedian	33.5	43.8	55.0	67.3	81.8	89.4
Speed range	29.1–38.0	38.4–49.2	49.6–60.4	60.8–73.8	74.2–89.4	>89.4

, https://www.spc.noaa.gov/faq/tornado/ef-scale.html, accessed on 16 March 2022) and explain these changes from the point of view of a developed statistical approach.

Even with ever-increasing societal interest in tornado activities engendering catastrophes of loss of life and property damages [1,2], a long-term change of tornado activities in the U.S. over the past six decades is still puzzling or not robustly proved scientifically. To date, research on trends for U.S. tornado activity (counts and days) can be grouped into the following categories (see

Table 2. Comparisons of trend analysis for U. S. tornado activity.

Authors	Entire or Regions of U.S.	Tornado Counts and/or Days	Trend Analysis Methods	Individual (E)F Scale or (E)F+	Data Coverage (Years)
Verbout et al. (2007)	Entire	Counts	Linear Regression	(E)F0–(E)F5 (E)F1+…(E)F4+	1954–2003 (50 years)
Kunkel et al. (2013)	Entire	Counts	Eye-ball Slope Line	(E)F0 (E)F1+	1950–2011 (62 years)
Moore (2017)	Entire	Days	MK test with pre-whitening	(E)F1+	1974–2015 (42 years)
Guo et al. (2016)	Regions	Counts	Linear Regression	(E)F1+	1950–2013 (64 years)
Gensini and Brooks (2018)	Regions	Counts	Linear Regression	(E)F0+	1979–2017 (39 years)
Moore (2018)	Regions	Counts	MK test with pre-whitening	(E)F1+	1954–2016 (63 years)
Our work	Entire	Counts and Days	(1) Linear Regression (2) MK test with trend-free pre-whitening (3) Monte Carlo simulation	(E)F1+ (E)F1…(E)F5	1954–2018 (65 years)

for detailed comparisons): (1) estimating trends for tornado counts over the entire U.S. by plotting the time series [8] and/or slope lines [9], implying that the occurrence of (E)F1+ tornadoes shows no trend since 1954 [8], and both (E)F2+ and (E)F3+ tornadoes have a declining trend over 1954–2003 [9]; (2) performing trend analyses for tornado days over the entire U.S. using the MK (Mann and Kendall) test [10,11] with the pre-whitening procedure [12]; and (3) performing trend analyses for tornado counts over different regions of the U.S. using a linear regression method [13,14] or using the MK test with the pre-whitening procedure [15]. The research [12] used 41 years (1974–2014) of data on U.S. tornado days and a MK test [10,11] with pre-whitening but without applying trend-free procedures. The results showed that since the mid-1970s, the number of a tornado days decreased but the mid-frequency tornado days (10–29 tornadoes) had no trend and the number of high-frequency tornado days (30+ tornadoes) had increased [12]. It should be noted that when only pre-whitening is used, part of the true trend can be removed [16,17]. Use of the same method for trend analysis of tornado counts (1954–2016) over six regions of the U.S. found that the number of tornadoes per year declined in the West, North Great Plains, South Great Plains, and Midwest regions, but increased in the Southeast [15]. However, utilizing Kendall’s $\tau$ linear correlation coefficient to assess the spatial trends of all (E)F-scale tornadoes reported over the U.S. since 1979 identified positive trends over the portions of the Midwest [14]. Contradictory to the previous two studies, the Great Plains was documented with a significantly increasing trend using a conventional linear regression method to examine the spatial trends over different regions of the U.S. tornado number since 1950 [13]. As a result, trend analyses of U.S. tornado activity are sensitive to the length of tornado data used and methods employed, for example, using different methods even for same geographic regions.

In this work, we perform trend analyses for both tornado counts and days over the entire U.S. using three independent methods [17,18,19,20] such as the MK test with trend-free pre-whitening, the Monte Carlo simulation, and linear regression (see methods), and the reported tornado data from 1954 to 2018 with the weakest tornadoes (E)F0 were removed to eliminate the reporting bias in (E)F0 tornadoes [13,21] (see methods). As a result, we can directly compare our results with those in categories (1) and (2) to find out what are the differences between our work and previous studies, and what is new in our work. By using three independent methods, we will show for the first time that over the last six decades, individual (E)F1 tornado counts have a statistically significant upward trend; individual (E)F2, (E)F3, and (E)F4 tornadoes have a significant downward trend; and (E)F5 has a downward trend but not statistically significant. In addition, we will demonstrate a downward but not significant trend for the total of the tornadoes from the (E)F1 to (E)F5 scales, which provides scientific evidence to support assumptions made in the previous studies [8,9]. To gain further insights into these identified trends, we will develop a statistical framework to elucidate how long-term climate trends (e.g., U.S. tornadoes) can be quantitatively evaluated by internal climate variabilities through their two-way interconnections, especially concerning why the (E)F1 tornado has an upward trend, i.e., why (E)F1 tornadoes occur more frequently.

We also show that a tornado day has a statistically significant negative trend, consistent with the previous study [12]. On the other hand, we demonstrate that for tornado days with high frequency (≥10+ tornadoes per day), both the total of (E)F1–(E)F5 tornado days and (E)F1 tornado days have statistically significant positive trends, indicating that the increase of high frequency tornado days is mainly driven by (E)F1 tornadoes, which is absent in the previous study [12].

Before performing trend analyses of tornado activities, it is important to build up a solid foundation on the basis of the tornado data we use. To ensure that the tornado data are reliable and trustable, we employ a method of sample generation by replacement (SGR) to examine the tornado data used in this study. The SGR method was first developed and applied for phycology research dealing with potential fake data [22,23,24]. The idea of this method is to construct a new dataset, called a fake dataset, by manipulating an original dataset based on a replacement probability distribution [22,23,24]. To our knowledge, this method has not been applied for tornado climate analysis so far.

2. Data and Methods

2.1. Observational Data

The U.S. tornado data provided by NOAA’s (National Oceanic and Atmospheric Administration) National Weather Service Storm Prediction Center (SPC) are the same as those used in other studies [6,7,19,23]. These data are maintained by the Storm Prediction Center of NOAA, reviewed by the U.S. National Climatic Data Center, and updated on a yearly basis [13] (http://www.spc.noaa.gov/wcm/#data, accessed on 16 March 2022). This tornado dataset consists of one single text file where, for each reported U.S. tornado, information such as the date, time, start location, end location, EF-scale, and fatality are recorded.

There are uncertainties in the U.S. tornado databases [14,25,26,27,28,29,30] but recent studies show that the time series of (E)F1 or greater tornadoes is more stationary over time [9] and is representative of actual tornado activity since the 1950s [6,7]. Similar to others [4,6,7,31], we eliminated (E)F0 tornadoes from consideration and excluded 4-year (1950–1953) data from the tornado time series because of incomplete efforts to archive them [31].

2.2. Sample Generation by Replacement (SGR)

Essentially, SGR is a perturbation method to perturb a dataset through a two-step procedure: (1) using Monte Carlo simulation to generate a dataset and (2) then using an ad hoc probabilistic model to replace the original dataset with the generated data in step (1). The data replacement is satisfied with a probabilistic distribution, such as a uniform support fake-good distribution [22,23,24], as follows:

$$\mathrm{P}(f_i=k|d_i=h,{\theta }_p)=\{\begin{array}{c}1, k=h=max\\ {\theta }_{p, } k>h\\ 1-{\theta }_p, k=h\\ 0, k<h=min\end{array},$$

(1)

where P is a probability of fake (perturbation) data f_i replacing observed (reported) data d_i (i = 1,…,n) with prior knowledge about the distribution of faking or empirically based knowledge about the process of faking, such as the direction of faking, e.g., fake-good vs. fake-bad. Here, the fake-good refers to all values of k > h assumed to be equally likely in a process of replacing observed h while not allowing for replacing the original observed (reported) value h with lower-value k. ${\theta }_p$ is a percent of the replacement by the perturbations. This process mimics a possible inflation of tornado-number increase due to non-meteorological reasons.

2.3. Trend-Free Approach for Removing AR(1) Process

A lag-one autoregressive process, AR(1), in a time series increases the probability of the detection of a significant trend by the MK test. To detect a deterministic trend (systematic changes over the mean), not a stochastic trend introduced by red noise represented by an AR(1) process, AR(1) should be removed [17], if it is statistically significant, from the time series through a pre-whitening procedure. However, this pre-whitening also removes a portion of the true trend [16]. Hence, we have applied a trend-free pre-whitening approach [17] prior to applying the statistical MK test. As a result, the true trend is preserved and no longer influenced by the effects of AR(1) autocorrelation.

To detect a trend properly, we use a trend-free pre-whitening approach [17] through the following four steps.

(1)

A non-zero slope β of a trend in a time series {$X_t$, t = 1,2,…,n} is estimated by linear regression, and the sample data are detrended.

$$X_t^{\prime }=X_t-T_t=X_t-\beta t$$

(2)

(2)

A lag-one serial correlation coefficient ${\rho }_1$ of the detrended series $X_t^{\prime }$ is calculated and the AR(1), if it is statistically significant, is removed from $X_t^{\prime }$.

$$Y_t^{\prime }=X_t^{\prime }-{\rho }_1{X}_{t-1}^{\prime }$$

(3)

(3)

The identified trend $T_t$ from step (1) and the residual $Y_t^{\prime }$ are blended.

$$Y_t=Y_t^{\prime }+T_t$$

(4)

(4)

The MK test is then applied to the blended time series to assess the significance of the trend. The blended series $Y_t$ preserves the true trend and is no longer contaminated by the effects of autocorrelation.

2.4. Non-parametric MK Statistical Test for Trend Detection

Under the null hypothesis H_o (i.e., no trend in a time series) that a sample of data {$X_i$, i = 1,2,…,n} is independent and identically distributed, the MK test statistic $S$ is defined as:

$$S={\sum }_{i=1}^{n-1}{\sum }_{j=i+1}^{n}sgn\left(X_j-X_i\right)$$

(5)

where

$$sgn\left(\mathsf{\theta }\right)=\{\begin{array}{c}1 if \mathsf{\theta }>0\\ 0 if \mathsf{\theta }=0\\ -1 if \mathsf{\theta }<0\end{array}.$$

(6)

The MK showed that when n ≥ 8, the statistic $S$ is approximately normally distributed with the mean and variance as follows:

$$E\left(s\right)=0$$

(7)

$${\sigma }_s^2=\frac{1}{18}[n\left(n-1\right)\left(2n+5\right)-{{\displaystyle \sum }}_{i=1}^m{t}_i\left(t_i-1\right)\left(2t_i+5\right)],$$

(8)

where m is the number of tied (i.e., equal values) groups and t_i is the number of data points in the ith tied group. Under the null hypothesis, the standardized MK statistic $Z$ follows the standard normal distribution with a mean of zero and variance of one:

$$Z=\left\{\begin{array}{cc}(S-1)/{\sigma }_S& ifS>0\\ 0& ifS=0\\ (S+1)/{\sigma }_S& ifS<0\end{array}\right.$$

(9)

If |Z| > Z_1−α/2, the null hypothesis H_o is rejected and a trend is statistically significant at a level of 1 − α/2.

2.5. Monte Carlo Simulation for Confirming Detected Trends

The original time series of tornado frequency is first shuffled using a random number generator from a standard intrinsic function in Fortran 90 (http://www.nsc.liu.se/~boein/f77to90/a5.html, accessed on 16 March 2022). A linear regression slope is then computed. After repeating these processes 1000 times, we sort out these slopes in an ascending order and obtain the 990th test statistic (the absolute value of the slope for a two-tailed test). If the slope with the original time series is greater than the 990th test statistic, for example, then the detected trend is statistically significant at the 99% level [20].

3. Results

3.1. Trend Analyses for U.S. Tornadoes (E)F1–(E)F5 over the Last Six Decades

Before carrying out trend analyses of tornado activities, it is important to ensure that the tornado data are reliable and trustable. Hence, we employed a method of sample generation by replacement (SGR) to examine the tornado data. After operating the SGR procedure, we performed a statistic t-test [17,31] at a statistically significant level of 95% to see if there is any statistically significant difference between the original data and the data perturbed through the SGR procedure (see Section 2.2). As shown in

Table 3. T test for mean differences between original and SGR perturbed tornado time series.

F Scale	(E)F1	(E)F2	(E)F3	(E)F4
t > t0.05 (=1.96)	2.15	2.10	2.10	2.19
θp	0.25	0.45	0.40	0.35

, in order to make a statistically significant difference at a level of 95%, one needs to replace 25% to 45% (i.e., θ_p) of the original tornado (E)F1–(E)F4 data through the SGR process. Any perturbations below 25–45% replacement do not lead to a significant difference between the original data and the data with the perturbations. This indicates that the original tornado (E)F1–(E)F4 data tolerate high levels of uncertainty while they still keep accurate statistical features such as a mean. With 25%, 45%, 40%, and 35% replacements for (E)F1, (E)F2, (E)F3, and (E)F4 tornado data, respectively, these tornado trends still keep the same sign as the ones without the replacements. Therefore, the tornado (E)F1–(E)F4 data are accurate and trustable to serve our trend analyses. Note that the t-test of the (E)F5 tornado is not shown because the sample size of the (E)F5 tornado is much smaller than that of the (E)F1–(E)F4 tornado and the (E)F5 tornado has no significant trend, as demonstrated in the next few paragraphs.

At a first glance, the time series of the tornadoes shows an upward trend in (E)F1 tornadoes (Figure 1a) and a downward trend in (E)F2 to (E)F5 tornadoes (Figure 1b), resulting in a slightly downward trend for the total number of tornadoes from the (E)F1 to (E)F5 scales (Figure 1c). This needs to be verified through statistical tests for the individual time series of (E)F1 to (E)F5 tornadoes.

We used three independent methods [17,18,19,20] to identify any significant trends for (E)F1 to (E)F5 tornadoes (see methods). First, we employed the nonparametric MK test after removing AR(1) and applied a trend-free pre-whitening procedure to the (E)F1 to (E)F5 tornado time series [17,19,20]. Second, we applied the conventional linear regression method to the original tornado time series [13,17]. Third, we used Monte Carlo simulation.

For MK tests, the upward trend of (E)F1 tornadoes is at least statistically significant at the level of 99%; the downward trends of (E)F2, (E)F3, and (E)F4 tornadoes are significant at the level of 99%; and the downward trend for (E)F5 tornadoes is not statistically significant (

Table 4. Mann-Kendall test for U.S. tornado (E)F1 to (E)F5 frequency trends.

F Scale	(E)F1	(E)F2	(E)F3	(E)F4	(E)F5
Z statistic	3.66	−5.83	−4.26	−3.06	−1.06
$Z_{1-\alpha /2}$	2.58	2.58	2.58	2.58	1.65
α	0.01	0.01	0.01	0.01	0.10
Trend (tornado/year)	1.99	−2.00	−0.45	−0.10	−0.01

). After removing the AR(1) processes, the trends for (E)F1 to (E)F4 tornadoes are 1.99, −2.00, −0.45, −0.10 tornadoes per year, respectively. These trends were also independently identified by the conventional linear regression approach at a statistical significance level of 99% for tornadoes (E)F1 to (E)F4. The Monte Carlo method further confirmed the detected trends of (E)F1 to (E)F4 tornadoes at the 99% significance level. Hence, the detected trends for (E)F1 to (E)F4 tornadoes are robust. Furthermore, when we add the time series of (E)F1 to (E)F4 tornadoes, these new time series have a downward but not a statistically significant (at a level of 95%) trend (Figure 2b), while (E)F1 tornadoes show a statistically significant upward trend (Figure 2a). This is also true if we further include (E)F5 tornadoes in this time series. In short, our results demonstrate that over the last six decades, (E)F1 tornadoes have a significant upward trend whereas (E)F2, (E)F3, and (E)F4 tornadoes have a significant downward trend, and (E)F5 has a downward trend but not significant. Besides, the identification of the downward but not significant trend for the total of tornadoes from the (E)F1 to (E)F5 scales provides scientific evidence to support assumptions in the previous study [8,9].

To further robust the above findings, we carried out two types of numerical experiments to test (1) if tornado trends are sensitive to population density related ratios of (E)F1 to (E)F2 and (2) if tornado trends are sensitive to major tornado outbreaks, such as onesoccurred during 25–28 April 2011 and 3–4 April 1974. It has been pointed out that (E)F2–(E)F5 tornado reports would vary less with population density compared with (E)F0–(E)F1 tornado reports [25]. Recent research [26] found that the (E)F1 and stronger tornadoes related population bias disappears more quickly over time and that the (E)F0 tornadoes have a less pronounced population bias. It is hypothesized [25] that some bias in rural areas may be due to artificially decreasing the number of (E)F2–(E)F5 tornadoes and increasing the number of (E)F0–(E)F1 tornadoes, indicating that the (E)F2 numbers need to increase while the (E)F1 numbers need to decrease. In contrast, a recent study [27] suggested increasing (E)F1 tornado numbers and decreasing (E)F2 tornado numbers by about 28% for the period of 1954–1973 using the tornado count ratio of (E)F1 to (E)F2 over 1974–1994, assuming that the latter period is more reliable. Given the uncertainty of the ratio of (E)F1 to (E)F2, we performed numerical experiments by increasing (decreasing) (E)F1 tornado numbers and decreasing (increasing) (E)F2 tornado numbers about ±14% for the period of 1954–1973. The second situation of decreasing (E)F1 tornado numbers and increasing (E)F2 tornado numbers was to consider the proposed suggestion [25]. As a result, both experiments have demonstrated an upward trend for the (E)F1 tornado and a downward trend for (E)F2 tornadoes and (E)F1–(E)F5 tornadoes (Figure 3 and Figure 4), reaching the same conclusion as before. This indicates that population bias has little effect on the sign of (E)F1 and (E)F1–(E)F5 tornado trends but might have some effect on their magnitude.

To examine if a single major tornado outbreak, such as the one during 25–28 April 2011, affected the results of the trend analysis, we removed about 223 (E)F1–(E)F5 tornadoes that occurred during this time period and found that (E)F1 tornadoes still have an upward trend at a statistically significant level of at least 99%; (E)F2–(E)F5 tornadoes have a downward trend at a statistically significant level of at least 99%; and (E)F1–(E)F5 tornadoes have a downward trend but not statistically significant (Figure 5). Furthermore, we performed a sensitivity test to remove a tornado outbreak near the beginning of the entire time series, i.e., 3–4 April 1974, with 135 (E)F1–(E)F5 tornadoes reported during this tornado outbreak. Similarly, we found that both (E)F1 and (E)F2–(E)F5 tornadoes have an upward and downward trend, respectively, at a statistically significant level of at least 99% and (E)F1–(E)F5 tornadoes have a downward trend but not statistically significant (Figure 6). As a result, we draw the same conclusions as before.

In addition to U.S. tornado count trends, we also examined U.S. tornado trends from the tornado day [6,29,30,32] perspective, ranging from low frequency (e.g., a tornado day with one tornado at least on that day) to high frequency (e.g., 10+, 20+, 30+, 40+, and 50+ tornado days).

For a tornado day, (E)F1–(E)F5 tornado days have a negative trend (Figure 7a), which is consistent with Figure 2a of the recent study [12]. In addition, the negative trend of a tornado day for total (E)F1–(E)F5 is consistently attributed to different intensities/scales of tornadoes. These include contributions by (E)F1, (E)F2, (E)F3, and (E)F4 tornado days, in which all negative trends are statistically significant at the level of >99% (Figure 7b), except for (E)F5 tornado day which have a negative trend but not statistically significant (not shown). The study [12] did not show the contributions to the negative trend of a tornado day from different intensities/scales of tornadoes rather than presenting the contributions to 1–9 tornado days. The important difference is that in Figure 5a of the study [12], the (E)F1 tornado day showed a positive trend (at the 90% level of statistical significance), whereas in our Figure 7b, the (E)F1 tornado day shows a negative trend (at the >99% level of statistical significance). These differences may be due to the different lengths of tornado data used (65 vs. 42 years) and methods employed in the two studies. In our work, we show the tornado day trend and the contributions to a tornado day from (E)F1–(E)F5 tornado days, while the study [12] showed the tornado day trend and the contributions to the 1–9 tornado days from (E)F1–(E)F3 tornado days.

For tornado days with a high frequency of occurrence (e.g., 10+, 20+, 30+, 40+, and 50+ tornadoes per day), (E)F1–(E)F5 tornado days have positive trends, in which 10+ and 50+ tornado days pass a statistical test at the level of 95%, and 20+, 30+, and 40+ tornado days pass a statistical test at the level of 99%. For 20+ tornado days, for example, (E)F1–(E)F5 tornado days have a positive trend (Figure 8a) contributed by a positive trend of the (E)F1 tornado day (Figure 8b), and both positive trends are statistically significant at the level of >99%. Our results show that (E)F2, (E)F3, (E)F4, and (E)F5 tornado days are not statistically significant. Similar results were obtained for 10+, 30+, 40+, and 50+ tornado days, with positive trends for both (E)F1–(E)F5 and (E)F1 tornado days at the statistically significant level of >95%. This indicates that the increase of tornado days is toward the high frequency phase, which is mainly contributed by the (E)F1 tornado day. On the other hand, the research in [12] pointed out without statistical tests that the high-frequency tornado days tend to have a higher proportion of stronger tornadoes, such as (E)F3.

3.2. Two-Way Interconnections between Long-Term Climate Trends and Internal Variabilities

To understand why the (E)F1 tornado occurrence frequency increases while the (E)F2–(E)F4 tornado occurrence frequency decreases in recent decades, we developed a statistical framework to quantitatively evaluate two-way interconnections between long-term climate trends and internal variabilities.

A linear trend b in a regression equation ${\widehat{y}}_i=a+b{x}_i$, where $x_i$ is a predictor variable of time, ${\widehat{y}}_i$ is a predictand, and $a$ is a constant, can be expressed as:

$$b^2=\frac{{{\displaystyle \sum }}_{i=1}^n{\left({\widehat{y}}_i-\overline{y}\right)}^2}{{{\displaystyle \sum }}_{i=1}^n{\left(x_i-\overline{x}\right)}^2},$$

(10)

where ${{\displaystyle \sum }}_{i=1}^n{\left({\widehat{y}}_i-\overline{y}\right)}^2$ is the sum of the squared differences between the regression predictions and the sample mean $\overline{y}$, and it can further be formulated as:

$$\sum _{i=1}^n{\left({\widehat{y}}_i-\overline{y}\right)}^2=\sum _{i=1}^n{\left(y_i-\overline{y}\right)}^2-\sum _{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2,$$

(11)

where the first and second terms on the right-hand side of Equation (11) are, respectively, the sum of the squared deviation of the sample $y_i$ value around the mean value $\overline{y}$, and the sum of squared residuals, respectively. Substitute Equation (11) into (10), we have

$$b^2=\frac{{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\overline{y}\right)}^2-{{\displaystyle \sum }}_{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2}{{{\displaystyle \sum }}_{i=1}^n{\left(x_i-\overline{x}\right)}^2}.$$

(12)

The sum of the squared residuals is a function of the standard deviation of the residuals σ:

$${\sum }_{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2=(n-2){\sigma }^2.$$

(13)

Substitute Equation (13) into (12), we yield

$$b^2=\frac{{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\overline{y}\right)}^2-\left(n-2\right){\sigma }^2}{{{\displaystyle \sum }}_{i=1}^n{\left(x_i-\overline{x}\right)}^2},$$

(14)

where σ is the standard deviation of the internal variability in the context of climate change. If an AR(1) process exists in the time series of the residuals, σ needs to be multiplied by a factor of γ to account for the bias in the sample standard deviation due to persistence [33]:

$$\gamma ={\left[\frac{n-2}{n\left(\frac{1-r_1}{1+r_1}\right)-2}\right]}^{1/2},$$

(15)

where $r_1$ is the AR(1) of the residuals. To eliminate dimension and magnitude differences among climate variables, we introduce a concept of relative internal variability RIV defined as

$$\mathrm{RIV}=\frac{(n-2){\sigma }^2}{{\sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2}$$

(16)

so that RIV ranges are within [0, 1] and can be used to measure and compare the internal variability for different climate variables. Substitute Equation (16) into (14), we gain

$$b^2=\frac{1-\mathrm{RIV}}{\frac{{\sum }_{i=1}^n{\left(x_i-\overline{x}\right)}^2}{{\sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2}}$$

(17)

On the other hand, by rewriting Equation (17), we have impacts of trend b on RIV:

$$\mathrm{RIV}=1-b^2\frac{{\sum }_{i=1}^n{\left(x_i-\overline{x}\right)}^2}{{\sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2}$$

(18)

Equations (17) and (18) indicate that the large (small) internal variability statistically results in the small (large) of the trend, and vice versa. Mathematically speaking, for a given dataset the sum of the total squared deviation is constant, the less (more) the sum of the squared regressions (proportional to the trend), the more (less) the sum of the squared residuals (proportional to the internal variability). Graphically speaking, the less (more) data clamp onto a trend line, the more (less) variability of data points is.

To be able to perform proper computations based on Equation (17), one needs to verify if the residuals (i.e., the deviations about a linear trend b) obey a normal distribution. To confirm this, we plotted quantile–quantile (Q–Q) plots [28] for (E)F1–(E)F4 tornadoes. As displayed in Figure 9, there are excellent linear relationships between actual quantiles of (E)F1–(E)F4 tornadoes and normal theoretical quantiles. The correlation coefficients are 0.97, 0.96, 0.95, and 0.91 for tornadoes (E)F1, (E)F2, (E)F3, and (E)F4, respectively, and all are statistically significant at a level of at least 99.5%. Therefore, the residuals for the tornadoes (E)F1–(E)F4 obey the normal distributions and applications of Equation (17) based computations to tornadoes (E)F1–(E)F4 are well justified.

In fact, the derived Equation (17) or (18) can be applied to any variables and disciplines at any spatial and time scales as long as the residuals (i.e., the deviations about a linear trend b) obey a normal distribution.

As Equation (17) shows, the increase (decrease) of tornado RIV statistically results in the decrease (increase) of its occurrence frequency trend. The (E)F1 tornado RIV 0.58 counts year⁻¹ corresponds with its trend b = 1.99 counts year⁻¹ (b² = 3.96;

Table 5. Tornado trend square b² and tornado RIV (# year⁻¹).

F Scale	(E)F1	(E)F2	(E)F3	(E)F4
b2	3.9601	4.0000	0.2025	0.0100
RIV	0.5770	0.5414	0.7669	0.8881

), whereas the RIV for (E)F4, being 0.89 counts year⁻¹ corresponds with its trend b = −0.10 counts year⁻¹ (b² = 0.01; Table 5). All tornadoes (E)F1 to (E)F4 after removing AR(1) processes very well obey the theory derived from Equation (17) that the large RIV statistically leads to the small trend (Table 5), which is also true for tornadoes (E)F1 to (E)F4 before removing AR(1) processes.

4. Discussion and Conclusions

Using three independent statistical methods, we show for the first time that over the last six decades, U.S. (E)F1 tornado counts have a significant upward trend while (E)F2, (E)F3, and (E)F4 tornado counts have a significant downward trend, and (E)F5 tornado counts have a downward trend but not significant. Based on our numerical experiments, the trends of tornado counts are not sensitive to both population density related ratios of (E)F1 to (E)F2 and to the removal of tornadoes occurred in the major tornado outbreaks of 3–4 April 1974 and 25–28 April 2011.

The effect of population density on tornado intensity estimation is also discussed in [34]. It seems that due to the increasing number of storm chasers, the population bias is less pronounced for (E)F0 tornadoes, and disappears more quickly over time for (E)F1–(E)F5 tornadoes [35].

To assess if the overrating of strong tornadoes (≥F2) before the Fujita-scale implementation in 1974 [9,27] impacts our trend analyses for individual (E)F-scale tornadoes, we have carried out (1) SGR analyses for each (E)F-scale tornado dataset and (2) Monto Carlo simulations for each (E)F-scale tornado after the trend analyses. The SGR assessments show that even with 45%, 40%, and 35% random replacements for (E)F2–(E)F4 tornado data, the (E)F2–(E)F4 tornado trends still keep the same sign. Therefore, these (E)F2–(E)F4 tornado data are trustable to serve our trend analyses. Furthermore, the Monto Carlo evaluations by randomly shuffling the original time series of (E)F2–(E)F4 tornadoes for 1000 times confirm that the detected trends for (E)F2–(E)F4 tornadoes are statistically significant at 99%. It should be pointed out that the SGR approach and Monto Carlo simulation have already taken into consideration other types of heterogeneity, such as the implementation of the Weather Surveillance Radar-1988 in 1997 [27]. In fact, those two techniques can apply to any time series with inhomogeneity.

Consistently, a tornado day has a statistically significant negative trend that is contributed by statistically significant negative trends of (E)F1, (E)F2, (E)F3, and (E)F4 tornado days, except for the (E)F5 tornado day having a negative trend but not statistically significant. For tornado days with high frequency (≥10+ tornadoes per day), both the total of (E)F1–(E)F5 tornado days and (E)F1 tornado days have statistically significant positive trends, indicating that the increase of tornado days toward the high frequency is mainly contributed by (E)F1 tornadoes. In short, the number of days with at least one tornado decreases while the number of days with many tornadoes increases. This is consistent with previous studies [12,29,36,37,38].

Furthermore, we show the two-way interconnections between climate trends and internal variability based on our newly derived model. Using this novel model, we explained why over the last six decades, U.S. (E)F1 tornado counts have a significant upward trend while (E)F2 to (E)F4 tornado counts have a significant downward trend, as detected by three independent statistical methods. Our work offers a basis for understanding and predicting two-way interconnections between trends and internal variabilities that can be used for other disciplines and variables, and could make a contribution to seasonal forecasting of tornado activities.