Journal topic
Clim. Past, 14, 1035–1049, 2018
https://doi.org/10.5194/cp-14-1035-2018
Clim. Past, 14, 1035–1049, 2018
https://doi.org/10.5194/cp-14-1035-2018

Research article 11 Jul 2018

Research article | 11 Jul 2018

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

Dynamics of sediment flux to a bathyal continental margin section through the Paleocene–Eocene Thermal Maximum
Tom Dunkley Jones1, Hayley R. Manners2,3, Murray Hoggett1, Sandra Kirtland Turner4, Thomas Westerhold5, Melanie J. Leng6, Richard D. Pancost7, Andy Ridgwell4, Laia Alegret8, Rob Duller9, and Stephen T. Grimes2 Tom Dunkley Jones et al.
• 1School of Geography, Earth and Environmental Sciences, University of Birmingham, Edgbaston, Birmingham, B15 2TT, UK
• 2School of Geography, Earth & Environmental Sciences, Plymouth University, Drake Circus, Plymouth, Devon, PL4 8AA, UK
• 3School of Ocean and Earth Sciences, National Oceanography Centre, University of Southampton, Southampton, SO14 3ZH, UK
• 4Department of Earth Sciences, University of California, Riverside, CA 92521, USA
• 5MARUM – Center for Marine Environmental Sciences, University of Bremen, Leobener Strasse, 28359 Bremen, Germany
• 6NERC Isotope Geosciences Facilities, British Geological Survey, Nottingham NG12 5GG, UK & Centre for Environmental Geochemistry, University of Nottingham, Nottingham, NG7 2RD, UK
• 7Organic Geochemistry Unit, The Cabot Institute, School of Chemistry, University of Bristol, Bristol, BS8 1TS, UK
• 8Departamento de Ciencias de la Tierra & Instituto Universitario de Ciencias Ambientales, Universidad de Zaragoza, 50009 Zaragoza, Spain
• 9Department of Earth, Ocean and Ecological Sciences, University of Liverpool, Liverpool, L69 3GP, UK

Correspondence: Tom Dunkley Jones (t.dunkleyjones@bham.ac.uk)

Abstract

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 (δ13C) 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 CO2 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 Corg 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; 043184.5′′ N, 0021531.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 1Location 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 CaCO3 using a calibration set of 25 discrete calcium carbonate measurements. Weight percent CaCO3 contents were determined at UoB, based on inorganic carbon contents measured by CO2 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−3. Organic carbon MARs were calculated following Willmott et al. (2009).

A subset of 129 thumbnail samples were also analysed for the stable isotope composition (δ13Ccarb and δ18Ocarb) of bulk carbonates at the NERC Isotope Geosciences Laboratory (NIGL, Keyworth, Nottingham). Sample material was powdered to yield CO2 from the equivalent of 10 mg of 100 % calcite and reacted with anhydrous phosphoric acid in vacuo overnight at a constant 25 C. Liberated CO2 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 δ13C 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 δ13C values concerned are all within the background variations in δ13C seen in the marl–limestone couplets of the latest Paleocene (e.g. from −4.0 to −1.0 m). For example, the carbonate δ13C 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 δ13C 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 CaCO3 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 δ13C shift between 0.25 m (0.01 ‰) and 0.30 m (−0.28 ‰) is the start of the consistent negative deviation in δ13C 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 δ13C – 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 2New 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) SiO2 weight percent (green) and Si  Fe ratio (cyan) based on weight percentages of oxides (SiO2Fe2O3) from bulk sediment XRF analyses; pre-PETM and post-PETM precession cycles numbered in black; (d) weight percentage of CaO and CaCO3 from XRF; (e) new bulk carbonate δ13C (this study). Panels (f) to (l) show data against age (Westerhold et al., 2018): (f) Polecat Bench core (PCB) soil carbonate nodule δ13C (green; Bown et al. 2015); (g) deep-ocean benthic foraminiferal δ13C (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 δ13Ccarb data diverge from the typical PETM CIE profile between  0.7 and 2.8 m. Here, δ13Ccarb 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 δ13Ccarb 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 δ13C 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 SiO2 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−1 ( 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 δ13C 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 3New sediment geochemistry data for the Zumaia section: (a) lithological log produced during sample collection for this study (2011); (b) bulk carbonate δ13C; (c) weight percentages of CaO and CaCO3 from bulk sediment XRF analyses; (d) Si  Fe ratio based on weight percentages of oxides (SiO2Fe2O3) 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 4Multitaper 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 CaCO3, 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 CaCO3, 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 (CaCO3 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 5New 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 δ13C; (c) weight percentage of CaCO3 derived from XRF CaCO calibrated to discrete CaCO3 analyses (see methods); (d) CaCO3 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−2 ka−1, 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 δ13C 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−2 ka−1. 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−2 ka−1. Coincident with the CIE recovery phase, between 55.83 and 55.81 Ma, TOC MARs again increase to over  0.02 g cm−2 ka−1.

Figure 6Comparison 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 δ13C for deep-ocean sites; (b) bulk carbonate δ13C for Zumaia; (c) soil nodule δ13C from the Big Horn Basin; (d) CaCO3 records from Zumaia and the deep-ocean sites; (e) detrital (red) and CaCO3 (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 CO2 set to 834 ppm with a δ13C of −4.9 ‰. A further 250 kyr “open” system spin-up with temperature-dependent weathering feedbacks enabled and input of CO2 via volcanic outgassing (Colbourn et al., 2013) was used to bring the geological 13C cycle to steady state. In this, the specified non-carbonate detrital flux to marine sediments (0.18 g cm−2 ka−1), the rain ratio of carbonate to particulate organic carbon (0.2), and the global average concentrations of Mg2+ and Ca2+ (29.89 and 18.22 mmol kg−1, respectively) and alkalinity (1975 µ mol kg−1) 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 (CO2) inputs to the atmosphere such that the δ13C of surface ocean dissolved inorganic carbon (DIC) followed a prescribed inversion target taken from the Zumaia δ13C record, adjusted to the initial cGENIE average surface ocean DIC δ13C 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 δ13C 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 δ13C 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−1) 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 δ13C 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 δ13C 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 δ13C 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 CO2. 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 1Summary 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 CO2? 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 pCO2 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 δ13C${}_{{\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. δ13C and CaCO3 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 7Results of inverse modelling with the cGENIE Earth system model: (a) the bulk carbonate δ13C for Zumaia which drives model inversions; (b) soil nodule δ13C record of the Big Horn Basin, as representative of atmospheric CO2 δ13C for comparison to model results; (c) modelled δ13C composition of atmospheric carbon dioxide; (d) modelled and actual weight percentage of CaCO3 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 pCO2. The model scenarios are as follows: (1) carbon input at δ13C of −22 ‰ during the onset, no carbon removal (green); (2) carbon input and removal with δ13C of −22 ‰ to follow Zumaia target CIE (red); (3) carbon input with δ13C of −60 ‰ and removal of carbon with δ13C 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 CO2 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 ‰ δ13C 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−2 ka−1) 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 CO2 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 δ13C 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.

Data availability
Data availability.

Correspondence and requests for materials should be addressed to Tom Dunkley Jones.

Supplement
Supplement.

Author contributions
Author contributions.

TDJ and STG conceived and co-ordinated this study; TDJ, HRM and LA logged and sampled the section; MH ran XRF and carbonate content analyses; SKT and AR ran and interpreted GENIE model simulations; TW undertook spectral analyses of the data and correlations with other PETM sections; MJL ran stable isotope analyses; HRM ran organic carbon content analyses; RDP advised on the interpretation of isotope and organic carbon data; RD provided assistance with the regional sedimentology; TDJ wrote the manuscript with the assistance of all co-authors.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

Tom Dunkley Jones acknowledges funding from a Royal Society Dorothy Hodgkin Fellowship, Laia Alegret project CGL2014-58794-P (Spanish Ministry of Economy and Competitiveness and FEDER funds), and Richard D. Pancost acknowledges the Royal Society Wolfson Research Merit Award. This research was supported by the Deutsche Forschungsgemeinschaft (DFG) to Thomas Westerhold. We would also like to thank Victoriano Pujalte for guidance in accessing sampling sites, Enrique Cantero and Georgia Hole for assistance in sampling, and Chris Kendrick for assistance with TOC analysis. Finally, we appreciate the discussion and the helpful comments and suggestions, which improved this study, provided by Victoriano Pujalte, Michael Clare and the anonymous referee.

Edited by: Laurie Menviel
Reviewed by: Victoriano Pujalte, Michael Clare, and one anonymous referee

References

Abdul Aziz, H., Hilgen, F. J., Van Luijk, G. M., Sluijs, A., Kraus, M. J., Pares, J. M., and Gingerich, P. D.: Astronomical climate control on paleosol stacking patterns in the upper Paleocene-lower Eocene Willwood Formation, Bighorn Basin, Wyoming, Geology, 36, 531–534, 2008.

Alegret, L., Ortiz, S., Orue-Etxebarria, X., Bernaola, G., Baceta, J. I., Monechi, S., Apellaniz, E., and Pujalte, V.: The Paleocene-Eocene Thermal Maximum: New data on microfossil turnover at the Zumaia Section, Spain, Palaios, 24, 318–328, 2009.

Alegret, L., Matías, R., and Pérez Manuel, V.: Environmental instability during the latest Paleocene at Zumaia (Basque-Cantabric Basin): The bellwether of the Paleocene-Eocene Thermal Maximum, Palaeogeogr. Palaeocl., 497, 186–200, 2018.

Baceta, J. I., Pujalte, V., Dinarès-Turell, J., Payros, A., Orue-Etxebarria, X., and Bernaola, G.: The Paleocene/Eocene boundary interval in teh Zumaia Section (Gipuzkoa, Basque Basin): Magnetostratigraphy, and high-resolution lithostratigraphy, Revista de la Sociedad Geológica de España, 13, 375–391, 2000.

Baceta, J. I., Pujalte, V., and Caballero, F.: Paleocene and Early Eocene Facies and Events: a Basin-Platform-Coastal Plain Transect (South-Central and Western Pyrenees), Climate and Biota of the Early Paleogene 2006, Post Conference Field Trip Guidebook, Bilbao, 2006.

Bains, S., Corfield, R., and Norris, R.: Mechanisms of climate warming at the end of the Paleocene, Science, 285, 724–727, 1999.

Bernaola, G., Baceta, J. I., Orue-Etxebarria, X., Alegret, L., Martin-Rubio, M., Arostegui, J., and Dinares-Turell, J.: Evidence of an abrupt environmental disruption during the mid-Paleocene biotic event (Zumaia section, western Pyrenees), Geol. Soc. Am. Bull., 119, 785–795, 2007.

Bowen, G. J. and Zachos, J. C.: Rapid carbon sequestration at the termination of the Palaeocene-Eocene Thermal Maximum, Nat. Geosci., 3, 866–869, 2010.

Bowen, G. J., Maibauer, B. J., Kraus, M. J., Rohl, U., Westerhold, T., Steimke, A., Gingerich, P. D., Wing, S. L., and Clyde, W. C.: Two massive, rapid releases of carbon during le onset of the Palaeocene-Eocene thermal maximum, Nat. Geosci., 8, 44–47, 2015.

Canudo, J. I., Keller, G., Molina, E., and Oritz, N.: Planktic foraminiferal turnover and d13C isotopes across the Paleocene-Eocene transition at Caravaca and Zumaya, Spain, Palaeogeogr. Palaeocl., 114, 75–100, 1995.

Carmichael, M. J., Lunt, D. J., Huber, M., Heinemann, M., Kiehl, J., LeGrande, A., Loptson, C. A., Roberts, C. D., Sagoo, N., Shields, C., Valdes, P. J., Winguth, A., Winguth, C., and Pancost, R. D.: A model-model and data-model comparison for the early Eocene hydrological cycle, Clim. Past, 12, 455–481, https://doi.org/10.5194/cp-12-455-2016, 2016.

Carmichael, M. J., Inglis, G. N., Badger, M. P. S., Naafs, B. D. A., Behrooz, L., Remmelzwaal, S., Monteiro, F. M., Rohrssen, M., Farnsworth, A., Buss, H. L., Dickson, A. J., Valdes, P. J., Lunt, D. J., and Pancost, R. D.: Hydrological and associated biogeochemical consequences of rapid global warming during the Paleocene-Eocene Thermal Maximum, Global Planet. Change, 157, 114–138, 2017.

Clare, M. A., Talling, P. J., and Hunt, J. E.: Implications of reduced turbidity current and landslide activity for the Initial Eocene Thermal Maximum - evidence from two distal, deep-water sites: Earth Planet. Sc. Lett., 420, 102–115, 2015.

Colbourn, G., Ridgwell, A., and Lenton, T. M.: The Rock Geochemical Model (RokGeM) v0.9, Geosci. Model Dev., 6, 1543–1573, https://doi.org/10.5194/gmd-6-1543-2013, 2013.

Cotton, J. M., Sheldon, N. D., Hren, M. T., and Gallagher, T. M.: Positive feedback drives carbon release from soils to atmosphere during Paleocene/Eocene warming, Am. J. Sci., 315, 337–361, 2015.

Cui, Y., Kump, L. R., Ridgwell, A. J., Charles, A. J., Junium, C. K., Diefendorf, A. F., Freeman, K. H., Urban, N. M., and Harding, I. C.: Slow release of fossil carbon during the Palaeocene-Eocene Thermal Maximum, Nat. Geosci., 4, 481–485, 2011.

Dickens, G. R., O'Neil, J. R., Rea, D. K., and Owen, R. M.: Dissociation of oceanic methane hydrate as a cause of the carbon isotope excursion at the end of the Paleocene, Paleoceanography, 10, 965–971, 1995.

Dinarès-Turell, J., Baceta, J. I., Pujalte, V., Orue-Etxebarria, X., and Bernaola, G.: Magnetostratigraphic and cyclostratigraphic calibration of a prospective Palaeocene/Eocene stratotype at Zumaia (Basque Basin, northern Spain), Terra Nova, 14, 371–378, 2002.

Dinarès-Turell, J., Baceta, J. I., Pujalte, V., Orue-Etxebarria, X., Bernaola, G., and Lorito, S.: Untangling the Palaeocene climatic rhythm: an astronimically calibrated Early Palaeocene magnetostratigraphy and biostratigraphy at Zumaia (Basque basin, northern Spain), Earth Planet. Sc. Lett., 216, 483–500, 2003.

Dinarès-Turell, J., Baceta, J., Bernaola, G., Orue-Etxebarria, X., and Pujalte, V.: Closing the Mid-Palaeocene gap: Toward a complete astronomically tuned Palaeocene Epoch and Selandian and Thanetian GSSPs at Zumaia (Basque Basin, W Pyrenees), Earth Planet. Sc. Lett., 262, 450–467, 2007.

Dunkley Jones, T., Lunt, D. J., Schmidt, D. N., Ridgwell, A., Sluijs, A., Valdes, P. J., and Maslin, M.: Climate model and proxy data constraints on ocean warming across the Paleocene/Eocene Thermal Maximum, Earth-Sci. Rev., 125, 123–145, 2013.

Gavrilov, Y. O., Kodina, L., Lubchenko, I., and Muzylev, N.: The late Paleocene anoxic event in epicontinental seas of the Peri-Tethys and formation of the sapropelite unit, Sedimentology and geochemistry, Lithol. Miner. Resour., 32, 427–450, 1997.

Giusberti, L., Boscolo Galazzo, F., and Thomas, E.: Variability in climate and productivity during the Paleocene-Eocene Thermal Maximum in the western Tethys (Forada section), Clim. Past, 12, 213–240, https://doi.org/10.5194/cp-12-213-2016, 2016.

Gutjahr, M., Ridgwell, A., Sexton, P. F., Anagnostou, E., Pearson, P. N., Pälike, H., Norris, R. D., Thomas, E., and Foster, G. L.: Very large release of mostly volcanic carbon during the Palaeocene-Eocene Thermal Maximum, Nature, 548, 573–577, 2017.

John, C. M., Bohaty, S. M., Zachos, J. C., Sluijs, A., Gibbs, S., Brinkhuis, H., and Bralower, T. J.: North American continental margin records of the Paleocene-Eocene thermal maximum, Implications for global carbon and hydrological cycling, Paleoceanography, 23, PA2217, https://doi.org/10.1029/2007PA001465, 2008.

Kirtland Turner, S. and Ridgwell, A.: Development of a novel empirical framework for interpreting geological carbon isotope excursions, with implications for the rate of carbon injection across the PETM: Earth Planet. Sc. Lett., 435, 1–13, 2016.

Manners, H. R., Grimes, S., Sutton, P. A., Domingo, L., Leng, M. J., Twitchett, R. T., Hart, M. B., Dunkley Jones, T., Pancost, R. D., Duller, R., and Lopez-Martinez, N.: Magnitude and profile of organic carbon isotope records from the Palaeocene-Eocene Thermal Maximum: evidence from northern Spain, Earth Planet. Sc. Lett., 377, 220–230, 2013.

McCarren, H., Thomas, E., Hasegawa, T., Röhl, U., and Zachos, J. C.: Depth dependency of the Paleocene-Eocene carbon isotope excursion: Paired benthic and terrestrial biomarker records (Ocean Drilling Program Leg 208, Walvis Ridge), Geochem. Geophy. Geosy., 9, Q10008, https://doi.org/10.1029/2008GC002116, 2008.

McInerney, F. A. and Wing, S. L.: The Palaeocene-Eocene Thermal Maximum: A perturbation of carbon cycle, climate, and biosphere with implications for the future, Annu. Rev. Earth Pl. Sc., 39, 489–516, 2011.

Meissner, K. J., Bralower, T., Alexander, K., Dunkley Jones, T., Sijp, W., and Ward, M.: The Paleocene-Eocene Thermal Maximum: How much carbon is enough?, Paleoceanography, 29, 946–963, 2014.

Meyers, S. R.: Astrochron: An R Package for Astrochronology, available at: https://CRAN.R-project.org/package=astrochron (last access: 28 June 2018), 2014.

Murphy, B. H., Farley, K. A., and Zachos, J. C.: An extraterrestrial 3He-based timescale for the Paleocene-Eocene thermal maximum (PETM) from the Walvis Ridge, IODP Site 1266, Geochem. Cosmochem. Acta, 74, 5098–5108, 2010.

Panchuk, K., Ridgwell, A., and Kump, L.: Sedimentary response to Palaeocene-Eocene Thermal Maximum carbon release: A model-data comparison, Geology, 36, 315–318, 2008.

Penman, D. E.: Silicate weathering and North Atlantic silica burial during the Paleocene-Eocene Thermal Maximum, Geology, 44, 731–734, 2016.

Penman, D. E., Turner, S. K., Sexton, P. F., Norris, R. D., Dickson, A. J., Boulila, S., Ridgwell, A., Zeebe, R. E., Zachos, J. C., Cameron, A., Westerhold, T., and Röhl, U.: An abyssal carbonate compensation depth overshoot in the aftermath of the Palaeocene-Eocene Thermal Maximum, Nat. Geosci., 9, 575–580, 2016.

Pujalte, V., Baceta, J. I., and Schmitz, B.: A massive input of coarse-grained siliciclastics in the Pyrenean Basin during the PETM: the missing ingredient in a coeval abrupt change in hydrological regime, Clim. Past, 11, 1653–1672, https://doi.org/10.5194/cp-11-1653-2015, 2015.

Ridgwell, A. and Schmidt, D. N.: Past constraints on the vulnerability of marine calcifiers to massive carbon dioxide release, Nat. Geosci., 3, 1–5, 2010.

Röhl, U., Westerhold, T., Bralower, T., and Zachos, J.: On the duration of the Palaeocene-Eocene thermal maximum (PETM), Geochem. Geophy. Geosy., 8, 1–13, 2007.

Schmitz, B. and Pujalte, V.: Abrupt increase in seasonal extreme precipitation at the Palaeocene-Eocene boundary, Geology, 35, 215–218, 2007.

Schmitz, B., Asaro, F., Molina, E., Monechi, S., Salis, K. V., and Speijer, R. P.: High-resolution iridium, d13C, d18O, foraminifera and nannofossil profiles across the latest Paleocene benthic extinction event at Zumaya, Palaeogeogr. Palaeocl., 133, 49–68, 1997.

Schmitz, B., Pujalte, V., and Nunez-Betelu, K.: Climate and sea-level perturbations during the Initial Eocene Thermal Maximum: evidence from siliciclastic units in the Basque Basin (Ermua, Zumaia and Trabakua Pass), northern Spain, Palaeogeogr. Palaeocl., 165, 299–320, 2001.

Short, D. A., Mengel, J. G., Crowley, T. J., Hyde, W., and North, G. R.: Filtering of Milankovitch Cycles by Earth's Geography, Quaternary Res., 35, 157–173, 1991.

Sluijs, A., Bowen, G. J., Brinkhuis, H., Lourens, L. J., and Thomas, E.: The Palaeocene-Eocene Thermal Maximum super greenhouse: biotic and geochemical signatures, age models and mechanisms of global change, in: Deep-time perspectives on climate change: Marrying the signal from computer models and biological proxies, edited by: Williams, M., Haywood, A. M., Gregory, J., and Schmidt, D. N., The Micropalaeontological Society Special Publication, London, 323–350, 2007.

Sluijs, A., van Roij, L., Harrington, G. J., Schouten, S., Sessa, J. A., LeVay, L. J., Reichart, G.-J., and Slomp, C. P.: Warming, euxinia and sea level rise during the Paleocene-Eocene Thermal Maximum on the Gulf Coastal Plain: implications for ocean oxygenation and nutrient cycling, Clim. Past, 10, 1421–1439, https://doi.org/10.5194/cp-10-1421-2014, 2014.

Storme, J.-Y., Devleeschouwer, X., Schnyder, J., Cambier, G., Baceta, J. I., Pujalte, V., Di Matteo, A., Iacumin, P., and Yans, J.: The Palaeocene/Eocene boundary section at Zumaia (Basque-Cantabric Basin) revisited: new insights from high-resolution magnetic susceptibility and carbon isotope chemostratigraphy on organic matter (δ13Corg), Terra Nova, 24, 310–317, 2012.

Svensen, H., Planke, S., Malthe-Sorenssen, A., Jamtveit, B., Myklebust, R., Eidem, T., and Rey, S.: Release of methane from a volcanic basin as a mechanism for initial Eocene global warming, Nature, 429, 542–545, 2004.

Westerhold, T., Röhl, U., Laskar, J., Raffi, I., Bowles, J., Lourens, L., and Zachos, J.: On the duration of magnetochrons C24r and C25n and the timing of early Eocene global warming events: Implications from the Ocean Drilling Program Leg 208 Walvis Ridge depth transect, Paleoceanography, 22, PA2201, https://doi.org/10.1029/2006PA001322, 2007.

Westerhold, T., Röhl, U., Raffi, I., Fornaciari, E., Monechi, S., Reale, V., Bowles, J., and Evans, H. F.:, Astronomical calibration of the Paleocene time, Palaeogeogr. Palaeocl., 257, 377–403, 2008.

Westerhold, T., Röhl, U., Pälike, H., Wilkens, R., Wilson, P. A., and Acton, G.: Orbitally tuned timescale and astronomical forcing in the middle Eocene to early Oligocene, Clim. Past, 10, 955–973, https://doi.org/10.5194/cp-10-955-2014, 2014.

Westerhold, T., Röhl, U., Wilkens, R. H., Gingerich, P. D., Clyde, W. C., Wing, S. L., Bowen, G. J., and Kraus, M. J.: Synchronizing early Eocene deep-sea and continental records – cyclostratigraphic age models for the Bighorn Basin Coring Project drill cores, Clim. Past, 14, 303–319, https://doi.org/10.5194/cp-14-303-2018, 2018.

Willmott, V., Rampen, S., Domack, E., Canals, M., Sinninghe Damste, J. S., and Schouten, S.: Holocene changes in Proboscia diatom productivity in shelf waters of the north-western Antarctic Peninsula, Antarct. Sci., 22, 3–10, 2009.

Zachos, J. C., Röhl, U., Schellenberg, S. A., Sluijs, A., Hodell, D. A., Kelly, D. C., Thomas, E., Nicolo, M., Raffi, I., Lourens, L. J., McCarren, H., and Kroon, D.: Rapid acidification of the ocean during the Paleocene-Eocene thermal maximum, Science, 308, 1611–1615, 2005.

Zachos, J. C., Mccarren, H., Murphy, B., Röhl, U., and Westerhold, T.: Tempo and scale of late Paleocene and early Eocene carbon isotope cycles: Implications for the origin of hyperthermals, Earth Planet. Sc. Lett., 299, 242–249, 2010.

Zeebe, R. E., Zachos, J. C., and Dickens, G. R.: Carbon dioxide forcing alone insufficient to explain Palaeocene-Eocene Thermal Maximum warming, Nat. Geosci., 2, 576–580, 2009.

Zeebe, R. E., Ridgwell, A., and Zachos, J. C.: Anthropogenic carbon release rate unprecedented during the past 66 million years, Nat. Geosci., 9, 325–329, 2016.