Silicate weathering in the semi-arid Southern Pyrenees during the PETM: lithium isotope evidence

The Palaeocene-Eocene Thermal Maximum (PETM), a hyperthermal event ∼ 56 Ma ago, allows the Earth system response to abrupt climate change to be explored. Recent investigations link the PETM with a negative lithium isotope (δ⁷Li) excursion, interpreted as an increase in continental silicate weathering fluxes, which can regulate Earth's surface temperature over geological timescales. However, the silicate weathering response under different climatic regimes has yet to be constrained. Here we aim to address the chemical weathering response to the PETM in the semi-arid Southern Pyrenees, and to explore how different archives (i.e. clays and carbonate nodules) record the weathering changes.

We investigated two continental sections in the Southern Pyrenees. In the Esplugafreda section, we measured δ⁷Li values as a silicate weathering proxy and ε_Nd values as a provenance proxy in the clay minerals. In the Rin section, we characterised the PETM locally by analysing δ¹³C values in organic matter and examined the clay mineralogy in the paleosols, as well as measuring δ⁷Li values in clays and carbonate nodules to trace silicate weathering. In the Esplugafreda section, we observe temporally stable ε_Nd values, while the δ⁷Li_clays record shows two small positive excursions, one during the Pre-Onset Excursion (∼ 0.7 ‰) and a second during the body of the PETM (∼ 0.8 ‰). In the Rin section, the PETM is characterised by a negative carbon isotope excursion of 2.8 ‰. The clays consist mostly of illite/smectite, illite, kaolinite, and chlorite, consistent with a seasonal climate in the region, and we find a positive δ⁷Li_clays excursion of ∼ 0.8 ‰.

The combined δ⁷Li_clays and ε_Nd records indicate increased clay formation and increased silicate weathering fluxes in the semi-arid Pyrenees, while the sediment provenance was stable. The δ⁷Li values in the carbonate nodules indicate more variability, potentially due to clay contamination. Constrained by the bedrock type of dominantly reworked sediments and the seasonal precipitation regime, the initially low weathering fluxes, despite a comparatively high weathering intensity, evolved towards higher weathering fluxes with enhanced erosion during the PETM.

1 Introduction

Continental silicate weathering is a critical feedback mechanism that stabilises Earth's climate over geological timescales by regulating atmospheric CO₂ through the long-term carbon cycle (Maher and von Blanckenburg, 2023; Raymo and Ruddiman, 1992; Walker et al., 1981). Through the breakdown of silicate minerals, the transport of cations in river systems, and the precipitation and burial of carbonates in the ocean, silicate weathering sequesters atmospheric CO₂, acting as a natural climate thermostat. Understanding how this process responds to abrupt climate change is essential for evaluating its capacity to modulate carbon fluxes under a range of future warming scenarios.

The Palaeocene-Eocene Thermal Maximum (PETM), a hyperthermal event ∼ 56 Ma ago, resulted from the rapid release of greenhouse gases that triggered a 5–8 °C global temperature increase over a geologically brief interval (Dickens et al., 1995; Kennett and Stott, 1991; McInerney and Wing, 2011; Westerhold et al., 2009; Zachos et al., 2003, 2008). Global records from the PETM suggest increases in silicate weathering fluxes (e.g. Hessler et al., 2017; Pogge von Strandmann et al., 2021a; Jaimes-Gutierrez et al., 2025; Rush et al., 2025), while some local records have been interpreted to show increased weathering intensity (e.g. Ramos et al., 2022; Chen et al., 2023), underscoring the potential for weathering to buffer atmospheric CO₂ during extreme warming events. In the context of modern anthropogenic warming, these insights are crucial for understanding the capacity of natural systems to mitigate rising CO₂ levels (Zeebe et al., 2016; Carmichael et al., 2017 and references therein).

In mid-latitude records, a range of sedimentological, geochemical, and mineralogical proxies suggest that the PETM resulted in a hydrological perturbation with episodic extreme rainfall events, increased seasonality, and aridification, leading to a loss of vegetation, extreme flooding, and enhanced channel mobility (Barefoot et al., 2022; Carmichael et al., 2017; Chen et al., 2018; Rush et al., 2021; Schmitz et al., 2001; Schmitz and Pujalte, 2007; Vimpere et al., 2023). These changes were particularly pronounced in semi-arid regions such as the Southern Pyrenees (∼ 35° N paleolatitude, Fig. 1), where sedimentary records document hydrological seasonality, enhanced erosion, and increased sediment transport (Chen et al., 2018; Jaimes-Gutierrez et al., 2024; Prieur et al., 2024, 2025; Pujalte et al., 2015; Rush et al., 2021; Schmitz and Pujalte, 2007).

Figure 1 (A) Location map. (B) Palaeogeography of the Tremp-Graus Basin during the late Palaeocene, modified after Jaimes-Gutierrez et al. (2024) and references therein. (C) Sediment routing system during the late Palaeocene, with the floodplain Esplugafreda section and the more coastal Rin section, at the marine-continent transition, modified from Prieur et al. (2025). (D) Stratigraphy and correlation between the Esplugafreda and Rin terrestrial sections. CC, Claret Conglomerate. Alv. Ls., Alveolina Limestone.

The Southern Pyrenees (Fig. 1) offer an exceptional setting for investigating climate-driven weathering dynamics. This region experienced tectonic quiescence during the PETM (Rosenbaum et al., 2002), allowing for the isolation of the effects of climate and hydrology on weathering. Sedimentary records indicate enhanced hydrological seasonality and increased runoff, consistent with amplified denudation rates during this interval (Pujalte et al., 2015; Rush et al., 2021; Schmitz and Pujalte, 2007). In this study, we use lithium isotopes (δ⁷Li) as a proxy for silicate weathering and neodymium isotopes (ε_Nd) as a tracer for sediment provenance, in order to quantify the weathering responses in the Southern Pyrenees and to assess their regional contribution to CO₂ regulation during the PETM.

We focus on two continental floodplain sections to answer two primary questions: (i) What was the chemical weathering response to the PETM in the semi-arid Southern Pyrenees? (ii) How do different sedimentary archives, such as clays and carbonate nodules, record the weathering changes? In the Esplugafreda section, we measured δ⁷Li values in clay minerals as a weathering proxy, together with ε_Nd values in two clay size fractions to determine sediment provenance. In the Rin section, we characterised the PETM locally through δ¹³C measurements in organic matter, and analysed the clay mineralogy of paleosols, and δ⁷Li values in both clays and carbonate nodules. These geochemical and mineralogical datasets allow us to reconstruct weathering dynamics in the region and to assess how they compare with existing globally-distributed records of PETM weathering (Pogge von Strandmann et al., 2021a; Ramos et al., 2022; Chen et al., 2023; Jaimes-Gutierrez et al., 2025; Rush et al., 2025).

1.1 Silicate weathering as Earth's surface thermostat

Silicate weathering rates are influenced by climate (Dessert et al., 2003; West et al., 2005), vegetation (Moulton et al., 2000; Porder, 2019), lithology (Caves et al., 2016; Dessert et al., 2003; Murray and Jagoutz, 2024), and regolith properties (Caves Rugenstein et al., 2019; Kump and Arthur, 1997). Weathering is driven by the availability of fresh mineral surfaces, reactive fluids, and dissolution kinetics (Bufe et al., 2021; Maher and von Blanckenburg, 2023; Riebe et al., 2004). Denudation (D= erosion rate [E] + silicate weathering [W]) links surface processes to the carbon cycle because erosion supplies fresh minerals, enhancing CO₂ sequestration through chemical weathering (Anderson et al., 2007; Gaillardet et al., 1999; Hilton, 2023; Riebe et al., 2004; West et al., 2005).

Two end-member regimes can be used to describe chemical weathering dynamics. In supply-limited regimes, mature soils dominated by secondary clays shield bedrock, limiting fresh mineral exposure and resulting in low weathering rates (Goddéris et al., 2008). In kinetically-limited regimes, which are typical of high-relief areas with thin soils, weathering rates are controlled by mineral dissolution kinetics (Kump et al., 2000; Riebe et al., 2004; West et al., 2005). Investigating how climate and erosion interact to shape these regimes under hyperthermal events such as the PETM is thus essential for understanding the weathering mechanisms and rates underpinning Earth's carbon cycle feedbacks in a warming climate.

1.2 Lithium isotopes as a chemical weathering tracer

Secondary clay minerals, which form as a by-product of primary silicate rock dissolution, preferentially incorporate ⁶Li over ⁷Li, resulting in isotopically light clays and isotopically heavy waters (e.g. Pogge von Strandmann et al., 2020). As weathering progresses, both dissolved lithium and the clays that precipitate from it become isotopically heavier, with the δ⁷Li value of the water and soil being linked by an approximately constant fractionation factor (Pogge von Strandmann et al., 2021b). Measuring δ⁷Li values in detrital and carbonate archives therefore allows past weathering regimes to be reconstructed. Because carbonate weathering has a minimal influence on riverine lithium budgets, δ⁷Li variations primarily reflect silicate weathering processes (Kısakűrek et al., 2005). Consequently, lithium isotopes have become widely applied as a proxy for tracking clay mineral formation, thereby tracing silicate weathering intensity changes, both in modern systems (e.g. Dellinger et al., 2015, 2017; Pogge von Strandmann et al., 2023) and during past geological events (e.g. Misra and Froelich, 2012; Pogge von Strandmann et al., 2013, 2021a; Ramos et al., 2022; Jones et al., 2023).

Weathering congruency, which represents the balance between primary mineral dissolution and secondary clay mineral formation, determines the δ⁷Li composition of river waters and sediments (Dellinger et al., 2015; Zhang et al., 2022 and references therein). In rapidly eroding regions with low $W/D$, congruent weathering results in minimal isotopic fractionation, because clay formation is relatively low (Fig. 2). In contrast, incongruent weathering in soil-mantled environments with moderate $W/D$, such as floodplains with high clay formation, yields both clays and waters with higher δ⁷Li values (Figs. 2 and 3). Finally, in supply-limited regimes with high $W/D$, such as rainforests, there is no remaining primary rock material to weather, so pre-formed clays are re-dissolved, which drives solutions to low δ⁷Li values, but with a very low weathering flux (e.g. Dellinger et al., 2015). In modern rivers, clays take up their Li from solution with an approximately constant fractionation factor (Pistiner and Henderson, 2003; Pogge von Strandmann et al., 2023; Ramos et al., 2024), so their composition also mimics this boomerang curve (Pogge von Strandmann et al., 2023; Wei et al., 2025; Winnick et al., 2022).

Figure 2 Large-river dissolved lithium isotope composition (δ⁷Li_dissolved) corrected for bedrock composition (Δ⁷Li_{dissolved-bedrock}) plotted versus weathering intensity ($W/D$). Modified from Dellinger et al. (2015) and references therein.

In detrital sediment archives, only part of this boomerang trend is typically observed because of mixing of the neoformed clays with primary silicate material, especially at low $W/D$ conditions (Dellinger et al., 2017). Therefore, continental and marine detrital records may need to be interpreted differently (e.g. Pogge von Strandmann et al., 2021a; Ramos et al., 2022, 2024; Jones et al., 2023; Jaimes-Gutierrez et al., 2025; Rush et al., 2025; Wei et al., 2025). Because finer sediment fractions tend to be preferentially transported further offshore due to hydrodynamic sorting during river to marine transport, clay-sized records may be more clearly expressed in some marine sedimentary records (e.g. Gibbs, 1977; Liu et al., 2023). Such biases resulting from mixing with primary silicate grains in bulk sediment samples can potentially be reduced by analysing the clay size fraction (<2 µm), although this fraction can still also contain some primary minerals.

Finally, lithium isotopes can also be fractionated by direct climatic fluctuations. For example, temperature (Li and West, 2014; Vigier et al., 2008) and hydrological controls (Zhang et al., 2022) have both been found to influence the δ⁷Li composition of river water, and consequently the composition of the sedimentary archives that form in equilibrium with them (Pogge von Strandmann et al., 2023). In particular, riverine dissolved δ⁷Li values have been shown to have a negative correlation with runoff, because it controls the water-rock interaction time that affects clay formation, with the dry season exhibiting enhanced clay formation and higher δ⁷Li values than the wet season (Wilson et al., 2021; Zhang et al., 2022).

A detailed discussion on the lithium isotope interpretative framework. including the roles of grain size, hydrodynamic sorting, and lithology, is provided in Jaimes-Gutierrez et al. (2025, Supplement).

Figure 3 Processes determining the lithium isotope composition of bulk sediments (δ⁷Li_sed) and dissolved lithium flux along a sediment routing system from source to sink in relation to denudation. The weathering intensity ($W/D$) expresses the relative share of weathering (W) over denudation (D), where $D=W+E$ (erosion). Modified from Tofelde et al. (2021), Pogge von Strandmann et al. (2021b), and Bufe et al. (2024).

During the PETM, Pogge von Strandmann et al. (2021a) documented a ∼ 3 ‰ negative δ⁷Li excursion in several marine carbonate sections, indicating globally enhanced weathering fluxes (50 %–60 %) and erosion rates (2–3×), and a shift to an overall lower weathering intensity regime. At a continental scale, detrital lithium isotope records from North America show coherent negative δ⁷Li excursions in both floodplain and deep-marine settings, indicating rapid propagation of erosion- and weathering-related signals through sediment-routing systems under intensified hydrological conditions, despite largely stable sediment provenance (Ramos et al., 2022; Jaimes-Gutierrez et al., 2025; Rush et al., 2025). On a regional scale, Ramos et al. (2022) reported a rapid, sustained increase in silicate weathering intensity in the Bighorn Basin floodplains that was attributed to seasonal hydrological variability. Similarly, Chen et al. (2023) identified a ∼ 100 % increase in silicate weathering intensity in the Nanyang Basin, East Asia. These studies highlight the roles of local hydrology, lithology, and erosion in shaping regional weathering responses and the associated δ⁷Li changes. However, they also reveal significant gaps in our understanding of how regional processes integrate into driving global δ⁷Li records and carbon cycle feedbacks. Notably, discrepancies between the proposed increases in weathering intensity at a regional scale (Chen et al., 2023; Ramos et al., 2022) and the inferred decrease at a global scale (Pogge von Strandmann et al., 2021a; Jaimes-Gutierrez et al., 2025; Rush et al., 2025) require further assessment of how regional climatic and geological controls translate into weathering responses.

To address these gaps, we focused on the silicate weathering response to the PETM climatic perturbation in two sections of the Southern Pyrenees (Espugafreda and Rin, Fig. 1). With its semi-arid climate, seasonal precipitation, and relatively unreactive lithologies comprising reworked sediments and significant carbonate content (e.g. Eichenseer, 1988; Eichenseer and Luterbacher, 1992; Gómez-Gras et al., 2016), this setting represents a contrasting regional weathering regime to previous PETM studies. Our results contribute to understanding how floodplain paleosols, which are often overlooked in global weathering studies, respond to climatic perturbations, with broader implications for the recovery of Earth's climate system after significant warming events.

2 Geological context

The Pyrenees formed as a result of convergence between the Iberian and European plates, a process that initiated in the Late Cretaceous and continued into the Miocene (Mattauer and Henry, 1974; Roest and Srivastava, 1991; Rosenbaum et al., 2002; Roure et al., 1989). The orogenic evolution began with the mid-Cretaceous hyper-extension of the Iberian margins, followed by the late Cretaceous subduction and collision with the European plate (Teixell et al., 2016). Foreland basins formed on both sides of the fold and thrust belt (Puigdefàbregas and Souquet, 1986; Muñoz, 1992; Gómez-Gras et al., 2016). The Southern Pyrenean foreland basin was active between the Late Cretaceous and the Oligocene, and contains well-preserved sedimentary archives of continental and marine environments.

The Tremp-Graus Basin is located in the South-Central Pyrenean Foreland Basin (Spain), delimited by the Boixols Thrust to the north and the Montsec Thrust to the south (Fig. 1C). During the Palaeocene, the Tremp-Graus Basin was dominated by continental sedimentation sourced from the Pyrenees (Gómez-Gras et al., 2016). The continental deposits of the Thanetian Esplugafreda Formation (Fm.) predominantly represent floodplain sediment accumulation and consist of clay, silt, carbonate nodules, and Microcodium grains, with some isolated sandy to conglomeratic channels (Dreyer, 1993; Puigdefàbregas and Souquet, 1986; Schmitz and Pujalte, 2003, 2007).

2.1 Esplugafreda section

The Esplugafreda section (42°14^′50^′′ N; 0°45^′13^′′ E, Fig. 1B) has been widely studied for its well-preserved Palaeocene-Eocene sedimentary record (Baceta et al., 2005; Basilici et al., 2022; Jaimes-Gutierrez et al., 2024; Khozyem, 2013; Schmitz and Pujalte, 2003; Tremblin et al., 2022) (Fig. 1D). The Upper Thanetian sediments belong to the Esplugafreda Fm. in the Tremp Group of the Tremp-Graus Basin (Dreyer, 1993). This formation consists of coarse-grained stream deposits intercalated with red floodplain sediments that are rich in carbonate nodules and characterised by mature paleosols. The PETM sediments have been classified into five stratigraphic members (Basilici et al., 2022; Colombera et al., 2017; Pujalte et al., 2014; Pujalte and Schmitz, 2005). Member 1 belongs to the Claret Fm. and consists of a fining-upwards sequence of conglomerates and cross-laminated sandstones, known as the Incised Valley Fill (IVF) sediments. During this interval, a first negative carbon isotope excursion (CIE) marks the Pre-Onset Excursion (POE) (Khozyem, 2013; Tremblin et al., 2022). Member 2 at the onset of the Ypresian is represented by the Claret Conglomerate (Pujalte and Schmitz, 2005), a 3–5 m thick conglomeratic unit, corresponding to a braid plain which has been interpreted as the proximal part of a megafan (Schmitz and Pujalte, 2007). Member 3, the Yellowish soils, consists of yellow mudstone with purple mottling, and the main body of the CIE is recorded during this interval (Pujalte and Schmitz, 2005). Member 4, consisting of red soil with gypsum, and Member 5, comprising light red mudstones with scarce carbonate nodules, correspond to the recovery interval of the PETM in this locality (Baceta et al., 2011; Basilici et al., 2022; Khozyem, 2013; Pujalte et al., 2014; Pujalte and Schmitz, 2005; Tremblin et al., 2022).

Figure 4 Rin section between the upper Thanetian and lower Ypresian. Members 1–5 described in the literature for the Esplugafreda section can be identified in the Rin section, except that Member 4 from the recovery phase of the PETM has not been preserved. E. Fm.: Esplugafreda Formation.

2.2 Rin section

The Rin section (42°19^′42.01^′′ N; 0°32^′42.16^′′ E, Figs. 1B and 4) is a Palaeocene-Eocene sequence comprising mudstone-dominated alluvial deposits and very shallow marine carbonate alternations, indicating episodes of transgression and regression on the coastal plain (Schmitz and Pujalte, 2007) (Fig. 1). The upper Esplugafreda Fm. soils are characterised by grey mottling and sparse iron nodules, with preserved pedogenic features such as peds. Member 1, the IVF, consists of 4 m-thick reddish-yellow soils that are rich in carbonate nodules. Member 2, the Claret Conglomerate, outcrops as a 3 m-thick calcareous conglomerate with pale red clay pockets, and has sparse carbonate nodules and charophyte occurrences. Member 3, the Yellowish soils, consists of 13 m-thick reddish-yellow clays and silts. The base of Member 3 records sparse occurrences of lignite and carbonate nodules. Member 4 is not preserved in the Rin section, and the upper 3 m of the sequence consists of Member 5, which has light grey to reddish yellow soils with grey mottling, before the overlying Alveolina Limestone.

3 Material and methods

3.1 Size fraction separation

Standard protocols (e.g. Adatte et al., 1996; Bauer et al., 2016) were followed for decarbonation and size fraction separation at the Institute of Earth Sciences clay laboratory at the University of Lausanne (ISTE-UNIL). Samples (∼ 5 g) were leached with 10 % HCl for 30 min in a bubble bath, including 3 min in an ultrasonic bath, to disaggregate sediments and dissolve calcite. Distilled water was used to remove the acid until a neutral pH was obtained. Subsequently, the <2 µm fraction was separated by settling and enhanced with a centrifuge. Settling and extraction were repeated three times.

3.2 Clay mineralogy

The clay minerals were identified on air-dried and ethylene glycol-solvated samples at ISTE-UNIL following the protocol described in Adatte et al. (1996). An aliquot of the separated size fractions was pipetted on glass slides and dried at room temperature. The air-dried samples were further analysed with a Thermo Scientific ARL X'TRA powder diffractometer equipped with a Cu anode, operated at 45 kV and 40 mA. The step size was 0.02°, with a scan rate of 0.5–1.2° min⁻¹. Samples were glycolated to identify smectite (Moore and Reynolds, 1997). Diffractograms were analysed using the XRDWin software, where the background was removed, and a deconvolution was performed for overlapping peaks (e.g. K002 and Ch004).

3.3 Nodule purification

Carbonate nodules were washed with running distilled water until visible clay clumps were removed. They were then placed in a beaker with distilled water and in an ultrasonic bath to remove the remaining clay particles. A second round in the ultrasonic bath was then carried out with some drops of 10 M HCl in order to remove the outermost layer. The nodules were later washed in running distilled water, dried at 40 °C, and ground.

3.4 Rock-Eval pyrolysis

Organic matter analyses were performed on powdered bulk rock samples using a Rock-Eval 6 at ISTE-UNIL, following standard methodology (Behar et al., 2001; Espitalie et al., 1985). For calibration, the IFP 160000 standard was used. The Rock-Eval pyrolysis parameters measured were hydrogen index (HI, mg HC/g TOC, HC = hydrocarbons), oxygen index (OI, mg CO₂/g TOC), Tmax (°C), and total organic carbon content (TOC, wt %). The HI, OI, and Tmax values give an overall measurement of the type and degree of maturation of the organic matter (e.g. Espitalie et al., 1985).

3.5 Isotope geochemistry

3.5.1 Organic matter carbon isotopes

The carbon isotope composition of the decarbonated bulk rock samples was determined at the Institute of Earth Surface Dynamics at the University of Lausanne (IDYST-UNIL) by elemental analysis/isotope ratio mass spectrometry (EA/IRMS). The EA/IRMS system consisted of 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 Thermo Fisher Scientific, Bremen, Germany) operated under continuous helium (He) flow (Spangenberg, 2006; Spangenberg and Zufferey, 2019). The carbon isotope compositions were reported in the delta (δ) notation as permil (‰) variations of the molar ratio of the heavy to light isotope (¹³C $/$ ¹²C) relative to the international standard Vienna Pee Dee Belemnite limestone (VPDB). For calibration and normalisation of the measured δ¹³C values to the Vienna Pee Dee Belemnite limestone (VPDB) standard, a four-point calibration was used with international reference materials and in-house standards (Spangenberg and Zufferey, 2019). The used standards included UNIL-Glycine (δ¹³C $=-\mathrm{26.10}\pm \mathrm{0.05}$ ‰), UNIL-Urea-1 (δ¹³C $=-\mathrm{43.00}\pm \mathrm{0.04}$ ‰), UNIL-Pyridine (δ¹³C $=-\mathrm{29.25}\pm \mathrm{0.06}\mathrm{‰}$), and the RM USGS24 graphite (δ¹³C $=-\mathrm{16.05}\pm \mathrm{0.04}$ ‰). Analyses were done in duplicates. The accuracy of the analyses was checked periodically through the analysis of international RM standards not used for calibration. The reproducibility and precision of the EA/IRMS δ¹³C analyses were determined by the standard deviation of separately replicated analyses and were better than 0.1 ‰.

3.5.2 Lithium isotopes

Sample digestion, column chemistry, and mass spectrometry were conducted in the London Geochemistry and Isotope Centre (LOGIC) laboratories at University College London (UCL) and Birkbeck, University of London. Clay samples were subjected to bulk digestion using concentrated HF, HNO₃, and HClO₄ in Teflon beakers on a hot plate at 130 °C, followed by steps in concentrated HNO₃ and 6 M HCl. The carbonate nodules were subject to leaching to separate the carbonate and detrital fractions. The carbonate fraction was extracted by leaching ∼ 100 mg of sample in 8 mL 0.1 M HCl for 1 h (Pogge von Strandmann et al., 2013; Wilson et al., 2021), allowing a maximum of ∼ 40 mg of calcium carbonate to be dissolved.

A standard method of elution was applied for lithium isotope separation in 0.2 M HCl. Two-column passes were applied through AG50W-X12 resin to ensure matrix removal (Pogge von Strandmann et al., 2013). Given that lithium isotopes are fractionated during ion chromatography, sample splits were collected before and after the lithium collection interval to assess column yields. For example, a 1 % loss in yield at UCL has been assessed to lead to an offset of 1.7 ‰ (Wilson et al., 2021). Here, yields between two column passes were 99.8 %–100 %, indicating excellent recovery.

Lithium isotope measurements were performed on a Nu Plasma 3 MC-ICP-MS at UCL, using a Cetac Aridus 2 desolvation system, “super-lithium” cones, and standard-sample bracketing with the IRMM-016 Li standard (Pogge von Strandmann et al., 2019). Samples were measured at least three times within an analytical session, with each measurement integrating ∼ 50 s, and the reported values are the mean and standard deviation (2 SD) of these values, given in permil (‰) relative to the IRMM-016 standard. Accuracy and external reproducibility were assessed using seawater and USGS standard BCR-2, which gave δ⁷Li values of $+\mathrm{31.3}\pm \mathrm{0.6}\mathrm{‰}$ (2 SD, n=28 ) and $+\mathrm{2.5}\pm \mathrm{0.3}\mathrm{‰}$ (n=5), respectively.

3.5.3 Neodymium isotopes and rare earth element concentrations

After decarbonation using 10 % HCl for 30 min, clays were separated from decarbonated sediments into <0.5 and 0.5–2 µm fractions (analytical protocol for size fraction separation reported in Jaimes-Gutierrez et al., 2024). A total of 18 samples (8 in the <0.5 µm size fraction and 10 in the 0.5–2 µm size fraction) were analysed for their neodymium (Nd) isotopic composition and their Nd and samarium (Sm) concentrations. Aliquots of about 1.5 mg of each clay fraction followed a sequential leaching procedure to remove Fe-Mn oxides and organic matter, based on the protocol of Bayon et al. (2002) and Gutjahr et al. (2007), slightly adapted. The Fe-Mn oxides were removed using a solution of 0.5 M hydroxylamine hydrochloride in 20 % $v/v$ acetic acid for 48 h. Then, the organic matter was removed with a 5 % H₂O₂ solution for 48 h.

The leached samples were dried and digested by alkaline fusion following the protocol of Bayon et al. (2009), along with certified standards (BHVO-2, BRC-2) from the United States Geological Survey (USGS). Approximately 50 mg of each sample underwent alkaline fusion in a carbon crucible with 0.6 g of NaOH and 1.2 g of Na₂O₂ heated at 650 °C for 12 min in a furnace, before adding ultrapure water in which Fe-hydroxides precipitated, concentrating rare earth elements. After centrifugation, the samples were dissolved in 3 mL 4 M HCl.

From this solution, an aliquot of 0.3 mL was extracted for analyses of Nd and Sm concentrations. Part of the samples were measured for their Nd and Sm concentrations on an Agilent 7500 quadrupole ICP-MS spectrometer in the Laboratoire Magmas et Volcans (LMV) in Clermont-Ferrand (France), and quantified using standard bracketing with a solution of BHVO-2 during the session. Accuracy and reproducibility were assessed using two BHVO-2 and one BCR-2 samples among the samples. Deviations of Nd and Sm concentrations from these standards were below 11 %. The other part of the samples was measured for their Nd and Sm concentrations on a Thermo Scientific X-Series II^® at the Pole Spectrométrie Océan in Brest (France), and quantified using multi-element calibration standards prepared from single element standards purchased from SCP science (Baie d'Urfé, Québec, Canada). Accuracy and reproducibility were assessed using one BHVO-2 and one BCR-2 sample, which were analysed among the samples. Deviations of Nd and Sm concentrations from these standards were below 7 %.

Purified neodymium fractions were isolated from the mother solution by ion chromatography following the protocol described in Gaitan et al. (2023) for the low-pressure, automated column chromatography PrepFAST-MC^® system device, using AG50W-X8 (200–400 mesh) resin for rare earth element separation and Ln Spec (50–100 µm) resin for Nd separation. Part of the neodymium isotopic measurements was performed on a MC-ICP-MS Neptune Plus (Thermo Scientific) at the Laboratoire Magmas et Volcans in Clermont-Ferrand (France). Ratios were corrected for mass bias using an exponential law and a ¹⁴³Nd $/$ ¹⁴⁴Nd ratio of 0.7219. Mass-bias-corrected ¹⁴³Nd $/$ ¹⁴⁴Nd were normalised to a JNdi-1 value of 0.512115 (Tanaka et al., 2000). Repeated measurements of JNdi-1 throughout the session gave an external reproducibility of ±0.000009 (2σ, n=15), corresponding to ±0.18 in the standard ε_Nd(0) notation. Analyses of two BHVO-2 reference materials yielded a ¹⁴³Nd $/$ ¹⁴⁴Nd ratio of 0.512991±0.000005 for each, in excellent agreement with the published value of 0.512990±0.000010 (Weis et al., 2005). The other part of the samples was analysed on an MC-ICP-MS Neptune Plus (Thermo Scientific) at the ENS of Lyon (France). Ratios were corrected for mass bias using an exponential law and a ¹⁴³Nd $/$ ¹⁴⁴Nd ratio of 0.7219. Mass-bias-corrected ¹⁴³Nd $/$ ¹⁴⁴Nd values were normalised to a JNdi-1 value of 0.512115 (Tanaka et al., 2000). Repeated measurements of JNdi-1 throughout the session gave an external reproducibility of ±0.000018 (2σ, n=16), corresponding to ±0.34 in the standard ε_Nd(0) notation. Analyses of four BHVO-2 reference materials gave an average ¹⁴³Nd $/$ ¹⁴⁴Nd ratio of 0.512985±0.000009 for each, in agreement with the published value of 0.512990±0.000010 (Weis et al., 2005).

The data are reported in the standard epsilon notation ε_Nd = [((¹⁴³Nd $/$ ¹⁴⁴Nd)_sample/(¹⁴³Nd $/$ ¹⁴⁴Nd)_CHUR) −1] ×10⁴), corrected for the radioactive decay of ¹⁴⁷Sm to ¹⁴³Nd based on the Nd and Sm concentrations measured for each sample (¹⁴⁷Sm $/$ ¹⁴⁴Nd = Sm/Nd ×0.6049), an age of 55.8 Ma, and the ¹⁴⁷Sm radioactive decay constant λ ($\mathrm{6.54}\times {\mathrm{10}}^{-\mathrm{12}}$ yr⁻¹; Lugmair and Marti, 1978). The CHUR (CHondritic Uniform Reservoir) ¹⁴³Nd $/$ ¹⁴⁴Nd ratio was also corrected using a ¹⁴⁷Sm $/$ ¹⁴⁴Nd ratio of 0.1960 and a present-day value of 0.512630 (Bouvier et al., 2008).

Figure 5 Rin section isotopes and clay mineralogy. (A) The δ¹³C_OM record shows the negative Carbon Isotope Excursion (CIE), with an onset before the Claret Conglomerate and sustained negative values until the Alveolina Limestone. (B) The CIE was accompanied by a positive excursion in lithium isotopes of the clays (δ⁷Li_clays), reaching a 0.9 ‰ excursion in the Yellowish soils member. (C) Lithium isotopes in the carbonate nodules showed high variability and a less conclusive trend. (D) Rin section clay mineralogy. The pre-PETM and body intervals were determined based on the δ¹³C_OM record and the stratigraphy. Dashed lines in panel (A) represent average values for δ¹³C_OM in the pre-PETM and syn-PETM intervals. Grey bars in panel (B) outline the analytical uncertainty (2 SD) of δ⁷Li_clays.

4 Results

4.1 Clay mineralogy in the Rin section

The clay mineralogy in the Rin section (Fig. 5D, Table S1 in the Supplement) comprises mixed-layer illite-smectite (I/S), with a mean abundance of 34±11 (1σ) wt %; illite, 29±9 wt %; kaolinite, 30±11 wt %; and minor chlorite, 7±5 wt %. The pre-PETM samples have a mixed-layer I/S abundance of 42±10 wt %, decreasing to 32±10 wt % during the body of the PETM. The kaolinite abundance increases from 22±14 wt % to 32±9 wt % during the PETM body, while the illite and chlorite abundances remain stable.

4.2 Organic matter carbon isotopes in the Rin section

Throughout the section, the mean δ¹³C_OM value is −25.7 ‰, with a standard deviation (1σ) of 1.2 ‰ (Fig. 5A, Table S2). The pre-PETM samples, between 5.8 and 9.4 m, have a mean value of $-\mathrm{23.5}\pm \mathrm{0.5}\mathrm{‰}$. A negative excursion begins in samples at 10.7 and 11.4 m, with values decreasing to $-\mathrm{24.6}\pm \mathrm{0.1}\mathrm{‰}$. The most depleted values occur between 12.8 and 38.9 m, with a mean of $-\mathrm{26.3}\pm \mathrm{0.5}\mathrm{‰}$. The final sample at 39.5 m, below the Alveolina Limestone, suggests a return to pre-PETM levels, with a value of −23.7 ‰.

4.3 Lithium isotopes in the Rin and Esplugafreda sections

The clays of the Rin section have a mean lithium isotope composition of $-\mathrm{2.9}\pm \mathrm{0.5}\mathrm{‰}$ (1σ) (Fig. 5B, Table S2). Between 5.8 and 18.2 m, the mean composition is $-\mathrm{3.4}\pm \mathrm{0.2}\mathrm{‰}$. Above this, from 20.3 to 39.5 m, the mean composition is $-\mathrm{2.6}\pm \mathrm{0.2}\mathrm{‰}$, which corresponds to a shift towards more positive values of ∼ 0.8 ‰. The minimum value of −3.7 ‰ is seen before the Claret Conglomerate, and the maximum value of −2.2 ‰ occurs immediately after the Claret Conglomerate, indicating a total range of up to ∼ 1.5 ‰. The δ⁷Li values measured on carbonate nodules have a maximum value of 8.6 ‰, a minimum value of −3.0 ‰, and a mean composition of 0.9±4.0 ‰ (Fig. 5C). No clear temporal trend is observed in the δ⁷Li_nodules record. No systematic correlation is observed between δ⁷Li_clays and the relative abundance of individual clay minerals within analytical uncertainty (Fig. S1 in the Supplement).

Figure 6 Esplugafreda section isotopes and clay mineralogy. (A) The δ¹³C record from microcrystalline carbonate nodules (from Khozyem, 2013) shows a negative carbon isotope excursion during the Pre-Onset Excursion (POE) and during the main body of the PETM. (B) Positive excursions in lithium isotopes of the clays (δ⁷Li_clays) during both the POE and the main CIE of the PETM. (C) Neodymium isotopes (ε_Nd (t=55.8 Ma)) show no variation throughout the section, indicating constant provenance. (D) Esplugafreda section clay mineralogy (modified from Jaimes-Gutierrez et al., 2024).

At the Esplugafreda section, the clays have a mean lithium isotope composition of $-\mathrm{3.7}\pm \mathrm{0.7}\mathrm{‰}$ (Fig. 6B, Table S3). The pre-PETM samples (0–10 and 21–28 m, Jaimes-Gutierrez et al., 2024 and references therein) have a mean composition of $-\mathrm{3.8}\pm \mathrm{0.2}\mathrm{‰}$. The POE samples (∼ 15–21 m) have a composition of $-\mathrm{3.2}\pm \mathrm{0.2}\mathrm{‰}$, and the syn-PETM sediments have a composition of $-\mathrm{3.0}\pm \mathrm{0.2}\mathrm{‰}$, indicating a PETM shift towards more positive values of ∼ 0.8 ‰. The post-PETM sediments have a composition of $-\mathrm{4.5}\pm \mathrm{0.4}\mathrm{‰}$. Samples from the pre-PETM, Pre-Onset Excursion (POE), and syn-PETM intervals display relatively higher δ⁷Li_clays values, whereas recovery-phase samples form a distinct cluster characterised by lower δ⁷Li_clays (Fig. S2).

4.4 Neodymium isotopes in the Esplugafreda sections

Throughout the Esplugafreda section, the 0.5–2 µm clays have a mean ε_Nd (t=55.8 Ma) composition of $-\mathrm{10.94}\pm \mathrm{0.16}$ (2σ), while the <0.5 µm fraction has a mean composition of $-\mathrm{10.88}\pm \mathrm{0.11}$ (Fig. 6C, Table S4). In comparison, the typical analytical uncertainty on any individual sample measurement was 0.21 (2σ). Hence, these values are considered constant through time, with no deviation significantly outside the analytical uncertainty. The neodymium isotope measurements on the <0.5 µm fraction are also indistinguishable from those on the 0.5–2 µm fraction (Fig. 6C).

5 Discussion

5.1 The PETM in the Rin section

The PETM sediments in the Rin section represent an archive of the climatic perturbation in a coastal terrestrial setting (Prieur et al., 2025; Pujalte et al., 2014). This locality records a negative δ¹³C_OM excursion of −2.8 ‰ from pre- to syn-PETM (Fig. 5), in agreement with the CIE excursion of 3 ‰–5 ‰ identified in other southern Pyrenean sections and other global settings (e.g. Schmitz et al., 2001; Schmitz and Pujalte, 2007; McInerney and Wing, 2011; Pujalte et al., 2015). The slightly reduced magnitude compared to the global record is consistent with observed systematic differences in the CIE across different types of terrestrial archives, with paleosol carbonates typically recording a 1 ‰–2 ‰ larger CIE than paleosol organic matter (Bowen et al., 2004; Cotton et al., 2015; Gallagher et al., 2019). We do not identify the POE in the Rin section, and we suggest that it may have been missed due to its occurrence further down in the section. Likewise, the recovery to pre-PETM values is also largely absent, with just one sample below the Alevolina Limestone showing less depleted δ¹³C_OM values.

The five members recognised in the Claret Fm. show an evolution from the eastern terrestrial setting into the western marine domain. At Esplugafreda (Fig. 6), the five members are recorded (Basilici et al., 2022 and references therein), including Member 1, IVF (pre-PETM); Member 2, the Claret Conglomerate (Onset at Esplugafreda); Member 3, Yellowish soils (syn-PETM, or body of the PETM); Member 4, red paleosols with gypsum; and Member 5, consisting of red mudstones with carbonate nodules (e.g. Schmitz and Pujalte, 2007; Baceta et al., 2011; Pujalte et al., 2014; Colombera et al., 2017; Basilici et al., 2022). However, Member 4, the gypsum-rich member, only occurs in the eastern part of the basin (Pujalte et al., 2014). Given the coastal position of the Rin section, at the marine-continental transition and only ∼ 20 km east of the Serraduy section, representing the westernmost expression of the interfingering between continental deposits from the Esplugafreda Fm. and marine carbonates (Prieur et al., 2025), the absence of Member 4 supports a further downstream position of the Rin section relative to the Esplugafreda floodplain section.

Duller et al. (2019) estimated a lag time of approximately 16.5±7.5 kyr between the CIE and the onset of coarse-grained deposition at terrestrial sites in the Pyrenees. While sections such as Tendrui, Claret, and Campo (Domingo et al., 2009; Pujalte et al., 2009) display a stratigraphic offset consistent with this lag, the Esplugafreda section does not show such an offset (Duller et al., 2019). In the Rin section discussed here, we observe a clear offset between the onset of the CIE and the arrival of the Claret Conglomerate (Fig. 5). To correctly position the lag time and explore the missing POE, future work could focus on high-resolution δ¹³C characterisation of the section.

The clay mineralogy at Rin further suggests a potential signal propagation effect. A shift from smectite-dominated clays during the pre-PETM interval to an increase in kaolinite during the syn-PETM interval (Fig. 5D) could indicate a transition to more hydrolysing conditions and to an increase in weathering intensity, or enhanced erosion of former sedimentary formations rich in kaolinite. However, this trend also corresponds to a downstream transition from authigenic smectite-rich paleosols at Esplugafreda (Fig. 6D) (e.g. Khozyem, 2013; Basilici et al., 2022; Jaimes-Gutierrez et al., 2024) to kaolinite-dominated sediments in the Zumaia deep-marine section (Bolle and Adatte, 2001; Gawenda et al., 1999; Schmitz et al., 2001). This mineralogical gradient from Esplugafreda to Rin and Zumaia underscores the system connectivity across the basin (Pujalte et al., 2014). However, such variations in clay composition may also reflect differences in sediment provenance (although temporal changes are not observed in the ε_Nd record from Esplugafreda); differential mineral transport; enhanced floodplain weathering; an increased proportion of eroded sedimentary formations downstream, bringing reworked kaolinite (Pujalte et al., 2015); or a larger catchment area feeding the marine system (Chamley, 1989 and references therein).

An influence from Marine Authigenic Aluminosilicate Clay (MAAC) formation is feasible in the coastal setting of the Rin section, but is unlikely to have been a significant driver of changes in the lithium isotope record given the dominance of continental inputs. In marine sediments, clay mineral assemblages dominated by kaolinite, illite, smectite, and mixed-layer illite/smectite are generally interpreted as mainly detrital in origin, reflecting continental weathering and fluvial transport rather than in situ marine precipitation (Fagel, 2007; Thiry, 2000; Velde and Meunier, 2008). Although the authigenic formation of these clay minerals is known to occur in marine environments, it is typically restricted to specific conditions or processes, such as the alteration of volcanic ash, or evaporitic or hydrothermal settings, and is therefore not expected to be significant in shallow, nearshore depositional systems with high terrigenous sedimentation rates (Meunier, 2005; Środoń, 2001; Wise et al., 2001). In contrast to glauconite or other green clays, which may form authigenically under low sedimentation rates, kaolinite and illite in coastal marine settings are widely regarded as inherited from continental sources (Bernhardt et al., 2020; Fagel, 2007; Meunier, 2005; Presti and Michalopoulos, 2008; Thiry, 2000). Given the dominantly continental depositional setting of the Rin section and the large detrital clay input inferred for the nearshore environment, any MAAC contribution is expected to have been minor relative to the terrigenous signal. MAAC form from isotopically heavy seawater or porewaters and are therefore expected to have δ⁷Li values substantially higher than detrital clays (Pogge von Strandmann et al., 2021b). A simple mass balance indicates that even under extreme assumptions (i.e. 2 wt % MAAC, 10 ‰ δ⁷Li values), such a contribution would shift bulk δ⁷Li_clays values by <0.3 ‰, far smaller than the observed PETM excursion.

Figure 7 Comparison of isotopic results from the Rin and Esplugafreda sections. Grey solid lines show the raw data, while dashed black and red lines indicate typical values for different intervals, where red represents the syn-PETM. The pre-PETM and post-PETM shifts are indicated with the black arrows, and the magnitudes of the δ⁷Li_clays excursions are reported.

5.2 Evolution of weathering intensity in the continental realm of the Southern Pyrenees

Our δ⁷Li_clays records from Rin and Esplugafreda both show a positive (∼ 1 ‰) lithium isotope excursion in the continental Southern Pyrenees during the onset and body of the PETM (Fig. 7). The δ⁷Li values from carbonate nodules at Rin show greater variability (Fig. 5C), but remain inconclusive due to potential clay contamination or cation exchange between clays and carbonates (e.g. Pogge von Strandmann et al., 2019). Given the high Li content in silicate minerals, even a minor clay particle content in the nodules could contaminate the carbonate signature. Critically, the invariant ε_Nd composition of both size fractions throughout the Esplugafreda record (Fig. 6C) supports a constant provenance of the sediments, which suggests that the δ⁷Li_clays records can be reliably interpreted as a reflection of weathering regime changes in response to the climatic perturbation.

In many other PETM records, both marine and terrestrial δ⁷Li values show a negative excursion from pre-PETM to syn-PETM conditions. Pogge von Strandmann et al. (2021a) documented a negative δ⁷Li excursion of ∼ 3 ‰ during the PETM in several marine carbonate sections and in detrital shales, indicating intensified global erosion rates (by 2–3×) and a 50 %–60 % increase in silicate weathering fluxes, which was proposed to have contributed to climate stabilisation. Ramos et al. (2022) also found a negative δ⁷Li_clays excursion, albeit with a smaller magnitude of ∼ 1.5 ‰, in fine sediments of the Bighorn Basin, North America, during the PETM, and this change was sustained during the recovery phase. Marine records off the North American margin also reflect a negative excursion in the Gulf of Mexico (0.4 ‰–1.5 ‰) and the Mid-Atlantic Coast (∼ 3 ‰) (Jaimes-Gutierrez et al., 2025; Rush et al., 2025). Consistent with these results, Chen et al. (2023) found a negative δ⁷Li_clays excursion of ∼ 3 ‰ in the Nanyang Basin, East Asia, which, together with the negative δ⁷Li excursion in the lacustrine carbonates, was interpreted as recording a doubling of the regional silicate weathering intensity during the PETM.

Why, then, do the Southern Pyrenees floodplains record a positive δ⁷Li_clays excursion during the PETM? We interpret the positive δ⁷Li_clays excursion during the PETM as a shift towards increased incongruent weathering (moderate weathering regime), characterised by enhanced clay formation in the floodplain deposits. This regime would be characterised by increased chemical weathering, but relatively greater increases in physical erosion and sediment transport (e.g. Pujalte et al., 2015; Chen et al., 2018; Prieur et al., 2024), due to enhanced runoff causing short water residence times and rapid sediment export (i.e. lower $W/D$; Fig. 8). Hydrological changes during the PETM are widely documented in the Southern Pyrenees (Chen et al., 2018; Jaimes-Gutierrez et al., 2024; Prieur et al., 2024, 2025; Pujalte et al., 2015; Rush et al., 2021; Schmitz and Pujalte, 2007) and likely played a central role in driving those weathering changes. A shift towards more intense, episodic rainfall, without an increase in mean annual precipitation (Rush et al., 2021), could have reduced infiltration, increased runoff, and shortened water–mineral interaction times. These processes would lead to decreased clay formation in the uplands by shortening water-rock interaction times (Kump et al., 2000; Riebe et al., 2004), while shifting clay production and accumulation towards lowland environments, where longer sediment residence times promote authigenic clay formation, consistent with modern river floodplain processes (e.g. Dellinger et al., 2015; Maffre et al., 2020). We also note that weathering processes are highly heterogeneous, and therefore global variability during the PETM can be expected (e.g. Frings, 2019). Critically, a decrease in $W/D$ can lead both to positive or negative δ⁷Li excursions, depending on the starting weathering regime (Krause et al., 2023).

Increasing evaporation, as recorded by gypsum lenses in Esplugafreda (Baceta et al., 2011; Khozyem, 2013; Jaimes-Gutierrez et al., 2024 and references therein), could also result in oversaturated pore waters, favouring clay formation. Experimental and field-based studies show that enhanced evaporation and reduced water availability increase the dissolved δ⁷Li values in soil and pore waters, leading to higher δ⁷Li compositions of clays forming in equilibrium with these fluids (Pogge von Strandmann et al., 2023; Xu et al., 2022). Hydrological controls can also influence the δ⁷Li_clays signatures by regulating water-rock interactions, sediment transport, and secondary mineral formation (Fig. 8). In supply-limited weathering regimes, increased runoff shortens water–rock interaction times and lowers dissolved δ⁷Li values in river waters (Zhang et al., 2022), while enhanced sediment transfer promotes sediment storage and prolonged water–sediment interaction in floodplains, where continued clay alteration and isotopic re-equilibration can yield relatively higher δ⁷Li values in the clay fraction preserved in lowland deposits. Overall, rapid sediment transport limits basin-scale weathering intensity, despite localised clay formation in floodplain environments.

Figure 8 Weathering regime change from pre-PETM to syn-PETM based on δ⁷Li_clays in the continental deposits of the Southern Pyrenees. The floodplain records were characterised by a decrease in weathering intensity and an increase in clay formation. Enhanced transport efficiency resulted in a major increase in physical erosion, with reaction kinetics limiting chemical weathering. Modified from Dellinger et al. (2015, 2017). At low $W/D$, the empirical river sediment mixing line includes the effects of mixing with primary minerals.

In the Esplugafreda section, the recovery phase sees a shift towards more negative values than the pre-PETM conditions. During this interval, we observe a coeval increase in illite, kaolinite, and chlorite abundances in the Esplugafreda section (Jaimes-Gutierrez et al., 2024). This shift in the clay mineral assemblage strongly suggests an increase in detrital input. For this reason, we consider that the observed excursion cannot be straightforwardly interpreted as a change in weathering intensity, but rather reflects the clay mineral assemblage (Fig. S2).

Temperature effects on lithium isotope fractionation are minor in the Southern Pyrenees. While lithium isotope fractionation during clay formation is temperature-dependent (Li and West, 2014; Vigier et al., 2008), the fractionation factor (α) for incorporation in smectite is nearly constant across typical surface weathering temperatures (Vigier et al., 2008). Using the ∼ 3 °C warming estimated for the continental Pyrenees across the PETM (Jaimes-Gutierrez et al., 2024), the maximum temperature-driven change in clay δ⁷Li values is expected to be only a few tenths of a per mil (Li and West, 2014). This effect would shift clay δ⁷Li values slightly towards heavier values and could therefore contribute marginally to the observed excursion, but it is insufficient to explain the full magnitude of the ∼ 1 ‰ positive δ⁷Li_clays shift observed in the records. Hence, climatic and hydrological processes, rather than direct temperature effects, must dominate the δ⁷Li signal.

The differences in δ⁷Li values between river water and bedrock are controlled by the balance between lithium release by mineral dissolution and lithium removal by secondary mineral formation (Bouchez et al., 2013). The Southern Pyrenees during the PETM was a relatively high-erosion regime, such that physical erosion dominated, increasing the sediment supply and exposing fresh minerals (Chen et al., 2018; Prieur et al., 2024, 2025; Pujalte et al., 2015, 2016; Schmitz and Pujalte, 2007). When erosion exposes fresh minerals, weathering rates increase with total denudation, albeit less strongly than the increases in erosion, consistent with shared controls on chemical weathering and physical denudation rates (Riebe et al., 2004; West et al., 2005). Even though high-relief regions produce weakly weathered sediments, their high sediment yields and moderate clay formation rates result in elevated weathering fluxes (Gaillardet et al., 1999). Therefore, we propose that the Southern Pyrenean floodplains record a shift from a high-weathering intensity regime to a moderate-weathering intensity regime during the PETM (Fig. 8). The pre-PETM conditions were characterised by a low reactivity of the parent lithology (Caves Rugenstein et al., 2019; Kump and Arthur, 1997), associated with the carbonate-rich, reworked sediments in the floodplain deposits, and hence low total weathering fluxes. The above scenario is also consistent with the “system-clearing” event documented in western North America (Foreman et al., 2012), where sediment transport surged in response to rapid climatic forcing, as well as with other Eocene warming events such as the Mid-Eocene Climatic Optimum, which saw a shift towards enhanced clay formation and a lower weathering intensity (Krause et al., 2023).

The progressive increase in kaolinite content from Esplugafreda to Rin may reflect an evolving weathering signal during sediment transport from the hinterland towards the coastal plains, which potentially extended into the marine realm. This scenario supports a basin-wide connectivity between climate-driven terrestrial processes and marine sedimentary records. In addition, the increase in kaolinite content from Esplugafreda to Rin supports a shift in clay formation processes during the PETM. Kaolinite is typically associated with intense leaching and more advanced weathering, often forming under warm, humid, and periodically saturated conditions (Chamley, 1989; Velde and Meunier, 2008). Hence, its enrichment suggests either intensified in-situ clay formation in the floodplains or increased transport of weathered material from the uplands to the lowlands. In either case, this shift implies greater clay mineral production, consistent with a more incongruent weathering regime driving the observed positive δ⁷Li_clays excursion. Alternatively, enhanced kaolinite supply from the erosion-driven exhumation of older sediments cannot be ruled out based on the current evidence.

Despite these insights, key questions remain unresolved. In particular, a comprehensive study of the provenance and evolution of clay mineralogy is still needed to determine to what extent the observed patterns along the sediment routing system reflect changes in weathering intensity, differential mineral transport, or sediment reworking. Equally important is the need to constrain the precise age of the clay formation in relation to the timing of the different phases of the PETM, which is critical for reconstructing the temporal dynamics of the weathering regime in the Southern Pyrenees. Addressing these gaps will be crucial for better understanding how continental weathering systems responded to extreme climatic perturbations in the past and how they may behave under future global warming scenarios.

6 Conclusions

We explored the silicate weathering response to the PETM in two terrestrial sections from the Tremp-Graus Basin of the Southern Pyrenees. These floodplain records show a positive δ⁷Li_clays excursion, contrasting with the commonly observed global negative δ⁷Li excursion in clays and carbonates. We interpret this excursion as reflecting a shift towards a moderate-intensity, incongruent weathering regime from an initial high-intensity, supply-limited regime. The high erosion rates associated with increased extreme rainfall events and channel mobility may have been the central factor influencing sediment residence times, with rapid sediment transport limiting the extent of chemical weathering. Nevertheless, the elevated denudation rates would have led to higher sediment and dissolved cation fluxes to the ocean, thereby enhancing regional CO₂ drawdown.

We explored two potential archives for recording continental weathering processes using lithium isotopes. The clay records show a distinct response, reflected in positive δ⁷Li_clays excursions synchronous with the negative CIE. However, the δ⁷Li_nodules signal recorded in the carbonate nodules is less conclusive, and we interpret the high temporal variability as a sign of potential contamination by clays in the nodules. Given that Li concentrations in silicate minerals are higher than in carbonates by several orders of magnitude, even minor amounts of clays could have resulted in a mixed response in the nodules. Future studies should explore weaker leaching approaches on such nodules and seek to validate such data with major and trace element analyses.

Provided coeval formation, the increase in kaolinite content from Esplugafreda to Rin provides mineralogical support for more hydrolysing conditions and clay formation during the PETM in the Tremp-Graus Basin, reinforcing the interpretation of more incongruent weathering under altered hydroclimatic conditions. Notably, the parent material in these floodplain paleosols is carbonate-rich and relatively unreactive. These results highlight the critical role of hydrological controls, especially rainfall intensity, runoff dynamics, and sediment residence time, in shaping continental weathering responses during extreme climate events.

Finally, we propose that to fully quantify weathering dynamics during the PETM in the Southern Pyrenees, further work is needed to: (1) constrain the chronology of clay formation; (2) trace the evolution of clay mineralogy and provenance from source to sink; and (3) integrate continental and marine weathering records across the sediment routing system. Together, these steps will be essential for refining our understanding of weathering behaviour and the associated climate feedbacks under rapid climatic perturbations, and for improving predictions of Earth's surface processes in semi-arid floodplain systems in future global warming scenarios.