{'fraction_non_alphanumeric': 0.034093959731543624, 'fraction_numerical': 0.028944478340451496, 'mean_word_length': 4.675346260387812, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Very large release of mostly volcanic carbon during the Paleocene-Eocene Thermal Maximum

Global warming during the Palaeocene-Eocene Thermal Maximum1,2 (PETM, ~56 Ma) is commonly interpreted as being primarily driven by the destabilization of carbon from surficial sedimentary reservoirs such as methane hydrates3. However, the source(s) of carbon remain controversial1,3–5. Resolving this is key to understanding the proximal cause, as well as quantifying the roles of triggers versus feedbacks in driving the event. Here we present new boron isotope data – a proxy for seawater pH – that demonstrate the occurrence of persistently suppressed surface ocean pH across the PETM. Our pH data, alongside a paired carbon isotope record, are assimilated in an Earth system model to reconstruct the unfolding carbon cycle dynamics across the event6,7. We find strong evidence for a much larger (>10,000 PgC) and on average isotopically heavier carbon source than considered previously8,9. This leads us to identify volcanism associated with the North Atlantic Igneous Province, rather than carbon from a surficial reservoir, as the main driver of the PETM10,11. We also find that, although amplifying organic carbon feedbacks with climate likely played only a subordinate role in driving the event, enhanced organic matter burial was important in ultimately sequestering the released carbon and accelerating the recovery of the Earth system12.

suppressed surface ocean pH across the PETM. Our pH data, alongside a paired carbon isotope record, are assimilated in an Earth system model to reconstruct the unfolding carbon cycle dynamics across the event6,7. We find strong evidence for a much larger (>10,000 PgC) and on average isotopically heavier carbon source than considered previously8, 9. This leads us to identify volcanism associated with the North Atlantic Igneous Province, rather than carbon from a surficial reservoir, as the main driver of the PETM10, 11. We also find that, although amplifying organic carbon feedbacks with climate likely played only a subordinate role in driving the event, enhanced organic matter burial was important in ultimately sequestering the released carbon and accelerating the recovery of the Earth system12.
Aside from climate13 and ecological sensitivities14, arguably the greatest uncertainties surrounding the response of the Earth system to massive carbon release concern the role of carbon-cycle feedbacks15. A past event with considerable potential to evaluate such feedbacks is the Palaeocene-Eocene Thermal Maximum (PETM)1 -a 4-5°C transient surface warming2 associated with ecological disruption occurring around 55.8 million years ago16. Estimates of total carbon release vary from ~3,000 PgC to over 10,000 PgC7,8, spanning the range of present-day fossil fuel reserves17 but equally reflecting considerable uncertainty in current understanding. The source(s) of carbon is also highly uncertain, and has been proposed to involve methane hydrates3, permafrost4 and marine sedimentary5 organic matter. To further complicate the matter, proposed triggers for the PETM include orbital variations4 and an extraterrestrial impact18. Massive flood basalts and sill emplacement occurring around the time of the PETM and associated with the North Atlantic Igneous Province (NAIP) 10,11,19, constitute an additional potential source of carbon, but one not linked to a feedback with climate. If we are to fully understand the paleo-record, as well as exploit it to improve our understanding of the longer-term consequences of anthropogenic carbon emissions, we must resolve the balance of carbon source(s) that gave rise to the PETM, and thereby deconvolve the role(s) of triggers versus feedbacks. To provide new insight into the amount and source of carbon involved in PETM warming, we present new, paired, surface ocean boron (a well-established proxy for ambient surface seawater pH20,21) and carbon isotope data, and simultaneously use these to constrain the time-varying sources and sinks of carbon across the PETM in a novel data assimilation approach in an Earth System model (ESM).
We generated near-continuous boron, oxygen and carbon isotope records from NE Atlantic DSDP Site 401, using the surface ocean mixed-layer dwelling foraminifer Morozovella subbotinae (Fig. 1). We sampled the sediment sequence over an interval corresponding to lowest δ 13 C and δ 11 B values are both observed about ~25 ka after the onset of the CIE in our preferred age model, giving an inferred duration of the onset phase of the CIE in good agreement with an independently dated record from Spitsbergen16.
Because of uncertainties in early Cenozoic seawater boron isotopic composition (δ 11 B SW ), we tie our initial, pre-CIE boron isotope derived pH to mean ocean pH (7.75) as simulated by the 'GENIE' Earth System Model (ESM)6 and following the approach of a previous PETM model-data pH study20. Our δ 11 B measurements then dictate the timing and magnitude of how ocean pH deviated from this value across the PETM. In our pH reconstruction, we calculate an uncertainty envelope accounting for uncertainties in surface ocean temperature and salinity plus δ 11 B measurement errors, and test two contrasting endmember δ 11 B-pH calibrations for the extinct foraminifer M. subbotinae (see Methods). We focus on the δ 11 B foram = δ 11 B borate calibration, giving an estimated δ 11 B SW (38.9 ± 0.4‰) consistent with a recent reconstruction of Eocene δ 11 B SW based on δ 11 B21.
Evolution of ocean pH across the PETM is characterized by a negative excursion of 0.27 (range: 0.18-0.41) or 0.36 (0.21-0.56) pH units, depending on which δ 11 B-pH calibration is used ( Fig. 2 and Extended Data Fig. 3a), and in general agreement with a recently published PETM δ 11 B record20 (Fig. 2). The wide geographic distribution, but close correspondence in magnitude of all PETM δ 11 B-pH records (Pacific, S. Atlantic and N. Atlantic) gives us confidence that a global surface pH excursion signal is captured at DSDP Site 401. The fact that ocean surface pH responds relatively uniformly in models14 supports the evidence from multiple δ 11 B records (Fig. 2) that a single open ocean site can be representative of the global trend (see Methods).
To reconstruct PETM carbon release and its average isotopic composition, we devised a novel data assimilation methodology. We build on previous work7 in which a single δ 13 C record was assimilated ('inverted') to constrain the time-varying addition of carbon, but here exploit a more direct indicator of carbon addition -ocean surface pH (Fig. 2). This allows our δ 13 C record to simultaneously provide a second, independent constraint on the isotopic composition of the carbon emissions in a transient, 500 kyr duration assimilation of both records (see Methods). We explored a wide range of different model parameterizations and proxy assumptions (Extended Data Table 1a) but focus here on the results of the data assimilation of the smoothed record.
With our preferred age model ('R07sm', Extended Data Table 1a) we diagnose a cumulative PETM carbon release reaching ~10,200 PgC with almost all emissions occurring in the first 50 kyr (Fig. 3d). This estimate is largely independent of the choice of age model (Extended Data Table 1), which primarily affects the cumulative carbon emissions associated with the onset interval itself (defined as: from the first trace of the δ 13 C decline in our records up to peak CIE values) rather than with total emissions associated with the event as a whole. We demonstrate this in idealized model experiments (Extended Data Fig. 5 and Extended Data Table 1b) in which we find total carbon emissions over 50 kyr essentially independent of the assumed duration of the onset interval, and varying by only ±20% at the 20 kyr horizon (Extended Data Fig. 5 and Extended Data Table 1b). Thus, it is the extended duration of low pH across the PETM as a whole and the existence of the so-called carbon isotope 'plateau'2, rather than the duration of the onset interval alone, that lead to our diagnosis of total PETM emissions on the order of 10,000 PgC. Uncertainty in the duration of low pH equates to ~100 PgC kyr -1 at the 50 kyr horizon (Extended Data Fig. 5), consistent with the ~12,000 PgC total emissions deduced for our alternative age model with an extended duration of low pH (Extended Data Fig. 3c).
In response to carbon emissions, atmospheric pCO 2 in the model increases from ~866 to a peak PETM value of 2176 +1904/-669 μatm, consistent with independent atmospheric pCO 2 constraints based on variable terrestrial and marine δ 13 C gradients over the PETM22. The corresponding projected annual mean sea surface temperature (SST) increase is 3.6°Cclose to the observation-based global mean warming estimate of 4-5°C2. Also in response to carbon emissions (and surface ocean pH suppression), there is a shoaling of the carbonate compensation depth (CCD) in the model -the depth horizon below which calcium carbonate (CaCO 3 ) is not preserved23 (Extended Data Fig. 7). In previous global carbon cycle model analyses of the PETM, the CCD has been used as a data constraint, with the conclusion that carbon emissions on the order of 10,000 PgC are too high8. In contrast, here, the relatively long (>50 kyr) duration of low ocean pH conditions (Fig. 3) in conjunction with weathering feedbacks, leads to a partial decoupling of pH and ocean carbonate saturation24, hence a relatively muted response of the CCD despite the large emissions (Extended Data Fig. 7 and Methods).
Diagnosed carbon emission rates peak at 0.58 PgC yr -1 ( Fig. 3c; Extended Data Table 1a), although we assign rather less confidence to these, because their value is sensitive to the duration of the onset of the PETM and hence the specific age model (Extended Data Table  1a). To put this in perspective, for carbon input rates to approach those of current fossil fuel emissions (~10 PgC yr -1 17), the PETM onset would have to occur within 200-500 yr -a duration not supported by any independent age model7,16,24,25. The much lower than modern carbon emissions rate we diagnose here then implies reduced PETM ocean acidification impacts (especially in carbonate saturation) compared to the future6,15. However, we cannot rule out multiple, short-lived pulses of carbon release >0.58 PgC yr -1 having occurred throughout an extended (e.g. 20 kyr) onset24.
In addition to the emissions diagnosed by matching the pH decline, using the δ 13 C data as an independent constraint leads us to deduce a flux-weighted mean δ 13 C of released carbon of -11‰ (Fig. 3f, n). However, the smoothed δ 13 C record (-2.6‰ excursion) on which we focus, very likely underestimates the isotopic magnitude of the event. For instance, if the 'true' PETM CIE was as large as -4.0‰7,24 and we simply proportionally scale δ 13 C input , diagnosed on the basis of a -2.6‰ excursion, we obtain a more depleted mean source of -17‰. Uncertainty in our ocean pH reconstruction also affects the diagnosed carbon source composition. Our minimum pH decrease of 0.18 pH units requires only 5,700 PgC, with a mean δ 13 C input of -19‰. However, the comparatively muted surface warming seen in this 'minimal pH change' model experiment (2.25°C, Extended Data Table 1a -experiment 'R07am_HI') is difficult to reconcile with an observed warming of 4-5°C2. Conversely, the upper end of our measured pH increase would require emission of considerably more carbon (19,960 PgC) with a correspondingly heavier carbon isotopic composition of -6.6‰ (Extended Data Table 1a).
Our diagnosed carbon input over the event likely reflects a combination of carbon source(s) -for instance, a mean of -11‰ could reflect a 75% contribution of mantle-derived carbon (δ 13 C source ~-6‰) plus 25% from permafrost (~-26‰), or 90% mantle-derived plus 10% methane hydrates (~-60‰). In such scenarios, volcanism triggered the PETM, and thawing permafrost in Antarctica4 or destabilization of methane hydrates provided amplifying feedback. Assuming a -4‰ magnitude excursion and mean δ 13 C input of -17‰ still requires a substantial CO 2 contribution from volcanism10, but would allow for the possibility of a greater role for organic carbon feedbacks -almost 60% for organic matter or ~20% for methane hydrates.
To date, the PETM has predominantly been viewed as an event dominated by feedbacks between climate and reservoirs of carbon3. Yet, there is abundant evidence of an intimate link in time with the opening of the North Atlantic11, with volcanism and ash deposition occurring from immediately prior to PETM onset, as also recorded by declining 187 Os/ 188 Os in sediments19. Radiometric dating places the PETM coincident with a ~1 Myr interval of massive flood basalt volcanism11 and the emplacement of magmatic sills26, both of which represent large carbon sources. Degassing CO 2 from magma yields an estimated 3,600-6,000 gC m -3 27 and combining this with the estimated volume of the NAIP as a whole (5×10 6 km 3 to 10×10 6 km 3 11,27), equates to a potential carbon source of 18,000-60,000 PgC. The interaction of magmatism with organic rich sediments could enhance carbon release via thermogenic methane production10,11, which is estimated to range from 3,000-6,000 PgC28 to as high as 15,000 PgC10. Available carbon reservoirs are thus more than sufficient to provide the 10,200-12,200 PgC required by our data assimilation and we further note that an all-volcanic carbon driver scenario for the PETM is possible if thermogenic methane10,11 provided the isotopically lighter end-member. On the other hand, NAIP magmatic activity took place over several million years10,11 and how carbon emissions were distributed with time over this interval is currently unknown. Dating, biostratigraphy, and seismic constraints do however: (1) place an interval of volcanism in East Greenland11 and sill emplacement in the Vøring Basin (offshore Norway)26, both coeval with PETM onset, (2) identify 100s of degassing structures consistent with thermogenic carbon release as forming close to the P-E boundary10 and with one structure constrained to have been active during the body of the PETM itself9. Release of a disproportionate amount of NAIP carbon associated with the PETM is hence consistent with available geological evidence as well as our data-inferred carbon source and total release. More work dating further specific volcanic episodes and refining carbon reservoir estimates is however clearly needed.
Our paired δ 11 B-δ 13 C data also provide insights into climate system recovery from PETM warming. Once carbon emissions ceased (ca. ~55 kyr after PETM initiation - Fig. 3c), elevated global temperatures (Fig. 3a) and enhanced rates of silicate weathering (Fig. 3c) in cGENIE (see Methods) drive a trend of increasing ocean surface pH that closely follows the observed surface ocean pH recovery (Fig. 3b). However, we find a model-data misfit of up to ~1‰ in δ 13 C during the recovery phase (Fig. 3e). We therefore performed an additional set of experiments in which, after peak CIE, organic carbon (C org ) is removed from the ocean surface29 and assumed buried whenever modelled mean ocean surface δ 13 C registered lower values than the observed trend (see Methods). These final experiments provide close agreement with the recovery trend in the δ 13 C data (Fig. 3m), with cumulative C org burial ( Fig. 3l, blue bars) of 2,500 PgC (at an average modelled marine value of -30.5‰), in agreement with other estimates (~2,000 PgC)12 of the role of enhanced organic matter burial in PETM recovery29 as well as the ensuing reduction in deep-sea oxygenation30.
These findings collectively lead us to a view of the PETM as having been on the smaller end of a spectrum of severe perturbations of climate and carbon cycling during the Cretaceous and Jurassic (Ocean Anoxic Events -OAEs30,31), despite it having been by far the largest end-member in a series of Paleocene-Eocene 'hyperthermal' events32. Our pH reconstruction, in conjunction with the observed δ 13 C decline, constrains the dominant carbon source during the PETM onset to have had a comparatively heavy carbon isotope ratio, strongly implicating volcanism as having been dominant in triggering and driving the event. Our inferred mean δ 13 C source of -11 to -17‰ is consistent with the isotopically relatively heavy source (ca. -15‰33) inferred for the end-Permian event, suggesting mechanistic similarities between the two events27. The implied important role for organic carbon deposition in the recovery from peak warming12 represents another diagnostic feature of OAEs31 (and end-Permian). Further quantifying and understanding the precise role of feedbacks -both those amplifying initial CO 2 release, and those aiding recovery from global warming -is arguably where the PETM is of greatest value in helping reduce uncertainties surrounding the response of the global carbon cycle and climate system to perturbation.

Site and sample selection
The open northeast Atlantic DSDP Site 401 (47° 25.65' N, 08° 48.62' W, 2495 m) was selected for this study. Its depth during the PETM was approximately 2000 m34. Around 2 mg of the 250-300 μm size fraction of mixed-layer dweller Morozovella subbotinae were picked for the carbon, oxygen and boron isotopic analyses. Furthermore, over the studied interval, very high-resolution δ 18 O and δ 13 C analyses of bulk carbonate were conducted to establish a revised age model for Site 401. Planktic foraminifera are extremely well preserved at Site 40135, free from infilling and, particularly from the onset of the CIE upwards, are semi-glassy in appearance36.

Sample treatment
Using a binocular microscope, picked foraminifera were cracked open under glass plates, the sample then homogenised, before splitting into a fraction for stable isotope (δ 18 O and δ 13 C) analysis and another for the boron isotopic and elemental analyses (with a ratio of ca. 10:90). Purification and measurement of the boron fraction followed established protocols37,38. Samples were thoroughly cleaned to remove any adhering clays and samples were oxidatively cleaned using buffered peroxide in a warm water bath closely following39. Boron isotopic and elemental analyses were carried out on a Thermo Scientific Neptune MC-ICPMS and Element XR ICPMS, respectively, at the University of Southampton. Sample purification and handling was done in low-boron clean labs at the University of Southampton. The average boron total procedural blank was on the order of 30 to 50 pg (n>10) and is hence negligible given our typical sample size (~5 to 15 ng of B). Boron isotopic uncertainties are reported at the 2 sigma level calculated using repeats of in-house carbonate standards40. Boron isotopic and elemental aliquots were measured using additional ammonia gas for better sample washout between samples and strictly monitored during every analytical session37. Prior to analysis for boron isotopic composition, samples were screened for chemical consistency by checking various elemental ratios (B/Ca, Mg/Ca, Al/Ca etc.) (Extended Data Fig. 1). While few samples had elevated Al/Ca (up to ~ 3400 μmol/mol) this feature did not translate into altered δ 11 B (Extended Data Fig. 1).
Carbon and oxygen isotope aliquots were measured on a Thermo Finnigan MAT252 stable isotope mass spectrometer at the GEOMAR Helmholtz Centre for Ocean Research Kiel, Germany. Additionally, some foraminifera-based δ 18 O and δ 13 C analyses as well as all bulk carbonate stable isotope measurements were carried out at the MARUM Bremen, Germany on a Finnigan 251 gas isotope ratio mass spectrometer, coupled to a Kiel I automated carbonate preparation device. All produced isotope records are shown in Extended Data Fig.  2 plotted against depth in core. The carbon isotope excursion seen in our record is 3.4‰, significantly expanded relative to the benthic carbon isotope excursion presented by Nunes and Norris41 that only reported an excursion on the order of 1.8‰. This discrepancy arises from the lower resolution data this earlier study41 and the fact that samples were not taken through the core interval of the CIE at Site 401 (202.55 to 202.41 mcd) in this earlier study. We note that Bornemann et al. 35 reproduced a very similar magnitude of change in δ 13 C to us; their δ 13 C data obtained from the same species (Morozovella subbotina) registered a shift from 4.87‰ at 202.58 mcd to 1.47‰ at 202.46 mcd (an identical excursion magnitude of 3.4‰). The core containing the PETM (core 401-14) shows some rotary-drilling induced core deformation across the CIE. Such deformation commonly occurs across abrupt changes in lithology, but there is no obvious coring gap35.

Effect of δ 11 B-pH calibration used on resulting pH excursion
Using the appropriate δ 11 B-pH calibration in order to convert calcite δ 11 B into ambient seawater pH is essential for any paleo-pH reconstruction. For late Neogene studies using extant foraminifer species, the species used are typically calibrated for their δ 11 B calcite to pH dependency using culture or field studies42,43 in order to assess the magnitude of δ 11 B-vital effects that relate to foraminiferal physiology44-46. However, the species used here is extinct, making such calibrations impossible.
In order to bracket the likely magnitude of vital effects, and following ref. 21, we present two calibrations, one using the δ 11 B to pH relationship of aqueous borate47 and the other using the T. sacculifer calibration43. While the aqueous borate calibration is used for pH trends shown in Figs. 2 and 3, Extended Data Fig. 3a also present the alternative outcome. As noted previously20,46, when pre-PETM pH is fixed (as is the case here), the choice of δ 11 B-pH calibration has little impact on the reconstructed pH curve. We note that the aqueous borate ion calibration is more conservative and is our preferred option. This is for the following reasons: (i) not all modern species show a reduced sensitivity to pH relative to aqueous borate48; (ii) previous studies have argued for a reduced magnitude of δ 11 B vital effects in Eocene foraminifera.

Determination of δ 11 B sw
Boron in seawater has a residence time of between ~11 to 20 Ma52,53 and to date the δ 11 B sw is not well constrained for the PETM. In order to create a self-consistent model-data setup we therefore used the output of GENIE ESM in the pre-CIE configuration which for the open NE Atlantic provides a pH of 7.756. Using this pH information and employing the generic borate ion calibration47 for the pH-dependent incorporation of boron into the studied foraminifera Morozovella subbotinae resulted in a δ 11 B sw of 38.94 ± 0.41‰. The uncertainty in deriving this bulk seawater δ 11 B is based on 10,000 realizations of a borate ion to pH conversion using the commonly used experimentally derived boron fractionation factor47, varying the given δ 11 B randomly within its 2 sigma measurement uncertainty, and also varying salinity by ±1.5 psu and temperature by ±1.5°C. Utilising the T. sacculifer δ 11 B-pH calibration43, but following the same approach, gives a δ 11 B sw = 37.6 ± 0.5‰.

Chronology for Site 401
A new and detailed age model was established for Site 401 by aligning our new ultra-high resolution (1 cm-spacing) bulk carbonate δ 18 O and δ 13 C records with equivalent bulk carbonate isotope records from Site 690 using the 'Analyseries' software54. Most stratigraphic correlation tie points (vertical lines in Extended Data Fig. 4) were made using the δ 18 O records, which gave excellent agreement between the sites. The bulk δ 18 O record from Site 401 shows high structural similarity to the δ 18 O of the mixed layer-dwelling planktic foraminifer M. subbotinae from this same site (Extended Data Fig. 2), and also to the δ 18 O of thermocline-dwelling S. patagonica35, suggesting that bulk sediment δ 18 O at Site 401 provides a reliable record of the basic trends in upper ocean warming and cooling across the PETM. A dominant control by temperature on the bulk δ 18 O signal makes sense, given the scale of global surface ocean warming across the PETM2 (4-5°C). The fidelity of the bulk δ 13 C record from Site 401 is supported by the fact that it shows high structural similarity to the δ 13 C of mixed layer-dwelling M. subbotinae (Extended Data Fig. 2), and also to the δ 13 C of thermocline-dwelling S. patagonica35. It is also consistent with bulk δ 13 C from another nearby location (the Forada section in northern Italy) that also shows an unusually early recovery to higher δ 13 C following the initial excursion to lowest δ 13 C at the PETM's onset55. The Forada section is considered to be complete, because the CIE interval covers the maximum number of precession cycles25. Site 690 currently has two detailed age models. By detailed correlations to Site 401, we were thus able to transpose both the astronomically calibrated chronology25,56 and an extra-terrestrial He-based chronology57 onto Site 401. Extended Data Figs. 3b and c compares our pH record from Site 401 on both chronologies. These uncertainties relating to choice of age model, and their impact on the calculated duration of the onset phase, have been evaluated via modelling sensitivity experiments (Extended Data Fig. 5) and have no impact on our main findings as discussed in the main text.
A very different timescale for PETM carbon release during the CIE was suggested in an earlier study, arguing for an onset of the PETM CIE within only 13 years58. The proposal of a CIE onset within such a short timescale has proven controversial59-63. In particular, geochemical modelling constraints59 as well as drilling disturbance of the core creating the impression of annual layering have together cast significant doubt on the suggested very rapid (~13 year) CIE onset. Indeed, Further Earth system model based analysis of the carbon and oxygen isotope records, leads to an estimate of 4 kyr or longer for PETM onset64. Given previously presented age constraints for the duration of the PETM CIE based on cyclostratigraphy25 and a 3 He-based age model from ODP Site 69057 in addition to absolute and cyclostratigraphic age constraints from Spitsbergen16, we regard an age model that leads to a multi-millennia-scale CIE onset as more plausible. However, as analysed (Extended Data Fig. 5) and discussed earlier, assumptions regarding the duration of PETM onset interval itself are not critical to our conclusions. We introduce three separate model innovations here. The first builds on previous work7,71 'inverting' an observed δ 13 C record to recover the underlying time-history of carbon release. In this, cGENIE adjusts mean atmospheric or surface ocean δ 13 C to match a (proxy data) target at each model time-step (~1 week). If the current mean model value lies above the data value (observed data is automatically linearly interpolated to the model time-step), a pulse of carbon is released to the atmosphere (or ocean). If the model lies below the data value, depending on the experimental setup, carbon is either removed from the atmosphere, or nothing is done (cf. Fig. 3). The magnitude of the carbon pulse emitted at each time-step is prescribed and chosen such that the fastest observed change in δ 13 C can be closely tracked, but without creating excessive overshoots in modelled δ 13 C. Here, we allow a maximum rate of carbon emissions to the atmosphere of 10 PgC yr -1 and hence a magnitude of an individual pulse of ~0.21 PgC, corresponding to an instantaneous increase in atmospheric pCO 2 of about 0.1 ppm.

Earth system modelling -configuration and data inversion methodology
We diverge from an earlier approach7,71 in that rather than utilizing a record of δ 13 C as our model target to assimilate, we instead employ our Site 401 reconstructed surface ocean pH record. The methodology is inherently the same, but rather than comparing mean model and observed δ 13 C each time-step, we contrast (model and data) pH, diagnosing the required carbon flux to the atmosphere in order that surface pH in the model tracks the data. The model-data comparison is done on the basis of a mean global surface ocean pH value calculated in cGENIE because utilizing a single (Site 401) surface ocean grid point in cGENIE creates artefacts in the diagnosed carbon emissions because there is seasonality in pH in the model but not in the data. We justify the assumption that proxy reconstructed surface ocean pH at Site 401 can be representative of the global mean, firstly on the basis of the relatively close degree of correspondence (visually) between the globally distributed pH records available, as show in Fig. 2. Secondly, ocean surface pH, both today and during the Paleocene-Eocene, is relatively uniform in the model (and supported by observations and proxies, respectively), with maximum surface gradients between upwelling regions and subpolar regions of no more than 0.1 pH units for modern, and considerably less than this in the late Paleogene (likely primarily due to the non-linear nature of the pH scale) (Extended Data Fig. 6). Furthermore, these muted patterns are retained largely unaltered in response to CO 2 emissions. For instance, when we calculate the annual mean surface ocean pH anomaly at different times across the PETM (experiment ID 'R07sm_Corg') as compared to the pre-PETM pattern, we find a generally uniform (to within ±0.02 pH units) pattern in pH change (Extended Data Fig. 6). If we contrast the evolution of global and annual mean surface ocean pH across the PETM ('R07sm_Corg') with the annual mean surface pH at the location of Site 401 for the time points available (Extended Data Fig. 6, top), we also find Site 401 pH is globally representative (and vice versa). All this goes to illustrate that there is unlikely to be any substantive artefact in our assumption of treating our pH record at Site 401 as a surrogate for the global mean in the model inversion experiment. Finally, and for comparison, a similar analysis for the modern ocean under a future ocean acidification scenario (here, chosen to follow RCP6.072) is shown in Extended Data Fig. 6 and demonstrates a comparably spatially uniform pattern of pH change.
The second innovation involves the determination of the δ 13 C of the carbon emitted to the atmosphere. Previously7,71, the δ 13 C of the carbon was treated as an unknown and a range of different possible values (and hence carbon sources and reservoirs) tested in turn. However, since observed pH constrains the magnitude of carbon emissions, we can now simultaneously employ our observed δ 13 C record to determine the source of carbon. The way in which the 'double inversion' methodology then works is that on each model timestep, following the assessment of whether or not a pulse of carbon is emitted to the atmosphere (based on the model-data pH difference), mean global model and observed Site 401 δ 13 C values are compared. If the current mean model surface ocean δ 13 C value lies above the current data value, the carbon emitted is assigned a carbon isotopic value of -100‰. If however, the mean model value lies below the data value, an isotopic value of 0‰ is assigned to the carbon values. By binning the emission fluxes in time and calculating a flux-weighted average δ 13 C, as per in Fig. 3, intermediate (between -100 and 0‰) δ 13 C values are achieved. We emphasize that we are not assuming a source that could be -100‰ per se -this choice of extremely depleted value simply gives the model greater flexibility in tracking the trend in δ 13 C emissions -isotopically intermediate mean annual carbon emissions arise by varying proportions of individual 0‰ and -100‰ carbon pulses. We could have used any value just as long as it is as least as light as the lightest potential source (e.g. -60‰). Gutjahr

Europe PMC Funders Author Manuscripts
Finally, in the situation that the mean model surface ocean δ 13 C value becomes lower than the observed Site 401 value, we also test the importance of marine organic carbon (C org ) burial. This works identically to the negative emissions diagnosed in previous studies7,71 (when carbon is removed from the system to force δ 13 C more positive) but rather than prescribing the δ 13 C value, we calculate it according to a simple phytoplankton organic matter fractionation scheme67,73.
For all our experiments, we first spun up the model under late Paleocene boundary conditions24,64, here choosing an open system run time of 200 kyr in order to fully bring the long-term δ 13 C cycle into balance (and following on from an initially closed system spin-up of 20 kyr used to established the basic climate and ocean circulation state). We then carried out a range of experiments as summarized in Extended Data Table 1a. We tested combinations (not all are reported here) of: (i) age model -orbital cyclostratigraphy ('R07') vs. 3 He-based age model ('FE'), uncertainty in the pH reconstruction -mean vs. the 2.5% and 97.5% confidence limits ('LO' and 'HI', respectively), whether or not the data is smoothed ('sm') or raw ('rw'), whether or not climate-dependent weathering feedback was allowed, or weathering was fixed ('noW'), and whether or not C org burial was enabled to recover δ 13 C to more positive (and data tracking) values (Corg when carbon burial was enabled). These experiments were run for 500 kyr, with the exception of the carbon burial Corg series of experiments (Extended Data Table 1a), which were run for an initial interval of 72.6 kyr and up until the peak of the CIE with no organic carbon burial allowed, and then a further 227.4 kyr with carbon burial allowed when needed (for a total of 300 kyr of simulation). Model results are plotted relative to the observed data point defining PETM onset.

Earth system modelling -additional sensitivity experiments and analysis
We also carried out a range of sensitivity experiments to explore the importance (or otherwise) of the assumed duration of the CIE onset -in other words, whether there is a strong age model dependence of diagnosed total carbon emissions. In this series of experiments, the CIE onset phase was assumed to occur as a simultaneous linear decline in both δ 13 C (by -3.5‰) and pH (by -0.3 pH units). We varied the duration of this decline from 100 to 20,000 yr. Once the minimum in δ 13 C and pH was reached, these values were held constant up until the end of the experiment (a total of 50 kyr). The exact same double inversion methodology was employed and starting from the same spin-up state as the main experiments. The results of these sensitivity experiments are plotted in Extended Data Fig. 5 and summarised in Extended Data Table 1b.
Further details of the model and its paleo configuration, plus comprehensive discussion of model uncertainties, can be found in the supplementary information file SI 1. Additional assessments of the evolution of model-projected global mean as well as spatial patterns of sedimentary wt% CaCO 3 and sea-surface temperature are illustrated in Extended Data Figs. 7 and 8, respectively (and described in SI). Site-specific model-data comparisons are shown in Extended Data Fig. 9 (and again discussed in full in SI 1).

Europe PMC Funders Author Manuscripts
Earth system modelling -model code and supporting file availability The source code of the cGENIE Earth system model used to generate the results presented in this paper, together with specific experiment configuration, boundary conditions, and dataforcing files, is available for download. A brief overview of and directions to: obtaining the code and configuring the cGENIE Earth system model, basic usage of the cGENIE Earth system model and required software, plus details of and how to execute and analyse the published model experiments, is given here: http://www.seao2.info/cgenie/pubs/gutjahretal. 2017.txt Further specific details of e.g. using the provided plotting functions to process the model results as per in the paper, configurations for the experiments presented in Extended Data and/or described in the SI, or the raw model output, can be obtained directly from A.R. (andy@seao2.org).

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 δ 13 C (by 3.5‰) and pH (by 0.3 pH units) that varies from 100 to 20,000 yr, with the target pH and δ 13 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 δ 13 C on the right hand y-axis. The lower rows of panels show: (b) maximum emission rate per time interval,

Europe PMC Funders Author Manuscripts
(c) cumulative carbon emission for respective onset phase in EgC (1 Eg = 10 18 g) and (d) average emitted δ 13 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.072. 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 δ 11 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  down-core model-data comparison is carried out (Extended Data Fig. 9). (c-f) Model projected spatial pattern of the surface sedimentary wt% CaCO 3 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.

Europe PMC Funders Author Manuscripts
Europe PMC Funders Author Manuscripts 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)  atmospheric δ 13 C of CO 2 in the model. For wt% CaCO 3 and δ 13 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)  Extended Data Table 1 Key results from individual model runs.
(a) Summary of the main double inversion experiments carried out. The terminology "R07" refers to configurations tying the Site 401 records to the chronostratigraphy of ref. 25, the notation "FE" refers to the 3 He-based age model of ref. 57). Annotation "sm" refers to inversion of analytically smoothed δ 13 C and pH data sets, "rw" to usage of original sample data for double inversions. "HI" and "LO" represent potentially extreme configurations taking into account the boron proxy uncertainty at 95% confidence level. "noW" has silicate (and carbonate) weathering feedbacks disabled. "Corg" denote model configurations that allow removal of excess organic carbon from the surface ocean. Grey shading highlights experiments focussed upon in the main text and plotted in Figure   Extended Data Table 1 Key results from individual model runs.
(a) Summary of the main double inversion experiments carried out. The terminology "R07" refers to configurations tying the Site 401 records to the chronostratigraphy of ref. 25, the notation "FE" refers to the 3 He-based age model of ref. 57). Annotation "sm" refers to inversion of analytically smoothed δ 13 C and pH data sets, "rw" to usage of original sample data for double inversions. "HI" and "LO" represent potentially extreme configurations taking into account the boron proxy uncertainty at 95% confidence level. "noW" has silicate (and carbonate) weathering feedbacks disabled. "Corg" denote model configurations that allow removal of excess organic carbon from the surface ocean. Grey shading highlights experiments focussed upon in the main text and plotted in Figure 3 ("R07sm" in Fig. 3a-f and "R07sm_Corg" in 3i-n.). Note: (1)