Very large release of mostly volcanic carbon during the Palaeocene-Eocene Thermal Maximum

Abstract

The Palaeocene-Eocene Thermal Maximum (PETM) was a global warming event that occurred about 56 million years ago, and is commonly thought to have been driven primarily by the destabilization of carbon from surface sedimentary reservoirs such as methane hydrates. However, it remains controversial whether such reservoirs were indeed the source of the carbon that drove the warming. Resolving this issue is key to understanding the proximal cause of the warming, and to quantifying the roles of triggers versus feedbacks. Here we present boron isotope data-a proxy for seawater pH-that show that the ocean surface pH was persistently low during the PETM. We combine our pH data with a paired carbon isotope record in an Earth system model in order to reconstruct the unfolding carbon-cycle dynamics during the event. We find strong evidence for a much larger (more than 10,000 petagrams)-and, on average, isotopically heavier-carbon source than considered previously. This leads us to identify volcanism associated with the North Atlantic Igneous Province, rather than carbon from a surface reservoir, as the main driver of the PETM. This finding implies that climate-driven amplification of organic carbon feedbacks probably played only a minor part in driving the event. However, we find that enhanced burial of organic matter seems to have been important in eventually sequestering the released carbon and accelerating the recovery of the Earth system.

Extended Data Fig. 1. Elemental and stable isotope cross-plots for M. subbotinae measured in this study.

Extended Data Fig. 2. Foraminifera- and bulk carbonate stable isotope data plotted against depth in core.

Foraminifera-based stable isotope compositions were generated from identical samples after splitting of δ¹³C / δ¹⁸O fraction from the δ¹¹B fraction.

Extended Data Fig. 3. Illustration of δ¹¹B to pH conversion as well as age model differences.

(a) Comparison of pH evolution at Site 401 over the PETM CIE using either the borate ion (red) or alternatively the T. sacculifer (green) calibration. Age scale used is following Röhl et al.. (b) Direct comparison of our two age models, showing the reconstructed pH evolution of Site 401 plotted using either the age model of Farley and Eltgroth or our preferred age model of Röhl et al.. (c) Expanded view of (b).

Extended Data Fig. 4. Selection of age model tie points.

Bulk carbonate δ¹³C and δ¹⁸O comparison between Site 401 and Site 690 presented in Röhl et al.. Vertical lines highlight age tie points used to derive the age model relative to the PETM carbon isotope excursion (see methods for discussion).

Extended Data Fig. 5. Key results of sensitivity experiments.

Illustrating the influence of uncertainties in the CIE onset duration on diagnosed total carbon release. In these idealized experiments, the CIE onset phase is assumed to occur linearly, with a duration of the decline in δ¹³C (by 3.5‰) and pH (by 0.3 pH units) that varies from 100 to 20,000 yr, with the target pH and δ¹³C values thereafter held constant until the end of the experiment (50,000 yr). The evolution with time of these target ocean surface variables is shown in the uppermost panels (a), with pH on the left hand y-axis, and δ¹³C on the right hand y-axis. The lower rows of panels show: (b) maximum emission rate per time interval, (c) cumulative carbon emission for respective onset phase in EgC (1 Eg = 10¹⁸ g) and (d) average emitted δ¹³C per time interval.

Extended Data Fig. 6. Spatial and temporal evolution of mean annual surface ocean pH in cGENIE.

Illustrated both across the PETM and for comparison, modern pH patterns projected from preindustrial and into the future under RCP 6.0. Shown are: (a) Global and annual mean surface ocean pH (black solid line) across the PETM from experiment ‘R07sm_Corg’ (our central pH estimate, using the inorganic borate ion calibration and the RH07 age model, and including an assumption of organic carbon burial post peak PETM). Red circles represent the annual mean pH values at the location of Site 401 in the model (see location in panel b) taken at times in the model simulation that have a corresponding δ¹¹B derived pH data points (cf. Fig. 3b) (but note that we do not utilize all of the observed data points). (b) Model projected spatial pattern of annual mean surface ocean pH at time zero (i.e. PETM onset). (c-f) Model projected spatial pattern of the annual mean surface ocean pH anomaly compared to time zero, for the highlighted time-points in (a) – 5.0, 31.6, 58.2, and 71.5 kyr following onset. (g) Model projected spatial pattern of annual mean surface ocean pH in the modern ocean under pre-industrial atmospheric CO₂ (278 ppm). The model is configured as per described in Cao et al. and driven with a CO₂ emissions scenario calculated consistent with RCP 6.0. (h-i) Model projected spatial pattern of the annual mean surface ocean pH anomaly compared to 1765, at year 2010 and 2050. The scale is chosen to be the same as per (c-f).

Extended Data Fig. 7. Spatial and temporal evolution of surface sedimentary carbonate content in cGENIE across the PETM.

(a) Global mean surface sedimentary wt% CaCO₃ (black solid line) across the PETM from experiment ‘R07sm_Corg’. White circles represent the times from PETM onset onwards that correspond to the δ¹¹B derived pH data points as per in Fig. 3b and Extended Data Fig. 6. Note that the white circles do not represent ‘values’ and are plotted simply as markers of specific time-points (see Extended Data Fig. 6). (b) Model projected spatial pattern of surface sedimentary wt% CaCO₃ at time zero (i.e. PETM onset). Shown are the locations of sites for which surface ocean pH has been reconstructed (see Fig. 2) and at which detailed down-core model-data comparison is carried out (Extended Data Fig. 9). (c-f) Model projected spatial pattern of the surface sedimentary wt% CaCO₃ anomaly compared to time zero, for the highlighted time-points in (a) – 5.0, 31.6, 58.2, and 71.5 kyr following onset. (g) For reference – the assumed seafloor bathymetry in the model (together with the locations of the four data-rich sites focussed on in the SI analysis).

Extended Data Fig. 8. Spatial and temporal evolution of sea surface temperature in cGENIE across the PETM.

(a) Global and annual mean sea surface temperature (SST) (black solid line) across the PETM from experiment ‘R07sm_Corg’. Yellow circles represent the annual mean SST values at the location of Site 401 in the model at the times from PETM onset onwards that correspond to the δ¹¹B derived pH data points (cf. Fig. 3b). Orange and blue filled circles represent Mg/Ca and δ¹⁸O derived, respectively, SST estimates. (b) Model projected spatial pattern of annual mean SST at time zero. The location of Site 401 in the model is highlighted by a star. (c-f) Model projected spatial pattern of the annual mean SST anomaly compared to time zero, for the highlighted time-points in (a) (yellow circles) – 5.0, 31.6, 58.2, and 71.5 kyr following onset.

Extended Data Fig. 9. Down-core model-data evaluation at four data-rich sites.

Shown are comparisons for four ocean drilling sites for which surface ocean pH has been reconstructed across the PETM (Fig. 2) – 401, 865, 1209, and 1263 (this study and ref. 20). Their paleo locations in the cGENIE Earth system model are shown to the side (panel q). Model-data comparisons are made for: (i) wt% CaCO₃ (far LH panel for each site), (ii) δ¹³C of bulk carbonate (second-from-left series of panels), and (iii) surface ocean pH (third-from-left series of panels). To provide an orientation in time with regard to the evolution across the PETM event, the farthest-right series of panels shows the projected evolution of atmospheric δ¹³C of CO₂ in the model. For wt% CaCO₃ and δ¹³C of bulk carbonate, model points (resolved at 1 cm resolution) are plotted as filled yellow circles. Model-projected pH (global and annual mean, as per shown in Fig. 3j and Extended Data Fig. 6a) and atmospheric δ¹³C of CO₂ are shown as continuous red lines. In all cases, observed data values are shown as stars (*). The age models for Sites 865, 1209 and 1263 employing original relative age model constraints used to convert from model-simulated sediment depth (resolved at 1 cm intervals) at each location in the cGENIE Earth system model, are calculated using a constant detrital flux accumulation rate. The observed data are plotting on their respective site 690-derived age models. Both model and data age scales are synchronized to age zero at PETM onset (horizontal line). See SI for details.

Fig. 1. New DSDP Site 401 stable isotope data.

Foraminifera (M. subbotinae) (a) and bulk carbonate δ¹³C (b), δ¹¹B (c) and δ¹⁸O (d and e) records plotted relative to the onset of the PETM carbon isotope excursion (CIE) from DSDP Site 401 (47° 25.65’ N, 08° 48.62’ W, 2495 m) using our preferred age model (see Methods).

Fig. 2. M. subbotinae based δ¹³C and boron isotope based pH reconstructions of Site 401.

Panels A and B show the entire record, while C and D focus on the CIE interval. Also shown are data of ref. on the original age model with pH values recalculated using a laboratory offset such that pre-PETM pH calculated using our Monte Carlo approach at Site 1209 = 7.74 given the distribution of seawater δ¹¹B determined at Site 401 (38.9 ± 0.4‰). This resulted in a mean correction of the literature data of -0.32‰.

Fig. 3. Results of Earth system model data assimilation.

The right hand panels also account for organic carbon burial during PETM recovery. (a,i) Atmospheric pCO₂ (red, LH axis) and mean global SST (blue, RH axis). (b,j) Modelled mean global ocean surface pH (observed smoothed surface ocean pH data as yellow symbols). (c,k) Model diagnosed rates of CO₂ release (red) and excess CO₂ consumption due to silicate weathering (green) from PETM onset onwards. (d,l) Cumulative CO₂ release (red) and organic carbon burial (blue). (e,m) Modelled mean global ocean surface δ¹³C (observations as yellow symbols). (f,n) Model diagnosed δ¹³C of the CO₂ release (red) and isotopic composition of buried carbon (blue). Shaded bands (a,b,e,i,j,m) and empty bars (c,d,f,k,l,n) reflect 95% uncertainty limits. Bars reflect 2 kyr averaging (c,f,k,n) or integration (d,l) bins. All model results and related data are plotted from -50 to +150 kyr relative to the onset of the CIE, on our preferred orbital age model.