Anthropogenic forcing changes coastal tropical cyclone frequency

Abstract

It remains a mystery if and how anthropogenic climate change has altered the global tropical cyclone (TC) activities, mainly due to short reliable TC observations and climate internal variabilities. Here we show with large-ensemble TC-permitting simulations that the observed increases in TC frequency since 1980 near the US Atlantic coast and Hawaii are likely related to the aerosol and greenhouse gases (GHG) effects, respectively. The observed decrease in the South China Sea after 1980 could be associated with GHG emissions alone, whereas the observed increase near Japan and Korea during this period would be related to the aerosol and GHG combined effects. These changes in coastal TC frequency are explained by the responses of large-scale environmental conditions to anthropogenic forcing.

Introduction

Anthropogenic activities have been affecting our climate since the industrial revolution¹. The effect of anthropogenic climate change on tropical cyclones (TCs) is of particular concern because of their catastrophic damage at landfall and relatively frequent occurrence compared to other natural hazards. However, it is notoriously challenging to attribute a regionally observed tropical cyclone frequency (TCF) change to anthropogenic climate change. The reasons for this difficulty are twofold. First, climate internal variability can considerably modulate regional TC activities^2,3,4,5. Second, reliable TC observations are only available since the satellite era, and therefore, not long enough to separate the potential effect of anthropogenic climate change from climate internal variability.

With the rapid development of high-performance computing, large-ensemble climate simulations have become a powerful tool to overcome this obstacle. Since the simulated climate internal variability in models can be, by design, out of phase among ensemble members through the initial condition large-ensemble simulations, the model-estimated effect of anthropogenic climate change on TCF can emerge after averaging out the internal variability by taking the ensemble mean⁶. Following this path, it was reported recently that anthropogenic climate change might have played an important role in the observed global TCF change since 1980, including the substantial increase in the North Atlantic and decrease in the western North Pacific^6,7. Other multiple attempts^{8,9,10,11,12,13} have also been made by scientists on the detection and attribution of the anthropogenic impact on TCF, but mostly from a basin-mean perspective. However, considering the actual TC threat posed to our society, the coastal TCF change under anthropogenic forcing is a central and perhaps ultimate concern.

Based solely on TC observations, it has recently been found that TC activities have been migrating towards coastal regions since 1980, related to a westward shift of TCF in most of the global basins¹⁴. However, it is still unknown whether these observed coastal changes are related to anthropogenic climate change or climate internal variability. In this study, we will explore this vital question by investigating the effect of anthropogenic forcing on the coastal TCF change globally with large-ensemble TC-permitting simulations (see Methods). This is an important step forward to understanding better the TC threat as the substantial impact of human activity on our climate system continues now and will do in the future.

Results

SVD analysis

This study’s TCF is defined as the annual total number of TC records observed in each $5^\circ \times 5^\circ$ grid covering the globe. The target period of TCF change is confined from 1980 to 2020 when the global TC observation is more reliable than in the pre-satellite era¹⁵. Figure 1a shows the observed epochal TCF difference (2001–2020 minus 1980–2000). For coastal regions, we have been seeing significantly increased TCF along the US Atlantic coast, around Hawaii, near the Northeast Asian coast (i.e., the Japan-Korea region), along the western coast of the Arabian sea, and to the north of Madagascar. Meanwhile, significant decreases in the observed TCF have happened near Mexico/Central America west coast, in the South China Sea, and in proximity to the Australian coast. With the listed areas above, we first try to identify the coastal regions where the TCF change is less likely to be related to climate internal variability with the singular value decomposition (SVD) analysis.

Fig. 1: Observed TCF change for 1980–2020 and regional SVD modes.

a TCF epochal change after and before the year 2000. The white “x” denotes statistically significant change at a 90% confidence interval based on a two-sided bootstrap test. b–g Regional SVD analysis using SST and TCF. b Composite of the SVD mode representing climate internal variability characterized by the IPO. c As in (b), but for the SVD mode representing the global mean SST change. d Time series of the expansion coefficient (EC) for SST of the natural variability mode superimposed on the observed IPO index with a flipped sign. e As in (d), but for the global mean SST mode. f The squared covariance fraction (SCF) map of the SVD mode representing natural variability (only the grids with SCF ≥ 50% are colored). g As in (f), but for the global mean SST mode. The magenta boxes in (a) highlight six coastal regions where the SCF of global warming mode is no less than 50% in (g) and the observed TCF change in (a) is statistically significant.

Two kinds of SVD analysis are performed (global and regional SVD analyses, see Methods). In the global SVD analysis (Supplementary Fig. 1), the spatial pattern of the first SVD sea surface temperature (SST) shows high similarity to the SST configuration of the Interdecadal Pacific Oscillation^16,17 (IPO). Supplementary Table 1 confirms that the expansion coefficient for the first SST mode shows the highest correlation with the IPO index (r² = 0.89, p-value < 0.05) among all four natural variability indices examined. The second SVD mode represents the global mean SST change (Supplementary Fig. 1d–f), showing a high covariance between the expansion coefficient and the long-term global mean SST increase (r² = 0.87, p-value < 0.05). Consistent global SVD results can be reproduced with the TC records in the North Atlantic, eastern North pacific and western North Pacific with an extended period of 1960–2020 (Supplementary Fig. 2 and Supplementary Table 2) so that the results of global SVD analysis appear to be independent on the length of analyzed period.

Based on the physical representation of the first and second modes in the global SVD analysis, we next conduct the regional SVD analysis to identify the coastal regions where TCF is highly correlated with the global mean SST mode but less with the natural variability mode (hereafter referring explicitly to the IPO). A regional SVD on one map grid utilizes the same global SST field but different regional TCF information on the central and eight surrounding grids. The regional SVD is repeated until the global TCF field is covered (see Methods). Figure 1b–e shows the composite of SVD SST and expansion coefficient in the modes that are highly correlated with the natural variability and global mean SST change, respectively. On average, the natural variability mode appears as the first mode with a global mean SST mode as the second mode, which is consistent with the global SVD outcomes.

Figure 1f, g shows the regional co-variability of SST and TCF in the natural variability and global mean SST modes, respectively, measured by the squared covariance fraction (SCF) with all the regional SVD results combined. Here we define a group in the variation of regional TCF at one grid is less likely dominated by natural variability when the SCF score of the natural variability mode is less than 50% in Fig. 1f, or equivalently, more than 50% in Fig. 1g representing the grids dominated by the mode of global mean SST change. Superposing the colored regions in Fig. 1g and observed significant change of TCF in Fig. 1a, we can identify several coastal regions where the observed TCF change in the past four decades are more likely to be related to the global mean SST warming, indicating the potential influence of anthropogenic forcing on coastal TCF changes. These coastal regions include the US Atlantic coast, Hawaii region, Northeast Asian coast near Japan and Korea, South China Sea, western coast of the Arabian Sea, and Madagascar region (all highlighted with magenta boxes in Fig. 1a). Although SVD analyses do not guarantee perfect separation between natural variability and anthropogenic modes, we assume that the global mean SST warming mode based on observations (Fig. 1) is potentially related to anthropogenic climate change. Next, we will leverage large-ensemble simulations to explore further if the effects of anthropogenic climate change on simulated TCF changes over these coastal regions are also predominant over the same period.

Large ensemble simulations

The effect of anthropogenic climate change can be broken down into anthropogenic aerosol effect (e.g., related to sulfate, organic carbon and black carbon emissions), and anthropogenic GHG effect (e.g., due to carbon dioxide, methane, and other GHG emissions). Four suites of large-ensemble simulations are performed with both aerosol and GHG forcing (named as AllForc), only GHG forcing (named as AllForc_NoAE), only natural forcing (named as NatForc, including, e.g., volcanic eruption, solar insolation variability), and pre-industrial forcing representative of 1850 (named as CNTL). Supplementary Fig. 3 shows that the large-ensemble simulations generate a consistent climatology of annual global TC counts as in observation, with a median and one standard deviation of 85 and 9 TCs, respectively. By applying the regional SVD analysis to all the ensemble members of AllForc, we find the potentially dominant role of anthropogenic forcing on coastal TCF changes around the US Atlantic coast, Hawaii region, Northeast Asian coast near Japan and Korea, South China Sea, western coast of the Arabian Sea, and Madagascar region (Supplementary Fig. 4g), which showed significant trends in TCF by the observations (Fig. 1g). We next identify the single forcing effects of anthropogenic aerosol, GHG and natural forcing on the potential coastal TCF changes, by calculating the ensemble mean TCF trend difference in pairs of the four simulation suites (see Methods for details).

Compared to the observed coastal TCF changes in Fig. 1a, the TCF trend difference in Fig. 2 shows consistent changes—with statistical significance at a 10% level—along the US Atlantic coast and around Madagascar under the aerosol effect (Fig. 2a), around Hawaii and in the South China Sea under the GHG effect (Fig. 2b), and near the Japan-Korea region under a combined effect of aerosol and GHG (Fig. 2d). These coastal regions are where we speculate that the observed TCF change is less likely to be related to the natural variability but more likely related to at least part of anthropogenic climate change. The observed significant increase in TCF to the western coast of the Arabian Sea was not simulated under any forcing in the model. Natural forcing may be important to the observed increase of TCF around Madagascar but does not show any consistent sign of changes with the observed TCF changes over the other targeted coastal regions.

Fig. 2: Differences in simulated ensemble-mean TCF trends for 1980–2020.

The trends are calculated based on annual TCF numbers, due to (a) the anthropogenic aerosol effect (AllForc minus AllForc_NoAE), b The anthropogenic greenhouse gas effect (AllForc_NoAE minus NatForc), (c) natural forcing (NatForc minus CNTL) and (d) the combined effect of anthropogenic aerosol and GHG effects (AllForc minus NatForc). The white dot shows the grid where the ensemble mean TCF trends in two experiment suites are statistically different at a 90% confidence interval based on a two-sided t-test approach. The magenta boxes are the same as in Fig. 1a.

TCF origen analysis

TCF changes could be related to TC track and genesis changes. We next utilize a previously developed TCF origen analysis⁹ (see Methods) to further attribute the differences in TCF trends to those in TC genesis and/or TC tracks. Along the US Atlantic coast, the positive TCF trend under the aerosol effect (Fig. 3a) is primarily contributed by the TC genesis increase in the TC main development region (Fig. 3b, highlighted with the green box). Consistent with a previous study⁹, the positive TCF trend around Hawaii under the GHG effect (Fig. 3e) is also largely related to the TC genesis increase to the south of Hawaii (Fig. 3f). Meanwhile, TC track changes over the eastern North Pacific in the deep tropics (Fig. 3g) also favor the TCF increase near Hawaii (Fig. 3e). For the negative TCF trend in the South China Sea under the GHG effect (Fig. 3i), both TC genesis reduction and unfavorable TC track change to the east of the Philippines are related (Fig. 3j, k). For the TCF increase around the Japan-Korea region (Fig. 3m), the TC track changes to the south of this region are found to be the most important (Fig. 3o), followed by the increase in TC genesis (Fig. 3n), under the combined effect of aerosol and GHG emissions. Figure 3k, o suggest a general track direction shift for the TCs over the Philippine Sea from entering the South China Sea to moving toward the Japan-Korea region, which agrees with previous findings focusing on the future projection^18,19,20. For the positive TCF trend around the northern tip of Madagascar under the aerosol effect, both the local TC genesis increase and TC track change are important (Fig. 3r, s).

Fig. 3: TCF origen analysis applied to the five coastal regions.

a–d The US Atlantic coast under the anthropogenic aerosol effect. a Local TCF trend change under the aerosol effect as shown in Fig. 2. b–d The contribution of each grid to the total TCF change ($\varDelta {df}$) in the magenta box due to genesis change ($\varDelta d{g}_{B,{A}_{o}}$), TC track change ($\varDelta d{t}_{B,{A}_{o}}$), and nonlinear combination of genesis and track changes ($\varDelta d{n}_{B,{A}_{o}}$). e–h As in (a–d), but for the Hawaii region under the anthropogenic GHG effect. (i–l) As in (a–d), but for the South China Sea under the anthropogenic GHG effect. m–p As in (a–d), but for the Japan-Korea region under the anthropogenic aerosol and GHG effects. q–t As in (a–d), but for the Madagascar region under the anthropogenic aerosol effect. The magenta boxes show the same regions as in Fig. 1a. The green boxes highlight the key regions with the genesis and track contributions.

DGPI Analysis

We next explain the causes for the TC genesis and track changes in the key regions highlighted by the green boxes in Fig. 3, by analyzing two environment-based proxies—the genesis potential index (GPI) and steering flow. We chose a recently developed dynamic GPI^21,22 (DGPI) to represent the TC genesis potential. This choice was made mainly considering a higher consistency of the simulated TC genesis frequency (Supplementary Fig. 5) with the DGPI (Fig. 4a–c, spatial correlation r² = 0.47) than that with a conventional choice—the Emanuel-Nolan GPI²³ (Supplementary Fig. 6, spatial correlation r² = 0.34), particularly in the Southern Hemisphere.

Fig. 4: Trend difference in TC dynamic genesis potential index (DGPI) and steering flows.

The ensemble-mean trends are based on the mean change in each year of 1980–2020 during the main TC seasons only. Northern Hemisphere: July-October; Southern Hemisphere: December (of the previous year)-March. DGPI and 300–850 hPa mass-weighted steering flows are calculated with monthly-mean data before taking the main TC season mean. The white dot in (a–c) shows the statistically significant change in DGPI. The dark black arrow in (d–f) denotes significant changes in zonal and/or meridional steering components. The white and gray backgrounds in (d–f) show easterly and westerly zonal wind anomalies, respectively. A statistical significance test is conducted with a 90% confidence interval based on a two-sided t-test approach. The green boxes highlight the key regions with the genesis and track contributions.

The simulated DGPI trend difference shows evident and significant increases in the Atlantic TC main development region and around the northern tip of Madagascar under the aerosol effect (Fig. 4a), an increase to the south of Hawaii and a decrease to the east of the Philippines under the GHG effect (Fig. 4b), and an increase to the south of Japan-Korea region under the combined effect of aerosol and GHG emissions (Fig. 4c). It is encouraging to see that the DGPI trends in these key regions (green boxes in Fig. 4a–c) agree well with the genesis contribution analyzed in Fig. 3. This allows us to further break down the DGPI index to investigate which environmental factor contributes the most to the DGPI trend change with different anthropogenic forcing (Table 1, Supplementary Fig. 7 and see Methods). Around Madagascar, although the vertical motion shows the highest contribution of 119% to the DGPI increase (Table 1), we decided not to investigate this region further in this study, considering the DGPI trend difference there in Fig. 4a is highly localized with limited statistical significance.

Table 1 Decomposition of DGPI trend change by region.

For the DGPI increase in the Atlantic TC main development region, Table 1 shows that the vertical wind shear change contributes the most (by 40%) under the aerosol effect, mainly due to a zonal shear reduction (Supplementary Fig. 8a–c). This reduced shear is dominated by a significant easterly anomaly at 200 hPa (Supplementary Fig. 8d), which indicates the weakening of subtropical jet with the anthropogenic aerosol forcing.

For the DGPI increase to the south of Hawaii, Table 1 also shows that the vertical wind shear change, dominated by the variation in zonal winds (Supplementary Fig. 9a–c), is the largest contributor (by 56%). Supplementary Fig. 10d shows that under the GHG effect, an SST warming tongue extends from midlatitudes to the central Pacific. Such a central Pacific warming can perturb the climatological Walker circulation, weakens 200-hPa westerly (Supplementary Fig. 9d) and 850 hPa easterly (Supplementary Fig. 9e) simultaneously, and therefore, reduces the vertical wind shear and increases DGPI to the south of Hawaii.

For the DGPI decrease to the east of the Philippines, the vertical motion weakening in the mid-troposphere appears as the most important factor (by 70%) under the GHG effect (Table 1). This suppressed vertical motion, located over the Pacific warm pool, connects to an easterly anomaly at 200 hPa (Supplementary Fig. 9d) and a westerly anomaly at 850 hPa (Supplementary Fig. 9e) to the west of the dateline, and an anomalously upward motion above the equatorial central Pacific (Supplementary Fig. 10b) where the SST is anomalously warmed (Supplementary Fig. 10d). Thus, the central Pacific warming under the GHG effect may be a primary driver for the DGPI reduction to the east of the Philippines by weakening the Walker circulation.

For the DGPI increase to the south of Japan, enhanced meridional shear vorticity at 500 hPa is a dominant factor under a combined effect of aerosol and GHG emissions (55%, Table 1). Figures S11a–c show an easterly anomaly around 30^oN in the western North Pacific under aerosol and/or GHG effect. Strong easterly anomalies cover the northern part of the key genesis region to the south of Japan (green box in Supplementary Fig. 11c), producing a positive change in meridional shear vorticity and causing an increase in the DGPI there.

What is the physical mechanism for the aforementioned 500 hPa easterly anomaly at 30^oN? Further analysis reveals that the easterly anomaly at 500 hPa is related to a weakening of the subtropical jet at 200 hPa (Supplementary Figs. 8d, 9d). We can also see an enhancement of 200 hPa zonal wind to the north of 50^oN (Supplementary Figs. 8d and 9d), indicating a poleward shift of the subtropical jet over the western North Pacific and Northeast Asia. As a response due to thermal wind balance, the subtropical jet may be shifted poleward if the meridional temperature gradient is weakened²⁴. Indeed, Supplementary Fig. 12a–c confirm the mid-latitude warming in the lower troposphere under aerosol and GHG effects. The mid-latitude warming decreases the meridional temperature gradient, weakens the subtropical jet via a thermal wind adjustment, and therefore, generates 500-hPa easterly anomalies at 30^oN as shown in Supplementary Fig. 11.

Steering flow analysis

For the steering change in the western North Pacific, a westward steering anomaly around 30^oN (Fig. 4d–f) is found to be tied to the weakened subtropical jet aloft. This westward steering anomaly, covering the northern part of the box in Fig. 4f, lets TCs stay longer in the Japan-Korea region under a combined effect of aerosol and GHG emissions. In the southern part of the box in Fig. 4f, we can see a cyclonic circulation over the Philippine Sea, which should be mainly due to the GHG effect (Fig. 4e). This cyclonic circulation generates eastward TC steering to the east of the Philippines (Fig. 4e), creating an unfavorable steering condition that reduces the probability for TCs entering the South China Sea under the GHG effect. Meanwhile, it can also increase the chance for TCs to move meridionally into higher latitudes where westward steering anomaly can push TCs toward the Japan-Korea region under the combined effect of aerosol and GHG emissions (Fig. 4f).

What is the physical driver behind the cyclonic circulation over the Philippine Sea in Fig. 4e? We show that under the GHG effect, there is an equatorial central Pacific warming (Supplementary Fig. 10d), accompanied by a central Pacific ascending anomaly (Supplementary Fig. 10b), and 200 hPa easterly anomaly (Supplementary Fig. 9d) and 850 hPa westerly anomaly (Supplementary Fig. 9e) between the Philippines and dateline. Together with the cyclonic circulation over the Philippine Sea, all these features dramatically resemble a classic Gill-type circulation pattern²⁵. We, therefore, speculate that the cyclonic circulation over the Philippine Sea is very likely to be a Gill-type response to the central Pacific warming under the GHG effect.

The Gill-type circulation pattern may also be related to the steering change affecting the TCF near Hawaii. Figure 4e shows that under the GHG effect, there is a weak westward steering enhancement in the eastern part of the Hawaii-related box—favoring TCs moving toward the central North Pacific and Hawaii. The enhanced westward steering is part of the equatorial easterly anomaly in the East Pacific, which vanishes in the central Pacific where the SST shows anomalous warming (Supplementary Fig. 10d). This equatorial easterly anomaly is also consistent with the zonal-vertical circulation of the Gill-type response to the central Pacific warming due to GHG emissions.

Discussion

The regionally oceanic and atmospheric warming may be the primary cause of the coastal TCF change under anthropogenic forcing. We find that two kinds of warming are related.

First, the anomalous lower-tropospheric warming at midlatitudes under the aerosol effect weakens the subtropical jet over the North Atlantic, which reduces the vertical wind shear over the Atlantic Main development region, increases the cyclogenesis probability, and thus increase the TC activity along the US Atlantic coast. Over Northeast Asia and the western North Pacific, aerosol and GHG effects can also generate lower-tropospheric warming at midlatitudes, which weakens the subtropical jet, therefore increasing the chance for TCs approaching the Japan-Korea region.

It was reported that the observed anomalous tropospheric warming at midlatitudes for 1979–2005 reduces the meridional temperature gradient, providing evidence that the subtropical jets have been shifting poleward²⁴. The causality between the midlatitude warming in the lower troposphere and the poleward shift of the subtropical jet was also examined and established in another simulation study²⁶. Since the jet stream can be regarded as the Hadley cell terminus, the subtropical jet shifts poleward in tandem with a widening of the Hadley cell²⁷. Increases in GHG concentration²⁸ can expand the tropical belt, which is in line with our analysis showing a poleward shift of subtropical jet under the GHG effect. Still, more research inputs are warranted to identify the key mechanisms linking the anthropogenic forcing and the observed tropical expansion in the preceding decades²⁷.

Second, the equatorial central Pacific warming under the GHG effect stimulates a classical Gill-type circulation pattern in the northern tropical Pacific, which could be related to the increase in TCF near Hawaii and reduction in the South China Sea by varying the local TC steering and genesis conditions. This equatorial central Pacific warming also appears in the SST pattern of the global mean SST mode in our global SVD analysis (Supplementary Fig. 1d), in line with the SVD result in a previous study⁶. In the North Pacific, the SST warming pattern under the GHG effect (Supplementary Fig. 10d) shows substantial similarity to the warming configuration during an El Niño Modoki II event²⁹, the intensity and relative frequency of which have been increasing since the late twentieth century³⁰. A similar SST warming trend in the central Pacific can also be seen in a reconstructed SST dataset³¹ for 1854-2019. More analyses are still needed to confirm further if the equatorial central Pacific warming has been happening since 1980. However, our analysis, at least, indicates that this potentially warming pattern may be associated with anthropogenic GHG emissions.

It should be noted that there is only one climate model used in this study and the results may be model dependent. Despite this caveat, this study shows that the observed TCF change in several coastal regions over the past 40 years is potentially associated with anthropogenic climate change, mainly through regional warming and its remote influence on circulations that changes the TC genesis and steering conditions. Many studies have focused on the global- and basin-mean change of TC activity^6,8,32. However, considering the real threat posed by TCs, the TCF changes in coastal regions are perhaps our ultimate concern. From this perspective, our study reveals that the inhabited coastal regions worldwide may face further uneven changes, with increases in landfalling TC threats in some regions and decreases in others, as anthropogenic climate change continues.

Methods

Data

The TC best track data was taken from the International Best Track Archive for Climate Stewardship³³ (IBTrACS, Version 4.0). Six ocean basins are considered, including the North Atlantic (NA), eastern North Pacific (EP), western North Pacific (WP), North Indian Ocean (NI), South Indian Ocean (SI), and South Pacific (SP). The borders of the basins are shown in Fig. 1a. The best track data for the NA and EP were provided by the National Hurricane Center, and the Joint Typhoon Warning Center’s best track was used for the WP, NI, SI, and SP. We utilized the best track records at 00, 06, 12, and 18 Coordinated Universal Time with a maximum sustained near-surface wind speed of at least 34 kts (tropical-storm strength) observed over oceans within the 60^oN-60^oS latitude band.

The monthly sea surface temperature (SST) data were taken from the Hadley Centre Sea Ice and Sea Surface Temperature data set³⁴ (HadISST). The ERA5 monthly reanalysis data from the European Centre for Medium-Range Weather Forecasts³⁵ were also used for historical climatology analysis.

In this study, all the TC frequency-related analysis is based on the annually accumulated TC count (i.e., in Figs. 1–3, S1–S5, S13–S15, S17). We only consider environmental condition changes during the main TC seasons, that is, July-October for the Northern Hemisphere, and December (of the previous year)-March in the Southern Hemisphere (i.e., in Fig. 4 and S6-S12). Note that the NI shows different TC seasonality and we do not explain the TCF change in the NI with the environmental condition changes presented in this study.

Before further analysis, TC best tracks and environmental variables were linearly interpolated onto a 5^o-latitude by 5^o-longitude map. The results shown in this study are not sensitive to the choice of grid resolution for interpolation.

The IPO index is calculated as the standardized second principal component of the empirical orthogonal function (EOF) analysis for the 13-year low-pass-filtered global SST¹⁶. This calculation is applied to both observations and simulations.

Models

SPEAR³⁶ was used in this study, which comprises 50 km atmospheric and land mesh and 100 km mesh for the ocean and sea-ice components. Consistent with the pre-processing of the observed best tracks, model-generated TCs were identified and tracked every 6 h using a previously developed TC tracking scheme³⁷ with a maximum wind speed of at least 16.5 m s⁻¹, well-established warm core, and minimum duration of 36 h³⁸. The SPEAR model can well simulate the seasonality of TC number in all the basins compared to the observation (Supplementary Figs. 13, 14), and reasonably reproduce the IPO-like variability⁶ and its relation to TCF (Supplementary Fig. 15).

Large-ensemble experiments

We used SPEAR to perform initial condition large-ensemble experiments. Five kinds of simulations were performed, that is, CNTL, AllForc, AllForc_NoAE, NatForc, and AllForc_Nudge_SST, which are explained in detail as follows.

CNTL

A 3000-year pre-industrial climate simulation was generated. We then resampled thirty times for 41 consecutive years (same length as from 1980 to 2020) within the last 2000 years of the simulation to form 30 ensemble members. These 30 members represent the climate scenario with pre-industrial forcing.

AllForc

The historical external forcing was prescribed for 1921–2014, followed by the Shared Socioeconomic Pathway 5-85 (SSP5-85)³⁹ after 2015. The prescribed external forcing includes both anthropogenic and natural forcing. Anthropogenic forcing includes greenhouse gases, anthropogenic aerosols, and ozone, whereas natural forcing includes solar constant, volcanic eruptions, and natural aerosols such as dust. Volcanic forcing was only prescribed before 2006. There are 30 ensemble members.

AllForc_NoAE

The external forcing is identical to the AllForc experiment, but with the emissions of anthropogenic aerosols, such as sulfates, organic carbon, and black carbon, fixed at the level of 1921. There are 12 ensemble members. Supplementary Fig. 16 shows that these 12 members are enough to largely remove the natural internal variability after taking the ensemble mean.

NatForc

The natural forcing is identical to the AllForc experiment, but with the fixed anthropogenic forcing at the level of 1921. There are 30 ensemble members.

AllForc_Nudge_SST

As in AllForc, but with SSTs nudged to HadISST at 5-day time scale for 1971–2020. There is only 1 ensemble member. This simulation is used as part of the validation of SPEAR TC climatology against the observation for 1980–2020 (Supplementary Fig. 3).

Separation of anthropogenic effects on TC activity

We consider four kinds of forcing on TC activity, including the climate internal variability (e.g., the IPO), natural forcing (e.g., volcanic eruption), anthropogenic GHG forcing (e.g., CO₂ emission), and anthropogenic aerosol forcing (e.g., sulphate emission).

The initial condition of the ensemble members is based on a 2000-year “preindustrial” control run with atmospheric composition fixed at levels representative of calendar year 1850. The ensemble simulations were started from the atmospheric, ice and ocean restart files from the control simulation at years 101, 121, and every 20 years thereafter combined with the 1921 land restart file (due to limitations imposed by land vegetation transitions) of the first historical member of the ensemble. In this way, the fully coupled SPEAR can simulate largely uncorrelated natural internal variability among ensemble members (e.g., as shown in Supplementary Fig. 16). By taking the ensemble mean, climate internal variabilities can be largely averaged out in SPEAR since they are out of phase among ensemble members by design⁶. In this way, the effect of anthropogenic climate change and natural forcing on TC activity can be separated. For example, the first global SVD mode in the ensemble mean of SPEAR simulation represents a clear global warming mode (Supplementary Fig. 17).

We consider four kinds of climate change effects on TC activity under the following forcing:

• Natural forcing: the ensemble-mean difference between experiments NatForc minus CNTL,

• Anthropogenic GHG forcing: the ensemble-mean difference between experiments AllForc_NoAE minus NatForc,

• Anthropogenic aerosol forcing: the ensemble-mean difference between experiments AllForc minus AllForc_NoAE, and

• Anthropogenic forcing (combined effect of GHG and aerosol): the ensemble-mean difference between experiments AllForc minus NatForc.

To explore the effect of each forcing on TCF, we calculated the difference of the ensemble-mean TCF trend for 1980–2020 between two suites of experiments in the related pair.

SVD analysis

SVD analysis is a multivariate statistical method that can identify the synchronized spatial and temporal correlations between two 3-dimensional variables, which, in this study, refer to as the SST and TCF matrixes. The SVD outputs are multiple paired SST and TCF patterns, each of which corresponds to an orthonormal mode. For each pair of the SST and TCF patterns (i.e., each mode), a pair of expansion coefficient (EC) time series were also generated for SST and TCF, respectively. A squared covariance fraction (SCF, in a unit of %) was calculated for each mode, describing the fraction of co-variability between the SST and TCF patterns in that mode relative to the total co-variability.

Two types of SVD analyses were performed in this study. For the first kind (hereafter, the global SVD), the global SST and TCF fields were utilized as inputs (Supplementary Figs. 1 and 2), which was also applied in a previous study⁶ by Murakami et al. Some consistencies between regional TC activity and SST pattern can be observed. For example, Supplementary Fig. 1a, b show a large reduction of TCF in the open ocean of the western North Pacific with weak increase in TCF along the East Asia coast, which is in line with a westward shift of storm activity under a La Niña type SST pattern².

For the second kind (hereafter, the regional SVD), the global SST field and regional TCF information were employed as inputs (Fig. 1b–g). Specifically, for one grid point ${A}_{i,j}$, nine grids of ${A}_{i\pm 1,j\pm 1}$ were used as the TCF input for the regional SVD analysis at ${A}_{i,j}$. In this way, each regional SVD analysis can generate 9 modes. Most of the grids yielded the first two leading modes as global SST mode and IPO mode. By calculating the Pearson correlation coefficient of the global mean SST and IPO time series, respectively, with the EC time series of each mode, we could identify two modes that best represent the global mean SST change and IPO, respectively. The regional SVD was repeated to cover all the global grids. Two modes selected from each regional SVD analysis were used for the composite results shown in Fig. 1b–g.

Origin analysis for TCF trend

The TCF origen analysis⁹ is a method that can identify the locations of the contribution of TC genesis, TC track change, and the non-linear combination of the two to the mean TCF trend change in a region of interest.

The TCF at an individual grid $A$ in one year can be written as

$$f={\iint }_{C}g\times t\,d{A}_{o},$$

(1)

where $f=f\left(A\right)$ is the TCF at grid A in this year, $g=g\left({A}_{o}\right)$ is the frequency of TC genesis in a remote grid cell ${A}_{o}$, $t=t\left(A,{A}_{o}\right)$ is the probability that a TC generated in the grid ${A}_{o}$ can propagate to grid $A$, and $C$ is the region over which the integration is calculated (e.g., over one ocean basin).

Now we decompose $f$ in one year with the climatological mean ($\bar{f}$) and the anomaly (${f}^{\mbox{'}}$) in this year, i.e., $f=\bar{f}+{f}^{\mbox{'}}$, and Eq. (1) can be rewritten as

$$\bar{f}+{f}^{{\prime} }={\iint }_{C}\left(\bar{g}+{g}^{{\prime} }\right)\times \left(\bar{t}+{t}^{{\prime} }\right)\,d{A}_{o}.$$

(2)

By removing the climatological mean relationship ($\bar{f}={\iint }_{C\,}\bar{g}\times \bar{t}{d}{A}_{o}$), Eq. (2) can be simplified as

$${f}^{{\prime} }={\iint }_{C}{g}^{{\rm{\text{'}}}}\times \bar{t}\,d{A}_{o}+{\iint }_{C}\bar{g}\times {t}^{{\prime} }\,d{A}_{o}+{\iint }_{C}{g}^{{\prime} }\times {t}^{{\prime} }\,d{A}_{o}.$$

(3)

The map of annual TCF trend, denoted as $\frac{\partial }{\partial t}\left(\bar{f}+{f}^{\mbox{'}}\right)=\frac{\partial }{\partial t}\left({f}^{\mbox{'}}\right)$, can be then decomposed with Eq. (3) as

$$\frac{\partial }{\partial t}\left({f}^{{\prime} }\right)=\frac{\partial }{\partial t}\left({\iint }_{C}{g}^{{\prime} }\times \bar{t}\,d{A}_{o}\right)+\frac{\partial }{\partial t}\left({\iint }_{C}\bar{g}\times {t}^{{\prime} }\,d{A}_{o}\right)+\frac{\partial }{\partial t}\left({\iint }_{C}{g}^{{\prime} }\times {t}^{{\prime} }\,d{A}_{o}\right).$$

(4)

To establish the origen analysis for TCF trends, we rewrite Eq. (4) as

$$\underbrace{\frac{\partial }{\partial t}({f^{\prime} }_{B,{A}_{o}})}_{d{f}_{B,{A}_{o}}}=\underbrace{\frac{\partial }{\partial t}({\iint }_{B}{g^{\prime} }_{{A}_{o}}\times {\bar{t}}_{{A}_{o}}\,dA)}_{d{t}_{B,{A}_{o}}}+\underbrace{\frac{\partial }{\partial t}({\iint }_{B}{\bar{g}}_{{A}_{o}}\times {t^{\prime} }_{{A}_{o}}\,dA)}_{d{g}_{B,{A}_{o}}}+\underbrace{\frac{\partial }{\partial t}({\iint }_{B}{g^{\prime} }_{{A}_{o}}\times {t^{\prime} }_{{A}_{o}}\,dA)}_{d{n}_{B,{A}_{o}}},$$

(5)

where $A$ presents all the grids in the region of interest, $B$ (e.g., the US Atlantic coast). The three terms on the r.h.s. shows the effects—i.e., genesis change, track change and nonlinear combination of the two, respectively—of a remote grid cell ${A}_{o}$ on the mean TCF change in the region of $B$. The term on the l.h.s. represents the overall effect of the changes in grid ${A}_{o}$ on the TCF mean change in the region of $B$.

Next, we apply Eq. (5) to an experiment 1 that serves as a reference (e.g., AllForc_NoAE), which can be written as

$$\frac{\partial }{\partial t}\left({{f}^{{\prime} }}_{B,{A}_{o}(1)}\right)=\frac{\partial }{\partial t}\left({\iint }_{B}{{g}^{{\prime} }}_{{A}_{o}(1)}\times {\bar{t}}_{{A}_{o}(1)}\,{dA}\right)+\frac{\partial }{\partial t}\left({\iint }_{B}{\bar{g}}_{{A}_{o}(1)}\times {{t}^{{\prime} }}_{{A}_{o}(1)}\,{dA}\right)+\frac{\partial }{\partial t}\left({\iint }_{B}{{g}^{{\prime} }}_{{A}_{o}(1)}\times {{t}^{{\prime} }}_{{A}_{o}(1)}\,{dA}\right).$$

(6)

Similarly, we can also apply Eq. (5) to an experiment 2 that can be regarded as a sensitivity test (e.g., AllForc), which can be written as

$$\frac{\partial }{\partial t}\left({{f}^{{\prime} }}_{B,{A}_{o}(2)}\right)=\frac{\partial }{\partial t}\left({\iint }_{B}{{g}^{{\prime} }}_{{A}_{o}(2)}\times \left({\bar{t}}_{{A}_{o}(1)}+{\bar{t}}_{{A}_{o}(\varDelta )}\right)\,{dA}\right)+\frac{\partial }{\partial t}\left({\iint }_{B}\left({\bar{g}}_{{A}_{o}(1)}+{\bar{g}}_{{A}_{o}(\varDelta )}\right)\times {{t}^{{\prime} }}_{{A}_{o}(2)}\,{dA}\right)+\frac{\partial }{\partial t}\left({\iint }_{B}{{g}^{{\prime} }}_{{A}_{o}(2)}\times {{t}^{{\prime} }}_{{A}_{o}(2)}\,{dA}\right),$$

(7)

where ${\bar{t}}_{{A}_{o}(\varDelta )}={\bar{t}}_{{A}_{o}(2)}-{\bar{t}}_{{A}_{o}(1)}$ and ${\bar{g}}_{{A}_{o}(\varDelta )}={\bar{g}}_{{A}_{o}(2)}-{\bar{g}}_{{A}_{o}(1)}$, which represent the difference of climatological mean of track probability and genesis frequency, respectively, in experiment 2 relative to experiment 1.

By taking the difference of Eqs. (7) and (6) we have

$$\frac{\partial }{\partial t}\left({{f}^{{\prime} }}_{B,{A}_{o}\left(2\right)}\right)-\frac{\partial }{\partial t}\left({{f}^{{\prime} }}_{B,{A}_{o}\left(1\right)}\right)=\frac{\partial }{\partial t}\left({\iint }_{B}{{g}^{{\prime} }}_{{A}_{o}\left(2\right)}\times {\bar{t}}_{{A}_{o}\left(1\right)}\,{dA}-{\iint }_{B}{{g}^{{\prime} }}_{{A}_{o}\left(1\right)}\times {\bar{t}}_{{A}_{o}\left(1\right)}\,{dA}\right)+\frac{\partial }{\partial t}\left({\iint }_{B}{\bar{g}}_{{A}_{o}\left(1\right)}\times {{t}^{{\prime} }}_{{A}_{o}\left(2\right)}\,{dA}-{\iint }_{B}{\bar{g}}_{{A}_{o}\left(1\right)}\times {{t}^{{\prime} }}_{{A}_{o}\left(1\right)}\,{dA}\right)+\frac{\partial }{\partial t}\left({\iint }_{B}{{g}^{{\prime} }}_{{A}_{o}\left(2\right)}\times {{t}^{{\prime} }}_{{A}_{o}\left(2\right)}\,{dA}-{\iint }_{B}{{g}^{{\prime} }}_{{A}_{o}\left(1\right)}\times {{t}^{{\prime} }}_{{A}_{o}\left(1\right)}\,{dA}+{\iint }_{B}{{g}^{{\prime} }}_{{A}_{o}\left(2\right)}\times {\bar{t}}_{{A}_{o}\left(\varDelta \right)}\,{dA}+{\iint }_{B}{\bar{g}}_{{A}_{o}\left(\varDelta \right)}\times {{t}^{{\prime} }}_{{A}_{o}\left(2\right)}\,{dA}\right),$$

(8)

where the term on the l.h.s., denoted now as $d{f}_{B,{A}_{o}}$, shows the total effect of a remote grid ${A}_{o}$ on the TCF trend difference between the two experiments in the region of $B$. The three terms on the r.h.s., denoted as $d{f}_{B,{A}_{o}}$, $d{g}_{B,{A}_{o}}$ and $d{n}_{B,{A}_{o}}$, respectively, show the breakdown effects of genesis, track, and nonlinear combination of the two from a remote grid ${A}_{o}$ on the mean TCF change in the region of $B$ (e.g., the anthropogenic aerosol effect on the mean TCF change along the US Atlantic coast in Fig. 3b–d). Compared to the individual influence of TC genesis and track, the nonlinear combination of the two (the last column of Fig. 3) shows a much weaker impact on the simulated TCF trend difference in all five regions.

Genesis Potential Index (GPI)

GPI is an index that describes the likelihood to have a TC genesis given the local environmental conditions. The dynamic GPI (DGPI) can be defined as²²

$$DGPI=\mathop{\underbrace{{(2.0+0.1\times {V}_{s})}^{-1.7}}}\limits_{V-SHEAR}\mathop{\underbrace{{(5.5-\frac{d{u}_{500}}{dy}\times {10}^{5})}^{2.3}}}\limits_{M-SHEAR}\mathop{\underbrace{{(5.0-20\times {\omega }_{500})}^{3.4}}}\limits_{OMEGA}\mathop{\underbrace{{(5.5+|{\zeta }_{a850}\times {10}^{5}|)}^{2.4}}}\limits_{VORT}{e}^{-11.8}-1.0,$$

(9)

where ${V}_{s}$ is the vertical wind shear magnitude (m s^-1) between 200 and 850 hPa, ${u}_{500}$ is the zonal wind (m s^-1) at 500 hPa, $\left(-\frac{d{u}_{500}}{{dy}}\right)$ denotes the meridional shear vorticity (s^-1) at 500 hPa, ${\omega }_{500}$ is the vertical velocity in a pressure coordinate (Pa s^-1), ${\zeta }_{a850}$ is the absolute vorticity (s^-1) at 850 hPa, and $e$ is the natural base. The calculated DGPI value was set back to zero between 5^oN and 5^oS or where the relative SST becomes negative. These two criteria were justified against observations in previous studies^21,22. Note that the DGPI includes thermodynamic factors indirectly via the OMEGA term. ω₅₀₀ highly correlates with the 600-hPa relative humidity and local SST anomaly relative to the tropical mean²¹. These high correlations reflect a physical linkage among the three: locally high SST may foster 500-hPa ascent-related moisture convergence that tends to moisten the lower troposphere by increasing the relative humidity at 600 hPa²¹.

The definition of DGPI, i.e., Eq. (9), can be linearized by taking the logarithm as

$$\log ({DGPI}+1)=\log (V-{SHEAR})+\log (M-{SHEAR})+\log ({OMEGA})+\log ({VORT})-11.8.$$

(10)

Equation (10) is used to calculate the relative contribution of different terms to the mean DGPI change in a region of interest (Table 1 and Supplementary Fig. 7).

For comparison, the conventional GPI, i.e., Emanuel-Nolan GPI (ENGPI) is also calculated, which can be defined as²³

$${ENGPI}={\left|{10}^{5}\eta \right|}^{3/2}{\left(\frac{{RH}}{50}\right)}^{3}{\left(\frac{{V}_{{pot}}}{70}\right)}^{3}{\left(1+0.1{V}_{s}\right)}^{-2},$$

(11)

here ${V}_{{pot}}$ is the maximum potential intensity (m s^-1) defined by the local thermodynamic profile⁴⁰.

Statistical significance test

The statistical significance of this study is defined with 90% confidence intervals. We calculated the statistical significance between two linear trends with a t-test (Figs. 2, 4, S5, and S8-S12). For the epochal TCF change (Fig. 1a), a bootstrap approach¹⁴ was utilized. Two tested populations were resampled in pairs for 10,000 times. The difference in the means of each pair was then calculated to form a new distribution with 10,000 samples, from which the 90% confidence intervals were obtained.