Observational trends

Thus far, TC intensification trend analysis has exclusively focused on International Best-Track Archive for Climate Stewardship17 (IBTrACS) data (Methods), which likely adds systematic biases to published results. IBTrACS is composed of TC best-track data from different operational agencies whose TC intensity observations became more reliable when global satellite coverage was introduced in the early 1980s (except for the Indian Ocean which did not gain continuous satellite coverage until 1998)18. Over much of the globe, satellite data quality has improved through increasing resolution and decreasing view angles, which enhances the accuracy of the most recent intensity estimates. The recent addition of more in situ measurements in the Atlantic and the East Pacific leads to lower errors in these basins, which results in additional spatial and temporal inhomogeneities in the observations19.

In order to gauge the impact of these inhomogeneities on TC intensification trends, we first compare the probability of 24-h wind speed changes in IBTrACS and a more homogeneous record of TC intensity, the Advanced Dvorak Technique-Hurricane Satellite-B1 (ADT-HURSAT)12 (Methods). Kossin et al.12 developed ADT-HURSAT in order to maintain the same protocol to determine TC intensities throughout all ocean basins during the period 1982–2009. ADT-HURSAT incorporates artificially degraded intensity estimates in certain regions so its data quality is consistent both spatially and temporally, which is better suited for trend analysis. However, this increased data consistency means the best technology and analysis techniques are not utilized in ADT-HURSAT, which leads to individual intensity estimates with higher average errors than those in IBTrACS12,20. However, for the scope of our study, which involves assessing rates of change throughout a time series, temporal consistency is preferable to accuracy in individual estimates.

Figure 1 shows the empirical probability density plots calculated from ADT-HURSAT and IBTrACS 24-h intensity changes, using a logarithm transformation on the probabilities to better illustrate differences in the tails of the distributions. Supplementary Figure 1 shows the more traditional probability density plots. For all of our analyses, we impose criteria involving storm longevity, location, and intensity to ensure only intensity changes derived from well-defined storms are included in our results (Methods). For the global results, we entirely exclude the north and south Indian Oceans from the sample because of the well-known absence of quality satellite data in the region until 199812. Nevertheless, our conclusions do not significantly change when these basins are introduced into the analysis.

Fig. 1 Observed 24-intensity change probability densities. a, b Common logarithm of the probability densities calculated from IBTrACS (black) and ADT-HURSAT (blue) 24-h intensity changes. a Global and b Atlantic basin results for the period 1982–2009 are plotted. Data are binned in 20 knot increments between −110 and 110 knots. The dashed lines indicate the 90% confidence interval (between 5th and 95th percentiles of the data) (Methods). All distributions are bounded below by 10−5 Full size image

As expected, the differences in the ADT-HURSAT and IBTrACS intensity change distributions are larger for the global data, while the two datasets show more agreement in the Atlantic basin. For the global results, TCs in ADT-HURSAT maintain their intensity for a larger percentage of cases than in IBTrACS but also have significantly more instances of intensification rates between 50 and 90 knots (resulting in a “shelf” in this part of the distribution). IBTrACS curves appear almost linear in both plots which imply that probability densities are exponentially distributed. Thus, IBTrACS suggests that TC intensification and dissipation are likely controlled by statistically random environmental and internal processes21. Meanwhile, ADT-HURSAT’s shelf-like pattern suggests that RI is a special process that could preferentially augment the probability of the highest intensification rates. The former interpretation is much more likely considering the ADT-HURSAT algorithm exhibits well-documented problems during TC eye formation (Methods).

As discussed in Bhatia et al.6, the discrepancies among basins and datasets cannot be attributed to specific physical processes or data deficiencies. A summary of the potential explanations for this behavior are included in the Methods section. We here argue that the agreement among the observational datasets in the Atlantic basin provides some level of confidence in the trend analysis in this basin. On the other hand, both observational datasets have larger uncertainties in global intensity changes. Therefore, any conclusions stemming from the global trend analysis must be treated with caution.

Figure 2 shows Atlantic basin and global trends in 24-h intensity changes between 1982 and 2009. The 28-year period is selected for analysis because it is the longest continuous sample where ADT-HURSAT produces relatively reliable global intensity data. Quantile regression is calculated for every 5% quantile between 5% and 95% (Methods). IBTrACS shows significant trends in all quantiles except the ones between 35% and 45%. The 95th percentile shows the largest trend of about +4 knots decade−1. The negative slopes for the lower percentiles and positive slopes for the upper quantiles suggests a broadening of the intensity change distribution, implying less TCs maintaining a steady intensity and more TCs exhibiting high-intensity fluctuations. ADT-HURSAT shares a similar pattern for the slopes of the quantiles but the trends are more muted. Despite being relatively small in magnitude, ADT-HURSAT quantile slopes greater than the 50% quantile are statistically significant and positive, with the 90th percentile having the largest slope of approximately 1 knot decade−1.

Fig. 2 Quantile regression of 24-h intensity changes. a–d Slope of the quantiles for 24-h intensity changes during the period 1982–2009. Slopes are shown for IBTrACS (a, c) and ADT-HURSAT (b, d) globally (a, b) and in the Atlantic basin (c, d). The black dots represent the slope derived from least squares regression of intensity change as a function of year for each quantile from 0.05 to 0.95 in steps of 0.05. Shading represents the 5th and 95th percentiles of the regressions with randomly perturbed observational data (Methods). The red solid line shows the (constant value) trend in the mean as measured by ordinary least squares regression, and the red dotted lines show the 90% confidence interval Full size image

For the global data, ADT-HURSAT trends are clearly much smaller than those in IBTrACS, yet IBTrACS trends are likely more susceptible to data homogeneity problems. Therefore, it is our interpretation that the weaker global trends shown by ADT-HURSAT are more likely to be correct. This conclusion is particularly notable because multiple studies4,13,14,15 have assumed that IBTrACS is reliable during the time frame considered in our study when in reality, there is likely a spurious trend masked by the data. It is possible that resolution and algorithm issues with ADT-HURSAT prevent it from resolving the intensity changes and have slightly dampened largest trends. However, it is much more likely that the influx of new satellites has enabled better detection of intensity changes in the later years of the time series and thus, artificially enhanced the slope of the extreme quantiles in IBTrACS.

In the Atlantic basin, the quantile regression for ADT-HURSAT and IBTrACS yields similar results for extreme intensity changes. For both observational datasets, the 95th percentile shows the largest trend. This quantile increases by over +4 knots decade−1 in ADT-HURSAT and +3 knots decade−1 in IBTrACS. The upward trend in TC activity during the 1990s is typically attributed to the Atlantic Multidecadal Oscillation (AMO), which features multidecadal fluctuations in SSTs in the Atlantic basin22. In the late 1990s, the phase of the AMO shifted from negative to positive which coincides with warming in parts of the Atlantic basin where TCs develop14. Therefore, the AMO phase favored RI at the end of the time series. In addition to the AMO behavior, the reliability of the observational data in the Atlantic basin suggests the positive trend in intensification rates during 1982–2009 is likely robust. IBTrACS has its lowest errors in the Atlantic basin because of the plentiful in situ observations that are incorporated into the best-track analysis. ADT-HURSAT relies on best-track positions to accurately assign a scene type and find an eye, which suggests that the additional reconnaissance in this basin could also help with ADT-HURSAT’s intensity estimates.

As a complement to quantile regression, we also analyze how the annual rapid intensification ratio (RI ratio) has changed from 1982 to 2009. RI ratio is defined as the number of 24-h intensity changes greater than 30 knots divided by the total number of intensity changes6, and it reflects changes in the probability of the highest TC intensification rates. Figure 3 illustrates the Atlantic basin and global trends in RI ratio for ADT-HURSAT and IBTrACS. For both observational datasets, slopes are positive and significant in the global and Atlantic basin data. Much like Fig. 2, the Atlantic basin trends for the observational datasets are very similar. Over the 28-year period, the percentage of 24-h intensity changes that exceed an intensification rate of 30 knots approximately triples in IBTrACS and ADT-HURSAT. IBTrACS shows a similar rise in RI ratio for the global data but ADT-HURSAT has a much smaller slope.

Fig. 3 Rapid intensification ratio trends. a, b Observed trends in the rapid intensification (RI) ratio of ADT-HURSAT (black) and IBTrACS (blue) over the 28-year period 1982–2009 using a global and b Atlantic data. RI ratio is defined as the number of 24-h intensity changes above 30 knots divided by the total number of 24-h intensity changes. Trends in the time series of the annual mean RI ratio are denoted by dashed lines. The slopes of the trend lines as well as their 90% confidence intervals are provided. The slopes and confidence intervals are calculated using 1000 randomly perturbed samples of the observational data. Shading represents the 5th and 95th percentiles of the 1000 regressions with these randomly perturbed observational data (Methods) Full size image

Supplementary Figure 2 shows the annual RI ratio of the two observational datasets plotted against each other. In the Atlantic basin, 35.7% of the variance in IBTrACS RI ratio is explained by ADT-HURSAT RI ratio. For the global data, the percentage explained drops to 1.8%, further highlighting the lack of agreement between the two observational datasets at global scales. The variations in the slope of RI ratio among the global observations as well as the lack of a published theory23 that explains the source of a significant positive slope in global TC intensification metrics highlight the uncertainty in the global results. Conversely, the recent uptick in intense TCs and RI magnitude in the Atlantic basin appears more plausible because of the AMO shifting phase.

Ordinary trend analysis over a relatively short period (1982–2009) may mistake multidecadal natural variability for a long-term trend. For a more refined estimate of the possible influence of natural variability on observed trends and to quantitatively assess whether the recent increase in Atlantic RI ratio can be explained by natural variability alone, we compared the observed trends to those in a TC-permitting GCM that can simulate multidecadal climate variability in the Atlantic basin. Global results are also presented with the caveat that the ADT-HURSAT data are likely more realistic because of the homogeneous time series.

Climate model simulations

HiFLOR is a high-resolution coupled GCM that can recover many aspects of the highest TC intensification rates observed in nature and capture the connection between low-frequency climate oscillations and TC behavior6,24,25,26. A recent study24 showed that HiFLOR, when nudged toward observed SSTs, can skillfully reproduce the observed year-by-year variations of the frequency of Category 4 and 5 hurricanes (maximum wind speed ≥ 113 kts) in the Atlantic basin (r ≈ 0.64) and other basins to a lesser extent. HiFLOR’s seasonal predictions of major hurricanes were also skillful in the Atlantic basin, and the high positive correlations between HiFLOR predictions and observations are still among the highest documented for a dynamical climate model25. The confirmed presence of important modes of internal climate variability within HiFLOR as well as the resolved relationship between TCs and the large-scale climate instills some confidence that HiFLOR can diagnose whether observed changes are unusual compared to expected natural variability.

For an initial exploration of whether anthropogenic climate may be contributing to the increase in TC intensification, we follow the methodology of Murakami et al.27 and examine TC variations from a suite of HiFLOR multicentury experiments with different levels of anthropogenic forcing. Four HiFLOR experiments are run using anthropogenic forcing (e.g., CO 2 , aerosols, and ozone) and natural forcing (e.g., volcanic aerosol loading and solar insolation) representative of the years 1860, 1940, 1990, and 2015 (1860CTL, 1940CTL, 1990CTL, and 2015CTL; Methods).

To compare the control simulations to observations, Supplementary Figure 3 shows the logarithm of the 24-h intensity change probabilities observed between 1982 and 1998 in IBTrACS and ADT-HURSAT and those simulated in the 1990CTL. As a technique to account for the typical error associated with measuring TC intensity, the probability of each bin was calculated using 1000 subsamples that incorporated random error (Methods). To verify the quality of the control simulations, the 17-year period between 1982 and 1998 was selected for analysis because it is the longest available period with observational data that is centered on a control simulation. The 1990CTL data matches the shape of the IBTrACS curve but the distribution is not as broad as either of the observational datasets, indicating that HiFLOR underestimates the frequency of rapid weakening and intensification. A discrepancy in the probability of the highest intensity changes was also documented by Bhatia et al.6 but the differences appear larger in Supplementary Figure 3, possibly due to the lack of SST nudging or the different years sampled.

HiFLOR underestimates the climatological rates of intensification, which will impact estimates of the amplitude of internal variability from this model. In order to mitigate the effects of systematic biases in HiFLOR intensity changes on estimating the internal variability of intensification rates, we utilize a statistical downscaling technique known as quantile mapping28. Specifically, we test two methods: quantile delta mapping29 (QDM) and bias correction quantile mapping30 (BCQM) (Methods). All subsequent analysis only incorporates the QDM algorithm because it is more adept at preserving relative changes in distributions of meteorological variables from different climate change simulations29. Figure 4 shows the effects of both quantile mapping algorithms. In this case, random error is not included in the probability estimates because we aim to nudge HiFLOR to the best guess of the intensity distribution. Clearly, the bias corrections are aligning the tails of the 1990CTL more closely with the tails of the observations. When we explore the potential influence of anthropogenic forcing on TC intensification rates, it is critical that HiFLOR provides a realistic picture of natural climate variability. Bias-corrections increase the probability of extreme intensity changes in the HiFLOR control simulations, which broadens the intensity change distribution. Thus, the bias-corrections enhance the annual variations in RI ratio and the range of the RI ratio slopes in the HiFLOR simulations, which translates into more stringent statistical tests for establishing significant observational trends.

Fig. 4 The effects of quantile mapping on 24-intensity change probability densities. a, b Common logarithm of the probability densities calculated from IBTrACS (black), ADT-HURSAT (blue), HiFLOR 1990CTL (ORIG), HiFLOR QDM-corrected 1990CTL, and HiFLOR BCQM-corrected 1990CTL 24-h intensity changes. a Global and b Atlantic basin results are plotted. QDM and BCQM are quantile mapping designed to bias correct 1990CTL to more closely represent the distribution of IBTrACS. ADT-HURSAT and IBTrACS probabilities are computed using the 17-year period centered around 1990 (1982–1998). HiFLOR probability density curves are generated using 250 years of intensity changes from the 1990CTL. Data is binned in 20 knot increments between −110 and 110 knots, and each bin entry is plotted as a dot on a curve. All distributions are bounded below by 10−5 Full size image

To test for observed RI trends that are outside of natural, internal climate variability, overlapping (not independent) 28-year RI ratio slopes in the bias-corrected HiFLOR control simulations are compared to the ADT-HURSAT and IBTrACS RI ratio slopes between 1982 and 2009. Supplementary Figure 4 shows how the probability distribution of RI ratio slopes changes when selecting 28-year segments from HiFLOR 1860CTL that are independent and not independent. The plotted bins show the independent sample has a slightly broader distribution but in general, the probabilities of the bins are similar to those in the overlapping bins.

Figure 5 contains a box and whisker plot that shows the distribution of RI ratio slopes for the QDM-corrected 1860CTL. The observed slopes during 1982–2009 are overlaid on each plot. Supplementary Figure 5 is similar to Fig. 5 but shows the box and whisker plots for all of the HiFLOR control simulations. Supplementary Figure 6 portrays the corresponding cumulative distribution functions (CDF) of the RI ratio slopes for all of the control simulations as well as the observed RI ratio slopes in ADT-HURSAT and IBTrACS. In the Atlantic basin, the slope of the RI ratio for ADT-HURSAT and IBTrACS are both above the 99th percentile of the slopes of any of the bias-corrected HiFLOR control simulations. Therefore, the large positive slope of RI ratio in both observational datasets is outside HiFLOR’s estimate of expected internal climate variability, which suggests the model’s depiction of climate oscillations like the AMO cannot explain the observed trend.

Fig. 5 Observed trends in RI ratio vs. 1860CTL natural variability. a, b Box and whisker plot represents the distribution of slopes of RI ratio in the QDM-corrected 1860 HiFLOR control simulation. a Global and b Atlantic basin results are plotted. Each slope is calculated by applying least squares regression analysis to annual RI ratio values in overlapping 28-year periods. Thus, the number of slopes for a control simulation is the number of available years subtracted by 28 (i.e., 1860CTL has 1422 slopes). The red line in each box indicates the median of the slopes. The box is bounded by the 25th and 75th percentiles of the data, and the whiskers bracket approximately 99% of the data. Red plus signs indicate outliers whose values are outside of whiskers’ range. IBTrACS and ADT-HURSAT trends in annual mean RI ratio between 1982 and 2009 are respectively represented by blue and green dotted lines and the corresponding p values are listed below each line Full size image

For the global data, IBTrACS data indicates a trend that is well above the bounds of the observed natural variability of RI ratio in any HiFLOR control simulation, but the noted temporal and spatial heterogeneities leads us to question the validity of this result. The slope of RI ratio for ADT-HURSAT is approximately the 80th percentile of the HiFLOR 1860CTL, indicating the observed trend is not significant compared to modelled climate variability. Importantly, ADT-HURSAT observed data also yields an RI slope that is still relatively unusual compared to the HiFLOR preindustrial natural variability. This emerging trend is consistent with an increasing trend in the global intensity of TCs12. In the near future, ADT-HURSAT will be extended to 2016, and it is therefore vital for future analysis to determine whether the global trend also becomes significant.

In order to explore whether the observed changes between 1982 and 2009 could be attributable to anthropogenic forcing, we compare RI ratio in the three climate change simulations to the preindustrial 1860CTL. Figure 6 shows the percent difference in RI ratio between 1860CTL and the HiFLOR simulations with stronger anthropogenic forcing (1940CTL, 1990CTL, and 2015CTL). The analysis only involves HiFLOR simulations so bias corrections and random error are not applied to the data (Methods). Grid boxes that are not statistically significant are demarcated with a white “X”. The prevalence of red boxes throughout the maps conveys that a larger percentage of TCs are undergoing RI in 1940CTL, 1990CTL, and 2015CTL than in 1860CTL. There are no significant blue grid boxes in the 1990CTL and 2015CTL maps, and the number of gridboxes that record significant increases in RI ratio is 41 (9.7% of all shaded gridboxes) in 1940CTL, 156 (34.7% of all shaded gridboxes) in 1990 CTL, and 164 (38.2% of all shaded gridboxes) in 2015CTL. Historical radiative forcing changes applied to HiFLOR clearly lead to more frequent large intensification rates, but the available simulations do not allow us to strictly attribute the identified growth in RI during the period 1982–2009 to anthropogenic forcing. An ensemble of experiments that simulate anthropogenically forced climate change during this specific period, including changes in RI ratio, would be better suited for more comprehensive attribution analysis. Thus, for now, we can only conclude that anthropogenic forcing significantly increases extreme TC intensification rates in the HiFLOR model compared to preindustrial (1860CTL) conditions.

Fig. 6 Anthropogenic forcing’s effects on RI ratio in HiFLOR. a–c Simulated changes in RI ratio by the 1940CTL (a), 1990CTL (b), and 2015CTL (c) relative to the 1860CTL. Percent difference in RI ratio between HiFLOR 1860CTL and each climate change simulation is plotted in each 5° × 5° grid box. Data is only plotted in a grid box if at least one TC passes through the grid box every 50 years in the two experiments used to calculate percent difference. Red (blue) squares indicate grid boxes where a larger (smaller) percentage of 24-h intensity changes exceed 30 knots in the climate change simulations than in the 1860CTL. Grid boxes that achieve a p value of 0.05 using a binomial proportion test are considered statistically significant. White “Xs” are located in grid boxes that are not statistically significant Full size image

In the two most reliable long-term observational records available for TC intensity changes, the proportion of 24-h TC intensification rates greater than 30 knots significantly increases in the Atlantic basin between 1982 and 2009. Results averaged over all basins show a significant increase in TC intensification rates in IBTrACS but not in ADT-HURSAT. By itself, a 28-year upward trend in a TC intensification metric does not necessarily reflect the effects of anthropogenic climate forcing because of the intrinsic natural variability of the climate system. Here, we use bias-corrected TC intensification rates simulated by HiFLOR to demonstrate that natural variability cannot explain the magnitude of the observed upward trend in the Atlantic basin. These conclusions are possible because HiFLOR is a unique climate model that can successfully simulate the most intense TCs and highest intensification rates in multicentury simulations.

This study is limited by the ability of a climate model to accurately represent natural variability as well as the uncertainty around the trends in relatively short observational records. However, this study represents a crucial first step in quantifying the precise roles of stochastic processes, anthropogenic warming, and natural variability when assessing changes in TC intensification rates. Further analysis with additional high-resolution climate models and a longer and more reliable observational record is required to confirm these conclusions. Regardless, these trends provide another example of the potentially serious repercussions of anthropogenic warming for TCs.