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 15 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 20 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 (


Introduction
The Middle Eocene Climatic Optimum (MECO) is a global warming event that occurred during the  Ma) and briefly reversed the longer-term cooling trend of the middle to upper Eocene ( Fig. 1; Arimoto et al., 2020;Bijl et al., 35 2010;Bohaty and Zachos, 2003;Bohaty et al., 2009;Bosboom et al., 2014;Galazzo et al., 2014;Henehan et al., 2020;Edgar et al., 2010Edgar et al., , 2020Giorgioni 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., 2020). Marine bulk and benthic oxygen isotope compositions ( 18 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 ( 13 C values), unlike earlier hyperthermals of the Cenozoic 40 (e.g., Palaeocene-Eocene Thermal Maximum PETM, Eocene Thermal Maximum ETM 2 among others), 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 (pCO2) during the warming peak (Henehan et al., 2020;Bijl et al., 2010), and numerous potential CO2 sources have been proposed. Among them, a prolonged pulse of metamorphic decarbonization possibly linked 45 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. 2020), or lower continental weatherability (van der Ploeg et al. 2018). However, the pCO2 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 CO2 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, 50 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. 55 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). 60 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). In some locations, the warmer conditions reduced nutrient availability, decreasing the benthic productivity (Arimoto et al., 2020;Bijl et al., 2010, Galazzo et al., 2014Moebius et al., 2015). In contrast, the Southern Ocean (Moebius et al., 2014) or the Neo-Tethys Ocean (Galazzo et al., 2013) record 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. 75 (2015) suggested a boost of precipitation in the North American plateau derived from low δ 18 O values, while Bosboom et al. (2014) documented a shift towards arid conditions in the Tarim Basin with a reduction in fern palynomorphs. Such 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 nutrients 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 80 response of weathering, erosion, and sediment transport in terrestrial systems to global warming, as it has also been posed for other hyperthermals recently (e.g., Chen et al., 2018;Foreman et al., 2012Foreman et al., , 2017Honegger 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 85 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 geochemical and isotopic evolution in the marine domain. We focus on two separated 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 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;130 Puigdefàbregas, 1975;Puigdefàbregas and Souquet, 1986;Millán et al., 1994Millán et al., , 2000Castelltort 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 135 Arguís marls and the Belsué sandstones, while the Gabardiella and Pico del Águila anticlines were growing. Different authors studied the influence in the stratigraphy of local tectonic movement in 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 paleobathymetry in the BS area 140 (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 150 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.

Field and sampling 155
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 (BS-E) was sampled at a resolution of 0.5-1.0 m and the upper one (BS-O) every 3-6 m. In the northern margin, the higher average 160 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 (YB-HR) and each 9-15 m in the other intervals (YB-C and 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 165 geochemical analyses in the laboratories of the University of Geneva and the University of Lausanne. The sample surface was cleaned with deionized water, the weathered material was removed, and then dried at 45ºC for 2-3 days. 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 170 because of the fine-grained nature of the marls (e.g., Quaternary cover).

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 (1983Kübler ( , 1987 and Adatte et al. (1996). 175 Ground chips were mixed with deionized water (pH 7-8) and agitated. The carbonate fraction was removed by treatment with 10% HCl at room temperature and then for 20 minutes or more until all carbonate was dissolved. The insoluble residue was disaggregated (ultrasonication, 3 min), washed and centrifuged (8 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. XRD analyses of oriented clay samples were 180 made after air drying at room temperature at 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 µm and 2-16 µm).

Major and trace element compositions
Major and trace element concentrations from 24 representative samples per section were determined by X-ray fluorescence 185 (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 one night to obtain the loss on ignition (LOI) value. Then, 1.2000 ± 0.0005 g of the calcinated sample was mix with 6.0000 ± 0.0005 g lithium tetraborate (Li2B4O7) to prepare the fused glass disc using an automated glass bead-casting machine (Pearl-X'3) at the University of Geneva. 190 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.0005g of non-calcinated sample powder. The mixture was homogenised and pressed (1 ton hydraulic press, 20 seconds). 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.

Rock-Eval pyrolysis 195
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 160000 standard. Aliquots of samples were placed in an oven, heated at 300 o C under an inert atmosphere, and then gradually pyrolyzed 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 CO2 and hydrocarbon (HC) concentration during the entire process. The calculated parameters included total organic 200 carbon content (TOC, wt.%), hydrogen index (HI in mg HC/g TOC), oxygen index (OI in mg CO2/g TOC), and Tmax (ºC) according to Espitalié et al. (1985) and Behar et al. (2001).

Carbonate carbon and oxygen stable isotopes
Carbonate carbon and oxygen stable isotope ratios ( 13 Ccarb and 18 Ocarb) of whole rock powders containing > 10 wt.% CaCO3 (n = 237) were determined at the laboratories of the Institute of Earth Surface Dynamics of the University of Lausanne (IDYST-205 UNIL) using a Thermo Fisher Scientific Gas Bench II carbonate preparation device connected to a Delta Plus XL isotope ratio mass spectrometer. The CO2 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 δ = (Rsample -Rstandard)/Rstandard × 1000 and R = 13 C/ 12 C or 18 O/ 16 O. The δ 13 Ccarb and δ 18 Ocarb values were standardized relative to the international VPDB scale by calibration of the reference gases and working standards with the 210 international reference materials NBS 18 (carbonatite, δ 13 C = -5.04 ‰, δ 18 O = -23.00 ‰) and NBS 19 (limestone, δ 13 C = +1.95 ‰, δ 18 O = -2.19 ‰). Analytical uncertainty (1 sigma), monitored by replicate analyses of the international calcite standard NBS 19 and the laboratory standard Carrara Marble (δ 13 C = +2.05 ‰, δ 18 O = -1.7 ‰), was not greater than ±0.05‰ for δ 13 C and ±0.1 ‰ for δ 18 O.

Organic carbon stable isotopes 215
The organic carbon stable isotope ratios ( 13 Corg values in ‰ vs. VPDB) were determined in 155 samples, which were previously decarbonated by treatment with 10% v/v HCl, thoroughly washed with deionized water and dried (40 °C, 48 h).
The 13 Corg measurements were performed at the IDYST-UNIL by elemental analysis/isotope ratio mass spectrometry, using a Carlo Erba 1108 (Fisons Instruments, Milan, Italy) elemental analyzer 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 220 helium flow (Spangenberg and Herlec, 2006). The calibration and normalization of the measured 13 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 13 Corg.

Stratigraphy and sedimentology 225
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 oxidized OM fragments. Sandstone beds (Belsué Fm.) are interlayered within the marls, forming two major coarsening and thickening upwards sequences that consist of medium sandstone beds (5-10 m thick) with sharp erosion base, parallel stratification, undifferentiated ripples, and glauconite rich horizons (Fig. 4). The Arguís marls are interpreted as prodelta 230 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 paleocurrents show a corrected east/south-east 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 235 (Castelltort et al., 2003).  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 upwards sequences. The sandstone beds present planar and through cross-stratification erosion scours, sigmoidal beds, and flasser-wavy stratification. The upper boundary (ca. 245 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é et al., 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 prodelatic marls deposition (Lafont, 1994;Puigdefàbregas, 1975). 250 As mentioned in section 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 approximately at 100-120 m 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, results of lack of sufficient exposure within the marls and the presence 255 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;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 260 bounded by the two shallowest facies (flooding surface, FS;van Wagoner et al., 1988van 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 summarized in Fig. 4, where parasequences thickness from Belsué and Yebra de Basa vary from a few to tens of meters (5-50 m), and its stacking pattern defines two main P-R cycles in both sections.

Clay mineralogy 265
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 270 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, 2000). Kaolinite content positively correlates with the deltaic progradation, indicating that kaolinite could be mainly transported.

Major and trace elements 280
The major and trace elements have been normalized to aluminium (Al) to limit the dilution effect caused by different proportions of terrigenous sediment components (van der Weijden et al., 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 extrabasinal 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 285 interval (depicted in red in Figure 6) coincides with a relative decrease of the major Si, Ca, Ti, and K, and the trace Sr, Zr, and Sn, normalized to Al. In contrast, the OM-related elements increase, such as the V/Cr ratio, which is related to OM-rich and suboxic/anoxic conditions (van der Weijden et al., 2006), or the Ni/Co ratio, which is related to biogenic production (Tribovillard et al., 2006).

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 295 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 300 this significant excursion, most TOC values are close to 0.3% wt.% and 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.%. The Tmax 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 Tmax values of the samples range between 422 and 445 o C (Fig. 8), with some 310 exceptions reaching almost 460 o 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 with 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 CO2/g TOC 315 (average 82 mg CO2/g TOC), but 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).

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 ẟ 18 Ocarb shift of ∼1 ‰ before the sandstone unit, just in the chron transition C19r-C18n.2n (Fig. 9). The ẟ 18 Ocarb 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 355 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 ẟ 18 Ocarb values in both sections return closely to those before the excursion (Fig. 9). 360 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., 2010Edgar et al., , 2020Spofforth et al., 2010;Jovane et al., 2007;Giorgioni et al., 2019;Galazzo et al., 2014). They all show the same trend towards more 370 negative ẟ 18 Ocarb 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 ẟ 18 Ocarb values by ∼1 ‰, as reported in other sections worldwide (Bohaty et al., 2009;Galazzo et al., 2014;Edgar et al., 2010Edgar et al., , 2020Giorgioni et al., 2019;Spofforth et al., 2010).
Contrarily to the agreement between our new data and most of the available oxygen isotope records, the ẟ 13 Ccarb results do not show a clear correlation with global curves (Fig. 7). On one hand, the Belsué section records a positive ẟ 13 Ccarb excursion, with 375 a delay respect the ẟ 18 Ocarb minimum. The Yebra de Basa section shows a prominent positive ẟ 13 Ccarb 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 Myr; Westerhold and Röhl, 2013;Bohaty et al., 2009;hemispheres (Henehan et al., 2020;Giorgioni et al., 2019). This CIE, like in other northern hemisphere sites (Sppofforth et 380 al., 2010;Giorgioni et al., 2019), is not well represented in our data. Therefore, our ẟ 13 Ccarb data seem to confirm that the MECO is not associated with an extensive input of depleted 13 C in the environment, as suggested in a 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 ẟ 13 Ccarb 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, 385 2015). Here, given the proximity of the continent in the shelf environment of the studies section, we cannot rule out that spatial and temporal variation 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 δ 13 C and δ 18 O values (e.g., Marshall, 1992;Saltzman and Thomas, 2012;Wendler, 2013Läuchli et al., 2021. This freshwater input could produce carbonate depleted in 13 C and 18 O (see detailed discussion in section 5.2). The fact that the Belsué and Yebra de Basa isotopic records preserve the MECO 390 excursion in ẟ 18 Ocarb suggests that the ẟ 13 Ccarb values could also record a global signal. However, this is difficult to appreciate because of the small ẟ 13 Ccarb 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 ẟ 18 Ocarb records (Fig. 9). In summary, our results record a decoupling between oxygen and carbon isotopes. The 18 Ocarb values seem to follow the global trend, while the variation of the 13 Ccarb values remains ambiguous. The particular carbon isotope record during the MECO, with variable and opposite 395 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.

Primary versus diagenetic signals 400
The carbonate primary carbon and oxygen isotope compositions may be affected by postdepositional processes, including the neoformation of authigenic and diagenetic phases. Therefore, before the paleoenvironmental interpretation of 13 Ccarb and ẟ 18 Ocarb records from shallow marine environments, it is necessary to determine primary versus diagenetic signal components.
This discrimination requires understanding 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). 405 Oxygen isotopes in carbonates are controlled by the temperature of formation, the 18 O value of the carbonate-precipitating fluid ( 18 Ow), the mineralogy (e.g., higher 18 O in dolomite vs. calcite), and any environmental parameter (e.g., pH, salinity) affecting the rate of carbonate precipitation (Swart, 2015). The effect of diagenetic alteration is more pronounced in the case of oxygen isotopes than carbon isotopes due to the high amount of oxygen relative to carbon present in postdepositional fluids and their variable ẟ 18 O values (e.g., Marshall, 1992;Schrag et al., 1995;Fio et al., 2010). Carbonate with low ẟ 18 O values can 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, 2015). However, 13 C values are also controlled by kinetic effects, mineralogy, and mainly by the 13 C value from the DIC (Wendler, 2013). The primary diagenetic process that affects 415 the ẟ 13 C values of the DIC is the oxidation of the organic matter, which produce CO2 (and DIC species) depleted in 13 C (low ẟ 13 C values). Therefore, the ẟ 13 C values of the DIC and derived carbonates indicate the source of carbon, including the type of degraded/oxidized organic matter (OM) of different types, original seawater carbon, skeletal and non-skeletal carbonate sources (e.g., Swart, 2015). In proximal depositional environments, however, the ẟ 13 C values could be modified by (1) OM source, productivity, and burial rate, (2) extrabasinal carbonate input, (3) water circulation/stratification and evaporation, (4) 420 terrestrial runoff and weathering (Saltzman andThomas, 2012, Läuchli et al., 2021). Considering this, 13 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 CO2 input; Wendler 2013 and references therein). The degree of diagenetic alteration was assessed through three different approaches. First, was evaluated the relationship between 13 C and 18 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 statistical significant correlation (r < 0.3) was found between the 18 Ocarb and 13 Ccarb values. This lack of 430 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 3 rd diagenetic zone of temperature (Tmax) reached during the Rock-Eval Pyrolysis (S2), which marks the maturity of the OM. The Tmax values obtained in samples with relatively high OM content (TOC > 0.5 wt.%; S2 > 0.2) were < 440 o C (Fig. 8), which corresponds to the beginning of the oil window (ca. 60 o 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 13 C and 18 O values, illite crystallinity, and thermal maturation of the 440 organic matter (Tmax), 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 paleoenvironmental conditions and be compared to global key isotopic curves during the MECO event.

The organic matter peak in Yebra de Basa
In 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 445 a negative isotope excursion of -1.5‰ for ẟ 18 Ocarb, -2.0‰ for ẟ 13 Corg and -0.8‰ for ẟ 13 Ccarb values ( Fig. 7 and 9). The OMrich 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 et al., 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., 450 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 (Tribovillard et al., 455 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 main terrestrial components of the organic matter (low HI-OI) and no low nutrient availability increase (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 dark-marl beds westwards suggests it was not a unique episode but a series of recurrent events (Boya, 2018). In 460 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 ẟ 13 C than terrestrial OM (Popp et al., 1989). Thus, an alternative explanation for the negative ẟ 13 Corg and ẟ 13 Ccarb excursion may be an increased input of organic matter released from soils containing bacterial biomass with low ẟ 13 Corg values (Fio et al., 2010). This agrees with the 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, being the first evidence of the main MECO excursion in the region. 470

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., 475 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 et al., 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;Puigdefabregas et al., 1975, Vinyoles et al., 2021. Finally, in the deeper sink environments of the Jaca and 480 Pamplona basins, the correlation with the turbiditic systems is still debated and needs further research. Previous works (Puigdefàbregas, 1975;Lafont, 1994) interpreted these deltaic sequences as eustatic fluctuations of the relative sea level, which can relate 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 icehouse conditions could have significant implications during the MECO (Edgar et al., 2007;Huyghe et al., 2012;Baatsen et al., 2020). However, 485 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 Paleogene hyperthermals are often associated with an enhanced flux of terrigenous material interpreted as a boost of the hydrological cycle and higher seasonality (Schmitz et al., 2001;Chen et al., 2018;Foreman et al., 2017;Pujalte et al., 490 2015). Although the MECO is not an abrupt event like other hyperthermals, instead a more extended period of gradual warming (ca. 500 kyr; Bohaty et al., 2009), we also observe this progradation 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 originated upstream, for 495 instance, linked to enhanced sediment remobilization (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 point towards an increase in soil 505 remobilization, erosion, and transport in continental environments during the MECO event.

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 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 ẟ 18 Ocarb values (Mulch et al., 2015). In contrast, in the Tarim 510 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., 2019Benyamovskiy et al., 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 515 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 ẟ 13 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). 520 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 525 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., 2012Foreman et al., , 2014Pujalte et al., 2015;Foreman and Straub, 2017;Honegger et al., 2020). Hence, despite its important differences with the early Eocene hyperthermals, the MECO shares several attributes with them around the warming peak. In summary, our results point to a 530 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.

Conclusions
In the South-Pyrenean Foreland Basin, an important progradation affected the entire sediment routing system from fluvial to deltaic environments at times 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 isotopes 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, for the first time reported in the Pyrenean region. In Yebra de Basa, an organic-rich interval of terrestrial origin is found before the main 540 deltaic progradation. It is associated with a negative excursion in oxygen and carbon isotopes. The correlation between the MECO and the basin-wide progradation, and the new geochemical results present 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. 545 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, quantify sediment and organic export, and its relationship with carbon burial and silicate weathering. 550 Our results support the view that high-accommodation settings in foreland basins are important recorders of paleoenvironmental 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 highresolution but more condensed deep-sea paleoclimatic records. In particular, during high-CO2 globally warm episodes of the Earth's history when the carbonate-rich oceanic records may undergo intervals of non-deposition or dissolution. 555

Data availability
All the data (stable isotopes, clay minerals, organic matter, major and trace elements) can be found in the supplementary material.

Authors contribution
SPC led fieldwork, sampling, sample preparation, data interpretation, and writing. LV contributed to the fieldwork, data 560 interpretation, discussion, and writing. JES performed stable isotope analyses, data interpretation, discussion, and writing. TA performed XRD analyses, data interpretation, and discussion. JV and AV contributed to the field work preparation, sampling, and discussion. MT, SW and NS contributed to the discussion, and writing. CP contributed to fieldwork, discussion, and writing. AV and MG helped with magnetostratigraphic data interpretation and discussion. SC supervised the project, funding, interpretation, and writing. 565 of the Northeastern Peri-Tethys, Austrian Journal of Earth Science, V. 105/1. P. 117-128, 2012. Bijl, P. K., Houben, A. J., Schouten, S., Bohaty, S. M., Sluijs, A., Reichart, G. J., ... and Brinkhuis, H