Middle Eocene Climatic Optimum (MECO) and its imprint in the continental Escanilla Formation, Spain

The Middle Eocene Climatic Optimum (MECO) is a global warming event mainly recognized in the marine domain and described less in the terrestrial environment. Here we present a comprehensive geochemical record of the MECO from the Escanilla Formation, a fluvial sedimentary succession in the southern Pyrenees, Spain, based on a suite of sampled paleosols, fluvial stromatolites, and pedogenic nodules. Our multiproxy approach involves using carbon and oxygen stable isotope compositions to identify the regional preservation of the MECO, calculate chemical weathering intensity and mean annual precipitation, perform clumped isotopes on carbonates, and identify clay mineralogy assemblages of paleosols. Results indicate that the Middle Eocene interval in the southern Pyrenees was characterized by low weathering rates under warm and arid climatic conditions. This is further supported by the presence of smectite, palygorskite, illite, and chlorite, which suggest seasonal rainfall but under generally dry conditions resulting in weak chemical weathering. Importantly, a negative organic carbon isotopic excursion indicates the regional, terrestrial impact of the MECO, highlighting that fluvial sedimentary successions even in active foreland basins can represent particularly interesting terrestrial archives of past changes in global climate.

1 Introduction

The Middle Eocene Climatic Optimum (MECO) is a transient global warming event that occurred during the Bartonian (40.0 Ma) and represents a significant reversal in the long-term cooling trend of the Middle to Upper Eocene (Bohaty et al., 2009). It is characterized by a distinct negative oxygen isotope excursion (OIE) of about 1 ‰ interpreted as a 3–6 °C rise in global ocean temperatures, and a subdued and spatially non-uniform negative carbon isotope excursion (CIE) in marine carbonates (e.g., Bohaty et al., 2009; Bijl et al., 2010; Henehan et al., 2020). This lack of a distinctly negative and robust CIE, a key diagnostic feature of earlier Cenozoic hyperthermals such as the Paleocene Eocene Thermal Maximum (PETM), and its comparatively long duration (∼ 500 kyr), along with elevated pCO₂ levels sustained over a long period by a diminished silicate weathering feedback (van der Ploeg et al., 2018), points towards an instability in the long-term carbon cycle during the MECO (Sluijs et al., 2013). In principle, elevated global temperatures and CO₂ levels are expected to be reflected in continental environments in the form of increased silicate weathering of exposed rocks. This process helps restore the climate through a drawdown of atmospheric CO₂ typically associated with an intensified hydrological cycle (e.g., Sluijs et al., 2013; Methner et al., 2016; van der Ploeg et al., 2018).

While the MECO has been extensively identified in both onshore and offshore marine sedimentary sections (e.g., Jovane et al., 2007; Bohaty et al., 2009; Spofforth et al., 2010; Gandolfi et al., 2023), only a few studies have identified it in terrestrial successions (e.g., Bosboom et al., 2014; Mulch et al., 2015; Methner et al., 2016). Currently available data from the terrestrial realm suggest important regional differences in response to the MECO. For instance, Bosboom et al. (2014) discuss the onset of region-wide arid conditions in the Xining Basin, China, while Methner et al. (2016) predicted large temperature seasonality and an intensification of the hydrologic cycle in the continental interiors of the northwestern USA. These regional differences highlight the need to further document records of the MECO in continental depositional systems to provide more constraints for understanding the dynamics of the climate and earth's surface perturbations during that period.

Given the need for further research, we document a new record in the Escanilla Formation (Fm), a fluvial sedimentary succession in the southern Pyrenees at Olsón (Spain), using carbon stable isotope compositions (δ¹³C) of paleosol organic matter to highlight the regional preservation of the MECO. We also measured the carbon and oxygen stable isotope ratios of carbonates in stromatolites, pedogenic nodules, and bulk paleosols to further examine the terrestrial MECO-related changes in the carbon and water cycles. We calculated weathering indices to quantify weathering intensity during the MECO and reconstructed climatic conditions using mean annual precipitation estimates, preliminary carbonate clumped isotope thermometry data, and clay mineralogy.

2 Study locality and age constraints

2.1 The Escanilla sediment routing system

The Escanilla sediment routing system, situated in the Pyrenees, represents a well-preserved mid- to late Eocene-age (ca. 41–34 Ma) fluviatile succession. Catchment areas of the high Pyrenees were linked to the southern Pyrenean foreland basin through the Sis and Gurb paleovalleys (Fig. 1a) (Bentham et al., 1993; Labourdette and Jones, 2007; Labourdette, 2011; Michael et al., 2013). Extensive paleocurrent data suggest that these paleovalley systems predominantly derived sediments from the axial zone of the Pyrenees, converging in the Viacamp area (Fig. 1a). From this point, sediments were carried downstream to the west through the Ainsa Basin (Vincent, 2001; Whittaker et al., 2011; Parsons et al., 2012; Michael et al., 2013) into the shallow marine Jaca Basin (Fig. 1a) (Peris Cabré et al., 2023).

Figure 1 Panel (a) is a map of the Escanilla paleo-sediment routing system in the southern Pyrenees, Spain, also showing the main tectonic structures. Red arrows mark the water discharge and sediment transport direction of the Escanilla system away from the source regions of Sis and Gurb paleovalleys in the axial zone of the Pyrenees. (Figure panel modified after Michael et al., 2014.) Also displayed is an inset map of Spain, indicating the study area near the village of Olsón. (Map modified from Labourdette and Jones, 2007.) Panel (b) is a field image depicting the sampled Escanilla Fm. Panel (c) shows that the lithostratigraphic framework of the Escanilla Fm at Olsón consists of two main members – the Mondot and the Olsón members – with the Olsón conglomerate (OC; red line) at the transition between the two members. Panel (d) is a geological map of the Escanilla Fm around Olsón. (This map was prepared using QGIS Desktop 3.22.8 (https://qgis.org/en/site/, last access: 17 June 2022) Panel (e) shows a sampled composite section, with the position of each collected sample, of the Escanilla Fm together with the local magnetostratigraphic interpretation by Vinyoles et al. (2020) correlated to the Geomagnetic Polarity Time Scale (GPTS 2020) (Ogg, 2020). The thickest normal magnetozone C18n in the local magnetostratigraphic interpretation includes C18n.1n, C18n.1r, and C18n.2n.

The entire Escanilla sediment routing system has been meticulously documented within a comprehensive source-to-sink framework based on provenance tools, including clast lithologies (Mesozoic carbonates, Upper Carboniferous to Triassic clastic and igneous rocks, Hercynian granites, and Paleozoic basement), heavy minerals, U–Pb geochronology of detrital zircons, apatite fission track analysis, paleocurrent analysis, as well as magneto- and biostratigraphy. A detailed explanation of this system can be found in the works of Michael et al. (2013, 2014) and references therein.

2.2 Escanilla Fm

The Escanilla Fm at Olsón has a maximum thickness of 1000 m and is subdivided into the Mondot and Olsón members (Fig. 1c) (e.g., Labourdette and Jones, 2007). At the transition of these two members lies a basin-wide extending conglomeratic channel–complex, hereafter referred to as the Olsón Conglomerate (OC), which is the interpreted stratigraphic expression of peak MECO warming in the Escanilla Fm based on available age constraints (Fig. 1b–e). The Escanilla Fm predominantly consists of conglomeratic to sandy channel-fill deposits and their adjacent fine-grained floodplain/overbank deposits. These floodplain fines have been previously characterized as Entisols, displaying a low degree of pedogenesis (Dreyer et al., 1992). In the field, paleosols were identified based on pedogenic features such as horizonation, coloring, root traces, and burrowing. In addition, we document the presence of three pedogenic carbonate nodule horizons (Fig. 2c) and nine fluvial stromatolites preserved within the Escanilla floodplains (Fig. 2d–f). Although fluvial stromatolites are common in the upper Paleocene to Eocene fluvial deposits of Spain (Zamarreño et al., 1997), to our knowledge, they have not been identified in the Escanilla Fm before. These stromatolites have been preserved as elongated domes with an asymmetrical shape and have varying lengths from a few centimeters to a meter, and range in diameter from 15 to 60 cm. They have been interpreted as having formed as overbank deposits of fluvial channels with a preferential elongation presumably parallel to the flow direction (Zamarreño et al., 1997). These stromatolites are also similar in appearance to Eocene stromatolites that grew on tree trunks from the Green River Formation, Colorado, USA (Awramik and Buchheim, 2015), although further research would be required to investigate this.

Figure 2 Panel (a) is an outcrop panorama depicting the Olsón Conglomerate (OC) and the thick floodplain interval sampled for paleosols. Panel (b) shows a typical sampled paleosol. Panel (c) shows pedogenic carbonate nodules sampled below the weathering depth. Panel (d) shows long, tubular fluvial stromatolites (white arrows) preserved in the floodplain. Panels (e, f) show a cross-section of the stromatolite samples NS7 and NS1, respectively, with well-preserved individual layering. The white stars indicate micro-drilling sites.

Within this context, a composite section of about 120 m was sampled for a suite of paleosol (ca. 200 g of fine-grained and fresh (unweathered) rock) and carbonate (pedogenic nodules and stromatolites) samples, collected below the weathering depth (at least 20 cm below modern surface) to avoid diagenetic and grain size bias (e.g., Lupker et al., 2011). Locations of sampled sites and corresponding raw data from different analyses have been provided as Supplement.

2.3 Age constraints

The age model used in this study is based on the magnetostratigraphic framework, involving high- and intermediate-quality samples, of the Escanilla Fm by Vinyoles et al. (2020) and agrees well with the magnetostratigraphic interpretation of Bentham et al. (1993). Age constraints by Vinyoles et al. (2020) were also used by Peris Cabré et al. (2023) for the identification of the MECO at Belsue and Yebra de Basa in the Jaca Basin.

To compare our geochemical data relative to the target MECO isotopic excursion from ODP site 738 (Bohaty et al., 2009) and the Geomagnetic Polarity Time Scale (GPTS 2020) (Ogg, 2020), we rely on a preliminary first-order linear scaling by matching the base (at 40 m) and top (at 210 m) of Chron C18n in the magnetostratigraphic data of Vinyoles et al. (2020) to the base of Chron C18n.2n and the top of Chron C18n.1n on the GPTS 2020 (Fig. 1e). All data are presented relative to the thickness of the sampled Olsón section accompanied by the local magnetostratigraphic interpretation of Vinyoles et al. (2020), linearly scaled to the GPTS 2020.

3 Materials and methods

Analyses were carried out on a suite of sampled paleosols (n= 45), stromatolites (n= 9), and pedogenic nodules (n= 9 from three sample horizons) (Fig. 1d). Pedogenic nodules were small and could not be cut open to be examined by cathodoluminescent microscopy or transmitted light. Hence, entire nodules were crushed and their homogenized powders were used for analysis.

Powders of bulk paleosol samples were prepared and analyzed for geochemical indicators, including total organic carbon (TOC) content, Rock-Eval parameters, organic carbon isotope compositions (δ¹³C_org), carbonate carbon and oxygen isotopes (δ¹³C_carb and δ¹⁸O_carb), major elements, and clay mineral assemblages. Stromatolite and pedogenic nodule powders were analyzed for δ¹³C_carb and δ¹⁸O_carb, as well as clumped isotope compositions (δ⁴⁷ and Δ₄₇).

3.1 Rock-Eval analysis

Rock-Eval pyrolysis was performed at the Institute of Earth Sciences of the University of Lausanne (ISTE-UNIL) to assess the total organic carbon (TOC) content in the paleosols. The analyses were performed using a Rock-Eval 6 instrument, following the procedure described by Behar et al. (2001).

3.2 Organic carbon isotopes

The carbon isotope compositions of organic matter in paleosol samples were analyzed in the stable isotope laboratories of the Institute of Earth Surface Dynamics, University of Lausanne (IDYST-UNIL). Samples first underwent de-carbonatation with 10 % $v/v$ HCl and then were thoroughly washed with deionized water and dried at 40 °C for 48 h. The δ¹³C_org measurements were made using a Carlo Erba 1100 (Fisons Instruments, Milan, Italy) elemental analyzer connected to a Thermo Fisher Scientific Delta V Plus isotope ratio mass spectrometer, both operated under continuous helium flow. Measured δ¹³C values were calibrated and normalized using international reference materials and in-house standards (Spangenberg, 2006, 2016) and reported in permil (‰) vs. the Vienna Pee Dee Belemnite limestone standard (VPDB). The precision of the δ¹³C_org values was better than 0.1 ‰.

3.3 Carbonate isotopes

Carbon and oxygen isotope compositions were determined at the IDYST-UNIL laboratories. Bulk paleosol samples containing > 10 wt% CaCO₃, including carbonate samples, were analyzed using a Thermo Fisher Scientific Gas Bench II carbonate preparation device connected to a Delta V Plus isotope ratio mass spectrometer according to a method adapted after Spötl and Vennemann (2003). CO₂ gas was produced by reaction with 99 % orthophosphoric acid at 70 °C. The δ¹³C_carb and δ¹⁸O_carb values are reported in permil vs. VPDB. Replicate measurements of the international calcite standard NBS 19 (limestone; δ¹³C $=+$ 1.95 ‰, δ¹⁸O = −2.19 ‰) and an in-house standard (Carrara marble; δ¹³C $=+$ 2.05 ‰, δ¹⁸O = −1.7 ‰) yielded an analytical precision (1σ) of ±0.05 ‰ for δ¹³C_carb and ±0.1 ‰ for δ¹⁸O_carb.

3.4 Major elements

SiO₂, Al₂O₃, Fe₂O₃, MnO, MgO, CaO, Na₂O, K₂O, P₂O₅, Cr₂O₃, NiO, and loss on ignition (LOI) were measured in the powdered bulk paleosol samples by X-ray fluorescence (XRF; Phillips PANalytical PW2400 spectrometer) at the ISTE-UNIL laboratories. The analysis was performed on fused glass discs prepared with 1.2000±0.0005 g ignited sample powder and 6.0000±0.0005 g of lithium tetraborate (Li₂B₄O₇). The concentrations of the major elements were expressed as weight percent of oxides. The analytical precision (1σ) assessed by replicate analysis of international reference materials is 0.4 %.

3.5 Weathering indices

The chemical index of alteration (CIA in %; Eq. 1), proposed by Nesbitt and Young (1982), was used to quantify the degree of weathering by using the molar ratio of immobile Al₂O₃ and the mobile oxides CaO, Na₂O, and K₂O in the silicate fraction (e.g., Deng et al., 2022). The CIA is commonly used to estimate the intensity of alteration and test for environmental factors (e.g., temperature, precipitation, elevation, and slope) affecting silicate weathering:

$$\begin{array}{}\text{(1)}& \mathrm{CIA}={\displaystyle \frac{{\mathrm{Al}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{3}}}{{\mathrm{Al}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{3}}+{\mathrm{Na}}_{\mathrm{2}}\mathrm{O}+{\mathrm{CaO}}^{*}+{\mathrm{K}}_{\mathrm{2}}\mathrm{O}}}\times \mathrm{100},\end{array}$$

where CaO^* is the CaO incorporated in the silicate fraction and is calculated as:

$$\begin{array}{}\text{(2)}& {\mathrm{CaO}}^{*}=\mathrm{CaO}-\left[\left(\mathrm{10}/\mathrm{3}\right)\times {\mathrm{P}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{5}}\right].\end{array}$$

Intense weathering removes mobile oxides, concentrating Al₂O₃ to a maximum value of almost 100 wt %, whereas weak weathering would give lower CIA values (< 50 %) because of the dominance of mobile oxides.

For comparison with the CIA, we used the chemical index of weathering (CIW in %; Eq. 3) proposed by Harnois (1988), which is a modified version of the CIA, to account for potassium metasomatism in paleosols during weathering (see Fig. S7 in the Supplement) (Stein et al., 2021):

$$\begin{array}{}\text{(3)}& \mathrm{CIW}={\displaystyle \frac{{\mathrm{Al}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{3}}}{{\mathrm{Al}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{3}}+{\mathrm{Na}}_{\mathrm{2}}\mathrm{O}+{\mathrm{CaO}}^{*}}}\times \mathrm{100},\end{array}$$

where CaO^* is the CaO incorporated in the silicate fraction (Eq. 2).

3.6 Mean annual precipitation

Mean annual precipitation (MAP) was estimated from the CIW (Eq. 4; standard error $=\pm \mathrm{182}$ mm yr⁻¹) using the equation proposed by Sheldon et al. (2002) as

$$\begin{array}{}\text{(4)}& \mathrm{MAP}=\mathrm{221}{e}^{\mathrm{0.0197}\left(\mathrm{CIW}\right)}.\end{array}$$

3.7 Clumped (Δ₄₇) isotopes

Clumped isotope compositions were measured at the IDYST-UNIL laboratories using a Nu Perspective dual-inlet mass spectrometer with an automated NuCarb sample preparation device. Carbonate digestion, the CO₂ purification procedure, and measurements followed Anderson et al. (2021). For Δ₄₇ analyses, 3.8–4.2 mg of carbonate were reacted at 70 °C with 110 µL 105 wt % phosphoric acid (H₃PO₄), and the liberated CO₂ was purified after digestion in a series of temperature-controlled, liquid-nitrogen-cooled cold fingers, including an adsorption trap (packed with Porapak Q $\mathrm{50}/\mathrm{80}$ mesh) held at −30 °C. We applied the acid fractionation factor of +0.088 ‰ (Anderson et al., 2021) to the data in order to compare with the ETH standards reacted at 90 °C (Bernasconi et al., 2021).

Three standards (ETH-1, ETH-2, and ETH-3) were used to calculate Δ₄₇ values of unknown samples using linear regression with the respective Δ₄₇ values of 0.205, 0.209, and 0.613 (Bernasconi et al., 2021). Temperatures were calculated using the calibration of Anderson et al. (2021), where Δ₄₇ is in weight percent and T in degrees Celsius, as

$$\begin{array}{}\text{(5)}& {\mathrm{\Delta }}_{\mathrm{47}}=(\mathrm{0.0391}\pm \mathrm{0.0004})\times \left({\mathrm{10}}^{\mathrm{6}}/T^{\mathrm{2}}\right)+(\mathrm{0.154}\pm \mathrm{0.004}).\end{array}$$

Multiple analyses were carried out on the drilled powders of each stromatolite (n= 2–4) and pedogenic nodule (n= 4–8) samples to obtain better statistics (see data in the Supplement).

3.8 Clay mineralogy

Clay mineralogical assemblages in paleosol samples were determined by X-ray diffractometry (XRD) at the ISTE-UNIL laboratories. Samples were prepared following the procedure described in Adatte et al. (1996). Analyses were made using a Thermo Fisher Scientific ARL X-TRA diffractometer and the intensities of the XRD peaks characteristic of each mineral were used to estimate the relative percentage in bulk rock and < 2 and 2–16 µm clay size fractions.

3.9 Uncertainty in reported data

All data reported in this study are associated with uncertainties in the form of standard error of the mean (SE) calculated as $\mathrm{SE}=\frac{\mathrm{SD}}{\sqrt{n}}$, where SD is the standard deviation and n is the number of replicates analyzed. Uncertainty propagation was done using the uncertainties package of Python (Spyder 4.0.1), which is an open-source and cross-platform program that handles calculations with numbers involving uncertainties.

4 Results and discussion

4.1 TOC content and δ¹³C of bulk paleosol organic matter

TOC content in paleosol samples varies from 0.01±0.01 to 0.57±0.01 wt % with an average value of 0.07±0.01 wt % (N= 45) (Fig. 3). Low TOC may be indicative of low primary productivity, in this case “vegetation” including grasses and higher plants, and (cyno)bacteria or low preservation of organic matter due to an oxidizing (oxygenated) environment (Tyson, 1995).

Figure 3 Total organic carbon (TOC) and dispersed organic carbon isotope compositions (δ¹³C_org), as well as associated standard error, in paleosol samples (circles) presented with a 7-point moving average. Also displayed is the change in benthic foraminiferal (genus Cibicidoides) carbon isotope ratios (δ¹³C_cib) from ODP site 738.

The δ¹³C_org values of the paleosols have a range between $-\mathrm{26.0}\pm \mathrm{0.2}$ ‰ and $-\mathrm{20.4}\pm \mathrm{0.2}$ ‰ with an average value of $-\mathrm{23.2}\pm \mathrm{0.2}$ ‰ (N= 45) (Fig. 3). A negative CIE is marked by a 3 ‰ shift from the base of the section (from 0 to 30 m), where the onset begins, followed by a plateau of low values (from 30 to 50 m) that gradually return to higher values 60 m upwards. This negative CIE is most likely coeval to the 0.5 ‰ negative excursion observed in the benthic foraminifera δ¹³C_cib values from ODP sites 738 (Fig. 3) (Bohaty et al., 2009), indicating general agreement in the change of δ¹³C values, even though absolute differences in the magnitude of excursions exist.

A similar large-magnitude negative CIE was previously identified for the PETM within the intermontane Piceance Creek Basin of western Colorado (USA), where a negative CIE of about 3 ‰ was reported (Foreman et al., 2012). The δ¹³C_org record from the Middle Eocene Alano di Piave section deposited in the marginal Tethys Ocean recorded a negative CIE of about 1 ‰ (Spofforth et al., 2010), while the coeval shallow water Sealza section from Italy recorded a negative CIE of 2 ‰ (Gandolfi et al., 2023).

The δ¹³C_org values can also be used as indicators of paleoecology and paleoclimate (Kohn, 2010). C₃ plants which include trees, most shrubs, and cool-season grasses, have δ¹³C values between −37 ‰ and −20 ‰, and have dominated the history of terrestrial vegetation (Kohn, 2010). This wide range in δ¹³C values of plants is dependent on several factors such as temperature, altitude, latitude, and MAP (Schulze et al., 1996; Kohn, 2010). Non-water-stressed C₃ plants are enriched in ¹²C and hence have more negative δ¹³C values, typically lower than −26 ‰. Higher δ¹³C values (> −26 ‰) are associated with plants growing under water-deficient conditions and low soil transpiration rates (MAP < 500 mm yr⁻¹) (e.g., Cerling and Quade, 1993; Kohn, 2010; Methner et al., 2016). Measured δ¹³C_org values suggest a predominance of C₃ vegetation consistent with an Eocene ecosystem (Cerling and Quade 1993; Methner et al., 2016). A significant proportion of measured values have relatively high δ¹³C values (> −23 ‰) that are characteristic of dry environments with MAP < 500 mm yr⁻¹ (Kohn, 2010). Low primary productivity and low organic matter preservation complemented by elevated δ¹³C_org values likely indicates sparse vegetation in a dry and arid ecosystem.

Vegetation plays a key role in influencing landscape response to terrestrial hydroclimates during global warming events such as the PETM (Foreman et al., 2012). For instance, de-vegetated fluvial banks enhance sediment erodibility, bank erosion, and lateral channel mobility, resulting in peak discharge and bedload sediment downstream flux even under a negligible increase in precipitation (Gran and Paola, 2001; Barefoot et al., 2021). Sparse vegetation in the fluvial Escanilla Fm reflects the response of the Escanilla rivers to peak discharge and sediment flux events, which modified the stratigraphic architecture described as alternating sequences of high amalgamation (HA) and low amalgamation (LA) intervals (Sharma et al., 2023). Such large flux in the Escanilla Fm at Olsón during the MECO is also consistent with the clastic progradation event observed in the Jaca Basin (Peris Cabré et al., 2023), indicating the transmission and preservation of sedimentary flux signals in the sedimentary record from source-to-sink. Lastly, enhanced channel mobility during the MECO could also have led to the export of large quantities of fine-grained (floodplain) sediment to the marine domain and potentially affected the carbonate environments developing at that time on the southern Margin of the Jaca–Pamplona trough (e.g., Garcès et al., 2023).

4.2 δ¹³C and δ¹⁸O of paleosol bulk carbonates, stromatolites, and pedogenic nodules

The δ¹³C_carb values in paleosol bulk carbonates have a range of $-\mathrm{2.6}\pm \mathrm{0.01}$ ‰ to $-\mathrm{1.0}\pm \mathrm{0.01}$ ‰ with an average of $-\mathrm{1.6}\pm \mathrm{0.01}$ ‰ (N= 45) (Fig. 4). A transient decrease in δ¹³C_carb values is observed at 40 m where a negative CIE of about 1.0 ‰ magnitude is considered to represent the MECO negative CIE in the Escanilla Fm, followed by an increase in δ¹³C_carb values to −1.5 ‰ towards the top of the section (Fig. 4). These values are most likely synchronous to the negative CIE recorded in paleosol organic matter.

Figure 4 Carbon isotope compositions (δ¹³C_carb) and associated standard error from paleosol bulk carbonates (circles) with a 7-point moving average, and from stromatolite (triangles) and pedogenic nodules (squares) for the Olsón section. Also marked is the MECO onset and peak interval based on δ¹³C_org values from this study.

δ¹³C_carb values in stromatolites range from $-\mathrm{3.9}\pm \mathrm{0.1}$ ‰ to $-\mathrm{1.3}\pm \mathrm{0.1}$ ‰ with an average value of $-\mathrm{2.5}\pm \mathrm{0.1}$ ‰ (N= 63; 5–9 replicate measurements from nine sample horizons), while δ¹³C_carb values in pedogenic nodules range from $-\mathrm{4.9}\pm \mathrm{0.3}$ ‰ to $-\mathrm{1.1}\pm \mathrm{0.3}$ ‰ with an average value of $-\mathrm{3.2}\pm \mathrm{0.3}$ ‰ (N= 29; 9–10 replicate measurements from three sample horizons). Since sample size is limited stratigraphically, it does not permit a direct evaluation of the isotopic signal relative to the MECO. However, δ¹³C_carb values from stromatolites and pedogenic nodules show a consistent 1 ‰–2 ‰ negative offset when compared with values from paleosol bulk carbonates, which is possibly due to the presence of detrital Mesozoic carbonates (δ¹³C_carb= 0 ‰; Zamarreño et al., 1997) in the bulk sediments from the source area.

The δ¹⁸O_carb values in paleosol bulk carbonates have a range between $-\mathrm{6.7}\pm \mathrm{0.1}$ ‰ and $-\mathrm{4.2}\pm \mathrm{0.1}$ ‰ with an average of $-\mathrm{5.8}\pm \mathrm{0.1}$ ‰ (N= 45) (Fig. 5). A positive OIE of ca. 0.5 ‰ magnitude at 40–60 m suggests an increase in freshwater ¹⁸O content and perhaps represents peak MECO conditions in the Escanilla Fm. Peak warming would correspond to the OC where the highest discharge and flux estimates have been predicted by Sharma et al. (2023). Following the positive OIE, δ¹⁸O_carb values decline and return to relatively stable values of around −6.0 ‰, 60 m onwards until the top of the Olsón section and may represent the post-MECO cooling phase (Fig. 5).

Figure 5 Oxygen isotope compositions (δ¹⁸O_carb) from paleosol bulk carbonates (circles) with a 7-point moving average and associated standard error, as well as carbonate samples (stromatolites (triangles) and pedogenic nodules (squares)) for the Olsón section. Also marked are the MECO onset and peak interval based on δ¹³C_org values from this study.

δ¹⁸O_carb values in stromatolites range from $-\mathrm{8.1}\pm \mathrm{0.1}$ ‰ to $-\mathrm{4.7}\pm \mathrm{0.1}$ ‰ with an average value of $-\mathrm{7.1}\pm \mathrm{0.1}$ ‰ (N= 63), while δ¹⁸O_carb values in pedogenic nodules range from $-\mathrm{7.6}\pm \mathrm{0.1}$ ‰ to $-\mathrm{6.5}\pm \mathrm{0.1}$ ‰ with an average value of $-\mathrm{6.9}\pm \mathrm{0.1}$ ‰ (N= 29); and crudely match the δ¹⁸O trend in paleosol bulk carbonates (Fig. 5).

In summary, irrespective of the presence of authigenic and detrital carbonates in paleosol samples, the negative CIE in paleosol organic matter suggests that the MECO can be regionally recognized in the Escanilla Fm. Stable isotope data from the Escanilla Fm at Olsón are also compatible with climate perturbations through excursions similar to the isotope excursions in the marine records, even though there are differences in the magnitude of excursions. These excursions have also been identified downstream in the time-equivalent marine sediments in the Jaca Basin, Spain (Peris Cabré et al., 2023), indicating the preservation of MECO climate perturbation signals in the source-to-sink Escanilla sediment routing system.

4.3 Primary vs. diagenetic signals

Carbon and oxygen isotope composition of bulk paleosol carbonates may be affected by the diagenetic alteration of mineral phases. It is therefore important to evaluate the potential diagenetic overprint on primary geochemical signatures (e.g., Marshall, 1992).

The degree of alteration was assessed through the relationship between δ¹³C_org and δ¹³C_carb, as well as between δ¹⁸O_carb and δ¹³C_carb values. The Pearson correlation coefficient, r < 0.6, indicates a statistically non-significant relationship and indicates that a diagenetic overprint on the primary signal can be excluded (e.g., Fio et al., 2010). In both correlation plots (Fig. 6), no statistically significant correlation was found (δ¹³C_carb vs. δ¹³C_org: r= 0.21 (P= 0.16, N= 45); δ¹³C_carb vs. δ¹⁸O_carb: r= 0.05 (P= 0.74, N= 45)) indicating almost none or very minor diagenetic modification of the primary signal. Also, no correlation trend was observed between TOC and δ¹³C_org (Fig. S1).

Figure 6 Panel (a) is a scatter plot of paleosol δ¹³C_carb vs. δ¹³C_org values and (b) is a scatter plot of paleosol δ¹³C_carb vs. δ¹⁸O_carb values. For both plots, the Pearson correlation coefficient (r) and regression line are shown.

Maximum Temperature (T_max) from Rock-Eval analysis was used as a second approach to assess diagenetic alteration. T_max obtained in samples with high TOC (> 0.5 wt %) was < 440 °C which is the beginning of the oil window and indicates immature organic content (ca. 60 °C; Espitalié et al., 1985).

As a third approach, paleosol samples (S9, S12, S17, and S30) were analyzed using scanning electron microscopy (SEM). SEM images show the presence of authigenic calcite (Figs. S2–S4), the presence of authigenic clay minerals palygorskite and smectite (Fig. S5), as well as detrital illite and chlorite (Fig. S6). Collectively, the three approaches suggest that the primary signal is largely preserved in the Escanilla Fm at Olsón.

4.4 Weathering conditions

CIA values ranged from 9±1 % to 43±1 % with an average value of 19±1 % (N= 45). A slight peak in CIA values may represent peak weathering conditions during the MECO (Fig. 7). CIA values have further been compared with CIW values (Fig. S7), which take into account potassium metasomatism and range from 10±1 % to 48±1 % with an average value of 20.5±1 % (N= 45). Such low values are broadly indicative of weak chemical weathering, and if true, suggest a prevalence of physical weathering and erosion during the Middle Eocene. This contrasts with the assumption of intense weathering during global warming and recently published results for earlier hyperthermals (e.g., Tanaka et al., 2022). For instance, during the PETM, CIA values have been estimated to be in the 75 %–85 % range (Stokke et al., 2021). Our relatively low CIA values at Olsón could be related to the long-term trend of low silicate weathering in response to elevated pCO₂ and warming levels during the Middle Eocene as indicated by osmium isotopes (van der Ploeg et al., 2018) and more recently by lithium isotope data available from the marine environment (Krause et al., 2023). Current understanding of the earth's carbon cycle suggests strengthening of the negative silicate weathering feedback in response to rising atmospheric pCO₂ (Colbourn et al., 2015; Penman et al., 2020). However, the strength of the feedback depends primarily on several local/regional environmental variables such as temperature and precipitation, in addition to geomorphology and lithology (e.g., Richey et al., 2020; Deng et al., 2022). Our estimated weathering intensities are in line with a predicted dry and arid ecosystem (based on δ¹³C_org values) and therefore further reinforce our interpretation of a local dry and arid climate with a weakened hydrologic system in northern continental Spain during the Middle Eocene. Such low weathering rates, if confirmed at a scale relevant to the global carbon cycle, could also explain the sustained elevated carbon levels for a longer duration than during the previous hyperthermals, highlighting the different dynamics of the MECO (Sluijs et al., 2013; van der Ploeg et al., 2018).

Figure 7 The chemical index of alteration (CIA) and associated standard error values from paleosols (circles) displayed with a 7-point moving average to quantify chemical weathering in the Escanilla Fm at Olsón.

4.5 Regional climate record during the Middle Eocene

4.5.1 Mean annual precipitation estimates

MAP values in the Olsón section range from 270±10 to 570±10 mm yr⁻¹ with an average of 330±10 mm yr⁻¹. Values stay constant at 300±10 mm yr⁻¹ until 40 m followed by a 20 % increase in precipitation, up to 370±10 mm yr⁻¹, which most likely corresponds to the OC (peak MECO conditions). Above the OC, MAP values return to an average value of 340±10 mm yr⁻¹ until the top of the section (Fig. 8). Overall, these values predict arid to semi-arid climate in this area of the southern Pyrenees during the Middle Eocene and are coherent with the high δ¹³C_org values (water-stressed environments) as well as the low CIA and CIW (diminished chemical weathering) values in the studied section. At Igualada in the Ebro Basin, 200 km away from Olsón, palynological, pollen taxa, and floral diversity studies suggest warm climate and humid vegetation, with preservation of mangrove swamp vegetation along the coast (Cavagnetto and Anadón, 1996; Haseldonckx, 1972). The absence of humid climate in Olsón could be due to its location being higher in elevation and away from the coastline as compared with Igualada. Such regional differences in climate could also be the result of a climate transition phase during the Middle Eocene, oscillating from a warm tropical Early Eocene to a cold and arid Early Oligocene, expressed differently in different regions and at possibly different sampled intervals.

Figure 8 Mean annual precipitation (MAP) values (white circles) and the associated standard error from the Olsón section are presented using a 7-point moving average, followed by mean clumped isotope (Δ₄₇) temperatures and the standard error from replicate measurements of stromatolites (triangles) and pedogenic nodules (squares). We compare our terrestrial temperature estimates to sea surface temperature (SST) proxies TEX₈₆ (gray circles) and $U_{\mathrm{37}}^K$ (black circles) (Bijl et al., 2010).

MAP estimates based on well-dated megafloras from the Weisselster and Lausitz basins (both in northeast Germany), consisting of shallow marine and continental deposits, are in the range of 1100–1400 mm yr⁻¹ (Mosbrugger et al., 2005). Other proxy data from southern France indicate a MAP less than 500 mm yr⁻¹ in the Bartonian (Kocsis et al., 2014) and is similar to our calculation from Spain. In conclusion, the values reported here should be regarded as being representative of a local signal, most likely influenced by rain-shadow effects imposed by the Pyrenean topography at that time, which rose to 2000 m between 49 and 41 Ma (Huyghe et al., 2012), thus possibly inducing orographic effects as observed in the modern situation (Vacherat et al., 2017; Huyghe et al., 2018).

4.5.2 Carbonate clumped isotope thermometry

Δ₄₇ values of carbonate in stromatolites have a range from 0.634 ‰ to 0.704 ‰, translating into temperatures of 18–43 °C (Fig. 8). Each stromatolite was analyzed two to four times for replicate measurements with the standard error ranging from 0.006 ‰ to 0.011 ‰. Δ₄₇ values of carbonate in pedogenic nodules have a range from 0.629 ‰ to 0.704 ‰, translating into temperatures of 18–45 °C (Fig. 8). Each pedogenic nodule was analyzed four to eight times for replicate measurements having a standard error ranging from 0.006 ‰ to 0.009 ‰.

While measurement reproducibility was good, care must be taken while interpreting Δ₄₇ data as a number of potentially significant uncertainties are associated with them. For instance, Δ₄₇ temperatures may not necessarily reflect primary formation temperatures but could instead be the result of a combination of primary formation temperatures and secondary effects, such as potential diagenetic temperatures, that bias primary compositions, although secondary overprinting is unlikely to produce cooler temperatures (Hren and Sheldon, 2020). Second, organic contaminants could cause Δ₄₇ values to be variable, although any volatile part would be discarded when the samples were dried at high temperatures overnight (70–80 °C). Third, more replicate measurements from the same homogenized powder are required to assess interference by contaminants in the drilled powders and to better constrain the spread in data. For preliminary data, such as those presented here, the mean value could be considered as a good temperature estimate. Finally, significant diagenetic alteration could cause Δ₄₇ values to be variable, although a similar temperature range in both stromatolites and pedogenic nodules further suggests that analyzed samples most likely did not undergo significant diagenetic alteration after their formation. This would, however, need to be verified using petrography and/or cathodoluminescence.

Figure 9 Identified clay mineral assemblages (smectite, palygorskite, illite, and chlorite) and smectite/illite ratios are presented with respective standard error values and a 7-point moving average across the studied Escanilla Fm at Olsón. These minerals constitute up to 98 % of the total clay mineralogy and suggest the presence of an arid to semi-arid climate under weak chemical weathering conditions during the MECO in the southern Pyrenees, Spain.

Mean temperatures vary from 32.1±1.8 °C to 38.6±0.7 °C in the lower half of the section until 50 m followed by a peak mean temperature of 42 °C, just 25 m above the OC, without any observed change in lithology. Above 90 m until the top of the section, values return to an average value of 30 °C. Based on the available age constraints, our results suggest a potential lag between marine and terrestrial MECO climate signals. Our results also suggest a land–sea temperature gradient of 5–10 °C when compared with sea surface temperature records from ODP site 1172 (Tasman shelf (∼ 65° S paleolatitude, southwestern Pacific Ocean; Bijl et al., 2010), as well as IODP sites U1408 and U1410 (northwest Atlantic Ocean; van der Ploeg et al., 2023), most likely indicating an amplifying effect due to continentality. Similar continental temperature sensitivity during the Middle Eocene has also been previously identified through clumped temperatures of pedogenic carbonates in the continental interiors of southwest Montana, USA (Methner et al., 2016). Further research and sample analysis would, however, be required to investigate this in more detail.

Δ₄₇ temperatures were further used to calculate the δ¹⁸O of fluids in equilibrium with carbonates using the temperature-dependent fractionation factor of Epstein et al. (1953). Δ₄₇ temperatures and δ¹⁸O_carb values in stromatolites give water δ¹⁸O values in the range of −5.9 ‰ to −2.1 ‰ (average of −3.3 ‰), while values from pedogenic nodules give water δ¹⁸O values ranging from −6.5 ‰ to −0.8 ‰ (average of −3.6 ‰) (see data in the Supplement). For comparison, we also used the approach of Kim and O'Neil (1997), which gave water δ¹⁸O values from stromatolites to be in the range of −6.2 ‰ to −2.3 ‰ (average of −3.5 ‰) and −6.7 ‰ to −0.9 ‰ (average of −3.8 ‰) from pedogenic nodules (see the Supplement). These calculated water isotope values are consistent with meteoric water isotope composition at low latitude but signify that ¹⁸O-isotope enrichment was most likely due to excess evapotranspiration under arid climatic conditions, consistent with estimates from other proxies used in this study.

4.5.3 Clay mineralogy

Clay mineral assemblages in paleosols are also important paleoclimatic indicators and reflect detrital mineral input, composition of the source area lithology, and type of weathering of the source rocks, to provide integrated records of the overall climate (Singer, 1984; Franke and Ehrmann, 2010; Rego et al., 2018).

Smectite, palygorskite, illite, and chlorite form up to 98 % of the identified mineral assemblages in the studied section (Fig. 9). Smectite, commonly derived from alteration in volcanic rocks, forms under seasonal rainfall conditions with a pronounced dry season (e.g., Singer, 1984; Tabor and Myers, 2014) and constitutes on average 17 % of the total identified clay mineral assemblage, while individual values reach up to 58 %. Palygorskite, an authigenic mineral (Fig. S5), is indicative of an arid to semi-arid environment where evapotranspiration exceeds precipitation (e.g., Birkeland, 1984; Singer, 2002; Meunier, 2005), and constitutes up to 25 % (average of 10 %) of the total clay mineral assemblage. Illite content is between 11 % and 53 % (average of 36 %), while chlorite content is between 7 % and 64 % (average of 35 %) in the analyzed paleosols. High amounts of illite and chlorite, both detrital minerals (Fig. S6), are typically found in sediments formed by physical erosion of low-grade metamorphic rocks and are thus indicative of weak, incipient chemical weathering (Tabor and Myers, 2014; Rego et al., 2018), and support our calculated CIA values. Except for palygorskite, clay mineral variations can be inferred to reflect weathering conditions in the source area of the Escanilla sediment routing system.

An increase in authigenic palygorskite content above the OC indicates that environmental conditions on land became drier, which is in agreement with previous studies showing an increase in aridification following the MECO interval in northwest China (Bosboom et al., 2014) and in the Neo-Tethys (Baskil section, eastern Turkey; Rego et al., 2018). This is also in agreement with our decreasing smectite/illite ratios above the OC (Fig. 9) and is consistent with values reported by Rego et al. (2018). Overall, interpretation of clay mineral assemblages corroborates well with environmental conditions deduced from organic carbon stable isotope data, weathering indices, and mean annual precipitation.

5 Conclusions

Detailed geochemical and mineralogical analysis of paleosols, stromatolites, and pedogenic nodules provides new insights into terrestrial records of the MECO in the Ainsa Basin of the southern Pyrenees, Spain. A negative CIE measured on organic matter indicates the local preservation and identification of the MECO in the fluvial Escanilla Fm, demonstrating that continental sedimentary successions can serve as important climate archives and highlighting stable isotope proxies as a powerful dating and correlation tool in notably difficult-to-date fluvial successions. Low CIA values in the Escanilla Fm suggest poor silicate weathering feedback in response to elevated pCO₂ levels and a prevalence of physical erosion during the Middle Eocene. This is compatible with an arid to semi-arid climate with a locally diminished hydrological cycle, as supported by low MAP estimates and identified clay mineral assemblages in the fluvial sedimentary succession. Carbonate clumped isotope thermometry suggests high temperatures up to 42 °C and a possible amplifying effect of 10–15 °C on continents, compared with temperature records from the deep sea.