# Assessing the importance of thermogenic degassing from the Karoo Large Igneous Province (LIP) in driving Toarcian carbon cycle perturbations

Thea Heimdal, Yves Goddéris, Morgan Jones, Henrik Svensen

► To cite this version:

Thea Heimdal, Yves Goddéris, Morgan Jones, Henrik Svensen. Assessing the importance of thermogenic degassing from the Karoo Large Igneous Province (LIP) in driving Toarcian carbon cycle perturbations. Nature Communications, 2021, 12 (1), 10.1038/s41467-021-26467-6. http://hal.science/hal-03412325

HAL Id: hal-03412325

https://hal.science/hal-03412325

Submitted on 5 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ARTICLE

https://doi.org/10.1038/s41467-021-26467-6

OPEN

<figure><img src="image_2.png" /></figure>

Check for updates

# Assessing the importance of thermogenic degassing from the Karoo Large Igneous Province (LIP) in driving Toarcian carbon cycle perturbations

Thea H. Heimdal<fnref n="1" />✉, Yves Goddéris², Morgan T. Jones<fnref n="1" /> & Henrik H. Svensen<fnref n="1" />

The emplacement of the Karoo Large Igneous Province (LIP) occurred synchronously with the Toarcian crisis (ca. 183 Ma), which is characterized by major carbon cycle perturbations. A marked increase in the atmospheric concentration of CO₂ (pCO₂) attests to significant input of carbon, while negative carbon isotope excursions (CIEs) in marine and terrestrial records suggest the involvement of a ¹²C-enriched source. Here we explore the effects of pulsed carbon release from the Karoo LIP on atmospheric pCO₂ and δ¹³C of marine sediments, using the GEOCLIM carbon cycle model. We show that a total of 20,500 Gt C replicates the Toarcian pCO₂ and δ¹³C proxy data, and that thermogenic carbon (δ¹³C of -36 ‰) represents a plausible source for the observed negative CIEs. Importantly, an extremely isotopically depleted carbon source, such as methane clathrates, is not required in order to replicate the negative CIEs. Although exact values of individual degassing pulses represent estimates, we consider our emission scenario realistic as it incorporates the available geological knowledge of the Karoo LIP and a representative framework for Earth system processes during the Toarcian.

<footnote n="1">¹Centre for Earth Evolution and Dynamics (CEED), University of Oslo, Oslo, Norway. ²Géosciences-Environnement Toulouse, CNRS-Université Paul Sabatier, Toulouse, France. ✉email: t.h.heimdal@geo.uio.no</footnote>
The Toarcian crisis (ca. 183 Ma) is characterized by extinctions in both the marine and terrestrial realm, localized ocean anoxia (The Toarcian Ocean Anoxic Event; T-OAE), global warming and major carbon cycle perturbations¹⁻¹⁵. The carbon cycle disturbances are evidenced by increases in the atmospheric concentration of CO₂ (pCO₂;¹⁻³) and negative carbon isotope excursions (CIEs) recorded in both carbonate and organic matter in marine and terrestrial archives⁴⁻¹⁴. A prominent negative CIE marks the onset of the T-OAE, and varies in magnitude in both carbonate and organic matter (δ¹³Ccarb ≈ -0.6 to -6.0%; δ¹³Corg ≈ -0.8 to -8.6%; see review in ref. ¹⁵). High-resolution δ¹³C curves from multiple sections show a stepwise character of the excursion, with several abrupt negative shifts making up the full body of the CIE³,⁴,⁶⁻⁹,¹¹,¹²,¹⁴⁻¹⁷.

High-precision U-Pb geochronology demonstrates that the Karoo LIP in South Africa and Lesotho (Fig. 1) was emplaced synchronously with the T-OAE negative CIE (Fig. 2; ref. ¹⁸⁻²¹). This has led to the hypothesis that Karoo-derived emissions were responsible for the carbon cycle perturbations (e.g., ¹,⁷,¹⁸⁻²¹); however, the negative shift in δ¹³C is of such magnitude that mantle-derived carbon alone is unlikely to have been the source²². Alternative suggestions include marine methane clathrate dissociation (e.g., ⁴,⁸,²²), permafrost melting²³, terrestrial organic matter²⁴,²⁵, or thermogenic carbon release generated by sill emplacement into organic-rich sedimentary rocks of the Karoo Basin (e.g., ¹,¹⁹,²²,²⁶).

Volcanic basins, such as the Karoo, are characterized by a pulsed emplacement of sills (e.g.,²⁷), leading to pulsed contact metamorphism and fluid release. Notably, a stepwise decrease in δ¹³C is suggestive of isotopically light carbon supplied in pulses to the ocean-atmosphere system²⁸. Field investigations, borehole studies, and numerical modeling suggest thermogenic carbon generation on the order of several thousand gigatons within the Karoo basin²⁶,²⁹⁻³¹. The presence of thousands of explosion pipes and hydrothermal vent complexes rooted in the contact metamorphic aureoles represent a plausible release mechanism, and attests to rapid emissions of the thermogenic gases²⁶,³². Available high-precision U-Pb ages of the sills and thermogenic emission estimates from the Karoo Basin allows to identify detailed and realistic emission scenarios from the Karoo LIP; however, these have not yet been tested by carbon cycle modeling.

Here, we test such scenarios using the GEOCLIM model³³,³⁴, in order to explore whether pulsed carbon release from the Karoo LIP could explain the full extent of observed Toarcian pCO₂ and δ¹³C changes. The emission scenarios are grounded by U-Pb geochronology, detailed evolution of the Karoo LIP, and estimates of both the mantle-derived and thermogenic carbon from volcanic and contact metamorphic degassing, respectively.

## Results

In order to investigate the consequences of carbon release from the Karoo LIP, it is essential to constrain the timing between the negative CIEs, pCO₂ excursion, and Karoo activity. Following the correlations used here (see Supplementary Note 1;¹⁸,¹⁹,³⁵), the T-OAE lasted for 500 kyr, and overlaps with the full extent of the carbon isotope perturbations (Fig. 2). This corresponds to the onset of the negative CIEs, as well as peak negative values, and the subsequent interval of δ¹³C increase eventually reaching background levels. Reconstructed pCO₂ curves from the Danish Basin show an increase in pCO₂ of up to ~2000 ppm during the T-OAE (including uncertainties; ref. ¹). Following our correlations, the pCO₂ peak coincides with the peak negative δ¹³C values (Fig. 2). The pCO₂ records do not include absolute time constraints, so the correlation of pCO₂ data with the Karoo activity and various

carbon isotope records was based on the extent of the T-OAE (see Supplementary Note 1). Available U-Pb ages of Karoo intrusives demonstrate that the Karoo igneous activity overlaps with the full extent of the Toarcian CIEs and by inference the pCO₂ peak (Fig. 2; ref. ¹⁸⁻²¹).

The Karoo intrusives postdate a ~1-2‰ negative CIE that occurs around the Pliensbachian-Toarcian (PI-T) boundary in some sections (e.g.;⁵,³⁶, Fig. 2), excluding a temporal link between Karoo igneous activity and this earlier disruption of the carbon cycle³⁷. As this CIE is not systematically recorded, and there is no evidence for pCO₂ increase or Karoo activity at this time, it is not considered in this study. Furthermore, we stress that there are alternative correlations and interpretations for the Toarcian carbon cycle perturbations that differ from those used here, particularly with respect to the duration of the full extent of the T-OAE CIE, which varies from 120 to 2400 kyr (Supplementary Note 1;⁴,⁷,⁸,¹⁸,¹⁹,³⁵,³⁸⁻⁴³).

The total model run is set to 500 kyr, which corresponds to the total duration of the T-OAE (following the correlations used here; Fig. 2; Supplementary Note 1). Both volcanic and thermogenic processes associated with the Karoo LIP activity are accounted for in the emission scenario. Seven thermogenic carbon pulses (δ¹³C of -36‰; Supplementary Note 1) are released between model time (t) = 30 and 262 kyr (pulses #2-8; Table 1; Supplementary Fig. 1), which is constrained by the timing of negative CIEs. The thermogenic carbon pulses vary in magnitude, between 300 and 3800 Gt. The total magnitude of released mantle-derived carbon (with δ¹³C of -5‰) is set to 8000 Gt, following average estimates for the Karoo LIP (e.g., ref. ²²; Supplementary Note 1). 4500 Gt mantle-derived carbon is released at t = 160 (pulse #6; Table 1), while the remaining 3500 Gt C is released as a long-lasting and continuous pulse between t = 0 and 262 kyr (pulse #1; Table 1; Supplementary Fig. 1; Supplementary Note 1). Note that carbon pulse #6 represents a mix of both mantle-derived and thermogenic carbon, with a total magnitude of 7000 Gt and a δ¹³C value of -16‰ (Table 1; Supplementary Note 1). Except for the continuous mantle-derived carbon pulse (pulse #1), the carbon releases have durations of 3-22 kyr, which is constrained by the duration of observed carbon cycle perturbations targeted, and in agreement with modeling studies suggesting that the generation of thermogenic fluids occurs on timescales of 100's to 1000's years²⁹,⁴⁴.

The pulsed release of a total of 20,500 Gt C leads to a stepwise increase in pCO₂ from ~650 ppm up to a maximum of ~1700 ppm. The modeled pCO₂ curve plots generally within the range of observed data, with the exception of the fourth peak (at t = 111 kyr), which overestimates the pCO₂ (Fig. 3). The modeled δ¹³Ccarb curve decreases in seven steps (between ~0.3 and ~1.8‰), with a total negative shift (δ¹³Cmax - δ¹³CpreT-OAE) of ~-4‰. Following the δ¹³Ccarb data, the modeled δ¹³Corg curve also decreases in seven steps (between ~0.5 and ~3.7‰), but with larger magnitudes, as well as a larger total negative shift of ~-6‰. Both the modeled δ¹³Ccarb and δ¹³Corg curves generally plot within the range of observed data considering both the magnitude and shape of the negative shift (Fig. 3). There are however some exceptions: 1) the modeled δ¹³Ccarb curve predicts slightly lower δ¹³C values at the end of the model run during the recovery interval of δ¹³C increase, and 2) the modeled δ¹³Corg curve does not replicate all of the short-lived negative CIEs making up the prolonged interval of negative δ¹³C values following peak negative values. Note that all observed δ¹³C curves are from shallow marine stratigraphic sections, and is accordingly compared to model δ¹³C data representing a shallow ocean epicontinental sea or shelf environment (i.e., epicontinental surface reservoir boxes; see Methods).
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 1 Location and cross section of the Karoo Basin.** Top: a plate reconstruction at 180 Ma (modified from ref. ⁶⁰; Fig. 12.2) showing the known and inferred extent of the Karoo Large Igneous Province (LIP)²¹,⁶¹. Bottom: a detailed map of the Karoo Basin stratigraphy in South Africa and Lesotho, with a schematic cross section across the province (modified from ref. ⁴⁷). Photo A shows deep laminar sills of the Ecca Group in western South Africa, while photo B shows the continental flood basalt lavas exposed in Lesotho (positions are labelled as A and B on the map and cross section).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Correlation between Toarcian carbon cycle perturbations and the timing of Karoo intrusions. Timeline of Karoo activity, the Toarcian Ocean Anoxic Event (T-OAE; green box), carbon isotope curves from shallow marine stratigraphic sections and pCO₂ data (including uncertainties) after ref. ¹⁸,¹⁹,³⁵. The age of the Pliensbachian-Toarcian (Pl-To) boundary (orange line) is from ref. ¹⁹. U-Pb ages of Karoo intrusives are from ref. ¹⁸⁻²¹. U-Pb color codes refer to sill emplacement levels in the Karoo Basin, and corresponds to the colors used in the map and cross section in Fig. 1. Numbers in boxes represent the total negative shift (δ¹³Cmax − δ¹³CpreT-OAE) in δ¹³Corg and δ¹³Carb.</figcaption></figure>

<table><caption>**Table 1** Overview and input values for Karoo emission scenario.</caption><thead><tr><th>Carbon pulse no.</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th></tr></thead><tbody><tr><td>Magnitude (Gt)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Volcanic C</td><td>3500</td><td>N/A</td><td>N/A</td><td>N/A</td><td>N/A</td><td>4500</td><td>N/A</td><td>N/A</td></tr><tr><td>Thermogenic C</td><td>N/A</td><td>300</td><td>500</td><td>2000</td><td>3800</td><td>2500</td><td>1700</td><td>1700</td></tr><tr><td>Total</td><td>3500</td><td>300</td><td>500</td><td>2000</td><td>3800</td><td>7000</td><td>1700</td><td>1700</td></tr><tr><td>δ¹³C (‰)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Volcanic C</td><td>-5</td><td>N/A</td><td>N/A</td><td>N/A</td><td>N/A</td><td>-5</td><td>N/A</td><td>N/A</td></tr><tr><td>Thermogenic C</td><td>N/A</td><td>-36</td><td>-36</td><td>-36</td><td>-36</td><td>-36</td><td>-36</td><td>-36</td></tr><tr><td>Total</td><td>-5</td><td>-36</td><td>-36</td><td>-36</td><td>-36</td><td>-16</td><td>-36</td><td>-36</td></tr><tr><td>Model time t (ky)</td><td>0</td><td>30</td><td>62</td><td>82</td><td>111</td><td>160</td><td>210</td><td>254</td></tr><tr><td>Duration (ky)</td><td>262</td><td>3</td><td>3</td><td>8</td><td>19</td><td>22</td><td>11</td><td>8</td></tr><tr><td>Flux (Gt C/y)</td><td>N/A</td><td>0,1</td><td>0,2</td><td>0,3</td><td>0,2</td><td>0,3</td><td>0,2</td><td>0,2</td></tr></tbody></table>

N/A, not applicable

Discussion

The emission scenario presented here replicates an overall global trend in δ¹³C, including several negative steps, and a total negative shift (δ¹³Cmax − δ¹³CpreT-OAE) of ~−4‰ for δ¹³Ccarb and ~−6‰ for δ¹³Corg. This is comparable to observed Toarcian carbonate δ¹³C records showing generally lower CIE magnitudes compared to those of organic matter, with mean total CIE values of ~−3 vs. ~−5‰, respectively¹⁵. Globally, the total magnitude of the T-OAE CIE ranges considerably (up to ~−6 for δ¹³Ccarb and ~−9‰ for δ¹³Corg; ref. ¹⁵), however, the largest excursions are likely exaggerated due to local effects¹⁵,⁴³. Furthermore, all high-resolution Toarcian carbon isotope curves show a stepwise character of the T-OAE excursion, and a series of at least five prominent negative shifts is recognized in numerous carbonate and organic matter Toarcian sections (e.g., Kaszewy, Mochras, Yorkshire, Peniche and Luxembourg;¹²,⁴⁵). δ¹³C curves showing one main long-lived T-OAE CIE, or a limited number of negative shifts, have likely failed to record the full extent of the stepped descent due to low sample or stratigraphic resolution. We chose not to take the observed prolonged interval of negative δ¹³Corg values following peak negative values into account, because it is uncertain whether or not this represents a global trend. It has been suggested that, due to improved preservation of marine organic matter under oxygen-depleted conditions, artificial Toarcian δ¹³Corg excursions have been produced⁴⁰. This could potentially explain why there is an apparent prolonged interval of δ¹³Corg CIEs, while δ¹³Carb values instead return to background

levels (which is also consistent with the shape of the pCO₂ curve). Although Toarcian δ¹³C curves vary with respect to the magnitude, number, timing, and durations of the negative CIEs, we consider our modeled δ¹³Ccarb and δ¹³Corg curves to be representative for the overall global trend of Toarcian carbon cycle changes.

While the modeled δ¹³C data is compared to several observed δ¹³C curves, available reconstructed Toarcian pCO₂ curves are scarce. As seen in Fig. 3, the modeled pCO₂ curve generally plots within the range of observed Toarcian pCO₂ data from ref. ¹, but potential target pCO₂ values represent quite a wide range due to data points with large error bars. Consequently, the magnitude of released carbon could vary greatly and still replicate the observed pCO₂ data within the range of uncertainties. However, by respecting both observed pCO₂ and δ¹³C curves simultaneously, the possible magnitude of carbon release is narrowed down considerably. Furthermore, based on Toarcian seawater pH (ref. ²) and Δ13C data (δ¹³Cterrigenic − δ¹³Cmarine; ref. ³), a pCO₂ increase from ~800–850 ppm to ~1750–1800 ppm has been calculated (when assuming a release of thermogenic carbon with δ¹³C of −40 to −30‰ for the latter; ref. ³). As seen in Fig. 3, an increase of this magnitude (i.e., 800–850 to 1750–1800 ppm) is in very close agreement with the modeled pCO₂ data. The overestimation of pCO₂ at t = 111 kyr shown by the modeled curve (Fig. 3) is likely the result of the low sample resolution of the observed data, as it is very unlikely that a shift in δ¹³C of up to 4‰ is not reflected by a change in pCO₂. Therefore, we consider the modeled pCO₂ curve presented here to
<figure><img src="image_6.png" /><figcaption>**Fig. 3 Modeled pCO₂ and δ¹³C curves compared to observed Toarcian data.** Model response (pink dashed lines) of δ¹³Corg and δ¹³CCarb of the shallow ocean (GEOCLIM epicontinental surface reservoir boxes) and atmospheric pCO₂ to a Karoo emission scenario including eight pulses of carbon release (see Table 1; Supplementary Fig. 1). Gray and blue outlines represent the range of observed Toarcian δ¹³CCarb<sup>4,5,9</sup>, δ¹³Corg<sup>8,16,36</sup>, and pCO₂ data¹, which corresponds to the curves presented in Fig. 2.</figcaption></figure>

be realistic as it reflects Toarcian pCO₂ estimates constrained by
several methods.

According to our emission scenario, a total release of 20,500 Gt carbon replicates the observed Toarcian pCO₂ and δ¹³C data. This number includes 12,500 Gt thermogenic carbon, which is within the range of thermogenic carbon generation estimations from the Karoo Basin (~7500–23,000 Gt C; ref. ²⁶,²⁹–³¹). The ratio of trapped versus released thermogenic gases from volcanic basins is poorly constrained, although 100% degassing from intrusive LIP components is considered unrealistic (e.g., ⁴⁶). Consequently, the magnitude of released carbon must be lower than the magnitude of carbon generated. The lowermost thermogenic carbon generation estimates are lower than the release value constrained by our modeling. However, these numbers are likely significantly underestimated as they have been upscaled based on a much smaller area of affected host rock compared to that of ref. ³¹, and do not include devolatilization of Ecca coals or Beaufort, Dwyka and Stormberg organic-bearing rocks in the Karoo strata (see Fig. 1). According to ref. ³¹, an estimated 72 % of the generated carbon gases in the Karoo Basin were released, which corresponds to ~17,000 Gt C. We therefore consider 12,500 Gt to represent a realistic value of cumulative thermogenic carbon release from the Karoo Basin.

Considering the TOC content and sill distribution in the Karoo Basin, sill emplacement into Ecca and Beaufort shales volatilized significantly more carbon-bearing gases compared to sill emplacement at other stratigraphic levels. Within the Ecca and Beaufort series, sills are abundant and can be traced almost

continuously with depth, as well as between the basin margins (Fig. 1; ref. ⁴⁷). In contrast, sills are rare in the Stormberg Group, while devolatilization of Dwyka sedimentary rocks was likely restricted to the upper parts of the series based on the distribution of the sills. The Ecca shales have the highest TOC contents, reaching locally as much as 18 wt.%, while maximum values within the Stormberg, Beaufort and Dwyka Group sedimentary rocks range from 3 to 6 wt.%⁴⁸⁻⁵³. Although the TOC contents are overall lower compared to the Ecca Group, the presence of closely spaced large sills (up to 270 m thick; ref. ⁴⁷) in the Beaufort Group with TOC up to ~6 wt.%⁵² suggests significant gas generation potential at these stratigraphic levels as well. Sill thicknesses and vertical spacing are important factors influencing the gas generation potential; gas generation can increase by up to ~35 % for more closely spaced sills, compared to separate sills emplaced into the same host-rock²⁷.

According to our emission scenario, the magnitudes of the individual thermogenic carbon pulses vary throughout the model run (Table 1), which is constrained by the magnitudes of the observed carbon cycle perturbations targeted. The available U-Pb data cannot be used to directly assess the magnitude of the individual degassing pulses as 1) there is no systematic correlation of sill age with stratigraphic emplacement level, and 2) they lack the sufficient resolution (large error bars) to match that of the stepped δ¹³C curves. Importantly, however, the U-Pb geochronology links sill emplacements within organic-rich sedimentary rocks (i.e., Ecca and Beaufort) to the entire interval of Toarcian carbon cycle
perturbations (Fig. 2). Furthermore, the proposed carbon emission fluxes of 0.2 to 0.3 Gt C yr⁻¹ (Table 1) corresponds to those estimated for other LIPs (e.g., NAIP and CAMP; 0.1–0.5 Gt C yr⁻¹; ref. 54,55). We therefore consider the magnitudes of the individual pulsed thermogenic carbon releases to represent realistic values.

Realistic temporal and volumetric estimates of Karoo LIP volatile emissions can explain the carbon cycle perturbations observed in the Toarcian. Pulsed sill intrusions and contact metamorphism of organic matter-bearing sedimentary rocks in the Karoo Basin represents a significant source of ¹³C-depleted carbon, which can explain the observed negative CIEs. As the lowest δ¹³C value of the released thermogenic carbon in our model is −36‰, an extremely isotopically depleted carbon source (e.g., methane clathrates; −60‰; ref. 22) is not a required component (see also Supplementary Fig. 2). High-precision U-Pb geochronology links sill emplacements within the organic-rich strata of the Karoo Basin to the negative CIEs and the abrupt pCO₂ peak. Although the exact values for the different degassing pulses represent estimates, we consider our emission scenario realistic as it incorporates the available geological knowledge of the Karoo LIP and a representative framework for Earth system processes during the Toarcian.

## Methods

**Model description.** GEOCLIM couples the 10-box model COMBINE⁵⁶, describing different geochemical cycles (e.g., C, P, O₂, alkalinity, Sr) and their associated isotopic cycles (δ¹³C, ⁸⁷Sr/⁸⁶Sr), to the 3D FOAM general circulation climate model⁵⁷. COMBINE consists of one atmospheric box and nine ocean boxes, including 1) an epicontinental surface reservoir divided in two boxes (surface box up to 100 m depth; deep box from 100 to 200 m depth), 2) three mid-latitude ocean boxes (surface, thermocline, and deep), and 3) two polar oceans both divided in two boxes (surface and deep). A full description of the GEOCLIM model can be found in refs. ³⁴ and ⁵⁶.

Boundary conditions for the simulations presented here include a Jurassic (180 Ma) paleogeography, which is derived from a synthesis of paleomagnetic data, hot spot tracks and geological constraints³⁴,⁵⁸. The background volcanic degassing rate was set to $4 \times 10^{12}$ mol/yr so that the pCO₂ equals ~ 650 ppm at t = 0 (prior to any Karoo degassing), which represents an average value of observed pCO₂ data (including uncertainties; ref. 1).

GEOCLIM calculates the δ¹³C of the various carbon species (e.g., H₂CO₃, POC; particulate organic carbon) in each reservoir/box (e.g., epicontinental surface box). The isotopic fractionation between H₂CO₃* (H₂CO₃ + CO₂; Supplementary Note 1) and POC is set to depend on the calculated H₂CO₃* concentration at each timestep for each oceanic reservoir. As shown in Supplementary Fig. 3a, the fractionation between H₂CO₃* and the POC increases during the interval of carbon injections, but then decreases towards “background values”. A decrease in fractionation reduces the difference between the δ¹³C of carbonate and POC, which allows the δ¹³Cₚc to increase faster compared to δ¹³Cₚc,⁵⁹. This could potentially explain the apparent decoupling between the δ¹³Cₚc and δ¹³Cₚcarb modeled curves (Fig. 3; Supplementary Fig. 3b); while the modeled δ¹³Cₚc curve recovers faster and reach observed values during the recovery interval at the end of the model run, the modeled δ¹³Cₚcarb curve predicts slightly lower δ¹³C values compared to observed values.

## Data availability

All data used for this paper is available in the main text and in the Supplementary Information (Supplementary Note 1, Supplementary Figs. 1–3 and Supplementary References).

## Code availability

The GEOCLIM source code is available upon request to the first author.

Received: 7 April 2021; Accepted: 24 September 2021;
Published online: 28 October 2021

## References

1. McElwain, J. C., Wade-Murphy, J. & Hesselbo, S. P. Changes in carbon dioxide during an oceanic anoxic event linked to intrusion into Gondwana coals. *Nature* **435**, 479–482 (2005).

2. Müller, T. et al. Ocean acidification during the early Toarcian extinction event: Evidence from boron isotopes in brachiopods. *Geology* **48**, 1184–1188 (2020).

3. Rubesam, W., Rolid, M. & Schwark, L. δ¹³C of terrestrial vegetation records Toarcian CO₂ and climate gradients. *Sci. Rep.* **10**, 117 (2020).

4. Hesselbo, S. P. et al. Massive dissociation of gas hydrate during a Jurassic oceanic anoxic event. *Nature* **406**, 392–395 (2000).

5. Hesselbo, S. P., Jenkyns, H., Duarte, L. & Oliveira, L. C. Carbon-isotope record of the Early Jurassic (Toarcian) Oceanic Anoxic Event from fossil wood and marine carbonate (Lusitanian Basin, Portugal). *Earth Planet. Sci. Lett.* **253**, 455–470 (2007).

6. Cohen, A. S., Coe, A. L., Harding, S. M. & Schwark, L. Osmium isotope evidence for the regulation of atmospheric CO₂ by continental weathering. *Geology* **32**, 157–160 (2004).

7. Cohen, A. S., Coe, A. L. & Kemp, D. B. The Late Palaeocene-Early Eocene and Toarcian (Early Jurassic) carbon isotope excursions: a comparison of their time scales, associated environmental changes, causes and consequences. *J. Geol. Soc. Sci.* **164**, 1093–1108 (2007).

8. Kemp, D. B., Coe, A. L., Cohen, A. S. & Schwark, L. Astronomical pacing of methane release in the Early Jurassic period. *Nature* **437**, 396–399 (2005).

9. Hermoso, M., Callonnec, L. L., Minoletti, F., Renard, M. & Hesselbo, S. P. Expression of the Early Toarcian negative carbon-isotope excursion in separated carbonate microfractions (Jurassic, Paris Basin). *Earth Planet. Sci. Lett.* **277**, 194–203 (2009).

10. Hermoso, M. et al. Dynamics of a stepped carbon-isotope excursion: ultra high-resolution study of Early Toarcian environmental change. *Earth Planet. Sci. Lett.* **319–320**, 45–54 (2012).

11. Sabatino, N. et al. Carbon-isotope records of the Early Jurassic (Toarcian) oceanic anoxic event from the Valdoria (Umbria–Marche Apennines) and Monte Mangart (Ju-lian Alps) sections: palaeoceanographic and stratigraphic implications. *Sedimentology* **56**, 1307–1328 (2009).

12. Hesselbo, S. P. & Pieńkowski, G. Stepwise atmospheric carbon-isotope excursion during the Toarcian Oceanic Anoxic Event (Early Jurassic, Polish Basin). *Earth Planet. Sci. Lett.* **301**, 365–372 (2011).

13. Léonide, P. et al. Drowning of a carbonate platform as a precursor stage of the Early Toarcian global anoxic event (Southern Province sub-Basin, South-east France). *Sedimentology* **59**, 156–184 (2012).

14. Pittet, B., Suan, G., Lenoir, F., Duarte, L. V. & Mattioli, E. Carbon isotope evidence for sedimentary discontinuities in the lower Toarcian of the Lusitanian Basin (Portugal): sea level change at the onset of the Oceanic Anoxic Event. *Sediment. Geol.* **303**, 1–14 (2014).

15. Remírez, M. N. & Algeo, T. J. Carbon-cycle changes during the Toarcian (Early Jurassic) and implications for regional versus global drivers of the Toarcian oceanic anoxic event. *Earth Sci. Rev.* **209**, 103283 (2020).

16. Xu, W. et al. Evolution of the Toarcian (Early Jurassic) carbon-cycle and global climatic controls on local sedimentary processes (Cardigan Bay Basin, UK). *Earth Planet. Sci. Lett.* **484**, 396–411 (2018).

17. Fantasia, A. et al. Global versus local processes during the Pliensbachian–Toarcian transition at the Peniche GSSP, Portugal: a multi-proxy record. *Earth-Sci. Rev.* **198**, 102932 (2019).

18. Sell, B. et al. Evaluating the temporal link between the Karoo LIP and climatic–biologic events of the Toarcian Stage with high-precision U–Pb geochronology. *Earth Planet. Sci. Lett.* **408**, 48–56 (2014).

19. Burgess, S. D., Bowring, S. A., Fleming, T. H. & Elliot, D. H. High-precision geochronology links the Ferrar large igneous province with early-Jurassic ocean anoxia and biotic crisis. *Earth Planet. Sci. Lett.* **415**, 90–99 (2015).

20. Corfu, F., Svensen, H. & Mazzini, A. Comment to paper: Evaluating the temporal link between the Karoo LIP and climatic–biologic events of the Toarcian Stage with high-precision U-Pb geochronology by Bryan Sell. *Earth Planetary Sci. Lett.* **434**, 349–352 (2016). Maria Ovtcharova, Jean Guex, Annachiara Bartolini, Fred Jourdan, Jorge E. Spangenberg, Jean-Claude Vicente, Urs Schaltegger in Earth and Planetary Science Letters **408** (2014) 48–56.

21. Greber, N. D., et al. New high precision U-Pb ages and Hf isotope data from the Karoo large igneous province: implications for pulsed magmatism and early Toarcian environmental perturbations. *Results Geochem. 1*, 100005 (2020).

22. Beerling, D. J. & Brentnall, S. J. Numerical evaluation of mechanisms driving Early Jurassic changes in global carbon cycling. *Geology* **35**, 247–250 (2007).

23. Rubesam, W., Mayer, B. & Schwark, L. Cryosphere carbon dynamics control early Toarcian global warming and sea level evolution. *Glob. Planet. Change* **172**, 440–453 (2019).

24. Pieńkowski, G., Hodbod, M. & Ullmann, C. V. Fungal decomposition of terrestrial organic matter accelerated Early Jurassic climate warming. *Sci. Rep.* **6**, 31930 (2016).

25. Them, T. R. II et al. High-resolution carbon isotope records of the Toarcian Oceanic Anoxic Event (Early Jurassic) from North America and implications for the global drivers of the Toarcian carbon cycle. *Earth Planet. Sci. Lett.* **459**, 118–126 (2017).

26. Svensen, H. et al. Hydrothermal venting of greenhouse gases triggering Early Jurassic global warming. *Earth Planet. Sci. Lett.* **256**, 554–566 (2007).

27. Aarnes, I., Svensen, H., Polteau, S. & Planke, S. Contact meatmorphic devolatilization of shales in the Karoo Basin, South Africa, and the effects of multiple sill intrusions. *Chem. Geol.* **281**, 181–194 (2011).
28. Jenkyns, H. C. Geochemistry of oceanic anoxic events. *Geochemistry, Geophysics, Geosystems* **11**, GC002788 (2010).

29. Aarnes, I., Svensen, H., Connolly, J. A. & Podladchikov, Y. Y. How contact metamorphism can trigger global climate changes: Modeling gas generation around igneous sills in sedimentary basins. *Geochimica et. Cosmochimica Acta* **74**, 7179–7195 (2010).

30. Iyer, K., Rüpke, L. & Galerne, C. Y. Modeling fluid flow in sedimentary basins with sill intrusions: Implications for hydrothermal venting and climate change. *Geochem., Geophysics, Geosystems* **14**, 5244–5263 (2013).

31. Galerne, C. Y. & Hasenclever, J. Distinct degassing pulses during magma invasion in the stratified Karoo Basin – new insights from hydrothermal fluid flow modeling. *Geochem. Geophysics, Geosystems* **20**, 2955–2984 (2019).

32. Svensen, H. H. et al. Understanding thermogenic degassing in Large Igneous Provinces: Inferences from the geological and statistical characteristics of breccia pipes in the western parts of the Karoo Basin. Geological Society of America Books, Special Paper, 544 (2020).

33. Donnadieu, Y. et al. A mechanism for brief glacial episodes in the Mesozoic greenhouse. *Paleoceanogr. Paleoclimatol.* **26**, PA002100 (2011).

34. Goddéris, Y., Donnadieu, Y., Le Hir, G., Lefebvre, V. & Nardin, E. The role of palaeogeography in the Phanerozoic history of atmospheric CO2 and climate. *Earth-Sci. Rev.* **128**, 122–138 (2014).

35. Bouilila, S. et al. Astronomical calibration of the Toarcian Stage: implications for sequence stratigraphy and duration of the early Toarcian OAE. *Earth Planet. Sci. Lett.* **386**, 98–111 (2014).

36. Caruthers, A. H., Grocke, D. R. & Smith, P. L. The significance of an Early Jurassic (Toarcian) carbon-isotope excursion in Haida Gwaii (Queen Charlotte Islands), British Columbia, Canada. *Earth Planet. Sci. Lett.* **307**, 19–26 (2011).

37. De Lena, L. F. et al. The driving mechanisms of the carbon cycle perturbations in the late Pliensbachian (Early Jurassic). *Sci. Rep.* **9**, 1–12 (2019).

38. McArthur, J. M., Donovan, T. D., Thirlwall, M. F., Fouke, B. M. & Mattey, D. Strontium isotope profile of the early Toarcian (Jurassic) oceanic anoxic event, the duration the ammonite biozones, and belemnite palaeotemperatures. *Earth Planet. Sci. Lett.* **179**, 269–285 (2000).

39. Kemp, D. B., Coe, A. L., Cohen, A. S. & Weedon, G. P. Astronomical forcing and chronology of the early Toarcian (Early Jurassic) oceanic anoxic event in Yorkshire, UK. *Paleoceanography* **26**, 1–17 (2011).

40. Clémence, M. E. et al. Micropalaeontologic, geochemical and cyclostratigraphic approach for the timing of the early Toarcian oceanic anoxic event in the Paris Basin (GPF-Sancerre borehole). *Volumina Jurass.* **4**, 154–156 (2006).

41. Krencker, F.-N., Bodin, S., Suan, G., Kabiri, L. & Immenhauser, A. Assessing the duration and possible causes of the earliest Toarcian carbon isotopic excursion. *EGU General Assembly Conference Abstracts.* 4970 (2013).

42. Suan, G. et al. Duration of the Early Toarcian carbon isotope excursion deduced from spectral analysis: consequence for its possible causes. *Earth Planet. Sci. Lett.* **267**, 666–679 (2008).

43. Suan, G., van de Schootbrugge, B., Adatte, T., Fiebig, J. & Oschmann, W. Calibrating the magnitude of the Toarcian carbon cycle perturbation. *Paleoceanography* **30**, 495–509 (2015).

44. Stordal, F., Svensen, H. H., Aarnes, A. & Roscher, M. Global temperature response to century-scale degassing from the Siberian Traps Large igneous province. *Palaeogeogr., Palaeoclimatol., Palaeoecol.* **471**, 96–107 (2017).

45. Pienkowski, G., Hesselbo, S., Barbacka, M. & Leng, M. J. Non-marine carbon-isotope stratigraphy of the Triassic-Jurassic transition in the Polish Basin and its relationship to organic carbon preservation, pCO2 and palaeotemperature. *Earth-Sci. Rev.* **210**, 103383 (2020).

46. Jones, M. T., Jerram, D. A., Svensen, H. H. & Grove, C. The effects of large igneous provinces on the global carbon and sulphur cycles. *Palaeogeogr., Palaeoclimatol., Palaeoecol.* **441**, 4–21 (2016).

47. Svensen, H. H., Polteau, S., Cawthorn, G. & Planke, S. Sub-volcanic intrusions in the Karoo Basin, South Africa in physical geology of shallow magmatic systems (eds Breitkreuz, C., & Rocchi, S.) (Advances in Volcanology, Springer, 2018).

48. Yemane, K., Siegenthaler, C. & Kelts, K. Lacustrine environment during low Beaufort (upper Permian) Karoo deposition in Northern Malawi. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **70**, 165–178 (1989).

49. Faure, K. & Cole, D. Geochemical evidence for lacustrine microbial blooms in the vast Permian Main Karoo, Paraná, Falkland Islands and Huab basins of southwestern Gondwana. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **152**, 189–213 (1999).

50. Scheffler, K., Buehmann, D. & Schwark, L. Analysis of late Palaeozoic glacial to postglacial sedimentary successions in South Africa by geochemical proxies – Response to climate evolution and sedimentary environment. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **240**, 184–203 (2006).

51. Branch, T., Ritter, O., Weckmann, U., Sachsenhofer, R. F. & Schilling, F. The Whitehill Formation – a high conductivity marker horizon in the Karoo Basin. *South Afr. J. Geol.* **110**, 465–476 (2007).

52. Svensen, H. et al. Nitrogen geochemistry as a tracer of fluid flow in a hydrothermal vent complex in the Karoo Basin, South Africa. *Geochimica et. Cosmochimica Acta* **72**, 4929–4947 (2008).

53. Götz, A. E., Ruckwied, K. & Wheeler, A. Marine flooding surfaces recorded in Permian black shales and coal deposits of the Main Karoo Basin (South Africa): Implicatipns for basin dynamics and cross-basin correlation. *Int. J. Coal Geol.* **190**, 178–190 (2018).

54. Jones, S. M., Hogget, M., Greene, S. E. & Dunkley Jones, T. Large Igneous Province thermogenic greenhouse gas flux could have initiated Paleocene-Eocene Thermal Maximum climate change. *Nat. Commun.* **10**, 5547 (2019).

55. Heimdal, T. H., Jones, M. T. & Svensen, H. Thermogenic carbon release from the Central Atlantic magmatic province caused major end-Triassic carbon cycle perturbations. *Proc. Natl Acad. Sci. USA* **117**, 11968–11974 (2020).

56. Goddéris, Y. & Joachimski, M. M. Global change in the Late Devonian: modelling the Frasnian-Famennian short-term carbon isotope excursions. *Palaeogeogr., Palaeoclimatol., Palaeoecol.* **202**, 309–329 (2004).

57. Jacob, R. *Low frequency variability in a simulated atmosphere ocean system*. Ph.D. Thesis, University of Wisconsin, Madison, USA (1997).

58. Besse, J. & Courtillot, V. Apparent and true polar wander and the geometry of the geomagnetic field over the last 200 Myr. *J. Geophys. Res.: Solid Earth.* **107**, EPM 6-1-EPM 6-31 (2002).

59. Freeman, K. H. & Hayes, J. M. Fractionation of carbon isotopes by phytoplankton and estimates of ancient CO2 levels. *Glob. Biogeochemical Cycles* **6**, 629–644 (1992).

60. Torsvik, T. H. & Cocks, L. R. M. *Earth history and palaeogeography*. (Cambridge University Press, Cambridge, United Kingdom, 2016).

61. Luttinen, A. V. Bilateral geochemical asymmetry in the Karoo large igneous province. *Sci. Rep.* **8**, 5223 (2018).

## Acknowledgements

This work was supported by the Research Council of Norway through its Centers of Excellence funding scheme (project number 223272), which included a field excursion to South Africa and Lesotho in 2018. M.T.J. was funded by the Research Council of Norway with project number 263000.

## Author contributions

T.H.H. and H.H.S. initiated the project that led to this manuscript. T.H.H. and Y.G. performed the GEOCLIM modelling. M.T.J. created Fig. 1. T.H.H. created the remaining figures and wrote the manuscript, but all authors contributed to discussions of the results and commented on the manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

Supplementary information The online version contains supplementary material available at https://doi.org/10.1038/s41467-021-26467-6.

Correspondence and requests for materials should be addressed to Thea H. Heimdal.

Peer review information *Nature Communications* thanks Guillaume Paris and Mariano Remírez for their contribution to the peer review of this work.

Reprints and permission information is available at http://www.nature.com/reprints

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_7.png" /></figure>

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.