Modeling hyperthermal events in the Mesozoic-Paleogene periods: a review

Introduction

Hyperthermal events, characterized by rapid and extreme global warming, have significant impacts on the Earth’s surface environment and biosphere through multiple factors such as ocean acidification, anoxia, acid rain, and aridity, ultimately resulting in mass extinctions (Montoya-Pino et al., 2010; Clarkson et al., 2015; Penn et al., 2018; Shen et al., 2019; Rush et al., 2021). At the present day, we are living through a >1°C global warming event caused by anthropogenic carbon emissions (IPCC, 2022). However, the mechanisms and environmental effects of past global warming remain unclear, hindering our understanding of the mechanisms of global warming and its bio-environmental effects. Climate warming extremes have occurred numerous times, particularly during the Mesozoic to Paleogene eras, which saw the Permian-Triassic Boundary event (PTB, ~252 Ma), early Toarcian Oceanic Anoxic Event (TOAE, ~183 Ma), early Aptian oceanic anoxic event (OAE 1a, ~120 Ma), Cenomanian-Turonian oceanic anoxic event (OAE 2, ~94Ma), and the Paleocene-Eocene Thermal Maximum (PETM, ~56 Ma) (see Figure 1) (Hu et al., 2020).

Figure 1 Hyperthermal events, their possible related large volcanism activities, and their typical seawater carbon isotope pattern from the Mesozoic to Paleogene. Note: the seawater carbon isotope patterns are from Kuroda et al. (2007), Di Lucia et al. (2012), Owens et al. (2013), Shen et al. (2013), Suan et al. (2015), and Minshull et al. (2016).

These past hyperthermal events offer valuable geological windows for studying the mechanisms of global warming and its bio-environmental effects. For instance, the PTB event at the end of the Permian (~252 Ma) was characterized by a rapid increase in atmospheric carbon dioxide (CO₂), leading to global warming of 6–11°C (Penn et al., 2018; Shen et al., 2019), ocean acidification (Payne et al., 2010; Clarkson et al., 2015; Silva-Tamayo et al., 2018), and marine anoxia (Song et al., 2012). The PTB event resulted in the extinction of 81% of marine species (Stanley, 2016) and 70% of terrestrial species (Sahney and Benton, 2008). The carbon emission responsible for the PTB hyperthermal event has been linked to the Siberian Traps Large Igneous Province (LIP) based on their strong temporal association (Joachimski et al., 2020). Similarly, other hyperthermal events (TOAE, OAE 1a, OAE 2, and PETM) are also believed to be causally linked to LIPs that share a temporal association with them. These associated LIPs are Karoo-Ferrar Province (Percival et al., 2016), Ontong-Java Plateau (Bauer et al., 2017), Caribbean-Madagascar flood basalts (Turgeon and Creaser, 2008), and North Atlantic Igneous Province (Gutjahr et al., 2017), respectively (see Figure 1).

Numerical modeling has proven to be an effective approach for quantitatively investigating the driving mechanisms for hyperthermal events, including the identification of potential carbon sources, the amount of carbon that drove warming (e.g., Cui et al., 2021b; Wu et al., 2023), and the resulting environmental impacts. These models span a range of complexity and spatial resolutions: from simple box models and Earth system box models to 3D Earth system models. Simple box models represent the major geochemical or geophysical processes in the geological cycles of the elements or isotopes of interest, by treating each reservoir as a box and calculating the fluxes among these reservoirs. ‘Earth system box models’, which incorporate the interactions among the primary chemical cycles in the atmosphere, ocean, land, and biosphere, represent these processes with a higher level of complexity than simple box models for single chemical species.

One of the most widely used Earth system box models is the Carbon-Oxygen-Phosphorus-Sulfur-Evolution model (COPSE) (Bergman et al., 2004; Lenton et al., 2018) which involves the surficial C, O, P, and S cycles among the atmosphere, ocean, land, sediments, and biosphere. For example, Ullmann et al. (2020) used the COPSE model to predict seawater C isotopes, temperature, weathering, nutrient dynamics, and ocean anoxia intensity to compare with geological evidence for the TOAE event, suggesting that the C cycle perturbation of the TOAE was caused by the injection of carbon derived from the interaction of LIPs with sedimentary organic matter and methane.

In contrast to the simple/Earth system box models, 3D Earth system models integrate the interactions of the atmosphere, ocean, land, ice, and biosphere across a spatially-defined grid, focusing on the expressions of physiochemical processes and energy exchange, allowing estimation of local physiochemical processes and states of each grid. Like the box models, 3D models use a defined number of boxes or containers to store information about concentrations of species of interest, but in these models, the boxes are assigned specific spatial locations, thus they are not usually termed to be ‘box models’. An example of a low-resolution 3D Earth system model (intermediate complexity) is the Carbon centric-Grid Enabled Integrated Earth system model (cGENIE; Ridgwell et al., 2007), while an example of a high-resolution 3D Earth System Model is the Community Earth System Model (CESM, available at https://www.cesm.ucar.edu). Even the low-resolution 3D models have substantially longer computation times than the simple box models or Earth system box models, meaning they typically produce ‘snapshots’ of climate and biogeochemistry for a certain time in the past or geologically short (<100 kyr) transient simulations.

This review explores recent modeling studies of hyperthermal events and their implications for warming mechanisms and bio-environmental influences. As hyperthermal events typically last more than 100 kyr, we will particularly focus on simple box models and Earth system box models, which can be applied for these timescales, and will discuss how these models are being extended to include 3D representations of steady-state climate (Goddéris et al., 2014; Mills et al., 2021), bringing their functionality closer in line with the 3D Earth system models.

Previous modeling work on climate extremes

In this section, we briefly introduce each hyperthermal event and summarize the modeling work on these events chronologically: PTB, TOAE, OAE 1a, OAE 2, and PETM.

The PTB event

As previously mentioned, the PTB event led to the extinction of most life on Earth, both on land and in the ocean, due to a combination of factors such as ocean acidification, anoxia, and acid rain (Clarkson et al., 2015; Penn et al., 2018; Dal Corso et al., 2022; Wang et al., 2023). Although the interplay of kill mechanisms during this event is still debated, it is widely agreed to have been ultimately triggered by massive amounts of greenhouse gases such as CO₂ released into the atmosphere during the emplacement of the Siberian Traps LIP (Figure 2). This process led to a global temperature increase of approximately 6–11°C within a few thousand years, as evidenced by oxygen isotope records (Penn et al., 2018; Shen et al., 2019). The PTB event is characterized by a negative seawater C isotope excursion of 3–8‰ (Shen et al., 2013; Joachimski et al., 2020) (see Figure 1 and Table 1).

Figure 2 Driving mechanisms and environmental impacts of the PTB hyperthermal event (modified from Hu et al., 2020).

Table 1 Comparison of hyperthermal events showing the global warming magnitude, carbon isotope excursion, duration, and proposed possible carbon sources.

Using calculated annual air temperature and runoff from a one-dimensional energy balance model, Grard et al. (2005) applied a box model of biogeochemical carbon and alkalinity cycles to investigate the negative δ¹³C excursion (NCE) driver. By applying time-dependent forcing functions for the eruptive phase and mass extinctions, the model generated a seawater NCE of approximately −3‰. The authors suggested that volcanic eruptions with ~19,200 Pg mantle carbon with a −5‰ δ¹³C value (Table 1) could generate the negative C isotope anomaly, while a reduction in primary productivity that decreased the removal of ¹³C-depleted organic carbon from the ocean also contributed. However, when combined with proxy data for atmospheric CO₂, the carbon source was determined to have significant contributions of ¹³C-depleted carbon (<−5‰) (Table 1) (Payne et al., 2010; Wu et al., 2021).

Observations of negative calcium (Ca) isotope excursion during the PTB, when combined with a box model that coupled the C and Ca cycles, suggested that the average δ¹³C of CO₂ released was more likely near −15‰ (Payne et al., 2010). This conclusion was further supported by simulations with the Long-term Ocean-atmosphere Sediment CArbon cycle Reservoir Model (LOSCAR; a type of box model, Zeebe, 2012), where a combination of mantle sources (~−5‰) and organic matter-derived CO₂ (~−30‰) (Shen et al., 2022) or volcanic CO₂ sources ranging from −8 to −12‰ (Wang et al., 2023) was required to match the observed NCE pattern (Table 1). Similarly, a massive (~36,000 Pg C) and rapid emission of C source (~−15‰) (Cui et al., 2021b) or two-phased C release (1. ~5,000 Pg thermogenic C of −60 to −20‰, 2. ~21,000 Pg volcanic C of −20 to 0‰) (Wu et al., 2023) was required to drive the observed NCE pattern (Table 1) in the cGENIE Earth system model with intermediate complexity (Ridgwell et al., 2007).

In addition to thermogenic or volcanic CO₂ release, the collapse of vegetation on land has also been suggested to impact the NCE pattern significantly. This is particularly true for explaining the sudden negative shift superimposed on the relatively gentle negative excursion (Figure 1). Such a shift has been explained through the use of box modeling of coupled C, P, and mercury (Hg) cycles (Dal Corso et al., 2020) and the strontium (Sr) cycle (Song et al., 2015).

Combinations of other geochemical proxies with numerical modeling techniques for elemental cycles have led to further hypotheses for environmental change during the PTB event. By utilizing a simple Sr ocean box model to simulate the accelerated rise of seawater ⁸⁷Sr/⁸⁶Sr during the PTB event, it was proposed that the continental weathering rate in the Early Triassic was enhanced by more than 1.9 times compared to the late Permian (Table 2), due to the intensification of runoff and vegetation die-off (Song et al., 2015). Similarly, a simple lithium (Li) ocean box model and light seawater Li isotopic signatures during the P–T transition were interpreted as evidence for enhanced continental weathering (Sun et al., 2018).

Table 2 Comparison of hyperthermal events showing modeling weathering enhancement factors and related proxy records.

In contrast to the consensus on the increased continental weathering, primary productivity during the P–T transition is debated. High nutrient levels and therefore high primary productivity was suggested using the cGENIE model to explain the gradient of C isotopes versus water depth (Meyer et al., 2011), while the nitrogen isotope record indicated lower nutrient levels and low productivity caused by suppressed upwelling (Grasby et al., 2016). Due to uncertain changes in primary productivity, its variations during the P–T transition were not included in Figure 2.

Ocean acidification is believed to have played a major role in the P–T mass extinction, as evidenced by seawater boron (B) and Ca isotopes and models of those systems (Payne et al., 2010; Clarkson et al., 2015; Silva-Tamayo et al., 2018). The observed negative isotopic excursions of C and Ca, along with a forward box model of the C and Ca cycles, suggested a CO₂-driven ocean acidification scenario (Figure 2) accompanied by intensified weathering of Ca-bearing rocks during the PTB event (Payne et al., 2010; Silva-Tamayo et al., 2018). With further constraints on seawater pH variations from B isotopes, a box model of the C cycle was utilized to simulate seawater C and B isotope fluctuations during the PTB, suggesting that the acidifying effects of CO₂ released during the formation of the Siberian Traps were buffered for approximately 150 kyr until global ocean acidification occurred (Clarkson et al., 2015). In addition to ocean acidification, global warming could also lead to increased ocean stratification (Figure 2) and decreased ocean circulation due to lower oxygen solubility and increased runoff into high-latitude oceans, as demonstrated by the CESM (Penn et al., 2018). These factors could help account for the extinction of marine species during the PTB event (Penn et al., 2018).

There has been relatively little modeling work on the demise of the PTB event. Existing work suggests the termination of the warming was related to the weakening of volcanism (Dal Corso et al., 2020), enhanced continental weathering (Song et al., 2015; Sun et al., 2018), and processes such as organic carbon burial and carbonate formation (Silva-Tamayo et al., 2018) which played important roles in reducing atmospheric CO₂ levels back to pre-event levels.

TOAE

TOAE, also called the ‘Jenkyns event’, occurred during the Early Jurassic period (~183 Ma) and was characterized by the deposition of organic-rich sediments on a global scale (Suan et al., 2013). The event is widely attributed to a sudden increase in atmospheric CO₂ levels, which is thought to have led to a greenhouse warming effect of over 5°C (Figure 3 and Table 1) (McElwain et al., 2005; Suan et al., 2015). Isotope analysis reveals that the TOAE is associated with an NCE in both marine and terrestrial environments, with a magnitude of −3‰ to −8‰ (Suan et al., 2015; Them et al., 2017a) (Figure 1). This event has been linked to the early magmatism of the Karoo-Ferrar LIP (Figure 1) (Percival et al., 2016).

Figure 3 Driving mechanisms and environmental impacts of TOAE (modified from Hu et al., 2020).

The perturbation in the carbon cycle was initially attributed to the sudden release of light C from either magma-intruded organic-rich rocks (Kump and Arthur, 1999) or marine gas hydrates (Hesselbo et al., 2000). McElwain et al. (2005) estimated the marine gas hydrates reservoir and conducted a simple carbon cycle calculation, using CO₂ proxy data and the NCE pattern to support the carbon source from magma-intruded organic-rich rock (~2,600–4,400 Pg C) (Figure 3 and Table 1). Subsequently, negative Ca isotope excursions found during the TOAE, together with an ocean box model of the Ca cycle, were used to constrain the continental weathering input to the ocean and the causes of the NCE (Brazier et al., 2015). This study suggested that injections of carbon with a range of −6‰ to −14‰ were the most likely cause for the NCE, supporting the carbon source of thermogenic CO₂ (McElwain et al., 2005). However, recent COPSE modeling of the NCE and temperature variations during TOAE suggests that the interactions of the Karoo-Ferrar large igneous province with organic C (~−7‰) and methene (~−60‰) both contributed to the global warming and the NCE (Ullmann et al., 2020) (Table 1).

Regarding the impact of the TOAE warming event, Brazier et al. (2015) utilized an ocean box model of the Ca cycle to study the causes of the negative δ⁴⁴/⁴⁰Ca excursions during the TOAE, and proposed that a transient 90% decrease of carbonate accumulation due to ocean acidification followed by a 5-fold increase in continental weathering rates (Table 2) were required to match the Ca isotope record. However, this estimated enhance factor for weathering was later adjusted to 2.15–5.30 times (Table 2) compared to the pre-event baseline by Them et al. (2017b) based on a positive osmium (Os) isotope excursion and box modeling of the Os cycle. Finally, the COPSE model was used to match the NCE and temperature variations, and it was found that the enhanced factor of continental weathering required was approximately 1.67 times (Table 2) (Ullmann et al., 2020).

Continental weathering also played a vital role in ending the TOAE by enhancing CO₂ sequestration as a negative feedback mechanism (Them et al., 2017b). Moreover, the increased fluxes from continental weathering would have brought more nutrients into the ocean, increasing organic carbon burial (Ullmann et al., 2020). This, in turn, likely contributed to the termination of the TOAE. The increased organic carbon burial was suggested to be caused by stratification and anoxia resulting from a sluggish ocean circulation (Figure 3), which was simulated in general circulation models such as the Fast Ocean-Atmosphere Model (FOAM) (Dera and Donnadieu, 2012; Baroni et al., 2018) and Princeton Ocean Model (POM) (Bjerrum et al., 2001).

The sulfur (S) cycle during the TOAE also drew the attention of geologists, especially in reconstructing seawater sulfate concentrations. Forward box modeling of the S cycle was used to simulate the extended recovery of seawater sulfate S isotopes after the TOAE. This suggested a Toarcian marine sulfate concentration of 4–8 mM (Gill et al., 2011). This low sulfate concentration was further supported by additional seawater sulfate S isotope data and the maximum sulfate concentrations (Han et al., 2022) calculated using the ‘rate method’, which infers the sulfate concentration based on the rate of change of sulfur isotopes (Algeo et al., 2015).

OAE 1a

OAE 1a, which occurred in the early Aptian (~120 Ma) during the mid-Cretaceous, was a period of rapid and extreme global warming (Figure 4), which has been estimated to be between 2 and 6°C (Ando et al., 2008). This event was also associated with widespread anoxia in the oceans and high atmospheric CO₂ concentrations (Jenkyns, 2010; Percival et al., 2021), and had significant impacts on both marine and terrestrial ecosystems, causing changes in nutrient cycles, shifts in primary productivity, and biological turnover (e.g., Keller et al., 2011). The event is marked by widespread preservation of organic matter in marine settings (Jenkyns, 2010) and a significant NCE at its onset, followed by a broad positive carbon isotope excursion (PCE) with a magnitude of 2–5‰ (Figure 1) (Di Lucia et al., 2012). This event has been causally linked to the eruption of the Ontong Java Plateau LIP (Figure 1) (Tejada et al., 2009; Keller et al., 2011; Kashiwagi, 2016; Bauer et al., 2017).

Figure 4 Driving mechanisms and environmental impacts of OAE 1a (modified from Hu et al., 2020).

Based on large un-radiogenic Os isotope shifts during OAE 1a, box modeling of the Os isotope excursions suggests a rapid 6–10 times increase in volcanic activity during the event (Bauer et al., 2017) and/or weathering of a huge amount (~30–60%) of the Ontong Jave Plateau (Tejada et al., 2009). Additionally, the LOSCAR C-cycle box model (Zeebe, 2012) was used to investigate the impact of increased volcanic emissions on the carbon cycle, suggesting that a ~4 times rise in mantle carbon emissions (δ¹³C: ~−6‰) (Table 1) could explain the abrupt negative carbon isotope excursion (NCE) during the OAE 1a alone (Bauer et al., 2017).

Subsequently, the cGENIE model was employed to assimilate an atmospheric CO₂ reconstruction and the carbon isotope record of OAE 1a, suggesting a shift from a ¹³C-depleted reservoir (e.g., thermo- and/or biogenic methane from sill intrusions: −70 to −30‰) towards a dominantly volcanic source (e.g., LIPs: −20 to 0‰) (Table 1) (Adloff et al., 2020). Recently, hydrothermal activity has been proposed as the primary cause of the un-radiogenic Os isotope shift based on additional Os isotope data (Matsumoto et al., 2022), which may support the hypothesis of a shifting carbon source (Adloff et al., 2020).

To explain the negative Ca and Sr isotope shifts, an oceanic box model with coupled Sr and Ca isotope cycles was used to show that a threefold increase in continental weathering (Table 2) and a 20% increase in hydrothermal activity for 5 Myr could drive these changes (Blättler et al., 2011). The relative contributions of hydrothermal activity and continental weathering were later revised to a 1.7-fold increase in hydrothermal activity with a 2.5-fold enhancement in continental weathering (Table 2) for 6 Myr to match the Li, Sr, Os, and Ca isotope variations during OAE 1a in an ocean box model of coupled elemental cycles (Lechler et al., 2015). To match the S isotope record under low seawater sulfate concentrations, a 4.2-fold increase in hydrothermal inputs and a threefold increase in weathering fluxes (Table 2) over 500 kyr, combined with increased evaporite and pyrite burial, was necessary within an ocean box model of the S cycle (Bauer et al., 2022). An S-cycle box model previously showed that increased weathering inputs would increase marine sulfate concentration from a low background value, which was responsible for decreased seawater δ³⁴S isotopes (Gomes et al., 2016).

As for the termination of OAE 1a, a carbon cycle model based on the GEOCARB model (Berner and Kothavala, 2001) was applied to study the atmospheric CO₂ variations (Kashiwagi, 2016). The study suggested that the end of the OAE 1a event was associated with decreases in CO₂, which was likely due to a combination of factors, including reduced volcanic activity, enhanced continental weathering, carbonate precipitation, and organic C burial (Kashiwagi, 2016; Bauer et al., 2017; Adloff et al., 2020).

OAE 2

Cretaceous OAE 2, also known as the Cenomanian-Turonian extinction and the Bonarelli Event, occurred around 94 million years ago and lasted up to ~900 kyr (Clarkson et al., 2018; Hu et al., 2020) (Table 1). This event was characterized by a rapid increase in atmospheric CO₂ levels, possibly caused by Madagascar and/or Caribbean LIP (Figure 1) (Turgeon and Creaser, 2008), which led to global warming of around 2–10°C (Forster et al., 2007; Huber et al., 2018) and marine anoxia (Clarkson et al., 2018) (Figure 5). This event ultimately caused widespread devastation of marine ecosystems (Coccioni, 2004) and widespread deposition of organic matter (Jenkyns, 2010). The carbon cycle perturbations during the event were globally recorded as a broad positive carbon isotope excursion with a magnitude of around 2–5‰, possibly with an abrupt ~3‰ NCE superimposed on the onset of the PCE (Figure 1) (Kuroda et al., 2007; Jenkyns, 2010; Owens et al., 2013).

Figure 5 Driving mechanisms and environmental impacts of OAE 2 (modified from Hu et al., 2020).

In 2011, a modeling study by Flögel et al. utilized a combination of biogeochemical box models (the Cretaceous Ocean Box Model with 27 ocean boxes, COBM) and an Earth system model (known as GENESIS) to investigate the carbon cycle perturbations during the OAE 2 event. The study confirmed that under a volcanic LIP degassing scenario with a δ¹³C value of −5‰ (Wallmann, 2001), high productivity and its following organic carbon burial due to increased nutrient supply from continental weathering, could produce a ~2.5‰ PCE (Figure 5), and lead to >90% global deep-ocean anoxia. In addition, Flögel et al. (2011) proposed that C:P ratios probably increased under high CO₂ conditions in a well-ventilated ocean, contributing to ocean anoxia. With a light carbon source with a δ¹³C value of −5‰, ~19,000–33,000 Pg C (Table 1) was required for an abrupt ~3‰ NCE, possibly superimposed on the onset of the PCE (Figure 1) in a simple one-box model of the atmosphere–ocean C cycle (Kuroda et al., 2007). The different C isotope patterns produced in these two models probably highlight a significant role of primary productivity (and also varied C:P ratios) in the C cycle of OAE 2 (Flögel et al., 2011).

During OAE 2, the varying C:P ratios were likely caused by increased P supply from land and P recycling back into seawater under low oxygen conditions (Nederbragt et al., 2004; Tsandev and Slomp, 2009; Flögel et al., 2011). By applying a box model of oceanic P, C, and O cycles (Slomp and Van Cappellen, 2007), Tsandev and Slomp (2009) argued that oceanic anoxia could be triggered by enhanced continental weathering, which provided more P supply into the ocean. They also suggested that the C-P-O cycle is particularly sensitive to oceanic mixing. Aside from increased weathering P inputs, the oceanic P concentrations can be increased by enhanced benthic (seafloor) regeneration of organic P under anoxic conditions. This can promote primary productivity, and organic C burial, as shown in the cGENIE model (Monteiro et al., 2012), and then generate a rapid 2‰ PCE in coupled C-P models (Nederbragt et al., 2004). Additionally, ocean circulation can influence local primary productivity and organic C burial by changing oxygen availability and nutrient concentrations between ocean basins. Studies using the CESM (Trabucho Alexandre et al., 2010) and a regional ocean circulation model (Topper et al., 2011) have shown this to be the case.

Increases in weathering fluxes were also required in the box modeling of other element cycles, such as Ca, Sr, Os, Li, and S (Blättler et al., 2011; Pogge Von Strandmann et al., 2013; Gomes et al., 2016; Nana Yobo et al., 2021). For example, a threefold increase in weathering flux (Table 2) was required in an oceanic box model of coupled Sr and Ca cycles to explain their behavior during OAE 2 (Blättler et al., 2011). Later, the increased factor was revised to be approximately 1–3 times (Table 2) when the Os and Li isotopic variations during the OAE 2 were also considered in a box model of coupled Li and Os cycles (Pogge Von Strandmann et al., 2013). More recently, box modeling of the increases in seawater ⁸⁷Sr/⁸⁶Sr suggested a rise of only 1.8 times in the weathering flux of Sr (Table 2), but with a more significant role of enhanced hydrothermal activity in the event (Nana Yobo et al., 2021).

Simple box modeling has also been used to investigate the extent of marine anoxia during OAE 2, particularly through the uranium (U) cycle (Montoya-Pino et al., 2010; Flögel et al., 2011; Clarkson et al., 2018; McDonald et al., 2022). U isotope compositions of seawater (δ²³⁸U) were represented in a simple mass-balance equation of the U cycle, to calculate the anoxic area of the ocean sedimentary environment. Here, an increase from a modern value of 0.3% to ~1–2% of the anoxic seafloor was suggested to have occurred during OAE 2 (Montoya-Pino et al., 2010). However, subsequent studies have suggested different extents of marine anoxia during this event. For example, the GENESIS-COBM model proposed an anoxic extent of >90% of the deep ocean (Flögel et al., 2011), while the cGENIE model suggested an extent of >50% of the ocean volume could have been anoxic (Monteiro et al., 2012) (Table 3).

Table 3 Summary of assumed oceanic anoxic extent during OAE 2.

These high estimates were questioned by Owens et al. (2013), who used a forward geochemical box model of C and S cycles and suggested that only ~2.5–5% of global seafloor area (Table 3) was overlain by euxinic waters during the OAE 2 interval. In contrast, a dynamic box model of coupled C, P, and U cycles proposed that ~5–15% anoxia area (Table 3) was required to reproduce the negative δ²³⁸U excursion with the PCE during OAE 2 (Clarkson et al., 2018). Recently, McDonald et al. (2022) revisited the negative U isotope excursion and found a larger magnitude of ~0.9‰, and re-modeled it in a dynamic U box model, which suggested that the anoxic area of the seafloor was only ~21% during the event (Table 3). The differences in the U isotope records and indicated anoxic areas indicate that anoxia may differ between ocean regions (Table 3). This is supported by the inconsistent I/Ca ratio records on a global scale and high oxygen availabilities of some regional deep waters simulated in the cGENIE model (Zhou et al., 2015).

The periodicity of redox conditions during the Cretaceous, as suggested by the in-turn development of black shales and red beds during the Cretaceous (Wang et al., 2011), has also been studied using modeling approaches. Handoh and Lenton (2003) noticed periodicities of ~5 Myr in the C, P cycles and seawater redox conditions, and hypothesized it to be a self-sustaining oscillation in response to an increase in the P supply from the land above a critical threshold. By applying a model of coupled C, N, P, and O cycles (a predecessor of the COPSE model), they suggested that the oscillations were primarily maintained by positive feedback between P concentration, biological productivity, and anoxia in the ocean (Handoh and Lenton, 2003). In addition to the long-term oscillation, a distinct short-term redox cyclicity (~1 Myr) was recognized during OAE 2 in iron speciation data (Poulton et al., 2015). In their study, new S isotope data and a simple S isotope box model indicated that fluctuations in the weathering influxes of S and reactive iron due to the dynamic hydrological cycle likely drove this short-term redox cyclicity.

Therefore, phosphorus cycling during OAE 2 likely played a critical role in both the initiation and cessation of ocean anoxia. Organic carbon burial (e.g., Monteiro et al., 2012), as well as enhanced continent weathering (e.g., Pogge Von Strandmann et al., 2013), were both essential processes that contributed to the reduction of atmospheric CO₂ levels (Kashiwagi, 2016), ultimately leading to the termination of the OAE 2 event.

The PETM event

The PETM event (Figure 6), which occurred approximately 56 million years ago and lasted for 150–220 thousand years, saw a global average temperature increase of around 5–8°C (Röhl et al., 2007; McInerney and Wing, 2011; Frieling et al., 2017; Haynes and Hönisch, 2020). The North Atlantic Igneous Province (Figure 1) (Kraus and Riggins, 2007; Gutjahr et al., 2017; Jones et al., 2019) and/or methane clathrate release (Dickens et al., 1995; Dickens, 2011) have been proposed as the causes of the onset of this event, leading to extreme changes in the C cycle. The perturbation of the carbon cycle was recorded as significant NCEs in both marine and terrestrial carbon reservoirs (Koch et al., 1992; Cui et al., 2021a), with magnitudes ranging from −2‰ to −7‰ (Figure 1) (Koch et al., 1992; McInerney and Wing, 2011; Minshull et al., 2016). This event was also associated with profound biological consequences, including the largest mammalian faunal turnover of the Cenozoic (Koch et al., 1992).

Figure 6 Driving mechanisms and environmental impacts of the PETM event (modified from Hu et al., 2020).

The origin of the ¹³C-depleted carbon source for the PETM event is still a matter of debate. In 1995, a carbon mass balance calculation suggested that a 2–3‰ NCE during the PETM event was likely caused by the thermal dissociation of oceanic CH₄ hydrate with a δ¹³C value of −60‰ (Table 1) (Dickens et al., 1995). Later, the rapid burning of accumulated Paleocene terrestrial organic carbon (Table 1) was proposed as an alternative to releasing gas hydrates to explain the decoupled C and S isotope variations during this interval with a simple box model of the C and S cycles (Kurtz et al., 2003). When seafloor CaCO₃ dissolution data during the PETM became available, the cGENIE model was used to match the data and suggested a −22‰ carbon source (~6,800 Pg) (Table 1), possibly a combination of methane, volcanic CO₂, organic matter-derived CO₂ (Panchuk et al., 2008).

However, the methane hypothesis (~1,500–3,000 Pg carbon with a δ¹³C value of −50‰) (Table 1) was shown to also reproduce the seafloor CaCO₃ dissolution curve and a ~−4‰ NCE in the LOSCAR model (Zeebe et al., 2009; Frieling et al., 2016). Similarly, a combined ~900–1,400 Pg of carbon as CO₂ and CH₄ (Table 1) was shown to drive the observed δ¹³C excursion and global warming of 3°C in a global C cycle box model coupled with an atmospheric CH₄ box (Carozza et al., 2011). However, the cGENIE model (Cui et al., 2011) required at least double the amount of CH₄ (Table 1) for isotopic mass balance since most of the released methane remains in the sediments, as shown in a thermohydraulic model (Minshull et al., 2016).

Seawater pH variations indicated from B isotope data were then used together with carbon isotopes to constrain the carbon source in the cGENIE model (Gutjahr et al., 2017), showing that a volcanic-dominant C source of −11‰ to −17‰ (isotopically heavier than previously considered) (Table 1) was likely related to the North Atlantic Igneous Province, although records of Hg enrichment in sediments suggest that volcanism did not provide all of the carbon input (Kender et al., 2021). An abrupt input of an isotopically-light carbon source (−35‰) combined with elevated volcanic outgassing (−6‰) could also satisfy the constraint of C isotopes, pH, and temperatures in the cGENIE model (Kirtland Turner, 2018). Recent box modeling of the C cycle suggested a carbon source of −22‰ (Table 1) based on new C isotope data from biomarkers (Cui et al., 2021a). In conclusion, the C composition and input amounts during the PETM event are still sources of debate (Table 1).

Simulations using various box models and 3D Earth system models have suggested possible environmental outcomes of the PETM event. Examples are listed here in order of modeling complexity from more simple to more complex. A simple box model of the Li cycle was applied to study the observed negative δ⁷Li excursion, suggesting enhanced continent weathering during the PETM event (Pogge Von Strandmann et al., 2021). Earth system models with intermediate complexity (represented by the cGENIE model) have proposed globally enhanced calcite burial (or deepening of the calcite compensation depth) (Penman et al., 2016), enhanced organic burial (Gutjahr et al., 2017), enhanced ocean stratification (Kirtland Turner et al., 2017), and ocean deoxygenation (Remmelzwaal et al., 2019) (Figure 6). Earth system models with higher complexity, such as CESM, have also been used to study the climate impact of the PETM event, revealing significant changes in hydroclimate such as extreme precipitation (Rush et al., 2021; Shields et al., 2021) and increased surface seawater temperature (Frieling et al., 2017) due to increased CO₂. Hydroclimate responses simulated in the Hadley Centre Atmosphere-Ocean General Circulation Model (HadCM3L), were suggested to be heterogeneous, meaning precipitation-evaporation increased in some regions but decreased in others (Carmichael et al., 2017). These impacts on climate or environment, particularly enhanced organic carbon and carbonate burial, acted as negative feedbacks on the increased CO₂ and temperature, ultimately contributing to the termination of the PETM event (e.g., Kurtz et al., 2003; Penman et al., 2016).

Short summaries and comparisons of hyperthermal events

Plentiful modeling studies, as well as proxy records, have suggested that the Mesozoic-Paleogene hyperthermal events were caused by vast emissions of ¹³C-depleted carbon into the atmosphere and ocean. The proposed composition of carbon emissions that triggered these hyperthermal events varies significantly for a specific hyperthermal event or among these hyperthermal events: δ¹³C values are in the range of −60 to −5‰ for the PTB event, of −60 to −6‰ for TOAE, of −70 to 0‰ for OAE 1a, −5‰ for OAE 2, of −60 to −11‰ for the PETM event, respectively (Table 1). Whereas their amount of carbon emissions is in the range of 900–450,000 Pg C (Table 1). We find that the light carbon source that drove the PETM event is believed to have huge contributions from organic carbon or methane (Dickens et al., 1995; Panchuk et al., 2008; Zeebe et al., 2009; Carozza et al., 2011; Cui et al., 2011; Gutjahr et al., 2017; Cui et al., 2021a), and the carbon source for OAE 2 is probably from a volcanic CO₂ source (Kuroda et al., 2007; Flögel et al., 2011). Besides the PETM and OAE 2, the assumed carbon sources for other hyperthermal events range from volcanic sources to organic matter and methane-derived carbon (e.g., Adloff et al., 2020). The shortest onset duration (8–23 kyr) for warming 5–8°C during the PETM event (Hu et al., 2020) likely explains the necessity of involving a carbon source with a δ¹³C value lighter than volcanic sources (< −5‰) for the PETM event. More detailed duration and warming magnitudes comparisons refer to Hu et al. (2020).

In addition, both models and proxy records have shown the biogeochemical consequences of warming, including but not limited to intensified continental weathering, ocean anoxia, and acidification. However, their exact magnitudes are still subject to debate. Enhanced continent weathering is an important and necessary precondition (Table 2) to respond to high CO₂ conditions and to replicate the isotope records (i.e., negative δ⁴⁴/⁴⁰Ca and δ⁷Li excursions) (e.g., Lechler et al., 2015; Sun et al., 2018), although the radiogenic isotopes (Sr and Os) might show different behaviors among hyperthermal events likely due to the different level of eruptions of LIPs and their weathering and different hypothermal intensity (Pogge Von Strandmann et al., 2013; Them et al., 2017b). Given different box models of different combined element cycles applied, the proposed enhanced weathering factor ranges from 1.5 to 18 times for these hyperthermal events (Table 2). Considering that these simple box models generally lack the constraints from the carbon cycle, combing them with Earth system box models or 3D Earth system models is promising to further quantify the weathering process during each hyperthermal event.

Modeling studies also suggest global warming could lead to ocean anoxia due to lower oxygen solubility, weaker ocean circulation (e.g., Penn et al., 2018), and likely higher surface productivity that produces more organic matter consuming oxygen in deeper seawater (e.g., Monteiro et al., 2012). However, quantifying the extent of ocean anoxia is still subject to modeling and proxy efforts since the exact extent of ocean anoxia is still unknown for all hyperthermal events. For instance, the assumed anoxic extent ranges significantly from ~1–2% to >90% of the global seafloor area (Montoya-Pino et al., 2010; Flögel et al., 2011; Owens et al., 2013) during OAE 2 (Table 3). Meanwhile, compared with OAE 2, the quantitative estimates of ocean anoxia for other hyperthermal events lack, which requires more proxy and modeling efforts.

Previous studies have classified these hyperthermal events into two types based on seawater carbon isotope characteristics (Hu et al., 2020). One type of hyperthermal event (PTB, TOAE, and PETM) generally shows an NCE pattern, and the other type (OAE 1a and OAE 2) exhibits PCE possibly with an abrupt NCE at its onset (Figure 1). However, the C isotope pattern for these hyperthermal events requires further global stratigraphic correlation work to verify. For example, whether an abrupt NCE exists at the onset of OAE 2 is uncertain. In contrast to an NCE record at the initial part of the PCE of OAE 2 from Bonarelli of Western Tethys (Kuroda et al., 2007), the C isotope profile from the Gongzha of Southern Tibet lacks an NCE at the onset of PCE (Hu et al., 2020). In addition, other seawater carbon variations prior to or after the hyperthermal event may also interfere with the identification of the C isotope variation pattern for each hyperthermal event, such as a PCE prior to the NCE of PETM (Hu et al., 2020; Wang et al., 2022). Future stratigraphic correlation works are expected to help address these uncertainties in identifying seawater C isotope patterns.

The hyperthermal events of PCE have been previously considered to be triggered by volcanism associated with LIPs in the deep-sea environment, while the NCE type was initiated by LIPs on land (Hu et al., 2020). However, this theoretical hypothesis requires further modeling verifications. As far as the current modeling work suggests, vast emissions of ¹³C-depleted carbon into the atmosphere and ocean would drive a (short-term) NCE in seawater (e.g., Kuroda et al., 2007), while the enhanced weathering input of nutrients and regeneration of P from the seafloor promoted surface productivity and organic carbon burial in anoxic ocean, ultimately resulting in a (relatively broad) PCE record (e.g., Flögel et al., 2011; Chen et al., 2022). Therefore, we propose that the different responses of seawater carbon isotope patterns among hyperthermal events likely result from interactions between the light carbon input and increased organic carbon burial. In other words, interactions between the light carbon emission pattern (i.e., carbon source, carbon emission rate, and total carbon emission) and the extent of primary productivity and organic carbon burial largely determine the seawater C isotope response pattern to the warming event.

To summarize, the development of simple box models, Earth system box models, and 3D Earth system models progressed in tandem with the increasing availability of proxy constraints, with a focus on the composition and amount of carbon emission that triggered these hyperthermal events (Table 1), and on the biogeochemical consequences of warming (Tables 2, 3). However, the exact amplitude of environmental impact is still debated and requires further modeling and proxy efforts. In addition, the mechanism that terminated the warming has been neglected in the previous modeling works.

Simulating hyperthermal events in long-term climate-biogeochemical models

To date, the models used to study hyperthermal events span from simple single-element box models, and Earth system box models (e.g., LOSCAR, COPSE), to complex 3D spatially-resolved models. These models have provided valuable insights into the drivers and potential impacts of these hyperthermal events, but they have limitations and uncertainties. The simple box models or Earth system box models do not capture the spatial resolution and full complexity of the Earth system and its feedback mechanisms, meaning that geological processes, particularly weathering processes, are very poorly represented within them (e.g., Goddéris et al., 2014). The weathering of LIPs during hyperthermal events, which is crucial in reducing the increased atmospheric CO₂ (e.g., Pogge Von Strandmann et al., 2013; Them et al., 2017b), has been chiefly modeled in box models, but they do not account for the spatial variations in temperature, runoff, and topography, which are essential factors influencing the weathering process.

On the other hand, the 3D Earth system models are computationally expensive and have limitations in terms of the simulated timescale – they cannot be integrated across the whole duration of most of the events discussed. A compromise between these approaches can be made with so-called ‘semi-spatial’ Earth system models. These models use data structures of steady-state 3D climate model simulations to represent continental surface processes in 3D while the ocean is represented by a box model. The first of these models to be developed was GEOCLIM (Donnadieu et al., 2006), which can run up to around 10 Myr – limited by the changing continental configuration. It has been applied to hyperthermal events by producing predictions over long >100 kyr timescales for the Deccan Traps warming at the K–Pg boundary, combining the LIP warming with the bolide impact event (Le Hir et al., 2020).

A new semi-spatial model is the Spatial-Continuous Integration ‘SCION’ model (Mills et al., 2021; Zhang et al., 2023). It uses the same climate model data structure as GEOCLIM but interpolates this over geological time to run continuously over the 539 million years of the Phanerozoic Eon. This allows for extensive model validation against multiple environmental datasets (atmospheric CO₂, atmospheric O₂, marine SO₄) and isotope proxies (δ¹³C, δ³⁴S, ⁸⁷Sr/⁸⁶Sr), and therefore provides a relatively robust background model setup for simulating a hyperthermal event (e.g., CO₂ level, marine phosphate level, marine productivity, background CO₂ degassing rate) that closely aligns with the Phanerozoic-scale geological record.

Models like SCION and GEOCLIM are currently under-utilized for exploring hyperthermal events and could help tie together the long-term but more poorly-constrained experiments in simple box models or Earth system box models, with the detailed but shorter-term experiments in 3D models. For example, SCION has been used to test the interplay between CO₂ degassing and LIP weathering in controlling climate over the Columbia River Basalt emplacement (Longman et al., 2022). These types of simulations that include long-term continental and marine processes may help us better understand these events – particularly concerning what terminates a hyperthermal event. Event termination is currently under-explored compared to the initiation, given that 3D Earth system models typically cannot run to the full duration of the event, which could last up to millions of years (e.g., for the PTB-Early Triassic hyperthermal; Sun et al., 2018; Wu et al., 2021).

Conclusions

Numerous simulation studies have been conducted to investigate the mechanisms of hyperthermal events. When summarizing the existing modeling work on the Mesozoic-Paleogene hyperthermal events, it can be observed that: 1) the developments of simple box models, Earth system box models, and 3D Earth system models have progressed in tandem with the increasing availability of proxy constraints, 2) the biogeochemical consequences of warming, including but not limited to intensified continental weathering, ocean anoxia, and acidification, have been widely acknowledged in both models and proxy records, but their exact magnitudes and drivers are still subject to debate and require further modeling and proxy efforts, 3) determining the composition and amount of carbon emissions that triggered these hyperthermal events also demands additional modeling and proxy constraints. Semi-spatial models like SCION and GEOCLIM may help to bridge the gap between complex 3D representations of continental processes and long timescales of geological records, potentially offering fresh insights and perspectives on the drivers of, and biogeochemical responses to, hyperthermal events.

Author contributions

YZ: Conceptualization; Investigation; Writing – original draft; Writing – review & editing; BM: Conceptualization; Funding acquisition; Supervision; Writing – review & editing; TH: Conceptualization; Supervision; Writing – review & editing; XH: Supervision; Writing – review & editing; MZ: Conceptualization; Funding acquisition; Supervision. All authors contributed to the article and approved the submitted version.

Funding

This research was financially supported by the National Natural Science Foundation of China (41888101 and 41921002), the China Postdoctoral Science Foundation (2023M733592), and the Jiangsu Funding Program for Excellent Postdoctoral Talent. BM’s contribution was supported by the UK Natural Environment Research Council (NE/S009663/1 and NE/X011208/1).

Acknowledgments

This manuscript contributes to the Integrated Understanding of the Early Jurassic Earth System and Timescale (JET) project and the IGCP 739 project.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer SJ declared a shared affiliation with the author TH to the handling editor at time of review.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References