TC Observations
We used the International Best Track Archive for Climate Stewardship (IBTrACS18), v04r00, and the Advance Dvorak Technique-Hurricane Satellite-B1 (ADT-HURSAT17,19) for the period 1982–2017. The ADT-HURSAT data record was recently expanded by 8 years to span the 39-year period between 1979-2017. We omit the first three years of the record where limited geostationary data results in missing storms and focus on the 36-year period between 1982-2017. For our IBTrACS analysis, we only consider best-track data from the National Hurricane Center for the Atlantic and east Pacific and the Joint Typhoon Warning Center for the remainder of the globe18. One of the benefits of only using data from these U.S. agencies is they follow the same definition of maximum winds: the highest 1-min average at 10 m height over a smooth surface37. Best-track data start as operational estimates of the intensity and track of a TC and are refined at the end of a TC season with a combination of in situ (e.g., dropsondes, scatterometers, buoys), radar, and satellite measurements. Best-track intensity and position estimates are available every six hours at the four synoptic times (0000, 0600, 1200, and 1800 UTC) and are recorded to the nearest 5 knots (1 kt = 0.5144 m s−1) and 0.1° latitude/longitude38.
The creation of ADT-HURSAT consists of three main steps. Geostationary satellite imagery is first analyzed from International Satellite Cloud Climatology Project (ISCCP)-B1 data39-41. Then, the data is centered on IBTrACS TCs and subsampled to be both spatially and temporally homogeneous. Finally, a simplified version of the advanced Dvorak technique42 is used to evaluate the data and determine a maximum TC wind speed. ADT-HURSAT data are produced every three hours based on satellite data that has been uniformly subsampled to a horizontal resolution of 8 km, and wind speeds are recorded to the nearest tenth of a Dvorak “T-number” (depending on the current intensity, between 1-3 knots).
Defining Storm-Local and Tropical-Mean Environments
Data for the analysis of cyclone environments were taken from ERA527 and MERRA-228,43. Reanalyses combine a forecast model with continuous assimilation of observational data to construct a global representation of historical atmospheric variability. The native spectral resolution of ERA5 is T639 (nominally 31 km) and the native resolution of MERRA-2 is 0.5º ´ 0.625º (nominally 50 km). We focus on these reanalyses because they capture realistic moisture, temperature, and wind values in the lower atmosphere44,45. For storm-local environment calculations, we identified and tracked tropical cyclones in ERA5 and MERRA-2 objectively using the Lagrangian feature-tracking algorithm—TRACK—of Hodges31. This methodology was documented in greater detail by Hodges et al.46.
Vertically averaged relative humidity (at 850, 700 and 600 hPa), vertical wind shear (computed as the square-root of the sum of the squared differences in the zonal and meridional winds between 850 and 200 hPa), and sea-surface temperature were spectrally filtered to T11 resolution to remove cyclonic circulation features and retain only the large-scale, background environmental fields. Sensitivity analysis of ERA5, the highest-resolution reanalysis available, shows that 95 % of the T639-resolution cyclonic circulation is removed at a T11 truncation (Supplementary Fig. 3). Along-track sampling of mean values from the spectrally filtered fields within a 5º storm-centered radius (geodesic) was performed. To compute potential intensity along cyclone tracks, reanalysis data (sea-surface temperature, mean sea-level pressure, air temperature, and specific humidity) were regridded first to 1º resolution. Regridding from the native reanalysis resolutions to 1º has a negligible effect on potential intensity and no spectral filtering was performed. Potential intensity was computed by taking input fields from the grid cell nearest to the storm center at each timestep, following Bister and Emanuel47 and using published code48,49 (Gilford et al., 2017; 2019). Vertical soundings of temperature and specific humidity were constructed from reanalysis data on 14 isobaric levels (1000, 925, 850, 700, 600, 500, 400, 300, 250, 200, 150, 100, 70, and 50 hPa). By default, the code allows reversible ascent and dissipative heating; the ratio of the exchange coefficients of enthalpy and momentum flux is 0.9; and output velocity is scaled by 0.8 to reflect surface drag. Further discussion of these constants is given in Emanuel50 and Gilford et al.48,49 (2017; 2019).
The tropical-mean of these fields are also calculated but with no spectral filtering. Annual averages are comprised of peak-TC-season means in each hemisphere, August-October (ASO) in the Northern Hemisphere and February-April (FMA) in the Southern Hemisphere and then averaged between the two Hemispheres (area-weighted). June-October (JJASO) in the northern hemisphere and December-February (DJFMA) were also tested but yielded comparable results. In the northern hemisphere, the tropical-mean is an area-average of environmental values over the ocean between 10°-30°N and 40°E-20°W. In the southern hemisphere, the bounds for averaging are 10°-30°S and 30°E-150°W.
HiFLOR experiments
HiFLOR control simulations introduced in Murakami et al.51 and Bhatia et al.6 were used here to represent natural (unforced) climate variability and provide the framework for exploring anthropogenic effects on TC intensification. We focused on the control simulation that used anthropogenic forcing fixed at 1860 (1860CTL) levels because it has the longest simulation length: 1500 years. The first 50 years of the simulation were disregarded to mitigate effects of model drift. Additional details on data preparation and methodology for the HiFLOR simulation can be found in Bhatia et al.6. The approach developed by Harris et al.52 is used to track TCs in HiFLOR and is applied using the parameter values of Zhang et al.53 and Murakami et al.15.
CMIP6 Experiments
We examine linear trends of the four tropical-mean environmental fields from all available CMIP6 simulations over the period 1982-2014. The fields are defined in an identical way to the observed fields in MERRA-2 and ERA5. Models and the number of ensemble members for each model that contain the relevant environmental variables are listed in Supplementary Table 3.
Storm Criteria
For consistency, intensity change values in HiFLOR, and the observational datasets are rounded to the nearest five knots. We only consider TCs that are active for at least 72 hours and exceed wind speeds of 34 knots for at least 36 hours. We restrict our analysis sample to only consider cases where the TC center is located over the ocean, the starting and ending TC position are below 40 degrees of latitude, and the TC intensity stays above 34 knots. The warm core criteria discussed in Murakami et al.15 are also applied to the HiFLOR data before analysis. For the analysis involving storm-local environments, all storm fixes within 0.5 degrees of land in any direction are removed to reduce the number of spurious PI readings.
Uncertainty Quantification
For Fig. 1, we use Monte Carlo techniques to create random noise before analyzing the discretized data. Random noise prevents multiple data points from having the same value and provides an estimate of the typical error associated with measuring TC intensity. 1000 subsamples were produced by adding random noise from a uniform distribution on the interval ± knots to each intensity change value. The magnitude of this random noise is derived by adding 5 knots of error in quadrature (propagation of errors stemming from the intensity change calculation), which is a conservative estimate for the typical error associated with each TC intensity observation. Fig. 1 involved the calculation of the 5th and 95th percentile in each of the 1000 subsamples for each year. 1000 slopes of each percentile were calculated, and the mean of the slopes was considered the best estimate of the 1982-2017 slope of the quantile. The 5th and 95th percentiles of the 1000 slopes were shaded as the uncertainty bounds.