Dynamics of sediment flux to a bathyal continental margin section through the Paleocene–Eocene Thermal Maximum

The response of the Earth system to greenhouse-gas-driven warming is of critical importance for the future trajectory of our planetary environment. Hyperthermal events – past climate transients with global-scale warming significantly above background climate variability – can provide insights into the nature and magnitude of these responses. The largest hyperthermal of the Cenozoic was the Paleocene–Eocene Thermal Maximum (PETM ∼ 56 Ma). Here we present new high-resolution bulk sediment stable isotope and major element data for the classic PETM section at Zumaia, Spain. With these data we provide a new detailed stratigraphic correlation to other key deep-ocean and terrestrial PETM reference sections. With this new correlation and age model we are able to demonstrate that detrital sediment accumulation rates within the Zumaia continental margin section increased more than 4-fold during the PETM, representing a radical change in regional hydrology that drove dramatic increases in terrestrial-to-marine sediment flux. Most remarkable is that detrital accumulation rates remain high throughout the body of the PETM, and even reach peak values during the recovery phase of the characteristic PETM carbon isotope excursion (CIE). Using a series of Earth system model inversions, driven by the new Zumaia carbon isotope record, we demonstrate that the silicate weathering feedback alone is insufficient to recover the PETM CIE, and that active organic carbon burial is required to match the observed dynamics of the CIE. Further, we demonstrate that the period of maximum organic carbon sequestration coincides with the peak in detrital accumulation rates observed at Zumaia. Based on these results, we hypothesise that orbital-scale variations in subtropical hydro-climates, and their subsequent impact on sediment dynamics, may contribute to the rapid climate and CIE recovery from peak-PETM conditions.

1 Introduction

The PETM is a marked interval of climate warming (Dunkley Jones et al., 2013), shoaling of oceanic carbonate saturation horizons (Zachos et al., 2005) and a large global negative carbon isotope excursion (CIE) in all rapidly exchangeable marine and terrestrial carbon reservoirs (McInerney and Wing, 2011). The most consistent explanation for these coupled perturbations is the release of carbon from a large shallow lithospheric reservoir, with a depleted carbon isotopic (δ¹³C) signature (∼ −20 to −60 ‰), on a multi-millennial timescale (Bowen et al., 2015; Dickens et al., 1995; Gutjahr et al., 2017; Kirtland Turner and Ridgwell, 2016; Zeebe et al., 2016). Although there is still no confidence on the identity of such a large (> 4000 Pg C) and unstable carbon reservoir, its release and oxidation within the ocean–atmosphere system caused rising atmospheric CO₂ concentrations, warming and a range of Earth system perturbations associated with pronounced global warming (Sluijs et al., 2007).

Although considerable attention has been paid to constraining the rates of carbon release, based on deep-ocean carbonate dissolution (Panchuk et al., 2008; Zachos et al., 2005; Zeebe et al., 2009), rates of warming (Meissner et al., 2014; Zeebe et al., 2016), carbon isotope profiles (Bowen et al., 2015; Kirtland Turner and Ridgwell, 2016) and surface ocean pH (Gutjahr et al., 2017) the mechanisms responsible for both the climatic and isotope recovery at the end of this transient event are still not well constrained (Bowen and Zachos, 2010). The timescales of silicate weathering and carbonate burial (∼ 100–200 kyr) are suggested to be too long to drive the main phase of CIE recovery, which the best records available to date indicate is an order of magnitude faster (∼ 10–20 kyr) (Bowen and Zachos, 2010). However, observed patterns of enhanced carbonate and possibly biogenic silica deposition in the deep ocean through the PETM recovery phase are interpreted as the signature of a silicate weathering-driven recovery (Penman, 2016; Penman et al., 2016). The alternative mechanism for carbon removal is the substantial burial of organic carbon within either terrestrial or marine sedimentary systems (Bowen and Zachos, 2010). One means of increasing global organic carbon burial rates is through enhanced sedimentation rates on continental margins, which might contribute significantly to C_org burial rates with little or no rise in measured total weight percent organic carbon content of sediments (Sluijs et al., 2014). There is also direct evidence for substantial increases in the organic carbon content of syn-PETM marine sediments in the eastern Tethys, associated with increased terrestrial to marine weathering fluxes and reduced seawater oxygenation (Carmichael et al., 2017; Gavrilov et al., 1997).

If organic carbon burial does drive recovery from the PETM, there is no explanation for the temporal offset between increased sedimentation rates and organic carbon burial during the body of the PETM (John et al., 2008) – when models indicate no or little carbon removal from the systems – and the major phase of carbon burial required for the CIE recovery at the end of the event. Models that include a continued leakage of carbon to maintain a prolonged CIE provide a mechanism to delay recovery (Penman et al., 2016; Zeebe et al., 2009), but this delay only exacerbates the mismatch between the observed fast recovery and the timescales of silicate weathering-driven carbon drawdown (Penman et al., 2016). As yet, there is no satisfactory model to explain the timing of the large-scale carbon drawdown required to drive the rapid CIE recovery phase ∼ 100 kyr after the onset of the PETM (Bowen and Zachos, 2010). Understanding the ultimate fate of thousands of petagrams of “excess” carbon injected into the ocean–atmosphere system during the PETM has profound implications for understanding the behaviour of the Earth's climate future over the coming millennia. Here, we present new high-resolution bulk carbonate isotopic analyses, calcium carbonate and organic carbon contents, and bulk sediment elemental data from the Iberian upper continental slope (∼ 1000 m paleo-water depth) to address the linkage between climatic controls on sedimentation rates and carbon drawdown. We reference this data within a new age model, constrained by detailed stratigraphic correlations to the most complete terrestrial and marine PETM reference successions.

2 Location and methods

Our studied PETM section is located at the Itzurun beachfront of the town of Zumaia, northern Spain (Fig. 1; 043^∘18^′4.5^′′ N, 002^∘15^′31.2^′′ W). Situated within the deep basin of the Pyrenean Gulf, it is the deep-water end of a sediment distribution system that has been mapped in detail over the past 2 decades (see references and summary within Pujalte et al., 2015). Logged sections span terrestrial alluvial fans and plains, shallow shelf carbonates, incised channel distribution systems and a deep-water basin (Pujalte et al., 2015, and references therein). Within this, the Zumaia section is regarded as the most complete and representative end-member of the deep-water depositional system, recording the distal deposition of fluvial-derived finely grained sediment plumes (Pujalte et al., 2015) at middle–lower bathyal (1000 m) paleo-depths (Alegret et al., 2009).

Figure 1 Location of Zumaia section and relationship to the terrestrial and shallow shelf carbonate depositional environments, after Pujalte et al. (2015).

The Zumaia section is a classic PETM locality and the subject of a series of bio-, magneto-, chemo- and lithostratigraphic studies (Alegret et al., 2018, 2009; Baceta et al., 2000; Bernaola et al., 2007; Canudo et al., 1995; Dinarès-Turell et al., 2007, 2002; Manners et al., 2013; Schmitz et al., 1997, 2001; Storme et al., 2012). The aim of this study is to build on the existing detailed depositional framework developed for the region, and to use the Zumaia section as a key recorder of the temporal patterns of climate-driven sediment flux from terrestrial systems to the deep ocean during the PETM. To do this we sought to substantially increase the resolution of existing bulk carbonate isotope records and supplement these with high-resolution bulk sediment compositional analyses using X-ray fluorescence (XRF). This approach is consistent with a progressive shift in PETM studies, from broad characterisation of the CIE as a marker of the event within wider Paleogene stratigraphy, to a focus on the millennial-scale dynamics of the carbon cycle, climate and other Earth system properties within the PETM itself (Bowen et al., 2015; Zeebe et al., 2016). For this study we took a total of 248 small “thumbnail” 1cc bulk sediment samples at ∼ 3 to 5 cm resolution from ∼ 4.3 below to ∼ 8.2 m above the onset of the characteristic “siliciclastic unit” (SU) that marks the body of the PETM at this location. The zero reference point for all sampling is taken at the top of the prominent limestone bed, located ∼ 0.3 m below the base of the SU. In between the limestone and the SU is a 0.3 m thick greenish-grey marl unit, which is distinct from both the limestone–marl units found above and below the SU and the typically red-brown clays and marly clays of the SU itself. The section was logged in detail and integrated with the existing stratigraphic framework for the site (Fig. 2).

All samples were analysed for bulk sediment elemental composition with XRF using dried, powered sediment, mounted in wax pellets. Pellets were analysed with a Bruker S8 TIGER XRF spectrometer with 8 min analysis time, at the University of Birmingham (UoB). CaO concentrations from XRF analyses were converted into high-resolution estimates of weight percent CaCO₃ using a calibration set of 25 discrete calcium carbonate measurements. Weight percent CaCO₃ contents were determined at UoB, based on inorganic carbon contents measured by CO₂ generation with phosphoric acid reaction within a 200 ^∘C furnace using a Shimadzu TOC-V carbon analyser. Detrital mass accumulation rates (MARs) were estimated assuming that both detrital and carbonate fractions had an estimated dry bulk density of 1.8 g cm⁻³. Organic carbon MARs were calculated following Willmott et al. (2009).

A subset of 129 thumbnail samples were also analysed for the stable isotope composition (δ¹³C_carb and δ¹⁸O_carb) of bulk carbonates at the NERC Isotope Geosciences Laboratory (NIGL, Keyworth, Nottingham). Sample material was powdered to yield CO₂ from the equivalent of 10 mg of 100 % calcite and reacted with anhydrous phosphoric acid in vacuo overnight at a constant 25 ^∘C. Liberated CO₂ was separated from water vapour under vacuum and collected for analysis. Measurements were made on a VG Optima mass spectrometer (standard reproducibility of < ± 0.1 ‰). All data are on the Vienna Pee Dee Belemnite scale. An additional 67 bulk samples were analysed for total organic carbon content at NIGL using a Carlo Erba 1500 elemental analyser with acetanilide used as the calibration standard. Replicate analyses indicated a precision of ±0.1 wt % in well-mixed samples (1 standard deviation, SD). All data are available in the Supplement.

3 Results

3.1 Carbon isotope stratigraphy

The new high-resolution bulk carbonate isotope stratigraphy presented here is consistent with previous studies of the Zumaia section, in particular the isotopic data of Schmitz et al. (1997). The capacity for the analyses of low-weight carbonate samples at NIGL has, however, provided substantial additional data within the critical phase of most rapid change at the PETM onset, including the base of the SU, between ∼ 0.3 and 0.7 m. These new data provide information about subtle lead–lag relationships between lithological change at the base of the SU, sediment carbonate content, the CIE, and benthic foraminifera community change and extinction (Alegret et al., 2018). Although there is a clear negative δ¹³C shift of ∼ 1 ‰ between the prominent limestone unit (−0.7 to 0 m) and the 0.3 m thick marl unit (0 to 0.3 m) that immediately underlies the SU, this change and the absolute δ¹³C values concerned are all within the background variations in δ¹³C seen in the marl–limestone couplets of the latest Paleocene (e.g. from −4.0 to −1.0 m). For example, the carbonate δ¹³C values of ∼ 0 ‰ within the precursor marl unit (0 to 0.3 m) are close to the average value for the latest Paleocene succession analysed in this study. On this basis, we take a conservative approach and do not identify this δ¹³C shift as related to the onset of the main phase of the PETM. This is not to say that there are not significant environmental perturbations through the limestone and precursor marl intervals. There is evidence for an early shift in both benthic (Alegret et al., 2018) and planktonic foraminiferal assemblages within this marl unit (Schmitz et al., 1997). Further, the distinct precursor limestone bed is itself an unusual feature within the typically thinner and lower CaCO₃ limestone–marl couplets of the latest Paleocene. The nature of any such precursor environmental changes is beyond the scope of this study, but it does underline the potentially great importance and utility of the Zumaia section because, unlike many deep-ocean sites, it appears to be more continuous and expanded across the onset of the PETM.

Here, we place the base of the CIE at 0.30 m. The δ¹³C shift between 0.25 m (0.01 ‰) and 0.30 m (−0.28 ‰) is the start of the consistent negative deviation in δ¹³C values that continues through the base of the SU. This level is also the point of lithological change between the precursor marl and the base of the SU. The largest negative step in δ¹³C – from −0.8 to −3.6 ‰ – is recorded between samples 0.34 and 0.38 m with a further decline to values of −4.0 ‰ or lower above 0.5 m (Fig. 2). The total CIE is more than ∼ 3.7 ‰ (Fig. 2). We note that there is a slight lag between the onset of the CIE and the corresponding decline in calcium carbonate contents. Carbonate contents are stable between 0.25 and 0.30 m, at ∼ 44 wt %, but then drop rapidly across successive samples to 26 wt % (0.34 m) and then 7 wt % (0.38 m). The placement of the onset of the CIE at 0.30 m also agrees well with the main phase of the benthic foraminifera extinction event (BEE), which was placed at ∼ 0.3 m in the original study of Schmitz et al. (1997) as well as in detailed analyses of foraminiferal turnover (Alegret et al., 2018, 2009). A gradual but rapid disappearance of benthic foraminifera is recorded through the precursor marls (0 to 0.30 m), but the main phase of the BEE coincides with the onset of the negative CIE at ∼ 0.3 m (Alegret et al., 2018).

Figure 2 New sediment geochemistry data for the Zumaia section and correlation to existing deep-ocean and terrestrial records (inter-site correlations and age model of Westerhold et al. 2018). Panels (a) to (e) show data against depth in the Zumaia section: (a) lithological log after Storme et al. (2012); (b) lithological log produced during sample collection for this study (2011); (c) SiO₂ weight percent (green) and Si ∕ Fe ratio (cyan) based on weight percentages of oxides (SiO₂ ∕ Fe₂O₃) from bulk sediment XRF analyses; pre-PETM and post-PETM precession cycles numbered in black; (d) weight percentage of CaO and CaCO₃ from XRF; (e) new bulk carbonate δ¹³C (this study). Panels (f) to (l) show data against age (Westerhold et al., 2018): (f) Polecat Bench core (PCB) soil carbonate nodule δ¹³C (green; Bown et al. 2015); (g) deep-ocean benthic foraminiferal δ¹³C (Bains et al., 1999; McCarren et al., 2008); (h) Gaussian filter of the longer 8.2 m cycle (precession) of the Fe data from PCB-B (red; Westerhold et al., 2018); (i) ODP Site 690 XRF core scanning Fe data (red) and ODP Site 1262 A* colour reflectance data (black); (j) composite core section photos from deep-ocean sites ODP 1262, 1267, 1266, 1265, 1263 and 690; (k) terrestrial PCB sites A and B, with the latter correlated to (l) nearby outcrop lithologies (Westerhold et al., 2018).

Although an increase in sedimentation rates at the start of the PETM in the Zumaia section limits the effects of any carbonate burn-down on the earliest stages of the CIE profile, the measured δ¹³C_carb data diverge from the typical PETM CIE profile between ∼ 0.7 and 2.8 m. Here, δ¹³C_carb values return to higher values, of ∼ −3 ‰, which lie between those of the peak CIE and pre-excursion values (Fig. 2). We attribute this to the increased relative contribution of reworked detrital carbonates, as evidenced by the presence of Cretaceous nannofossils (Tom Dunkley Jones, personal observation, 2017), during periods of very low autochthonous carbonate flux. This is supported by the return of δ¹³C_carb to near peak-excursion values once carbonate contents increase above 2.8 m. Above this level, there is a CIE recovery phase similar to other high-resolution records from the terrestrial section and the deep ocean (e.g. Bowen et al., 2015; Zachos et al., 2005; Fig. 2).

In order to develop a refined age model for the Zumaia section, we used our new high-resolution stable isotope stratigraphy, together with the established cyclostratigraphy for the Zumaia section, based on precession-paced lithological cycles (Dinarès-Turell et al., 2003; Westerhold et al., 2008). This combination of carbon isotope and cyclostratigraphy provided multiple direct stratigraphic ties to existing deep-ocean and terrestrial reference sections (Bowen et al., 2015; Röhl et al., 2007; Westerhold et al., 2018). These tie points include precession cycles −9 to −14 below the onset of the PETM, and cycles 6 to 14 above the onset of the PETM, as well as the CIE onset and marked inflexion point in δ¹³C during the recovery phase (∼ 55.81 Ma). In the latest Paleocene, distinct marl–limestone couplets are lost with the deposition of more continuous higher carbonate lithologies between ∼ −1.0 and 0 m, including the distinct marker limestone bed. This high carbonate interval may be due to clastic sediment starvation and has previously been associated with the base of a transgressive systems tract (Storme et al., 2012). The new bulk sediment chemistry presented here is also unable to clearly identify cycles within this interval, but interpolations to Chron 25n and biostratigraphic correlations between Zumaia and deep-ocean sections would suggest that this high carbonate interval is somewhat condensed relative to the rest of the succession. In this study, we associate the first clear limestone–marl couplet beneath this high carbonate interval (marl between −1.4 and −1.8 m) with the ninth precession cycle below the onset of the CIE. Above the SU, two prominent thicker carbonate beds are correlated with the 100 kyr eccentricity minima of Westerhold et al. (2007; precession cycles 6/7 and 11), which is supported by the placement of the last rare occurrence of the calcareous nannofossil genus Fasciculithus at 8.7 m (Baceta et al., 2006; confirmed by Tom Dunkley Jones, personal observation), between cycles 15 and 16, which is consistent with the Forada section and ODP Site 690 (Röhl et al., 2007).

3.2 Sedimentary cycles through the PETM

Precession-paced lithological cycles are well established within the Zumaia succession, both above and below the main body of the PETM (Fig. 3; Dinarès-Turell et al., 2003; Westerhold et al., 2008). This same prominent cyclicity is also resolved, both above and below the SU, in the new bulk sediment chemistry data presented here, including in calcium carbonate and SiO₂ contents and Si ∕ Fe ratios (Fig. 3). In addition to these already established precession cycles, we note, for the first time, a prominent variation in Si ∕ Fe ratios through the SU itself. Spectral analysis, using the “astrochron” software package (Meyers, 2014), of the new high-resolution Zumaia Si ∕ Fe records resolve significant power at frequencies of ∼ 2.2 cycles m⁻¹ (∼ 0.45 m period), above the 99 % confidence limit (Fig. 4). Evolutive power spectra through the whole record show that these cycles are most strongly developed within the body of the PETM. In total there are 10 of these cycles in Si ∕ Fe between the CIE onset and the δ¹³C inflection point F of Röhl et al. (2007; Fig. 3). With a total CIE duration of 90–167 kyr (Murphy et al., 2010; Röhl et al., 2007) this variability is on the half-precession timescale ($\mathrm{10}\cdot \mathrm{10.5}=\mathrm{105}$ kyr) and bears a striking resemblance to variations observed in terrestrial paleosol successions from the Bighorn Basin (Abdul Aziz et al., 2008; Westerhold et al., 2018). In the Bighorn Basin these cycles are related to changes in catchment hydrology (e.g. Abdul Aziz et al., 2008), resulting from the interference of both Northern Hemisphere and Southern Hemisphere precession signals in tropical hydro-climates (Short et al., 1991). A similar mechanism for the origin of the Zumaia intra-PETM cycles is consistent with a strong hydroclimate control on sediment erosion and transport in the Tremp-Graus Basin during the PETM (Pujalte et al., 2015; Schmitz and Pujalte, 2007).

Figure 3 New sediment geochemistry data for the Zumaia section: (a) lithological log produced during sample collection for this study (2011); (b) bulk carbonate δ¹³C; (c) weight percentages of CaO and CaCO₃ from bulk sediment XRF analyses; (d) Si ∕ Fe ratio based on weight percentages of oxides (SiO₂ ∕ Fe₂O₃) from XRF analyses with pre-PETM and post-PETM precession cycles numbered in black, intra-PETM cycles numbered in red; (e) Si ∕ Fe ratios as in (d) with Gaussian filters of the dominant cycles 0.45 m (red) and 1 m (thin blue line); (f) evolutive spectral analysis for Si ∕ Fe.

In our preferred analysis, we use these 10 intra-PETM Si ∕ Fe cycles to provide additional age constraints through the PETM at Zumaia. The use of these cycles does not impact the wider detailed correlations to other locations based on carbon isotope stratigraphy and cyclostratigraphy outside of the body of the PETM as presented above. We believe these correlations are robust and provide a new context for interpreting the Zumaia record. The presence of distinct variability in Si ∕ Fe ratios though the body of the SU at Zumaia provides additional, fine-scale lithological and bulk rock geochemical data for an otherwise poorly resolved interval. The use of XRF-derived Fe contents of bulk sediments is established practice for the development of cyclostratigraphic models within Paleogene successions (e.g. Röhl et al., 2007; Westerhold et al., 2007, 2008, 2014, 2018; Zachos et al., 2010). Although other sections within the Basque Basin show distinct lithological variability through the PETM, including turbidite deposition at the up-dip and shallower-water section of Ermua (Schmitz et al., 2001) and the development of large deep-water distribution channels (Pujalte et al., 2015), no evidence for either turbidite deposition or any significant unconformity has been reported for the section at Zumaia. Indeed, the Zumaia section is regarded as the most complete and representative PETM succession within the Basque Basin (Pujalte et al., 2015). Throughout the main phase of the PETM the Zumaia section is dominated by finely grained sediments, with the < 28 µm fraction typically contributing 99 % or more of sediment volume during the SU (Schmitz et al., 2001). This is consistent with the absence of any recognised turbidite deposition during the PETM phase at Zumaia (Clare et al., 2015). In summary, all available evidence indicates that deposition is continuous and dominated by finely grained hemi-pelagic clays through the PETM interval at Zumaia. It is this continuity that makes the Zumaia section an important location for the detailed study of biotic and geochemical perturbations through the PETM.

Figure 4 Multitaper method (MTM) power spectra for Si ∕ Fe ratios; highlighted are spectral peaks equivalent to cycles with a period > 1 and ∼ 0.45 m. Autoregressive spectral analysis undertaken using the mtmML96 routine of the “astrochron” software package (Meyers, 2014); dotted red lines represent the 90, 95 and 99 % confidence levels (CL).

3.3 Calcium carbonate and detrital sediment mass accumulation rates

Based on our age model and estimated weight percentage of CaCO₃, we estimate calcium carbonate, total organic carbon (TOC) and detrital MAR throughout the study interval (Figs. 5 and 6). Carbonate MAR closely follows the weight percentage of CaCO₃, with a marked decline during the body of the SU. These two records of carbonate content do, however, diverge in two key intervals: first during the distinct limestone bed below the onset of the PETM, between ∼ −0.8 and 0 m (∼ 56.03 to 55.95 Ma), and then in the later stage of the CIE recovery between ∼ 4.4 and 4.7 m (55.82 to 55.81 Ma) (Figs. 5 and 6). In the first of these, the distinct high carbonate content of the limestone unit is weakly represented in carbonate MAR, indicating that this is the result of clastic sediment starvation, rather than increased carbonate production. The second horizon occurs during the CIE recovery phase, when carbonate MARs nearly treble compared to pre-PETM values before stabilising at levels ∼ 50 % higher than pre-PETM conditions during the post-PETM phase above 4.7 m (∼ 55.81 Ma). This persistent increase in carbonate MARs above the PETM is reflected in consistently higher carbonate contents (CaCO₃ wt %) of sediments above the PETM than those in the latest Paleocene. This is clear evidence for significantly increased carbonate accumulation rates post-PETM, even in relatively shallow (∼ 1000 m) supra-lysoclinal depths on the continental margins. The rate and timing of this recovery in carbonate MARs is remarkable, with the recovery to, and then exceedance of, background levels all within the ∼ 20 kyr of the main phase of the CIE recovery (Fig. 6).

Figure 5 New carbonate, detrital and organic carbon mass accumulation and sedimentation rates for the Zumaia section: (a) lithological log produced during sample collection for this study (2011); (b) bulk carbonate δ¹³C; (c) weight percentage of CaCO₃ derived from XRF CaCO calibrated to discrete CaCO₃ analyses (see methods); (d) CaCO₃ mass accumulation rate; (e) detrital mass accumulation rate; (f) organic carbon accumulation rate; (g) linear sedimentation rates derived from the new cyclostratigraphic framework.

The MAR of non-carbonate components is here considered equivalent to the detrital MAR, in the absence of any significant contribution of biogenic silica. Before the CIE, the detrital MARs are relatively low, fluctuating around 2 g cm⁻² ka⁻¹, with a further decline during the precursor limestone unit (−0.8 to 0 m). During the onset of the PETM (0.3–1.5 m) detrital fluxes increase dramatically, to between 3 and 4 times pre-PETM levels. Within the peak PETM (1.5–1.7 m; ∼ 55.88 to 55.85 Ma), detrital flux declines slightly and shows variability on similar timescales to Si ∕ Fe ratios. It then increases markedly towards the end of the peak-PETM interval, above 3.7 m (∼ 55.85 Ma), and reaches peak values for the whole PETM during the main phase of the δ¹³C recovery interval (3.7–4.2 m; ∼ 55.83 to 55.82 Ma). Both detrital and carbonate MAR increase to peak values for the whole PETM, during the main interval of CIE recovery, with detrital MARs peaking first, followed by carbonate MARs. We note that these high detrital MARs during the body of the PETM do not coincide with a noticeable increase in grain size, with sediments dominated by clay-grade material throughout. This is consistent with the marked absence of more coarsely grained turbidite deposition during the PETM within the Zumaia section (Clare et al., 2015), although increased coarsely grained siliciclastic deposition during the PETM is recorded within deep-water channels within the basin (Pujalte et al., 2015). Total organic carbon MARs are relatively stable before the onset of the event, at ∼ 0.01 g cm⁻² ka⁻¹. At the onset of the PETM there is an increase in TOC MARs, but this lags the CIE onset by ∼ 20 kyr (Fig. 6). During the body of the PETM there are distinct peaks in TOC MARs at ∼ 55.90 and ∼ 55.85 Ma, but in general values remain between 0.015 and 0.02 g cm⁻² ka⁻¹. Coincident with the CIE recovery phase, between 55.83 and 55.81 Ma, TOC MARs again increase to over ∼ 0.02 g cm⁻² ka⁻¹.

Figure 6 Comparison of the Zumaia records against terrestrial (Big Horn Basin) and deep-ocean sites (Walvis Ridge ODP Sites 1262, 1263, 1265, 1266, 1267; Maud Rise ODP Site 690) in the same cyclostratigraphic age model: (a) bulk carbonate δ¹³C for deep-ocean sites; (b) bulk carbonate δ¹³C for Zumaia; (c) soil nodule δ¹³C from the Big Horn Basin; (d) CaCO₃ records from Zumaia and the deep-ocean sites; (e) detrital (red) and CaCO₃ (blue) MARs from Zumaia; (f) organic carbon MAR from Zumaia.

4 Inverse modelling the Zumaia CIE and the PETM carbon source

To constrain the rates of carbon input and sequestration through the onset and termination of the PETM, based on the CIE profile and cyclostratigraphy from Zumaia, we undertook a series of PETM simulations using the late Paleocene–early Eocene configuration of the cGENIE Earth system model (Kirtland Turner and Ridgwell, 2016). This includes a 3-D dynamic ocean model, 2-D energy balance atmospheric model and a representation of marine biogeochemical cycling plus the preservation of carbonates in deep-sea sediments (Ridgwell and Schmidt, 2010). We also include a 0-D (global average) temperature-dependent terrestrial weathering feedback as described in Colbourn et al. (2013). We follow a two-part spin-up procedure, first for 20 kyr as a “closed” system (i.e. weathering exactly balancing sedimentation) and with atmospheric CO₂ set to 834 ppm with a δ¹³C of −4.9 ‰. A further 250 kyr “open” system spin-up with temperature-dependent weathering feedbacks enabled and input of CO₂ via volcanic outgassing (Colbourn et al., 2013) was used to bring the geological ¹³C cycle to steady state. In this, the specified non-carbonate detrital flux to marine sediments (0.18 g cm⁻² ka⁻¹), the rain ratio of carbonate to particulate organic carbon (0.2), and the global average concentrations of Mg²⁺ and Ca²⁺ (29.89 and 18.22 mmol kg⁻¹, respectively) and alkalinity (1975 µ mol kg⁻¹) were selected to give an average global mean sedimentary carbonate preservation of 46.75 %, comparable with late Paleocene estimates (Panchuk et al., 2008).

Four inversion experiments were undertaken, all forced with carbon (CO₂) inputs to the atmosphere such that the δ¹³C of surface ocean dissolved inorganic carbon (DIC) followed a prescribed inversion target taken from the Zumaia δ¹³C record, adjusted to the initial cGENIE average surface ocean DIC δ¹³C of +2.7 ‰. Within the very low carbonate interval, in which the primary signals appear to be lost due to the disproportionate influence of non-autochthonous carbonates, the δ¹³C target curve was linearly interpolated between the data points at ∼ 19 and 78 kyr after the onset of the PETM. Experiments included variation in the isotope composition of the source (−60 and −22 ‰) and CIE recovery driven by silicate weathering only and with active carbon removal to match the CIE profile (Table 1). In order to directly compare modelled and observed carbonate contents, all experiments also include the temporally varying observed changes in detrital fluxes within the sediment model at the model grid box corresponding to the paleo-location of Zumaia only.

Net input fluxes depend on the δ¹³C composition of the source (Table 1), with smaller carbon inputs required for a source with a more negative isotope composition. Total carbon inputs are also slightly larger where no removal fluxes are allowed (Experiment 1). In all experiments, the majority of C input occurs in the first 20 kyr with a low rate of C input (< 0.05 Pg C yr⁻¹) required to maintain low isotope values until 78 ka. Total modelled C input is broadly comparable to previous model-based estimates constrained by deep-ocean carbonate dissolution that exceed ∼ 4000 Pg C (Bowen et al., 2015; Cui et al., 2011; Gutjahr et al., 2017; Panchuk et al., 2008). In this case the lower bound of carbon input (4154 Pg C) is generated by an inversion of the CIE isotopic profile assuming a pure biogenic methane carbon source with isotope composition of −60 ‰. As an independent test of the reliability of these simulations, in terms of the estimated CIE magnitude and temporal dynamics, the modelled atmospheric δ¹³C excursion is compared to the best available terrestrial record of the CIE, from soil carbonates of the Bighorn Basin (Fig. 7; Bowen et al., 2015). There is a strong correspondence between the two, in the shape of the onset and recovery, the magnitude of the CIE, and the overall duration of the event. Although these model runs have no independent second constraint on the carbon cycle perturbation through the PETM – for example a measure of surface ocean pH change (Gutjahr et al., 2017) or deep-ocean carbonate saturation state (Panchuk et al., 2008; Zeebe et al., 2009) – model simulations with a mixed source input with δ¹³C of −22 ‰ are close in total mass input across the PETM onset (11 316 Pg C) to the preferred simulations of Gutjahr et al. (2017; ∼ 10 200 Pg C), which are constrained by surface ocean pH changes. The combination of these two records – pH from DSDP Site 401 and the CIE dynamics from Zumaia – is, however, informative. Total carbon input masses, constrained by the most detailed surface ocean pH records of the PETM available to date, are modelled to be > 10 000 Pg C (Gutjahr et al., 2017), but the dynamics of this release must be considerably quicker than assumed within the deep-ocean records recovered from DSDP Site 401 in order to satisfy the rate of CIE onset observed at Zumaia, which is ∼ 3.3 ‰ over less than 5 kyr. We argue that this CIE onset rate is robust, given that it is observed in both terrestrial (Bowen et al., 2015) and marine (this study) sections in which the onset of the CIE is expanded and set within a robust intra-PETM cyclostratigraphic age model. The CIE recorded at Zumaia, and within the Big Horn Basin (Bowen et al., 2015), of 4 ‰ or more, is also greater than that assumed in the favoured model simulations of Gutjahr et al. (2017) of 2.6 ‰, and is closer to their alternate run assuming a CIE of ∼ 4 ‰, which reconstructs a carbon source with a δ¹³C composition of ∼ −17 ‰. Both these aspects of the Zumaia record – a CIE onset of less than 5 kyr, and a CIE with a magnitude of ∼ 4 ‰ or more – requires the rapid release of carbon from a source that is more isotopically depleted than volcanic CO₂. Whether thermogenic methane generated by sill intrusion during the onset phase of the North Atlantic Igneous Province (NAIP) (Svensen et al., 2004) can meet the emission rates and carbon isotopic compositions required to generate the early stages of the CIE is yet to be resolved.

Table 1 Summary details of cGENIE simulations detailing isotopic compositions of source carbon input during the onset of the PETM to match the Zumaia CIE profile – and the nature of modelled CIE recovery, whether by silicate weathering only (Experiment 1), or by active carbon removal. The isotopic composition of this carbon removal (“Sink”) is also listed. Net C fluxes are calculated by the model to match the CIE profile and do not include carbon removed through silicate weathering.

5 Discussion: sediment accumulation rates and carbon removal during the PETM recovery

Two fundamental unknowns remain about the recovery phase of the PETM: first, what is the ultimate sink of thousands of petagrams of carbon when added to the ocean–atmosphere system as CO₂? Second, what controls the temporal dynamics of carbon flux into this sink, and thence controls the timing of the relatively rapid CIE and PETM climate recovery? The combination of the age-constrained CIE records from Zumaia, the inverse Earth system modelling of this CIE and the generally excellent agreement between the atmospheric CIE generated in these simulations with independent records from the terrestrial realm (Fig. 7; Bowen et al., 2015) provide new insights into these key questions. In simulations without active carbon removal, the only process restoring atmospheric pCO₂ to pre-CIE levels is the silicate weathering feedback, in which net removal depends on the size of the carbon input and associated temperature change. This “weathering only” experiment fails to match the observed CIE recovery, with minimal recovery of atmospheric δ¹³C${}_{{\mathrm{CO}}_{\mathrm{2}}}$ even when a significant mass of carbon has been removed from the atmosphere (Fig. 7). This is due to the isotopic imbalance between the ultimate carbon sink in this case – marine carbonates (−1 to 4 ‰) – compared to the carbon source (−22 or −60 ‰). With the weathering feedback only, the rate of removal of carbon from the atmosphere (Fig. 7), and associated climate recovery, also occurs on timescales significantly longer (> 100 kyr) than indicated by climate proxy data, where temperatures are closely coupled to the recovery of the CIE (Dunkley Jones et al., 2013). The configuration of the silicate weathering feedback within the model is thus an insufficient driver of PETM recovery on two counts: first its inability to drive recovery in the carbon isotopic composition of the exogenic carbon reservoir, and second, in its inability to match the rate of recovery of the climate system at the end of the PETM. The Zumaia record also provides evidence for a rapid climate transition at the end of the PETM, independent of carbon cycle system parameters (e.g. δ¹³C and CaCO₃ accumulation), in the observed rapid switch back to pre-PETM detrital fluxes during the recovery phase. This implies a rapid climate recovery – back to a hydrological and erosive system similar to pre-event conditions – that is very closely coupled to the CIE recovery (Fig. 6).

Figure 7 Results of inverse modelling with the cGENIE Earth system model: (a) the bulk carbonate δ¹³C for Zumaia which drives model inversions; (b) soil nodule δ¹³C record of the Big Horn Basin, as representative of atmospheric CO₂ δ¹³C for comparison to model results; (c) modelled δ¹³C composition of atmospheric carbon dioxide; (d) modelled and actual weight percentage of CaCO₃ content of the sediment column at Zumaia; (e) carbon emission rate in petagrams of carbon per year; (f) cumulative carbon emissions (cumulative input minus removal); (g) modelled atmospheric pCO₂. The model scenarios are as follows: (1) carbon input at δ¹³C of −22 ‰ during the onset, no carbon removal (green); (2) carbon input and removal with δ¹³C of −22 ‰ to follow Zumaia target CIE (red); (3) carbon input with δ¹³C of −60 ‰ and removal of carbon with δ¹³C of −22 ‰ to follow Zumaia target CIE (blue). All simulations have active temperature-dependent silicate weathering feedbacks operative throughout.

To correctly simulate both the CIE and the climate system (temperature, hydrology) recovery rates, the active removal of large masses of isotopically light carbon is required during the PETM recovery phase. Simple isotope mass balance requirements require a greater mass of carbon to be removed than was first released, if the isotopic composition of the sink is more positive (organic carbon −22 %) than that of the source (biogenic methane −60 %). In these methane input and organic carbon burial simulations (Experiment 4), the required excess carbon removal (removal mass minus input mass) causes a remnant post-event climate cooling, with an ∼ 60 % decline in post-PETM atmospheric CO₂ levels relative to pre-PETM conditions (Fig. 7), and a global mean surface temperature reduction of ∼ 4 ^∘C. At present there is no evidence for such a post-PETM cooling, which would suggest that the source and sink terms are relatively close in their isotopic compositions and precludes a biogenic methane input and organic carbon burial scenario.

The simulations that best match the constraints of the CIE and PETM recovery are those with the active removal of carbon, focused on the main phase of the CIE recovery, and with an isotopic composition close to that of the original carbon source. The failure of the silicate weathering feedback alone to match the rate and magnitude of CIE recovery is similar to the recent inverse modelling by Gutjahr et al. (2017), in which 2500 Pg of carbon removal, with isotopic composition of −30.5 ‰, was required to match the CIE. The larger mass of carbon removal required in the simulations presented here is due to the smaller-magnitude CIE recorded at DSDP Site 401 and the isotopically lighter composition of the carbon sink used (−30.5 ‰) by Gutjahr et al. (2017) compared to this study (−22‰).

The question remains as to what is the driver of enhanced carbon burial during the CIE recovery interval? In the Zumaia records it is notable that exactly this interval – the CIE recovery phase – is directly coupled with maxima in total sedimentation rates and detrital MARs (Fig. 6). There is also a close coupling between these peak detrital fluxes and a transient peak in carbonate MAR, which is exactly the coupling predicted if increased detrital accumulation is accompanied by the widespread burial of organic carbon, with a resultant rapid increase in ocean pH and carbonate saturation state (Bowen and Zachos, 2010). Inverse Earth system modelling of the CIE at this location demonstrates that this interval of peak MARs (∼ 80 to 120 kyr after onset) is coincident with ∼ 8700 Pg of carbon drawdown from the ocean–atmosphere system in the −22 ‰ δ¹³C input and removal scenario, equivalent to ∼ 87 % of total removal fluxes. The magnitude of increase in detrital accumulation rates at the end of the event (from < 6 to > 8 g cm⁻² ka⁻¹) is also of the same order as the ∼ 40 % increase in organic carbon burial rates estimated for the CIE recovery (Bowen and Zachos, 2010). Although Zumaia is a single PETM record from a single location, it has striking similarities to other Tethyan deep-water successions, such as Forada in Italy (Giusberti et al., 2016), and is likely representative of regional-scale dynamics of climate and the associated response of continental margin sedimentary systems. Within the main phase of the PETM, variability within the Forada section is interpreted to be driven by precession-paced wet and arid phases associated with a monsoonal-type subtropical climate (Giusberti et al., 2016). Measured TOC MARs through the Zumaia section are elevated during the recovery interval, ∼ 55.83 to 55.81 Ma (Fig. 6), although these data are not able to discriminate between the burial fluxes of old (pre-PETM) reworked carbon and those for newly produced (intra-PETM) organic carbon – only the latter will contribute to the removal of carbon from the exogenic carbon cycle. The Zumaia records cannot yet be used to trace the rates of burial of newly formed organic carbon through the PETM recovery phase; however they do provide a feasible mechanism for the timing of the recovery phase of the PETM, namely a strong precession forcing of tropical and subtropical hydrological regimes and associated sediment erosion, flux and accumulation. Alongside the observed precession-forced variations in sediment accumulation, there is likely a complex influence of temperature and CO₂ on the balance between rates of net primary production and organic carbon remineralization during the PETM (Cotton et al., 2015). This production–remineralization balance will control the availability of newly formed sedimentary organic carbon for long-term burial in the sedimentary carbon sink (Cotton et al., 2015). During conditions of peak PETM warmth, enhanced organic carbon remineralization in terrestrial systems may have drastically limited the flux to continental margin systems, but with even a gradual recovery from peak-warmth conditions – forced, for example, by known enhanced carbon burial in the anoxic marine basins of the eastern Tethys during the peak PETM (Carmichael et al., 2017) – a subtle shift in the production to respiration balance towards the availability of organic carbon for long-term burial, combined with a precession-forced peak in sediment flux and accumulation rates across terrestrial and continental margin systems, could account for the magnitude and timing of the sudden switch to widespread carbon sequestration required during the CIE recovery phase. The greatly enhanced detrital accumulation rates observed in continental margin settings such as Zumaia could also effectively mask patterns of enhanced burial of newly formed organic matter (Sluijs et al., 2014), unless the MARs of this newly formed organic matter can be constrained. Future studies should focus on discriminating the relative contributions of this newly formed (intra-PETM) and reworked fossil (pre-PETM) organic carbon within these continental margin successions, as only the former will contribute to carbon sequestration.

6 Conclusions

New high-resolution stable isotopic and sediment chemistry records across the PETM interval of the Zumaia section provide new insights into the dynamics of sediment accumulation as well as their relationship to local climate dynamics and global carbon cycle perturbations. In comparison to deep-ocean locations, the Zumaia section has enhanced sediment accumulation rates across the critical phase of PETM onset. This provides a high-fidelity and high-resolution record of the PETM CIE that is comparable to the most complete terrestrial and marine records. The magnitude and rate of CIE onset observed at Zumaia is in close agreement with soil carbonate nodule data from the Big Horn Basin and imply PETM onset times of less than 5000 years, with a CIE magnitude of at least ∼ 4 ‰. Detailed correlations of existing cyclostratigraphy from the Zumaia section, together with carbon isotope stratigraphy of the CIE, to other key marine and terrestrial locations with comparable records, provides a new age framework for the Zumaia PETM section. This age model is then used to determine detrital and carbonate mass accumulation rates through the PETM interval.

During the deposition of the characteristic PETM SU of Zumaia, sedimentation rates double and the mass accumulation rates of detrital sediments increase by up to 4 times relative to pre-event conditions. Dramatic increases in sedimentation rates during the main phase of the PETM are interpreted to represent significant changes in regional hydrology, most likely an increase in the magnitude and frequency of extreme rainfall and run-off events, potentially within an overall lower mean annual precipitation regime (Carmichael et al., 2016, 2017). These changes in sediment accumulation rates – and inferred perturbations to hydrology – persist throughout the PETM interval, only returning towards pre-event levels during the rapid phase of CIE recovery. This supports a close coupling between the CIE and the climate system. Peak sedimentation rates and detrital accumulation rates are observed through the recovery phase of the CIE and not the onset. We suggest that this peak in sedimentation rates is coupled to the precession forcing of subtropical terrestrial hydrology and may provide a key to understanding the timing of the recovery phase of the CIE.

Earth system model (cGENIE) inversion experiments, designed to replicate the observed CIE dynamics at Zumaia, reproduce an atmospheric CIE that is close to the best available terrestrial records of the CIE (Bowen et al., 2015), in timing and magnitude. Depending on the δ¹³C composition of the source of PETM carbon, estimated total carbon inputs during the onset of the PETM are ∼ 3400 Pg C (methane source, −60 ‰) or ∼ 11 300 Pg C (organic carbon source, −22 ‰). Simulations with the silicate weathering feedback as the only mechanism for carbon cycle recovery fail to reproduce the observed record in two key ways: first, they are unable to recover the coupled CIE–climate system at a rate fast enough to match the observations (< 20 kyr), and, second, the final sink of carbon within this feedback – marine carbonates – is isotopically too heavy to provide any significant recovery of the CIE. Simulations that match the patterns of CIE recovery require active organic or methane carbon removal of a magnitude similar to the original carbon input. Simulations with an isotopic mismatch between source carbon (methane −60 ‰) and sink carbon (organic matter −22 ‰) require significantly more carbon removal at the end of the event than originally released during the onset. This causes a marked post-event cooling of ∼ 4 ^∘C relative to pre-event conditions that is inconsistent with current climate records.

The results presented here demonstrate that the Earth system was capable of the sequestration of several thousand petagrams of carbon at the end of the PETM, but only after nearly 100 kyr within a profoundly altered climate state. The mechanism of carbon burial also appears to rely on a radical perturbation to global hydro-climates, sufficient to cause erosion and dramatically increased sediment flux rates to continental margins between 3 and 4 times pre-perturbation levels.