Alluvial record of an early Eocene hyperthermal within the Castissent Formation, the Pyrenees, Spain

The late Palaeocene to the middle Eocene (57.5 to 46.5 Ma) recorded a total of 39 hyperthermals – periods of rapid global warming documented by prominent negative carbon isotope excursions (CIEs) as well as peaks in iron content – have been recognized in marine cores. Documenting how the Earth system responded to rapid climatic shifts during hyperthermals provides fundamental information to constrain climatic models. However, while hyperthermals have been well documented in the marine sedimentary record, only a few have been recognized and described in continental deposits, thereby limiting our ability to understand the effect and record of global warming on terrestrial systems. Hyperthermals in the continental record could be a powerful correlation tool to help connect marine and continental deposits, addressing issues of environmental signal propagation from land to sea. In this study, we generate new stable carbon isotope data (δ¹³C values) across the well-exposed and time-constrained fluvial sedimentary succession of the early Eocene Castissent Formation in the south central Pyrenees (Spain). The δ¹³C values of pedogenic carbonate reveal – similarly to the global records – stepped CIEs, culminating in a minimum δ¹³C value that we correlate with the hyperthermal event “U” at ca. 50 Ma. This general trend towards more negative values is most probably linked to higher primary productivity leading to an overall higher respiration of soil organic matter during these climatic events. The relative enrichment in immobile elements (Zr, Ti, Al) and higher estimates of mean annual precipitation together with the occurrence of small iron oxide and iron hydroxide nodules during the CIEs suggest intensification of chemical weathering and/or longer exposure of soils in a highly seasonal climate. The results show that even relatively small-scale hyperthermals compared with their prominent counterparts, such as PETM, ETM2, and ETM3, can leave a recognizable signature in the terrestrial stratigraphic record, providing insights into the dynamics of the carbon cycle in continental environments during these events.

1 Introduction

From the end of the Palaeocene, a period of global warming reached its climax during the Early Eocene Climatic Optimum (EECO) (Westerhold and Röhl, 2009; Hyland and Sheldon, 2013). The EECO started ca. 53 Ma ago and lasted until ca. 49 Ma ago (Westerhold et al., 2018), after which the climate began to cool (∼ Eocene–Oligocene transition; Zachos et al., 2001, 2008). Superimposed on, and coeval to, this globally warm epoch, brief periods of pronounced global warming known as “hyperthermals” standout as anomalies outside of background climate variability (Kirtland-Turner et al., 2014; Dunkley Jones et al., 2018). The Palaeocene–Eocene Thermal Maximum (PETM; ∼56 Ma) was the first of these events to be identified globally because of its exceptional magnitude and preservation in both marine and continental deposits (Koch et al., 1992). To date, for the late Palaeocene – early Eocene period, a total of 39 hyperthermal events of lesser magnitude have been identified from marine cores (Lourens et al., 2005; Sexton et al., 2011; Kirtland-Turner et al., 2014; Lauretano et al., 2015, 2016; and Westerhold et al., 2018), among which the most prominent and studied are the Eocene Thermal Maximum (ETM) 2 and 3, H2, and I1 and I2 events (Cramer et al., 2003; Lourens et al., 2005; Nicolo et al., 2007; Lunt et al., 2011; Deconto et al., 2012; Kirtland-Turner et al., 2014; Lauretano et al., 2016; and Westerhold et al., 2017) (Fig. 1). In the marine stratigraphic record, these events are primarily characterized by paired negative excursion in carbon and oxygen isotope data exceeding background variability (Cramer et al., 2003; Nicolo et al., 2007; Zachos et al., 2008; Sluijs and Dickens, 2012; and Lauretano et al., 2016), i.e. typically with amplitude greater than the standard deviation (SD) of pre-hyperthermal background values.

Figure 1 Late Palaeocene and early Eocene benthic carbon isotope record from Sites 1209, 1258, 1262, and 1263. Top of Chron C22r and top of T. orthostylus zone from site 1263 from Westerhold et al. (2017). Hyperthermal nomenclature from Cramer et al. (2003), Lauretano et al. (2016), and Westerhold et al. (2017). Castissent Fm. extension in green.

In deep marine settings, the carbon isotope excursions (CIE) are typically paired with an increase in iron concentration and decrease in carbonate content, indicating ocean acidification potentially linked with high atmospheric CO₂ concentrations (Nicolo et al., 2007; Slotnick et al., 2012; and Westerhold et al., 2018). In coastal marine sections, early Eocene hyperthermal events are generally associated with an enhanced flux of terrigenous material, interpreted as linked to accelerated hydrological cycle and higher seasonality (Schmitz et al., 2001; Bowen et al., 2004; Nicolo et al., 2007; Slotnick et al., 2012; Payros et al., 2015; and Dunkley Jones et al., 2018), although several studies document a spatially heterogeneous hydrological climatic response during the PETM (Bolle and Adatte, 2001; Kraus and Riggins, 2007; Giusberti et al., 2016; and Carmichael et al., 2017). In fluvial systems, the abrupt warming of the PETM was found to be associated with expansion and coarsening of alluvial facies combined with an increase in the magnitude of flood discharge (Foreman et al., 2012; Pujalte et al., 2015; and Chen et al., 2018), as well as enhanced pedogenesis (Abels et al., 2012). Yet, how continental systems reacted to the other, smaller-magnitude hyperthermals of the early Eocene remains to be documented. In particular, because of the subaerial nature and lateral preservation dynamics of alluvial systems (e.g. Foreman and Straub, 2017; Straub and Foreman, 2018), the extent to which fluvial successions can provide complete and faithful archives of past climatic events, especially those with the smallest magnitudes, is still largely unknown (Foreman and Straub, 2017; Trampush et al., 2017; and Straub and Foreman, 2018). Addressing this question is particularly critical for studies focussing on environmental signal propagation in source-to-sink systems (e.g. Castelltort and Van Den Driessche, 2003; Duller et al., 2019; Romans et al., 2016; and Schlunegger and Castelltort, 2016), which require high-resolution continental-marine correlations such as those provided by the PETM (e.g. Duller et al., 2019) or by other hyperthermals of the early Eocene.

Figure 2 Simplified situation and geological map of the study area with main depositional palaeo-environments (e.g. Nijman, 1998). The Castissent Fm. is a prominent fluvial unit particularly well exposed in the Noguera Ribagorçana and Isáabena river valleys. (1) Chiriveta section (2) Mas de Faro (3) La Roca section. Main palaeoflow directions indicated in orange (from Nijman and Puigdefabregas, 1978). Regional map after Teixell (1998).

To address these issues, we explored the geochemical signature (carbon and oxygen stable isotopes; major and trace elements) and the sedimentology of the fluvial deposits of the Ypresian Age Castissent Fm. (south central Pyrenees, Spain, Fig. 2). First, we generated a new carbon isotope profile from a palaeosol succession rich in carbonate nodules across the Castissent Fm. in order to compare these results with a global δ¹³C record. The data suggest that this fluvial succession preserves a record of hyperthermal “U” event at ca. 50 Ma, providing important constraints to its depositional the age. Second, we used the major and trace element composition of bulk floodplain material in order to explore the climatic impact of such a hyperthermal, including empirical reconstructions of mean annual precipitation, allowing us to discuss soil dynamics during global warming. This study identifies for the first time in a continental succession an event so far only recorded in marine sediments, thereby demonstrating the global breadth of these climatic events and the complementarity of oceanographic and terrestrial archives.

Figure 3 (a) Time constraints on the Castissent Fm. MP zone from the continental section from Checa Soler (2004) and Payros et al. (2009). SBZ and NP in the Campo section from (Schaub, 1966, 1981; Kapellos and Schaub, 1973; and Tosquella, 1995), magnetostratigraphy from Bentham and Burbank (1996). SBZ in El Pueyo section from Payros et al. (2009). Magnetostratigraphy in El Pueyo from Poyatos-Moré (2014). (b) Extended map of the study area. For map legend and references, see Fig. 2.

2 Geological setting

The Castissent Formation comprises fluvial deposit of Ypresian age cropping out in the Tremp-Graus Basin (South Pyrenean foreland basin; Marzo et al., 1988, Fig. 2). The Castissent Fm. is defined by its prominent overall sand-rich character and is composed by three coarse-grained channel complexes (labelled as Members A, B, and C) separated by four marine incursions (M0 to M3) inferred from the observation of marginal coastal bioclast-rich horizons developed up into the upper deltaic plain and correlative with finer dark-grey mudstones and calcretes in the fluvial segment of the Castissent (Marzo et al., 1988). This major fluvial progradation is correlated westwards with deep-water turbidite sequences of the Arro and Fosado Formations in the Ainsa Basin (Fig. 3, Mutti et al., 1988; Nijman and Nio, 1975; Nijman and Puigdefabregas, 1978; and Pickering and Bayliss, 2009). In the upstream, eastern counterparts of the Castissent Fm., the channel complexes are intercalated with yellow to red coloured palaeosols. Sub-spherical to slightly elongated carbonate nodules with diameters ranging from 1 mm to 4 cm are omnipresent in the palaeosols (Fig. S1 in the Supplement). Studies of the Castissent Fm. tentatively attributed their occurrence to an important pulse of exhumation and thrust activity in the hinterland at ca. 50 Ma, in possible combination with a late Ypresian sea-level fall (Puigdefabregas et al., 1986; Marzo et al., 1988; Whitchurch et al., 2011; and Castelltort et al., 2017), both resulting in reduced available accommodation space enhancing progradation and amalgamation (Chanvry et al., 2018).

The Chiriveta section, encompassing the Castissent Fm., is situated in a continental palaeogeographic position prone to pedogenesis and slightly off-axis from the more “in-axis” amalgamated sand-rich-type section of Mas de Faro (Fig. 2); for palaeo-position and correlation see also Figs. 10 and 12 in Marzo et al. (1988).

In the Chiriveta location, stratigraphic constraints are limited to the identification of European Mammals zone MP10 (Badiola et al., 2009), which provides an age range of 50.73 to 47.4 Ma (GTS2012). This age span is refined by bio- and magnetostratigraphic data from the Castissent Fm. outcrops of the Campo location, about 40 km further west (Kapellos and Schaub, 1973; Tosquella, 1995; Bentham and Burbank, 1996; Tosquella et al., 1998; and Payros et al., 2009) (Fig. 3). Because of its outcropping extent, the Castissent Fm. has been mapped from west to east across these sections (Nijman and Nio, 1975; Nijman, 1998; Poyatos-Moré, 2014; and Chanvry et al., 2018). The low slope of the Castissent Fm. (ca. $\mathrm{2.3}\times {\mathrm{10}}^{-\mathrm{4}}$ m m⁻¹; see Supplement Table S1) indicate an elevation drop of ca. 1 m between the Chiriveta section and the Campo section. Given average flow depths of 3.75 m in the Castissent channels based on measurement in the Chiriveta and La Roca sections, we thus assume no significant time lag of deposition between both sections. In the Campo section, Kapellos and Schaub (1973) find the transition between the D. lodoensis and the T. orthostylus nannoplankton (NP) zones at ca. 200 m below the base of the Castissent Fm. and the transition between the T. orthostylus and the D. sublodoensis NP zones in the transgression ca. 100 m above the uppermost member of the Castissent Fm. This indicates that the Castissent Fm. was deposited during NP13. Magnetostratigraphic data of the same section by Bentham and Burbank (1996) place the transition between the C22r and C22n magnetozones closely above the top of the Castissent Fm. We thus used the recent astrochronological age models of Westerhold et al. (2017), which obtain numerical ages of 50.777±0.01 and 49.695±0.043 Ma for the base and top of C22r, respectively, and obtain a numerical age of 50.534±0.025 Ma for the base of NP13 based on the Ocean Drilling Program site 1263. Considering the data available and their resolution, we suggest a depositional age span between 50.5 and 49.7 Ma for the Castissent Fm. (reported in green on Fig. 1). According to global isotopic records (Fig. 1), this period was marked by four hyperthermals labelled S/C22rH3, T/C22rH4, U/C22rH5, and V/C22nH1 (Cramer et al., 2003; Lauretano et al., 2016; and Westerhold et al., 2017).

3 Material and methods

3.1 Sampling

A total of 74 samples were collected from the lower Eocene Chiriveta section for geochemical studies. All samples consist of floodplain material and were taken below the weathering depth (∼50 cm), with an average resolution of 1 m. Resolution was increased by a factor of 2 in specific horizons such as red beds. When important sandbodies occurred, lateral equivalent floodplain material or intercalated palaeosol horizons were sampled. Each sample was split in two aliquots, one for major and trace element analysis and the other for carbon and oxygen stable isotope analysis on pedogenic carbonate nodules. The carbonate nodules were extracted from the bulk palaeosol material by sieving and then cleaned by repeated washes with deionized water in an ultrasound bath. From each cleaned nodules set, subsamples of one to four nodules were taken, leading to a total of 149 subsamples of pedogenic carbonate nodules.

3.2 Carbon and oxygen stable isotopes

Pedogenic carbonate nodules were crushed and powdered in an agate mortar and analysed for stable carbon and oxygen isotope composition at the Institute of Earth Surface Dynamics of the University of Lausanne (Switzerland) using a Thermo Fisher Scientific (Bremen, Germany) carbonate-preparation device and Gas Bench II connected to a Thermo Fisher Delta Plus XL isotope ratio mass spectrometer. The carbon and oxygen isotope compositions are reported in the delta (δ) notation as the per mil (‰) isotope ratio variations relative to the Vienna Pee Dee Belemnite standard (VPDB). The analytical reproducibility estimated from replicate analyses of the international calcite standard NBS-19 and the laboratory standard Carrara marble was better than ±0.05 ‰ (1σ) for δ¹³C and ±0.1 ‰ (1σ) for δ¹⁸O.

3.3 Major and trace element composition

Fifty-two bulk palaeosol samples were analysed for major and trace elements using X-ray fluorescence (XRF) spectrometry. Crushed bulk powders (< 80 µm) were mounted in a plastic cup covered by a thin polypropylene film (4 µm thick) and analysed in the laboratory with a Thermo Fisher Niton XL3t^® portable XRF analyser fixed on a test stand. Analyses were performed with a beam diameter of 8 mm, to determine the concentrations of 34 major and trace elements (from Mg to Au). Each measurement took 120 s, consisting of two 60 s cycles on four different filters (15 s on low, main, high, and light ranges), operating the X-ray tube at different voltages to optimize the fluorescence and peak/background ratios of the different elements. The limits of detection were of tens of parts per million for most elements, except for Mg, Si, and Al which are at wt % level. Sodium is too light to be detected. The acquired spectra were transferred to a computer using NDT software version 8.2.1. (Thermo Fisher Scientific, Waltham, Ma, USA). The same material has been analysed for 23 major and trace elements on fused and pressed discs, respectively, using a PANalytical PW2400 XRF spectrometer with a copper (Cu) tube at the University of Lausanne to cross-calibrate the compositions measured with the Niton XL3t^® portable XRF analyser.

3.4 Mean annual precipitation

The mean annual precipitation estimate (MAP) used in this study was estimated from the empirical relationship between MAP and CaO∕Al₂O₃ ratio for Mollisols from a national survey of North American soils according to the following equation: MAP (mm) $=-\mathrm{130.9}\times \ln \left(\mathrm{CaO}/{\mathrm{Al}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{3}}\right)+\mathrm{467}$ (Sheldon et al., 2002). CaO and Al₂O₃ concentrations were measured on bulk palaeosol material. Climate linked to the MAP estimate was classified based on the following boundaries: arid to semiarid at 250 mm and semiarid to subhumid at 500 mm (Bull, 1991).

3.5 Grain-size estimation

The relative grain-size variation in the sediment samples was estimated from their major element compositions. Si, Ti, and Zr are more concentrated in the coarse fraction of the sediment as they are found in larger mineral grains, whereas Al is more concentrated in the finer fraction of the sediment because it is mostly linked to clay minerals (Lupker et al., 2011, 2012; Croudace and Rothwell, 2015). Grain size variation throughout the section was estimated using Si∕Al, Ti∕Al, and Zr∕Al ratios, therefore, an increase in these ratios suggests a relative increase in the proportion of coarser material in the sample.

3.6 Correlation with target curves

The measured δ¹³C dataset was compared with a time-equivalent ODP 1263 global δ¹³C record reported by Westerhold et al. (2017) using the AnalySeries software (Paillard et al., 1996). The δ¹³C record of site 1263 was favoured over those of ODP 1209 and 1258 covering the Castissent Fm. time period, because it is continuous and has a higher resolution. Correlations between the δ¹³C record of site 1263 and the δ¹³C record of the Chiriveta section were performed in order to optimize the Pearson correlation coefficient (r) and to minimize abrupt variations in sedimentation rates. Well-defined peaks in both δ¹³C records were used as tie points for the correlation and the number of tie points was kept to a minimum (< 10) so as to not force the correlations.

Figure 4 Field images of the Chiriveta section (42^∘7^′56.57^′′ N, 0^∘41^′19.45^′′ E). (a) Outcrop view of Members A and B of the Castissent Formation. (b) Close-up view of the upper part of Castissent A Member. Fluvial channel-fill deposits, intercalated in reddish floodplain and overbank deposits and regional marine incursions (M1). (c) M0, first marine incursion at the base of the Castissent Fm. described by Marzo et al. (1988) expressed in the Chiriveta section by a tidal-influenced coarse sandstone with herringbone cross-stratification. (d) Yellow mottled palaeosol between CIE C and D. (e) Red floodplain interval equivalent of the CIE C. (f) 2 m thick grey interval interpreted as poorly drained brackish water facies and equivalent to the marine incursion M1. (g) An ∼6 m thick laterally extensive Castissent B sandbody incised in the underlying floodplain deposits. (h) Mottled silt, interpreted as pedogenetic fluvial channel overbank deposits.

4 Results

4.1 Sedimentology of the Castissent Formation at Chiriveta

We describe here the section logged and sampled in this work (Fig. 4). At Chiriveta, the Castissent Fm. is a palaeosol-rich succession, which shows greyish-yellow to red-brown mottled floodplain palaeosols (Fig. 4a–b), corresponding laterally to thick, medium to coarse-grained quartz-rich channel-fill deposits (width/depth ratio =20–50; Marzo et al., 1988) and overbank deposits flowing parallel to the main structures of the growing Pyrenean orogeny (Marzo et al., 1988). At the base of the section, the first marine incursion M0 is situated at the top of a 20 m thick coarse-grained tidal bar deposit with herringbone cross-stratifications and oyster shells (Fig. 4c). In the Chiriveta section, the Castissent Member A is a 48 m thick interval comprising two main medium-grained sandbodies of light colouration of 5.40 and 1.5 m in thickness respectively. Bedforms observed in the first sandbody have a mean height of 24 cm (n=9). The second marine incursion M1 is located at 48 m just below the Castissent B Member and consists of a 2 m thick grey interval interpreted by Marzo et al. (1988) as brackish–lagoonal water facies (Fig. 4b–f). The Castissent B Member (Fig. 4g) is a 12 m thick and laterally extensive (width/depth ratio ≥250; Marzo et al., 1988) amalgamated sandbody with a micro-conglomeratic erosive base. Grain size is overall larger than in Member A and ranges from fine sand to large pebbles. Sandbody tops show a fining-upward trend and are capped by mottled siltstone packages. Mottled siltstone layers are interpreted as pedogenized overbank deposits based on roots traces and their capping relationship with underlying sandbody deposits (observed at 26, 76, 89, and 96 m in Figs. 5 and 4h). More regular and sheet-like sandbodies interbedded with mottled siltstone layers are observed upwards. The section ends with a 23 m thick, medium to very coarse, tidally influenced sandstone deposit interpreted as the equivalent M3 marine incursion by Marzo et al. (1988). Although Castissent Member C was not interpreted by Marzo et al. (1988) in this section, a 2 m thick fine-grained sandbody at ca 80 m in our section could be the condensed lateral equivalent of it (Fig. 5).

Figure 5 Isotopic and geochemical data from the Chiriveta section. For the isotope dataset, the curves passes through the mean values at each sample position. Samples with minimum in δ¹³C values below 1 and 2 standard deviations are labelled A to F. Mean annual precipitation (MAP) was estimated from the empirical relationship between MAP and CaO to Al₂O₃ ratio (Sheldon et al., 2002).

4.2 Stable isotopic record

Carbon and oxygen isotope ratios from the carbonate nodules are presented in Fig. 5. The δ¹³C values vary between −10.9 ‰ and −1.9 ‰ with a mean value and 1 SD of $-\mathrm{7.7}\pm \mathrm{1.6}$ ‰. Six CIEs (named A to F in Fig. 5 and colour coded in Fig. 6) are more negative than −9.3 ‰ (i.e. the mean value – 1 standard deviation) amongst which one (CIE D) is below 2 SDs. The values are −9.6 ‰, −9.8 ‰, −9.9 ‰, −10.9 ‰, −9.9 ‰, and −9.4 ‰ for CIEs A to F respectively. At the bottom of the section, CIE A is followed by a relatively constant interval of mean δ¹³C values. CIE B, situated in the first red bed, marks the beginning of a stepped δ¹³C trend (around ±1 SD) leading to the minimum CIE D. The second part of the section shows two more CIEs separated by the highest δ¹³C value at 74 m. CIE F is the least prominent of all CIEs. The δ¹⁸O values vary between −7.0 ‰ and −5.0 ‰ with a mean value of $-\mathrm{6.0}\pm \mathrm{0.4}$ ‰, which makes them less dispersed than the δ¹³C record. Nine negative oxygen isotope excursions are more negative than the mean value minus 1 SD, amongst which one is below 2 SD reaching a minimum value of −6.8 ‰ at 19 m. The oxygen isotope excursions do not correspond with CIEs described above.

4.3 Major and trace elements

Titanium (Ti), aluminium (Al), and zirconium (Zr) concentrations measured on bulk palaeosols are plotted in Fig. 5. These elements are commonly considered as immobile and are expected to concentrate in more weathered soils. Ti values vary between 0.18 % and 0.52 % with a mean value of 0.34 % and a standard deviation of 0.08. Al values vary between 3.03 % and 9.35 % with a mean value of 5.85 % and a standard deviation of 1.53. Zr values vary between 67 and 204 ppm with a mean value of 128 ppm and a standard deviation of 35. Mean annual precipitation (MAP) estimates values vary between 185 and 754 mm yr⁻¹ with a mean value of 376 mm yr⁻¹ and a standard deviation of 111. Ti, Al, Zr, and MAP show a similar trend starting from the base of the section with a global increase for all values toward CIE C and a decrease afterwards. All CIEs show higher values of Ti, Al, Zr, and MAP except CIE F. Based on Bull (1991), an average value of 387 mm yr⁻¹ for the MAP in the Chiriveta section represents a semi-arid climate (Fig. 5). All CIEs show an increase in MAP.

5 Discussion

5.1 Carbon and oxygen isotopic record

5.1.1 Identifying the CIE

In continental successions, the carbon isotope composition of pedogenic carbonate nodules – which consists of calcareous concretions between 1 mm and 4 cm diameter formed in situ in the floodplain – have been shown to be sensitive to environmental conditions during their formation (e.g. Millière et al., 2011a, b) and are therefore a promising tool to track how environments respond to carbon cycle perturbation. The carbon isotope composition of the soil carbonate nodules depend on the δ¹³C value of the atmospheric CO₂ and soil CO₂, which in turn is a function of the δ¹³C of the atmospheric CO₂ and the overlying plants, as well as the soil respiration flux and the partial pressure of atmospheric CO₂ (Cerling, 1984; Bowen et al., 2004; Abels et al., 2012; and Caves et al., 2016).

The δ¹³C vs. δ¹⁸O diagram for the pedogenic carbonate nodules from the Chiriveta section ($r=-\mathrm{0.26}$, n=149) suggests a good preservation of the primary isotopic signal (Fig. 6), with an average value of δ¹³C $=-\mathrm{7.7}\pm \mathrm{1.6}$ ‰ similar to mid-latitude late Palaeocene to Eocene continental δ¹³C values (excluding the PETM samples) observed elsewhere (e.g. McInerney and Wing, 2011; and references therein) and a spread comparable with δ¹³C values from carbonate nodules analysed for the same period in the Bighorn Basin (Bowen et al., 2001). Figure 6 emphasizes that lower Eocene carbonate nodules display overall more negative δ¹³C values than the Holocene nodules, which is consistent with a large compilation of data from eastern Eurasia (Caves Rugenstein and Chamberlain, 2018). Pre-PETM δ¹⁸O values from carbonate nodules from the same area ($-\mathrm{4.5}\pm \mathrm{0.4}\mathrm{‰}$) (Hunger, 2018) show values similar to measurements from the Chiriveta section ($-\mathrm{6.0}\pm \mathrm{0.4}\mathrm{‰}$). Oxygen and carbon isotopes are not coupled during hyperthermal events in continental record as already observed by Schmitz and Pujalte (2003) Bowen et al. (2001) for the PETM isotopic excursion. Though the precise mechanisms that produce stable δ¹⁸O during CIE are still debated, mid-latitude precipitation δ¹⁸O appears to be relatively insensitive to changes in atmospheric pCO₂ and warming, particularly in greenhouse climates (Winnick et al., 2015). In contrast, the stable δ¹⁸O values of soil carbonates from the Pyrenean foreland basin ($-\mathrm{5.5}\pm \mathrm{0.9}\mathrm{‰}$) are likely additionally stabilized by its position close to the coast (Cerling, 1984; Kukla et al., 2019) compared for example to those of the Bighorn Basin ($-\mathrm{9.0}\pm \mathrm{0.6}\mathrm{‰}$). This is in line with a more continental palaeogeographical position of the Bighorn Basin compared to the Tremp-Graus Basin at the time (Seeland, 1998).

Figure 6 Continental δ¹³C and δ¹⁸O values from the early Eocene Castissent Fm. in the Chiriveta section (this study) plotted with pre- and syn-PETM δ¹³C and δ¹⁸O values from the same area (Khozyem Saleh, 2013; Hunger, 2018) and pre-, syn-, and post-PETM values from the Bighorn Basin (Bowen et al., 2001) as well as recent pedogenic carbonate isotopic values (Cerling and Quade, 1993; Gallagher and Sheldon, 2016).

A hyperthermal event recorded in marine sediments is defined by paired negative carbon and oxygen stable isotope excursions that are more negative than the mean value minus 1 SD (Kirtland-Turner et al., 2014). This definition may not be applicable to continental deposits, because continental systems respond differently than marine systems to the carbon cycle perturbations. Though the marine δ¹³C record is thought to record the global CO₂ δ¹³C, the δ¹³C value of the marine dissolved inorganic carbon is also influenced by dissolution of carbonates at depth (McInerney and Wing, 2011). In contrast, δ¹³C in pedogenic nodules varies with soil properties, atmospheric and soil pCO₂ and δ¹³C, and the rate and nature of carbon input and/or output by soil respiration (Bowen et al., 2004; Sheldon and Tabor, 2009). These processes create complexities in estimating CIEs in soil carbonate nodules and in marine carbonates (McInerney and Wing, 2011). Nevertheless, we used the hyperthermal definition from Kirtland-Turner et al. (2014) as a starting point to filter the high-resolution variations in the Chiriveta section. We identify 16 samples with CIE values more negative than the mean minus 1 SD. Among these 16 samples, we recognized six discrete CIEs (named A–F in Figs. 5 and 7). Both marine incursion M1 and M2 show an abrupt shift from −9 ‰ to −10 ‰ in continental δ¹³C values towards more (positive) marine values of −4 ‰ to −2 ‰; this points to a progressive higher contribution of seawater to the formation of the carbonate nodules.

Figure 7 (a) Scaling of the Chiriveta isotopic section with the time equivalent interval of site 1263 (Westerhold et al., 2017). The correlation was calculated using the AnalySeries software (Paillard et al., 1996) and centred on CIE D and hyperthermal U. Mean, minus 1 and 2 SD lines on the global record were calculated sets over the selected time period. The correlation coefficient (r) between the two curves is 0.65. (b) Hyperthermal U amplitude in palaeosol carbonate and benthic foraminifera (inset B after Abels et al., 2016).

Six correlation options with the global record were explored in the time window of the Castissent Fm. (Figs. S2 and S3). The correlation presented in Fig. 7a was favoured as it shows (i) reasonable sedimentation rates variations, (ii) a similar amplitude to the CIE in the global record, and (iii) the highest correlation coefficient (r=0.65, n=71). Moreover, it plots along the same trend regarding hyperthermal CIE amplitudes in marine and continental environments, suggesting a common mechanism of global climatic change with events I1, I2, H2, and ETM2 (Fig. 7b). Based on these observations and obtained correlation, we suggest that only hyperthermal U is preserved in the Chiriveta section and that it is correlated with CIE D. Sedimentation rate obtained with the favoured correlation (Fig. 7) varies between 0.1–0.29 mm yr⁻¹, consistent with sedimentation rates reported for other Eocene floodplain successions (Kraus and Aslan, 1993). The correlation coefficient of r=0.65 suggests an overall good signal preservation in the studied continental section for a 40 kyr climatic event.

5.1.2 Mechanisms causing the CIE

An increase in temperature could potentially release a significant amount of CO₂ into the atmosphere (Trumbore et al., 1996; Melillo et al., 2014). The amplitude and duration of Eocene CIEs are approximately 30 % of the one recorded for the PETM; we hypothesize that the climatic effects of smaller-scale hyperthermals can be linearly scaled to the PETM. Based on this assumption and in order to get a rough approximation without considering a nonlinear sensitivity response, a smaller-scale hyperthermal would imply a release of approximately 500 to 1500 Gt of carbon to the ocean and atmosphere reservoir and a global temperature rise of about 1.5–2.5 ^∘C. This estimation corresponds to the 1500–4500 Gt of carbon released during the PETM, causing a rise of 5–8 ^∘C (Bowen et al., 2006), and is in line with previous estimations of ∼3 and ∼2 ^∘C warming for ETM2/H1 and H2 events respectively (Stap et al., 2010). A release of 500 to 1500 Gt of carbon in the form of methane would imply a marine CIE of 0.8 ‰ to 2.3 ‰ or 0.3 ‰ to 0.9 ‰ if the carbon origin is dissolved organic carbon (DOC) (Sexton et al., 2011). The latter seems more plausible regarding the observed amplitude of ∼1 ‰ measured in the marine record for hyperthermal U (Westerhold et al., 2017) and the supposed origin linked to the oxygenation of deep-marine DOC of post-PETM hyperthermals (Sexton et al., 2011). A global shift of −1 ‰ in δ¹³C can however not fully explain the 3 ‰ shift in δ¹³C observed in this study.

Figure 8 Components influencing the δ¹³C values of pedogenic carbonate nodules. Mean early Eocene bulk marine carbonate and small-scale hyperthermal (all except PETM) are from Westerhold et al. (2018). Fractionation value between organic matter and carbonate nodules are based on Sheldon and Tabor (2009). All other fractionation values are based on Koch et al. (1995). Mean carbonate nodule values come from this study.

The δ¹³C mean value in the Chiriveta section is $-\mathrm{7.7}\pm \mathrm{1.6}$ ‰. This value reflects an overall equilibrium with a mean atmospheric CO₂ of −7 ‰ (Koch et al., 1995) and is coherent with pre-PETM δ¹³C values of $-\mathrm{7.1}\pm \mathrm{0.9}\mathrm{‰}$ found in the same area (Hunger, 2018; Fig. 6). It is possible to calculate from the (small-scale) hyperthermal δ¹³C excursions in the marine environment the shift to be expected in soil carbonate nodules by using known fractionation coefficients (Koch et al., 1995, 2003); the expected δ¹³C value in carbonate nodules, only considering the respiration of organic matter, is −11 ‰ (Fig. 8). This value is within the range of those measured from the Chiriveta section, where some nodules reach values as low as −10.9 ‰. We suggest that the bacterial respiration of organic matter, enhanced by warmer temperatures (e.g. Davidson and Janssens, 2006; Trumbore et al., 1996), may also have contributed to the lower δ¹³C values of nodules during the CIEs (Fig. 8). On geological timescales, soil organic carbon can be considered at steady state with equal organic carbon inputs and outputs from the soil (Koven et al., 2017). Respiration (carbon output after mineralization as CO₂) is thought to be more sensitive to global warming than gross primary productivity (organic carbon input as organic matter), leading to a depletion of the total soil carbon pool with time during transient global warming events; although the precise sensitivity of gross primary productivity remains poorly constrained (Davidson and Janssens, 2006). Large uncertainties remain about carbon dynamics and their timescale in the soils during climate changes. Parameters such as the vegetation type (Klemmedson, 1989), temperatures (Koven et al., 2017), soil geochemistry (Torn et al., 1997; Doetterl et al., 2015), and soil water content (Davidson et al., 2000) have been shown to be important controlling factors within historical timescales.

Considering these caveats, we estimate the maximum possible contribution of enhanced soil carbon respiration to negative δ¹³C excursions during the CIEs. Using typical values for the organic carbon reservoir comprising fast and slow cycling carbon in soils in arid to semi-arid ecosystems of 5.6–19.2 kg C m⁻² (Klemmedson, 1989; Raich and Schlesinger, 1992), respiration fluxes starting at a steady-state value of 0.5 kg C yr⁻¹ and a respiration rate sensitivity ca. 5 % per degree (Raich and Schlesinger, 1992) (Q₁₀=1.5), we estimate that all of the organic carbon in soils would be consumed within 250 to 850 years, given an increase of 1 ^∘C and without changing the carbon input rate. Though there are a number of assumptions in this first-order estimate, the timescale of soil carbon depletion is substantially shorter than our estimate of the timescale of the CIE (∼36 kyr) (Fig. 7). As evidenced by this calculation, an increase in soil respiration triggered by warmer temperatures cannot be the sole mechanism driving the CIE shift over multi-millennial timescales. Instead, we suggest that during these transient warmings, this mechanism is associated with a high primary productivity – resulting in a greater input of carbon to the soil – leading to an overall higher soil respiration of organic matter. Coupled with lower atmospheric δ¹³C during hyperthermals, this mechanism caused a pronounced CIE in soil carbonate nodules.

5.2 Geochemical signature of hyperthermal events

Major and trace elements compositions of floodplain sediments is a function of river dynamics, climate, and sediment grain-size (Lupker et al., 2012; Turner et al., 2015). Based on the CIEs, we defined six intervals showing a relative enrichment (up 10 % to 30 % compared to the average value) in immobile elements such as Ti, Al, and Zr (Fig. 5). To ensure that major and trace concentrations are not grain-size biased, we plotted grain-size proxies Si∕Al, Ti∕Al, and Zr∕Al (Lupker et al., 2012; Turner et al., 2015), which all exhibit a relatively stable trend, not connected with the immobile element concentrations (Fig. S4). The enrichments in Ti, Al, and Zr suggest mature palaeosols with potential intense weathering due to enhanced humid climatic conditions; but this may also correspond to a longer exposure time on a stable floodplain, allowing leaching of mobile elements and relative enrichment of immobile elements (Sheldon and Tabor, 2009). Pedogenic nodules are frequent in well-drained soil profiles associated with a climate regime where the potential evapotranspiration is greater than the mean annual precipitation rate (Slessarev et al., 2016) and with a mean annual precipitation < 800 mm yr⁻¹ (Cerling, 1984; Retallack, 1994; and Sheldon and Tabor, 2009). These conditions correspond to climate ranging from arid to subhumid conditions (Hasiotis, 2004; Prochnow et al., 2006; Hyland and Sheldon, 2013). This agrees with MAP values obtained for the palaeo-precipitation estimate (Fig. 5) and with a smectite/kaolinite > 1 assemblage dominating some of the studied soils (Nicolaides, 2017, Table S2); all suggestive of a semi-arid to subhumid climate with seasonal humidity (Arostegi et al., 2011). Associated with CIEs C and D in red bed deposits, sub-millimetric iron-oxide and iron-hydroxide nodules made of concentric hematite and goethite were found together with carbonate nodules (Fig. S1). This suggest a seasonal climate as hematite forms under more arid soil condition than goethite (Kraus and Riggins, 2007). Together, these observations are in line with an acceleration of the hydrological cycle and a higher seasonality, as has been observed during the PETM, H1, H2; I1 and I2 hyperthermals (Bowen et al., 2004; Nicolo et al., 2007; Slotnick et al., 2012; and Dunkley Jones et al., 2018). Therefore, combined with CIEs, we suggest that small-scale hyperthermals in continental records can be recognized by an increase in the weathering index (Hessler et al., 2017) and by an increase in the immobile element concentrations, both related to an increase in precipitation intensity.

5.3 High-resolution hyperthermal signal

The high-resolution isotopic and elemental record of the Chiriveta profile allow us to highlight the dynamics and variability in a hyperthermal event. We do not observe a unique peak in δ¹³C but rather a stepped isotopic signal suggesting, together with above-discussed geochemical data, a climatic oscillation alternating with variably intense precipitations and leaching conditions during a climax spanning ca. 150 kyr (interval CIE B to D). Such a climatic behaviour, was already described for the PETM, during the pre-onset excursion (Bowen et al., 2015) and in the core CIE of the PETM (Giusberti et al., 2016). Moreover, the δ¹³C climax (CIE D) does not correspond to the highest concentrations of immobile elements or maximum MAP estimates, which we estimate occur during CIE C, which predates the CIE D by ca. 50 kyr (Fig. 7) . The minimum δ¹³C value therefore does not seem to be coeval with the most extreme climatic response, suggesting a complex environmental response. However, because sedimentation in floodplain depositional settings is a function of the channel position and flood frequency, the relative concentration of elements may reflect the changes in river dynamics instead of climatic variability, which could explain the mismatch between minimum values in CIE and the climatic response. More high-resolution hyperthermal studies in coeval continental sections are needed to better understand the relationships between proxies.

5.4 Possible implication for the preservation potential of hyperthermals in continental sections

Major events such as the PETM have proven to be detectable in both marine and continental environments (e.g. Abels et al., 2016; Koch et al., 1992), but the signal and preservation potential of smaller-scale climatic events (e.g. hyperthermal events L to W in Lauretano et al., 2016) may be more difficult to detect (Foreman and Straub, 2017) because of the inherent highly dynamic nature of sedimentation in fluvial deposits. To address this issue in the present case study, we calculated the compensation timescale (T_c) of the Castissent Fm. T_c is a timescale characteristic of an alluvial basin below which stratigraphic signals with shorter durations may be of autogenic origin, thereby giving a scale below which allogenic forcing should be interpreted carefully (Wang et al., 2011; Foreman and Straub, 2017; and Trampush et al., 2017). In other words, an external forcing signal with a duration smaller than T_c will be challenging to identify from background variability; the external forcing must be therefore of a duration longer than T_c and optimally twice T_c (Foreman and Straub, 2017). The T_c max can be calculated by dividing the topographic roughness or maximum channel depth by the average subsidence or deposition rate (Wang et al., 2011). Using an average sedimentation rate of 0.17 mm yr⁻¹ and an average channel depth of 3.75 m, we obtained a mean T_c of 22 000 years, which means that hyperthermal events of 40 kyr duration (timescale of hyperthermal U and preceding CIE) have the potential to be recorded despite fluvial system dynamics. Our estimate of preservation potential assumes steady sedimentation rates throughout the section. But, sedimentation in terrestrial records is not uniform (steady) but rather highly variable, resulting in spatial and temporal changes in facies and deposition rates ranging from < 0.1 to 1–2 mm yr⁻¹ (Marriott and Wright, 1993; Bowen et al., 2015; Kraus et al., 2015). However, mean accumulation rates give a reasonable estimate approximating more realistic (i.e. variable) sedimentation rates as observed in the Bighorn Basin (Bowen et al., 2015). Additionally, we analyse the vertical movement of the nearby structures to evaluate their potential influence on disrupting deposition at Chiriveta during Castissent times. The Chiriveta section was deposited near or at the axis of the Tremp-Graus basin (Nijman, 1998), which is bounded by the Bóixols thrust in the north and the in the south (Marzo et al., 1988). The Tremp-Graus basin is transported as a piggyback basin on the Montsec thrust emerging at the time approximatively 4 km south of the studied section (Nijman, 1998). In the basin axis, subsidence is the highest with rates of 0.1 to 0.29 mm yr⁻¹ (this study and Marzo et al., 1988). Taking into account a vertical movement rate of the Montsec thrust of 0.03 to 0.1 mm yr⁻¹ during the Castissent time-interval (based on a horizontal displacement of 7 km, a period of activity lasting 26 Ma, and a thrust dip between 6 and 20^∘; Farrell et al., 1987; Nijman, 1998; Clevis et al., 2004; and Whitchurch et al., 2011), we estimate that the vertical displacement is no more than equal to sedimentation rates in the basin axis. This is consistent with the general absence of growth strata in the basin axis, although growth strata can indeed be observed closer to the Montsec (Nijman, 1998).

The rates of accumulation, distance to the main structures, and characteristic compensation timescale together suggest that hyperthermal events of ca. 40 kyr duration can be recorded in the Castissent Fm. These results confirm that, despite its highly dynamic nature, fluvial sedimentation may contain valuable record of high-frequency events, even in active tectonic contexts.

6 Conclusions

A new high-resolution isotopic record from the palaeosol-rich deposits at the Chiriveta section identified a prominent negative carbon isotope excursion (CIE) in continental settings. We suggest that the CIE recorded in fluvial succession of the early Eocene Castissent Formation is the “U” event, identified for the first time in continental deposits. This climatic event reaches δ¹³C values of 2σ (standard deviation) below the mean and is heralded and followed by several smaller-scale stepped CIEs, which are interpreted as moments of enhanced primary productivity, leading to an overall higher soil respiration. We show that all these CIEs are relatively enriched in immobile elements (i.e. Ti, Zr, and Al) and display an increase in MAP estimates. These observations coupled with the presence of iron oxide nodules on an overall weathered succession, suggest an increase in precipitation rates during these events. The data presented in this study suggest a period of ca. 150 kyr of contrasted climate alternating average and above background weathering conditions. Finally, the results of this demonstrate the importance of hyperthermal events in continental successions as well as in the preservation potential of such deposits.