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 – past climate transients with global-scale warming signiﬁcantly above background climate variability Here we present new high-resolution bulk sediment stable and major element data for the classic PETM With data we provide a

tive 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 (δ 13 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 2 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 Zeebe et al., 2009), but this delay only exacerbates the mis-match 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.

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).
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(Alegret et al., , 2009Baceta et al., 2000;Bernaola et al., 2007;Canudo et al., 1995;Dinarès-Turell et al., 2007, 2002Manners et al., 2013;Schmitz et al., 1997Schmitz et al., , 2001Storme 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 greenishgrey marl unit, which is distinct from both the limestonemarl 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 3 using a calibration set of 25 discrete calcium carbonate measurements. Weight percent CaCO 3 contents were determined at UoB, based on inorganic carbon contents measured by CO 2 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 (δ 13 C carb and δ 18 O carb ) of bulk carbonates at the NERC Isotope Geosciences Laboratory (NIGL, Keyworth, Nottingham). Sample material was powdered to yield CO 2 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 2 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.

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 www.clim-past.net/14/1035/2018/ Clim. Past, 14, 1035-1049, 2018 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 δ 13 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 δ 13 C values concerned are all within the background variations in δ 13 C seen in the marl-limestone couplets of the latest Paleocene (e.g. from −4.0 to −1.0 m). For example, the carbonate δ 13 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 δ 13 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 3 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 δ 13 C shift between 0.25 m (0.01 ‰) and 0.30 m (−0.28 ‰) is the start of the consistent negative deviation in δ 13 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 δ 13 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(Alegret et al., , 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).
Although an increase in sedimentation rates at the start of the PETM in the Zumaia section limits the effects of any car-bonate burn-down on the earliest stages of the CIE profile, the measured δ 13 C carb data diverge from the typical PETM CIE profile between ∼ 0.7 and 2.8 m. Here, δ 13 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 δ 13 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 deepocean 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 δ 13 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).

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 2 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 highresolution 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 δ 13 C inflection point F of Röhl et al. (2007;Fig . 3). With a total CIE duration of 90-167 kyr Röhl et al., 2007) this variability is on the half-precession timescale (10 · 10.5 = 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).
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 www.clim-past.net/14/1035/2018/ Clim. Past, 14, 1035-1049, 2018 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., 2007Westerhold et al., , 2008Westerhold et al., , 2014Westerhold et al., , 2018Zachos 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 dom-inated 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.   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 3 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).
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 ra-tios. 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 δ 13 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 deepwater 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 .

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 Paleoceneearly Eocene configuration of the cGENIE Earth system model (Kirtland . 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 2 set to 834 ppm with a δ 13 C of −4.9 ‰. A further 250 kyr "open" system spin-up with temperature-dependent weathering feedbacks enabled and input of CO 2 via volcanic outgassing (Colbourn et al., 2013) was used to bring the geological 13 C 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 (  vation of 46.75 %, comparable with late Paleocene estimates (Panchuk et al., 2008). Four inversion experiments were undertaken, all forced with carbon (CO 2 ) inputs to the atmosphere such that the δ 13 C of surface ocean dissolved inorganic carbon (DIC) followed a prescribed inversion target taken from the Zumaia δ 13 C record, adjusted to the initial cGENIE average surface ocean DIC δ 13 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 δ 13 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 δ 13 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 −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 δ 13 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 δ 13 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 δ 13 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 2 . 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. 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 2 ? 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 2 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 δ 13 C CO 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 . 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 cli-mate 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. δ 13 C and CaCO 3 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).
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 2 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 (1) carbon input at δ 13 C of −22 ‰ during the onset, no carbon removal (green); (2) carbon input and removal with δ 13 C of −22 ‰ to follow Zumaia target CIE (red); (3) carbon input with δ 13 C of −60 ‰ and removal of carbon with δ 13 C of −22 ‰ to follow Zumaia target CIE (blue). All simulations have active temperature-dependent silicate weathering feedbacks operative throughout.
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 oceanatmosphere system in the −22 ‰ δ 13 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 −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 carbononly 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 2 on the balance between rates of net primary production and organic carbon remineralization during the PETM (Cot-  , 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.

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 highresolution 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 in-terpreted 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(Carmichael et al., , 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 δ 13 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.