Elevated atmospheric CO₂ drove an increase in tropical cyclone intensity during the early Toarcian hyperthermal

During the early Toarcian hyperthermal, large-scale environmental conditions support two TC genesis potential centers across Tethys in our simulations. One is concentrated on northwestern Tethys; the other is distributed in the southeastern part (Fig. 2A). Storm generation over the two centers is relatively more active during July–October (defined as the storm season) in the Northern Hemisphere and during January–April in the Southern Hemisphere (Fig. 2B), which is similar to modern conditions (26). The regions favorable for TC production are generally within the areas bounded by the 26.5 °C isotherm and high potential intensity (>70 m/s) during the storm season (Fig. 2A). In addition, large-scale steering flow that controls the general pattern of storm trajectory (27, 28) shows an anticyclonic circulation over northern Tethys during the storm season (Fig. 2A), with an obvious easterly flow at low latitudes (~5 to 20°N). This favors a storm forming over northwestern Tethys to make landfall along the continental coasts of northern Gondwana. A westward steering flow is observed within the ~15 to 30°S band over southern Tethys (Fig. 2A), conducive to TC landfalls over the southwestern tropical Tethyan margin. Importantly, the modeled locations of active storm activity over Tethys during the early Toarcian hyperthermal are supported by geological occurrences of storm deposits at this time (Fig. 2A).

The enhanced storm activity during the early Toarcian hyperthermal was previously argued to be linked with the CO₂ rise and warming (20–22). To test this hypothesis, we set up a pair of sensitivity experiments and examine how TC may vary with the empirically determined approximate doubling of CO₂ (560 to 1,120 ppmv) during the early Toarcian hyperthermal, as determined from fossil plant stomatal density change (16) and fossil leaf-wax carbon-isotope composition (17) (see also Methods). Our results demonstrate that TC potential intensity is enhanced over almost the entire Tethys Ocean during the storm season in response to higher CO₂ concentration, accompanied by a uniform ocean warming (Fig. 3A). In particular, there is an increased possibility for stronger storms over Tethys during the early Toarcian hyperthermal, manifested by a shift in the tail of the probability distribution toward extremely high potential intensity (e.g., >75 m/s; Fig. 3 B and C). This rightward shift in the probability distribution is more profound over northern Tethys, around which the majority of sedimentary storm deposits during the early Toarcian hyperthermal have been found (Fig. 2A).

Atmospheric CO₂ concentrations during the early Toarcian hyperthermal may have exceeded 2,000 ppmv based on maximum likely values from fossil leaf stomatal density (16), so we further test the variation of potential intensity using a more extreme warming scenario (CO₂ = 2,800 ppmv). The results also point to a significant shift to extremely strong storms under such a high (but still geologically plausible) CO₂ scenario, and the higher the CO₂ concentration, the more extreme the storms (Fig. 3 B and C). This is supported by the findings based on modern continental configurations (29).

Despite an overall warming and higher potential intensity, there are obvious spatial differences in the variation of TC genesis potential across Tethys in response to higher CO₂ (Fig. 4). Genesis potential shows a general increase over northern Tethys during the storm season, with more profound change at the western part. This suggests a westward shift and increased favorability for storm generation. In contrast, there is a quasi-eastward shift in genesis potential over southern Tethys, with more favorable conditions for storm formation in the eastern part. Moreover, under a higher CO₂ scenario, an anticyclonic anomaly in steering flow appears over northwestern Tethys, with anomalous northward/northeastward flows around the continental coasts. This implies that coastal regions of northern Gondwana might be more susceptible to storm strikes and hence increased coastal erosion. Over the tropics of southern Tethys, there is an anomalous eastward flow at the low latitudes (~5 to 20°S) but a westward flow at the relatively higher latitudes, suggesting less favorable conditions for TC landfalls at low latitudes but more favorable conditions for landfalls at middle latitudes.

In addition to the inferred TC activity over Tethys, increases in genesis potential and mean potential intensity are observed over almost the entire tropical ocean in response to increased CO₂ concentration during the early Toarcian hyperthermal (SI Appendix, Fig. S1). In particular, there is an increased possibility for more intense storms across the tropical ocean, including Panthalassa, based on the estimated probability distribution (SI Appendix, Fig. S2). This is supported by evidence for storm deposits on the margin of western Panthalassa during the early Toarcian hyperthermal (24).

The shift in TC genesis potential during the early Toarcian hyperthermal is tied to the changes in large-scale environmental factors that arise from the adjustment of regional temperature structure and the associated variations of atmospheric circulations. In response to increased CO₂ concentration, the entire Tethys Ocean experiences higher potential intensity, but with relatively larger changes in the northwestern and southeastern parts (Figs. 3A and 5 A and E). This is largely attributable to the intensified thermodynamic disequilibrium over the two regions that results mainly from the reduced surface wind speed, with positive contributions from larger thermodynamic efficiency (Fig. 5 B and F and SI Appendix, Fig. S3). Furthermore, the uneven warming in the atmosphere leads to a zonal dipole pattern in meridional temperature gradient in northern Tethys, which results in enhanced (reduced) westerly to the east (west) of ~40°E based on thermal wind balance and hence gives rise to larger (smaller) vertical wind shear (Fig. 5C and SI Appendix, Fig. S4). Similarly, the meridional temperature gradient in southern Tethys exhibits an increase at the western part and a decrease to the east, which closely matches the pattern of vertical wind shear anomaly (Fig. 5G). The anomaly of moisture entropy deficit is characterized by negative and positive values in the western and eastern part of northern Tethys, respectively (Fig. 5D). The reduced moist entropy deficit over northern Tethys is dominated by increased midtropospheric relative humidity, whereas the decreased temperature contrast between the midtroposphere and the surface that weakens the strength of the surface heat fluxes largely contributes to the enhanced moist entropy deficit (Fig. 5D and SI Appendix, Fig. S5). In southern Tethys, the smaller moist entropy deficit in the eastern part is largely caused by increased midtropospheric relative humidity, while the larger moist entropy deficit in the western part is dominated by the dampened temperature contrast (Fig. 5H). Overall, the zonal shift of genesis potential is linked with the variations of potential intensity, vertical wind shear, and moist entropy deficit (Fig. 5), but the relative contributions of individual environmental factors vary across Tethys (SI Appendix, Fig. S6).

Through a set of coupled climate model experiments, our results suggest that the doubling of CO₂ concentrations that likely characterized the early Toarcian hyperthermal led to increased favorability for stronger TCs to develop over the northwestern/southeastern Tethys, as well as other ocean regions at tropical paleolatitudes. Our results confirm the hypothesized linkage between extreme warmth from CO₂ rise and intensified storm activity during the early Toarcian hyperthermal.

In addition to the early Toarcian hyperthermal, increased storm intensity associated with global warming has also been observed in other warmer-than-present time intervals in Earth history. For example, there is growing evidence for increased tempestite deposition (30–32) during the Paleocene Eocene Thermal Maximum (~55 Ma)—an event characterized by a global temperature rise of ~5 to 8 °C (33). An increase in global TC intensity has also been predicted for the Pliocene warm period (~5 to 3 Ma) (34, 35) when the planet warmed by ~3 °C with increased atmospheric CO₂ based on ref. 36. Extremely strong storms were documented over the North Atlantic (37–39) during the peaking warming of the Last Interglacial (~129 to 124 ka) in which global sea surface temperature was ~1 °C higher than present (40) owing to changes in Earth’s orbital configuration and hence insolation. Warmer sea surface temperatures also contributed to the recorded intense hurricanes over the North Atlantic during several intervals of the last two millennia (9). This geological evidence emphasizes the general expectation that increases in TC intensity accompany global warming. Interestingly, observations point to increased mean TC intensity and higher frequency of category 4 to 5 TCs under current warming (1, 3–5), though with differences across ocean basins. Thus, we should expect more intense TCs in a future warmer climate, which is consistent with the majority of future projections (1, 2).

There are a number of limitations in our work and findings to be considered. The model predicts significantly higher SSTs during the early Toarcian hyperthermal in response to higher CO₂ concentration (SI Appendix, Fig. S7), but the magnitude of warming may be underestimated compared to the limited proxy-based estimates that exist (14, 15). This may be partially attributable to the lower climate sensitivity of the climate model used here (41) or limitations with the reliability of proxy temperature estimates based on the chemistry of extinct fauna (i.e., belemnites) (14). However, this discrepancy has limited influence on the qualitative results reported here, as the Earth still experiences much stronger storms in our extreme warming scenario (CO₂ = 2,800 ppmv; Fig. 3 B and C) in which the modeled warming matches better with the proxy data (SI Appendix, Fig. S7). Additionally, potential intensity and genesis potential used here have been proven useful to depict the intensity and location of actual TCs on regional scales, but they provide no information on the behaviors of individual storms (e.g., storm tracks). This hampers the direct comparison of our results with the known storm-related sedimentary deposits. On the other hand, sedimentary storm deposits do not have the stratigraphic and temporal resolution needed to resolve spatial patterns of TC intensity, which in turn is insufficient to constrain the model results. Thus, combinations of high-resolution (e.g., 25-km resolution) climate model or other downscaling techniques (e.g., ref. 42) and additional proxies for TCs would be useful to further assess the linkage between extreme warmth and storm activity in the future. Nevertheless, our results and the available geological evidence now clearly point to increased storm intensity associated with global warming, which is what we observe in current warming. The evidence we present from the early Toarcian hyperthermal aligns with the emerging evidence from other past warm periods that experienced a marked increase in intense storms, which is also fully consistent with future projections determined by modeling (1, 2). Ultimately, therefore, this study highlights the value in studying how past warming events can be used to inform predictions of future extreme climate change.

Sedimentary indicators of storm events have been recorded during the early Toarcian hyperthermal from a number of localities worldwide (Fig. 1). In the western Tethys region, tempestites have been recorded in the Central High Atlas Basin (Dades Valley) of Morocco (21), the Lusitanian Basin of Portugal (21, 43), France (20, 44), and Germany (21, 45). Further east in Tethys, siliciclastic tempestites are recorded in carbonate-dominated successions at Nianduo and Wolong (Kioto Platform), Tibet (22). In Panthalassa, putative tempestites were observed in thin section from a marginal marine succession at Sakuraguchi-dani, Japan (24). Importantly, tempestites in all these locations are nearly always restricted to the early Toarcian hyperthermal (i.e., the carbon isotope excursion) where peak temperatures and maximum atmospheric pCO₂ occurred (Fig. 1). Sequence stratigraphic work has indicated that significant eustatic sea-level rise accompanied the early Toarcian hyperthermal (46), with this rise likely driven by thermal expansion of seawater and polar ice melting (13). As such, the sudden occurrence of storm deposits restricted to the early Toarcian hyperthermal cannot be attributed to a drop in sea level and consequent deposition of strata above storm wave base. On the contrary, rising sea level across the early Toarcian hyperthermal emphasizes that storm events likely increased in magnitude in order to impinge on the sea floor (21). In other locations, sedimentary deposits indicative of high-energy conditions and/or increased sediment fluxes (such as turbidites and hyperpycnites) have been recorded (Fig. 1). Locations include Japan (Sakuraguchi-dani) (24), the United Kingdom (Mochras Farm borehole and Yorkshire) (21, 25), the North Sea (Dutch Central Graben) (47), and Portugal (48). Although not strictly diagnostic of storm activity, these deposits are generally restricted to the early Toarcian hyperthermal and are suggestive of higher-energy conditions and/or increased sediment loading despite increased sea level.

The Community Earth System Model (CESM) version 1.2.2 (49) is a fully coupled global climate model developed by the National Center for Atmospheric Research. The model consists of components for the atmosphere, ocean, land, sea ice, and rivers, and geophysical fluxes across these components are exchanged via a central coupler. The resolution of the atmospheric component, Community Atmosphere Model version 4 (CAM4), is 3.75° in longitude by 3.75° in latitude, with 26 vertical levels. The land component, Community Land Model version 4 (CLM4), has the same horizontal resolution as CAM4. Here, the prognostic carbon–nitrogen model with dynamic vegetation is turned on to simulate unmanaged vegetation including tree, grass, and shrub plant functional types. The river transport model directs all runoff to oceans and is run with the default resolution (0.5° × 0.5°). The ocean component, Parallel Ocean Program version 2 (POP2), has a nominal 3° irregular horizontal grid (116 and 100 grid points in the meridional and zonal directions, respectively), with 60 vertical levels. The sea ice component (i.e., Community Sea Ice Model version 4) is run on the same horizontal grid as POP2. The CESM model has been validated for modern climates and has also been successfully implemented in a variety of paleoclimate studies spanning from the Precambrian (~540 Ma) (50, 51) to mid-Holocene (~6 ka) (52). Equally, CESM has relatively better performance in depicting the key features of past climates (53, 54). Thus, we used the fully coupled CESM version 1.2.2 to provide state-of-the-art computer simulations of past climate states.

The paleogeographic map for the early Toarcian (SI Appendix, Fig. S8) is from the paleo-digital elevation model (paleoDEM) of ref. 55, which is a digital representation of paleotopography and paleobathymetry reconstructed back in time. The paleoDEM describes distributions of deep oceans, shallow seas, lowlands, and mountainous regions, and it is an estimate of the elevation of the land surface and depth of the ocean basins. Some modifications are made for the Tethys following refs. 21 and 56, and we do not consider the isolated small islands due to the coarse resolution used here. No ice sheets are prescribed in the simulations. The paleogeography is produced at an original resolution of 1° × 1°, from which we generated model-resolution land–sea mask, topography, and bathymetry.

To investigate the sensitivity of TC activity to atmospheric CO₂ concentrations, three simulations are performed with 2×, 4×, and 10× the preindustrial concentration (280 ppmv). Estimates of early Toarcian hyperthermal CO₂ concentrations and changes are inherently uncertain given the limitations of the available proxies (e.g., ref. 17 and references therein). However, the 2× preindustrial concentration (560 ppmv) is close to the pre-early Toarcian hyperthermal estimate of atmospheric CO₂ determined by both fossil leaf stomatal density data (16) and fossil leaf wax carbon-isotope data (17) (~500 to 700 ppmv). The 4× preindustrial concentration (1,120 ppmv) is close to the estimated peak CO₂ concentration of the early Toarcian hyperthermal based on fossil leaf wax carbon-isotope data (17) (~1,000 ppmv) and close to the minimum concentration estimated by fossil leaf stomatal density data (16) (~1,250 ppmv). The 10× preindustrial concentration (2,800 ppmv) exceeds the maximum estimated peak CO₂ concentration of the early Toarcian hyperthermal based on fossil leaf wax carbon-isotope data (17) (~1,200 ppmv) but is closer to the maximum estimated concentration based specifically on fossil Ginkgoales stomatal density data (16) (~2,300 ppmv). The solar constant is set to 1341.4 W m⁻² for the early Toarcian, considering that solar radiation increases linearly at an approximate rate of 0.08% per 10 Myr (57) to the present-day value of 1360.89 W m⁻². Orbital parameters in the simulations are set the same as the present-day values.

Following the guidelines of the Deep-Time Model Intercomparison Project (58), the simulations are integrated for ≥2,000 model y to reach equilibrium states at which the net radiation at the top of the atmosphere averaged over the last 100 model y is less than 0.1 W m⁻² (SI Appendix, Table S1). Monthly outputs from the last 100 model y in each experiment are archived and used for the relevant statistical analysis (i.e., the calculation of potential intensity and genesis potential).

Potential intensity is the theoretical maximum intensity of a storm constrained by the thermodynamic state of the atmosphere and sea surface (59, 60). It has been proven useful to measure the climatology and variability of actual storm intensity (61–63). Potential intensity is defined as

where $C_{\mathrm{k}}$ and $C_{\mathrm{d}}$ are the enthalpy and momentum surface exchange coefficients, respectively, $T_{\mathrm{o}}$ is the mean outflow temperature, $h_0^{\ast }$ is the enthalpy of air saturated at the sea surface temperature and pressure, and $h$ is the enthalpy of an ambient boundary layer parcel. The term $\frac{SST-T_{\mathrm{o}}}{T_{\mathrm{o}}}$ depicts the thermodynamic efficiency, and the term $(h_0^{\ast }-h)$ measures thermodynamic disequilibrium (61). We estimate the contributions from the two factors by squaring and taking the natural logarithm of Eq. 1:

We employ a Genesis Potential Index (GPI) that summarizes large-scale environmental factors that are important to TC formation to predict the regions where storms may generate. GPI has high skills in predicting the spatial pattern of actual TC genesis (26, 63, 64) (e.g., SI Appendix, Fig. S9) and shows similar performance compared with downscaling method (65, 66), and hence are widely employed to study TC activity in past, present, and future climates (67–70). The GPI used here is defined as

where a is a normalizing coefficient, PI is the potential intensity, VS is the vertical wind shear, $\chi$ is the moist entropy deficit, and $\eta$ is the low-level absolute vorticity. Briefly, vertical wind shear is defined as the magnitude of the vector difference between the 200 and 850 hPa horizontal winds, as the most profound difference in TC environmental wind profiles generally appears between these two levels (71), and the vertical wind shear measured between 200 and 850 hPa generally has the strongest negative correlation with storm intensity (72). Vertical wind shear generally weakens storm formation and development via shearing the convective towers and ventilating the storm’s core with subsaturated air (73). Moist entropy deficit is adopted to assess the midtropospheric moisture content (at 600 hPa) and is determined by the ratio of the difference in moist entropy that moist convection must eliminate to the strength of the surface fluxes supplying the moisture (65–67):

where $s_0^{*}$ and $s^{*}$ are the saturation moist entropies of the sea surface and free troposphere (600 hPa), respectively, and $s_m$ represents the moist entropy of the middle troposphere (600 hPa). Larger moist entropy deficit suggests less favorable conditions for TC generation and vice versa. Low-level absolute vorticity is computed at 850 hPa and acts as a spin-up mechanism by inducing important synoptic convergence that favors cyclone formation (74). Overall, regions of high genesis potential identify the areas in which TC genesis is possible; low values are found only in regions that cannot support the deep convection that defines such systems. More information concerning individual genesis factors and GPI is given in ref. 69.

D.B.K., Z.Z., and Y.H. designed research; X.L. and J.G. performed research; Q.Y. analyzed data; and Q.Y. and D.B.K. wrote the paper.

The authors declare no competing interest.