# Paleoceanography and Paleoclimatology®

## RESEARCH ARTICLE

10.1029/2025PA005121

**Key Points:**

* $\delta^{18}\text{O}$ values of soil carbonates ($\delta^{18}\text{O}_c$) from Teruel, Spain are used to reconstruct changes in rainfall from 4.1 to 2.4 Ma

* $\delta^{18}\text{O}_c$ decreased with cooling during the Marine Isotope Stage M2 and during onset and intensification of Northern Hemisphere Glaciation

* Multi-proxy evidence suggests an establishment of seasonality and summer-specific drying in the Mediterranean Basin beginning at ~3.3 Ma

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

N. Czwakiel, nicole.czwakiel@austin.utexas.edu

**Citation:**

Czwakiel, N., Gallagher, T., Serach, L., Ludvigson, G., Gao, P., Nie, J., et al. (2025). The onset of summer drought on the Iberian Peninsula following Pliocene global cooling. *Paleoceanography and Paleoclimatology, 40*, e2025PA005121. https://doi.org/10.1029/2025PA005121

Received 25 JAN 2025
Accepted 12 SEP 2025

## The Onset of Summer Drought on the Iberian Peninsula Following Pliocene Global Cooling

<figure><img src="image_1.png" /></figure>

N. Czwakiel¹ ID, T. Gallagher¹,² ID, L. Serach¹ ID, G. Ludvigson³ ID, P. Gao⁴ ID, J. Nie⁴ ID, J.-P. Suc⁵ ID, J. Da¹,⁶ ID, and D. O. Breecker¹ ID

¹The Department of Earth and Planetary Sciences, University of Texas at Austin, Austin, TX, USA, ²Now at The Department of Earth Sciences, Kent State University, Kent, OH, USA, ³Kansas Geological Survey, University of Kansas, Lawrence, KS, USA, ⁴School of Earth Science, Lanzhou University, Lanzhou, China, ⁵Sorbonne Université, CY Cergy Paris Université, CNRS, Institut des Sciences de la Terre de Paris, ISTE, Paris, France, ⁶Now at The Department of Geosciences, Penn State University, University Park, TX, USA

**Abstract** The highly biodiverse Iberian Peninsula has winter-dominated precipitation and summer drought. The onset of this modern climate regime, termed a Mediterranean-Type Climate, began during the Late Pliocene–Early Pleistocene, concurrent with global cooling. To advance our understanding of regional hydroclimate changes in relation to past global changes, we present a stable oxygen isotope record measured on paleosol carbonates ($\delta^{18}\text{O}_c$) from Teruel, Spain spanning 4.1–2.4 Ma, when the regional ecosystems shifted from subtropical requiring high humidity to dominantly drought-tolerant. The $\delta^{18}\text{O}_c$ values show a negative excursion during Marine Isotope Stage M2 (~3.3 Ma) and a subsequent decrease during the onset and intensification of the Northern Hemisphere Glaciation (iNHG: ~2.73 Ma). While the $\delta^{18}\text{O}_c$ decrease during global cooling could have resulted from less high-$\delta^{18}\text{O}$ summer rainfall and/or intensified low-$\delta^{18}\text{O}$ winter rainfall, comparison with previously published proxy records point toward drying in the Mediterranean Basin during the Late Pliocene. We therefore propose that global cooling at this time resulted in lower Atlantic and Mediterranean sea surface temperatures during the iNHG which decreased moisture availability to the Mediterranean Basin, and together with enhanced and equatorward shifted Hadley circulation, reduced summer rainfall in the Mediterranean Basin when compared to the conditions of the Early Pliocene.

**Plain Language Summary** The modern Mediterranean Basin has a Mediterranean-Type Climate (MTC) where it receives most of its annual rainfall during the autumn and winter months, whereas the summers are warm and dry. However, environmental proxy records (such as pollen data) suggest that the region had a subtropical climate with no summer drought prior to ~3.3 Ma. An approach to understand regional climate changes in the basin is to look at soil carbonate $\delta^{18}\text{O}$ values, which largely record the $\delta^{18}\text{O}$ values of rainfall, through which we can then infer changes in the hydroclimate. This study examined $\delta^{18}\text{O}$ values of soil carbonates from Teruel, Spain from 4.1 to 2.4 Ma and compared them with other complementary paleoclimate and environmental proxy records retrieved from this region. The $\delta^{18}\text{O}_c$ values of our record show negative shifts during periods of known global cooling, which point to either more winter rainfall and/or less summer rainfall. Comparison with other records suggest a decrease of summer rainfall and the beginning of the MTC during periods of global cooling likely due to changes in oceanic and atmospheric systems.

## 1. Introduction

The Mediterranean Basin has a climate with hot summers and wet winters, dubbed a Mediterranean-Type Climate (MTC), that is typical of the western margins of continents at latitudes 30°–45° and often found adjacent to desert regions (Koppen, 1936; Lionello et al., 2006). Regions with these climates host high levels of biodiversity, with the continental Mediterranean, despite covering only ~2% of the global surface, containing ~20% of the existing vascular plant species (Médail & Quézel, 1997). The Iberian Peninsula (IP) is influenced by two different, and distinct, moisture sources: the North Atlantic Ocean and the Mediterranean Sea. Owing to the complex topography related to the Peninsula's major mountain chains and sharp orographic relief, rainfall across the region is both spatially and temporally variable (Figure 1) (Moreno et al., 2021; Serrano et al., 1999). Autumn, winter, and spring rainfall is largely sourced from the Atlantic Ocean through westerly winds and Atlantic fronts (Moreno et al., 2021; Serrano et al., 1999). This phenomenon can be further explained by a more equatorward and weaker subtropical high-pressure zone during the winter months compared to summer. This allows the IP to be affected

© 2025. The Author(s).
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Precipitation of the Iberian Peninsula as represented by the percentage of rainfall of the annual amount that occurs each season. Data set from Herrera et al. (2019).</figcaption></figure>

by westerly winds and baroclinic synoptic-scale perturbations sourced from the polar front (Cabos et al., 2020). Summers are typically dry, owing to a northward shifted Azores high-pressure system and a weakened Icelandic low-pressure system, allowing for the Azores high to deflect the westerly circulation and precipitation away from most of the peninsula, excluding the northern portion of the region (Cabos et al., 2020). In addition, the regional Hadley cell migrates northward during summer months, positioning its descending branch and associated subsidence largely over the central and eastern Mediterranean inducing drying (Rodwell & Hoskins, 1996). Summer rain is typically caused by local effects such as convective storms, backdoor cold fronts, high moisture content, and high atmosphere instability with the Mediterranean Sea as the dominant moisture source (Moreno et al., 2021; Serrano et al., 1999). On the northeastern IP, where the Teruel Basin is located, surrounding mountain ranges result in low winter rainfall with a majority of annual rainfall occurring during late spring through early autumn, associated with sea breezes that bring moisture from the Mediterranean Sea to the region during the warm season (Figure 1) (Moreno et al., 2021; Serrano et al., 1999). These complex dynamics present the opportunity for the climate regime to be influenced by both global and local changes.

Importantly, floral (Jiménez-Moreno et al., 2010), faunal (van Dam & Weltje, 1999), and palustrine carbonate records (Alonso-Zarza et al., 2012) from the IP suggest that, while rainfall seasonality existed during the Miocene —Early Pliocene, modern environmental conditions were not established. Pollen records older than 3.4 Ma are indicative of western Mediterranean terrestrial subtropical ecosystems requiring high humidity and mean annual temperatures between 20°C and 24°C and requiring high humidity, designated as subtropical taxa in this paper. Modern-like ecosystems established between 3.4 Ma and 3.25 Ma with the replacement of these subtropical-dominated taxa by high abundances of drought-tolerant plant species (Fauquette et al., 1999; Jiménez-Moreno et al., 2010; Suc, 1984; Suc et al., 2018). Concurrent with this climatic shift was the onset of Northern Hemisphere Glaciation (~3.3 Ma), followed by the intensification of the Northern Hemisphere Glaciation (iNHG: ~2.73 Ma) which peaked at 2.6 Ma (De Schepper et al., 2013; Hennissen et al., 2015; McClymont et al., 2023; Mudelsee & Raymo, 2005; Zachos et al., 2001). Because Mediterranean vegetation changed significantly throughout the Early-to-early Late Pliocene, it is important to consider the global climate evolution throughout this time period.

Subtropical regions during the Early-to-Mid Pliocene are considered to have been largely warmer and wetter than modern conditions. These background climate conditions are attributed to significantly reduced Polar ice sheets, warmer extratropical SSTs, which redistributed moisture north of 30°, and a weakened meridional SST gradient, generating a weaker Hadley circulation, and a weakened Intertropical Convergent Zone (ITCZ) where convection
decreased (Burls & Fedorov, 2017; Colleoni et al., 2015; Dowsett et al., 2011; Haywood & Valdes, 2004; Haywood et al., 2000; Salzmann et al., 2011). Prior to the onset of Northern Hemisphere Glaciation there were several failed initiations of the permanent glaciers in the northern hemisphere occurring at 4.0 Ma (MIS Gi22/Gi20), 3.6 Ma (MIS Gi2 and/or Gi4), and 3.3 Ma (MIS M2) (De Schepper et al., 2014). Ultimately, however, glacial growth did not occur in earnest until the onset of the Northern Hemisphere Glaciation (~3.3 Ma), featuring low-amplitude glacial-interglacial cycles, with amplification of ice growth occurring during the iNHG (Lisiecki & Raymo, 2005; Loubere & Moss, 1986; Raymo, 1994; Rohling et al., 2021).

Thus, to better understand the response of regional climates of the IP to past global climate changes, we measured the δ¹⁸O values of pedogenic carbonates collected from the Teruel Basin, northeastern Spain, spanning 4.1–2.4 Ma. We compared our record with existing proxy records from the Mediterranean Basin to advance our understanding of the climatic changes influencing the basin during the Late Pliocene to the beginning of the Pleistocene responsible for the shift to the modern Mediterranean-type environment. Ultimately, detailed records of terrestrial hydroclimate spanning the Late Pliocene provide opportunities to assess mechanistic linkages between these global-scale climate changes and regional environmental turnover in the Mediterranean Basin.

## 2. Background

### 2.1. Geographic/Geologic Setting

The study sites (VA 2016 and VA 2019) are located in the 100 km-long Teruel Basin, Spain which crosses the southeast compressional Alpine Iberian Chain (Figure 2). The basin is considered a half-graben bounded by en echelon faults trending NNE-SSW largely within the eastern portion of the basin (Alonso-Zarza & Calvo, 2000; Anadón & Moissenet, 1996). This half-graben was formed in relation to rifting in the western Mediterranean, associated with the Valencia Trough, during the Miocene that led to crustal extension (Alonso-Zarza & Calvo, 2000; Anadón et al., 1989; Ezquerro et al., 2019; Roca, 1994; Simón et al., 2019).

The basin is infilled with ~500 m thick Neogene sequences deposited by variably sized alluvial fans at its margins that transitioned to palustrine and lacustrine conditions in the central region (Alonso-Zarza & Calvo, 2000; Ezquerro, 2017; Ezquerro et al., 2012, 2014, 2016; Moissenet, 1989; Simón et al., 2019). The northern Teruel Basin can be stratigraphically divided into six megasequences limited by alluvial progradations (Ezquerro et al., 2012, 2014). Within these megasequences are interpreted periods of lake expansion and contraction of a closed shallow lake within the Teruel Basin throughout the Miocene-Pliocene that is typic of endorheic active basins (Ezquerro, 2017; Ezquerro et al., 2014). Extensional tectonics and climate influence the sedimentation of the basin with active (decreased) tectonism associated with lake retraction (expansion) and alluvial fan progradation (retrogradation) (Ezquerro, 2017; Ezquerro et al., 2016, 2019, 2022; Liesa et al., 2019).

The Villalba Alta section sampled in this study lies along the western border of the northern end of the Teruel Basin. The sediments of this 95 m section are characterized by two main depositional environments: alluvial (mudstones and sandstones rich in soil carbonates) and lacustrine (limestones) (Gao et al., 2022). The paleosols sampled in this study contain nodular carbonates ranging from a few mm to 2 cm in diameter embedded in a clay or sandy-clay matrix.

Laminated light to dark gray limestones from 0 to 12 m and limestone units transitioning to black mudstone from 24 to 35 m are indicative of a water-saturated depositional setting, which we interpret as a lacustrine environment (Figure 3). The mudstones in this part of the section may have been deposited and pedogenically modified in mudflat deposits proximal to shallow lakes. Within this lower part of the section (below 17 m), discrete paleosols are exposed. The most intact discrete paleosol (Figure S2 in Supporting Information S1) consists of a 10–20 cm thick, drab greenish-gray B horizon with subrounded, blocky structure which transition abruptly to a yellowish, nodular Bk horizon. The top 10 cm of the Bk horizon contains the highest density of nodules (25%). The Bk horizon transitions gradually to a red-brown Ck and then C horizon with prismatic structure. The entire soil is about 70 cm thick. Taken together, the absence of primary sedimentary structures, the distinct horization defined by both color and structure and the occurrence of nodules strongly indicate this is a fossil soil. The other paleosols sampled in the lower part of the section contain a similar Bk-Ck-C horization but do not contain the overlying greenish-gray B horizon.

The section is dominated by red and calcareous mudstones from 12 to 24 m. Red mudstones, sandy mudstones, sandstones, and conglomerates characterize the section from 35 to 97 m (Figure 3). Paleosols above 17 m tend to
<figure><img src="image_3.png" /><figcaption>Figure 2.</figcaption></figure>
consist of Bk horizons welded on top of one another. For instance, the mudstone from 80 to 85 m consists of a continuous nodule-bearing Bk horizon, probably the result of the rate of soil formation keeping pace with sedimentation. Above 35 m, the abundance of laterally discontinuous, cross-bedded conglomerates, and sandstones suggests intervening mudstones in which the soils developed are overbank deposits within an alluvial system. The change of facies throughout the section suggests shifts from shallow lacustrine to alluvial depositional settings, consistent with other interpretations in the basin (Alonso-Zarza & Arenas, 2004; Ezquerro et al., 2019; Gao et al., 2022).

## 2.2. Pedogenic Carbonates and $\delta^{18}O$

Pedogenic carbonates record $\delta^{18}O$ values of soil water during their formation (Cerling & Quade, 1993). Rainfall-$\delta^{18}O$ is the primary control of soil carbonate $\delta^{18}O$ values (Bellanca et al., 1992; Breecker et al., 2009; Cerling, 1984) whereas evaporation is typically of secondary, although geographically widespread, importance (Kelson et al., 2023). As winter and summer precipitation have distinct $\delta^{18}O$ values, seasonal changes in rainfall affect the $\delta^{18}O$ values of the soil carbonates (Da et al., 2023). Modern rainfall in the Teruel Basin exhibits distinct isotope values for winter and summer months, with a weighted average of -5.4‰ during the summer and -10.8‰ in winter months (IAEA/WMO (2024) Global Network of Isotopes in Precipitation. The GNIP Database). Pedogenic carbonates thus can record seasonal rainfall changes in the Mediterranean Basin during the Late Pliocene. It is of note that diagenesis and/or overprinting of the original nodular material by recrystallization can alter the original isotopic signal, which may produce a record not representative of the climate the nodule formed in (Cerling, 1984; Deutz et al., 2001). This possibility of post-depositional alteration requires care in sampling unaltered regions of paleosol carbonate nodules.

# 3. Materials and Methods

## 3.1. Stable Isotope Analysis

Expeditions to the Teruel Basin were made in 2016 and again in 2019 to sample sections about ~1 km apart. The section sampled as part of a pilot study in 2016 (also studied by Opdyke et al. (1997)) is located just north of the road that traverses the escarpment to the west of the town of Villalba Alta (Figure 2b). The section sampled in 2019 (also sampled by Ezquerro et al. (2019)) is ~90 m tall and located just south of the same road; the top of this section is adjacent to the hairpin turn (Figure 2b) and was chosen for sampling in 2019 due to its better exposure. The section sampled in 2019 was measured as height above the lowest outcrop whereas the section sampled in 2016 was measured as depth below the prominent cliff-forming conglomerate at the top of the section. The age models were developed by Gao et al. (2022) using paleomagnetic data gathered during the 2019 sampling campaign, and the magnetostratigraphy developed by Opdyke et al. (1997) for the 2016 data set (Figure S1 in Supporting Information S1).

During the 2016 pilot expedition, 52 nodules (2–6 nodules per stratigraphic level) were collected from 13 different depths. During the 2019 expedition, 3–14 nodules were collected from each stratigraphic level at which nodules were observed resulting in 1,257 discrete soil carbonate nodules from 148 different heights (Figure 3).

In total, 342 nodules from 56 different depths were selected for stable isotope analyses. In the laboratory, nodules were split in half using a mallet and then the micrite was sampled with a Dremel drill avoiding spar. The resulting powders were then reacted with $H_3PO_4$ at 50°C to liberate $CO_2$ before being analyzed on a Thermo 253 isotope ratio mass spectrometer operating in continuous flow mode and interfaced with a Thermo Gasbench II (Spötl & Vennemann, 2003) at the University of Texas at Austin. Drift corrections were applied to the soil carbonate $\delta$ values using an internal marble (UTM) standard while laboratory standards NBS-18, NBS-19 were used to

**Figure 2.** Geologic map of the Mediterranean Basin and Teruel Basin. (a) Approximate proxy sampling locations used in this study are indicated by the following: blue stars represent sampling locations for North African (ODP 659, Amarathunga et al., 2024; Grant et al., 2022/ODP 658, Leroy & Dupont, 1997/DSDP397, Stein, 1985) and Eastern Mediterranean records (ODP 967, Amarathunga et al., 2024; Grant et al., 2022; Heslop et al., 2023), black star represents Mediterranean sea surface temperature (SST) reconstruction sampling locations from the Capo Rossello composite section (Beltran et al., 2011, 2021; Herbert et al., 2015, 2018), green star represents Garraf 1 pollen data set location (Suc, 1984; Suc et al., 2018), North Atlantic SST reconstruction sampling location is represented by an orange star (Naafs et al., 2020), and the Teruel Basin of this study is denoted by the square. The approximate seasonally averaged location of the modern position of the core of the Midlatitude Westerly Jet is shown as a dashed line (Belmecheri et al., 2017). (b) Location of study area in the basin (Villalba Alta Coordinates 0 58' W, 40 37' N). Red and purple stars indicate sampling locations of the two sampling campaigns (2019 and 2016). Map units adapted from Ezquerro et al. (2019).
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Stratigraphy of the Villalba Alta section sampled in 2019. Due to figure size constraints, the 10 cm interval measurement of section, and close sampling collection of carbonate nodules, the stratigraphic height for some nodules is approximated.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** (a) Image of one sampled nodule used in spar versus micrite analyses. The height of the nodule in the 2019 VA stratigraphic column is 84.60 m. (b) Magnified image of targeted sparry area. Sparry materials are outlined in white.</figcaption></figure>

calibrate delta values to the PDB scale. Averaging the standard deviations of the NBS-18, NBS-19, and UTM across all runs produced a reproducibility of 0.16‰.

Additionally, spar and coexisting micrite from eight nodules were imaged (Figure 4), sampled separately, and analyzed using the same procedure to determine if they had $\delta^{18}\text{O}$ values distinct from each other.

## 3.2. Cathodoluminescence

Fourteen nodules were selected for CL imaging. Thin sections (30 microns thick) with a microprobe polish were created for each nodule by National Petrographic. Imaging was done at the Kansas Geological Survey Digital Cathodoluminescence Imaging Laboratory using a Relion Industries Reliotron III cold cathode cathodoluminescence chamber under both transmitted and reflected light. Macroscopic imaging, using a 16 Mpx Canon EOS SL1 DSLR camera, was completed on all 14 samples. Two samples of interest were selected for microscopic imaging using a 17 Mpx Olympus DP73 color camera on an Olympus BX41 compound microscope.

# 4. Results

## 4.1. $\delta^{18}\text{O}$ Data

$\delta^{18}\text{O}_c$ values are plotted versus depth for the 2016 sampling campaign (Figure 5a) and versus height for the 2019 data set (Figure 5b). When $\delta^{18}\text{O}$ values are plotted versus age, the data from the two sections (2016 and 2019) align well with each other (Figure 5c), and data from both sections are therefore considered together in this study. The composite $\delta^{18}\text{O}_c$ record can be divided into three time intervals based on its variations. $\delta^{18}\text{O}_c$ values fluctuate between -7.0 and -4.8‰ (mean = -5.8‰, SD = 0.4‰) during 4.1–3.35 Ma with no long-term trend (Figure 5c). This is followed by a brief negative excursion, most prominently seen in the 2019 data set, defined by nodules from two stratigraphic levels at ~3.28 Ma, reaching the lowest value of <-7‰. Finally, $\delta^{18}\text{O}_c$ values generally decrease from -4.88‰ at 3.2 Ma to ~-6.0‰ at 2.7 Ma.

## 4.2. Spar and Micrite $\delta^{18}\text{O}$ and $\delta^{13}\text{C}$ Data

Seven nodules were documented and analyzed for their spar and micrite $\delta^{18}\text{O}$ values (Figure 6).

Absolute values of differences between spar and micrite $\delta^{18}\text{O}$ and $\delta^{13}\text{C}$ values averaged at 0.4 and 0.3‰, respectively, suggesting there was unlikely a large influence of inadvertently sampled diagenetic spar on the $\delta^{18}\text{O}$ record reported here. To better estimate the magnitude of effect on $\delta^{18}\text{O}$ values from spar contamination of the drilled micrite powders, an isotope mass balance was considered:
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Paleosol carbonate δ¹⁸O values from Villalba Alta. (a) Average δ¹⁸O values versus depth for the section sampled in 2016. Each circle represents the mean of 1–3 nodules. (b) Average δ¹⁸O values versus height for the section sampled in 2019. Each circle represents the mean of 3–14 nodules. (c) Combined average δ¹⁸O values versus age for both the 2016 and 2019 sampling campaigns. Gray arrows highlight negative shifts in the δ¹⁸O data and the gray circle outlines the scattered data with no long-term trend of the oldest part of the record. Age values for the 2016 data set are resolved from Opdyke et al. (1997) and the 2019 data set is constrained by Gao et al. (2022). In all panels, error bars represent standard errors.</figcaption></figure>

$$ \delta^{18}\text{O}_{\text{total}} = (\delta^{18}\text{O}_{\text{micrite}} * F_{\text{micrite}}) + (\delta^{18}\text{O}_{\text{spar}} * F_{\text{spar}}) $$

where *F* is the fraction and the subscript *total* represents the powder generated from drilling micrite with possible contamination from some spar, the subscript *micrite* represents the actual δ¹⁸O value of micrite free from any spar contamination and the subscript *spar* represents the targeted sparry material. Spar content across the seven studied nodules were estimated between 5% and 40% (average 20%) of total area of the surface that was drilled into. We make a conservative calculation that assumes spar was not avoided, when in actuality we did our best to avoid the spar when targeting micrite. We solved the equation above for δ¹⁸O micrite:
<figure><img src="image_7.png" /><figcaption>Figure 6. Paleosol carbonate δ¹⁸O and δ¹³C values from Villalba Alta. Circles represent targeted sparry material while triangles represent the micrite drilled in an attempt to avoid all spar. Arrows point from micrite to spar values for each sample set. Sample depths are reported in m.</figcaption></figure>

$$ \delta^{18}\text{O}_{\text{micrite}} = \delta^{18}\text{O}_{\text{total}} - \left( \frac{\delta^{18}\text{O}_{\text{spar}} * (1 - F_{\text{micrite}})}{F_{\text{micrite}}} \right) $$

and propagated error with a Monte Carlo simulation. We assigned a uniform distribution to the uncertainties of all variables on the right-hand side, using the maximum and minimum values determined from our subset of nodules in which we targeted spar (n = 7). For Fmicrite we assigned a maximum value of 1 and a minimum of 0.6. Fspar was then calculated during each iteration as:

$$ F_{\text{spar}} = 1 - F_{\text{micrite}} $$

The distribution for δ¹⁸O<sub>spar</sub> was bounded between a minimum of -6.7‰ and a maximum value of -5.3‰. To ensure a paired spar and micrite value from each nodule is considered and avoiding comparing spar and micrite δ¹⁸O values between different nodules (i.e., considering the possible dependence between spar and micrite δ¹⁸O values), the δ¹⁸O<sub>total</sub> variable is defined as:

$$ \delta^{18}\text{O}_{\text{total}} = \delta^{18}\text{O}_{\text{spar}} + \Delta^{18}\text{O}_{\text{(total-spar)}} $$

The Δ¹⁸O<sub>(total-spar)</sub> values, calculated by subtracting the spar δ¹⁸O value from the micrite δ¹⁸O for each of the seven nodules, were bounded by a maximum difference of 0.7‰ and a minimum of -0.9‰. We performed 10,000 iterations with random values chosen from each distribution during each iteration. The mean difference between the calculated (pure) δ¹⁸O<sub>micrite</sub> and the measured δ¹⁸O<sub>total</sub> was 0.1 ± 0.2‰, which represents the error introduced by potential spar contamination. This error is only slightly larger than our analytical precision and similar to the mean of standard error among replicate (same stratigraphic level) nodules. Therefore, we are confident that the reported δ¹⁸O<sub>c</sub> reflects the original isotopic value at precipitation.

### 4.3. Cathodoluminescence Imaging

Cathodoluminescence (CL) is a useful tool to investigate not only the conditions soil carbonates formed under, particularly whether they were precipitated in reducing or oxidizing conditions, but additionally to better
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Selected CL images of nodules exhibiting the following luminescence types: concurrent orange and violet luminescence ((a): 25.20 m (3.58 Ma)), violet luminescence ((b): 85.15 m (2.75 Ma)), no luminescence ((c): 66.90 m (3.12 Ma)), and bright orange luminescence ((d): 15.30 m (3.91 Ma), (e): 16.30 m (3.87 Ma), and (f): 22.10 m (3.66 Ma)).</figcaption></figure>

determine the existence of diagenetic overprinting of original material (Aguirre Palafox et al., 2024; Barnaby & Rimstidt, 1989; Hemming et al., 1989; Richter et al., 2003). Reducing conditions during initial soil formation and/ or during diagenesis can produce bright orange and/or weak purple/blue luminescence depending on the ratio of Mn²⁺/Fe²⁺. Eight nodules exhibited certain luminescence, ranging from bright orange to weak purple/blue observed under CL (Figure 7). Generally, nodules from the highest part of section (43.50–85.15 m) exhibited little to no luminescence while nodules at the bottom of section (15.30–22.10 m) exhibited bright orange luminescence (Figure 7). The nodules with bright orange luminescence occur toward the bottom of the section at 15.30, 16.30, and 22.10 m (3.91, 3.87, and 3.66 Ma). A mixture of blue, orange, and violet luminescence was seen in one nodule from 25.50 m (3.58 Ma). Very faint purple/blue hues were noted in three nodules from 53.65, 78.10, and 85.15 m (3.26, 2.87, and 2.75 Ma), and were the weakest examples of luminescence. Six samples at 43.50, 51.30, 66.90, 77.70, 81.70, and 84.60 m (3.37, 3.28, 3.12, 2.88, 2.81, and 2.76 Ma) were found to have no luminescence (Figure 7).

# 5. Discussion

## 5.1. Regional Climactic Conditions of the Teruel Basin in Pliocene

To interpret our δ¹⁸O_c record we first contextualize the data set by examining the local setting in which the studied soil carbonates formed, beginning with the depositional environment of the Villalba Alta section during the Late Pliocene—Early Pleistocene. The lower portions of the section (~4.3 to ~3.5 Ma) occur during megasequences M-IIbc and M-III, representing a transition from alluvial progradation to alluvial retrogradation sourced from alluvial fans in the eastern margin of the basin (Ezquerro et al., 2019). The upper portion of the section (~3.5 to ~2.5 Ma) aligns with megasequence M-IV which featured alluvial progradation throughout (Ezquerro et al., 2019). Ultimately, the transition from largely lacustrine sediments from 0 to 12 m (~4.3 to ~4.0 Ma) and 24 to 35 m (~3.6 to ~3.5 Ma) to alluvial fan sediments from 12 to 24 m (~4.0 to ~3.6 Ma) and 35 to 97 m (~3.5 to ~2.5 Ma) (Figure 8i) indicates a shift toward more alluvial conditions, reflecting progressive drying and/or increased tectonic activity.
Cathodoluminescence (CL) imaging reveals orange luminescence of the nodules in the lower part of the 2019 Villalba Alta section until ~25 m (~3.6 Ma) (Figure 8j). The main controls on the CL intensity of calcite are manganese (Mn²⁺), responsible for the orange luminescence, and iron (Fe²⁺), which quenches luminescence (Hemming et al., 1989). The respective concentrations of manganese and iron, along with the oxidation state of either element, are important for the brightness, or lack thereof, of luminescence (Barnaby & Rimstidt, 1989). Soil carbonate with higher Mn²⁺ and/or Fe²⁺ concentrations is typically formed under water-saturated, reducing environments, which prompts the dissolution of iron- and manganese-bearing minerals. The release of these elements into the soil solution and their incorporation into growing calcites, leading to bright orange luminescence in the case of Mn²⁺. Previous research (Hemming et al., 1989; Richter et al., 2003) suggests that the presence of Mn²⁺ is a dominant feature in luminescence, with calcite emitting mild luminescence even with small amounts of Mn²⁺ and moderate amounts of Fe²⁺.

Pedogenic carbonate nodules may additionally be subjected to water saturated environments post-formation as they are buried, positioning the nodules in the phreatic zone below the water table and spurring orange luminescence through phreatic overprinting (Aguirre Palafox et al., 2024). However, the luminescent nodules of older (lower) portion of this study exhibit much less (and in some cases no) sparry infill when compared to the younger (upper) portion. Because diagenetic overprinting/infilling typically forms spar (Budd et al., 2002; Cerling, 1984; Wieder & Yaalon, 1982), we favor pedogenic reducing conditions at the time of formation. Thus, the orange luminescence in the oldest part of our record (14–25 m, ~4.1–~3.6 Ma) points to a water-saturated and reducing environment, while weak and/or no luminescence of soil nodules representative of the upper part of the section suggests a shift to formation in a better drained and thus more aerated soil environment (Layzell et al., 2023; Ross et al., 2017), consistent with the intercalation with coarser grained sediment in a fluvial depositional setting, related to a drier depositional environment and/or climate (Ezquerro et al., 2019, 2022; Liesa et al., 2019).

While climate and tectonic signals can be difficult to uncouple, several studies in the Teruel Basin have sought to determine regional climate signals across the Miocene—Pliocene (Alonso-Zarza et al., 2012; Ezquerro et al., 2014; Lafuente et al., 2010; van Dam et al., 2023). The IP in the Early-early Late Pliocene, prior to significant Northern Hemisphere glaciation at ~3.3 Ma, was characterized by generally warmer and wetter conditions with a lower degree of seasonality (Ezquerro et al., 2022; Fauquette et al., 1999; Gao et al., 2022; Jiménez-Moreno et al., 2010; Suc et al., 2018; van Dam et al., 2006). Aeolian dune fields (Rodríguez-López et al., 2012), stable isotope records from fossil horse teeth (van Dam & Reichart, 2009), and lacustrine carbonate data (Ezquerro et al., 2022) point toward increasingly arid conditions on the IP following the onset of Pliocene global cooling and peaking from ~2.9 to 2.5 Ma. Additionally, superimposed on this drying trend of the Teruel Basin were orbitally forced cyclical wet/dry conditions in the region (Colleoni et al., 2012; Ezquerro et al., 2022; Jiménez-Moreno et al., 2010; Ochoa et al., 2018; Rodríguez-López et al., 2012; van Dam & Reichart, 2009). Thus, while periods of active tectonism that occurred during this time period likely influenced the sedimentation in the Teruel Basin (Ezquerro et al., 2014, 2019, 2022; Liesa et al., 2019), climate proxy records suggest an overall trend toward increasing aridity from ~3.3 to ~2.5 Ma, punctuated by shorter-term orbitally forced wet/dry periods.

## 5.2. Pedogenic Carbonate δ¹⁸O Record and Drying in Mediterranean Basin

A change in rainfall seasonality and the establishment of a summer drought is the most likely explanation for the observed negative δ¹⁸O_c shifts, given that modern summer rainfall-δ¹⁸O values are on average ~5.4‰ higher than those of winter rainfall across the peninsula (IAEA/WMO (2024) Global Network of Isotopes in Precipitation. The GNIP Database). The shift from a subtropical to winter wet/summer dry Mediterranean climate during 3.3 to 2.7 Ma is consistent with the multiple proxy records examined in this study. Here we compare several important proxy records from the Mediterranean Basin to our δ¹⁸O_c to further support this contention (Figure 8).

First, the δ¹⁸O_L values of lacustrine carbonates in the Villalba Alta region from the Teruel Basin increased when nearby paleosol-carbonate-δ¹⁸O_c values decreased and almost reach a minimum (Figure 8e). The δ¹⁸O values of lacustrine carbonates are controlled by temperature and the isotopic composition of the lake, which for closed, shallow lakes, is largely controlled by evaporation (Horton et al., 2016; Leng & Marshall, 2004). Sedimentary and geochemical data from Ezquerro et al. (2022) suggests a closed lake in the Teruel Basin from ~3.6 to ~1.8 Ma. Thus, the higher δ¹⁸O_L values of lacustrine carbonates (Figure 8e) during the Late Pliocene are likely due to enhanced evaporation in the region as a result of intensified summer drought. Compared to closed-basin
<figure><img src="image_9.png" /><figcaption>Figure 8.</figcaption></figure>
lake waters (Horton et al., 2016; Leng & Marshall, 2004), the δ¹⁸O values of soil waters in the Bk horizon are not as sensitive to changes in evaporation in humid and semi-arid climates (Kelson et al., 2023) that characterize the Teruel Basin. Taking these considerations into account, a shift toward summer drought explains the relatively small δ¹⁸O_c decreases (more influenced by change in rainfall δ¹⁸O) of paleosol carbonates and the relatively large δ¹⁸O increase of lake carbonates (more influenced by evaporation) during the iNHG. The concurrent increase that aligns with the negative shift in δ¹⁸O_c values at ~3.1 Ma, and variability of δ¹⁸O_L, supports both an increase in evaporation in the region and influence of seasonal oscillations (Ezquerro et al., 2022).

Second, frequency-dependent magnetic susceptibility ($\chi_{fd}$) shows an increase beginning at 3.4 Ma pointing toward the onset of seasonal wet-dry cycles and thus to aridification of the study area through the establishment of a drought season (Figure 8g). Studies from the Chinese Loess Plateau suggest that magnetic enhancement (higher $\chi_{fd}$) results from the accumulation of nanoscale ferromagnetic materials during pedogenesis, which form under oxidation-reduction cycles (Ahmed & Maher, 2018). In the Teruel Basin, the seasonal cycles of wet winters and dry summers may result in these redox cycles with winter/spring rainfall filling pore space and forming reducing conditions and summer drying creating oxidizing conditions (Gao et al., 2022).

Third, the Garraf 1 pollen data set from offshore Barcelona (Suc, 1984) (Figure 2a) demonstrates a shift from subtropical taxa to drought-tolerant taxa (herbs, steppe plants, Mediterranean sclerophyllous plants) that began ~3.4 Ma (Figure 8h). This shift intensified after ~3.2 Ma, following the M2 event, and continued to ~2.6 Ma. The trend reversed for a short time, with an increase in the relative abundance of subtropical taxa from ~3.1 to ~2.8 Ma after which the original trend continued (Figure 8h). It is also of note that after the MIS M2 event, the δ¹⁸O_c values track these shifts in pollen assemblages fairly closely. The establishment of a summer drought, represented by decreases in δ¹⁸O_c values, would explain the observed shift from subtropical to drought-tolerant plants.

Fourth, the salinity of Mediterranean bottom water, a proxy for Mediterranean Outflow Water (MOW) intensity, showcases the increase in aridity in the Eastern Mediterranean Basin and North Africa over the Pliocene. Salinity increases in bottom waters are associated with reduced freshwater input into the Mediterranean Sea and thus increased aridity. The δ¹⁸O<sub>w</sub> values, which represent estimated δ¹⁸O values of bottom waters and thus salinity, reported by Khélifi et al. (2014) increased from ~3.6 Ma onwards, indicating reduced freshwater input and suggesting drying began earlier than other records indicate (Figure 8f). Additionally, the onset of more frequent, thinner sapropel formation after 3.2 Ma linked to a drier, less vegetated, and more erodible North African landscape with increased seasonality and interrupted humid periods (Amarathunga et al., 2024) further supports a trend toward overall drying from North Africa to the Mediterranean Basin.

### 5.3. Additional Factors That Control Pedogenic Carbonate δ¹⁸O Values

Comparing our record to the LR04 benthic δ¹⁸O (δ¹⁸O_b), we see that from 3.3 to 2.5 Ma, our record generally mirrored the long-term trends in the benthic stack (Figures 8a and 8d). Benthic foraminifera δ¹⁸O values record changes in ice volume and deep-sea temperature and are important indicators of global climate conditions, generally biased toward high latitudes (Rohling et al., 2021). Thus, the similarity between these two curves indicates a sensitivity of the paleosol-carbonate-δ¹⁸O values to global cooling.

δ¹⁸O_c decreased from 3.1 to 2.5 Ma during the onset of Northern Hemisphere Glaciation and iNHG, and although our record generally does not resolve changes as the glacial/interglacial timescales, two especially low δ¹⁸O_c values occurred during the MIS M2 event (Figures 8a and 8d). The lack of resolution on shorter time scales in our age model is a result of coarse sampling resolution of paleomagnetic cores and linear interpolation of six age tie points related to the three reversed subchrons captured, which in turn have their own associated uncertainties of

**Figure 8.** Comparison of mid-late Pliocene proxy records. (a) δ¹⁸O_b values of the LR04 benthic foraminifera stack (Lisiecki & Raymo, 2005). (b) North Atlantic sea surface temperature (SST) determined from $U^{K}_{37}'$ values (Naafs et al., 2020). (c) Mediterranean SST determined from composite $U^{K}_{37}'$ values (Beltran et al., 2021). (d) Average δ¹⁸O_c values of pedogenic carbonates measured in this study (3–14 nodules per stratigraphic level, mean and ± 1 s.e. shown). (e) δ¹⁸O_L values of lake carbonates from the Teruel Basin (Ezquerro et al., 2022). Data set is a composite of two sites: SB-5 (900 m north of San Blas neighborhood) and Or-3 (600 m northwest of Orrios) (Ezquerro, 2017). (f) Bottom water δ¹⁸O_w values of Mediterranean site 978 (Khélifi et al., 2014). (g) $\chi_{fd}$ records from the Villalba Alta section (Gao et al., 2022). (h) Garraf 1 pollen record from offshore Barcelona depicting the percentages of drought-tolerant taxa (including Mediterranean sclerophyllous taxa) vs subtropical taxa (Suc et al., 2018). (i) Stratigraphic column of the Villalba Alta section, scaled to the time domain to facilitate comparison of lithological changes with other records. (j) Cathodoluminescence images of selected soil carbonate nodules from the 2019 Villalba Alta section in the Teruel Basin.
±~5–10 thousand years (Gao et al., 2022; Gradstein et al., 2020). We therefore are cautious in interpreting the following trends: a negative excursion in δ¹⁸O_c values defined by data from two stratigraphic levels at ~2.9 Ma, which may correspond to the short glacial period MIS G14 (~2.89 to ~2.88 Ma) and a positive excursion at ~2.75 Ma that may correspond to the extending interglacial MIS G7 (Lisiecki & Raymo, 2005) (Figure 8). Adding confidence to our recording of MIS G7, we note that North Atlantic SSTs remained elevated for nearly 0.5 myr near the highest temperatures recorded after 3.1 Ma (Naafs et al., 2020; Figure 8b), meanwhile the long-term expansion of drought tolerant taxa and paleosol carbonate-inferred summer drying trends both reversed temporarily (Figure 8h). Thus, while we largely interpret this δ¹⁸O_c data set as a record of long-term climate averages, it is possible this record is sensitive to individual interglacial/glacial stages.

Although they align temporally, the shifts we observed are in the opposite direction expected from ice-volume-controlled changes in seawater δ¹⁸O values. Changes in δ¹⁸O values of sea surface water (moisture source regions) during this cooling probably also do not explain the paleosol carbonate δ¹⁸O variability. Benthic oxygen isotope records from ODP Site 659 in the North Atlantic across the Pliocene-Pleistocene record a gradual increase of 1‰ from 3.15 to 2.5 Ma, which is again in the opposite direction of our observed shift in our soil carbonate δ¹⁸O record (Tiedemann et al., 1994). This is also true for the Mediterranean Sea, where data from ODP Site 967 derived from planktonic foraminifera show a similar trend in increasing δ¹⁸O values during the beginning of the Late Pliocene (Colleoni et al., 2012; Lourens, 2004; Wang et al., 2010). Therefore, secular changes in seawater δ¹⁸O values cannot explain our observations.

Temperature affects δ¹⁸O_c values by influencing the δ¹⁸O values of rainfall as well as the isotope fractionation between calcite and water (Cerling & Quade, 1993; Coplen, 2007; Dansgaard, 1964; Kim & O'Neil, 1997). If we consider only the effect of temperature on fractionation of oxygen between calcite and water, where the isotopic fractionation between soil water and calcium carbonate decreases with increasing temperatures, global cooling would result in increased fractionation and higher calcite δ¹⁸O values. However, as temperature affects the δ¹⁸O value of rainfall in an opposite relationship where cooling leads to more negative δ¹⁸O values, the relationship between temperature and δ¹⁸O values of soil carbonates is complex. To investigate the effect of this we look at the terrestrial paleotemperature equation derived from the δ¹⁸O of pedogenic carbonates (Dworkin et al., 2005), which considers the effect of temperature on both rainfall δ¹⁸O values (based on the spatial and temporal relationship between δ¹⁸Oprecipitation and Mean Annual Temperature, MAT, across modern midlatitudes) and the magnitude of the calcite-water oxygen isotope fractionation factor. Dworkin et al. (2005) favors the simultaneous solution based on modern spatial distributions of δ¹⁸Oprecipitation, so we use that equation here:

$$-0.50T^3 + (\delta^{18}O_{cc(SMOW)} + 152.04)T^2 - 2.78 * 10^6 = 0$$

We calculated that ~2° of surface cooling, relative to modern conditions, most representative of the MIS M2 and iNHG (Iberian MAT, estimated from models by Haywood & Valdes (2004) and Dolan et al. (2015)) can explain only 0.5‰ decrease in soil carbonate δ¹⁸O values, which is much smaller than our observations (>1.5‰), indicating other processes at play. Therefore, cooling alone can explain only about 1/3 of the signal we observe, indicating some other process must account for the remainder unless Iberian cooling was much larger than currently thought. Additionally, the modern midlatitude locations (25°–60° regions) used to construct the MAT/ δ¹⁸Oprecipitation relationship governing Dworkin's equation largely have climates in which rainfall occurs uniformly throughout the year. As the IP is hypothesized to have transitioned from a climate with a weaker rainfall seasonality component during the Early Pliocene to its modern MTC, where a significant drought season exists, the failure of this equation to account for the large shift in the δ¹⁸O_c record suggests that a change in rainfall seasonality likely occurred.

## 5.4. Mechanisms for Summer Drying

There are several mechanisms that could explain summer drying in the Mediterranean Basin during the iNHG, with the favored being related to cooling of the North Atlantic. At ~3.3 Ma and again from ~3.1 to 2.4 Ma, North Atlantic SSTs decrease (Naafs et al., 2020; Figure 8b) while Mediterranean SSTs did not begin to decrease until later at ~2.6 Ma (Beltran et al., 2021; Figure 8c). Cooling of the North Atlantic has a two-fold effect on Iberian rainfall: decreased moisture availability and an enhanced meridional temperature gradient that generated a
southward displacement and enhancement of the Hadley circulation, ITCZ, and mid-latitude jet (Brierley et al., 2009; Burls & Fedorov, 2017).

First, decreased SSTs in both bodies of water would act to decrease available moisture for advection over the Mediterranean Region. Less available moisture from both precipitation sources that contribute to summer and winter rainfall would result in less rainfall to the region annually. This is supported by climate models across the Pliocene-Pleistocene transition that suggest a cooling-associated decrease in specific humidity north of 20° that contributed to a simulated decrease of rainfall in both winter and summer months in the IP (Colleoni et al., 2015).

Second, the effect of an enhanced meridional temperature gradient can be considered when comparing background Early Pliocene conditions with the Pliocene-Pleistocene transition. Climate simulations and SST reconstructions show a wet Mediterranean was developed during the warm Early Pliocene due to warmer extratropical SSTs and a weakened meridional SST gradient that generated a weaker Hadley circulation (Burls & Fedorov, 2017; Colleoni et al., 2015). As a result ITCZ convection decreased and, coupled with warmer extra-tropical SSTs, caused a redistribution of moisture north of 30°N, creating a wetter Mediterranean (Colleoni et al., 2015). Following these background wet Mediterranean conditions, global cooling created an enhanced meridional temperature gradient through lower North Atlantic SSTs that may have shifted the ITCZ southward toward the warmer hemisphere (Broccoli et al., 2006; Chiang et al., 2003). Furthermore, a contracted Hadley circulation and southward shifted ITCZ positioned the subsiding branch of the Hadley circulation southward near the Mediterranean Basin and northern Africa during the summer, bringing dry conditions to these regions (Colleoni et al., 2015; Haywood et al., 2000; Sun et al., 2013). In addition, the Hadley cell intensified and shrank, resulting in decreased moisture transport from the equator toward higher latitudes (Amarathunga et al., 2024; Colleoni et al., 2015).

The mechanistic linkages above are supported by several proxy and model based studies connecting Eastern Mediterranean and North African climate changes to the onset of Northern Hemisphere glaciation at ~3.3–3.2 Ma (Amarathunga et al., 2024; Grant et al., 2022; Heslop et al., 2023; Leroy & Dupont, 1997; Popescu et al., 2010; Stein, 1985). Pollen evidence from the Eastern Mediterranean (DSDP site 380) highlights an abrupt shift from mega-mesothermic plants to more drought-tolerant taxa at ~3.3 Ma, pointing to cooler and drier conditions (Biltekin, 2010; Popescu et al., 2010). As this region has an MTC and receives much of its winter rainfall from westerlies to the north and northwest, this clear shift in the pollen record to drier conditions suggest a concurrent shift in the westerlies (Biltekin, 2010; Popescu et al., 2010). Enhancement of tropical circulation is linked to clear changes seen in several climatic proxies from these regions. Pollen records off of the coast of northwest Africa (ODP Site 658) yield high percentages of drought-tolerant plant groups in the Late Pliocene, with a maximum arid phase reached at 3.26 Ma, linked to an enhancement of trade winds (Leroy & Dupont, 1997). Coarser grain size and increased accumulation of terrigenous matter followed by increased paleoproductivity, in a similar site off this coast (DSDP Site 397) at 3.2 Ma additionally points to enhanced trade winds and strengthened upwelling (Stein, 1985). Enhanced sapropel preservation, $\delta^{18}$O, Ba/Al, Ti/Al, and IRM records (ODP Sites 967 and 659) point toward an increase of fluvial-sourced terrigenous sediment into the Eastern Mediterranean beginning at ~3.3 Ma related to an enhanced Hadley circulation and more equatorward African monsoon front as a result of global cooling (Amarathunga et al., 2024; Grant et al., 2022; Heslop et al., 2023). Climate models of MIS M2 and early Late Pliocene glaciation additionally show enhanced summer subsidence in the IP when compared to background warmer and wetter Early—early Late Pliocene conditions (Amarathunga et al., 2024).

## 5.5. Mechanisms for Changes in Winter Rainfall

Consideration must also be given to the winter rainfall component across the late Pliocene and Pleistocene. The meridional temperature gradient affects not only the position of the Hadley circulation but also the midlatitude westerly jet (Ceppi et al., 2013). Modeling shows the mid-latitude westerly jet, which creates the storm track from its related baroclinic instability and lies just south of the track, shifted poleward by ~3.6° of its modern position and weakened during the Early Pliocene (Blackmon et al., 1977; Colleoni et al., 2015; Harvey et al., 2020; Li et al., 2015). In this respect, iNHG likely resulted in a stronger, equatorward shifted jet during the Late Pliocene-Early Pleistocene, which brought the North Atlantic Storm track closer to the Mediterranean region, leading to intensified winter rainfall compared to the early Late Pliocene glaciation (Figure 9). Climate models of the MIS M2 and mid-Pliocene glacials suggest this with strengthened westerly winds resulting in stronger winter westerlies (Amarathunga et al., 2024). Thus, while there is evidence that winter rainfall likely increased during the
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Schematic of possible mechanisms affecting the rainfall regime of the Iberian Peninsula during (a) the Early Pliocene and (b) across Mid-Pliocene glaciation. Purple and orange lines represent schematic shifts of the latitudes of the westerly jet and subsiding branch of the Hadley Cell, respectively, based on Pliocene model data (Li et al., 2015) and modern locations.</figcaption></figure>

iNHG due to an intensified storm track, given all other evidence for drying in the region locally we suggest summer drying was the dominant change in seasonality that influenced our $\delta^{18}O_c$ record.

**6. Conclusions**

We observe decreased $\delta^{18}O_c$ values from the IP during ~3.1 to ~2.5 Ma that corresponds with Pliocene global cooling. Through detailed comparison with other proxy records, we attribute this shift to a change in the seasonality of rainfall, specifically a decrease in the amount of $^{18}$O-enriched summer rainfall that led to regional drying. Such a change explains the concurrent increase in $\delta^{18}O_L$, shifts toward drought tolerant taxa, a longer-term increase in aridity observed in this region as reflected in increases in MOW, changes from lacustrine to palustrine to fluvial depositional environments in the Teruel Basin and changes in both magnetic susceptibility and the luminescence of carbonate nodules in the sediments of the Villalba Alta section. This change in rainfall seasonality is best explained by the cooling of Atlantic SSTs, intensified meridional SST gradient, and an equatorward shifted stronger subsiding branch of the Hadley cell from initial conditions of the “mid-Pliocene” warm interval, resulting in decreased moisture delivery to the basin. This study is a novel approach that provides climatic evidence that points to summer specific drying and related mechanisms for the poorly understood onset of the Mediterranean Basin's modern climate.

**Conflict of Interest**

The authors declare no conflicts of interest relevant to this study.

**Data Availability Statement**

Data is available through the PANGAEA Data Publisher under the title: Pedogenic carbonate nodule $\delta^{18}O$ and $\delta^{13}C$ record from the Pliocene-Pleistocene Villalba-Alta Section, Teruel (Czwakiel et al., 2025).
Acknowledgments

We are grateful to both Ana María Alonso-Zarza and Miguel Cisneros for their help in the field. We additionally thank Jeff Cullen for help in the laboratory. This work was supported by the National Science Foundation (PIRE Grant 1545859).

References

Aguirre Palafox, L. E., Möller, A., McLean, N. M., Ludvigson, G. A., Colombi, C. E., & Montañez, I. P. (2024). U-Pb geochronology of paleosol carbonate cements by LA-ICP-MS: A proof of concept and strategy for dating the terrestrial record. *Geochemistry, Geophysics, Geosystems*, 25(6), e2024GC011488. https://doi.org/10.1029/2024GC011488

Ahmed, I. A. M., & Maher, B. A. (2018). Identification and paleoclimatic significance of magnetite nanoparticles in soils. *Proceedings of the National Academy of Sciences of the United States of America*, 115(8), 1736–1741. https://doi.org/10.1073/pnas.1719186115

Alonso-Zarza, A. M., & Arenas, C. (2004). Cenozoic calcites from the Teruel Graben, Spain: Microstructure, stable isotope geochemistry and environmental significance. *Sedimentary Geology*, 167(1–2), 91–108. https://doi.org/10.1016/j.sedgeo.2004.02.001

Alonso-Zarza, A. M., & Calvo, J. P. (2000). Palustrine sedimentation in an episodically subsiding basin: The Miocene of the northern Teruel Graben (Spain). *Palaeogeography, Palaeoclimatology, Palaeoecology*, 160(1–2), 1–21. https://doi.org/10.1016/s0031-0182(00)00041-9

Alonso-Zarza, A. M., Meléndez, A., Martín-García, R., Herrero, M. J., & Martín-Pérez, A. (2012). Discriminating between tectonism and climate signatures in palustrine deposits: Lessons from the Miocene of the Teruel Graben, NE Spain. *Earth-Science Reviews*, 113, 141–160. https://doi.org/10.1016/j.earscirev.2012.03.011

Amarathunga, U., Rohling, E. J., Grant, K. M., Francke, A., Latimer, J., Klaebe, R. M., et al. (2024). Mid-Pliocene glaciation preceded by a 0.5-million-year North African humid period. *Nature Geoscience*, 17(7), 660–666. https://doi.org/10.1038/s41561-024-01472-8

Anadón, P., Cabrera, L., Julià, R., Roca, E., & Rosell, L. (1989). Lacustrine oil-shale basins in Tertiary Grabens from NE Spain (Western European rift system). *Palaeogeography, Palaeoclimatology, Palaeoecology*, 70(1–3), 7–28. https://doi.org/10.1016/0031-0182(89)90077-1

Anadón, P., & Moissenet, E. (1996). Neogene basins in the eastern Iberian range. Tertiary basins of Spain. In *The stratigraphic record of crustal kinematics* (pp. 68–76).

Barnaby, R. J., & Rimstidt, J. D. (1989). Redox conditions of calcite cementation interpreted from Mn and Fe contents of authigenic calcites. *Geological Society of America Bulletin*, 101(6), 795–804. https://doi.org/10.1130/0016-7606(1989)101<0795:RCOCCI>2.3.CO;2

Bellanca, A., Calvo, J. P., Censi, P., Neri, R., & Pozo, M. (1992). Recognition of lake-level changes in Miocene lacustrine units, Madrid Basin, Spain. Evidence from facies analysis, isotope geochemistry and clay mineralogy. *Sedimentary Geology*, 76(3–4), 135–153. https://doi.org/10.1016/0037-0738(92)90080-B

Belmecheri, S., Babst, F., Hudson, A. R., Betancourt, J., & Trouet, V. (2017). Northern Hemisphere jet stream position indices as diagnostic tools for climate and ecosystem dynamics. *Earth Interactions*, 21(8), 1–23. https://doi.org/10.1175/ei-d-16-0023.1

Beltran, C., Flores, J.-A., Sicre, M.-A., Baudin, F., Renard, M., & de Rafélis, M. (2011). Long chain alkenones in the Early Pliocene Sicilian sediments (Trubi Formation – Punta di Maiata section): Implications for the alkenone paleothermometry. *Palaeogeography, Palaeoclimatology, Palaeoecology*, 308(3–4), 253–263. https://doi.org/10.1016/j.palaeo.2011.03.017

Beltran, C., Sicre, M. A., Ohneiser, C., & Sainz, M. (2021). A composite Pliocene record of sea surface temperature in the central Mediterranean (Capo Rossello composite section – South Sicily). *Sedimentary Geology*, 420, 105921. https://doi.org/10.1016/j.sedgeo.2021.105921

Biltekin, D. (2010). *Vegetation and climate of North Anatolian and North Aegean region since 7 Ma according to pollen analysis*. Université Claude Bernard - Lyon I; İstanbul teknik üniversitesi.

Blackmon, M. L., Wallace, J. M., Lau, N., & Mullen, S. L. (1977). An observational study of the northern hemisphere wintertime circulation. *Journal of the Atmospheric Sciences*, 34, 1040–1053. https://doi.org/10.1175/1520-0469

Breecker, D. O., Sharp, Z. D., & McFadden, L. D., (2009). Seasonal bias in the formation and stable isotopic composition of pedogenic carbonate in modern soils from central New Mexico, USA. *GSA Bulletin*, 121(3–4), 630–640. https://doi.org/10.1130/B26413.1

Brierley, C. M., Fedorov, A. V., Liu, Z., Herbert, T. D., Lawrence, K. T., & LaRiviere, J. P. (2009). Greatly expanded tropical warm pool and weakened Hadley circulation in the early Pliocene. *Science*, 323(5922), 1714–1718. https://doi.org/10.1126/science.1167625

Broccoli, A. J., Dahl, K. A., & Stouffer, R. J. (2006). Response of the ITCZ to northern hemisphere cooling. *Geophysical Research Letters*, 33(1), 1–4. https://doi.org/10.1029/2005GL024546

Budd, D. A., Pack, S. M., & Fogel, M. L. (2002). The destruction of paleoclimatic isotopic signals in Pleistocene carbonate soil nodules of Western Australia. *Palaeogeography, Palaeoclimatology, Palaeoecology*, 188(3–4), 249–273. https://doi.org/10.1016/S0031-0182(02)00588-6

Burls, N. J., & Fedorov, A. V. (2017). Wetter subtropics in a warmer world: Contrasting past and future hydrological cycles. *Proceedings of the National Academy of Sciences of the United States of America*, 114(49), 12888–12893. https://doi.org/10.1073/pnas.1703421114

Cabos, W., de la Vará, A., Álvarez-García, F. J., Sánchez, E., Sieck, K., Pérez-Sanz, J. I., et al. (2020). Impact of ocean-atmosphere coupling on regional climate: The Iberian Peninsula case. *Climate Dynamics*, 54(9–10), 4441–4467. https://doi.org/10.1007/s00382-020-05238-x

Ceppi, P., Hwang, Y. T., Liu, X., Frierson, D. M. W., & Hartmann, D. L. (2013). The relationship between the ITCZ and the Southern Hemispheric eddy-driven jet. *Journal of Geophysical Research: Atmospheres*, 118(11), 5136–5146. https://doi.org/10.1002/jgrd.50461

Cerling, T. E. (1984). The stable isotopic composition of modern soil carbonate and its relationship to climate. *Earth and Planetary Science Letters*, 71(2), 229–240. https://doi.org/10.1016/0012-821X(84)90089-X

Cerling, T. E., & Quade, J. (1993). Stable carbon and oxygen isotopes in soil carbonates. In P. K. Swart, K. C. Lohmann, J. Mckenzie, & S. Savin (Eds.), *Climate change in continental isotopic records*. https://doi.org/10.1029/GM078p0217

Chiang, J. C. H., Biasutti, M., & Battisti, D. S. (2003). Sensitivity of the Atlantic Intertropical Convergence zone to last glacial maximum boundary conditions. *Paleoceanography*, 18(4), 1–18. https://doi.org/10.1029/2003PA000916

Colleoni, F., Cherchi, A., Masina, S., & Brierley, C. M. (2015). Impact of global SST gradients on the Mediterranean runoff changes across the Plio-Pleistocene transition. *Paleoceanography*, 30(6), 751–767. https://doi.org/10.1002/2015PA002780

Colleoni, F., Masina, S., Negri, A., & Marzocchi, A. (2012). Plio-Pleistocene high-low latitude climate interplay: A Mediterranean point of view. *Earth and Planetary Science Letters*, 319–320, 35–44. https://doi.org/10.1016/j.epsl.2011.12.020

Coplen, T. B. (2007). Calibration of the calcite-water oxygen-isotope geothermometer at Devils Hole, Nevada, a natural laboratory. *Geochimica et Cosmochimica Acta*, 71(16), 3948–3957. https://doi.org/10.1016/j.gca.2007.05.028

Czwakiel, N., Gallagher, T. M., Serach, L., Ludvigson, G., Gao, P., Nie, J., et al. (2025). Pedogenic carbonate nodule $\delta^{18}$O and $\delta^{13}$C record from the Pliocene-Pleistocene Villalba-Alta section, Teruel [Dataset]. PANGAEA. https://doi.org/10.1594/PANGAEA.979978

Da, J., Li, G. K., & Ji, J. (2023). Seasonal changes in the formation time of pedogenic carbonates on the Chinese Loess Plateau during Quaternary glacial cycles. *Quaternary Science Reviews*, 305, 108008. https://doi.org/10.1016/j.quascirev.2023.108008

Dansgaard, W. (1964). Stable isotopes in precipitation. *Tellus*, 16(4), 436–468. https://doi.org/10.1111/j.2153-3490.1964.tb00181.x

De Schepper, S., Gibbard, P. L., Salzmann, U., & Ehlers, J. (2014). A global synthesis of the marine and terrestrial evidence for glaciation during the Pliocene Epoch. *Earth-Science Reviews*, 135, 83–102. https://doi.org/10.1016/j.earscirev.2014.04.003

De Schepper, S., Groeneveld, J., Naafs, B. D. A., Van Renterghem, C., Hennissen, J., Head, M. J., et al. (2013). Northern Hemisphere Glaciation during the globally warm early late Pliocene. *PLoS One*, 8(12), e81508. https://doi.org/10.1371/journal.pone.0081508
Deutz, P., Montañez, I. P., Monger, H. C., & Morrison, J. (2001). Morphology and isotope heterogeneity of Late Quaternary pedogenic carbonates: Implications for paleosol carbonates as paleoenvironmental proxies. *Palaeogeography, Palaeoclimatology, Palaeoecology, 166*(3-4), 293-317. https://doi.org/10.1016/S0031-0182(00)00214-5

Dolan, A. M., Haywood, A. M., Hunter, S. J., Tindall, J. C., Dowsett, H. J., Hill, D. J., & Pickering, S. J. (2015). Modelling the enigmatic late Pliocene glacial event - Marine isotope stage M2. *Global and Planetary Change, 128*, 47–60. https://doi.org/10.1016/j.gloplacha.2015.02.001

Dowsett, H. J., Haywood, A. M., Valdes, P. J., Robinson, M. M., Lunt, D. J., Hill, D. J., et al. (2011). Sea surface temperatures of the mid-Piacenzian warm period: A comparison of PRISM3 and HadCM3. *Palaeogeography, Palaeoclimatology, Palaeoecology, 309*(1-2), 83-91. https://doi.org/10.1016/j.palaeo.2011.03.016

Dworkin, S. I., Nordt, L., & Atchley, S. (2005). Determining terrestrial paleotemperatures using the oxygen isotopic composition of Pedogenic carbonate. *Earth and Planetary Science Letters, 237*(1-2), 56-68. https://doi.org/10.1016/j.epsl.2005.06.054

Ezquerro, L. (2017). *El sector norte de la cuenca neógena de Teruel: tectónica, clima y sedimentación* (Vol. 496). Universidad de Zaragoza. Retrieved from http://zaguan.unizar.es/record/77098

Ezquerro, L., Luzón, A., Liesa, C. L., & Simón, J. L. (2012). Megasequential evolution of the Mio-Pliocene filling of the northern sector of the Teruel basin. In *VIII geological congress of Spain, Oviedo* (pp. 218-221).

Ezquerro, L., Luzón, A., Navarro, M., Liesa, C. L., & Simón, J. L. (2014). Climatic vs. tectonic signals in a continental extensional basin (Teruel, NE Spain) from stable isotope (δ¹⁸O) and sequence stratigraphical evolution. *Terra Nova, 26*(5), 337-346. https://doi.org/10.1111/ter.12101

Ezquerro, L., Luzón, A., Simón, J. L., & Liesa, C. L. (2019). Alluvial sedimentation and tectono-stratigraphic evolution in a narrow extensional zigzag basin margin (northern Teruel Basin, Spain). *Journal of Palaeogeography, 8*(1), 29. https://doi.org/10.1186/s42501-019-0044-4

Ezquerro, L., Moretti, M., Liesa, C. L., Luzón, A., Pueyo, E. L., & Simón, J. L. (2016). Controls on space-time distribution of soft-sediment deformation structures: Approaching the apparent recurrence period of paleoseisms at the Concud fault (eastern Spain). *Sedimentary Geology, 344*, 91-111. https://doi.org/10.1016/j.sedgeo.2016.06.007

Ezquerro, L., Muñoz, A., Liesa, C. L., Simón, J. L., & Luzón, A. (2022). Late Neogene to Early Quaternary climate evolution in southwestern Europe from a continental perspective. *Global and Planetary Change, 211* (January), 103788. https://doi.org/10.1016/j.gloplacha.2022.103788

Fauquette, S., Suc, J.-P., Guiot, J., Diniz, F., Feddi, N., Zheng, Z., et al. (1999). Climate and biomes in the west Mediterranean area during the Pliocene. *Palaeogeography, Palaeoclimatology, Palaeoecology, 152*(1-2), 15-36. https://doi.org/10.1016/S0031-0182(99)00031-0

Gao, P., Nie, J., Breecker, D. O., Gallagher, T., Serach, L., & Alonso-Zarza, A. M. (2022). Similar magnetic enhancement mechanisms between Chinese loess and alluvial sediments from the Teruel Basin, NE Spain, and paleoclimate implications. *Geophysical Research Letters, 49*(6), e2021GL096977. https://doi.org/10.1029/2021GL096977

Gradstein, F. M., Ogg, J. G., Schmitz, M. D., & Ogg, G. M. (2020). *Geologic time scale*. Elsevier.

Grant, K. M., Amarathunga, U., Amies, J. D., Hu, P., Qian, Y., Penny, T., et al. (2002). Organic carbon burial in Mediterranean sapropels intensified during Green Sahara Periods since 3.2 Myr ago. *Communications Earth & Environment, 3*(1), 11. https://doi.org/10.1038/s43247-021-00339-9

Harvey, B. J., Cook, P., Shaffrey, L. C., & Schiemann, R. (2020). The response of the northern hemisphere storm tracks and jet streams to climate change in the CMIP3, CMIP5, and CMIP6 climate models. *Journal of Geophysical Research: Atmospheres, 125*(23), e2020JD032701. https://doi.org/10.1029/2020JD032701

Haywood, A. M., Sellwood, B. W., & Valdes, P. J. (2000). Regional warming: Pliocene (3 Ma) paleoclimate of Europe and the Mediterranean. *Geology, 28*(12), 1063-1066. https://doi.org/10.1130/0091-7613(2000)28<1063:RWPMP0>2.0.CO;2

Haywood, A. M., & Valdes, P. J. (2004). Modelling Pliocene warmth: Contribution of atmosphere, oceans and cryosphere. *Earth and Planetary Science Letters, 218*(3-4), 363-377. https://doi.org/10.1016/S0012-821X(03)00685-X

Hemming, G. N., Meyers, W. J., & Grams, J. C. (1989). Cathodoluminescence in diagenetic calcites: The roles of Fe and Mn as deduced from electron probe and spectrophotometric measurements. *Journal of Sedimentary Petrology, 59*(3), 404-411.

Hennissen, J. A. I., Head, M. J., De Schepper, S., & Groeneveld, J. (2015). Increased seasonality during the intensification of northern hemisphere glaciation at the pliocene-pleistocene boundary ~2.6 Ma. *Quaternary Science Reviews, 129*, 321-332. https://doi.org/10.1016/j.quascirev.2015.10.010

Herbert, T. D., Gideong, N., & Peterson, L. C. (2015). Evolution of Mediterranean Sea surface temperatures 3.5-1.5 ma: Regional and hemispheric influences. *Earth and Planetary Science Letters, 409*, 307-318. https://doi.org/10.1016/j.epsl.2014.10.006

Herbert, T. D., Lawrence, K. T., Tzanova, A., Peterson, L. C., Caballero-Gill, R. P., & Kelly, C. S. (2018). (Table S2) SST estimates as a function of age, Mediterranean sites. *PANGAEA*. https://doi.org/10.1594/PANGAEA.885603

Herrera, S., Margarida Cardoso, R., Matos Soares, P., Espírito-Santo, F., Viterbo, P., & Gutiérrez, J. M. (2019). Iberia01: A new gridded dataset of daily precipitation and temperatures over Iberia. *Earth System Science Data, 11*(4), 1947-1956. https://doi.org/10.5194/essd-11-1947-2019

Heslop, D., Amarathunga, U., & Rohling, E. J. (2023). Estimating Plio-Pleistocene North African monsoon runoff into the Mediterranean Sea and temperature impacts. *Paleoceanography and Paleoclimatology, 38*(11), e2023PA004677. https://doi.org/10.1029/2023pa004677

Horton, T. W., Defliese, W. F., Tripati, A. K., & Oze, C. (2016). Evaporation induced ¹⁸O and ¹³C enrichment in lake systems: A global perspective on hydrologic balance effects. *Quaternary Science Reviews, 131*, 365-379. https://doi.org/10.1016/j.quascirev.2015.06.030

IAEA/WMO. (2024). *Global network of isotopes in precipitation*. International Atomic Energy Agency. Retrieved from https://nucleus.iaea.org/wiser/

Jiménez-Moreno, G., Fauquette, S., & Suc, J.-P. (2010). Miocene to Pliocene vegetation reconstruction and climate estimates in the Iberian Peninsula from pollen data. *Review of Palaeobotany and Palynology, 162*(3), 403-415. https://doi.org/10.1016/j.revpalbo.2009.08.001

Kelson, J. R., Huth, T. E., Passey, B. H., Levin, N. E., Petersen, S. V., Ballato, P., et al. (2023). Triple oxygen isotope compositions of globally distributed soil carbonates record widespread evaporation of soil waters. *Geochimica et Cosmochimica Acta, 355*, 138-160. https://doi.org/10.1016/j.gca.2023.06.034

Khélifi, N., Sarnthein, M., Frank, M., Andersen, N., & Garbe-Schönberg, D. (2014). Late Pliocene variations of the Mediterranean outflow. *Marine Geology, 357*, 182-194. https://doi.org/10.1016/j.margeo.2014.07.006

Kim, S. T., & O'Neil, J. R. (1997). Equilibrium and nonequilibrium oxygen isotope effects in synthetic carbonates. *Geochimica et Cosmochimica Acta, 61*(16), 3461-3475. https://doi.org/10.1016/S0016-7037(97)00169-5

Koppen, W. (1936). Das geographische System der Klimate. In W. Koppen & R. Geiger (Eds.), *Handbuch der Klimatologie 3* (p. 46). Gebrüder Borntraeger.

Lafuente, P., Arlegui, L. E., Liesa, C. L., & Simón, J. L. (2010). Paleoseismological analysis of an intraplate extensional structure: The Concud fault (Iberian Chain, eastern Spain). *International Journal of Earth Sciences, 100*(7), 1713-1732. https://doi.org/10.1007/s00531-010-0542-1

Layzell, A. L., Ludvigson, G. A., Smith, J. J., & Mandel, R. D. (2023). Using the factors of soil formation to assess stable carbon isotope disequilibrium in late Pleistocene (MIS 3) buried soils of the Great Plains, North America. *Palaeogeography, Palaeoclimatology, Palaeoecology, 620*, 111574. https://doi.org/10.1016/j.palaeo.2023.111574
Leng, M. J., & Marshall, J. D. (2004). Paleoclimate interpretation of stable isotope data from lake sediment archives. *Quaternary Science Reviews*, 23(7-8), 811–831. https://doi.org/10.1016/j.quascirev.2003.06.012

Leroy, S. A. G., & Dupont, L. M. (1997). Marine palynology of the ODP site 658 (NW Africa) and its contribution to the stratigraphy of late Pliocene. *Geobios*, 30(3), 351–359. https://doi.org/10.1016/s0016-6995(97)80194-5

Li, X. Y., Jiang, D. B., Zhang, Z. S., Zhang, R., Tian, Z. P., & Yan, Q. (2015). Mid-Pliocene westerlies from PlioMIP simulations. *Advances in Atmospheric Sciences*, 32(7), 909–923. https://doi.org/10.1007/s00376-014-4171-7

Liesa, C. L., Simon, J. L., Ezquerro, L., Arlegui, L. E., & Luzón, A. (2019). Stress evolution and structural inheritance controlling an intra-continental extensional basin: The central-northern sector of the Neogene Teruel Basin. *Journal of Structural Geology*, 118, 362–376. https://doi.org/10.1016/j.jsg.2018.11.011

Lionello, P., Malanotte-Rizzoli, P., Boscolo, R., Alpert, P., Artale, V., Li, L., et al. (2006). The Mediterranean climate: An overview of the main characteristics and issues. *Developments in earth and environmental sciences*, 4, 1–26. https://doi.org/10.1016/s1571-9197(06)80003-0

Lisiecki, L. E., & Raymo, M. E. (2005). A Pliocene-Pleistocene stack of 57 globally distributed benthic $\delta^{18}$O records. *Paleoceanography*, 20(1), 1–17. https://doi.org/10.1029/2004PA001071

Loubere, P., & Moss, K. (1986). Late Pliocene climatic change and the onset of Northern Hemisphere glaciation as recorded in the northern Atlantic Ocean. *Geological Society of America Bulletin*, 97(7), 818–828. https://doi.org/10.1130/0016-7606(1986)97<818:LPCCAT>2.0. CO;2

Lourens, L. J. (2004). Revised tuning of Ocean Drilling Program Site 964 and KC01B (Mediterranean) and implications for the $\delta^{18}$O, tephra, calcareous nannofossil, and geomagnetic reversal chronologies of the past 1.1 Myr. *Paleoceanography*, 19(3), 1–20. https://doi.org/10.1029/ 2003PA000997

McClymont, E. L., Ho, S. L., Ford, H. L., Bailey, I., Berke, M. A., Bolton, C. T., et al. (2023). Climate evolution through the onset and intensification of northern hemisphere glaciation. *Reviews of Geophysics*, 61(3), e2022RG000793. https://doi.org/10.1029/2022RG000793

Médail, F., & Quézel, P. (1997). Hot-spots analysis for conservation of plant biodiversity in the Mediterranean Basin. *Annals of the Missouri Botanical Garden*, 84(1), 112–127. https://doi.org/10.2307/239957

Moissenet, E. (1989). The Neogene basins of the Iberian chain: A comparative analysis. *Bulletin of the Geological Society of France*, 5, 919–926. https://doi.org/10.2113/gssgfbull.V.5.919

Moreno, A., Iglesias, M., Azorin-Molina, C., Pérez-Meiá, C., Bartolomé, M., Stoll, H., et al. (2021). Measurement report: Spatial variability of northern Iberian rainfall stable isotope values - Investigating atmospheric controls on daily and monthly timescales. *Atmospheric Chemistry and Physics*, 21(13), 10159–10177. https://doi.org/10.5194/acp-21-10159-2021

Mudelsee, M., & Raymo, M. E. (2005). Slow dynamics of the northern hemisphere glaciation. *Paleoceanography*, 20(4). https://doi.org/10.1029/ 2005PA001153

Naafs, B. D. A., Voelker, A. H. L., Karas, C., Andersen, N., & Sierro, F. J. (2020). Repeated near-collapse of the Pliocene sea surface temperature gradient in the North Atlantic. *Paleoceanography and Paleoclimatology*, 35(5), 1–16. https://doi.org/10.1029/2020PA003905

Ochoa, D., Sierro, F. J., Hilgen, F. J., Cortina, A., Lofi, J., Kouwenhoven, T., & Flores, J.-A. (2018). Origin and implications of orbital-induced sedimentary cyclicity in Pliocene well-logs of the Western Mediterranean. *Marine Geology*, 403, 150–164. https://doi.org/10.1016/j.margeo. 2018.05.009

Opdyke, N., Mein, P., Lindsay, E., Perez-Gonzales, A., Moissenet, E., & Norton, V. L. (1997). Continental deposits, magnetostratigraphy and vertebrate paleontology, late Neogene of eastern Spain. *Palaeogeography, Palaeoclimatology, Palaeoecology*, 133(3-4), 129–148. https://doi. org/10.1016/S0031-0182(97)00080-1

Popescu, S. M., Biltekin, D., Winter, H., Suc, J.-P., Melinte-Dobrinescu, M. C., Klotz, S., et al. (2010). Pliocene and lower Pleistocene vegetation and climate changes at the European scale: Long pollen records and climatostratigraphy. *Quaternary International*, 219(1-2), 152–167. https:// doi.org/10.1016/j.quaint.2010.03.013

Raymo, M. E. (1994). The initiation of Northern Hemisphere glaciation. *Annual Review of Earth and Planetary Sciences*, 22(1), 353–383. https:// doi.org/10.1146/annurev.ea.22.050194.002033

Richter, D. K., Götte, T., Götzé, J., & Neuser, R. D. (2003). Progress in application of cathodoluminescence (CL) in sedimentary petrology. *Mineralogy and Petrology*, 79(3-4), 127–166. https://doi.org/10.1007/s00710-003-0237-4

Roca, E. (1994). La evolución geodinámica de la Cuenca Catalano-Balear y áreas adyacentes desde el Mesozoico hasta la actualidad. *Acta Geológica Hispánica*, 3-25.

Rodríguez-López, J. P., Liesa, C. L., Van Dam, J., Lafuente, P., Arlegui, L., Ezquerro, L., & De Boer, P. L. (2012). Aeolian construction and alluvial dismantling of a fault-bounded intracontinental aeolian dune field (Teruel Basin, Spain); a continental perspective on Late Pliocene climate change and variability. *Sedimentology*, 59(5), 1536–1567. https://doi.org/10.1111/j.1365-3091.2011.01316.x

Rodwell, M. J., & Hoskins, B. J. (1996). Monsoons and the dynamics of deserts. *Quarterly Journal of the Royal Meteorological Society*, 122(534), 1385–1404. https://doi.org/10.1002/qj.49712253408

Rohling, E. J., Yu, J., Heslop, D., Foster, D. L., Opdyke, B., & Roberts, A. P. (2021). Sea level and deep-sea temperature reconstructions suggest quasi-stable states and critical transitions over the past 40 million years. *Science Advances*, 7(26), eabf5326. https://doi.org/10.1126/sciadv. abf5326

Ross, J. B., Ludvigson, G. A., Möller, A., González, L. A., & Walker, J. D. (2017). Stable isotope paleohydrology and chemostratigraphy of the Albian Wayan formation from the wedge-top depozone, North American western Interior basin. *Science China Earth Sciences*, 60(1), 44–57. https://doi.org/10.1007/s11430-016-0087-5

Salzmann, U., Williams, M., Haywood, A. M., Johnson, A. L., Kender, S., & Zalasiewicz, J. (2011). Climate and environment of a Pliocene warm world. *Palaeogeography, Palaeoclimatology, Palaeoecology*, 309(1-2), 1–8. https://doi.org/10.1016/j.palaeo.2011.05.044

Serrano, A., García, J. A., Mateos, V. L., Cancillo, M. L., & Garrido, J. (1999). Monthly modes of variation of precipitation over the Iberian Peninsula. *Journal of Climate*, 12(9), 2894–2919. https://doi.org/10.1175/1520-0442(1999)012<2894:MMOVOP>2.0.CO;2

Simón, J. L., Ezquerro, L., Arlegui, L. E., Liesa, C. L., Luzón, A., Medialdea, A., et al. (2019). Role of transverse structures in paleoseismicity and drainage rearrangement in rift systems: The case of the Valdecebro fault zone (Teruel Graben, eastern Spain). *International Journal of Earth Sciences*, 108(5), 1429–1449. https://doi.org/10.1007/s00531-019-01707-9

Spötl, C., & Vennemann, T. W. (2003). Continuous-flow isotope ratio mass spectrometric analysis of carbonate minerals. *Rapid Communications in Mass Spectrometry*, 17(9), 1004–1006. https://doi.org/10.1002/rcm.1010

Stein, R. (1985). Late Neogene changes of paleoclimate and paleoproduction off northwest Africa (DSDP site 397). *Palaeogeography, Palaeoclimatology, Palaeoecology*, 49(1-2), 47–59. https://doi.org/10.1016/0031-0182(85)90004-5

Suc, J.-P. (1984). Origin and evolution of the Mediterranean vegetation and climate in Europe. *Nature*, 307(5950), 429–432. https://doi.org/10. 1038/307429a0
Suc, J.-P., Popescu, S. M., Fauquette, S., Bessedik, M., Moreno, G. J., Taoufiq, N. B., et al. (2018). Reconstruction of Mediterranean flora, vegetation and climate for the last 23 million years based on an extensive pollen dataset. *Ecologia Mediterranea, 44*(2), 53–85. https://doi.org/10.3406/ecmed.2018.2044

Sun, Y., Ramstein, G., Contoux, C., & Zhou, T. (2013). A comparative study of large-scale atmospheric circulation in the context of a future scenario (RCP4.5) and past warmth (mid-Pliocene). *Climate of the Past, 9*(4), 1613–1627. https://doi.org/10.5194/cp-9-1613-2013

Tiedemann, R., Sarnthein, M., & Shackleton, N. J. (1994). Astronomic timescale for the Pliocene Atlantic $\delta^{18}$O and dust flux records of ocean drilling program site 659. *Paleoceanography, 9*(4), 619–638. https://doi.org/10.1029/94pa00208

van Dam, J. A., Abdul Aziz, H., Angeles Alvarez Sierra, M., Hilgen, F. J., van den Hoek Ostende, L. W., Lourens, L. J., et al. (2006). Long-period astronomical forcing of mammal turnover. *Nature, 443*(7112), 687–691. https://doi.org/10.1038/nature05163

van Dam, J. A., Mein, P., Garcés, M., van Balen, R. T., Furió, M., & Alcalá, L. (2023). Macroevolutionary and macroecological response of Iberian rodents to late Neogene climatic oscillations and events. *Global and Planetary Change, 227*, 104153. https://doi.org/10.1016/j.gloplacha.2023.104153

van Dam, J. A., & Reichart, G. J. (2009). Oxygen and carbon isotope signatures in late Neogene horse teeth from Spain and application as temperature and seasonality proxies. *Palaeogeography, Palaeoclimatology, Palaeoecology, 274*(1–2), 64–81. https://doi.org/10.1016/j.palaeo.2008.12.022

van Dam, J. A., & Weltje, G. J. (1999). Reconstruction of the late Miocene climate of Spain using rodent palaeocommunity successions: An application of end-member modelling. *Palaeogeography, Palaeoclimatology, Palaeoecology, 151*(4), 267–305. https://doi.org/10.1016/s0031-0182(99)00015-2

Wang, P., Tian, J., & Lourens, L. J. (2010). Obscuring of long eccentricity cyclicity in Pleistocene oceanic carbon isotope records. *Earth and Planetary Science Letters, 290*(3–4), 319–330. https://doi.org/10.1016/j.epsl.2009.12.028

Wieder, M., & Yaalon, D. H. (1982). Micromorphological fabrics and developmental stages of carbonate nodular forms related to soil characteristics. *Geoderma, 28*(3–4), 203–220. https://doi.org/10.1016/0016-7061(82)90003-9

Zachos, J., Pagani, H., Sloan, L., Thomas, E., & Billups, K. (2001). Trends, rhythms, and aberrations in global climate 65 Ma to present. *Science, 292*(5517), 686–693. https://doi.org/10.1126/science.1059412

## References From the Supporting Information

Ogg, J. G. (2020). Chapter 5 – Geomagnetic polarity time scale. In F. M. Gradstein, J. G. Ogg, M. D. Schmitz, & G. M. Ogg (Eds.), *Geologic timescale 2020* (pp. 159–192). Elsevier. https://doi.org/10.1016/B978-0-12-824360-2.00005-X