Millennial-timescale thermogenic CO_2 release preceding the Paleocene-Eocene Thermal Maximum

Domain: earth_science
Published: 2025-06-30
Editorial outcome: major_revision

Abstract
Geologic records support a short-lived carbon release, known as the pre-onset excursion (POE), shortly before the Paleocene-Eocene Thermal Maximum (PETM; ~56 Ma). However, the source and pace of the POE carbon release and its relationship to the PETM remain unresolved. Here we show a high-temporal-resolution stratigraphic record spanning the POE and PETM from the eastern Tethys Ocean that documents the evolution of surface ocean carbon cycle, redox and eutrophication, confirming the global nature of the POE. Biomarkers extracted from the sedimentary record indicate a smaller environmental perturbation during the POE than that during the PETM in the eastern Tethys Ocean. Earth system modeling constrained by observed δ13C and pH data indicates that the POE was driven by a largely thermogenic CO2 source, likely associated with sill intrusions prior to the main eruption phase of the North Atlantic Igneous Province and possibly biogeochemical feedbacks involving the release of biogenic methane.

Introduction
A holistic understanding of the carbon-climate dynamics of past warming events has important implications for CO2-induced anthropogenic climate change. The Paleocene-Eocene Thermal Maximum (PETM; ~56 Ma1–3) represents the largest disruption of the global carbon cycle in the Cenozoic4, which led to 5–6°C global warming5,6, ocean acidification7, ocean deoxygenation8–11, and intensified tropical cyclones12,13. The prominent 3–6‰ negative carbon isotope excursion (CIE) registered in both terrestrial and marine sections is consistent with major emissions (~2000 to >13,000 Pg C) of 13C-depleted carbon to the atmosphere and/or ocean on a timescale of a few to no more than ~20 kyr14,15. Recent work suggests that the North Atlantic Igneous Province (NAIP) and associated CO2 emissions may have triggered the PETM16–18, followed by carbon sequestration through organic carbon burial19 and silicate weathering20. The PETM was proceeded by a transient warming accompanied by a smaller CIE21—known as the pre-onset excursion (POE) and which is recorded in terrestrial records from the Wyoming Bighorn Basin21 together with only a few shallow marine sections (Atlantic coastal plain, southwest Pacific Ocean, the North Sea and the Pyrenean foreland basins)22–26. The POE is a short-lived warming event that occurred about 38 kyr to >100 kyr27 prior to the PETM onset with an estimated duration of no more than a few centuries22 to millennia27. As an environmental precursor to the PETM, the POE is absent in deep-sea sedimentary records because its short duration may have limited its preservation to surface and shallow water records22. Resolving a global POE signal could be further complicated by bioturbation, sediment mixing, and chemical burndown of deep-sea carbonates7,28, which could only be understood by studying shallow marine and terrestrial sections. The POE warming may represent an early warning signal on the instability of carbon reservoirs and set the stage for a climatic threshold crossing occurred during the PETM. Previous studies suggest that the PETM is modulated by astronomical forcing29–31, and linked with the POE via repeated, catastrophic CO2 release27, such as methane hydrate dissociation22, either as a direct response of the warming or via positive feedback mechanisms. Furthermore, the close timing between the initial stage of the NAIP and the POE suggests that volcanism and magmatism may also serve as a viable trigger32. However, the global extent of the POE, its relationship with the PETM and exact mechanisms that triggered the POE—whether from methane hydrate release, volcanic activity, or orbital drivers—remain debated.

Here we report ultra-high-resolution biogeochemical records from a recently discovered coastal shallow marine section in the eastern Tethys that span both the POE and the PETM (Fig. 1). The Kuzigongsu section (39°45’10”N, 75°17’29”E) is located in the western Xinjiang Uygur Autonomous Region of China, which was covered by the Turan Sea—an arm of the Tethys Ocean during the early Paleogene (Fig. 1). The eastern Tethys was a restricted shallow-water carbonate platform environment33, and a critical site for the formation of warm and saline intermediate water and the burial of organic matter34. Abundant calcareous nannofossils35 and well-preserved organic matter and oyster shells (Fig. S1) allow for an integrated sedimentological, biogeochemical, isotopic, organic geochemical, and global carbon cycle modeling approach to unravel the paleoenvironmental evolution of the eastern Tethys during the POE and PETM, thus filling a critical spatial data gap and advancing knowledge on forcing and recovery mechanisms of ancient hyperthermals.

Results
Astronomically tuned high-resolution PETM and POE records from the understudied eastern Tethys

The PETM has been identified at the Kuzigongsu section through calcareous nannofossil biostratigraphy (the NP9/NP10 boundary) and carbon isotope stratigraphy35. It occurs at 19.9 m (on a depth scale of 0 to 48 meters in Fig. 2) and corresponds to a ~6–8‰ negative carbon isotope excursion (CIE)—among the largest CIEs observed in shallow marine sites4. The CIE magnitude is ~6.3‰ in carbonate35, ~6.0‰ in organic matter, and somewhat amplified in long-chain n-alkanes (~7.8‰), which is likely a result of an enhanced hydrological cycle36 and elevated pCO22,37. The primary δ13Ccarb signal is likely well preserved, based on: (1) the strong covariation between δ13Ccarb and δ13Corg (r2=0.75, p<0.001; Fig. S3); (2) the fact that most δ13Ccarb and δ18Ocarb data plot within the area of primary carbonates38 (Fig. S3), and (3) the presence of only a weak correlation (r2=0.18, p<0.001, Fig. S3) between δ13Ccarb and Mn/Sr, as a strong correlation would indicate diageneticalteration39.

The POE is found at ~8.4 m below the PETM onset within lower nannofossil Zone NP9a35 and occurs in a 1.2-meter-thick interval (10.3 to 11.5 m) characterized by a –1 to –2.5‰ CIE (Fig. 2). Specifically, we observed CIEs of –2.5‰ in carbonate and –2.1‰ in organic matter, but in contrast to the PETM, only ~–1‰ in long-chain n-alkanes. The relatively smaller recorded magnitude in the n-alkane record is likely due to the lack of data at 10.8 m depth where δ13Ccarb and δ13Corg values reach their minima (Fig. 2).

Power spectrum analysis of the detrended magnetic susceptibility (MS) data series shows significant peaks in wavelength at 0.8, 1.2, 1.9, 3, 5, 6.5, and 9.8 m (see Methods, SI and Figs. S3–S5), with the filtered 1.2–1.9 m cycles interpreted as precession signal with an assumed 21 kyr duration and the filtered 5 to 9.8 m cycles as short eccentricity (~100 kyr). Spectral analysis revealed sedimentation rates averaging between 6.0 and 8.3 cm kyr–1 (Fig. S5) and suggests that the durations of the PETM and the POE at our study site are ~127 kyr and ~21 kyr, respectively (age model option 1; see SI and Supplementary Data 2 for details). The PETM and POE are separated by ~144 kyr (±21 kyr). The estimated PETM duration of 127 kyr is shorter than inferred from the deep sea sites (e.g., ~170 kyr from Röhl et al.2 and Zeebe and Lourens3), likely due to incomplete preservation of the entire PETM at Kuzigongsu with a change in lithology that truncates the recovery phase. The POE onset duration of ~7.0 kyr (age model option 1) is similar to, but slightly longer than the 2 to 5.5 kyr estimated by Bowen et al.21. An alternative age model (Option 2) that accounts for the significant drop in wt.% CaCO3 and a likely truncation assumes the filtered 6–10 m cycles represent ~20 kyr precessional signal. This age model option provides a duration of ~39 kyr for the PETM, ~4 kyr for the POE and ~54 kyr between the PETM and POE, which suggests the study site only preserves the PETM onset and the plateau, rather than the recovery (see SI for more discussion). However, due to the uncertainty in the astronomically tuned age model, we assume that the POE onset duration ranges from 500 to 7000 years to cover the full range of reported values in the literature22,27.

Paleoenvironment of the eastern Tethys during the POE and PETM

We use a multi-proxy approach to reconstruct the paleoenvironmental evolution of the eastern Tethys during the POE and PETM (Supplementary Data 1). Our records (Figs. 2, 3) include C/N ratios as indicators of organic matter source, weight percent (wt.%) CaCO3 as a proxy for ocean acidification and detrital dilution, trace element geochemistry for marine nutrient and chemical weathering proxies, organic biomarkers as proxies for marine microbial communities, and mercury content as a possible indicator of the NAIP activity. Together, our new data suggest that the shallow eastern Tethys experienced profound environmental changes, including extreme warmth, eutrophication, and biological turnover. Furthermore, the moderately high sedimentation rates (optimal sedimentation rate fluctuates between ~6–8 cm kyr–1; Fig. S3, S4) at this shallow site (estimated water depth is ~30–50 m based on microfacies analysis and foraminifera indicators33) yield highly expanded records that provide unique details on the relationship between the PETM and the POE. Such details are generally obscured in deep-sea sites because of lower sedimentation rates, dissolution, and bioturbation7.

The section is characterized by a rapid decrease in wt.% CaCO3 from >80 wt.% to near 0 wt.% at ~9 m—a shift which precedes the POE and PETM and may be attributed to significant reduction of carbonate production, detrital dilution, and/or shallow ocean acidification7,40. The sharp decrease (−2.5‰) in oxygen isotopes of marine carbonate values (δ18Ocarb) is consistent with an abrupt and significant warming during the POE and PETM (Fig. 2b), though the magnitude is likely an artifact of diagenetic overprinting. Additionally, a portion of the δ18Ocarb decrease may represent a decline in local salinity as the δ18Osw at epicontinental sites can be strongly influenced by freshwater input from surrounding continents5. Clumped isotope data from a well-preserved oyster specimen (at 29.8 m; Fig. S1) indicate that the eastern Tethys surface water temperature was around 32.5 ± 1.5°C (1σ) at the recovery phase of the PETM (Fig. 2). This estimate is similar to our independent temperature estimate of 30.6 ± 4.5°C (1σ) based on the TEX86H proxy41 for the sample at the same depth. However, the thermal maturity is relatively high for this section and the cyclized isoGDGTs abundance is low, preventing us from obtaining a high-resolution and precise TEX86H temperature record at the site (Fig. 2).

In the organic matter fraction, peak TOC and C/N ratios coincide with the lowest δ13Corg values during the PETM, suggesting increased terrestrial organic matter input at the study site, a likely consequence of intensified continental weathering and/or higher terrestrial primary production42. The inferred increase in terrestrial weathering is supported by the higher values of Ti/Al and K/Al ratios11. Elevated C29 hopane ββ/(αβ+βα+αβ) ratios (average=0.3) during the PETM indicate increased input of fresh organic matter either due to higher primary productivity or increased flux of fresh terrestrial organic matter into the basin (Fig. 3e). Lower C29 hopane ββ/(αβ+βα+αβ) ratios (average=0.1) in the pre- and post-PETM samples suggest relatively low primary production in the surface waters with background input of reworked and more mature organic matter from the surrounding continents42. Similarly, Crenarchaeol/(Crenarchaeol + isoGDGT-0) ratios range from 0.1 to 1.0, with a significant decrease during the POE and PETM. Crenarchaeol (with four cyclopentane rings and one cyclohexane ring) is considered as a biomarker for Thaumarchaeota43. The lower Cren/(Cren+isoGDGT-0) ratios during the POE and PETM therefore likely reflect a reduction in marine Thaumarchaeota, which may be attributed to warmer surface ocean temperature and lower dissolved oxygen concentration44. The occurrence of 2-methylhopanes (2-MeHop) in the PETM interval indicates a transient perturbation of surface ocean characteristics (Fig. 3). The C29 2-MeHop Index, calculated as 100 × (C29 2-MeHop)/(C29 2-MeHop+C29 Hop)45, ranges from ~0–38% with two prominent peaks, at 20.9 m and 28.5 m respectively, corresponding to the peak values of TOC and C/N ratios. Several studies reported that the occurrence of 2-MeHop in the sedimentary record can be viewed as indicators of stress responses to the capacity of microbial respiration under hypoxia46, nitrogen fixation47, increased productivity48, and changes in pH49, corroborating the interpretations of elevated primary productivity discussed above. Furthermore, the anomalously high C29 2-MeHop Index during the PETM may be attributed to marine nitrogen cycle perturbation as a result of biogeochemical changes. This is similar to observations of other major carbon cycle perturbations of the Phanerozoic, such as the end-Permian mass extinction event50, the end-Triassic extinction event51, and the Mesozoic Oceanic Anoxic Events52,53.

Ocean deoxygenation may have been enhanced by increased primary productivity from elevated nutrient input due to enhanced terrestrial weathering. This suggestion is supported by negative Mn* values (Eq. 1) from the POE to the PETM (Fig. 3), which are associated with more reducing conditions due to significant redox-related changes in the solubility of Fe and Mn54.

Mn*=log[(Mn sample/Mn shales)/(Fe sample/Fe shales)] (1)

The values used for the Mnshales and Feshales are 600 and 46,150 ppm, respectively55. Furthermore, the inferred surface ocean deoxygenation is consistent with elevated V/Al ratios over the same interval (Fig. 3) because V ions (+4 and +5 valence) are closely coupled with the redox cycle of Mn56. Widespread deoxygenation is well documented in many ocean basins across the globe during the PETM8,57, including the North Sea58, the Arctic Ocean59, the Atlantic and Caribbean60,61, and the northwestern Tethyan margins62. However, no significant changes in these redox indicators were observed across the POE63, suggesting relatively stable redox conditions in the eastern Tethys at this time.

Mercury content (or Hg concentration normalized as a ratio to organic carbon content—Hg/TOC) has been used as a signal of NAIP activity by several previous studies18,64. Our site exhibits two prominent Hg/TOC peaks that preceded the onset of the POE (~11 kyr) and the PETM (~26 kyr) (Fig. 2), supporting a possible link between Hg source and the 13C-depleted carbon source. However, because of the overall low Hg concentrations at the study site, establishing a direct link between the NAIP and the Hg peaks is not straightforward. Low Hg is likely due to dilution by carbonate and detrital input, the long distance of the site relative to Hg source, and/or Hg transport via oceanic waters rather than global atmospheric transport18. Increased Hg concentrations across the POE and PETM compared to background values suggest that multiple possible sources and processes may have been at play in addition to the NAIP activity. For example, variations in Hg concentrations in the sedimentary record can be caused by changes in river runoff, weathering, transport of terrestrial materials, primary productivity, source of organic matter, and post-depositional processes (e.g., diagenesis and dissolution)18, which could become more important at the study site because of its restricted carbonate platform setting65. Deoxygenation and changes in organic matter preservation and transport cannot fully account for the excess Hg as shown by the steeper Hg gradient to TOC within the PETM and POE interval at our site (Fig. S6). Moreover, Hg fluxes associated with wildfire (e.g., Arctic region66, northeastern US margin67, and England68) may have been far less than the Hg fluxes associated with a large igneous province source event69, and therefore cannot provide sufficient Hg into the study site. Principal component analysis (PCA) suggests that Hg is most closely related to C/N ratios (higher C/N ratio indicates more terrestrial organic source) and δ13Corg during the PETM, which reflect changes in source of organic matter and 13C-depleted CO2 emissions (Fig. S6). The C/N ratio exhibits no significant change across the POE, suggesting the increase in Hg and Hg/TOC ratio is unrelated to changes in source of organic matter. On the other hand, C/N ratio shows a large increase across the PETM, which indicates that changes in source of organic matter may have contributed to the increased Hg concentrations. These potential processes do not preclude volcanic involvement, however, especially via more complex pathways than simple atmospheric loading and deposition70. Despite these potential complex sources of Hg, we cannot completely exclude direct and indirect involvement of the NAIP in driving the Hg changes in the study section71. For example, the NAIP was active as early as 62 Ma71, and its peak activity may have encompassed both the POE and the PETM72–75. A negative shift in 187Os/188Os ratios has been observed prior to the PETM in several sites globally32,71,76,77, lending support to the occurrence of large igneous province activity prior to the PETM. Furthermore, hydrothermal vent complexes in the northeast Atlantic region78,79 further support that the NAIP activity can at least partially explain the observed Hg records.

Thermogenic CO2 emissions associated with NAIP activity during the POE

The PETM carbon emission history has been extensively modeled in the past, with estimated carbon emission rates ranging from 0.3 to 1.7 Pg C yr-1 for a CIE onset duration from ~3000 to ~20,000 years and cumulative amount of carbon added ranging from ~2500 to >13,000 Pg C14,15,21,79–81. Because the carbon emission history preceding the PETM has not been systematically quantified in an Earth system model and very little is yet known about the CO2 source during this time22, we then focused our model analysis on the POE (Table S1 and Fig. 4). Our new high-resolution geochemical data, together with an orbitally tuned astronomical age model, provide a unique opportunity to assess the effects of CO2 emissions during the POE.

We quantify carbon emissions over the POE using a data assimilation approach that considers paired δ13C-pH variation across the POE within an Earth system model of intermediate complexity cGENIE, following the approach detailed in Gutjahr et al.15. In this, changes with time in annual global mean surface ocean pH (derived from δ11B proxy data from the Mid-Atlantic Coastal Plain with a change of ~−0.1 to −0.3 pH units22) constrain the emission rate of CO2 to the atmosphere. Similarly, the change with time in observed δ13C of annual global mean surface ocean DIC (δ13CDIC) (reconstructed by applying an anomaly derived from the δ13C data of a global compilation; Fig. S7) refines the δ13C value of the (pH-constrained) CO2 emissions. The novelty of this approach is that it offers a unique solution of the mean δ13Csource without having to make a specific assumption about the carbon source (e.g., compare with Cui et al.14; see Methods and SI for detailed model results and sensitivity tests). To account for the uncertainty in the POE onset duration, we place our records on four different age models, including age model options 1 and 2 from this study, an age model from the Bighorn Basin based on Bowen et al.21 and an assumed age of 500 years based on Babila et al.22 (a summary of our model results and sensitivity analyses for the POE is listed in Table S1). (Fig. 4).

The flux-weighted δ13Csource values across the entire emission duration vary between −30.8 and −44.5‰ for the four age models used in our simulations with the minimum change in pH suggested by Babila et al.22, consistent with a thermogenic CO2 source82 (−30 to −65‰; Fig. 4a–d and Table S1). Longer POE duration (e.g., Age 1 associated with ~7000 year POE onset) necessitates lower flux-weighted δ13Csource values (−44.5‰) over the entire emission interval (Fig. 4a) at slower emission rate (~0.2 Pg C yr−1). We note that the δ13Csource values become progressively lower from the POE onset, likely resulting from a faster rate of change toward its minimum values in the δ13C forcing. This may represent a shift from thermogenic methane to biogenic methane (~−34 to <−70‰83,84) emissions during the development of the POE. The average carbon emission rate over the entire emission period ranges from 0.2 to 1.3 Pg C yr−1 (Fig. 4e–h), comparable to those estimated for the PETM from sill-degassed CO2 and thermogenic methane (0.2 to 0.5 Pg C yr−1 from Jones et al.16; 0.6 Pg C yr−1 from Frieling et al.79). Larger magnitude of pH changes (e.g., ΔpH = ~0.2 to 0.3) yield overall larger average peak CO2 emission flux (2.9 Pg C yr−1) and higher average δ13Csource values (−19.5‰) (Table S1), still consistent with largely thermogenic methane source. The pH change for the POE has been documented at only a single location, using a novel approach to measuring boron isotopes (δ11B), and therefore has a high degree of uncertainty22. Considering the smaller magnitude of δ13C excursion, the smaller degree of warming, its shorter duration, and the minor ecological responses, the changes in pH during the POE are unlikely to exceed those during the PETM (ΔpH = ~0.3)85. A higher average carbon emission rate is associated with shorter POE onset duration (Fig. 4e), which represents a combined impact of the imposed ΔpH forcing and age models used (Fig. S7).

The cumulative CO2 emission during the POE ranges from ~1030 to 1765 Pg C (Fig. 4i–l), with peak pCO2 reaching ~1180 to ~1220 ppm—a rise of ~350 to ~390 ppm above ~830 ppm (Fig. S8). The modeled cumulative carbon emitted during the POE falls within the range of the 400 to 1600 Pg C suggested by Babila et al.22 using similar ΔpH. However, if the actual ΔpH was at the lower end (lower than ~0.1), it is more likely that the carbon source was primarily biogenic methane. Associated with the diagnosed carbon emissions is a modelled global sea surface temperature rise (ΔT) of ~1.1 to ~1.3°C (Fig. S8). Although the paleotemperature history of the POE is currently poorly known, existing Mg/Ca ratios of planktonic foraminifera from the mid-Atlantic coastal plain suggest that the surface ocean temperature increase was ~ 2°C with an uncertainty of ±1°C due to salinity variations22, consistent with our modeled temperature changes within uncertainty. The POE warming may also help explain the observed increase in warm-water coccolithophore taxa in the eastern Tethys35.

Thermogenic CO2 related to the NAIP activities may have been the dominant carbon sources during the POE via contact metamorphism by intrusive activity through hydrothermal vent complexes78. It should be noted that mantle convection models suggest that a peak NAIP carbon emission flux at ~0.5 Pg C yr−1 could occur between 1 and 20 kyr16, comparable to those simulated in our inversion experiments, despite the geochronology of the NAIP continental flood basalt sequences being not very well constrained71. It is also important to note that a caveat of cGENIE in interpreting our results is the lack of terrestrial biosphere and potential changes in orbital forcing, which could impact the climate responses and lead to uncertainties in carbon emission estimates. Although this study provides a range of estimates on the carbon source and emission flux during the POE, more precise δ11B-based global surface pH records, detailed history of the sill intrusion of the NAIP, sea surface temperature records from across different latitudes, and better-constrained geochronology of the NAIP activity are clearly needed to reduce the uncertainty of the estimated thermogenic carbon emission fluxes from the NAIP.

The evolution of mean core-top carbonate (CaCO3) with time in the model exhibits a smaller magnitude of δ13C decrease for simulations with bioturbation turned on compared to those without bioturbation (Fig. S9). Similarly, core-top wt.% CaCO3 also exhibit smaller degree of dissolution for experiments with bioturbation on (Fig. S9). Longer experiment duration allows for a larger CIE magnitude regardless of whether bioturbation is on. This is due to the combined effects of bioturbation and dissolution as a result of the cumulative carbon emission (Fig. S9 and Fig. 4i–l), supported by a comparable Eocene hyperthermal event86. These experiments support the inference that short POE onset duration (less than millennial timescale) and bioturbation are the main causes of the lack of POE signal in the deep-sea sedimentary records.

In conclusion, we report astronomically tuned, ultrahigh-resolution PETM and POE stratigraphic records from a recently discovered site in the eastern Tethys. Geochemical proxies based on carbonate, bulk organic matter, and biomarkers suggest that the eastern Tethys experienced profound carbon cycle perturbations during the POE and PETM. Our integrated stratigraphic data and Earth system modeling together suggest that the millennial time-scale POE may be attributed to mainly thermogenic CO2 emission associated with sill intrusion prior to the main eruption phase of the NAIP, with contributions from amplifying feedbacks such as biogenic methane release. Furthermore, our findings predict substantial carbon fluxes driving the POE (averaging 0.2 to 1.3 Pg C yr-1), which could be tested by refined geochronological investigations of potential sources such as the NAIP. The POE may have set the stage for the ecosystem threshold crossing and the extreme carbon cycle disruption occurred during the PETM.

Methods
Cyclostratigraphy and astronomically tuned age model based on magnetic susceptibility measurements and time series analysis

A total of 480 samples at 10 cm intervals spanning both the POE and the PETM weighing 4 to 8 grams were measured for bulk mass-normalized magnetic susceptibility (MS or χ) using KLY-4S Kappabridge after being crushed in a copper rock hammer and placed in a 2×2×2 cm3 cubic plastic holder. The MS measurements were conducted at the Paleomagnetism and Environmental Magnetism Laboratory at the China University of Geoscience (Beijing). Measurements were made at room temperature with an applied field amplitude of 200 A/m and frequency of 976 Hz. Each measurement is corrected for the contribution of the plastic sample holder. Each sample was measured three times, with the average value corrected by mass to obtain χ in units of m3 kg-1. Relative standard deviations between the three runs were smaller than 0.5%.

Time-series analysis was conducted using MS data with the open-source software Acycle V2.487 because MS measures the magnetic mineral concentration, and is considered as a proxy for detrital fluxes from land to the ocean88. The MS data series was first detrended by subtracting a 40 m “loess” trend (locally estimated scatterplot smoothing, a non-parametric method for a series of data smoothing with a default window size of 35%) to remove non-periodic or high-amplitude long-term trends following the procedures described in Li et al.89. The multi-taper method (MTM)90 with 2π tapers was used to estimate the spectrum for the detrended MS series and confidence levels (mean, 90%, 95%, and 99%) were provided to test against robust first-order autoregressive model AR(1) red noise in order to reveal the MS series’ dominant wavelength. The evolutionary power spectra were calculated with “Evolutionary Spectral Analysis” function in Acycle with a sliding window of 10 m and a step of 0.1 m to identify any secular trend in dominant frequencies, which may be attributed to variations in sedimentation rates. The time scale optimization (TimeOpt; Meyers91) and correlation coefficient (COCO; Li et al.87) methods were used to identify the optimal sedimentation rate using Acycle’s “COCO” and “TimeOpt” functions, which use 2000 Monte Carlo statistical simulations to test the null hypothesis of no orbital forcing. The evolutionary versions of COCO and TimeOpt functions (i.e., eCOCO and eTimeOpt) were used to track changes in sedimentation rates. In addition, the “Spectral Moments” function was used to estimate variable sedimentation rates based on a periodogram with two spectral moments: evolutionary mean frequency (μf) and evolutionary bandwidth (B) (Fig. S3, S4). Subsequently, “Dynamic Filtering” function was used to apply dynamic filtering and isolate interpreted precession cycles from the MS data series. Since the power of long-term cycles (i.e., short eccentricity cycles) may have muted the manifestation of precession cycles in the evolutive harmonic analysis (EHA), we remove the > 4 m cycles that may be associated with eccentricity cycles to reveal precession-related cycles as the most prominent signal in the EHA spectrogram (Fig. S5). The significant power of the interpreted precession cycles in the EHA spectrogram allows us to effectively isolate this signal from EHA (Fig. S5). We then use the precession cycles to construct an astrochronological timescale for the study interval. Analyses of TimeOpt and COCO indicate alternation of optimal sedimentation rates (i.e., 6.0 cm kyr-1 and 8.3 cm kyr-1) (Fig. S3). Spectral Moments, eTimeOpt and eCOCO together suggest the estimated sedimentation rate ranges from 4.2 to 10.6 cm kyr-1 with increased sedimentation rate during the PETM body (Fig. S3–S5).

Stable carbon isotopes of bulk organic matter and wt.% CaCO3

HCl-treated carbonate-free powders were measured for total organic carbon (TOC) and total nitrogen (TN) concentrations on a Vario EL-III elemental analyzer, and the δ13Corg analyses were made using a thermo DELTA plus XL mass spectrometer at State Key Laboratory of Organic Geochemistry, Guangzhou Institute of Geochemistry, Chinese Academy of Sciences. Three reference materials were used to monitor the measurement of carbon isotopic ratio of bulk organic carbon, which included black carbon (−22.43‰), Urea#1 (−34.13‰), and Urea#2 (−8.02‰). Precision based on repeated measurement of these three standards were 0.12‰, 0.08‰, and 0.09‰, respectively. δ13Corg values were reported in VPDB and analytical precision was better than ±0.1‰ based on replicate analyses of the standards processed with each batch of samples. Weight percent (wt.%) CaCO3 was measured using a modified acid soluble weight-loss method92.

Carbonate clumped isotope geochemistry

The carbonate clumped isotope thermometer is based on the thermodynamic stability of C–O bonds at varying temperature, in which “clumping” of the rare, heavy isotopes of carbon and oxygen (13C and 18O) occurs more frequently at lower temperatures93. The excess occurrence of the 13C18O16O isotopologue of CO2 relative to a stochastic distribution of the heavy isotopes among all CO2 molecules is referred to as the mass47 anomaly and notated as Δ47, in which Δ47 = 47R × 1000 where 47R=[13C16O18O+12C17O18O+13C17O2]/[12C16O2] and * denotes a stochastic distribution of isotopes. Clumped isotope thermometry presents a significant innovation over oxygen isotope-based thermometry because the temperature estimate is independent of the bulk isotopic composition, and thus requires no assumptions about δ18Ocarb or δ18Owater. This mineral formation temperature can be used to calculate δ18O of ancient waters when paired with δ18Ocarb values of the same sample, which is measured concurrently with Δ47. Carbonate clumped isotope measurements of one Eocene fossil oyster (Crassostrea sp.) and one modern oyster specimen (Crassostrea hongkongensis) collected from northern South China Sea (21°42'7.89” N, 111°55'44.61” E) in 2022 were made at the Pennsylvania State University in April 2022 (see SI). Approximately 8 mg of pure carbonate powder was digested in a 105% phosphoric acid common acid bath at 90°C to yield CO2. Evolved CO2 was passed through a Protium Isotope Batch Extraction (IBEX) carbonate preparation line to purify the sample gas. The gas is passed through a cryogenic trap to separate CO2 from water, a silver wool-packed borosilicate column to trap sulfides, and a gas chromatography column packed with Poropak to separate CO2 from other compounds with a He carrier gas. The purified CO2 gas is once more frozen into a cryogenic trap before being frozen into a microvolume, and passed through a polished nickel capillary to the MAT 253 Plus bellows. Purified CO2 sample gas was analyzed on a Thermo MAT253 Plus dual inlet IRMS relative to an Oztech working gas.

Δ47 values versus the working gas were projected to the Intercarb Carbon Dioxide Equilibrium Scale94 (I-CDES) using a carbonate standard-based empirical transfer function. ETH 1, 2, 3, and 4 were measured to build the reference frame and for interlaboratory comparison, and IAEA-C2 and Carrara Marble were treated as unknowns. Individual replicates were averaged to create final sample Δ47 values and reported with a 95% confidence interval. Temperatures were calculated using the T-Δ47 calibration of Anderson et al.95. The average measured Δ47 value for the oyster fossil is 0.573 ± 0.011 (2σ), while the Δ47 value for the modern oyster specimen is 0.604 ± 0.028 (2σ). The calculated sea surface temperature in the eastern Tethys based on early Eocene oyster fossil is 32.5 ± 3.9°C (2σ). The calculated modern sea surface temperature based on modern oyster specimen is 21.6 ± 8.7°C (2σ), falling in the range of the observed average annual sea surface temperature (24.1 ± 5.6°C) in northern South China Sea in 2022.

Biomarker and stable carbon isotopes of long-chain n-alkanes

Around 11 grams of dried and powdered sample were extracted for their biomarker content using a microwave system (Milestone Ethos EX) and using 20 ml of a dichloromethane and methanol mixture (9:1). The total lipid extract was separated using silica flash chromatography and elution with hexane: DCM (9:1) for the apolar and DCM: MeOH (2:1) for the polar fraction. The apolar fractions were characterized on a Thermo Scientific ISQ single quadrupole mass spectrometer (MS) coupled to a gas chromatograph (GC). Compounds were separated using a fused silica column (50 m × 0.32 mm) with a ZB1 stationary phase and helium as the carrier gas. The GC was programmed for: injection at 70°C (1 min hold), ramp to 130°C at 20°C/min, followed by a ramp to 300°C at 4°C/min (20 min hold). The MS continuously scanned between m/z 650-50. The apolar fractions were subsequently analyzed using an Isoprime 100 combustion isotope ratio mass spectrometer, coupled to an Agilent GC, to determine the δ13C of the long-chain n-alkanes. We used the same type of column and temperature program as used for the GC-MS analyses. Samples were measured in duplicate on the GC-C-IRMS, and the average is reported here. An in-house CO2 reference gas was used to calculate compound specific δ13C values relative to Vienna Pee Dee Belemnite (VPDB). δ13C values of the C29 n-alkane are not reported here due to possible co-elution with other lipids. All biomarker and stable carbon isotopes of long-chain n-alkane analyses were performed at the University of Bristol.

Methods for GDGTs

Polar fractions were filtered through a 0.45 µm filter at the university of Bristol. The filtered polar fractions were redissolved in hexane: isopropanol (99:1) and analyzed using a high-pressure liquid chromatography atmospheric pressure chemical ionization mass spectrometer for their GDGT distribution. We used two ultra-high performance liquid chromatography silica columns to separate compounds, following Hopmans et al. (2016)96, and analyses were performed in selective ion monitoring (SIM) mode.

The thermal maturity of the organic matter in this section was estimated using the hopane isomerisation index: C29 ββ/(αβ+βα+ββ)97. The results indicate that the thermal maturity changes across the section, but the C29 ββ/(αβ+βα+ββ) ratio is consistently below 0.4 (Fig. 3). This is indicative for an elevated thermal maturity, but well below the oil window. Although this level of thermal maturity will not affect apolar compounds like hopanoids or n-alkanes, it is likely to impact more labile biomarkers such as glycerol dialkyl glycerol tetraethers (GDGTs)98. We determined the GDGT distribution in all samples. As expected with this level of thermal maturity, GDGT concentrations were low and, in most samples, branched (br) GDGTs were absent, as were isoprenoidal (iso)GDGTs containing cyclopentane rings. However, a few samples did have isoprenoidal (iso)GDGTs with cyclopentane rings. This includes the sample at depth 29.8 m that hosts the well-preserved oyster shell fossil. This sample has a TEX86 value of 0.76, which results in an SST of 30.6 ± 4.5°C using the TEX86H calibration99. Although we treat this estimate cautious as thermal maturity might have impacted the GDGTs distribution, this TEX86-based SST is consistent with the clumped isotope data from well-preserved oyster shell fossils from the same sample, adding confidence that we are able to constrain the SSTs at this site during the recovery phase of the PETM.

Earth system modeling

The carbon-centric Grid Enabled Integrated Earth system model (cGENIE) is an intermediate complexity climate model that couples a 3D ocean (36×36 grid, 16 levels) with a 2D atmosphere that has the capability to track biogeochemical cycling of elements, stable carbon isotopes, marine sediments, and continental weathering15,100. Bathymetry, paleogeography, planetary albedo, and wind fields are configured for the late Paleocene-early Eocene with the same initial and boundary conditions as Gutjahr et al.15. For example, the δ13C value of late Paleocene-early Eocene atmospheric CO2 (δ13CCO2) is set as ~‒5‰, and the atmospheric pCO2 is set as ~830 ppmv. The moderately high pCO2 allows for a small buildup of sea ice (0.5%) in the northern polar regions. We then run a number of ‘double inversion’ experiments in which δ13C of surface ocean dissolved inorganic carbon (δ13CDIC) and surface ocean pH22 are used as the two data assimilation constraints for the POE. The δ13CDIC forcing is based on the high-temporal-resolution δ13Ccarb data from the shallow Tethys Kuzigongsu section using astronomically tuned age models. For our inversion experiments, the model was first spun up for 20 kyr to establish the basic ocean circulation and climatic state under published late Paleocene-early Eocene boundary conditions, including paleogeography and paleobathymetry101,102. This is followed by an open-system spin-up of 200 kyr to allow the long-term δ13C cycle to reach balance. A range of inversion experiments were carried out (Table S1; Fig. S7–9). Although uncertainty exists for pre-PETM δ11B, the surface ocean pH at the end of the open-system spin up is 7.75, same as those used in Gutjahr et al. (2017)15, which is adapted as the initial surface ocean pH forcing in the “double inversion” experiment.

First, the “double-inversion” modeling takes the observed pH data, which constrain the flux and magnitude of CO2 emissions, and the observed δ13C values of the dissolved inorganic carbon of the surface ocean, which simultaneously determines the source of the emitted carbon by computing the δ13C values of the carbon source. At each model timestep, a pulse of CO2 is emitted to the atmosphere at a given rate if the δ13C value is lower than the previous timestep, and the modeled surface DIC δ13C values and the observed δ13C values at the Kuzigongsu section are compared. If the current modeled surface DIC δ13C value is higher than the data value, the δ13C value of the emitted CO2 is assigned a value of –100‰. In contrast, if the current modeled surface DIC δ13C value is lower than the data value, the δ13C value of the emitted CO2 is assigned a value of 0‰. δ13C values of the emitted CO2 between –100‰ and 0‰ can be achieved by binning the emission fluxes in time and averaging flux-weighted δ13C values. Justification for the choice of these end-member δ13C values of the emitted CO2 is provided in Gutjahr et al.15. During the experiments, cGENIE continually adjusts the rate and δ13C value of emitted CO2 into the atmosphere in order to simultaneously reproduce the two proxy records as a function of time. In these experiments, we assume that the POE onset occurred as a linear decline in both δ13C and pH simultaneously (Fig. 4a, b). We use the same “double-inversion” methodology in both the main experiments and the sensitivity experiments, both starting from the same open-system spin-up state (Table S1).

Sensitivity experiments and analyses

We carried out sensitivity experiments to explore the importance of the duration of the POE onset (~7000, ~1600, ~850, and ~500 years based on age model option 1, age model option 2, Bowen et al. (2016), and Babila et al. (2022), respectively) using a global compilation of marine carbonate δ13Ccarb records (Table S1; Fig. S7). We also tested the effect of larger pH decrease (i.e., ‒0.24 and ‒0.32 pH unit) in combination with each of the four assumed age model (Table S1). Additionally, we test the role of bioturbation on the carbon isotope excursion magnitude of core-top carbonates (Fig. S9).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The paper of Jiang et al. presented new data from a section in Xinjiang province in China across the POE and PETM. The dataset is extensive, new and really exciting, so is therefore worthy of publication in Nature Communications. However, there are a few areas where the text and discussion should be improved prior to publication, so I recommend moderate revisions.\nWith best regards, Morgan Jones\nImportant Comments:\nMercury signal.\nThe main issue I have with the paper is the interpretation of the mercury signal. All of the samples presented here have Hg <20 ppb and similarly low TOC contents, which means that analytical errors are going to be enhanced. On the machines I’ve used, Hg contents below 10 ppb have errors comparable to the concentrations, and TOC <0.2 wt% has been shown as the broad cut off for using Hg/TOC as a proxy (e.g. Grasby et al., 2019). Therefore, a lot more care needs to be taken before we can arrive at “originating from the NAIP” (lines 40.42). Hg/TOC values are almost all <30 ppb/wt%, which is well below the average for shales. Even the enriched sections during the POE and PETM are below this mean value from current data sets (Grasby et al., 2019). The authors need to address other possible avenues for disturbances to the Hg cycle. The authors mention and rule out wildfires as a possible source, but the main one that jumps out to me from their dataset is a possible change in organic matter source.\nThere are several indicators of enhanced terrestrial runoff in this dataset (e.g. line 124), including a negative δ18O excursion and elevated C/N ratios (what does the BIT index of the GDGTs look like?). Mercury uptake into terrestrial plant biomass is inherently different to marine uptake, leading to different initial Hg/TOC ratios in each organic matter pool. Given the low Hg and TOC contents, how can the authors be sure that the signal they are observing is indeed due to changes in NAIP activity? There seem to be two distinct Hg/TOC populations in Figure S7, do the population with 10-16 ppb Hg and 0.5-1.5 wt% TOC also have high C/N ratios?\nMoreover, the authors interpret the POE as being the main warming driven by intrusive activity. Volcanic and thermogenic emissions are near impossible to differentiate in the far field, but the limited near-field dataset that we do have suggests that the submarine emissions from contact metamorphism around intrusions leads to much of the Hg being scrubbed into the local water bodies, rather than released to the atmosphere (e.g. Jones et al., 2019a). Why then would the POE have the larger relative Hg/TOC excursion at Kuzigongsu? I’m not saying that the signal here isn’t due to NAIP activity, but more nuanced discussion is needed in the interpretation of this dataset.\nNAIP Emplacement.\nNorth Atlantic Igneous Province (NAIP) activity has been oversimplified in this manuscript. The province was active as early as 62 Ma, and while the main acme of activity did encompass the PETM, it is not (yet) constrained to as fine a window as 56-55.6 Ma (lines 162-163). This interval is based on the start of the crystallisation of the Skaergaard intrusion as it was buried by 5-6 km of lava (Larsen & Tegner, 2006) in East Greenland. The upper constraint is the eruption of the Gronau tuff, which is linked to Ash -17 in Denmark (Storey et al., 2007a). Once corrected to the most recent Fish Canyon Tuff age estimate, this layer is dated to 55.48 ±0.12 Ma, which gives us a ~0.5 Myr window where much of the East Greenland basalts were emplaced. However, this is just one part of a much bigger province, and the Skaergaard magma chamber was intruded into older basalts. Other parts of the province, particularly places like the Faroe Islands and the Vøring plateau, are not nearly so well constrained, and it is unclear at the moment whether this spike in activity in East Greenland is indicative of the province as a whole, or more localised activity.\nFurthermore, I strongly disagree with argument in lines 237-238. The study of Tian and Buck (2022) does not mention the NAIP once, and focuses on how the Deccan Traps and the Colombia river basalts were emplaced. LIP emplacements and the subsequent climatic impacts are different between individual LIPs (see Jones et al., 2016), and other emplacement models exist for other LIPs (e.g. Burgess et al., 2017 suggesting extrusive-intrusive-extrusive for the Siberian Traps). We know from field evidence that there is considerable NAIP volcanism prior to 56.5 Ma and after 55.6 Ma. The intrusive activity is less well temporally constrained, but the only two hydrothermal vent complexes that have been sampled and analysed contain the PETM onset within the vent infill (Frieling et al., 2016; Planke et al., 2022; Berndt et al., in review). There may have been a pulse in intrusive activity at the POE, but almost certainly also at the PETM. Many studies argue, quite convincingly, that magmatic carbon fluxes are insufficient to induce global warming on their own (e.g. Jones et al., 2019b). This section needs to be rewritten to be more nuanced in terms of what constitutes NAIP activity, how (if at all) it is possible to distinguish between volcanic and thermogenic emissions at such distance, and whether such definitive conclusions can be derived from a distal locality.\nMinor comments:\nTitle: The title can be improved, as the PETM is by definition a period of climate instability. The strength of this paper is the well preserved POE and its relation to the PETM, so I suggest rewording to draw attention to what is new and exciting here.\nAbstract: Mercury concentrations are not elevated, <20 ppb throughout the entire section is Hg-starved.\nLine 41: The NAIP was not one eruption. Reword to “originating from North Atlantic Igneous Province activity” or similar.\nLine 72 (& Line 27 in SI): The coordinates for the Kuzigongsu section point to a warehouse in Wuqia county on Google Earth. Please check the coordinates.\nLines 87-88: “The onset of the PETM is defined as the most rapid decline in δ13C values in all studied substrates at 19.9 m” reads like this is the biggest excursion in this section, where I suspect you mean it is the largest of all PETM sections. I suggest you reword to “the onset of the PETM at 19.9 m displays a rapid decline in δ13C values, among the largest CIEs in studied PETM substrates (McInerney & Wing, 2011)”.\nLines 94-95: Please explain why these are deemed “more likely” and “less likely”? Is it to do with the way the age models are set up, or is it that age model 2 predicting a 39 kyr PETM duration goes against a mountain of evidence from other localities?\nLines 116-117: Pleas add oxford commas.\nLines 162-163: 9 million km3 is the total volume of the NAIP emplaced between 62-54 Ma. The percentage that was emplaced between 56.0 and 55.6 Ma is going to be less than this. Suggest rewording.\nLines 185-186: The current best estimates for the onset of the PETM are either 55.93 Ma (Westerhold et al., 2017) or 56.01 Ma (Zeebe and Lourens, 2019). The Charles el al. (2011) paper provides an important ash marker horizon of 55.785 Ma, but this is during the recovery of the CIE. The shape of the CIE in Svalbard also makes it slightly difficult to pin down where in the PETM this marker horizon sits. Citations are needed for this sentence (e.g. Storey et al. 2007a; 2007b; Wilkinson et al., 2017).\nLine 226: How does these model results compare to the findings of Gutjahr et al.?\nLines 235 & 284: Activity is only classified volcanic at the surface. Intrusive activity is magmatism.\nLine 266: Also cite Papadomanolaki et al. (2022) here.\nLine 284: Reword to “as a result of intrusive and/or eruptive activity from a major large igneous province.”\nAuthor contributions: It is not clear what the last four authors contributed to the manuscript.\nFigure 5: The colour palate used for this figure is not very friendly towards the colour blind, suggest changing to colours with greater contrast (see Crameri et al., 2020, for useful tips).\nSupplementary lines 139-154: What were the BIT on analysed samples? Does this tell you anything about the amount of terrestrial organic matter input?\nSupplementary lines 248-251: A value of 2600 kg m3 seems a very high estimate for a density of shale. Are your samples that metamorphosed? We measured the density of the PETM clay in Denmark (paper in prep.) and got a value of 1400 kg m3. Admittedly this is an unlithified section, but if the value used here is correct then the section may have undergone significant compaction post-deposition. Some text addressing this, or better yet a measurement of the actual shale density, would be useful here.\nSupplementary line 254: Table S4 leaves out several recently studied sections that have much thicker complete sections than this, including Fur, Svalbard, the North Sea, and the Norwegian continental margin (e.g. Jin et al., 2022; Jones et al., 2019a; Planke et al., 2022). Perhaps tone down the language a little to say it is one of the thicker PETM sequences, or that it is one of the thickest carbonate-rich sections?\nSupplementary line 287: Correct typo to “performed”\nReferences mentioned:\nBurgess, S.D., Muirhead, J.D., Bowring, S.A., 2017. Initial pulse of Siberian Traps sills as the trigger of the end-Permian mass extinction. Nature Communications 8, 164. doi:10.1038/s41467-017-00083-9.\nCrameri, F., Shephard, G.E., Heron, P.J., 2020. The misuse of colour in science communication. Nature communications 11, 5444.\nFrieling, J., Svensen, H.H., Planke, S., Cramwinckel, M.J., Selnes, H., Sluijs, A., 2016. Thermogenic methane release as a cause for the long duration of the PETM. Proceedings of the National Academy of Sciences USA 113, 12059-12064.\nGrasby, S.E., Them II, T.R., Chen, Z., Yin, R., Ardakani, O.H., 2019. Mercury as a proxy for volcanic emissions in the geologic record. Earth Science Reviews 196, 102880.\nJin, S., Kemp, D.B., Jollet, D.W., Vieira, M., Zachos, J.C., Huang, C., Li, M., Chen, W., 2022. Large-scale, astronomically paced sediment input to the North Sea Basin during the Paleocene Eocene Thermal Maximum. Earth and Planetary Science Letters 579, 117340.\nJones, M.T., Jerram, D.A., Svensen, H.H., Grove, C., 2016. The effects of large igneous provinces on the global carbon and sulphur cycles. Palaeogeography Palaeoclimatology Palaeoecology 441, 4-21.\nJones, M.T., Percival, L.M.E., Stokke, E.W., Frieling, J., Mather, T.A., Riber, L., Schubert, B.A., Schultz, B., Tegner, C., Planke, S., Svensen, H.H., 2019a. Mercury anomalies across the Palaeocene–Eocene Thermal Maximum. Climate of the Past 15, 217-236.\nJones, S.M., Hoggett, M., Greene, S.E., Dunkley Jones, T., 2019b. Large Igneous Province thermogenic greenhouse gas flux could have initiated Paleocene-Eocene Thermal Maximum climate change. Nature Communications 10, 5547.\nLarsen, R.B., Tegner, C., 2006. Pressure conditions for the solidification of the Skaergaard intrusion: Eruption of East Greenland flood basalts in less than 300,000 years. Lithos 92, 181-197.\nPapadomanolaki, N.M., Sluijs, A., Slomp, C.P., 2022. Eutrophication and Deoxygenation Forcing of Marginal Marine Organic Carbon Burial During the PETM. Paleoceanography and Paleoclimatology 37, e2021PA004232.\nPlanke, S., Berndt, C., Alvarez Zarikian, C.A., Expedition 396 scientists, 2022. Expedition 396 Preliminary Report: Mid-Norwegian Continental Margin Magmatism. International Ocean Discovery Program.\nStorey, M., Duncan, R., Swisher III, C., 2007a. Paleocene-Eocene Thermal Maximum and the opening of the Northeast Atlantic. Science 316, 587-589.\nStorey, M., Duncan, R., Tegner, C., 2007b. Timing and duration of volcanism in the North Atlantic Igneous Province: Implications for geodynamics and links to the Iceland hotspot. Chemical Geology 241, 264-281.\nWesterhold, T., Röhl, U., Frederichs, T., Agnini, C., Raffi, I., Zachos, J.C., Wilkens, R.H., 2017. Astronomical calibration of the Ypresian timescale: implications for seafloor spreading rates and the chaotic behavior of the solar system? Climate of the Past 13, 1129-1152.\nWilkinson, C., Ganerød, M., Hendriks, B., Eide, E., 2017. Compilation and appraisal of geochronological data from the North Atlantic Igneous Province (NAIP)", "reviewer_in": "Péron-Pinvidic, G., Hopper, J.R., Stoker, M.S., Gaina, C., Doornenbal, J.C., Funck, T., Árting, U.E. (Eds.), The NE Atlantic Region: A Reappraisal of Crustal Structure, Tectonostratigraphy and Magmatic Evolution. Geological Society, London, Special Publications.\nZeebe, R.E., Lourens, L.J., 2019. Solar System chaos and the Paleocene–Eocene boundary age constrained by geology and astronomy. Science 365, 926-929.", "reviewer_2": "Jiang et al presents new records of the PETM and the precursor event named the POE from the eastern Tethys. There is a lot of new data presented in this manuscript - multiple proxies of carbon cycling and volcanism, ocean redox, and nutrient conditions (as the authors summarize in the abstract) as well as ocean temperature (clumped isotope measurements). In addition, the authors present results from simulations of the carbon injection across the PETM and POE using the cGENIE Earth system model.\nWhile the data are interesting, the new records raise many more questions than are addressed with this manuscript. My overall impression is that the manuscript fails to integrate the numerous data types and modeling in an effective or original manner. The content of the manuscript does not convincingly address the topic suggested by the title, which in any case is also not articulated in the text of the abstract that follows.\nThe title presents a very interesting question: is the climate state, presumably on millennial timescales, more or less stable than the ‘background’ Paleocene and Eocene climate? I say presumably because this is a key aspect of the question that is not elaborated - how do you define ‘climate instability?’ Over what timescale? I find only a few sentences in the main text that address the topic of the manuscript title (Lines 244-248) and argue that the use of a ‘sequential regime shift detection algorithm’ identifies ‘multiple regime shifts in the mean and variance of the time-series.’ None of this is sufficiently explained - neither the method nor the interpretation in either the main text or supplement. Out of all the records, why is this method only applied to the d13C and d18O records? Why use d18O at all if, as the authors caution, it may be ‘at least partially influenced by groundwater leaching and/or increased freshwater input’ (Line 123-124). How might this analysis be impacted by big changes in preservation across the section - in other words, is it valid to interpret noisier records in the body of the PETM as evidence for real climatic instability? Later in the Conclusions (Line 283) the authors suggest that ‘carbon cycle feedbacks…maintained the higher pCO2 during the body of the PETM, highlighting the climate instability.’ First, it isn’t clear that this study demonstrated the occurrence of carbon cycle feedbacks, and I’m not sure what evidence was meant to support this claim. Second, how is high CO2 alone evidence for climate instability?\nHowever, upon reading the abstract, the purpose of the manuscript instead seems to be to constrain the carbon source across the PETM and POE. The argument appears to be that volcanism was responsible for both events and organic carbon burial was an important recovery mechanism for the PETM. Neither of these arguments are novel, though much less is constrained about the POE compared to the PETM. However, the description of the findings is rather misleading because this study does not present any boron isotope data, which also has implications for the modeling conducted, as I describe below.\nThe authors use the double-inversion methodology of Gutjahr et al., 2017, in the cGENIE model. The main text inadequately describes what the authors have done here - they have taken d11B records from different sites and publications compared to their new d13C data (from Site 401 in the North Atlantic for the PETM and from a Maryland site for the POE). This has many unexplored consequences. For instance, have the authors considered the uncertainty in correlating these records? They’ve had to do this point-by-point to use records from two locations simultaneously as inversion targets. This correlation has direct impacts on the d13C of the source recovered in the model. Next, have they looked to see how well inversions of d13C from their new sites reproduce d13C curves in other locations, like equivalent to Site 401 from the Gutjahr paper? Each site is clearly impacted differently by dissolution and mixing and by combining the records used in the inversion the authors have inherited all these sources of uncertainty without doing anything to explore the implications for their results.\nBased on the modeling, the authors write, ‘shorter duration of the POE and less severe pH responses led us to conclude that the NAIP volcanic CO2 emissions may have led to the occurrence of the POE.’ Why does short duration and smaller pH response suggest volcanism? Instead, the diagnosed d13C source values from the double inversion are lower than for the PETM, and the authors appear to recognize this, writing that the results support ‘different triggering mechanisms for these two events (Line 234-235).’ I don’t understand what the authors mean by suggesting that both the POE and PETM are driven by NAIP volcanism but triggered by different mechanisms. Figure 5 summarizing the model results doesn’t provide any description in the caption for panels g through l for the POE inversion. There is no way that the records shown on the top in panel (g) were used to drive the entirety of the POE inversion given that no carbon emissions are diagnosed after 7k years.\nOverall, there is too much data and not enough rigorous interpretation in this manuscript. I’m left with many questions and without any clear takeaways. The title, abstract, and text of the manuscript are not well matched. I envision a study that could effectively use some (perhaps not all) of the presented data to drive model simulations of the PETM and POE and compare and contrast the results in terms of the diagnosed emissions, d13C, impact on CO2, and differential recovery mechanisms. Or, I can envision a study specifically focused on climate instability (as in the title), though this seems to require better constraints on temperature specifically."}, "author_response": "We appreciate reviewer #1’s positive feedback and thank you for your in-depth comments and suggestions, which have greatly improved the quality of the manuscript.\nImportant Comments:\nMercury signal.\nThe main issue I have with the paper is the interpretation of the mercury signal. All of the samples presented here have Hg <20 ppb and similarly low TOC contents, which means that analytical errors are going to be enhanced. On the machines I’ve used, Hg contents below 10 ppb have errors comparable to the concentrations, and TOC <0.2 wt% has been shown as the broad cut off for using Hg/TOC as a proxy (e.g. Grasby et al., 2019). Therefore, a lot more care needs to be taken before we can arrive at “originating from the NAIP” (lines 40.42). Hg/TOC values are almost all <30 ppb/wt%, which is well below the average for shales. Even the enriched sections during the POE and PETM are below this mean value from current data sets (Grasby et al., 2019). The authors need to address other possible avenues for disturbances to the Hg cycle. The authors mention and rule out wildfires as a possible source, but the main one that jumps out to me from their dataset is a possible change in organic matter source.\nThe TOC contents for our study section range from 0.05% to 3.8%, and therefore, we follow the reviewer’s suggestion to exert more caution when interpreting Hg/TOC as a proxy for volcanism. We add new discussions on the various processes that may have led to changes in the Hg/TOC ratios across the POE and the PETM:\n“However, because of the overall low Hg concentrations at the KZGS section, establishing a link between NAIP and the Hg peaks is not straightforward. The low Hg from the KZGS site is likely due to dilution by carbonate and detrital inputs, being more distant from the Hg source, and/or Hg transport via oceanic waters rather than global atmospheric transport12 (Fig. 3). Increased Hg concentrations across the POE and the PETM compared to background values suggest that multiple possible sources and processes may have been at play in addition to the NAIP volcanism. For example, changes in Hg concentrations in the sedimentary records can be caused by changes in river runoff, weathering, transport of terrestrial materials, primary productivity, source of organic matter, and post-depositional processes (e.g., diagenesis and dissolution)12, which could become more important at the KZGS section because of its restricted carbonate platform setting45. Ocean anoxia and changes in organic matter preservation and transport cannot fully account for the excess Hg as shown by the steeper Hg gradient to TOC within the PETM and POE interval at our site (Fig. 3).\nemissions (Fig. 3). C/N ratio exhibits no change across the POE, suggesting the increase in Hg and Hg/TOC ratio is unrelated to changes in source of organic matter. On the other hand, C/N ratio shows large increase across the PETM, which indicates that changes in source of organic matter may have contributed to the increased Hg concentrations. These potential processes do not preclude volcanic involvement, however, especially via more complex pathways than simple atmospheric loading and deposition47. Despite these potential complex sources of Hg, we cannot completely exclude direct and indirect involvement of the NAIP in driving the Hg changes in the study section48. For example, the NAIP was active as early as 62 Ma, and its peak activity may have encompassed the POE and the PETM 49-52. Enrichment in radiogenic Os isotopes suggests that the sediments from the Arctic and Peri-Tethys receive radiogenic materials from volcanic source53. Furthermore, hydrothermal vent complexes in the circum-Arctic and the Tethyan region suggest large emissions of Hg51,54,55, further supporting that the NAIP volcanic activities that encompassed the POE and the PETM can at least partially explain the observed Hg records. Geological evidence also supports increased activity of the NAIP, which encompassed both the POE and the PETM49-51 (PETM age: 56.01 Ma56).”\nRevised Figure S6. a) Cross plot of Hg concentration (ppb) with TOC (%) (b) PCA results that show the relative associations of sedimentary Hg with C/N ratios, δ13C org and other geochemical proxies.\nThere are several indicators of enhanced terrestrial runoff in this dataset (e.g. line 124), including a negative δ18O excursion and elevated C/N ratios (what does the BIT index of the GDGTs look like?). Mercury uptake into terrestrial plant biomass is inherently different to marine uptake, leading to different initial Hg/TOC ratios in each organic matter pool. Given the low Hg and TOC contents, how can the authors be sure that the signal they are observing is indeed due to changes in NAIP activity? There seem to be two distinct Hg/TOC populations in Figure S7, do the population with 10-16 ppb Hg and 0.5-1.5 wt% TOC also have high C/N ratios?\nMoreover, the authors interpret the POE as being the main warming driven by intrusive activity. Volcanic and thermogenic emissions are near impossible to differentiate in the far field, but the limited near-field dataset that we do have suggests that the submarine emissions from contact metamorphism around intrusions leads to much of the Hg being scrubbed into the local water bodies, rather than released to the atmosphere (e.g. Jones et al., 2019a). Why then would the POE have the larger relative Hg/TOC excursion at Kuzigongsu? I’m not saying that the signal here isn’t due to NAIP activity, but more nuanced discussion is needed in the interpretation of this dataset.\nThank you for these insightful suggestions. The GDGTs in this section are weak and sometimes completely absent, brGDGTs especially are absent in most samples, hence we focus on the two main isoGDGTs only: isoGDGT-0 and crenarchaeol. Where brGDGTs are present, BIT is low with an average value of 0.05 ± 0.18 across the section. Overall the brGDGTs are too weak in abundance to make any robust interpretations about changes in terrestrial run-off and hence we do not discuss this data in our manuscript. We added new discussions on why the POE have larger relative Hg/TOC excursion at KZGS (see above). We also revised Fig. 3 to show the relationship between Hg, TOC and C/N ratio (an indicator or sources of organic matter), which show that most data associated with population with 10-16 ppb Hg and 0.5-1.5 wt% TOC also have high C/N ratios (C/N >10), in particular for those with wt% TOC values greater than 1. PCA results also support C/N ratio and δ13C , org indicators of sources of organic matter and 13C-depleted CO emissions, exert primary 2 controls on the Hg concentrations.\nNAIP Emplacement.\nNorth Atlantic Igneous Province (NAIP) activity has been oversimplified in this manuscript. The province was active as early as 62 Ma, and while the main acme of activity did encompass the PETM, it is not (yet) constrained to as fine a window as 56-55.6 Ma (lines 162-163). This interval is based on the start of the crystallisation of the Skaergaard intrusion as it was buried by 5-6 km of lava (Larsen & Tegner, 2006) in East Greenland. The upper constraint is the eruption of the Gronau tuff, which is linked to Ash -17 in Denmark (Storey et al., 2007a). Once corrected to the most recent Fish Canyon Tuff age estimate, this layer is dated to 55.48 ±0.12 Ma, which gives us a ~0.5 Myr window where much of the East Greenland basalts were emplaced. However, this is just one part of a much bigger province, and the Skaergaard magma chamber was intruded into older basalts. Other parts of the province, particularly places like the Faroe Islands and the Vøring plateau, are not nearly so well constrained, and it is unclear at the moment whether this spike in activity in East Greenland is indicative of the province as a whole, or more localised activity.\nWe now follow the reviewer’s suggestion to add more discussions and related references on the NAIP activity. The revised text now reads as the following:\n“Despite these potential complex sources of Hg, we cannot completely exclude direct and indirect involvement of the NAIP in driving the Hg changes in the study section48. For example, the NAIP was active as early as 62 Ma, and its peak activity may have encompassed the POE and the PETM 49-52. Enrichment in radiogenic Os isotopes suggests that the sediments from the Arctic and Peri-Tethys receive radiogenic materials from volcanic source53. Furthermore, hydrothermal vent complexes in the circum-Arctic and the Tethyan region suggest large emissions of Hg51,54,55, further supporting that the NAIP volcanic activities that encompassed the POE and the PETM can at least partially explain the observed Hg records. Geological evidence also supports increased activity of the NAIP, which encompassed both the POE and the PETM49-51 (PETM age: 56.01 Ma56).”\nFurthermore, I strongly disagree with argument in lines 237-238. The study of Tian and Buck (2022) does not mention the NAIP once, and focuses on how the Deccan Traps and the Colombia river basalts were emplaced.\nThank you for pointing this out, the discussion related to Tian and Buck (2022) is now removed.\nLIP emplacements and the subsequent climatic impacts are different between individual LIPs (see Jones et al., 2016), and other emplacement models exist for other LIPs (e.g. Burgess et al., 2017 suggesting extrusive-intrusive-extrusive for the Siberian Traps). We know from field evidence that there is considerable NAIP volcanism prior to 56.5 Ma and after 55.6 Ma. The intrusive activity is less well temporally constrained, but the only two hydrothermal vent complexes that have been sampled and analysed contain the PETM onset within the vent infill (Frieling et al., 2016; Planke et al., 2022; Berndt et al., in review). There may have been a pulse in intrusive activity at the POE, but almost certainly also at the PETM. Many studies argue, quite convincingly, that magmatic carbon fluxes are insufficient to induce global warming on their own (e.g. Jones et al., 2019b). This section needs to be rewritten to be more nuanced in terms of what constitutes NAIP activity, how (if at all) it is possible to distinguish between volcanic and thermogenic emissions at such distance, and whether such definitive conclusions can be derived from a distal locality.\nThank you for these great suggestions. This section is now rewritten to include more discussions on what constitutes NAIP activity and how we estimate the likely contributions of volcanic versus thermogenic sources. The revised text now reads as the following:\n“Despite these potential complex sources of Hg, we cannot completely exclude direct and indirect involvement of the NAIP in driving the Hg changes in the study section48. For example, the NAIP was active as early as 62 Ma, and its peak activity may have encompassed the POE and the PETM 49-52. Enrichment in radiogenic Os isotopes suggests that the sediments from the Arctic and Peri-Tethys receive radiogenic materials from volcanic source53. Furthermore, hydrothermal vent complexes in the circum-Arctic and the Tethyan region suggest large emissions of Hg51,54,55, further supporting that the NAIP volcanic activities that encompassed the POE and the PETM can at least partially explain the observed Hg records. Geological evidence also supports increased activity of the NAIP, which encompassed both the POE and the PETM49-51 (PETM age: 56.01 Ma56).”\nMinor comments:\nTitle: The title can be improved, as the PETM is by definition a period of climate instability. The strength of this paper is the well preserved POE and its relation to the PETM, so I suggest rewording to draw attention to what is new and exciting here.\nThe relevant text is now revised to read: “Although the entire section is Hg-starved, peak Hg concentrations coincide with both the POE and the PETM, supporting direct and indirect involvement of the North Atlantic Igneous Province activity.”\nLine 41: The NAIP was not one eruption. Reword to “originating from North Atlantic Igneous Province activity” or similar.\nChanges are made following the reviewer’s suggestion: “… supporting direct and indirect involvement of the North Atlantic Igneous Province activity …”\nLine 72 (& Line 27 in SI): The coordinates for the Kuzigongsu section point to a warehouse in Wuqia county on Google Earth. Please check the coordinates.\nThe coordinates for the Kuzigongsu section are revised: “39°45'10\" N, 75°17'29\" E”.\nLines 87-88: “The onset of the PETM is defined as the most rapid decline in δ13C values in all studied substrates at 19.9 m” reads like this is the biggest excursion in this section, where I suspect you mean it is the largest of all PETM sections. I suggest you reword to “the onset of the PETM at 19.9 m displays a rapid decline in δ13C values, among the largest CIEs in studied PETM substrates (McInerney & Wing, 2011)”.\nChanges are made following the reviewer’s suggestion: “…the onset of the PETM at 19.9 m displays a rapid decline in the δ13C values, among the largest CIEs in studied PETM substrates1”\nLines 94-95: Please explain why these are deemed “more likely” and “less likely”? Is it to do with the way the age models are set up, or is it that age model 2 predicting a 39 kyr PETM duration goes against a mountain of evidence from other localities?\nFurther discussion is provided to explain why age model 2 is less likely. The revised text now reads: “Age model 2 is less likely because the predicted 39 kyr PETM duration goes against several lines of evidence from other localities suggesting the PETM duration ranges from 170 to 220 kyr27,28”\nLines 116-117: Pleas add oxford commas.\nOxford commas are added after “organic biomarker” and “nutrient”: “These proxies include C/N ratios as indications of organic matter source, wt.% CaCO as a proxy for 3 ocean acidification and detrital dilution, trace element geochemistry, organic biomarker, and mercury content as redox, nutrient, and volcanism indicators.”\nThis sentence is reworded following reviewer’s suggestion. The revised text now reads: “initiated at ca. 62 Ma and ended at ca. 54 Ma41 with a volume of ~9 million cubic km42”\nLines 185-186: The current best estimates for the onset of the PETM are either 55.93 Ma (Westerhold et al., 2017) or 56.01 Ma (Zeebe and Lourens, 2019). The Charles el al. (2011) paper provides an important ash marker horizon of 55.785 Ma, but this is during the recovery of the CIE. The shape of the CIE in Svalbard also makes it slightly difficult to pin down where in the PETM this marker horizon sits. Citations are needed for this sentence (e.g. Storey et al. 2007a; 2007b; Wilkinson et al., 2017).\nWe changed the age of the onset of the PETM to 56.01 Ma following Zeebe and Lourens (2019) and revised the text accordingly. The new text reads: “Geological evidence also supports increased activity of the NAIP, which encompassed both the POE and the PETM49-51 (PETM age: 56.01 Ma56).” These additional citations are now included in this sentence.\nLine 226: How does these model results compare to the findings of Gutjahr et al.?\nThe PETM model results are summarized in the new Table S5, and it is broadly similar to the results of Gutjahr et al. (2017). Because the PETM has been extensively modeled, we shifted our focus to the POE, and provided detailed sensitivity analyses on the carbon source, emission rate and total amount of the emitted carbon (revised Fig. 4 and Table S5).\nLines 235 & 284: Activity is only classified volcanic at the surface. Intrusive activity is magmatism.\nThank you for pointing this out. We now changed “intrusive volcanism” as “intrusive activity”.\nLine 266: Also cite Papadomanolaki et al. (2022) here.\nThis reference is now cited.\nLine 284: Reword to “as a result of intrusive and/or eruptive activity from a major large igneous province.”\nThis sentence is now reworded following reviewer’s suggestions: “highlighting the climate instability as a result of intrusive and/or eruptive activity from a major large igneous province”\nAuthor contributions: It is not clear what the last four authors contributed to the manuscript.\nFigure 5: The colour palate used for this figure is not very friendly towards the colour blind, suggest changing to colours with greater contrast (see Crameri et al., 2020, for useful tips).\nThank you for this suggestion. The colour palate is now changed to be more friendly towards the colour blind following Crameri et al. (2020). The revised Fig. 5 (now Fig. 4) is shown below:\nSimilar to the reply to comment above: The GDGTs in this section are weak and sometimes completely absent, brGDGTs especially are absent in most samples, hence we focus on the two main isoGDGTs only: isoGDGT-0 and crenarchaeol. Where brGDGTs are present, BIT is low with an average value of 0.05 ± 0.18 across the section. Overall the brGDGTs are too weak in abundance to make any robust interpretations about changes in terrestrial run-off and hence we do not discuss this data in our manuscript.\nSupplementary lines 248-251: A value of 2600 kg m3 seems a very high estimate for a density of shale. Are your samples that metamorphosed? We measured the density of the PETM clay in Denmark (paper in prep.) and got a value of 1400 kg m3. Admittedly this is an unlithified section, but if the value used here is correct then the section may have undergone significant compaction post-deposition. Some text addressing this, or better yet a measurement of the actual shale density, would be useful here.\nWe measured the density of the rock samples in the Kuzigongsu section, and got an average value of 2270 kg m-3 for mudstone samples (n = 5), which is higher than those of the PETM clay in Denmark as the reviewer mentioned. The calculation of organic carbon burial using the measured density is redone using these new density values.\nSupplementary line 254: Table S4 leaves out several recently studied sections that have much thicker complete sections than this, including Fur, Svalbard, the North Sea, and the Norwegian continental margin (e.g. Jin et al., 2022; Jones et al., 2019a; Planke et al., 2022). Perhaps tone down the language a little to say it is one of the thicker PETM sequences, or that it is one of the thickest carbonate-rich sections?\nWe now tone down the language to say it is one of the thickest carbonate-rich sections. The revised sentence now reads: “Given that the POE and the PETM span 1.0 m and 10.6 m, respectively, at the studied section, the KZGS section is considered as one of the thickest carbonate-rich sections”. These recent references are also added in Table S4.\nSupplementary line 287: Correct typo to “performed”\nThis typo is now corrected to “performed”."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The revised paper by Jiang and co-authors is a significant improvement on the first submission and is an excellent and exciting paper overall. The only area where I believe significant further investigation is required is the modelling results and the implications for a plausible source. After addressing this one one moderate comment and a few minor revisions, it is worth publishing in Nature Communications. With best regards, Morgan Jones\nMain comment:\nLines 266 to 272: I am a little worried about the validity of the modelling results that ascribe a purely NAIP volcanic carbon source for the POE, as 1 to 10 Pg C/yr is an enormous flux, and up to 12,300 Pg C is an enormous volume. To put these numbers into perspective, Stephen Jones et al. (2019)’s paper estimates that the total mantle-derived carbon reservoir for the entire NAIP is 14,100 Pg C, while our recent estimate puts it between 21,000 and 35,000 Pg C (M. Jones et al., 2023). These are estimates for the entire 62 to 54 Ma interval, and while much of this volume is believed to have erupted between 56 and 55.5 Ma, it would still need a phenomenal amount of magma/lava to degas in an extremely short amount of time to solely account for the POE or the PETM. Taking the approximation of 3.5 Tg C per km3 of magma as a potential degassing source (see Jones et al., 2016), a yearly flux of 1 to 10 Pg C/yr would require an eruptive / shallow intrusive flux of 286 to 2860 km3/yr of magma, which is absolutely colossal. The geochronology of the NAIP continental flood basalt sequences is not very well constrained, but the field evidence in localities such as East Greenland and the Faroe Islands shows individual lavas with evidence of hiatuses between that seem to contradict the hypothesis that large volumes of these lavas were erupted in a very short (e.g. 1–10 kyr) time periods. Therefore, the model results appear to generate unrealistic values for magmatic gas fluxes, and this discrepancy between the model findings and the NAIP magmatic carbon reservoir potential needs to be discussed in greater detail in the manuscript.\nI don’t think that the model needs to be rerun, but some discussion about the potential sources and differences between the model and field observations is warranted. Potential avenues include how the model setup may vary from natural conditions, and what NAIP sources may be active. Does the model cope well with an enhanced 12C sink (e.g. organic matter burial) that could be masking the degassing signature in atmospheric CO values, perhaps? The cGENIE model is not my expertise, but 2\nI think that the paper would be significantly enhanced with some discussion around what the limitations of the model are, what potential reasons can there be between the calculated fluxes from modelling and field observations, and how future studies could constrain these uncertainties.\nMinor comments:\nThe focus on the NAIP as a purely volcanic source is a little misleading, as numerous studies (e.g. Berndt et al., 2023) highlight the importance of sub-volcanic intrusions as a significant carbon flux across this key interval. The explosive ejection of thermogenic gases through these hydrothermal vent complexes would be a mix of magmatic and sedimentary sources. Stephen Jones et al. (2019) estimate that the partition would be would be 90% contact metamorphic and 10% mantle signature, and that peak emission fluxes of 0.2–0.5 Pg C/yr from thermogenic NAIP sources could have initiated the PETM. I personally think the 90:10 ratio is a little high, as they assume a thick metamorphic aureole in their calculations and do not account for magmatic degassing from deeper in the plumbing system. Therefore, the intrusive part of the NAIP can potentially be a substantial source of magmatic carbon, so I would err on the side of caution and recommend referring to it as “NAIP activity” or “emplacement of the NAIP” rather than explicitly “volcanism” throughout the text.\nLine 124: Add ‘a’ to “and mercury content as a possible indicator of volcanism”\nLine 128: Add ‘evidence of’ to “…represent evidence of an abrupt and significant warming…”\nLines 180–181: This explanation of Mn* would look better as an equation.\nLine 192: Elements at the start of sentences should be written out in full (i.e. change Hg to Mercury here)\nLine 203: For reasons stated above, I think it would be better to change “volcanism” to “activity” here.\nLine 216: As above. Either change to “The C/N ratio” or “Carbon/nitrogen ratio”\nLines 229–230: The locations of hydrothermal vent complexes here is wrong. The hydrothermal vent complexes are exclusively found on the continental margins around the NAIP (modern day Northeast Atlantic Ocean), caused by igneous intrusions into sedimentary sequences. See the recent paper by Berndt et al. (2023).\nLine 234: I would probably recommend not stating the estimated PETM onset age of 56.01 Ma (Zeebe and Lourens, 2019) here, as the U-Pb age of 55.785 ± 0.034 Ma from a bentonite layer that is within the PETM CIE body in Svalbard (Charles et al., 2011) does not match with the 127 kyr PETM duration presented in this paper. You can perhaps mention this discrepancy and the alternate 55.93 Ma PETM onset age (Westerhold et al., 2017) if you want to dig into what the cyclostratigraphy presented in this paper suggests for the validity of the various age models.\nLine 554 (Figure 1): Does the “CIE” marked on panel (d) refer to the PETM? Please clarify\nLines 635–636 and 756–758: References numbered 12 and 51 are a repeat, the latter is the correct version.\nReferences mentioned:\nBerndt, C., Planke, S., Alvarez Zarikian, C. A., Frieling, J., Jones, M. T., Millett, J. M., Brinkhuis, H., Bünz, S., Svensen, H. H., Longman, J., Scherer, R. P., Karstens, J., Manton, B., Nelissen, M., Reed, B., Faleide, J. I., Huismans, R. S., Agarwal, A., Andrews, G. D. M., Betlem, P., Bhattacharya, J., Chatterjee, S., Christopoulou, M., Clementi, V. J., Ferré, E. C., Filina, I. Y., Guo, P., Harper, D. T., Lambart, S., Mohn, G., Nakaoka R., Tegner, C., Varela, N., Wang, M., Xu, W., and Yager, S. L.:\nShallow-water hydrothermal venting linked to the Palaeocene–Eocene Thermal Maximum, Nat. Geosci.,\nhttps://doi.org/10.1038/s41561-023-01246-8, 2023.\nJones, S. M., Hoggett, M., Greene, S. E., and Dunkley Jones, T.: Large Igneous Province thermogenic greenhouse gas flux could have initiated Paleocene-Eocene Thermal Maximum climate change, Nat. Commun., 10, 5547,\nhttps://doi.org/10.1038/s41467-019-12957-1, 2019.\n(Currently reference No. 13 in this manuscript)\nJones, M. T., Jerram, D. A., Svensen, H. H., and Grove, C.: The effects of large igneous provinces on the global carbon and sulphur cycles, Palaeogeogr. Palaeocl., 441, 4–21, 2016.\n(Ref. 64)\nJones, M.T., Stokke, E.W., Rooney, A.D., Frieling, J., Pogge von Strandmann, P.A.E., Wilson, D.J., Svensen, H.H., Planke, S., Adatte, T., Thibault, N., Vickers, M.L., Mather, T.A., Tegner, C., Zuchuat, V., Schultz, B.P.: Tracing North Atlantic volcanism and seaway connectivity across the Paleocene–Eocene Thermal Maximum (PETM). Clim. Past 19, 1623-1652, 2023.\n(Ref. 55, updated from the preprint version now it is published)\nWesterhold, T., Röhl, U., Frederichs, T., Agnini, C., Raffi, I., Zachos, J. C., and Wilkens, R. H.: Astronomical calibration of the Ypresian timescale: implications for seafloor spreading rates and the chaotic behavior of the solar system?, Clim. Past, 13, 1129–1152", "reviewer_2": "I reviewed a previous version of this manuscript, and I agree with the authors’ decision to focus on the diagnosis of the likely carbon source for the POE. The additional experiments conducted as sensitivity analysis are welcome, but unfortunately I have many more unanswered questions raised by this version of the manuscript.\nOverall, I don’t think that a clear and succinct summary of the double inversion experiments has been presented. I do not think that the authors have clearly demonstrated that volcanic emissions can explain the POE without generating a d13C excursion that should have been observable in deep sea records. It is very challenging to interpret all experimental results by trying to combine the description of each experiment provided in a supplemental table with multiple figures with lines labeled ‘Experiment 1, 2, etc.’ and captions that are not re-written for each figure. It was difficult to remember what I was meant to be evaluating from each plot contrasting Figure 4 with the figures in SI. Key sensitivity analyses seem to be 1) the significance of varying the age model 2) importance of different assumptions about the pH constraint, and 3) importance of the alignment of d13C and pH targets. I guess the goal is to summarize the importance of 1) and 2) in Figure 4 and leave 3) to the SI. Overall, it would help for each subplot in each figure to identify the plotted lines in a way that clearly explains what is being varied (not ‘Experiment ID =’). Importantly, comparison of the significance of different age model assumptions is hindered a bit by the use of different axis limits for the right and left-side subplots showing carbon emissions. The reported diagnosed d13C of the source (-4 to -13 per mil) (Line 267-268) is not what it looks like from Figure 4 (which is also difficult to see because the axis extends to -60 per mil).\nA significant concern is the misrepresentation of the benthic d13C excursion generated in the experiments. In the main text, it is reported that the POE double inversions yield ‘no obvious CIE’ in the modeled deep ocean (Line 275). I don’t agree with that interpretation based on what is shown in Figure S11. The benthic CIE in Fig S11a is roughly 1.6 per mil and roughly 1 per mil in Fig S11b. In panels c-f benthic d13C is still decreasing because the duration of the experiment is so short that the deep ocean has either barely seen or has not seen the full magnitude of the excursion in the surface. It is very hard to argue that one would not expect to see a 1 to 1.6 per mil excursion in the deep ocean as predicted by these experiments - how do you explain the abundance of ‘hyperthermals’ identified of approximately this magnitude across the Early Eocene? It is a really interesting question whether inversions of the POE CIE and pH are compatible with the supposed absence of this event from the deep sea record. Viewing, for instance, the high resolution benthic foram d13C record from ODP Site 1262 (Littler et al., 2014) there does not appear to be evidence for a -1 per mil excursion ~100 kyr prior to the PETM. Does this suggest something about the amount of time between the two events and hence the likelihood that the POE could be erased by burndown? It is not the duration alone of an event that limits its appearance in the deep sea - it is the total size of the event that matters. (If only duration mattered, that would erroneously suggest that the impact of modern carbon emissions wouldn’t be recorded in the geologic record). Some of the inversions presented here show a really size-able event - CO2 more than doubling and masses of carbon predicted that are similar to some previously published estimates for the PETM itself - how are these ‘relatively small’ emissions estimates?\nThere are additional aspects of the analysis that are confusing to me (itemized by line below). The figures that present the age model development and the model results are difficult to read and have very small text. It is strange to leave the cyclostratigraphy to SI given the centrality of age model constraints to the interpretation of the POE via Earth system modeling. Text is duplicated in both Methods and SI but sometimes with slight differences. The SI should not repeat the Methods but provide additional explanation/discussion (for instance, an analysis of the sensitivity experiments).\nMoreover, the Methods (particularly regarding the cyclostratigraphy and modeling) are incomplete. Cyclostratigraphy uses magnetic susceptibility, but MS measurements are never described. What is the resolution of MS measurements? Exactly how is the record interpolated, detrended, smoothed or filtered? ‘Dynamical filter to isolate precession cycles’ needs elaboration - where is the evidence that there is significant precession power in the timeseries prior to filtering the record at this frequency and hence adding power? How is a ‘reasonable range of sedimentation rate’ set? What are the overall constraints? Multiple software packages are described (MATLAB scripts, AnalySeries) but then the final sentence states that all analyses are conducted with Astrochron unless otherwise noted. My overall impression is that this description of cyclostratigraphic methods doesn’t give the reader the ability to reproduce the results, which should be the goal.\nRegarding the description of Earth system modeling, the Methods first states that the isotopic inversion target is CO2 d13C (Line 503) and then later that the target is surface ocean DIC d13C (Line 511). What was the d13C record from the KZGS adjusted to in absolute values - CO2 d13C or DIC 13C? Line 496 is duplicated in Line 523 (details of the SPINUP) with configuration variously described as late Paleocene and early Eocene. In the SI (Line 315) there is mention of different pH scenarios for both the PETM and POE inversions but these appear to have been handled differently. I think the different pH scenarios for the PETM double inversion actually have adjusted pH targets in terms of absolute value whereas the different pH scenarios for the POE adjust only the magnitude of the pH change but this is really unclear. For the PETM double inversion, what do you mean by pH is ‘set to’? From Fig S8e, there is a difference in initial CO2 across the three different pH values (suggesting the absolute pH was offset initially in each inversion from the SPINUP). Why were these different approaches taken for the PETM and POE double inversions? (In Figure 4, the different pH options for the POE inversions don’t impact the initial value of CO2, but the data shown in grey versus red dots do have differences in the absolute values of pH in addition to the magnitude of the change in pH).\nAdditional comments below:\nLine 92 - provide the breakdown of CIE size across source for the POE as done above for the PETM?\nLine 95 - in the main text, could briefly summarize what the cycle assumption is for each of the two tuned age models?\nLine 144 - conceivably there could also be a significant decline in carbonate production (rather than reduced preservation alone)?\nLine 151 v Line 162 - duplicate description of elevated C29 hopane ratios during the PETM\nLine 174 - should be ‘perturbing’\nLine 172 - 177 - the wording of this sentence is unclear; why specifically related 2-MeHop to nitrogen cycle perturbation in this sentence when the previous sentence referred to a host of interpretations of which changes in N fixation was one option?\nLine 188 - ‘global extent’ is misleading; anoxia is documented in relatively restricted basins\nLine 195 - if there were a direct link between Hg source and 13C-depleted source, why would there be a lag of more than 10 kyr?\nLine 211-213 - is this meant to argue that wildfires in general cannot explain the Hg record or wildfires in the Arctic in particular?\nLine 305 - didn’t the double inversion provide constraints on the necessary magnitude of organic carbon burial during the POE? What is the rate and d13C of carbon removal in the simulated POE recovery?\nLine 309 - what does it mean for the system to ‘transition into the PETM’? If the argument is that the PETM was driven also by NAIP emissions, then does it matter whether or not organic carbon burial accelerated the recovery from the emissions that drove the POE? A more interesting question is to what extent CO2 has recovered after POE emissions but before the PETM and whether that requires additional drawdown fluxes.\nMethods - location of analyses is not mentioned for n-alkanes only\nFigure 2 - how is the choice made where to identify the PETM onset? In the supplement, the shading to indicate the PETM covers a different interval, with the onset aligned to the initial decline in d13C of CaCO3 rather than d13C of Corg. The overall description of the wt% CaCO3 record in the text doesn’t reflect the figure. The major decline in CaCO3 doesn’t appear to have anything to do with the timing of the POE or the PETM CIE, though an increase in CaCO3 is aligned with the end of the PETM.\nReference:\nLittler, K., Röhl, U., Westerhold, T. and Zachos, J.C., 2014. A high-resolution benthic stable-isotope record for the South Atlantic: Implications for orbital-scale changes in Late Paleocene–Early Eocene climate and carbon cycling. Earth and Planetary Science Letters, 401, pp.18-30.", "reviewer_author_response": "Reply: We wish to express our gratitude for reviewer’s constructive suggestions for improving the quality of the manuscript, and have made revisions according to reviewer’s comments. Line-by-line reply to reviewer’s comments is provided below.\nWe also further tighten up the text such that it exhibits a natural flow of speech and grammatical structure.\nMain comment:\nLines 266 to 272: I am a little worried about the validity of the modelling results that ascribe a purely NAIP volcanic carbon source for the POE, as 1 to 10 Pg C/yr is an enormous flux, and up to 12,300 Pg C is an enormous volume. To put these numbers into perspective, Stephen Jones et al. (2019)’s paper estimates that the total mantle-derived carbon reservoir for the entire NAIP is 14,100 Pg C, while our recent estimate puts it between 21,000 and 35,000 Pg C (M. Jones et al., 2023). These are estimates for the entire 62 to 54 Ma interval, and while much of this volume is believed to have erupted between 56 and 55.5 Ma, it would still need a phenomenal amount of magma/lava to degas in an extremely short amount of time to solely account for the POE or the PETM. Taking the approximation of 3.5 Tg C per km3 of magma as a potential degassing source (see Jones et al., 2016), a yearly flux of 1 to 10 Pg C/yr would require an eruptive / shallow intrusive flux of 286 to 2860 km3/yr of magma, which is absolutely colossal. The geochronology of the NAIP continental flood basalt sequences is not very well constrained, but the field evidence in localities such as East Greenland and the Faroe Islands shows individual lavas with evidence of hiatuses between that seem to contradict the hypothesis that large volumes of these lavas were erupted in a very short (e.g. 1–10 kyr) time periods. Therefore, the model results appear to generate unrealistic values for magmatic gas fluxes, and this discrepancy between the model findings and the NAIP magmatic carbon reservoir potential needs to be discussed in greater detail in the manuscript.\nReply: We now provide new discussions to compare the modeled carbon emission and field-based evidence of intrusive and extrusive activities. First, we provide a clear and succinct summary of the double-inversion experiments, with a focus on our preferred model scenario associated with a median ΔpH and age option 1, which provides an upper limit for the estimated carbon emission flux (lines 272-304). We also show model results of additional analyses turning on organic carbon burial feedback during the recovery phase of δ13C (lines 284-286; Methods and Table S5a). Greater details on the model findings and the NAIP magmatic carbon reservoir potential are now provided in lines 342-365, which justify the simulated carbon emission as an upper limit.\nI don’t think that the model needs to be rerun, but some discussion about the potential sources and differences between the model and field observations is warranted. Potential avenues include how the model setup may vary from natural conditions, and what NAIP sources may be active. Does the model cope well with an enhanced 12C sink (e.g. organic matter burial) that could be masking the degassing signature in atmospheric CO values, perhaps? The cGENIE model is not my expertise, but I think 2\nthat the paper would be significantly enhanced with some discussion around what the limitations of the model are, what potential reasons can there be between the calculated fluxes from modelling and field observations, and how future studies could constrain these uncertainties.\nReply: Regarding an enhanced 12C sink associated with organic carbon burial in cGENIE, we ran additional model simulations turning on organic carbon burial flux from ~8 kyr onward when δ13C starts to recover, which allows for the removal or excess organic carbon from the surface ocean (see Table S5a for model results with varying organic carbon burial scaling factor). This was done when the modelled mean global surface ocean δ13C value becomes lower than the observed value at the KZGS site, 13C-depleted carbon is removed from the system to force the δ13C to become more positive via negative emissions (e.g., Gutjahr et al., 2017). The new modeling result suggests ~ 800 to 1,900 Pg C cumulative removal of light carbon at average δ13C value of −30‰ and a rate of ~ −0.1 to −0.3 Pg C yr-1, which lead to a better match between modelled and observed recovery of both δ13C and pH. It is also DIC\nimportant to note that a caveat of cGENIE in interpreting our results is the lack of terrestrial biosphere and potential changes in orbital forcing, which could impact the climate responses and lead to uncertainties in carbon emission estimates. This caveat is acknowledged in lines 358-361. The revised text reads as the following:\n“Indeed, we find that applying a time-varying organic carbon burial calculated to nudge model δ13C closer to the data during recovery (~8 kyr onwards), we achieve DIC\nnot only a better match between the modelled and observed δ13C recovery but also DIC\nto the simulated pH recovery (Fig. 4a). The diagnosed organic carbon burial rate is approximately −0.1 to −0.3 Pg C yr-1 (Fig. 4c) with an average modelled δ13C value of −30‰, and leading to a cumulative C burial over the span of POE recovery of ~ 800 to 1,900 Pg C (Fig. 4d). …. It is also important to note that a caveat of cGENIE in interpreting our results is the lack of terrestrial biosphere and potential changes in orbital forcing, which could impact the climate responses and lead to uncertainties in carbon emission estimates.”\n“Although this study provides an upper constraint on the carbon emission flux during the POE, more precise δ11B-based global surface pH records, detailed history of the Iceland plume pulsing12, 77, sea surface temperature records from across different latitudes, and better-constrained geochronology of the NAIP activity are clearly needed to reduce the uncertainty of the estimated carbon emission fluxes.”\nMinor comments:\nThe focus on the NAIP as a purely volcanic source is a little misleading, as numerous studies (e.g. Berndt et al., 2023) highlight the importance of sub-volcanic intrusions as a significant carbon flux across this key interval. The explosive ejection of thermogenic gases through these hydrothermal vent complexes would be a mix of magmatic and sedimentary sources. Stephen Jones et al. (2019) estimate that the partition would be would be 90% contact metamorphic and 10% mantle signature, and that peak emission fluxes of 0.2–0.5 Pg C/yr from thermogenic NAIP sources could have initiated the PETM. I personally think the 90:10 ratio is a little high, as they assume a thick metamorphic aureole in their calculations and do not account for magmatic degassing from deeper in the plumbing system. Therefore, the intrusive part of the NAIP can potentially be a substantial source of magmatic carbon, so I would err on the side of caution and recommend referring to it as “NAIP activity” or “emplacement of the NAIP” rather than explicitly “volcanism” throughout the text.\nReply: Thank you for this great suggestion. We now reword “volcanism” as “NAIP activity” throughout the text.\nLine 124: Add ‘a’ to “and mercury content as a possible indicator of volcanism”\nReply: Done.\nLine 128: Add ‘evidence of’ to “…represent evidence of an abrupt and significant warming…”\nReply: Done.\nLines 180–181: This explanation of Mn* would look better as an equation.\nReply: Done.\nLine 192: Elements at the start of sentences should be written out in full (i.e. change Hg to Mercury here)\nReply: Done.\nLine 203: For reasons stated above, I think it would be better to change “volcanism” to “activity” here.\nReply: Done.\nLine 216: As above. Either change to “The C/N ratio” or “Carbon/nitrogen ratio”\nReply: Done.\nLines 229–230: The locations of hydrothermal vent complexes here is wrong. The hydrothermal vent complexes are exclusively found on the continental margins around the NAIP (modern day Northeast Atlantic Ocean), caused by igneous intrusions into sedimentary sequences. See the recent paper by Berndt et al. (2023).\nReply: We followed the reviewer’s suggestions to correct the locations of the hydrothermal vent complexes and cited relevant papers, including Berndt et al. (2023) and Frieling et al. (2016).\nLine 234: I would probably recommend not stating the estimated PETM onset age of 56.01 Ma (Zeebe and Lourens, 2019) here, as the U-Pb age of 55.785 ± 0.034 Ma from a bentonite layer that is within the PETM CIE body in Svalbard (Charles et al., 2011) does not match with the 127 kyr PETM duration presented in this paper. You can perhaps mention this discrepancy and the alternate 55.93 Ma PETM onset age (Westerhold et al., 2017) if you want to dig into what the cyclostratigraphy presented in this paper suggests for the validity of the various age models.\nReply: We followed the reviewer’s suggestion and removed the estimated PETM onset age of 56.01 Ma.\nLine 554 (Figure 1): Does the “CIE” marked on panel (d) refer to the PETM? Please clarify\nReply: The CIE marked on panel (d) does refer to the PETM. We now clarify this issue in the revised figure 1d.\nLines 635–636 and 756–758: References numbered 12 and 51 are a repeat, the latter is the correct version.\nReply: this is now corrected."}, "author_response": ""}, {"round": 3, "reviewer_comments": {"reviewer_1": "This is my third review of Jiang et al. for submission to Nature Communications. The manuscript is an improvement on the previous iteration and details an exciting new POE and PETM locality. I want to be supportive of this paper because of the high-resolution nature of the Kuzigongsu outcrop and the coverage of the POE. Moreover, I appreciate the efforts that the authors have put in to advance the paper and address my previous comments and suggestions. Sadly, I still find that the modelling part of this paper is its Achilles heel and not publishable in its current state. I am not a modeller by trade, so I cannot offer advice on where the issue lies, only that the output number of 230 km3 magma /year from the NAIP for 8000 years defies all our current understanding of LIP emplacements, mantle convection, and climate feedbacks. This magma flux estimate is probably around two orders of magnitude too high for the late Thanetian. I will detail a few of the reasons why here:\n1) Timing. If the POE is ~144 kyr before the PETM onset, then it is around 56.07 Ma (using Westerhold et al), 56.15 Ma using Zeebe and Lourens (2019). This predates much of the known NAIP lavas, even within uncertainty. In East Greenland, we have the precise age of the Skaergaard intrusion at 56.02 Ma (Wotzlaw et al., 2012), which is the anchor for the emplacement age of the flood basalts that buried the magma chamber as it crystallised (Larsen and Tegner, 2006). Skaergaard was intruded into the lower basalt sequence, but almost all the subsequent flood basalt activity postdates the intrusion. A similar picture emerges on the Norwegian margin. The lower series encountered on IODP 642 (see http://www-odp.tamu.edu/publications/104_SR/VOLUME/CHAPTERS/sr104_51.pdf) likely contains the PETM (IODP 396 scientists are working on this right now), but that means that the main volcanism is either syn- or post-PETM in age. The geochronology of the Faroes is less well constrained, but it is probably similar to these other localities. There are a few lava flows that might possibly be the right sort of age to coincide with the POE (see site U1566 for example; Planke et al., 2023), but these are volumetrically a very small component of the NAIP. It is certainly nowhere near “19% to 38% of the total NAIP volume” as stated in line 356. Even 0.1% of the total NAP volume in 8000 years would be an extreme scenario. Which brings us on to:\n2) Volume. 230 km3(magma)/year is an insanely high flux for a single eruption, let alone a flux from a LIP that continued for 8000 years. That is the equivalent of a Yellowstone sized-eruption’s worth of material emplaced in a decade, 800 times in a row. For comparison, the 1783-84 eruption Laki that is one of the largest known basaltic eruptions was 14 km3 volume. It is also very difficult to envision a mantle convection model that can emplace such a large volume in such a short space of time, go into almost a complete period of quiescence, then restart at similar flux levels for the PETM. Therefore, these estimated fluxes are well beyond the realms of possibility.\n3) Contemporaneous effects. Volcanoes don’t just emit carbon, and some co-emitted volatiles would have a significant climate impact. Emplacing 1.8 x 106 km3 magma in 8000 years would release considerable volumes of such volatiles. Take sulfur for instance, using the 2014-15 eruption of Holuhraun in Iceland as an analogue for some quick back of the envelope calculations to illustrate this point. The Holuhraun eruption was 1.6 km3 magma volume, with estimated emissions of 6.7 Tg SO (Carboni et al., 2019). This is roughly 3.35 Tg S, and 2.09 Tg S / km3 magma. Scaling this up to the proposed flux (230 2 km3/yr) gives a sulfur flux of 482 Tg S/yr. Comparing this to pre-industrial cycle would increase the total atmospheric S flux by more than an order of magnitude (sea spray = 13-36 Mt S/yr; volcanoes = 6.5 – 10.5 Mt S/yr; 1Mt = 1 Tg), and nearly a 50x increase on volcanic fluxes (see Jones et al., 2016 for discussion). Experiencing these conditions for 8000 years would have catastrophic effects on terrestrial and freshwater ecosystems in particular, through a sustained volcanic winter and increase in acid rain. To reuse the Laki example, that eruption caused crop failures in the Americas and Europe and mass famine in Iceland from one year of volcanic winter. There just isn’t the evidence in the geological record for this type of extinction-level disturbance during the POE, nor really during the PETM.\nThis leaves us with the conclusion that the model presented cannot be representative of the POE. This could be due to incorrect input parameters, issues with the simplified nature of the model itself, unidentified sources/sinks, or some combination thereof. This also calls into question the assertion that the PETM is largely driven by volcanic CO , if Gutjahr et al. (2017) used similar 2 model parameters as presented here. The main issue is that the only NAIP emissions that are likely to be capable of producing such large fluxes on sub-millennial timescales are thermogenic sources (cf. S. Jones et al., 2019). The only way that can be reconciled with the presented model is if the volatile flux from hydrothermal vents was magma-dominated over contact metamorphic-dominated and that many of the vents were active ~100 kyr before the first known vent activity. Even then, the fluxes required are colossal, which suggests we are missing other potential sources or processes that are not accounted for in the model set up. Unfortunately I don’t know what the best way forward is, as I do not have the expertise to suggest where the model is falling short. However, unless the model can produce a POE that is grounded in realistic NAIP fluxes, I don’t think it should be included in this paper.\nWith best regards, Morgan Jones\nMinor comments:\nLine 90 (and 119, ++): I don’t see the benefit of shortening Kuzigongsu to KZGS, it doesn’t significantly shorten the text and adds a barrier to readability. I suggest just using Kuzigongsu to describe the locality.\nLine 142: Is the ‘≥’ symbol correct? If the water depth is from around 30 m up to 50 m, then the symbol should be ‘≤’ instead.\nLines 158-162: This is a long sentence. I suggest you break in two to improve readability (I would break it into “…same depth. However, the thermal maturity…”\nLine 189: Suggest splitting this sentence too and fixing a slight grammar typo in doing so: “…as a result of biogeochemical changes. This is similar to observations of other major carbon cycle perturbations…”\nLine 244: Change “receive” to “received”.\nReferences mentioned:\nCarboni, E., Mather, T.A., Schmidt, A., Grainger, R.G., Pfeffer, M.A., Ialongo, I., Theys, N. (2019). Satellite-derived sulfur dioxide (SO2) emissions from the 2014–2015 Holuhraun eruption (Iceland). Atmospheric Chemistry and Physics 19 (7), 4851-4862.\nGutjahr, M., Ridgwell, A., Sexton, P. F., Anagnostou, E., Pearson, P. N., Pälike, H., Norris, R. D., Thomas, E., and Foster, G. L. (2017). Very large release of mostly volcanic carbon during the Palaeocene–Eocene Thermal Maximum, Nature, 548, 573–577.\nJones, S.M., Hoggett, M., Greene, S.E., and Dunkley Jones, T. (2019) Large Igneous Province thermogenic greenhouse gas flux could have initiated Paleocene-Eocene Thermal Maximum climate change, Nat. Commun., 10, 5547,\nhttps://doi.org/10.1038/s41467-019-12957-1.\nLarsen, R. B. and Tegner, C. (2006). Pressure conditions for the solidification of the Skaergaard intrusion: Eruption of East Greenland flood basalts in less than 300,000 years, Lithos, 92, 181–197.\nPlanke, S., Berndt, C., Alvarez Zarikian, C.A., Agarwal, A., Andrews, G.D.M., Betlem, P., Bhattacharya, J., Brinkhuis, H., Chatterjee, S., Christopoulou, M., Clementi, V.J., Ferré, E.C., Filina, I.Y., Frieling, J., Guo, P., Harper, D.T., Jones, M.T., Lambart, S., Longman, J., Millett, J.M., Mohn, G., Nakaoka, R., Scherer, R.P., Tegner, C., Varela, N., Wang, M., Xu, W., and Yager, S.L. (2023). Site U1566. In Planke, S., Berndt, C., Alvarez Zarikian, C.A., and the Expedition 396 Scientists, Mid-Norwegian Margin Magmatism and Paleoclimate Implications. Proceedings of the International Ocean Discovery Program", "reviewer_396": "College Station, TX (International Ocean Discovery Program). https://doi.org/10.14379/iodp.proc.396.104.2023.\nWotzlaw, J., Bindeman, I., Schaltegger, U., Brooks, C., and Naslund, H. (2012). High-resolution insights into episodes of crystallization, hydrothermal alteration and remelting in the Skaergaard intrusive complex, Earth Planet. Sc. Lett., 355, 199–212.", "reviewer_2": "I reviewed a previous version of this manuscript and I appreciate the authors responses to my comments, particularly the addition of more information about the age model construction. New high-resolution PETM records, including expression of the POE, are a welcome contribution, and additional age information about the duration of the POE and between the POE and PETM is especially significant. The problem is that I am ultimately not convinced that the way the double inversion modeling approach is applied and interpreted, and the central conclusions of the manuscript about carbon forcing that result, are valid for the POE. Central to the concept of the double inversion, as applied by Gutjahr et al., 2017, across the PETM, is the understanding that the PETM d13C excursion represents at whole ocean perturbation driven by carbon addition from a source external to the ocean-atmosphere system and that the record of d11B change across the PETM in the North Atlantic, due to its consistency with the existing Pacific record, is a good estimate of the magnitude of global surface ocean acidification. Many PETM records collected over decades substantiate the argument about the whole ocean d13C perturbation, and model-data comparison of simulated pH changes at the locations corresponding to the d11B records substantiate the application of those pH changes as global constraints. I don’t think that the evidence for the POE reaches the same bar, so I think too much weight is being given to interpreting the results from an experimental framework that I’m not convinced is valid. The POE has not been conclusively identified in deep-sea records, and despite the authors comparison to features in the Site 1209 and Site 1263 benthic foraminiferal records of the late Paleocene (not shown and no correlation attempted), it’s not even clear that the POE is identifiable in the deep sea. This was the basis of my previous major criticism of this paper, and I don’t think the revision has adequately addressed this problem. A modeled <1‰ d13C excursion is significant, and a -1.5‰ excursion doesn’t compare well to a -0.4‰ benthic excursion, with no quantification of how bioturbation should impact the modeled core-top d13C excursion magnitude (moreover, the issue of experiment duration in recording the core-top excursion magnitude has not been addressed). The authors modeling also imposes a pH change greater than that reconstructed across the PETM and consequently results in a CO2 increase that is larger than independent estimates for the PETM CO2 change (not those from boron isotopes, see Tierney et al., 2022).\nThe authors are using the d11B record from Babila et al. in their inversion, but that paper provides significantly more caveats in the interpretation of d11B as pH across the POE. They write, “the determination of the absolute magnitude of the pH excursion at the POE is limited by both the uncertainty of the initial environmental conditions and the larger uncertainty associated with the laser ablation analyses generated in this study…which makes determining the absolute magnitude of the pH excursion at the POE rather uncertain.” And also, “a full propagation of uncertainty…allows us to determine that the POE pH excursion at the Maryland sites was greater than -0.08, with an upper limit that is poorly constrained but overlaps with the magnitude estimated for the main CIE observed elsewhere.” Two important points I want to emphasize from these quotes: first, the authors here chose to use something close to the poorly constrained upper limit for their pH inversion target (-0.4) in the main results presented in this text but do not provide a justification for this choice. Second, the Babila paper correctly differentiates the pH change observed in the coastal North Atlantic to the pH decline for the CIE observed elsewhere (emphasis mine). I can think of a handful of reasons that the changes in pH near coastal sites and/or more restricted areas (as in this study location) would differ from the open ocean pH change. Changes in freshwater input, for instance, could account for some amount of local pH change. Hence using a -0.4 pH change in the double inversion seems to me as though the authors are choosing their forcing in order to arrive at a particular result. I appreciate that the authors have performed sensitivity analyses, so I think this entire discussion could be reworked largely using experiments already conducted, even though the authors did not test the minimum pH change estimate of Babila et al. (-0.08). I don’t think the “preferred scenario” is well-justified. Moreover, shouldn’t the discussion acknowledge and explain contrasts between their results and the modeling results presented by Babila et al., who evaluated a large range of pH changes and POE duration estimates, including much shorter durations than the minimum 1000 years tested here? Also, could the authors provide more discussion of why these age models differ from the age model developed for the Bighorn Basin terrestrial record from Bowen et al. 18? (The authors describe their estimated duration of 7 kyr as similar but slightly longer than” the Bowen estimate, but in reality their estimate is more than double the Bowen estimate, of “fewer than 2,000 years” (Bowen et al., 2018). Also, a more clear acknowledgement that Babila et al also suggested that “it is plausible that a short-lived C emission pulse preceded the main pulse of volcanic emissions” is warranted.\nAlso, regarding the modeled d13C target: could the authors include a plot in the supplement that shows a high-resolution image of all the d13C records across the POE alone for comparison to the target used in the inversion? Also, were the yellow dots in figure 4 actually the inversion target? Or was the inversion target a smoothed version of this curve?\nReference cited: Tierney, J.E., Zhu, J., Li, M., Ridgwell, A., Hakim, G.J., Poulsen, C.J., Whiteford, R.D., Rae, J.W. and Kump, L.R., 2022. Spatial patterns of climate change across the Paleocene–Eocene Thermal Maximum. Proceedings of the National Academy of Sciences, 119(42), p.e2205326119.\nLine 57 - add reference to Tierney et al., 2022.\nLine 96 – ‘The’ primary d13Ccarb signal…\nLine 113 – ‘average’ rather than ‘are averaging’\nLine 109: more about the age model – the 1.2 – 1.9 m peaks show the strongest power? These are well-identified throughout each of the POE and PETM? There is error on the duration between the POE and PETM of +/- 1 precession cycle but no error within the PETM? Presumably the 21-kyr duration estimates are consistent with identification of a single precession cycle, so please indicate that precession has an assumed 21 kyr duration on Line 112 (instead of ~20 kyr).\nLine 119 – how can you tell that it is the recovery phase in particular that is truncated? What about alignment based on features in the d13C records?\nAge model options 1 and 2 seem like extremes – what about 3 and 5 m cycles? What are these?\nLine 151 – to support the possibility of diagenetic overprinting or freshwater influence, what is the magnitude of the d18O excursion?\nLine 178 – say something about the significance of a reduction in Thaumarchaeota?\nLine 255 – perhaps add a sentence here specifying how the carbon emissions rates and sources are re-assessed? As in, substituting the new KZGS d13C record but using existing d11B constraints from Gutjahr?\nLine 263 – is this really a data assimilation approach?\nLine 274 – the POE change in pH is 0.4 units? As in, larger than the PETM pH change of 0.3?\nLine 279 – this could use a bit more explanation; specifically that you forced the model to follow the d13C excursion in DIC to the minimum but then not through the recovery, such that it was the modeled recovery mechanisms that controlled the mismatch between the modeled d13C of DIC during the recovery from the target curve…\nLine 302 – please provide the size of the simulated benthic d13C excursion\nLine 304 – so why not look at the downcore record?\nLine 379 – is this an argument against terrestrial organic carbon burial for the POE? Why would the PETM argument hold for the POE?\nLine 405-406 – are these two transient warming events? Based on the authors discussion, it seems like they are interpreting the entirety of the record as phases of NAIP carbon release.\nLine 412 – why round up the organic carbon burial estimate; use the same value given above (2500-2700 Pg C)"}, "author_response": "The authors have done an excellent job of addressing previous concerns I had about the presentation and interpretation of their exciting data set. The revised manuscript is much more polished and well presented. I only have a couple of small comments and a few grammar suggestions, otherwise I believe it is now ready for publication.\nReply: Thank you for the suggestions for improving the manuscript, which we have provided detailed responses below.\nResponse to reviewers: The authors did a great job of addressing my previous comments. One slight error is that they continue to advocate for a potential PETM age of 55.53 Ma from Westerhold et al. (2007) as possible evidence for the POE to be after the current dates we have for the big upswing in NAIP activity. Multiple lines of evidence, including an ash layer within the PETM dated to 55.785 Ma (Charles et al., 2011), rule out a PETM onset age of 55.5 Ma. This doesn’t make much difference to the manuscript as the focus has now shifted away from the precise timing of the POE, but I thought it was worth mentioning again.\nReply: We agree with the reviewer that although the focus is not on the precise timing of the POE, the age of the PETM was worth mentioning again. Therefore, we added three references that discuss the age of the PETM (~56 Ma) in Introduction. The three references are the following:\n1. Charles AJ, et al. Constraints on the numerical age of the Paleocene-Eocene boundary. Geochem Geophys Geosyst 12, Q0AA17 (2011).\n2. Röhl U, Westerhold T, Bralower TJ, Zachos JC. On the duration of the Paleocene-Eocene thermal maximum (PETM). Geochemistry, Geophysics, Geosystems 8, Q12002 (2007).\n3. Zeebe RE, Lourens LJ. Solar System chaos and the Paleocene–Eocene boundary age constrained by geology and astronomy. Science 365, 926-929 (2019).\nThe final piece of the puzzle that I think the manuscript needs prior to acceptance is that there is now a testable hypothesis from their findings. 0.2 to 1.3 Pg C/yr (line 304) is a sizeable flux of carbon to the ocean-atmosphere system, and if the NAIP is the dominant source, this should be discernible in the geological record. With more accurate dating of NAIP rocks, this may confirm the NAIP as a possible source. The discussion and/or conclusion could be improved with a sentence or two along the lines of “our findings predict substantial carbon fluxes driving the POE, and refining geochronological investigations of potential sources such as the NAIP could corroborate or refute this hypothesis.”\nReply: We now add the suggested sentence in the concluding paragraph, and the new text now reads as the following:\n“Furthermore, our findings predict substantial carbon fluxes driving the POE (averaging 0.2 to 1.3 Pg C yr-1), which could be tested by refined geochronological investigations of potential sources such as the NAIP.”\nGrammar edits (suggested additions in caps):\nLine 40: Add “(NAIP)”\nLine 45: Add “confirming THE global…”\nLine 50: Add “driven by A largely…”\nLine 51: Add “associated with sill intrusionS prior…”\nLine 52: Add “involving THE release…”\nLine 62: Suggest replacing “twenty millenia” with “20 kyr”\nLine 70: Suggest replacing “preceding” with “prior to”\nLine 71: Add “PETM ONSET with…”\nLine 73: Add “preservation to surface AND SHALLOW WATER records…”\nLine 74: Add “Resolving A global POE signal…”\nLine 83: Add “that volcanism AND MAGMATISM may also…”\nLine 84: Add “the global extent OF THE POE…”\nLine 104-105: Describe what the depth scale of the stratigraphic log is relative to when first referred to in the text.\nLine 181: Add “weathering AND/or higher terrestrial…”\nReply: All of these suggested changes are made in the revised text.\nLines 260-264: This needs rewording. Mantle inputs of osmium have low values of 187Os/188Os due to less rhenium in the source material that has decayed to 187Os. Low Os isotopes are therefore more unradiogenic, not the other way around.\nReply: The new text that addresses reviewer’s comment is replicated below:\n“A negative shift in 187Os/188Os ratios has been observed prior to the PETM in several sites globally32, 73, 78, 79, lending support to the occurrence of large igneous province activity prior to the PETM.”\nLine 345: Specify here whether you are referring to volcanic, thermogenic, or a mixture of these sources from the NAIP here.\nReply: This is now clarified, and the new text reads as the following:\n“Although this study provides a range of estimates on the carbon source and emission flux during the POE, more precise δ11B-based global surface pH records, detailed history of the sill intrusion of the NAIP, sea surface temperature records from across different latitudes, and better-constrained geochronology of the NAIP activity are clearly needed to reduce the uncertainty of the estimated thermogenic carbon emission fluxes from the NAIP.”"}]

Decision letter