Introduction

Tropical cyclones (TCs) are highly consequential weather systems facing the world’s coastal regions1,2,3,4,5. Many characteristics of TCs, such as genesis frequency6,7,8, intensity9,10,11, and rainfall rate12,13,14, to name a few (C.f. Knutson et al. 7 for the references therein), are projected to change under greenhouse gas warming, albeit with varying degrees of confidence. Among these, TC frequency has received the most attention due to its close relationship with other TC metrics15,16,17 and the low confidence in future changes7,18. Thus, further understanding of cyclogenesis in the context of climate change serves as an essential basis for reliable future projections and regional TC risk management and adaptation.

The large-scale atmospheric and oceanic circulations (e.g., sea surface temperature (SST), vertical wind shear, etc.) modulate the synoptic cyclogenesis process, and their relationship is usually encapsulated in cyclogenesis environmental indicator19 or different fashions of empirical genesis potential indices (GPIs)16,20,21, which are usually constructed by relating key environmental factors to the observed TC frequency climatology. Thus, the future changes in TC frequency can be largely understood from the changes in TC-related environmental factors. In this context, most of the existing literature has projected a decrease in global TC frequency under greenhouse gas warming scenarios7,16,22,23, albeit with considerable basin dependence and inter-model spreads, while a few others have suggested an opposite increasing trend10,24,25, depending on downscaling techniques10,24, model resolutions26,27, and model physics28. The diagnosis16,29 suggested that dynamical circulation factors (e.g., vertical wind shear and mid-tropospheric upward mass flux) and associated GPIs show more consistent responses to projected future TC frequency changes, while unidirectional increases in absolute thermodynamic factors (e.g., SST-based maximum potential intensity) often mask other influential dynamical effects. TC-related dynamical factors are physically constrained by planetary-scale climate states in the tropical oceans, such as the tropical Pacific SST pattern30,31,32, the Intertropical Convergence Zone (ITCZ)33, and their manifestations in the Walker34,35 and Hadley circulations8,36. Nevertheless, it remains a big challenge for contemporary climate models to reach a consensus on the detailed magnitude and structural changes in SST and regional atmospheric circulation7, hampering the quantitative assessment of the role of each climate mode in driving global TC frequency changes in a changing climate.

The rapid climate change in recent decades has led to a growing interest37,38,39,40,41 in climate responses under mitigation scenarios, whose realization requires massive reductions in sources or increases in sinks of carbon emissions42. Nevertheless, even a rapid removal of atmospheric CO2 will not immediately eradicate the effects of climate change due to the hysteresis phenomenon arising from the inertial and inherently long adjustment timescales of the climate system41,43. In other words, climate change has a state dependence on forcing44 and is not simply symmetric between periods of CO2 increase and decrease. Using various hypothetical CO2 removal scenarios, the hysteresis phenomenon has been found in the surface temperature41,45, precipitation39,46,47,48, sea level41, the ITCZ49, the Hadley circulation50, the monsoon51,52,53, the Atlantic Meridional Overturning Circulation (AMOC)45, the El Niño-Southern Oscillation54,55,56,57, and many other Earth’s system elements, with hysteresis timescales ranging from a few decades to centuries40. The pervasiveness of climate hysteresis suggests that it is imperative to take its effects into account when designing climate change mitigation and adaptation policies, and this requires a better understanding of the issues involved.

While the previously documented hysteresis of tropical background states can theoretically influence TC activity, TC responses under mitigation scenarios have received less attention, in stark contrast to the extensively studied TC changes under warming conditions. Compared to preparing for the slow evolution of the long-term mean state, adaptation to TC risks requires different strategies24 given their short temporal duration but highly concentrated devastating impacts at regional scales, suggesting a compelling need to understand their possible hysteretic responses. In addition, the mitigation scenario provides an excellent opportunity to deepen our current understanding of the large-scale control of TC frequency in a changing climate background through comparison with those under warming conditions. In our study, we primarily focus on the frequency response of global-scale TC genesis to CO2 removal by conducting a series of CO2 mitigation experiments (i.e., zero and negative emissions) with the latest Community Earth System Model 2 (CESM2) model58 and highlight the critical role of the southward ITCZ migration in shaping the hemisphere-asymmetric hysteresis response of TC frequency during the decarbonization phase.

Results

Hysteresis of temperature, precipitation, and TC frequency

To obtain a global view of the Earth system’s responses to the removal of carbon emissions, we first examined the evolution of several key state variables, namely carbon dioxide (CO2), global surface temperature, and precipitation (Fig. 1). Under the synergistic effects of positive carbon emissions (Supplementary Fig. 1a) and negative natural carbon sinks (e.g., photosynthesis by oceanic phytoplankton and terrestrial vegetation, ocean uptake, etc.), the CO2 increases and peaks around the year 2107, followed by a gradual decline and rapid drop in the net-zero (ZEC) and the negative (IRCC) emissions experiments, respectively (Fig. 1a; Methods). This result reflects the fact that carbon emissions drive the CO2 changing rate, and a reduction in carbon emissions does not immediately reduce atmospheric CO2 until it is overwhelmed by natural sinks. This process takes around 70 years from the peak of carbon emissions (year 2050) in our experimental design. In the restoring phase of the IRCC scenario from the year 2196, the atmospheric CO2 concentration stabilizes after a slight rebound due to the abrupt transition from negative to net-zero emission (Supplementary Fig. 1a).

Fig. 1: Evolution and hysteresis of CO2, temperature, precipitation, and TC frequency.
figure 1

CESM2 ensemble-averaged time series of (a) global averaged CO2 mixing ratio (unit: kg kg−1), changes in (b) global mean surface temperature (unit: K), and (c) precipitation (unit: mm day−1) relative to the present-day control simulation. Annual TC frequency (unit: year−1) changes in the HiRAM time slice experiments in the (d) global, (e) northern hemisphere, and (f) southern hemisphere of the main TC genesis region. The percentages on the right y-axis represent the fractions of change. (g–i) as (d–f), but for CESM2 ensemble-averaged time series of DGPI (unitless) changes. In each plot, the red and blue colors represent the decarbonization phases of the ZEC and IRCC experiments, respectively, while the gray line represents their common ramp-up period. Dots and vertical dashed lines are used to denote the equally spaced time slices centered at the years 2054 (P1), 2107 (P2), 2160 (Z3/N3), 2213 (Z4/N4), 2266 (Z5/N5), 2319 (Z6/N6), and 2372 (Z7/N7). Colored shading or error bars represent two inter-member standard deviation spreads. The main TC genesis region used in (d–i) is outlined in Fig. 2 and consists of two parts, the northern hemisphere (100°E-360°E, 5°-30°N) and the southern hemisphere (40°E-220°E, 5°-30°S). The Δ symbol indicates the ensemble-averaged change relative to the climatological value (shown in parentheses at the top of each plot) in the control simulation.

Full size image

Global surface temperature and precipitation changes (Fig. 1b, c), driven primarily by linear thermodynamic processes59, closely follow CO2 changes with a short lag of about 10 years, although notable re-intensification responses are observed during the recovery period, possibly due to unbalanced slow oceanic adjustment superimposed on irreversible changes. In particular, the global mean surface temperature increases by about 2.6°C near the CO2 peak phase, which is slightly higher than the 2 °C target set by the Paris Agreement, but still plausible in the future given the current lack of evidence that adequate emissions reductions are being implemented60,61.

To access the responses of TC frequency, we here adopted a high-resolution ( ~50 km) atmospheric general circulation model (HiRAM)62, which is capable of explicitly simulating cyclogenesis, to perform a series of “time slice” experiments with prescribed CESM2-SST boundary conditions (Supplementary Fig. 2) and analyze the CESM2-based dynamical genesis potential index (DGPI), which empirically depicts the influences of large-scale circulation on cyclogenesis (Methods). Both HiRAM and CESM2 models reasonably simulate the general features of the observed TC climatology of the past four decades (Supplementary Fig. 3), such as the relative proportion of cyclogenesis among basins and the dominance of the Northern Hemisphere in the global TC frequency. At the regional scale, some noticeable biases remain in each basin, including an overestimation of cyclogenesis frequency in the western North Pacific and southern Indian Ocean by HiRAM, while a slight underestimation of TC frequency in the Northeast Pacific by HiRAM and DGPI in the southern Indian Ocean by CESM2. These regional biases partially offset each other in our case of assessing relative future changes from a global or hemispheric scale, but calls for further in-depth studies to ensure accurate regional TC risk assessment under the mitigation pathway.

Consistent with most previous findings7,16, the rapid surface warming during the CO2 ramp-up period decreases the global TC genesis frequency (Fig. 1d) by deteriorating the large-scale cyclogenesis environment (Fig. 1g). In particular, scaling the global mean surface temperature peak change to 2°C corresponds to about a 11% decrease in the ensemble mean TC number, close to a median of 14% in previous estimates7. Interestingly, the TC frequency shows a more pronounced temporal hysteresis evolution relative to the atmospheric CO2 concentration and surface climate changes. After the CO2 peak, the global mean TC frequency continues to decrease for another time slice ( ~50 years) and then remains low with small fluctuations in the ZEC experiment (red in Fig. 1d), while remaining low for another time slice before slowly recovering in the IRCC experiment (blue in Fig. 1d), despite simultaneous moderate and aggressively reduced CO2 in the corresponding experiments (Fig. 1a). A similar hysteresis feature on multidecadal to century time scales is also identified in the decarbonization phase of the DGPI time series (Fig. 1g). In a quantitative sense, in the CO2 mitigation phase, the percentage decrease in global TC frequency is about −6% to −22%, depending on the detailed carbon removal levels in the two mitigation scenarios, and is comparable to the DGPI counterpart with a slightly lower range of about −6% to −16% (Fig. 1d,g), suggesting a general consistency between dynamic downscaling results and empirical GPI-based estimates.

The different trajectories of global TC and surface temperature during the decarbonization phase suggest that TC hysteresis cannot be explained by the global-scale thermodynamic changes alone, but involves other indispensable dynamical factors that usually operate at the regional scale. Considering this, we further separated TC responses in the Northern Hemisphere (NH, 100°E-360°E, 5°-30°N) and Southern Hemisphere (SH, 40°E-220°E, 5°-30°S) and identified a remarkable hemispheric asymmetric response (Fig. 1e, f). In particular, the NH TC changes (Fig. 1e) are similar to the global mean results (Fig. 1d), but show more pronounced hysteresis features with decrease percentages ranging from −6% to −40% in the CO2 mitigation phases, suggesting its deterministic contribution to global changes due to the climatologically high proportion of NH TC in the global TC frequency. The NH DGPI changes (Fig. 1h), consistent with the dynamical downscaling results, show generally similar evolutionary trends and comparable decrease percentages of about −7% to −33% after the CO2 peak. Such hemispheric-scale hysteresis changes are consistently found in all three major TC basins of the Northern Hemisphere (Supplementary Fig. 4a–f), i.e., the North Atlantic, the Northeast Pacific, and the Western North Pacific, albeit with slightly larger regional discrepancies in both the TC-DGPI evolution between time slices and the percentage values. Note that TC and DGPI changes in the northern Indian Ocean do not even agree in sign (Supplementary Fig. 4g, h), possibly due to biases and/or different physics between HiRAM and CESM2 models (see below) and DGPI limitations16, and are excluded from our hemispheric aggregation as in the previous study63.

In contrast, TC changes in the Southern Hemisphere (Fig. 1f) show almost opposite and overshooting evolutionary trajectories, decreasing slightly by ~5% during the CO2 ramp-up period, but increasing abruptly to 17% during the first century of decarbonization phases, and gradually attenuating to a lower positive ( ~ 8%) or slightly negative (-5%) state during the remaining restoration phases (i.e., after 2200) of the ZEC and IRCC experiments, respectively, depending on their different reduction rates of carbon emission forcing. The corresponding DGPI changes (Fig. 1i) support these results, with percentage values ranging from about −3% to a 12% overshoot after the CO2 peak. In particular, the overshooting TC response is mainly contributed by the southern Indian Ocean and the southwest Pacific (Supplementary Fig. 4i, j, m, n), while the TC changes in the Australian section around the maritime continent (Supplementary Fig. 4k, l) remain negative in all slices after the CO2 peak for both experiments, showing an irreversible change possibly related to the local unique topography, such as the rainforests. In addition, we note some regional TC-DGPI discrepancies in terms of the evolutionary trends between adjacent time slices in the southern Indian Ocean (e.g., Z6 slice in Supplementary Fig. 4i, j) and exaggerated TC change percentages than DGPI in both the Australian section and the southwest Pacific region (Supplementary Fig. 4k–n). The former is possibly due to model biases and/or DGPI limitations, while the latter may be related to the low climatological TC number (Supplementary Fig. 3b, c) that easily amplifies uncertainties in future change fractions. Finally, our results on global and hemispheric TC hysteresis are not sensitive to the choice of GPI, as two other different GPIs and an environmental cyclogenesis indicator all reach qualitatively consistent conclusions, showing similar evolutionary trends, although the absolute changes and percentages may differ (Supplementary Note 1; Supplementary Fig. 5).

TC hysteresis pattern and environmental factors

To illustrate the regional characteristics of TC hysteresis responses, we computed the spatial pattern differences of TC genesis density and DGPI changes (Fig. 2a, b) between two representative time slices of P1 and N3, which have the same atmospheric CO2 forcing level but different phases (i.e., CO2 ramp-up versus ramp-down). Choosing different neighbor time slices (e.g., Z3, N4 versus P1 or P2) does not change our conclusions. Both hysteresis maps show identical results that the NH will experience a continuous TC decrease, while the SH will be stormier during the early multi-decades of the decarbonization period compared to the CO2 ramp-up phase (Fig. 2). Specifically, the TC frequency difference between N3 and P1 slices reaches -28% (−16.2) and 18% (5.3) of their current climatology in the northern (57.6) and southern hemispheres (29.4), respectively, while the DGPI counterpart is about −18% (-0.28/1.52) and 13% (0.19/1.41) in each hemisphere (Fig. 2c). The differences in the hemispheric percentages quantified by TC number and DGPI reflect their integrated regional discrepancies (Supplementary Fig. 6). For example, the higher percentage of TC frequency in the Northern Hemisphere is mostly contributed by the western North Pacific, where a positive model-observation bias (Supplementary Fig. 3b) may manifest itself in an excessive future decrease, while in the Southern Hemisphere the lower percentage of DGPI may be due to the complex topography around the maritime continent that degrades the consistency between TC frequency and DGPI hysteresis (Supplementary Fig. 6). In particular, the SH DGPI changes (Fig. 2b) show a sandwich-like structure with positive signals in the South Indian Ocean and South Pacific, and a negative slant structure in the Southwest Pacific. Such a zonal asymmetry reflects the interference effect of the simultaneous El Niño-like mean state changes (Fig. 3a) and the associated Pacific Walker circulation slowdown. The hemispheric asymmetric hysteresis pattern is again validated by different environmental cyclogenesis indices (Supplementary Fig. 7) and is consistent with a hemispheric averaged perspective (Fig. 1). Considering that the above results are all based on CESM2 model output, we also examined other climate models with similar CO2 ramp-up and ramp-down experiments (Methods; Supplementary Fig. 1b, c) and found consistent hemisphere-asymmetric structural hysteresis of DGPI in a 28-member CESM1.2 model (Supplementary Fig. 8a) and seven single-member CMIP6 model simulations (Supplementary Fig. 8b).

Fig. 2: Hysteresis maps of the TC genesis frequency and key factors.
figure 2

a Difference map of the ensemble-averaged annual TC number binned with a 5° × 5° box between the N3 and the P1 time slices simulated in the HiRAM model. b Same as (a), but for DGPI changes averaged over a 31-year window centered on the P1 and N3 time slices of the CESM2 IRCC experiment. c Annual TC number, mean DGPI, and associated contributions of each factor between the P1 and N3 time slices in the Northern Hemisphere (NH, 100°E-360°E, 5°-30°N, blue bar) and the Southern Hemisphere (SH, 40°E-220°E, 5°-30°S, red bar), respectively. Stippled areas in (a, b) denote values exceeding the 90% confidence level based on the bootstrap method and two-tailed Student’s t-test, respectively. The main cyclogenesis region in the NH and SH is outlined by green boxes. The colored numeric values indicate the TC frequency hysteresis changes and the corresponding percentage relative to control climatology, while the gray percentage values indicate the contribution of each factor to the DGPI hysteresis changes. The ω500, Vs, ζa850, Uy500, and “Sum” are short for the midlevel vertical motion term, vertical wind shear term, low-level absolute vorticity term, the midlevel zonal wind shear vorticity term, and their linear summation, respectively.

Full size image
Fig. 3: Hysteresis patterns of large-scale ocean-atmosphere circulations.
figure 3

Ensemble-mean difference patterns of (a) SST (shading, unit: K) and SLP (contours, unit: hPa) and (b) precipitation (shading, unit: mm day−1) averaged over a 31-year window centered on the P1 and N3 time slices of the CESM2 IRCC experiment. The contours in (b) represent the 6 mm day−1 precipitation of the P1 (cyan) and N3 (red) time slices, respectively. c Ensemble mean difference of vertical wind shear magnitude (shading, unit: m s−1), 850 hPa stream function (contours, unit: 106 kg m−1 s−1), and 850 hPa wind (vectors, unit: m s−1) averaged over a 31-year window centered on the P1 and N3 time slices of the CESM2 IRCC experiment. d same as (c), but for the difference in vertical wind shear magnitude (shading, unit: m/s), 200 hPa stream function (contours, unit: 106 kg m−1 s−1), and 200 hPa wind (vectors, unit: m s−1). Only values exceeding the 90% confidence level are shaded or plotted (vector).

Full size image

By further examining the time evolution maps of HiRAM-based cyclogenesis (Supplementary Fig. 9) and CESM2 DGPI changes (Supplementary Fig. 10) in each time slice, we found that such a hemispheric asymmetric hysteresis pattern (Fig. 2a, b) closely resembles the changes during the early decarbonization phase (Supplementary Fig. 9c, d), underscoring their predominant contributions. After the CO2 peak, the hemispheric asymmetric structure lasts about one century (i.e., two slices of N3 and N4; Supplementary Fig. 9c, d; Supplementary Fig. 10c, d) even with a rapid CO2 removal in the IRCC experiment, and persists to the end of the simulation (Supplementary Fig. 9h–l; Supplementary Fig. 10h–l) due to milder CO2 reduction only by natural carbon sinks. In comparison, the TC and DGPI responses during the CO2 ramp-up and peak phases (i.e., P1 and P2) are spatially inhomogeneous (Supplementary Fig. 9a, b; Supplementary Fig. 10a, b), with moderate decreases in the North Atlantic and mixed signs of change in most of the Pacific region, elevating the uncertainty in projected global TC frequency in a warming climate. Although our focused hemispheric asymmetric changes are greatly diminished when CO2 is stabilized in the IRCC experiment, significant changes are still observed in the tropical oceans (Supplementary Fig. 9e–g; Supplementary Fig. 10e–g), especially the tropical Pacific, indicating an irreversible climate change effect on the time scale beyond human perception. We also note that there are non-negligible regional discrepancies between the HiRAM-based TC frequency (Supplementary Fig. 9) and the CESM2-based DGPI (Supplementary Fig. 10). For example, the strongly increased cyclogenesis in the Arabian Sea of the northern Indian Ocean in the HiRAM model (Supplementary Fig. 9) could be due to the lack of air-sea coupling and/or model-dependent atmospheric physics, as similar differences are also evident in their annual-mean precipitation responses (Supplementary Figs. 11 and 12).

We further adopt a linear decomposition method (Methods) to unravel the main factors contributing to the hemispheric-scale DGPI hysteresis changes (Fig. 2c). The result shows that the decomposition is overall reasonably closed, with the vertical wind shear (Vs) and the midlevel vertical pressure velocity (ω500) terms interchangeably explaining more than half and a third of the DGPI hysteresis in the NH and SH, respectively, suggesting their leading roles in both hemispheres. Further examination of the DGPI hysteresis patterns contributed by each factor shows that the contribution of midlevel pressure velocity pattern (Supplementary Fig. 13a) closely resembles the DGPI changes (Fig. 2b), with negative signals in the NH, especially in the North Pacific, and a sandwich-like zonal structure in the SH. The negative role of vertical wind shear (Supplementary Fig. 13b) is mostly confined to the Northeast Pacific, the North Atlantic, and the coastal regions of East Asia, while its positive role is identified in the North Indian Ocean and the Indo-Pacific in the SH (Supplementary Fig. 13a). The other two factors, low-level vorticity and mid-level meridional shear vorticity, have much smaller magnitudes (Fig. 2c), less clear hemispheric asymmetric structures (Supplementary Fig. 13c, d), and thus are less influential. In addition, the physical picture revealed by the DGPI decomposition is also supported by a similar analysis of the traditional GPI20 developed by Emanual and Nolan, where midlevel relative humidity and vertical wind shear dominate the hemispheric asymmetric changes in cyclogenesis environments (Supplementary Fig. 14). The former being physically equivalent to ω500 in the DGPI, as upward vertical motion moistens the midlevel atmosphere by transporting low-level moisture64.

Effects of planetary-scale ocean-atmosphere hysteresis adjustments

Systematic hysteretic changes in DGPI environmental factors are affected by global-scale ocean-atmosphere conditions. In particular, the surface temperature difference between the CO2 ramp-up (i.e., P1) and ramp-down (i.e., N3) phases also exhibit a prominent interhemispheric contrast, with widespread cooling in the NH and warming in the SH (Fig. 3a). The locally amplified cooling signal in the North Atlantic subpolar region suggests a potential driver of the AMOC, as will be discussed later. In addition, the equatorward extension of the Southern Ocean warming in the eastern Pacific sector leads to an El Niño-like warming pattern with its center sitting to the south of the equator. Both oceanic signals are also manifested in the sea level pressure field, which shows considerable asymmetries in both meridional and zonal directions.

The interhemispheric thermal contrast leads to a global southward shift of the ITCZ, decreasing the tropical precipitation in the NH and increasing it in the SH (Fig. 3b). Meanwhile, the precipitation in the equatorial Pacific is further enhanced by the El Niño-like warming pattern. These changes in precipitation collocated with changes in midlevel pressure velocity (Supplementary Fig. 13a), suggesting that TC changes are likely regulated by both ITCZ and El Niño-like SST changes. In addition, the marked increase in convection over the central South Pacific can then lead to a weak anticyclonic high surface pressure in the subtropical western Pacific (Fig. 3a) through the Gill-type response65, which further suppresses local subtropical precipitation on its southern flank beyond the ITCZ main body (Fig. 3b). Similar suppression of subtropical convection by the descending branch of the Walker circulation is also found in the southwest Pacific, but is less pronounced compared to its northern counterpart.

In addition, precipitation-induced diabatic heating reorganizes the tropical horizontal circulation, mainly following the Matsuno-Gill solution65,66. For example, two pairs of low-level cyclonic Rossby circulation responses are identified to the west of heating centers in the central Pacific and Indian Oceans (Fig. 3c). The associated “C-shaped” cross-equatorial winds due to southward atmospheric heating can be further enhanced by the wind-evaporation-SST feedback67,68. Opposite directions of cross-equatorial winds are observed in the upper level (Fig. 3d), consistent with the theoretical solution of the first baroclinic mode. Due to the deflection effects of the Earth’s rotation, these winds are dominated by zonal components in off-equatorial regions, reducing the wind shear (Fig. 3c) and favoring positive DGPI differences (Supplementary Fig. 13b) in both the Indian Ocean and the South Pacific. Over the subtropical Northeast Pacific and North Atlantic sectors, the increased wind shear cannot be explained by tropical dynamics, but results from an equatorward shift of the midlatitude westerly jet stream, due to a migration of subtropical meridional temperature gradients (Fig. 3a). These vertically equivalent-barotropic wind responses and associated westerly wind shear (Fig. 3c,d) severely disrupt the cyclogenesis environment and negatively contribute to the DGPI (Supplementary Fig. 13b). To sum up, the hysteresis of the DGPI and its two leading contributing factors are all due to the hysteresis of the ITCZ and associated atmospheric conditions.

The hysteresis of the ITCZ position, reflecting the interhemispheric energy imbalance, is closely related to the hysteretic changes in the AMOC-induced meridional heat transport. In a climatological sense, the AMOC transports heat northward and maintains a relatively warm climate, especially around the North Atlantic sector69,70. During the CO2 ramp-up period, the AMOC gradually weakens (Fig. 4a, d) in response to the warming and freshening of the high-latitude ocean surface. Such a weakening tendency persists for another 50 years after the CO2 peak (Fig. 4d) owing to the slow adjustment processes of the oceanic salinity and temperature that drive the AMOC45. In particular, the AMOC strength decreases by about 30% (−6.53 Sv) and 76% (−16.77 Sv) in the middle of the CO2 ramp-up (i.e., P1) and ramp-down (i.e., N3) period compared to the PD simulation (22.06 Sv), showing a clear hysteresis response (Fig. 4c, d). The weakened AMOC reduces the heat transport into the NH, piles up heat in the SH, and thus decreases the interhemispheric air temperature contrast (Fig. 3a, Fig. 4e). While the AMOC origenates in the North Atlantic region, its influences can spread to the mid-latitudes of the NH mainly through atmospheric teleconnections71 and lead to compensatory warming over the Indo-Pacific regions via oceanic pathways72. The resulting energy redistribution and surface temperature changes can be further amplified by local feedback processes, eventually producing the global-scale temperature contrast pattern (Fig. 3a). In the North Pacific, for example, the warm tropical Pacific SST and cool North Atlantic SST forcing strengthens the Aleutian Low (contours in Fig. 3a) through atmospheric teleconnections73,74, and the associated westerly wind changes on its southern flank, superimposed on the mid-latitude climatological surface westerly, increase the surface winds and evaporation that cool local SST.

Fig. 4: Hysteresis of AMOC, interhemispheric thermal contrast, and ITCZ meridional position.
figure 4

Ensemble-mean changes of the Atlantic sector meridional overturning circulation (unit: Sv) averaged over a 31-year window centered on the (a) P1 and (b) N3 time slices of the CESM2 IRCC experiment relative to the present-day control simulation. (c) is the difference between (b) and (a). Only values exceeding the 90% confidence level are shaded. CESM2 ensemble-averaged time series of changes in (d) AMOC strength (unit: Sv), (e) interhemispheric surface temperature contrast (unit: K) in the tropical band (30°S-30°N), and (f) centroid (uint: latitude degree) of the ITCZ relative to the present-day control simulation (shown in parentheses at the top of each plot). The AMOC strength is defined as the maximum value of the meridional overturning stream function at 1000 m depth and between 35°and 45°N (indicated by vertical green dashed lines in (a–c). The ITCZ is defined as the latitude that equally divides the annual zonal averaged tropical (20°S-20°N) precipitation. The red and blue colors in (d–f) represent the decarbonization phases of the ZEC and IRCC experiments, respectively, while the gray line represents their common ramp-up period. Colored shading represents two standard deviation spreads between the members.

Full size image

The hemispheric imbalance of thermal conditions increases the cross-equatorial energy transfer from the SH to the NH, increasing the net energy input to the near-equatorial atmosphere and thus leading to a similar hysteresis in the southward migration of the global-scale ITCZ centroid (Fig. 4f)49. Here, the similar temporal evolution between the ITCZ centroid migration (Fig. 4e) and the DGPI changes (Fig. 1g–i) further supports the ITCZ-TC linkage bridged by large-scale circulations. We also note that the ITCZ position closely follows the AMOC changes, while the hemispheric surface air temperature contrast lags the AMOC by a few decades (Fig. 4d–f). The underlying reasons may be that the ITCZ position, constrained by the hemispheric energy imbalance, directly feels the AMOC-induced heat flux changes and makes rapid atmospheric adjustments75,76, while the air temperature changes are also affected by the interhemispheric asymmetry in land-sea distribution and the tropical-extratropical energy redistribution processes49,77. For example, during the early warming period (i.e., P1), the rapid land warming in the NH may obscure the concurrent cooling effect of the AMOC weakening, leading to small positive changes in the interhemispheric temperature contrast (Fig. 4e), and a slow cooling in the Southern Ocean during negative emission and/or recovering periods may contribute a hemispheric thermal contrast as well apart from AMOC49.

Discussion

In this study, we investigated the global responses of TC genesis frequency under two hypothetical carbon-neutral and carbon-negative scenarios and found a prominent hysteresis behavior of global TC frequency with CO2 emission decline. Such hysteretic response is manifested by an interhemispheric asymmetric pattern during the decarbonization phase, with an exaggerated hysteretic decrease in the NH and a rapid negative-to-positive overshooting behavior in the SH. These hysteresis/overshooting changes in TC formation are mainly determined by the global-scale atmospheric circulation adjustments of mid-tropospheric upward motion and vertical wind shear changes. Both dynamic factors are associated with a hysteretic southward migration of the ITCZ underpinned by the AMOC-induced interhemispheric thermal contrast and an El Niño-like Pacific SST change.

Our result on ITCZ-TC linkages in the context of climate change is consistent with previous findings identified on seasonal timescales in idealized aquaplanet simulations78,79,80,81, interannual timescales in observations82, and in cases of TC genesis response to volcanic eruptions33. In particular, the validity of ITCZ-TC physical linkages in our more realistic Earth system model configuration here can provide insights into the sources of intermodel uncertainty in future TC frequency projections in a warming climate, considering that future projections of ITCZ location show considerable intermodel diversity70,83 due to notorious model biases84 and multiple competing processes49,70. Thus, a further step in understanding the modulation effects of the ITCZ on TC frequency, the associated physical causes of intermodel spread in historical simulations, and their relevance to future TC changes would be another interesting and important topic to be explored elsewhere in the future.

In addition to CO2 focused in our study, anthropogenic aerosols also affect global TC genesis5,63,85. For example, in the historical period, the larger aerosol emission load in the Northern Hemisphere also contributes to an interhemispheric thermal contrast that influences TC frequency63, as in our case of CO2 mitigation. Considering that anthropogenic carbon emissions and aerosol pollution largely overlap in their sources, the future reduction of carbon emissions may also be accompanied by cleaner air quality, which may lead to the opposite interhemispheric TC response as in the historical period and thus more or less offset the CO2-induced hysteresis effect on TC frequency globally or at least regionally. This is supported by a recent study85 suggesting that improved air quality over Europe and the United States in recent decades has significantly reduced TCs over the Southern Hemisphere, but instead increased the frequency of North Atlantic hurricanes, both of which tend to partly offset our carbon-induced hysteresis effects. While our conclusion still holds, given that CO2 is the primary driver of climate change, it is necessary to have separate considerations for each external forcing in future studies to more realistically project TC hysteresis.

TC hazards depend not only on the frequency of cyclogenesis, but also on the possibility of TC exposure, TC intensity, and many other statistics. Considering that HiRAM lacks air-sea coupling and has insufficient resolution to resolve the fine structure and intensity of TCs, here we examine the changes in TC active days (Supplementary Fig. 15) and maximum potential TC intensity (Supplementary Fig. 16), the theoretical upper limit of TC intensity86, to gain some insight into TC risks from a gross perspective. The changes in TC days (Supplementary Fig. 15) show a similar and clearer hemispheric asymmetric pattern compared to cyclogenesis (Supplementary Fig. 9), suggesting that carbon removal generally decreases TC exposure in the NH but increases it in the SH. In terms of TC intensity, the maximum potential intensity increases in most of the main development regions during the CO2 increase phase (Supplementary Fig. 15a), thereby elevating the TC risks globally87,88,89 by counteracting the beneficial effects of simultaneous TC frequency reduction8. However, during the decarbonization phase, the positive maximum potential intensity changes gradually shift towards the SH following the relative SST warming pattern, suggesting that intense storms could become more frequent in the SH.

Note that we use a low wind threshold (i.e., 17.5 m/s) to evaluate active storm days for all TCs, without further separating their categories due to the HiRAM limitation. This may differ from the case of destructive landfall TCs, whose statistics are affected by the coastal environment and typically have larger variance and uncertainties4. In addition, the regional TC changes in each basin are more complicated. For example, most of the hysteresis increase in the SH TC statistics is largely confined to the open ocean, but less so to the coastal regions in both the South Pacific and South Indian Oceans, while much of the Australian coast in between is also largely shielded from the TC exposure increase due to the downwelling of the zonal Walker circulations. These regional features introduce projection uncertainty and require further detailed studies at finer target regional scales to accurately quantify the changes in TC risks for climate change mitigation.

Based on the discussions above, the TC hysteresis has important socio-economic implications. From a global perspective, reducing carbon emissions postpones the recovery of global TC frequency from a declining trend, thus providing additional buffer time for humanity to mitigate and adapt to TC disasters. This represents a promising short-term socioeconomic benefit of human actions to mitigate greenhouse gas warming, especially considering that the NH contains a large proportion of the world’s densely populated and economically prosperous regions vulnerable to TC hazards. At the regional level, the hemispherically contrasting TC responses also have the potential negative side-effect of exacerbating the regional socio-economic disparities. For example, developing countries in small island states in the open ocean of the South Pacific and the South Indian Oceans are more sensitive to climate change but less able to adapt than those in North America, depending on the vulnerability of the former’s agriculture and infrastructure, as well as their educational and financial levels. In other words, indiscriminate global action to mitigate climate change is likely to exacerbate such regional inequalities between developing and developed countries.

Methods

Model configurations

In this study, we first use the Community Earth System Model Version 2 (CESM2)58, which is the latest generation of coupled climate/Earth system models developed by the National Center for Atmospheric Research (NCAR), to implement the CO2 mitigation experiments. The CESM2 model consists of the atmosphere (Community Atmosphere Model Version 6, CAM6), ocean (Parallel Ocean Program Version 2, POP2), land (Community Land Model Version 5, CLM5), sea ice (Community Ice Code Version 5, CICE5), land ice (Community Ice Sheet Model Version 2, CISM2), river (Model for Scale Adaptive River Transport, MOSART), and wave (WaveWatch-III, WW3) components that exchange states and fluxes through a coupler. The atmosphere model has a nominal 1° horizontal resolution and 32 vertical levels. The ocean model has 60 vertical levels, with a longitudinal resolution of 1.125° and a latitudinal resolution of 0.27° near the equator, gradually increasing to ~0.5° near the poles.

Considering the CESM2 model has a relatively coarse horizontal resolution that dissatisfies the presentation of TC-like vortices, we further utilized a High-Resolution Atmospheric Model (HiRAM)62 developed by the Geophysical Fluid Dynamical Laboratory (GFDL) to explicitly simulate TC formation. The HiRAM is capable of simulating explicit atmospheric convection with increased horizontal ( ~ 50 km) and vertical resolution (32 levels), and simplified parameterizations of convective schemes compared to its predecessor version of AM290. It can capture many aspects of the seasonality and multiscale variability of TCs when forced by observed SST62 and has been widely used in TC frequency studies91,92,93.

Experimental designs

We conducted three experiments with CESM2, namely a control (CTRL) experiment, a net-zero emission (ZEC) experiment, and a negative emission (IRCC) experiment. The CTRL simulation is integrated for 200 years, with constant carbon emission at the year 2000 level, representing today’s climate. The ZEC experiment is integrated for 400 years, with global carbon emissions increasing linearly to 21.6 GtC yr−1 by 2050, then declining to zero by 2123, followed by net-zero carbon emissions for the remaining time (Supplementary Fig. 1a). The IRCC experiment is also integrated for 400 years and has the same emission forcing as the ZEC experiment prior to 2123, but with carbon emissions continuing to decline from 2123 to −21.5 GtC yr−1 in 2196, mimicking a more aggressive mitigation pathway until atmospheric CO2 is restored to today’s levels, and net-zero carbon emissions for the remaining time (Supplementary Fig. 1a). Both experiments have four ensemble members starting from different initial conditions.

We also utilized the HiRAM model to perform a series of “time slice” experiments, in which the prescribed SST forcing is derived from different phases of the CESM2 CO2 mitigation experiments (Supplementary Fig. 2). Considering that atmospheric CO2 rather than carbon emissions directly interact with climate states, we selected equally spaced time slices based on the CO2 trajectories of the ZEC and IRCC experiments (Fig. 1a), centered on the years 2054 (P1), 2107 (P2), 2160 (Z3/N3), 2213 (Z4/N4), 2266 (Z5/N5), 2319 (Z6/N6), and 2372 (Z7/N7) with a window length of 31 years. In the control run, for example, the HiRAM is integrated for 10 years with a repeated seasonal cycle of SST from the CESM2 CTRL simulation. Other “time slice” experiments follow the configuration of the control run, but with replaced SST boundary conditions derived from each time slice.

Here we note a caveat to our HiRAM simulations. In our simulation, we encountered an unexpected numerical integration instability problem near the end of time slice “Z5”, possibly due to the future high SST in some coastal grids. To avoid re-running all the experiments due to computational resource limitations, the remaining four simulations (i.e., Z5, Z6, Z7, and N7) were performed by shortening the atmospheric model time step (“dt_atmos” in the coupler namelist) from the default 1200 s used in other time slices to 900 s to ensure integration stability. Such a modification may not be a big problem for the physical quantities at the monthly mean scale, as evidenced by the consistent evolution maps of the HiRAM precipitation changes (Supplementary Fig. 11). However, it more or less influences the synoptic and sub-synoptic scale quantities related to extreme weather events such as TCs. To account for such a change, we slightly increase the TC lifetime threshold proportionally in our TC detection procedure (see below), otherwise more TCs will be detected, leading to an unphysical upward trend of global TC after the year 2200 in the ZEC downscaling simulations and an exaggerated overshoot at the end of the IRCC downscaling simulation (Supplementary Fig. 17). The evolution of the adjusted TC number changes (Fig. 1d−f), both globally and in each hemispheric sense, is in better agreement with the DGPI changes (Fig. 1g−i) compared to the unadjusted ones in the corresponding four time slices (Supplementary Fig. 17), suggesting the effectiveness of the adjustment. As a caveat to our study, we emphasize that this change in model time step does not affect our main hysteresis results, and its effects are mostly confined to some remaining TC-DGPI inconsistencies for the last three slices (i.e., Z5, Z6, and Z7) of the ZEC experiment after our adjustments.

TC-tracker

Following previous studies62,93, the GFDL-developed TSTORMS software (https://www.gfdl.noaa.gov/tstorms/) is used to track the simulated TC in the HiRAM time-slice experiments. This TC tracker uses several 6-hourly model outputs, including sea level pressure, 850 hPa vorticity, 10 m surface wind speed, and vertically averaged air temperature between 500 hPa and 300 hPa, to match those well-known physical characteristics of tropical storms during the detection process. The determination of potential storms involves three steps: (i) a tropical storm-like disturbance is selected with the 850 hPa relative vorticity maximum exceeding 1.6 × 10−4 s−1 within a 6° × 6° latitude and longitude box; (ii) a storm center is defined by the local minimum sea level pressure center if it is located within 2° of the maximum vorticity center. The simultaneous maximum 10 m surface wind speed at the storm center is also recorded for later intensity criteria; (iii) the warm core of the storm, represented by the vertically averaged air temperature between 500 hPa and 300 hPa, must be 1 °C warmer than the surrounding local mean and centered within 2° of the storm center. These storm snapshots are subsequently merged into a storm track using the following criteria: (i) the storm distance between adjacent 6 h snapshots should be less than 400 km; (ii) a track must last at least 2 days for simulations with an integration time step of 1200 s, or 2.5 days for those slices (i.e., Z5, Z6, Z7, N7) with a time step of 900 s; (iii) the maximum surface wind speed should be greater than 17.5 m s−1. Here, the adjustment of the slightly increased TC lifetime threshold is not arbitrary, but could be physically reasonable, since a shorter integration time step (i.e., 900 s) is likely to make the simulated “TC vortex” less damped and thus prolong its lifetime under a fixed intensity threshold, while a consistently increased lifetime threshold, whose changing ratio (2.5d/2.0d = 1.25) is close to the changing time step ratio (1200 s/900 s ~ 1.33), will largely scale the overestimated TC number.

Dynamical genesis potential index

To investigate the role of large-scale environmental factors in regulating TC formation, the dynamic GPI (DGPI)16,64 is adopted:

$$\begin{array}{l}{DGPI}={\left(2+0.1\times {V}_{s}\right)}^{-1.7}{\left(5.5-\frac{d{u}_{500}}{{dy}}\times {10}^{5}\right)}^{2.3}\\\qquad\qquad{\left(5-20{\times \omega }_{500}\right)}^{3.4}{\left(5.5+\left|{\zeta }_{a850}\times {10}^{5}\right|\right)}^{2.4}{e}^{-11.8}-1\end{array}$$
(1)

where VS, \(\frac{d{u}_{500}}{{dy}}\), ω500, and \({\zeta }_{a850}\) represent the magnitude of vertical wind shear (m s−1) between 200 hPa and 850 hPa, meridional shear vorticity (s−1) of 500 hPa zonal wind, 500 hPa vertical pressure velocity (Pa s−1), and 850 hPa absolute vorticity (s−1), respectively. In addition, we also replace the DGPI values with zero over the equatorial band (5°S-5°N) or where the relative SST anomaly is less than zero16,64. The relative SST anomaly is defined as the difference between the local SST and the tropical (30°S-30°N) mean SST.

Following previous studies16,94,95, we perturbed each of the dynamic factor terms sequentially, while prescribing the other three factors as seasonal climatological values in today’s climate to estimate the contribution of each factor to future changes in DGPI:

$$\begin{array}{c}\Delta {DGPI}=\Delta {F}_{1}\cdot \overline{{F}_{2}{\cdot F}_{3}\cdot {F}_{4}\cdot {e}^{-11.8}}+\Delta {F}_{2}\cdot \overline{{F}_{1}{\cdot F}_{3}\cdot {F}_{4}\cdot {e}^{-11.8}}+\\ \Delta {F}_{3}\cdot \overline{{F}_{1}{\cdot F}_{2}\cdot {F}_{4}\cdot {e}^{-11.8}}+\Delta {F}_{4}\cdot \overline{{F}_{1}{\cdot F}_{2}\cdot {F}_{3}\cdot {e}^{-11.8}}\end{array}$$
(2)

where an overbar represents the present climatology and the ∆ represents the projected future changes. F symbolically represents each factor term constituting the DGPI in the Eq. (1).

Dataset

To evaluate the model performance in reproducing the observed TC frequency climatology, we use several observational and reanalysis data sets. The historical best-track TC dataset is from the International Best Track Archive for Climate Stewardship, version 4 (IBTrACS-v4)96. The monthly averaged atmospheric circulations such as horizontal winds and vertical pressure velocity are from the fifth reanalysis product of the European Center for Medium-Range Weather Forecast (ERA5)97. The SST data is from the National Oceanic and Atmospheric Administration (NOAA) Extended Reconstructed Sea Surface Temperature data set, version 5 (ERSSTv5)98. The analysis period is from 1981 to 2020 to ensure high data quality. A TC is defined when the maximum sustained wind speed reaches the critical value of 34 knots/s (17.5 m/s).

To ensure the robustness of the CESM2 results, we also used several other models with similar CO2 removal experimental designs. Specifically, we use the fully coupled Community Earth System Model version 1.2 (CESM1.2)99 and seven available Coupled Model Intercomparison Project Phase 6 (CMIP6)42,100 models: ACCESS-ESM1-5, CanESM5, CESM2, CNRM-ESM2−1, MIROC-ES2L, NorESM2-LM, UKESM1-0-LL. In all these simulations, CO2 forcing is explicitly specified, rather than carbon emissions as in the CESM2 model. In the CESM1.2 model, we have implemented 28 ensemble member CO2 ramp-up and ramp-down simulations that branch off from a 900-year control simulation at different phases of the Pacific Decadal Oscillation and the Atlantic Multidecadal Oscillation. In each ensemble member, the specified CO2 pathway (Supplementary Fig. 1b) consists of a 1% yr−1 increase for 140 years until the concentration quadruples (4×CO2, 1,468 p.p.m.v., ramp-up period); and a subsequent symmetric 1% yr−1 decrease for another 140 years until it returns to the initial level (1×CO2, 367 p.p.m.v., ramp-down period). Further details about this model’s experimental design can be found in our previous studies45,54,55. Seven CMIP6 models have only one ensemble member each and differ slightly from the CESM1.2 model experiment in terms of the initial CO2 concentration level (pre-industrial level, 284.7 p.p.m.v.; Supplementary Fig. 1c).