Fluvio-deltaic record of increased sediment transport during the Middle Eocene Climatic Optimum (MECO), Southern Pyrenees, Spain

The early Cenozoic marine sedimentary record is punctuated by several brief episodes (<200 kyr) of abrupt global warming, called hyperthermals, that have disturbed ocean life and water physicochemistry. Moreover, recent studies of fluvial–deltaic systems, for instance at the Palaeocene–Eocene Thermal Maximum, revealed that these hyperthermals also impacted the hydrologic cycle, triggering an increase in erosion and sediment transport at the Earth's surface. Contrary to the early Cenozoic hyperthermals, the Middle Eocene Climatic Optimum (MECO), lasting from 40.5 to 40.0 Ma, constitutes an event of gradual warming that left a highly variable carbon isotope signature and for which little data exist about its impact on Earth surface systems. In the South Pyrenean foreland basin (SPFB), an episode of prominent deltaic progradation (Belsué–Atarés and Escanilla formations) in the middle Bartonian has been usually associated with increased Pyrenean tectonic activity, but recent magnetostratigraphic data suggest a possible coincidence between the progradation and the MECO warming period. To test this hypothesis, we measured the stable-isotope composition of carbonates (δ¹³C_carb and δ¹⁸O_carb) and organic matter (δ¹³C_org) of 257 samples in two sections of SPFB fluvial–deltaic successions covering the different phases of the MECO and already dated with magnetostratigraphy. We find a negative shift in δ¹⁸O_carb and an unclear signal in δ¹³C_carb around the transition from magnetic chron C18r to chron C17r (middle Bartonian). These results allow, by correlation with reference sections in the Atlantic and Tethys, the MECO to be identified and its coincident relationship with the Belsué–Atarès fluvial–deltaic progradation to be documented. Despite its long duration and a more gradual temperature rise, the MECO in the South Pyrenean foreland basin may have led, like lower Cenozoic hyperthermals, to an increase in erosion and sediment transport that is manifested in the sedimentary record. The new data support the hypothesis of a more important hydrological response to the MECO than previously thought in mid-latitude environments, including those around the Tethys.

1 Introduction

The Middle Eocene Climatic Optimum (MECO) is a global warming event that occurred during the Bartonian (ca. 40.5–40.0 Ma) and briefly reversed the longer-term cooling trend of the middle to upper Eocene (Fig. 1; Arimoto et al., 2020; Bijl et al., 2010; Bohaty and Zachos, 2003; Bohaty et al., 2009; Bosboom et al., 2014; Galazzo et al., 2014; Henehan et al., 2020; Edgar et al., 2010, 2020; Giorgioni et al., 2019; Jovane et al., 2007; Mulch et al., 2015; Sluijs et al., 2013; Sppoforth et al., 2010; Pälike et al., 2012; van der Boon et al., 2021). Marine bulk and benthic oxygen isotope compositions (δ¹⁸O values) show a negative excursion of −1.5 ‰ over the event, which was interpreted as a gradual global warming of 3 to 6 ^∘C (Bohaty et al., 2009). In contrast, the evolution of carbon isotope composition (δ¹³C values), unlike earlier hyperthermals of the Cenozoic (e.g. Palaeocene–Eocene Thermal Maximum – PETM; Eocene Thermal Maximum – ETM 2), differs from site to site, showing opposite patterns between hemispheres and displaying a carbon isotope excursion (CIE) in some but not all marine records (Bohaty et al., 2009; Henehan et al., 2020; Westerhold and Röhl, 2013). This CIE suggests a rise in atmospheric partial pressure of carbon dioxide (pCO₂) during the warming peak (Henehan et al., 2020; Bijl et al., 2010), and numerous potential CO₂ sources have been proposed – among them, a prolonged pulse of metamorphic decarbonisation possibly linked with the Himalayan collision at that time (Bijl et al., 2010; Bohaty et al., 2009; Bouilhol et al., 2013; Sternai et al., 2020), an increase of volcanism (van der Boon et al., 2021), or lower continental weatherability (van der Ploeg et al., 2018). However, the pCO₂ record remains ambiguous and difficult to link in a straightforward way to a rapid injection of exogenous carbon during the MECO (e.g. Henehan et al., 2020). In addition, regardless of the CO₂ sources involved, the MECO coincides with a very long (2.4 Myr) eccentricity cycle, suggesting a possible orbital trigger (Westerhold and Röhl, 2013). Therefore, considering the unresolved MECO driving mechanism(s) and how the Earth system responded to this carbon cycle perturbation, the MECO poses a significant challenge to understanding carbon cycle variations on timescales of several hundreds of thousands of years (Sluijs et al., 2013; Henehan et al., 2020; Sternai et al., 2020). Addressing this challenge requires extensive documentation of the MECO in a range of environments and geodynamic contexts, as well as documentation of its effect on Earth surface dynamics.

Figure 1 (a) Cenozoic δ¹⁸O values compilation from the Pacific Ocean, compiled in Miller et al. (2020). The continuous blue bar represents permanent ice sheet presence, and the discontinuous blue bar represents ephemeral ice sheet, modified from Zachos et al. (2001). Main climate events: EECO – Early Eocene Climatic Optimum; PETM – Palaeocene–Eocene Thermal Maximum; MECO – Middle Eocene Climatic Optimum; E–O – Eocene–Oligocene transition. (b) Carbonate δ¹⁸O values of site 702 from Bohaty et al. (2009). General climatic context of the Middle Eocene Climatic Optimum. The MECO “event”, in yellow from ca. 40.5 to 40.0 Ma in the inset, is considered the last hyperthermal of the Eocene, immediately preceding the shift to genuine Antarctic glaciation and the ice-house world of the Oligocene.

Current data converge towards the view that, during the MECO, surface and deep oceanic waters experienced a gradual and uniform warming between 3 to 6 ^∘C for all latitudes (Arimoto et al., 2020; Bijl et al., 2010; Bohaty et al., 2009; Rivero-Cuesta et al., 2019). Moreover, deep-sea carbonates were affected either by significantly reduced carbonate accumulation rates or even by dissolution, suggesting broad acidification of sea-bottom waters, involving an estimated ca. 1 km shoaling of the carbonate compensation depth (CCD; Henehan et al., 2020; Cornaggia et al. 2020; Pälike et al., 2012; Arimoto et al., 2020). However, while the temperature increase in the oceans has been inferred in multiple sites, the MECO environmental perturbation affected differently the fauna communities (Arimoto et al., 2020; Cramwinckel et al., 2019). In some locations, the warmer conditions reduced nutrient availability, decreasing the benthic productivity (Arimoto et al., 2020; Bijl et al., 2010, Galazzo et al., 2014; Moebius et al., 2015; Cramwinckel et al., 2019). In contrast, the Southern Ocean (Moebius et al., 2014) or the Neo-Tethys Ocean (Galazzo et al., 2013) recorded increased productivity during the MECO.

In contrast to the oceanic realm, the expression of the MECO in non-marine records remains scarce and variable. Mulch et al. (2015) suggested a boost of precipitation in the North American plateau derived from low δ¹⁸O values, while Bosboom et al. (2014) documented a shift towards arid conditions in the Tarim basin with a reduction in fern palynomorphs. Such a drying trend in central Asia is opposite to the Neo-Tethys Ocean dynamic, where a greater burial of organic matter (OM) immediately following the MECO may have been caused by increased nutrient runoff due to an enhanced hydrological cycle during the warm period (Galazzo et al., 2014; Giorgioni et al., 2019; Spofforth et al., 2010). These studies raise the question of the response of weathering, erosion, and sediment transport in terrestrial systems to global warming, as has also been posed for other hyperthermals recently (e.g. Chen et al., 2018; Foreman et al., 2012; Foreman and Straub, 2017; Honegger et al., 2020). This issue highlights the need for further documentation of the clastic sedimentary successions that temporally cover single and long-term climate crises (i.e. Early Eocene Climatic Optimum, Palaeocene–Eocene Thermal Maximum, etc.; Fig. 1).

In this work, we aim to understand the effects of the MECO on surface systems by exploring the interface between ocean and continent. The shallow marine settings, very sensitive to sea level changes and sediment supply, potentially provide a unique perspective of the hydrological response to climate change in the continental domain and the geochemical and isotopic evolution in the marine domain. We focus on two separate deltaic successions in the southern (Belsué locality, B) and northern (Yebra de Basa locality, YB) margins of the Jaca basin in the South Pyrenean foreland basin (SPFB; Fig. 1). The successions are characterised by excellent exposure and have already been dated by high-resolution magnetostratigraphy. Both sections reveal progradation of deltaic and fluvial systems coeval with the magnetic reversal occurring at chrons C18r and C18n.2n, near or at the zenith of MECO warmth (Edgar et al., 2010, 2020; Garcés et al., 2014; Vinyoles et al., 2021). We generated new high-resolution profiles of δ¹³C_carb, δ¹⁸O_carb, major and trace elements, clay mineralogy, and Rock-Eval parameters across the chron C18r–C18n.2n reversal to identify geochemical changes associated with the MECO onset and its recovery. We tested the possible causative links between progradation and the MECO perturbation. Finally, we discuss the sedimentary evolution of both sections to understand landscape response during the MECO and to explore the significance of its identification in the SPFB and the impact of climate shifts in source-to-sink systems as recorded at the continent–ocean interface.

2 Geological setting

The Pyrenees are a nearly E–W-trending mountain belt formed by the collision of the Iberian and European plates from the Late Cretaceous (Santonian) to the Early Miocene (Muñoz, 1992; Roure et al., 1989; Teixell, 1998; Vergés et al., 2002). The south Pyrenean zone is composed of an imbricate system of synorogenic thrusted cover sheets propagating southwards, detached above the Triassic evaporites (Labaume and Teixell, 2018; Lagabrielle et al., 2010; Mochales et al., 2012; Pueyo et al., 2002; Teixell et al., 2016, 2018). Among them, the emplacement of the South-Central Unit (SCU) by the early Eocene resulted in the partitioning of the South Pyrenean basin (SPFB) and the development of an E–W elongated deep basin draining west towards the Atlantic Ocean (Mochales et al., 2012; Muñoz et al., 2018; Puigdefàbregas, 1975; Puigdefàbregas and Souquet, 1986; Seguret, 1972). Due to the westward propagation of deformation during the middle Eocene and the differential velocity of the thrust sheets, oblique thrusted anticlines developed at the southwestern termination of the SCU (Muñoz et al., 2013). These thrusts caused the fragmentation and piggy-back transport of a wider foreland region, which included, from east to west, the Tremp–Graus basin, the Ainsa basin, and eventually the Jaca basins (Muñoz et al., 2018; Fig. 2).

Figure 2 Geological and stratigraphic context of the study area (a) Synthetic geologic map of the Pyrenees and location of the study area in the Jaca basin (red square). Black lines represent the main tectonic structures. Modified from Teixell (1996) and Bosch et al. (2016). (b) Detailed geologic map of the Jaca basin and the study area (Belsué and Yebra de Basa). Modified from Puigdefàbregas (1975) and Boya (2018). (c) Chronostratigraphy of the Ainsa and Jaca basins, showing the westward progradation of all depositional systems. The names of the main lithostratigraphic units are represented in black. The studied sections (in white) are included in the work of Vinyoles et al. (2021), while the location of the geochemical analyses carried out in this paper is highlighted in red. The figure is modified from Vinyoles et al. (2021) and the age model GPTS 2012 from Gradstein et al. (2012).

The Tremp–Jaca basin (TJB) preserves an outstandingly exposed, complete source-to-sink system from during MECO times. In the middle Eocene, the alluvial–fluvial system of Sis–Escanilla flowed down the Tremp–Graus and Ainsa basins, draining the eroded sediments from the uplifting northeastern Pyrenees (Beamud et al., 2003; Roigé et al., 2016; Coll et al., 2020; Puigdefàbregas, 1975; Puigdefàbregas and Souquet, 1986). Sediments were transported westward into the Jaca basin, forming a mixed delta-carbonate ramp with tidal influence (Puigdefàbregas, 1975; Castelltort et al., 2003). Two main deltaic systems developed during the lower Bartonian in the southern (Belsué–Atarés) and northern (Sabiñánigo) margin of the Jaca basin that primarily fed the distal Hecho turbidites in the western sector (Mutti, 1977; Remacha and Fernández, 2003). Shallow marine environments are mainly dominated by marly facies, which, thanks to their high carbonate content and relatively deep depositional environment, are suitable for geochemical studies (Wendler, 2013). High-resolution magnetostratigraphy in the Pyrenean region correlates with different sections along the entire source-to-sink system (Vinyoles et al., 2021). We selected two lower-Bartonian sections, Belsué (BS) and Yebra de Basa (YB), because they present excellent exposition and are provided with magnetostratigraphic dating, aiming to unravel the geochemical history of the two main deltaic systems coeval to the MECO.

Figure 3 Detailed geological maps from Belsué (a) and Yebra de Basa (b), modified respectively from Puigdefàbregas (1975) and Boya (2018). Red lines show the location of the studied sites, and the darker grey lines in Belsué represent the depositional sequences defined by Millán et al. (1994). Small stratigraphic logs on the bottom right represent large-scale synthetic logs with the formation name, which are modified from Vinyoles et al. (2021) for Yebra de Basa and from Castelltort et al. (2003) for Belsué.

The 130 m thick BS section is located within the External Sierras, east of the “Pico del Águila” anticline (42.30^∘ N, 0.37^∘ W; Fig. 3). This section has been extensively studied as a perfect case study for the tectonic-sedimentation relationship (Fig. 2; Puigdefàbregas, 1975; Puigdefàbregas and Souquet, 1986; Millán et al., 1994, 2000; Castelltort et al., 2003; Huyghe et al., 2012; Garcés et al., 2014). The lower boundary corresponds to an encrusted and ferruginous surface on top of a shallow marine bioclastic limestone (Guara Fm.), locally overlain by sandy marls rich in glauconite (Millán et al., 1994). It has been interpreted as a drowning unconformity of the Guara carbonate platform (Puigdefàbregas, 1975; Silva-Casal et al., 2019), close to the Lutetian–Bartonian boundary (Rodriguez-Pintó et al., 2013). This major unconformity led to the syntectonic deposition of the Arguís marls and the Belsué sandstones, while the Gabardiella and Pico del Águila anticlines were growing. Different authors studied the influence on the stratigraphy of local tectonic movement in the Belsué-Arguís region (Lafont, 1994; Castelltort et al., 2003), concluding that local tectonics modify the stacking pattern and position of its genetic units along with different depositional environments. The entire section covers the lower-Bartonian interval (Garcés et al., 2014) up to a maximum flooding surface MFS-2 by Muñoz et al. (1994) (Fig. 3), which corresponds to the deepest palaeobathymetry in the BS area (ca. 150 m, Sztràkos and Castelltort, 2001).

The 800 m YB section is located between the Basa anticline and the Santa Orosia syncline (42.49^∘ N, 0.28^∘ W; Fig. 3) and comprises one of the best outcropping sections for the Sabiñánigo sandstone (Puigdefàbregas, 1975; Lafont, 1994; Boya, 2018). It is composed of three subsections covering the lower-Bartonian interval, between the upper part of the chron C18r and the chron C18n.1r (Vinyoles et al., 2021). The first subsection is located west (Yb-C), within the Larrés marls, whose top is correlated with the base of the Vinyoles et al. (2021) magnetostratigraphic profile. From here, the next two subsections were built following the Vinyoles et al. (2021) profile, which comprises the two deltaic levels of the Sabiñánigo sandstone and most of the overlaying Pamplona marls.

3 Materials and methods

3.1 Field and sampling

We performed a complete stratigraphic study and sampling of the lower-Bartonian sections from BS and YB. The stratigraphic thickness of these sections was measured using the Jacobs staff in the field and geometric calculations when direct measurements were impossible. New mapping based on fieldwork and orthophotos was performed to correlate the different subsections. The Belsué section in the southern margin is composed of two subsections (Fig. 3); the lower one (Belsué east section, BS-E) was sampled at a resolution of 0.5–1.0 m and the upper one (Belsué west section, BS-O) every 3–6 m. In the northern margin, the higher average sedimentation rate (SR) in Yebra de Basa (>80 cm kyr⁻¹; Vinyoles et al., 2021) motivated a high-resolution sampling of 1–3 m at the middle part of the section (Yebra de Basa high-resolution section, YB-HR) and every 9–15 m in the other intervals (lower Yebra de Basa section, YB-C; upper Yebra de Basa section, YB-sup).

A total of 101 samples in BS and 157 samples in YB were collected. Each was composed of ca. 200 g of fine-grained and fresh rock below the weathering depth to avoid alteration and grain size bias (e.g. Lupker et al., 2011). The samples were mostly marls, corresponding to rocks rich in carbonate, OM, and clays. These samples were prepared for mineralogical and geochemical analyses in the laboratories of the University of Geneva and the University of Lausanne. The sample surface was cleaned with deionised water, and the weathered material was removed and then dried at 45 ^∘C for 2–3 d. The dried samples were crushed with a hydraulic press and powdered using an agate mill. The exposure conditions were usually ideal for sampling in both sections. However, there were difficulties in four intervals, resulting in gaps in the data. The problems were due to a dominance in sandy facies at the outcrop, corresponding to moments of maximum deltaic progradation or to poor exposure because of the fine-grained nature of the marls (e.g. Quaternary cover).

3.2 Clay mineralogy

The clay mineral assemblages of 24 representative samples per section were determined by X-ray diffractometry. The used system was a Thermo Scientific ARL X-TRA diffractometer at the Institute of Earth Science of the University of Lausanne (ISTE-UNIL), following the methods described by Klug and Alexander (1974), Kübler (1983, 1987), and Adatte et al. (1996). Ground chips were mixed with deionised water (pH 7–8) and agitated. The carbonate fraction was removed by treatment with 10 % HCl at room temperature and then for 20 min or more until all carbonate was dissolved. The insoluble residue was disaggregated (ultrasonication, 3 min), washed, and centrifuged (eight times) until a neutral suspension was obtained (pH 7–8). Different grain size fractions (<2 to 16 µm) were separated by the time-settling method based on Stokes' law. The selected fraction was then pipetted onto a glass plate and air-dried at room temperature. X-ray diffraction (XRD) analyses of oriented clay samples were made after air drying at room temperature under ethylene-glycol-solvated conditions. The intensities of XRD peaks (2θ; Moore and Reynolds, 1997) characteristic of each clay mineral (e.g. chlorite, mica, kaolinite) were used for a semi-quantitative estimation of the relative percent of clay minerals present in two size fractions (<2 and 2–16 µm).

3.3 Major- and trace-element compositions

Major- and trace-element concentrations from 24 representative samples per section were determined by X-ray fluorescence (XRF) spectrometry using a PANalytical PW2400 spectrometer from ISTE-UNIL. The major elements were analysed on fused glass discs. First, 2.7 to 3 g of sample powder was heated in a crucible oven at 1050 ^∘C for 1 night to obtain the loss on ignition (LOI) value. Then, 1.2000±0.0005 g of the calcinated sample was mixed with 6.0000±0.0005 g lithium tetraborate (Li₂B₄O₇) to prepare the fused glass disc using an automated glass-bead-casting machine (Pearl-X'3) at the University of Geneva.

The trace elements were analysed using pressed powder discs prepared at the University of Geneva from 3.000±0.0005 g of wax and 12.000±0.0005 g of non-calcinated sample powder. The mixture was homogenised and pressed (1 t hydraulic press, 20 s). Accuracy of the XRF analyses, assessed by analyses of standard reference materials, was 0.4 wt % for the major elements and 1 to 3 ppm for the trace elements.

3.4 Rock-Eval pyrolysis

The quality and quantity of the organic matter (OM) were determined in 237 bulk rock powders using a Rock-Eval 6 instrument at ISTE-UNIL, following the method described by Behar et al. (2001) and using the IFP 160 000 standard. Aliquots of samples were placed in an oven, heated at 300^oC under an inert atmosphere, and then gradually pyrolysed up to 650 ^∘C. After the pyrolysis, the samples were transferred into another oven and gradually heated up to 850 ^∘C in the presence of air, analysing the CO₂ and hydrocarbon (HC) concentration during the entire process. The calculated parameters included total organic carbon content (TOC, wt %), hydrogen index (HI in mg HC g⁻¹ TOC), oxygen index (OI in mg CO₂ g⁻¹ TOC), and T_max (^∘C) according to Espitalié et al. (1985) and Behar et al. (2001).

3.5 Carbonate carbon and oxygen stable isotopes

Carbonate carbon and oxygen stable-isotope ratios (δ¹³C_carb and δ¹⁸O_carb) of whole-rock powders containing >10 wt % CaCO₃ (n=237) were determined at the laboratories of the Institute of Earth Surface Dynamics of the University of Lausanne (IDYST-UNIL) using a Thermo Fisher Scientific GasBench II carbonate preparation device connected to a Delta Plus XL isotope ratio mass spectrometer. The CO₂ extraction was done by reaction with phosphoric acid at 70 ^∘C. The carbon and oxygen stable-isotope ratios were reported in the delta (δ) notation as the per mil (‰) relative to the Vienna Pee Dee belemnite standard (VPDB), where $\mathit{\delta }=(R_{\mathrm{sample}}-R_{\mathrm{standard}})/R_{\mathrm{standard}}\times \mathrm{1000}$ and $R={}^{\mathrm{13}}\mathrm{C}/{}^{\mathrm{12}}$C or ${}^{\mathrm{18}}\mathrm{O}/{}^{\mathrm{16}}\mathrm{O}$. The δ¹³C_carb and δ¹⁸O_carb values were standardised relative to the international VPDB scale by calibration of the reference gases and working standards with the international reference materials NBS 18 (carbonatite, ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}=-\mathrm{5.04}$ ‰, ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}=-\mathrm{23.00}$ ‰) and NBS 19 (limestone, ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}=+\mathrm{1.95}$ ‰, ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}=-\mathrm{2.19}$ ‰). Analytical uncertainty (1σ), monitored by replicate analyses of the international calcite standard NBS 19 and the laboratory standard Carrara Marble (${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}=+\mathrm{2.05}$ ‰, ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}=-\mathrm{1.7}$ ‰), was not greater than ±0.05 ‰ for δ¹³C and ±0.1 ‰ for δ¹⁸O.

3.6 Organic carbon stable isotopes

The organic carbon stable-isotope ratios (δ¹³C_org values in ‰ vs. VPDB) were determined in 155 samples, which were previously decarbonated by treatment with 10 % v/v HCl, thoroughly washed with deionised water, and dried (40 ^∘C, 48 h). The δ¹³C_org measurements were performed at the IDYST-UNIL by elemental analysis and isotope-ratio mass spectrometry using a Carlo Erba 1108 (Fisons Instruments, Milan, Italy) elemental analyser connected to a Delta V Plus isotope-ratio mass spectrometer via a ConFlo III split interface (both of Thermo Fisher Scientific, Bremen, Germany) operated under continuous helium flow (Spangenberg and Herlec, 2006). The calibration and normalisation of the measured ¹³C to the VPDB scale was performed with international reference materials and UNIL in-house standards (Spangenberg and Herlec, 2006; Spangenberg, 2016). The repeatability and intermediate precision were better than 0.1 ‰ for δ¹³C_org.

4 Results

4.1 Stratigraphy and sedimentology

Belsué (BS) stratigraphic succession records the interfingering between prodelta (Arguís Fm.) and deltaic sediments (Belsué Fm.). The Arguís Fm. are highly bioturbated marls and silts, often rich in glauconite, with sparse bioclasts (e.g. bivalves) and oxidised OM fragments. Sandstone beds (Belsué Fm.) are interlayered within the marls, forming two major coarsening and thickening upward sequences that consist of medium sandstone beds (5–10 m thick) with sharp erosion bases, parallel stratification, undifferentiated ripples, and glauconite-rich horizons (Fig. 4). The Arguís marls are interpreted as prodelta deposits in a poorly circulated and relatively deep marine environment (Millán et al., 1994). The marls prelude deltaic mouth bars (Belsué Fm.) where the fluvial component predominates, although local effects from storms and tides are observed (Millán et al., 1994; Castelltort et al., 2003). Calculated palaeocurrents show a corrected east–southeast sediment supply source, in agreement with previous studies (Puigdefàbregas, 1975; Lafont, 1994; Millán et al., 1994; Garcés et al., 2014). Both formations are interpreted as a mixed delta-carbonate ramp system prograding westward, spanning from Bartonian to Priabonian (Castelltort et al., 2003).

Figure 4 Stratigraphic logs of Belsué – East (BS-E), Belsué – West (BS-O), and Yebra de Basa, with a more complete high-resolution log (Yebra de Basa HR). Red lines represent flooding surfaces (FS), and green lines represent the maximum flooding surfaces (MFS) correlation. Facies interpretation of the sedimentary logs is represented by grey bars, being more proximal the grey bar and more distal the white. Abbreviations used are as follows: Magneto. for magnetostratigraphy, Lithostrat. for lithostratigraphy, and Seq. Strat. for sequence stratigraphy. The poor-exposure zones in Belsué are covered with a semi-transparent white rectangle with a black cross.

On the northern margin of the basin, the Yebra de Basa (YB) section starts with laminated blue marls (Larrés Fm.) that are interlayered by sparse siltstone beds and two dark levels rich in OM (56–60 m in Yebra-HR). The Larrés Fm. transitions to the Sabiñánigo sandstone (SS), composed of two thickening and coarsening upward sequences. The sandstone beds present planar and, through cross-stratification, erosion scours, sigmoidal beds, and flasser-wavy stratification. The upper boundary (ca. 205 m in Yebra HR) is marked by a sharp contrast towards a highly bioturbated and fossiliferous horizon (hard ground), leading to the deposition of laminated grey-blue marls (Pamplona Fm.), less interlayered with siltstones beds but richer in fossiliferous horizons (Turritella sp. mainly). The system is interpreted as a deltaic siliciclastic shelf prograding W–SW (Roigé, 2018), defined as a fluvial delta with local tidal rework. The deltaic system ended abruptly with a major flooding that formed a hard ground and led to the Pamplona marls deposition (Lafont, 1994; Puigdefàbregas, 1975).

As mentioned in Sect. 3.1, although the exposure conditions were excellent for sampling, we had four intervals without continuous sampling, leading to local data gaps. Three of these data gaps were due to the dominance of sandy facies. In YB, the sandy intervals correspond to the Sabiñánigo sandstone deltaic bodies located at approximately 100–120 and 180–200 m (YB-HR section; Fig. 4). In BS, the Belsué sandstone interval is placed between 55 and 60 m (Belsué – E section; Fig. 4). The fourth data gap located at 85–115 m in Belsué – E is the result of a lack of sufficient exposure within the marls and the presence of a coarse-grained sandy interval (Fig. 4).

Facies associations were described by combining the observations in the field and information from previous studies in the Jaca basin (Millán et al., 1994, 2000; Castelltort et al., 2003; Lafont, 1994; Puigdefàbregas, 1975; Boya, 2018). Using the vertical variations of facies in the studied sections, we defined the depositional sequences that record the shoreline progradation and retrogradation (P–R) cycles. Here, we used the smallest correlatable sequences, which are termed parasequences when bounded by the two shallowest facies (flooding surface, FS; Van Wagoner, 1998; Van Wagoner et al., 1990), or genetic units when bounded by the two deepest facies (maximum flooding surface, MFS; Homewood et al., 1992). The sequence stratigraphic interpretation is summarised in Fig. 4, where parasequence thickness from Belsué and Yebra de Basa varies from a few to tens of metres (5–50 m), and its stacking pattern defines two main P–R cycles in both sections.

4.2 Clay mineralogy

In both sections, Belsué and Yebra de Basa, more than 90 % of the total recorded clay mineral assemblage correspond to the sum of chlorite, chlorite–smectite (CS), mica, or illite–smectite (IS) (Fig. 5). This association of clay minerals is characteristic of dominant physical erosion (Adatte et al., 2000). Mica is the most common clay in both sections (40 %–65 %), followed by chlorite (10 %–36 %). In contrast, the percentage of kaolinite is very low (<5 %). In Belsué, both progradations show different chlorite contents, being higher in the upper part. At Yebra de Basa, we observe an increase in mica that coincides with the OM peak. The absence of smectite indicates it has been transformed into CS or IS mixed layers during diagenesis. Its percentage, 20 %–30 % in the studied sections, can be used as a burial estimation (Kübler and Jaboyedoff, 2000). Kaolinite content positively correlates with the deltaic progradation, indicating that kaolinite could be mainly transported.

Figure 5 Schematic stratigraphic log of Belsué – E (a) and Yebra de Basa HR (b) with clay mineral assemblages, IS (illite–smectite), CS (chlorite–smectite), mica, chlorite, and kaolinite. Highlighted in pale grey is the OM-rich interval in Yebra de Basa, and in pale yellow are the main sandstone levels. The progradation–retrogradation cycles (P–R) are drawn with grey and white triangles. The dashed lines represent non-sampled intervals.

4.3 Geochemistry

4.3.1 Major and trace elements

The major and trace elements have been normalised to aluminium (Al) to limit the dilution effect caused by different proportions of terrigenous sediment components (van der Weijden, 2002; Fig. 6). At BS, two increasing pulses of detrital major (Si, Fe, K, and Ti) and trace elements (Mn and Sr) are concomitant with both deltaic progradations (Fig. 6). The similar trend of calcium (Ca) suggests an extra-basinal origin, likely from the eroded Mesozoic or Palaeocene carbonate platforms. Only K shows a negative trend compared to the detrital elements, likely related to clay abundance. At YB, the high TOC interval (depicted in red in Fig. 6) coincides with a relative decrease of the major Si, Ca, Ti, and K and the trace Sr, Zr, and Sn, normalised to Al. In contrast, the OM-related element increase, such as the $\mathrm{V}/\mathrm{Cr}$ ratio, which is related to OM-rich and suboxic–anoxic conditions (van der Weijden et al., 2006), or the $\mathrm{Ni}/\mathrm{Co}$ ratio, which is related to biogenic production (Tribovillard et al., 2006).

Figure 6 Stratigraphic profiles of Belsué – E and Yebra de Basa HR with the normalised major element concentrations (Si, Fe, Mg, Ti, and K), trace element concentrations (Sr and Zn), and trace element ratios ($\mathrm{Ni}/\mathrm{Co}$, $\mathrm{V}/\mathrm{Cr}$, and $\mathrm{U}/\mathrm{Th}$). The progradation–retrogradation cycles (P–R) are drawn with grey and white triangles. They are highlighted: in pale yellow, the sandstone levels, and in pale grey, the OM-rich interval in Yebra de Basa. The dashed lines represent non-sample intervals.

Figure 7 Stratigraphic logs of Belsué and Yebra de Basa including the results of total organic carbon (TOC), carbonate stable isotopes (δ¹⁸O_carb and δ¹³C_carb), and organic carbon stable isotopes (δ¹³C_org). The two progradation–retrogradation cycles referred to in the text are drawn next to the stratigraphy; note the scale change. Highlighted in grey is the MECO isotopic signature, and highlighted in pale red is the OM-rich interval coeval to the MECO thermal peak. The coloured broad lines correspond to the three-point moving average, whilst the central part of the Yebra de Basa section has a seven-point moving-average curve due to the high sampling resolution. Magnetostratigraphic logs from Vinyoles et al. (2021) and Garcés et al. (2014). The dashed lines represent non-sample intervals.

4.3.2 Organic matter content, type, and evaluation

The total organic carbon content (TOC) is low in both sections (average < 0.5 wt %; Fig. 7). At BS, the TOC values range from 0.06 % to 0.38 % (average 0.2±0.07 wt %). The lower one shows a decreasing trend of TOC that follows the deltaic progradation, showing that the OM concentrations decrease with increasing clastic material. The upper section also depicts this trend, with the higher and more stable TOC values (∼0.3 wt %) within marly prodelta deposits. Conversely, YB TOC values range from 0.14 % to 1.3 % (average 0.3±0.13 wt %). The most prominent feature is a dark-grey marl interval with a TOC spike of >1 wt %, which is associated with negative carbonate carbon and oxygen isotope excursions (Fig. 7). Apart from this significant excursion, most TOC values are close to 0.3 wt % and demonstrate no other significant variation along the section. In the high-resolution part of the section, there are small oscillations varying up to ±0.1 wt %.

Figure 8 (a) The scatter plot of hydrogen index (HI) vs. T_max values allows one to discriminate between the different kerogen types, (B), whereas the scatter plot of hydrogen (HI) vs. oxygen index (OI) allows one to assess the OM source and thermal maturity. Here we display the kerogen types (a) and the OM origin (b) from Yebra de Basa (orange) and Belsué (green). Samples with TOC values lower than 0.2 wt % have been excluded. Reference lines for kerogen types and maturity after Espitalié (1986).

The T_max values and the HI (hydrogen index) $/$ OI (oxygen index) ratios serve to classify the OM in terms of type (origin) and thermal maturity (Espitalié et al., 1985). The T_max values of the samples range between 422 and 445 ^∘C (Fig. 8), with some exceptions reaching almost 460 ^∘C. This indicates that the character of the preserved OM is generally immature or within the oil window (Fig. 8). The HI values in YB and BS are generally <100 mg HC g⁻¹ TOC (average 65 mg HC g⁻¹ TOC), corresponding to OM of types III and IV, indicative of a high input of terrestrial plants (Espitalié et al., 1985, Fig. 8). Some samples in Belsué (BS-W) have HI > 150 mg HC g⁻¹ TOC, which is probably related to a slightly different depositional condition between the sections (more distal in the W). The OI values show a similar trend to the HI values, having values below 100 mg CO₂ g⁻¹ TOC (average 82 mg CO₂ g⁻¹ TOC) but being more dispersed than HI values. In summary, the Rock-Eval parameters indicate a recycled source and/or terrestrial origin of the organic matter in both sections (YB and BS).

4.3.3 Carbonate carbon and oxygen stable isotopes

At BS, the δ¹³C_carb values range from −2.6 ‰ to −0.2 ‰ ($-\mathrm{1.5}\pm \mathrm{0.4}$ ‰) and the δ¹⁸O_carb values from −5.7 ‰ to −2.9 ‰ (average $-\mathrm{4.4}\pm \mathrm{0.6}$ ‰; Fig. 7). The δ¹⁸O_carb values show a gradual decrease (of −1.5 ‰) during the first two deltaic progradations, coinciding with a very gradual and small positive trend of the δ¹³C_carb values. After the second progradation, δ¹⁸O_carb rapidly returns to more positive values, which are maintained until the top of the section. This steadiness is not followed by the δ¹³C_carb values, which record a pronounced positive shift of +1 ‰ between 125 and 150 m height. The δ¹⁸O_carb values show a gradually decreasing trend during the first two deltaic progradations (−1.5 ‰), coinciding with a gradual and small positive trend of the δ¹³C_carb values. After the second progradation, δ¹⁸O_carb rapidly returns to more positive values, maintained until the top of the section. This steadiness is not followed by the δ¹³C_carb results that record a pronounced positive shift between 125 and 150 m height (of +1 ‰). At YB, the δ¹⁸O_carb values range from −6.4 ‰ to −4.2 ‰ (average $-\mathrm{4.9}\pm \mathrm{0.4}$ ‰), and the δ¹³C_carb values range from −1.8 ‰ to 0.6 ‰ ($-\mathrm{0.8}\pm \mathrm{0.4}$ ‰; Fig. 7). Small oscillations (±0.5 ‰) of the δ¹⁸O_carb dominate the lower part of the section and end with a significant negative shift at 285 m. There, the δ¹³C_carb values decrease by 0.8 ‰, and the δ¹⁸O_carb values decrease by 1.3 ‰. This level is OM rich (1.0 wt %–1.5 wt % TOC) and also shows a decrease of 2 ‰ in the δ¹³C_org values (see below). Above the OM-rich interval, the δ¹³C_carb and δ¹⁸O_carb values return to pre-event background values, with a shift to higher values towards the base of the Sabiñánigo sandstone, where the δ¹³C_carb values reach the maximum value of 0.6 ‰. A gradual decrease of the δ¹⁸O_carb values coincides with the recurrent progradation events evidenced by the Sabiñánigo sandstone. In contrast, the δ¹³C_carb follows a stable trend around −1 ‰ until the top of the section. Above the Sabiñánigo sandstone and within the Pamplona marls, the δ¹⁸O_carb values show no variance until the top of the section.

4.3.4 Organic carbon isotopes

At BS section, the δ¹³C_org values range from −26.3 ‰ to −22.6 ‰ ($-\mathrm{24.6}\pm \mathrm{0.5}$ ‰; Fig. 6), including two groups of samples with different δ¹³C_org values. The first group is formed by the samples before the first siliciclastic progradation (0–50 m), having relatively low TOC content (0.1 wt %–0.3 wt %) and relatively high δ¹³C_org values (−25 ‰ to −24 ‰). The second group is formed by samples between the deltaic propagations (65–85 m), which have higher TOC content (0.3 wt %–0.5 wt %) and lower δ¹³C_org values ($\sim -\mathrm{26}$ ‰). At YB, the δ¹³C_org values vary between −27.2 ‰ and −23.7 ‰ ($-\mathrm{24.9}\pm \mathrm{0.8}$ ‰; Fig. 6), with the lowest value of −27.2 ‰ having been measured within the OM-rich level at 285 m. The δ¹³C_org values increase upward till the base of the Sabiñánigo sandstone, where they show a negative spike, coinciding with the positive shift of the δ¹³C_carb and δ¹⁸O_carb values. Then the δ¹³C_org values first return to the pre-event background values and then show a negative excursion of 2 ‰ in the last two samples of the Sabiñánigo sandstone.

Figure 9 Oxygen isotope (δ¹⁸O_carb) correlation panel for the studied sections (Belsué and Yebra de Basa) with MECO target curves from Alano (Italy, Tethys Ocean; Spofforth et al., 2010), ODPS 1051 (N Atlantic Ocean; Edgar et al., 2010), ODPS 702 (S Atlantic Ocean; Bohaty et al., 2009), and ODPS 738 (S Indic Ocean; Bohaty et al., 2009). Data from the bulk- and fine-sediments fractions. Highlighted in red is the OM-rich interval (TOC peak) in Yebra de Basa. The two progradation–retrogradation cycles referred to in the text are drawn with grey and white triangles. The data are scaled according to magnetostratigraphic tie points, between C18r–18n.2n and C18n.2n–C18n.1r chrons.

5 Discussion

5.1 MECO isotopic record

In summary, at Belsué, the oxygen isotope record shows a general trend towards more negative values from the base to the middle sandstone units. This trend peaked with a negative δ¹⁸O_carb shift of ∼1 ‰ before the sandstone unit, just in the chron transition C19r–C18n.2n (Fig. 9). The δ¹⁸O_carb values in Yebra de Basa show a small-scale variability consistent with local effects. Then a negative shift of ∼1.2 ‰ occurs close to the deltaic progradation (Fig. 9). The MECO zenith around the magnetic reversal C19r–C18n.2n (Bohaty et al., 2009; Edgar et al., 2010; Henehan et al., 2020) is represented by the progradation of the deltaic facies over the prodelta, i.e. the deltaic facies of Belsué Fm. in Belsué and the Sabiñánigo sandstones in Yebra de Basa (Fig. 9). No isotopic data were obtained for this interval. Nevertheless, the onset of the main thermal event, just before the sandstone occurrence, is preserved as the negative excursion. After the sandstone progradation, the δ¹⁸O_carb values in both sections return closely to those before the excursion (Fig. 9).

The trend observed in the studied sections of the South Pyrenean basin (SPFB) is shared with most high-resolution offshore and nearshore isotopic records of the MECO (Bohaty et al., 2009; Bohaty and Zachos, 2003; Edgar et al., 2010, 2020; Spofforth et al., 2010; Jovane et al., 2007; Giorgioni et al., 2019; Galazzo et al., 2014). They all show the same trend towards more negative δ¹⁸O_carb values, which is intensified during the MECO peak (Fig. 9). The end of the event is marked in both sections by a rapid increase of the δ¹⁸O_carb values by ∼1 ‰, as reported in other sections worldwide (Bohaty et al., 2009; Galazzo et al., 2014; Edgar et al., 2010, 2020; Giorgioni et al., 2019; Spofforth et al., 2010).

Contrarily to the agreement between our new data and most of the available oxygen isotope records, the δ¹³C_carb results do not show a clear correlation with global curves (Fig. 7). On one hand, the Belsué section records a positive δ¹³C_carb excursion, with a delay with respect to the δ¹⁸O_carb minimum. The Yebra de Basa section shows a prominent positive δ¹³C_carb excursion just before the main deltaic progradation (320–340 m from YB). On the other hand, most of the oceanic geochemical records show a small carbon isotope excursion (CIE) at the MECO peak of warming (∼40 Ma; Westerhold and Röhl, 2013; Bohaty et al., 2009; Spofforth et al., 2010), but before and after the δ¹³C_carb values are highly variable, showing opposite trends between hemispheres (Henehan et al., 2020; Giorgioni et al., 2019). This CIE, like in other Northern Hemisphere sites (Sppofforth et al., 2010; Giorgioni et al., 2019), is not well represented in our data. Therefore, our δ¹³C_carb data seem to confirm that the MECO is not associated with an extensive input of depleted ¹³C in the environment, as suggested in another study (Henehan et al., 2020). However, an alternative explanation for this discrepancy would be that the sandstone progradation masked the CIE.

Instead of a global origin, however, the δ¹³C_carb variations could be caused by local changes in the carbon isotope composition of the dissolved inorganic carbon (DIC), pH, rate and temperature of carbonate precipitation, and its mineralogy (e.g. Swart, 1995). Here, given the proximity of the continent in the shelf environment of the studies section, we cannot rule out that spatial and temporal variations in freshwater input of different components (i.e. riverine, estuarine, and groundwater sources) could have altered the isotopic composition of the DIC and the carbonate δ¹³C and δ¹⁸O values (e.g. Marshall, 1992; Saltzman and Thomas, 2012; Wendler, 2013; Läuchli et al., 2021; Castelltort et al., 2017). This freshwater input could produce carbonate depleted in ¹³C and ¹⁸O (see detailed discussion in Sect. 5.2). The fact that the Belsué and Yebra de Basa isotopic records preserve the MECO excursion in δ¹⁸O_carb suggests that the δ¹³C_carb values could also record a global signal. However, this is difficult to appreciate because of the small δ¹³C_carb variations in the studied sections and the somewhat variable and peculiar published carbon isotope record during the MECO. Therefore, the correlations were based on the δ¹⁸O_carb records (Fig. 9). In summary, our results record a decoupling between oxygen and carbon isotopes. The δ¹⁸O_carb values seem to follow the global trend, while the variation of the δ¹³C_carb values remains ambiguous. The particular carbon isotope record during the MECO, with variable and opposite trends between hemispheres (Henehan et al., 2020) and a brief negative carbon isotope excursion recorded just at some sites (Westerhold and Röhl, 2013; Bohaty et al., 2009; Spofforth et al., 2010), could be an explanation. However, given the location of the studied sections on a continental shelf, it is important to check the possible diagenetic influence or alteration of the primary isotopic signal.

5.2 Primary versus diagenetic signals

The carbonate primary carbon and oxygen isotope compositions may be affected by post-depositional processes, including the neoformation of authigenic and diagenetic phases. Therefore, before the palaeoenvironmental interpretation of δ¹³C_carb and δ¹⁸O_carb records from shallow marine environments, it is necessary to determine primary versus diagenetic signal components. This discrimination requires an understanding of the factors controlling the primary marine isotopic composition and an evaluation of potential diagenetic overprints on the original geochemical signatures (e.g. Marshall, 1992; Schrag et al., 1995).

Oxygen isotopes in carbonates are controlled by the temperature of formation, the δ¹⁸O value of the carbonate-precipitating fluid (δ¹⁸O_w), the mineralogy (e.g. higher δ¹⁸O in dolomite vs. calcite), and any environmental parameter (e.g. pH, salinity) affecting the rate of carbonate precipitation (Swart, 1995). The effect of diagenetic alteration is more pronounced in the case of oxygen isotopes than in carbon isotopes due to the high amount of oxygen relative to carbon present in post-depositional fluids and their variable δ¹⁸O values (e.g. Marshall, 1992; Schrag et al., 1995; Fio et al., 2010). Carbonate with low δ¹⁸O values can be produced by increasing temperature, freshwater input, and meteoric diagenesis, whereas ¹⁸O enrichment could indicate either lower temperature or evaporation (e.g. Marshall, 1992; Patterson and Walter, 1994; Schrag et al., 1995). In contrast, carbon isotopes are not thought to be directly influenced by temperature and are generally more resistant to diagenetic processes (Patterson and Walter, 1994; Schrag et al., 1995; Swart, 1995). However, δ¹³C values are also controlled by kinetic effects, mineralogy, and mainly by the δ¹³C value from the DIC (Wendler, 2013). The primary diagenetic process that affects the δ¹³C values of the DIC is the oxidation of the organic matter, which produces CO₂ (and DIC species) depleted in ¹³C (low δ¹³C values). Therefore, the δ¹³C values of the DIC and derived carbonates indicate the source of carbon, including the type of degraded or oxidised organic matter (OM) of different types, original seawater carbon, and skeletal and non-skeletal carbonate sources (e.g. Swart, 1995). In proximal depositional environments, however, the δ¹³C values could be modified by (1) OM source, productivity, and burial rate; (2) extra-basinal carbonate input; (3) water circulation, stratification, and evaporation; and (4) terrestrial runoff and weathering (Saltzman and Thomas, 2012; Läuchli et al., 2021). Considering this, δ¹³C is usually used as a global correlation tool, since it can register eustatic sea level fluctuations, changes in weathering flux, or significant perturbations in the global carbon cycle (e.g. volcanic CO₂ input; Wendler, 2013, and references therein).

Figure 10 δ¹⁸O_carb–δ¹³C_carb scatterplot of all Belsué and Yebra de Basa samples. The size of the symbols is small for samples with TOC < 0.2 wt % and big for samples with TOC > 0.2 wt %. The Pearson correlation coefficients (r) and regression lines are shown.

The degree of diagenetic alteration was assessed through three different approaches. First, we evaluated the relationship between δ¹³C and δ¹⁸O values (Brasier et al., 1996). Statistically, a non-significant correlation (Pearson correlation coefficient r<0.6) indicates that a diagenetic overprint of the primary isotopic signature can be excluded (e.g. Fio et al., 2010). In both sections, no statistically significant correlation (r<0.3) was found between the δ¹⁸O_carb and δ¹³C_carb values. This lack of relationship suggests that no or minor diagenetic modifications affected the primary isotopic compositions (Fig. 10). The second approach used to assess the degree of alteration uses clay mineralogy. Kübler and Jaboyedoff (2000) defined four diagenetic zones by comparing illite crystallinity with mineral assemblages and organic-matter type. The Belsué and Yebra de Basa samples have 20 %–30 % smectite within the illite–smectite (IS) mixed layers and are within the third diagenetic zone of Kübler and Jabeyedoff (2000), i.e. shallow diagenesis (ca. 60–80 ^∘C). Another diagenetic indicator is the maximum temperature (T_max) reached during the Rock-Eval pyrolysis (S2), which marks the maturity of the OM. The T_max values obtained in samples with relatively high OM content (TOC > 0.5 wt %; S2 > 0.2) were <440 ^∘C (Fig. 8), which corresponds to the beginning of the oil window (ca. 60 ^∘C; Espitalié et al., 1985). This maturity level of the organic matter agrees with vitrinite reflectance and Raman measurements in the studied area (Labaume et al., 2016). In summary, the three approaches for assessment of the diagenetic degree – i.e. carbonate δ¹³C and δ¹⁸O values, illite crystallinity, and thermal maturation of the organic matter (T_max) – suggest that the diagenetic overprint in the studied Belsué and Yebra de Basa rocks is low. The primary isotopic signal is preserved largely in both sections. It can be safely used to study palaeoenvironmental conditions and can be compared to global key isotopic curves during the MECO event.

5.3 The organic matter peak in Yebra de Basa

In the Yebra de Basa (YB) section, an increase in TOC content at the 280–290 m interval (up to 1.5 wt % TOC) is associated with a negative isotope excursion of −1.5 ‰ for δ¹⁸O_carb, −2.0 ‰ for δ¹³C_org, and −0.8 ‰ for δ¹³C_carb values (Figs. 7 and 9). The OM-rich interval occurs 50 m below the main Sabiñánigo sandstone progradation in Yebra de Basa. It is not coincident with the prominent increase of detrital input, marked by an increase in grain size. This boost in OM burial is also observed in the Neo-Tethys region, like in Italy (Spofforth et al., 2010) and the Crimea–Caucasus (Benyamovsky, 2012), which may have played an important role in carbon drawdown and rapid cooling after the MECO event (Bohaty et al., 2009; Henehan et al., 2020).

Several possibilities could explain the presence of an OM-rich interval before a deltaic progradation. First, a significant freshwater input in a restricted basin can lead to water stratification where anoxic conditions are favoured, increasing OM preservation independently of its source. Nevertheless, the slight increase in redox-sensitive elements (V and Mo, Fig. 6) is too limited to support the development of water stratification and the resulting suboxic–anoxic conditions (Tribovillardet al., 2006). Second, the enhanced freshwater input could have increased nutrient availability and marine productivity. We, however, reject this hypothesis because the geochemical data point to the main terrestrial components of the organic matter (low HI–OI) and little to no nutrient availability increases (low Ni concentration; Tribovillard et al., 2006). Therefore, the most probable explanation is that the OM peak could be related to a significant increase in detrital input and terrestrial OM. The presence of several westward dark-marl beds suggests that it was not a unique episode but a series of recurrent events (Boya, 2018). In addition, the terrestrial origin is also supported by the strong correlation (r>0.7) observed between the siliciclastic elements (Al, Ti, Fe) and the TOC or all the OM-related trace elements (V, Mo, Ba, and Th; Tribovillard et al., 2006). Despite this, the isotopic results do not agree with this correlation because pre-Miocene marine OM had lower δ¹³C than terrestrial OM (Popp et al., 1989). Thus, an alternative explanation for the negative δ¹³C_org and δ¹³C_carb excursion may be an increased input of organic matter released from soils containing bacterial biomass with low δ¹³C_org values (Fio et al., 2010). This agrees with the Rock-Eval results, which point towards a terrestrial origin of this organic matter (Fig. 8).

As a result, the Sabiñánigo sandstone represents a singular deltaic event embedded in long-lasting prodelta conditions (Vinyoles et al., 2021) in which no evident organic events occur. Therefore, we interpret the occurrence of the OM-rich level just before the Sabiñánigo sandstone as the first indicator of a shift towards a setting with more fluvial conditions, this being the first evidence of the main MECO excursion in the region.

5.4 MECO response in the South Pyrenean foreland basin

The integration of available age constraints (Garcés et al., 2014; Vinyoles et al., 2021) and the new high-resolution isotopic record show that MECO's warming peak (∼40 Ma) is associated with isochronous progradation, which can be followed all along the SPFB source-to-sink system (Fig. 11; Vinyoles et al., 2021). In the Tremp–Graus basin, the Escanilla fluvial system was fed by the Sis–Gurp and Pobla alluvial systems, where a grain size increase is recorded at ca. 40 Ma (Whittaker et al., 2011). Downstream, in the time-equivalent sections in the Ainsa basin, an anomalous amalgamated Olsón sheet stands out from the landscape as a continuous and thick conglomeratic bed, interpreted as a stacking of several braided river channels (Fig. 11; Verité, 2019; Labourdette, 2011; Puigdefàbregas, 1975; Vinyoles et al., 2021). In the deltaic counterparts (Jaca basin), a significant progradation of deltaic deposits on top of slope marls is observed in the studied sections (BS and YB; Lafont, 1994; Puigdefàbregas, 1975; Vinyoles et al., 2021). Finally, in the deeper sink environments of the Jaca and Pamplona basins, the correlation with the turbiditic systems is still debated and needs further research.

Figure 11 Correlation panel between Belsué, Yebra de Basa, and Olsón section with the GPTS 2016 (Ogg et al., 2016). The stratigraphic sections are modified from Garcés et al. (2014) and Vinyoles et al. (2021). The oxygen isotopic (δ¹⁸O_carb) record from ODPS 738 and the MECO age constraints defined by yellow and blue bars are modified from Bohaty et al. (2009). The oxygen isotopic record (δ¹⁸O_carb) from Belsué and Yebra de Basa corresponds to our results, and the dashed lines represent non-sample intervals. The sedimentation rates (SRs) from Belsué, Yebra de Basa, and Olsón are the average SR between chron C18n.2n and C17r, as from the data of Vinyoles et al. (2021).

Previous works (Puigdefàbregas, 1975; Lafont, 1994) interpreted these deltaic sequences as eustatic fluctuations of the relative sea level, which can be related to different possibilities, such as thermal expansion or glacioeustasy. Ephemeral ice sheets in Antarctica during the Middle Eocene are likely, and it seems plausible that the progressive shift towards ice-house conditions could have significant implications during the MECO (Edgar et al., 2007; Huyghe et al., 2012; Baatsen et al., 2020). However, considering the temperature increase interpreted during the MECO zenith (+4 to 6 ^∘C; Bohaty et al., 2009), we should expect a sea level rise (ice caps melting and thermal expansion) instead of the observed regression and system progradation.

Alternatively, an abrupt sediment supply increase can also explain a progradation of deltaic systems. Several studies observed that the main Palaeogene hyperthermals are often associated with an enhanced flux of terrigenous material, interpreted as a boost of the hydrological cycle and higher seasonality (Schmitz and Pujalte, 2007; Chen et al., 2018; Foreman and Straub, 2017; Pujalte et al., 2015; Schmitz and Pujalte, 2007). Although the MECO is not an abrupt event like other hyperthermals, instead being a more extended period of gradual warming (ca. 500 kyr; Bohaty et al., 2009), we also observe this progradation being focused during the warming peak (ca. 40 Ma). Accordingly, an explanation for the progradation is that the MECO prolonged warming produced an enhanced hydrological cycle that favoured sediment production and transport, thus leading to an increase in sediment supply and favouring the system progradation at the peak of the event. The nature of a greater sediment provision (Qs) should be considered to have originated upstream, for instance, being linked to enhanced sediment remobilisation (e.g. floodplain) or accelerated hillslope processes (Foreman et al., 2012).

Therefore, the coincidence in time of a basin-wide progradation in the SPFB and the MECO might implicate a link between them. Our geochemistry analyses also suggest a terrestrial origin for this OM, which points towards an increase in soil remobilisation, erosion, and transport in continental environments during the MECO event.

5.5 Global implications and correlation

The global impact of the MECO event in continental settings remains poorly documented, with only a few studies in continental environments having been performed around the globe (e.g. Bosboom et al., 2014; Mulch et al., 2015). In the North American plateau, a boost of precipitation during the MECO is derived from lower δ¹⁸O_carb values (Mulch et al., 2015). In contrast, in the Tarim basin (China), a shift towards arid conditions has been interpreted from a reduction in fern palynomorphs (Bosboom et al., 2014). This aridification trend in central Asia differs from the documented Neo-Tethys Ocean dynamic, where marine records show an increase in organic-matter (OM) burial during the MECO peak and part of the post-MECO recovery (Spofforth et al., 2010; Giorgioni et al., 2019 Benyamovskiy, 2012). Increased sediment supply due to enhanced erosion and transport provides a mechanism for the more efficient burial of OM during this and other hyperthermals (Galy et al., 2007). If this enhanced OM burial is global or sufficiently widespread (it is absent in several sections, including Belsué in this study), it could represent an important mechanism to explain the carbonate δ¹³C increase that is recorded globally during the post-event recovery and the associated rapid return to pre-event conditions, maybe playing an essential role in the drawdown of atmospheric carbon (e.g. Bohaty et al., 2009; Henehan et al., 2020; Sluijs et al., 2013; Edgar et al., 2020; Giorgioni et al., 2019; Spofforth et al., 2010).

Considering the long duration of the MECO event (ca. 500 kyr; Bohaty et al., 2009), some of the most important effects in the ocean occur during its peak phase, e.g. ocean acidification (Bohaty et al., 2009; Henehan et al., 2020; Arimoto et al., 2020) or OM burial (Giorgioni et al., 2019; Spofforth et al., 2010). In the SPFB, the continental progradation also occurred at the end of the event, supported by the sedimentological and geochemical evidence that shows an increase of sediments delivered to the sea, including large amounts of organic matter of terrestrial origin. Hence, our work suggests a link between enhanced hydrological cycles and enhanced OM transport and burial, possibly accounting for the observations of enhanced OM burial around the Neo-Tethys region. This response in sediment delivery rate (OM burial in shallow and restricted basins) has been previously documented for other early Eocene hyperthermals (Chen et al., 2018; Foreman et al., 2012; Pujalte et al., 2015; Foreman and Straub, 2017; Honegger et al., 2020). Hence, despite its important differences in relation to the early Eocene hyperthermals, the MECO shares several attributes with them around the warming peak. In summary, our results point to a more intense hydrological cycle perturbing rainfall patterns in the Pyrenean region during the MECO peak and leading to increased sediment supply, expressed by a major progradation of sedimentary systems and, eventually, an increase in OM burial in the nearby oceanic basins.

6 Conclusions

In the South-Pyrenean foreland basin, an important progradation affected the entire sediment-routing system from fluvial to deltaic environments at the time of the Middle Eocene Climatic Optimum (MECO). Here we present a new high-resolution multiproxy dataset, including stable isotopes, Rock-Eval, XRF, and clay minerals, covering the different MECO phases from two well-dated key sections. The new stable-isotope records from Belsué (BS) and Yebra de Basa (YB) sections show a significant negative shift in the shallow marine sediments around the main warming peak of the MECO event, reported for the first time in the Pyrenean region. In Yebra de Basa, an organic-rich interval of terrestrial origin is found before the main deltaic progradation. It is associated with a negative excursion in oxygen and carbon isotopes. The correlation between the MECO and the basin-wide progradation, as well as the new geochemical results, presents compelling evidence for a climatic driver, suggesting an enhanced hydrological cycle in the Pyrenean region that caused a boost in sediment and carbon export. This is in agreement with previous studies from the Neo-Tethys Ocean that recorded an increase in organic-matter burial during the peak of the MECO and early post-MECO.

Although the duration of the MECO and its isotopic signature differ with respect to early Eocene hyperthermals (e.g. PETM), there are similarities around the warming peak that trigger a comparable response, including ocean acidification, organic-matter burial, or a boost in sediment supply export from land to sea. Nevertheless, further work is needed to understand the role of potential sediment supply increase from the proximal continental environments towards the deeper oceanic basins and, importantly, to quantify sediment and organic export and its relationship with carbon burial and silicate weathering.

Our results support the view that high-accommodation settings in foreland basins are important recorders of palaeoenvironmental signals, even in shallow marine environments. Although certainly noisy, the fact that climate signals are preserved in these settings provides a range of potentially expanded sections that can be an interesting complement to high-resolution but more condensed deep-sea palaeoclimatic records, particularly during high-CO₂ globally warm episodes of the Earth's history when the carbonate-rich oceanic records may undergo intervals of non-deposition or dissolution.