Temperature variability in southern Europe over the past 16 500 years constrained by speleothem fluid inclusion water isotopes

Abstract

In the Northern Hemisphere, the last 16.5 $\mathrm{kyr}$ was characterized by abrupt temperature transitions during stadials, interstadials, and the onset of the Holocene. These changes are closely linked to large-scale variations in the extent of continental ice sheets, greenhouse gas concentrations, and ocean circulation. Speleothems and their fluid inclusions serve as valuable proxies, offering high-resolution chronologies and quantitative records of past temperature changes for understanding global and regional climate mechanisms in the past. Here, we present a record based on five speleothems from two caves on the northeastern Iberian Peninsula (Ostolo and Mendukilo caves). Using hydrogen isotopic composition of fluid inclusions and rainfall samples, we developed a ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}/\mathrm{T}$ transfer function in order to reconstruct regional temperatures over the past 16.5 $\mathrm{kyr}$ (Ostolo–Mendukilo Fluid Inclusion Temperature record, OM-FIT). Our novel findings reveal abrupt temperature changes in SW Europe during the last deglaciation and Early Holocene, at millennial and centennial scales, anchored by a precise chronology. At the onset of Greenland Interstadial 1, the OM-FIT record shows an increase of 6.7 $\pm$ 2.8 $\mathrm{°}\mathrm{C}$ relative to the cold conditions of the preceding Greenland Stadial 2.1a. During the early phase of Greenland Stadial 1, OM-FIT records a temperature decline of 6.1 $\pm$ 2.8 $\mathrm{°}\mathrm{C}$. The end of this cold phase and the onset of the Holocene are marked by a rapid warming of about 5 $\mathrm{°}\mathrm{C}$, reaching a maximum at 11.66 $\pm$ 0.03 $\mathrm{kyr}\mathrm{BP}$. The OM-FIT record also exhibits abrupt events during the Holocene (e.g., the 8.2 $\mathrm{kyr}$ event), which are also reflected in the ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ values of the calcite. These abrupt temperature changes during the last deglaciation and the Holocene correspond to variations seen in paleotemperature records across Europe and in Greenland ice cores. This clearly and quantitatively illustrates the influence of changes in the Atlantic Meridional Overturning Circulation, driven by subarctic freshening, on the climate of southern Europe.

1 Introduction

The last deglaciation in the Northern Hemisphere (ca. 16.5–11.7 thousand years ($\mathrm{kyr}$) BP – before present; present $=$ 1950) was punctuated by a series of abrupt climatic changes driven by variations in the extent of large continental ice sheets, greenhouse gas concentrations, and deep-water ocean circulation (Clark et al., 2012). The Holocene was also characterized by variability in terms of temperature, precipitation seasonality, and glacier extent (Wanner et al., 2008, 2011), albeit at much smaller amplitudes compared to the Late Pleistocene. Reconstructing such paleoclimate changes quantitatively poses significant challenges due to the scarcity of quantitative techniques and the fact that proxy signals in archives may be influenced by more than one meteorological variable (e.g., temperature and precipitation), which complicates our understanding of past temperature variations (Heiri et al., 2014a; Moreno et al., 2014). These limitations greatly hinder the assessment of whether reconstructed paleotemperatures in different regions are reflecting climate variations or different methodologies. Therefore, it is crucial to obtain proxy data that accurately reflect quantitative changes in paleotemperature, independent of past changes in rainfall or humidity. Quantitative temperature reconstructions are needed to assess the ability of climate simulation models to predict regionally divergent trends in climate change and to better understand the mechanisms of global and regional climate variability (e.g., Affolter et al., 2019).

The last deglaciation in the Northern Hemisphere involved major climatic shifts associated with Greenland Stadials (GS-2.1a and GS-1) and Interstadials (GI-1 and the onset of the Holocene). The impact of these rapid climate changes on southwestern European environments is recorded by temperature proxies, e.g., pollen, speleothems, planktonic foraminifera, and chironomids (Català et al., 2019; Cheng et al., 2020; González-Sampériz et al., 2017; Heiri et al., 2014b; Millet et al., 2012; Tarrats et al., 2018). However, the available temperature reconstructions exhibit large regional climate differences across Europe (Affolter et al., 2019; Heiri et al., 2014b; Renssen and Isarin, 2001). For example, the chironomid study by Heiri et al. (2014b) revealed that temperature variations during the last deglaciation were more pronounced in western Europe than in southwestern, central, and southeastern Europe. Similar regional disparities are observed during the Holocene, where the long-term evolution of global and hemispheric temperature variations remains a subject of debate, with climate models and proxy records showing differing trends (Affolter et al., 2019; Marcott et al., 2013; Shakun, 2018). Given these uncertainties, quantitative studies using inorganic archives, such as fluid inclusions (FIs) in speleothems (Demény et al., 2016, 2021; Dublyansky and Spötl, 2009), are gaining increasing relevance (Affolter et al., 2019; Honiat et al., 2023; Wilcox et al., 2020) as a complement to existing studies largely based on biological archives. The strengths of this method are (a) the accurate and precise chronology provided by speleothems; (b) the well-established link between cave interior temperature and mean outside air temperature; and (c) the relationship between temperature and water isotopes, which is controlled by physical rather than biological processes. FI water isotopes can be measured using different analytical techniques (Affolter et al., 2014; Arienzo et al., 2013; Dublyansky and Spötl, 2009; Vonhof et al., 2006) and FI-based paleotemperature reconstruction methods (Demény et al., 2016, 2021; Uemura et al., 2020). One such approach is based on the ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ composition and uses the ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$–temperature relationship determined for a given study area (Affolter et al., 2019). The principal advantage of this method lies in its reliance on a relatively simple and robust analytical method. The ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$–temperature relationship is established using monitoring data, and the approach is most effective in settings where ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$ variability in rainfall is driven by surface temperature (Demény et al., 2021). For the FI water isotope thermometry method to yield reliable results, four aspects must be considered: (i) FIs must be of primary origin, well-sealed, and sufficiently abundant; (ii) the choice of the transfer function converting the hydrogen and/or oxygen isotope signal (${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$, ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{FI}}$) into temperature may bias temperature estimates; (iii) the relationship between ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ and ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{FI}}$ may have changed over time; and (iv) the ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{FI}}$ water isotope method assumes that speleothem calcite was deposited under isotopic equilibrium conditions, where isotopic non-equilibrium can be attributed to either recrystallization or kinetic isotope fractionation during calcite precipitation. The ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ values are therefore the most suitable temperature proxy as they witness surface air temperature avoiding possible alterations during carbonate precipitation (Affolter et al., 2019; Demény et al., 2021; Uemura et al., 2020).

Figure 1

Location of the study area in northern Spain. Yellow stars indicate the locations of the two studied caves, while the red star marks the site where the isotopic composition of rainfall was monitored (Las Güixas tourist cave in Villanúa).

[Figure omitted. See PDF]

Here, we assess the air temperature evolution in the northern Iberian Peninsula over the last 16.5 $\mathrm{kyr}$ using quantitative FI-based data from five well-dated stalagmites that overlap during the last deglaciation and Holocene, showing very similar stable isotope trends. This record (dubbed OM-FIT, Ostolo–Mendukilo Fluid Inclusion Temperature) in conjunction with other regional terrestrial proxy records allows us to better disentangle the effects of temperature and humidity reported by previous studies using calcite stable isotope data from caves in southwestern Europe (Bernal-Wormull et al., 2021, 2023). The paleotemperature data obtained from FIs in speleothems represent the first quantitative air temperature reconstruction for northeastern Iberia during the last deglaciation and provide a basis for future studies aiming to enhance our quantitative understanding of rapid regional climate changes.

2 Study sites

2.1 Ostolo and Mendukilo caves

Ostolo (43°11^′16^′′ N, 1°43^′56^′′ W; 248 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) and Mendukilo (42°58^′25^′′ N, 1°53^′45^′′ W; 750 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) caves are located in northern Iberia (Fig. ). Although only about 28 $\mathrm{km}$ apart, they exhibit different geological, geomorphological, and climatic settings. Ostolo cave is located in the Bidasoa river valley, formed within the Carboniferous limestones of the Cinco Villas Massif (Basque Mountains, Western Pyrenees). Mendukilo cave, on the other hand, is developed in Lower Cretaceous limestones (Urgonian, Albian-Aptian) along the eastern boundary of the Basque–Cantabrian basin. For additional details on the caves and the locations of the sampled stalagmites, see Bernal-Wormull et al. (2021, 2023).

The climate in the study region is dominated by the Atlantic Ocean, characterized by temperate summers, evenly distributed rainfall throughout the year, and no distinct dry season (Cfb of the Köppen–Geiger climate classification). Mediterranean fronts may also be secondarily responsible for rainfall. Mean annual air temperature (MAAT) and mean annual precipitation are higher in the Ostolo cave area (13.5 $\pm$ 0.8 $\mathrm{°}\mathrm{C}$; $>$ 2000 $\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$) compared to Mendukilo (12.2 $\pm$ 0.4 $\mathrm{°}\mathrm{C}$; $\sim$ 1365 $\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$). This temperature difference is even more pronounced inside the caves: the average annual cave air temperature in Ostolo is 13 $\mathrm{°}\mathrm{C}$, while in Mendukilo, it is 8.8 $\mathrm{°}\mathrm{C}$. The lower temperature inside Mendukilo is due to its more closed and hence less ventilated nature compared to Ostolo, which also contains a cave stream that helps stabilize its internal temperature (Bernal-Wormull et al., 2021). In contrast, the “cold-trap” behavior of Mendukilo is consistent with its more complex geometry, resulting in an anomalously low temperature (Bernal-Wormull et al., 2023). The vegetation around both caves is dominated by oak (Quercus robur and Quercus pyrenaica); alder (Alnus glutinosa); beech (Fagus sylvatica); and Atlantic-type polycultures, ferns, and heathers.

2.2 Isotopic composition of drip waters in Ostolo and Mendukilo

Quantitative reconstruction of past climate variability from speleothem isotope records relies on understanding the modern vadose karst flow regime (Lachniet, 2009). For Mendukilo, the ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ and ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$ values of drip waters feeding the stalagmites studied here remain relatively constant, with mean values of $-$7.7 $\pm$ 0.4 ‰ and $-$45.3 $\pm$ 2.9 ‰ Vienna Standard Mean Ocean Water (VSMOW), respectively (1$\mathit{\sigma }$ uncertainty), and lack of a seasonal pattern (Bernal-Wormull et al., 2023). The monitoring period in Mendukilo cave lasted nearly 3 years, with measurements taken every 2–3 months (2018–2021). In Ostolo, the ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ and ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$ values of drip water are also similarly stable, with mean values of $-$6.3 $\pm$ 0.2 ‰ and $-$37.8 $\pm$ 1.6 $\mathrm{‰}\mathrm{VSMOW}$, respectively, with carbonate precipitation throughout the year in only one gallery of the cave (Bernal-Wormull et al., 2021). The monitoring interval in Ostolo was 3–4 months over 1 year (2019–2020).

Figure 2

(a) ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$ and ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ values of precipitation events at Villanúa (black dots) along with the Local Meteoric Water Line (LMWL; black line). Samples that experienced evaporation prior to sampling and outliers were excluded (Giménez et al., 2021). The Global Meteoric Water Line (GMWL; dashed line; Rozanski et al., 1993) and the drip water lines of Mendukilo (MEN; red line) and Ostolo (OST; green line) are also represented. (b) ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ values of precipitation events and their respective temperature at Villanúa (Giménez et al., 2021). (c) ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$ values of precipitation events and their respective temperature at Villanúa (Giménez et al., 2021). The dashed line represents the linear regression of precipitation isotope and temperature data.

[Figure omitted. See PDF]

2.3 Isotopic composition of rainfall

The rainfall stable isotopic composition near the study sites was analyzed by Giménez et al. (2021) on an event basis above “Las Güixas” cave (Villanúa village), approximately 100 $\mathrm{km}$ east of the Ostolo and Mendukilo caves. This show cave, located in the Central South Pyrenees (Fig. ), experiences a transitional Mediterranean–Oceanic climate (Cfb of the Köppen–Geiger climate classification) with a MAAT of 11 $\mathrm{°}\mathrm{C}$ and around 1100 $\mathrm{mm}$ of annual precipitation. During the winter the westerly winds and Atlantic fronts are responsible for most rainfall, while the rest of the year is mixed between Mediterranean and Atlantic fronts (Giménez et al., 2021), similar to the conditions in the area of Ostolo and Mendukilo caves. Stable isotope data in precipitation and air temperature on an event scale spanning 2 years are available from this station (2017–2019, Giménez et al., 2021). The weighted mean values of ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ and ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$ are $-$8.6 ‰ and $-$60.1 ‰, respectively, with seasonal variations reaching total amplitudes of 23 ‰ and 174 ‰, respectively (Giménez et al., 2021). The Local Meteoric Water Line (LMWL) is defined as ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$ $=$ 7.56 $\cdot$ ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ $+$ 4.33 ($n$ $=$ 210; $R^{\mathrm{2}}$ $=$ 0.97). The slope of the LMWL is close to that of the Global Meteoric Water Line (GMWL; Rozanski et al., 1993) and aligns well with the water line defined by the drip waters of Mendukilo and Ostolo caves (Fig. a). In general, the isotopic composition of rainfall correlates with air temperature for the 2-year period ($n$ $=$ 210; $R^{\mathrm{2}}$ $=$ 0.44, Fig. b) and shows moderate correlation with relative humidity and a weaker correlation with rainfall amount at the event scale when performing Spearman's correlation ($r_{\mathrm{s}}$; $n$ $=$ 180; between rainfall amount and ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ [${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$]: $r_{\mathrm{s}}$ $=$ $-$0.27 [$-$0.25]; between temperature and ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ [${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$]: $r_{\mathrm{s}}$ $=$ 0.70 [0.69]; between relative humidity at the rainfall site and ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ [${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$]: $r_{\mathrm{s}}$ $=$ $-$0.46 [$-$0.41]) (Giménez et al., 2021).

3 Methods

3.1 Sampling and petrography

Stalagmites OST1, OST2, and OST3 were retrieved from a gallery in Ostolo cave, where active speleothem deposition was not observed. Stalagmites MEN-2 and MEN-5 were retrieved from a gallery in Mendukilo cave, where active calcite precipitation was only observed at the original dripping point of MEN-5. See Bernal-Wormull et al. (2021, 2023) for more details on these caves. All stalagmites were cut longitudinally, and the central slab was polished. Small blocks were cut along the growth axis for the preparation of doubly polished thin sections (about 200 $\mathrm{\mu }\mathrm{m}$). FIs were studied in these thin sections using a Nikon Eclipse transmitted-light microscope.

3.2 FI stable isotopic composition

A total of 344 carbonate subsamples (including duplicates) were crushed and analyzed for ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ (287 subsamples of Mendukilo stalagmites and 69 of Ostolo samples). Between 0.3 and 2.5 $\mathrm{g}$ of calcite was used to ensure a sufficiently high water yield (0.1–1 $\mathrm{\mu }\mathrm{L}$). Speleothem fluid inclusion water isotopes were analyzed at the University of Innsbruck using continuous-flow analysis of water via high-temperature reduction on glassy carbon in a thermal combustion/elemental analyzer (TC/EA) unit. The ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ measurements were performed using a Delta V Advantage isotope ratio mass spectrometer. For details about the method of fluid inclusion preparation and isotopic analysis, see Dublyansky and Spötl (2009). ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ values are reported in per mil relative to VSMOW. The average long-term precision (1$\mathit{\sigma }$ uncertainty) of replicate measurements of an in-house calcite standard is $\pm$ 1.5 ‰ for ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ for water amounts between 0.1 and 1 $\mathrm{\mu }\mathrm{L}$.

${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ is regarded as a more robust proxy of paleotemperature than ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{FI}}$, as it is less influenced by non-climatic parameters, with no other sources of hydrogen affecting the water trapped in the calcite (Demény et al., 2016, 2021; Affolter et al., 2019). In addition, ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{FI}}$ values obtained with the Innsbruck FI setup can be inaccurate for samples of low water content ($<$ 0.1 $\mathrm{\mu }\mathrm{L}$; Dublyansky and Spötl, 2009). Therefore, we only used ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ values in this study.

4 Results

4.1 Petrography

The Ostolo and Mendukilo stalagmites consist of coarse crystalline calcite and are macroscopically homogenous without any sign of recrystallization. The MEN-2 and MEN-5 stalagmites exhibit a columnar fabric, lack growth hiatuses, and do not show macroscopically visible laminae (Bernal-Wormull et al., 2023). In contrast, the Ostolo stalagmites shows a more porous columnar microcrystalline fabric that transitions into an elongated-columnar type (Bernal-Wormull et al., 2021). Two hiatuses are present in OST3, marked by organic inclusions and micrite layers.

Figure 3

FI assemblages in the studied stalagmites. (a) Primary FI throughout the growth layers in stalagmite MEN-5. (b) Inter-crystalline FI in stalagmite MEN-2. (c) Intra-crystalline FI in stalagmite MEN-5. (d) Primary intra- and inter-crystalline FI in stalagmite OST2, more frequently found in porous areas or associated with elongated (Ce) and/or microcrystalline (Cm) fabrics than with a tightly packed columnar fabric (Cc). (e) FIs in stalagmite OST2 are mostly intra-crystalline and do not necessarily align with the growth layers. (f) Pyriform and rounded intra-crystalline small FI in stalagmite OST2. Black arrows indicate the speleothem growth direction.

[Figure omitted. See PDF]

Primary FIs were observed in all stalagmite samples (Fig. ). The Mendukilo samples contain considerably more FIs compared to those from Ostolo, mainly concentrated along growth layers (Fig. a). In the Mendukilo stalagmites, primary inter-crystalline (10–30 $\mathrm{\mu }\mathrm{m}$; Fig. b) and intra-crystalline (10 to $>$ 100 $\mathrm{\mu }\mathrm{m}$; Fig. c) FIs are discernible. These intra-crystalline primary FIs are elongated and rounded or pyriform in shape (rounded at the base with a spike extending in the speleothem growth direction; Fig. c; Lopez-Elorza et al., 2021). In Ostolo, FIs are less prominent and are mostly intra-crystalline, located along or around white porous laminae and within the more elongated columnar or microcrystalline fabrics (Fig. d and e). The intra-crystalline FIs in Ostolo samples are, on average, smaller than those in the Mendukilo stalagmites (10–40 $\mathrm{\mu }\mathrm{m}$) and predominantly exhibit pyriform or rounded shapes (Fig. f). Petrographic observations confirm that the FIs in these samples are primary, well preserved, and suitable for their stable isotopic analysis.

Figure 4

(a) ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$ of FI water (${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$) and ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ of calcite (${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$) of Mendukilo (MEN-5 in orange and MEN-2 in black) and (b) Ostolo stalagmites (OST1 in blue, OST2 in red, and OST3 in green). ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ values are corrected for the ice-volume effect (Bintanja et al., 2005), with vertical error bars representing isotope measurements errors and 1$\mathit{\sigma }$ from repeated measurements. The yellow dashed line in the upper graphs of each panel indicates the annual mean ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$ value in drip water for each cave. Modeled $\mathrm{U}/\mathrm{Th}$ ages with 2$\mathit{\sigma }$ error bars for stalagmites from each cave are shown at the bottom. Heinrich event 1 (HE1) recorded in the Ostolo isotope record (Bernal-Wormull et al., 2021) is highlighted by a light-blue bar.

[Figure omitted. See PDF]

4.2 Last deglaciation and the Holocene ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ speleothem record

The chronology of the Ostolo stalagmites spans the last deglaciation between 16.5 and 11.7 $\mathrm{kyr}\mathrm{BP}$ with high precision due to their very high ${}^{\mathrm{238}}\mathrm{U}$ concentrations (10–80 $\mathrm{ppm}$). The carbonate ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ (${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$) profiles show consistency among the three stalagmites (Fig. ). OST1 and OST2 have more negative values ($-$5 ‰ to $-$8.9 ‰) during GS-1 and GS-2.1a and less negative values (up to $-$3.4 ‰) during GI-1 and the onset of the Holocene. OST3 did not grow during the intervals characterized by the most negative ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ values recorded by the other two stalagmites (Fig. ). On the other hand, the MEN stalagmites, despite having lower ${}^{\mathrm{238}}\mathrm{U}$ concentrations (100–350 $\mathrm{ppb}$), also have lower detrital ${}^{\mathrm{232}}\mathrm{Th}$ contents, enabling robust age models for both stalagmites. These models cover various intervals of the Holocene and GS-1 with good overlap (Fig. ), specifically as follows: (i) MEN-2 grew between 12.8 and 6.3 $\mathrm{kyr}\mathrm{BP}$, with ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ values that remain stable during GS-1, followed by an abrupt increase, reaching the highest values of the entire record at the GS-1/Holocene transition (from $-$5.2 ‰ in GS-1 to $-$4.3 ‰ at 11.6 $\mathrm{kyr}\mathrm{BP}$). (ii) MEN-5 spans the last 8.8 $\mathrm{kyr}$ and presents prominent negative values during certain short events (e.g., 8.2 $\mathrm{kyr}\mathrm{BP}$ with a value of $-$6.3 ‰, replicated by MEN-2), which are synchronous, within age uncertainties, with abrupt changes in the isotopic composition of North Atlantic surface waters (Kleiven et al., 2008; Carlson et al., 2008). More details on the chronology and isotopic data of these speleothems are provided by Bernal-Wormull et al. (2021, 2023).

4.3 FI isotopes

OST samples are characterized by variable water content, with replicates yielding a mean standard deviation of $\pm$ 2.7 ‰ for ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$. We assigned this value to individual measurements as an overall uncertainty estimate. Not all OST samples could be duplicated due to sometimes low water amounts and petrographically complex FI assembles in some samples (Fig. d and e), which restricted subsampling of some individual growth layers. All MEN measurements were duplicated, triplicated, or even quadruplicated. The ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$ values of sub-samples of MEN-2 and MEN-5 (ranging between $-$34 ‰ and $-$61 ‰) with water contents of 0.1 to 1 $\mathrm{\mu }\mathrm{L}$ replicated within 2.7 ‰.

${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ values for the Holocene and GI-1 (Tables  and ) are comparable to cave drip waters at Mendukilo and Ostolo caves (Fig. ). In contrast, values are more negative during GS-1 and GS-2.1a (Fig. ). GS-2.1a is represented by 8 OST subsamples with a mean ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ value of $-$50 ‰. One of these values, dated to 15.80 $\pm$ 0.05 $\mathrm{kyr}\mathrm{BP}$, is even more depleted ($-$59.7 $\pm$ 2.6 ‰). Values become less negative rapidly at 14.57 $\pm$ 0.05 $\mathrm{kyr}\mathrm{BP}$ (Fig. ; mean during GI-1: $-$37 ‰). This trend is interrupted in the three OST stalagmites at 14.13 $\pm$ 0.09 $\mathrm{kyr}\mathrm{BP}$, leading to more negative values (between $-$43 ‰ and $-$51 ‰). During GS-1, the ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ values decrease again (Fig. ), averaging $-$48 ‰ before showing a rapid increase at the onset of the Holocene ($-$37 ‰). The MEN-2 record also shows a mean of $-$48 ‰ during GS-1, though the transition to the Holocene is more gradual. Between 8.7 and 6.3 $\mathrm{kyr}\mathrm{BP}$, MEN-2 and MEN-5 ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ values show excellent correlation (Fig. ). There is no significant variation between the Greenlandian ($-$42 ‰), Northgrippian ($-$43 ‰), and Meghalayan ($-$42 ‰). Despite these relatively stable ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ values throughout the Holocene substages, a short negative shift is identified in two samples dated at 8.28 $\pm$ 0.08 ($-$53.5 $\pm$ 4.0 ‰) and 8.31 $\pm$ 0.06 ($-$56.3 $\pm$ 6.5 ‰) $\mathrm{kyr}\mathrm{BP}$.

5 Discussion

5.1 Interpretation of the ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ signal

At isotopic equilibrium, the ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ value is related to the ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ of the drip water and the cave temperature through its control on the equilibrium isotope fractionation between water and calcite (Lachniet, 2009). Additionally, variations in stalagmite ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ records may reflect changes in the ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ of surface ocean waters from the moisture source area as well as changes in atmospheric processes, which control the fractionation of oxygen isotopes in route to the site where rainfall occurs (Lachniet, 2009; McDermott, 2004). In the Ostolo stalagmites, the ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ signal is coherent with air temperature changes throughout the deglaciation period (Bernal-Wormull et al., 2021). The overall ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ pattern observed in these stalagmites is similar to that of speleothems from the Central Pyrenees (Bartolomé et al., 2015; Cheng et al., 2020) and the Alps (Li et al., 2020; Luetscher et al., 2015), which also predominantly receive Atlantic-derived moisture and where ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ primarily reflects atmospheric temperature. Superimposed on the temperature effect are changes in the isotopic composition of seawater (these changes are discussed in more detail in Sect. ), which may account for the negative excursion in the Ostolo ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ record during Heinrich event 1 (HE1) at 16.2–16.0 $\mathrm{kyr}\mathrm{BP}$, with values reaching as low as $-$8.9 ‰ (Bernal-Wormull et al., 2021; Fig. ).

Monitoring of the Mendukilo cave suggests that the MEN ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ record captured an annual signal, which is primarily influenced by rainfall amount (Bernal-Wormull et al., 2023). Conversely, the MEN ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ record captures a temperature signal that is obscured by the influence of rainfall amount, since temperature and humidity changes may have competing effects on the ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ signal (Bernal-Wormull et al., 2023). Additionally, during the earlier part of the record (13–8 $\mathrm{kyr}\mathrm{BP}$), changes in the oceanic isotopic composition associated with meltwater input (Català et al., 2019; Skinner and Shackleton, 2006) further affect the signal. A prominent feature of the MEN-2 and MEN-5 ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ records is a $-$0.7 ‰ anomaly (relative to the Holocene mean of $-$5.4 ‰) observed at 8.11 and 7.00 $\mathrm{kyr}\mathrm{BP}$ (Fig. ). These two events of anomalously low ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ values likely reflect rapid, short-lived decreases in temperature and in the ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ of the surface ocean water, rather than increased rainfall, as proposed in previous studies (e.g., Domínguez-Villar et al., 2009; García-Escárzaga et al., 2022; LeGrande and Schmidt, 2008; Matero et al., 2017; Stoll et al., 2022).

5.2 Isotope–temperature conversion

The composite paleotemperature records of the Ostolo and Mendukilo speleothems are based on 356 FI samples (and replicates), applying a regional water isotope–temperature relationship derived from monitoring data (isotopic data of drip water and outside temperature) of both caves (Bernal-Wormull et al., 2021, 2023) and the relationship between rainfall ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$ (${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{r}}$) and modern air temperature above Las Güixas cave. The latter provides a relationship between air temperature and the stable isotopic composition of rain (${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{r}}$ and ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{r}}$) observed from July 2017 to June 2019 ($n$ $=$ 210). The observed correlation between ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{r}}$, ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{r}}$, and air temperature is verified at biannual scale above Las Güixas cave, with significant correlation between MAAT and the values of ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{r}}$ and ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{r}}$, based on a multiple regression model using a univariate Spearman's correlation of different climatic parameters (Giménez et al., 2021). For example, the correlation between ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{r}}$ and air temperature at the time of precipitation (same data series) results in a $r_{\mathrm{s}}$ $=$ 0.69 ($p$ $\ll$ 0.01) by applying this multiple regression model (Giménez et al., 2021).

${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ and ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$ values of seawater vary on glacial–interglacial timescales due to the ice-volume effect: when surface water evaporates from the ocean, lighter stable isotopes are preferentially removed into the vapor phase, leading to increased ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ and ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$ values in the ocean water as more fresh water is stored as ice on continents (Lachniet, 2009). ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ values were corrected for the ice-volume effect during the deglaciation period covered by the MEN and OST speleothems. This correction used a gradient derived for ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ (Bintanja et al., 2005) converted to ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$ using a factor of 8. Paleotemperatures were then estimated using a linear ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}/\mathrm{T}$ transfer function anchored to the MAAT at both cave sites and the isotopic composition of drip water (${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{d}}$; Ostolo ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{d}}$ $=$ $-$37.8 ‰; Mendukilo ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{d}}$ $=$ $-$45.3 ‰), with corrections for the elevation of the Villanúa monitoring station (950 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$). The modern ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$ values were adjusted for the elevation difference between the rainfall sampling station and the studied caves, assuming a lapse rate of 0.2 ‰ per 100 $\mathrm{m}$ for ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$, i.e., 1.6 $\mathrm{‰}\mathrm{per}\mathrm{100}\mathrm{m}$ for ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{p}}$ (Poage, 2001). Uncertainties of the OM-FIT record reflect isotope measurement errors and the selected paleotemperature transfer function. The uncertainties associated with ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$, ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{d}}$, ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}/\mathrm{T}$, and MAAT at Las Güixas tourist cave, as well as the slope of the LMWL, were propagated through the calculation steps (for more details on the calculation of this error propagation, see Table  in the Appendix). Due to a lack of constraints on past seasonal changes in precipitation and effective infiltration, we assume constant annual infiltration over time.

Figure 5

(a) ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ records from Mendukilo and Ostolo stalagmites, compared to (b) the OM-FIT paleotemperature reconstruction (bottom). Heinrich event 1 (HE1) is highlighted by a light-blue bar. The MAAT outside the two caves is shown by dashed horizontal lines.

[Figure omitted. See PDF]

5.3 OM-FIT: paleothermometric record derived from FI stable isotope data

Our ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ values (Tables  and ; Fig. ) provide a robust record because (i) part of the record is well replicated by samples from two caves from different climatic settings (e.g., during the GS-1), (ii) stalagmites from the same cave are replicated (within their respective uncertainties), and (iii) a large proportion of the samples have multiple replications. We investigated the temperature dependence of the hydrogen (and oxygen) isotope composition of precipitation water in the study region, examining the modern-day ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}/\mathrm{T}$ and ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}/\mathrm{T}$ gradients. This relationship, which may change over time, was examined by Rozanski et al. (1992) for central Europe and applied by Affolter et al. (2019) to a 14 $\mathrm{kyr}$ record from Milandre cave (Switzerland). It was similarly applied to Last Interglacial records from Alpine caves (Honiat et al., 2023; Wilcox et al., 2020). The relationship between mean annual ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{r}}$ and MAAT (${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$/T) is 0.55 $\pm$ 0.05 $\mathrm{‰}\mathrm{°}{\mathrm{C}}^{-\mathrm{1}}$ (in the case of ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{r}}$, the ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}/\mathrm{T}$ is 4.00 $\pm$ 0.31 $\mathrm{‰}\mathrm{°}{\mathrm{C}}^{-\mathrm{1}}$) for the “Las Güixas” tourist cave in Villanúa (Giménez et al., 2021), which is consistent with the average European ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$/T gradient of 0.59 $\pm$ 0.08 $\mathrm{‰}\mathrm{°}{\mathrm{C}}^{-\mathrm{1}}$ (Rozanski et al., 1992). The OST and MEN FI isotope data overlap chronologically for the GS-1, allowing for their combination into a single temperature transfer function (OM-FIT) covering the last 16.7 $\mathrm{kyr}\mathrm{BP}$ (Fig. ). The OM-FIT is calculated using the corrected ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ values, ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{d}}$, MAAT ($T_{\text{modern}}$), and the modern-day ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}/\mathrm{T}$ gradient derived from the LMWL of rainfall isotopes:

1$$T_{\text{OM-FIT}}=T_{\text{modern}}-{\displaystyle \frac{{\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{d}}-{\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}\left(\mathrm{corrected}\right)}}{{\mathit{\delta }}^{\mathrm{2}}\mathrm{H}/T_{\text{gradient}}}}$$

Equation () expresses the fact that the paleotemperatures were lower than present-day temperatures; hence the calculated temperature difference should be subtracted from $T_{\text{modern}}$. The temperature reconstruction with Eq. () is based on the mean relationship of 4.00 $\mathrm{‰}\mathrm{°}{\mathrm{C}}^{-\mathrm{1}}$ (${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}/T_{\text{gradient}}$; Fig. c) for both caves. In this case, we assessed temperature variations using speleothems from two distinct caves. As a result, it is crucial to consider the values of ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{d}}$ and $T_{\text{modern}}$ (see Sect. ) that correspond to the ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ values of the speleothems from the same cave when using Eq. (). The final calculated uncertainty in the paleotemperature ranges from 1.6 to 2.6 $\mathrm{°}\mathrm{C}$ (Table ).

Figure 6

Paleotemperature reconstructions over the last 16.5 $\mathrm{kyr}\mathrm{BP}$, spanning Greenland to SW Europe, along with speleothem ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ records from the Iberian Peninsula. (a) NGRIP ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ (gray solid line; Rasmussen et al., 2014) and Greenland temperature reconstruction (black solid line; Kindler et al., 2014). (b) Milandre cave FI temperature record (MC-FIT) from NW Switzerland (Affolter et al., 2019). (c) July temperature inferred from chironomids at Basa de la Mora Lake (Tarrats et al., 2018). (d) Stacked and spliced chironomid-inferred July temperature record from SW Europe (Heiri et al., 2014b). (e) Ostolo and Mendukilo FI temperature record (OM-FIT; yellow star; this study). (f) Regional variations in freshwater levels identified in the NISA composite splice speleothem ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ record, compared to the benthic ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ record from the Iberian margin (${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{NISA}\text{-}\mathrm{BC}}$; light-blue line; Stoll et al., 2022), the sea-surface salinity (SSS) reconstruction derived from planktonic foraminifera (pink dashed line; Eynaud et al., 2012) and the ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ from Globigerina bulloides of marine core MD99-2334K (green dashed line; Skinner and Shackleton, 2006). (g) ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$c records from Mendukilo and Ostolo (Bernal-Wormull et al., 2021, 2023). (h) LV5 ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ record from Kaite cave (northern Iberia; Domínguez-Villar et al., 2017). (i) GAR-01 ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ record from La Garma cave (northern Iberia; (Baldini et al., 2019). (j) SIR-14 ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ record from El Soplao cave (northern Iberia; Kilhavn et al., 2022). Key abrupt climate events (Heinrich 1 (HE1), 9.3 and the 8.2 $\mathrm{kyr}$ events) and Greenland stadials (GS-1 and GS-2.1a) are highlighted by a light-blue bar. The gray envelope around the solid lines in (b)–(e) shows the uncertainties.

[Figure omitted. See PDF]

5.4 Temperature regime of northern Spain based on OM-FIT

5.4.1 Last deglaciation

The Ostolo cave ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ and ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ records (Fig. ) and the OM-FIT (Fig. ) show clear evidence of rapid temperature changes during GS-2.1a, GI-1, GS-1, and the onset of the Holocene. The timing and amplitude of these changes are in good agreement with other European oxygen isotope records from lake sediments (Van Raden et al., 2013; Von Grafenstein et al., 1999) and speleothems (Affolter et al., 2019; Cheng et al., 2020; Li et al., 2020; Luetscher et al., 2015). The strong similarity between these records and NGRIP ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ (Rasmussen et al., 2014) and temperature reconstructions (Kindler et al., 2014) (Fig. ) supports the idea of a common North Atlantic climate forcing during the last deglaciation on millennial to centennial timescales.

The OM-FIT record suggests that regional MAAT during GS-2.1a was slightly lower than during GS-1, characterized by a negative excursion at 15.80 $\pm$ 0.07 $\mathrm{kyr}\mathrm{BP}$ and a temperature decrease of approximately 2.5 $\mathrm{°}\mathrm{C}$ relative to the GS-2.1a average (Fig. ). This negative temperature excursion is supported by a single point in the ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ record from OST2 (Figs.  and ). The same ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ values were obtained in further analyses of the same stalagmite and OST1 but unfortunately with water contents $<$ 0.1 $\mathrm{\mu }\mathrm{L}$. Nevertheless, this anomaly coincides with the most negative ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ values in stalagmites OST1 and OST2 (Figs.  and ), which corroborate that during this short period of time in GS-2.1a there was a common catalyst in the variability of the isotopic values recorded in the speleothems of Ostolo cave (Bernal-Wormull et al., 2021). Indeed, this OM-FIT anomaly corresponds with the final phase of HE1, related to massive iceberg discharges from the Laurentide ice sheet, which collapsed around 16.2 $\pm$ 0.3 $\mathrm{kyr}\mathrm{BP}$ (Landais et al., 2018). Regionally, a significant glacier advance occurred at that time in the Pyrenees and other Iberian mountains (García-Ruiz et al., 2023), and speleothems from Meravelles cave (NE Iberia) record a notable ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ anomaly between 16.2 and 15.9 $\mathrm{kyr}\mathrm{BP}$ (Pérez-Mejías et al., 2021). Furthermore, this change during HE1 has also been observed in ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ records from the Atlantic Ocean (core MD99–2334; Skinner and Shackleton, 2006), sea-surface salinity (SSS) reconstructions in the Bay of Biscay (core MD95–2002; Eynaud et al., 2012) and already detected in the speleothem records from low altitude caves in the Cantabric range from N Iberia (Stoll et al., 2022) (Fig. ). The negative excursion in the Ostolo ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ and the OM-FIT precisely coincides with the maximum discharge of ice melting accordingly with these records sensitive to North Atlantic freshening and thus appears to reflect changes in the isotopic composition of the moisture source (Fig. ). This observation confirms that the OM-FIT captured not only temperature history on millennial scales but also abrupt climate events on a centennial scale.

A rapid temperature increase of 6.7 $\pm$ 2.8 $\mathrm{°}\mathrm{C}$ occurred at the onset of GI-1 (Fig. ). This sudden temperature shift, occurring between 14.78 $\pm$ 0.18 and 14.56 $\pm$ 0.05 $\mathrm{kyr}\mathrm{BP}$ in the OM-FIT record (Fig. ), aligns with data from Greenland ice cores marking the beginning of GI-1 (14.64 $\pm$ 0.19 $\mathrm{kyr}\mathrm{BP}$; Rasmussen et al., 2014). This indicates a nearly synchronous response to a significant climate warming due to a stronger Atlantic Meridional Overturning Circulation (AMOC), which enhanced the northward transport of ocean heat, thereby raising temperatures in the Northern Hemisphere (Clark et al., 2012; Muschitiello et al., 2019). This increase in the OM-FIT record coincides with an important glacier retreat in the Iberian mountains (García-Ruiz et al., 2023), an increase in chironomid-inferred July air temperatures (from ca. 11 to ca. 16 $\mathrm{°}\mathrm{C}$) from the West-Central Pyrenees, Millet et al., 2012), and an increase in MAAT (from ca. 12.2 to ca. 18.6 $\mathrm{°}\mathrm{C}$) recorded by branched glycerol dialkyl glycerol tetraethers in the Padul palaeolake record (Sierra Nevada, southern Iberian Peninsula; Rodrigo-Gámiz et al., 2022). The onset of GI-1 in the OM-FIT was recorded by ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ data from the OST1 and OST3 stalagmites. The amplitude of this abrupt warming is in agreement with other European temperature records, such as estimates based on ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ data from Alpine speleothems (Li et al., 2020; Luetscher et al., 2015). Von Grafenstein et al. (2013) used a combination of ostracod, mollusc, and charophyte data to estimate a rise of about 6 $\mathrm{°}\mathrm{C}$ in MAAT for this transition at the Gerzensee lake site. The Ammersee record, using a coefficient derived from a study of northern Switzerland stalagmites (0.48 $\mathrm{‰}\mathrm{°}{\mathrm{C}}^{-\mathrm{1}}$, Affolter et al., 2019), estimated a warming of about 5.5 $\mathrm{°}\mathrm{C}$ (4.1–8.4 $\mathrm{°}\mathrm{C}$) (Li et al., 2020) for this transition.

During GI-1, the ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ record is marked by higher ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$ values and similar temperatures in the OM-FIT record compared to the onset of the Holocene (Fig. ). As observed in the OST ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ record, ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ values follow a negative trend towards the end of GI-1. However, pollen and speleothem records of central and southern Europe suggest that temperate forests in southern and central Europe gradually expanded (Fletcher et al., 2010; Moreno et al., 2010), which contrasts with the temperature pattern in Greenland (Fig. ) and OST ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ values (Fig. ) that shows the warmest peak at the beginning of GI-1. This temperature rise is evident in most North Atlantic, European, and Greenland records and the OM-FIT (Fig. ), aligning with a more active AMOC at the start of the GI-1 (Clark et al., 2012; Heiri et al., 2014b; Kindler et al., 2014; Muschitiello et al., 2019; Rasmussen et al., 2014). Therefore, the contrasting patterns observed in the pollen and speleothem records compared to Greenland and OST ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ values could have been due to variations in moisture availability rather than temperature. Within this interstadial, a significant inflection point occurs, with a negative anomaly at 14.10 $\pm$ 0.09 $\mathrm{kyr}\mathrm{BP}$ in the OM-FIT record. This suggests that the OM-FIT minimum during GI-1, also registered at 14.10 $\pm$ 0.03 $\mathrm{kyr}\mathrm{BP}$ in the OST ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ record and equivalent to GI-1d in NGRIP (Rasmussen et al., 2014), involved the most pronounced cooling of GI-1 (between 3.0 and 4.7 $\pm$ 2.5 $\mathrm{°}\mathrm{C}$ in the OM-FIT record), occurring just after the GI-1e warm phase (Fig. ). This cooling event is contemporaneous with glacier expansions in the Pyrenees (García-Ruiz et al., 2023) and a centennial-scale cooling at Ech paleolake (Millet et al., 2012), Lake Estanya (Vegas-Vilarrúbia et al., 2013) and in the Portalet sedimentary sequence (González-Sampériz et al., 2006). Apparently, this relatively small decrease in temperature during GI-1d, as quantified by the OM-FIT record and chironomid-inferred July air temperatures (Millet et al., 2012) in this region, resulted in an important vegetation response (González-Sampériz et al., 2017), characterized by a decrease in juniper and an expansion of steppe herbs during this cold and dry event.

Between 13.0 and 12.5 $\mathrm{kyr}$ B.P., the ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ decrease (Fig. ) records a cooling of 6.1 $\pm$ 2.8 $\mathrm{°}\mathrm{C}$ in the OM-FIT record (Fig. ), marking the initial part of GS-1 Rasmussen et al. 2014). The OM-FIT is the first record of quantitative temperature changes captured by stalagmites from southwestern Europe and provides insights into the timing and intensity of temperature change during the GS-1 in this region due to rapid cooling associated with enhanced freshwater discharge into the North Atlantic or the Arctic Ocean, which led to an AMOC slowdown and triggered extreme cooling in the North Atlantic realm (Clark et al., 2012; Eynaud et al., 2012). Similar cooling magnitudes to those recorded by the OM-FIT were reported from a stalagmite in Seso cave, Central Pyrenees, based on ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ data (Bartolomé et al., 2015). On the other hand, this change appears slightly larger compared to cooling registered by summer air temperature records of the GI-1/GS-1 transition, such as those from lake sediments in NW Iberia (2–3 $\mathrm{°}\mathrm{C}$; Muñoz Sobrino et al., 2013) and the Central Pyrenees (1.5–2 $\mathrm{°}\mathrm{C}$; Millet et al., 2012). This important change in the OM-FIT record also agrees in magnitude with a rapid cooling recorded by (i) speleothems from the Alps (around 4–5 $\mathrm{°}\mathrm{C}$; Li et al., 2020) and the Jura Mountains (4.3 $\pm$ 0.8 $\mathrm{°}\mathrm{C}$; Affolter et al., 2019) and (ii) a drop in sea-surface temperatures of around 4 $\mathrm{°}\mathrm{C}$ off the Iberian coast at 12.9 $\mathrm{kyr}\mathrm{BP}$ (Martrat et al., 2014; Rodrigues et al., 2010).

The end of the GS-1 cold phase and the onset of the Holocene are marked by a rapid warming in the OM-FIT record of about $\sim$5 $\mathrm{°}\mathrm{C}$ (Fig. ), peaking at 11.67 $\pm$ 0.02 $\mathrm{kyr}\mathrm{BP}$. This warming aligns (within uncertainties) with the record from Greenland ice cores (14.65 $\pm$ 0.10 $\mathrm{kyr}\mathrm{BP}$; Rasmussen et al., 2014). This further reinforces OM-FIT's capability to detect temperature changes in the Northern Hemisphere associated with a stronger AMOC at the Holocene onset (Muschitiello et al., 2019). The variability of MEN ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ data during the GI-1/GS-1 and GS-1/Holocene onset transitions is less pronounced compared to OST ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$. This observation may be due to the proximity of Mendukilo cave to the Atlantic coast, with temperature and humidity changes having competing effects on ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$, as already reported in other speleothem records from this region (e.g., Baldini et al., 2019). In contrast, the ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ of speleothems from Pyrenean caves is predominantly controlled by temperature (Bartolomé et al., 2015; Bernal-Wormull et al., 2021; Cheng et al., 2020), resulting in a more subtle temperature signal (less amplitude of isotopic change) in the case of MEN ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ compared to the OST ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ record during GS-1, a cold and dry period (Fletcher et al., 2010). Nevertheless, the MEN ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ records capture important changes during the GI-1/GS-1 and GS-1/Holocene transitions and correlate quite well with the ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ data from OST (Fig. ).

5.4.2 Holocene

As mentioned above, the Holocene section of the OM-FIT record (Fig. ) is based on ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ values of the MEN stalagmites (Fig. ). This record not only captures variability in ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ composition influenced by temperature but also reflects past hydroclimatic conditions (Bernal-Wormull et al., 2023). This observation introduces a limitation in reconstructing periods of relatively stable temperature, such as the Holocene, which is represented by centennial-scale OM-FIT temperature variability that reaches up to 2 $\mathrm{°}\mathrm{C}$ in certain intervals. However, these variations are close to the uncertainty range of the OM-FIT record ($\pm$ 1.6 to $\pm$ 2.6 $\mathrm{°}\mathrm{C}$ for the Holocene). Therefore, these reconstructed quantitative temperature data for the Holocene must be viewed with caution. On millennial scales, the OM-FIT record indicates the highest temperatures at the onset of the Holocene (until $\sim$10 kyr BP), albeit with high variability. Between 11.7 and 10.9 $\mathrm{kyr}\mathrm{BP}$, OM-FIT records a mean temperature of 13.1 $\pm$ 1.9 $\mathrm{°}\mathrm{C}$, just 0.6 $\mathrm{°}\mathrm{C}$ higher than the mean temperature value for the GI-1 of 12.5 $\pm$ 1.9 $\mathrm{°}\mathrm{C}$ (Fig. ). This early rapid warming is also recorded by the hydroclimate-sensitive isotopic signal of the SIR-1 stalagmite from NW Iberia (Rossi et al., 2018) at the onset of the Greenlandian period. This observation underscores the value of obtaining a temperature-sensitive FI record in regions where the carbonate isotopic signal of speleothems is also influenced by the amount effect, such as the MEN ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ record.

The OM-FIT does not capture a clear cooling trend after the Holocene Thermal Maximum (HTM) compared to the ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ record from Greenland ice cores (Rasmussen et al., 2014) and the Milandre cave fluid inclusion temperature record (MC-FIT) from central Europe (Affolter et al., 2019) (Fig. ), instead suggesting stable temperatures. This Neoglacial cooling, widespread across the Northern Hemisphere, is well documented throughout Europe (e.g., Larocque-Tobler et al., 2010; Ilyashuk et al., 2011) and Iberia (Català et al., 2019; García-Ruiz et al., 2020; Leunda et al., 2019; Sancho et al., 2018) during the Northgrippian period. The absence of this cooling in the OM-FIT record is likely due to large centennial variability and large uncertainties. The temperature trends in MEN ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ and OM-FIT differ from those captured by chironomids in the Central Pyrenees (Tarrats et al., 2018), which indicate a millennial-scale cooling during the Middle Holocene compared to the HTM and the Late Holocene (Fig. ). This observation highlights the differences between temperature records derived from speleothems (OM-FIT, without seasonal bias) and chironomids (recording summer air temperature), as in the case for GS-1.

Despite the limited precision of OM-FIT, it can identify abrupt centennial events, some of which are also evident in the ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ values of MEN-2 and MEN-5 (Fig. ). For example, one of the lowest Holocene OM-FIT temperatures (9.6 $\mathrm{°}\mathrm{C}$) occurred at 11.50 $\pm$ 0.08 $\mathrm{kyr}\mathrm{BP}$ (mean temperature at the onset of the Holocene, 12.3 $\pm$ 1.8 $\mathrm{°}\mathrm{C}$), corresponding within age uncertainties to the Preboreal Oscillation (11.4 $\mathrm{kyr}$) recorded in Greenland ice cores (11.27 $\pm$ 0.03 $\mathrm{kyr}\mathrm{BP}$, based on the new ice-core chronology – Seierstad et al., 2014) and by MC-FIT in Switzerland (11.37 $\pm$ 0.15 $\mathrm{kyr}\mathrm{BP}$ – Affolter et al., 2019) (Fig. ). Another example is the 9.2 $\mathrm{kyr}$ event, documented across the Northern Hemisphere (e.g., Masson-Delmotte et al., 2005; Genty et al., 2006; Rasmussen et al., 2007; Fleitmann et al., 2008) and supported by terrestrial (Baldini et al., 2019; Carrión, 2002; Iriarte-Chiapusso, 2016; Mesa-Fernández et al., 2018; Vegas et al., 2010) and marine records from Spain (Nebout et al., 2009). This event is captured by a ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ value of $-$50.6 ‰ in MEN-2 and an OM-FIT temperature of 10.4 $\pm$ 2.1 $\mathrm{°}\mathrm{C}$ at 9.29 $\pm$ 0.08 $\mathrm{kyr}\mathrm{BP}$ (Fig. ). However, it is absent from the ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ record of MEN-2, and previous research suggests that the climate in northern Spain was likely considerably warmer and wetter $\sim$ 9 $\mathrm{ka}$ BP (Baldini et al., 2019; Morellón et al., 2018; Tarrats et al., 2018). This observation supports that the less variable ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ signal in Mendukilo cave is influenced by short-lived decreases in ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{sw}}$ and changes in humidity, while the OM-FIT shows the quantitative change in temperature during the 9.2 $\mathrm{kyr}$ event (Fig. ).

Catastrophic meltwater discharge during the “8.2 $\mathrm{kyr}$ event” from glacial lake Agassiz lowered the isotope composition of North Atlantic surface water by 0.4 ‰ (Carlson et al., 2008; Kleiven et al., 2008) and led to a widespread cooling across the circum-North Atlantic. The isotopic signal of this meltwater event was transported by the westerlies and left an imprint in the isotopic composition of precipitation in Iberia (Bernal-Wormull et al., 2023; LeGrande and Schmidt, 2008). The 8.2 $\mathrm{kyr}$ event overlapped a multi-centennial cool period from 8.29 to 8.10 $\pm$ 0.04 $\mathrm{kyr}\mathrm{BP}$ recorded by MEN ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$, characterized by an abrupt drop in temperature of about $\sim$ 3.0 $\mathrm{°}\mathrm{C}$ between 8.31 $\pm$ 0.06 and 8.28 $\pm$ 0.08 $\mathrm{kyr}\mathrm{BP}$ in the OM-FIT record (Fig. ). This cooling within an interglacial coincided with significant vegetation changes in the Iberian Peninsula (Allen et al., 1996; Carrión and Van Geel, 1999; González-Sampériz et al., 2006). This could be important for assessing future climate conditions in this region if changes in large parts of the climate system (climate tipping elements; Armstrong McKay et al., 2022) intensify beyond a warming threshold.

The cooling amplitude during the 8.2 $\mathrm{kyr}$ event recorded by OM-FIT appears at first glance more pronounced than in other Northern Hemisphere temperature and precipitation records, with proxy evidence across Europe indicating a cooling by $\sim$ 1–1.7 $\mathrm{°}\mathrm{C}$ during this event (Baldini et al., 2019; Davis et al., 2003; Morrill et al., 2013), but unfortunately the uncertainties involved in this records and the OM-FIT do not allow a valid comparison to be made on this occasion. Other terrestrial records in southwestern Europe offer important insights into the paleoenvironment during this event (e.g., Fletcher et al., 2013; González-Sampériz et al., 2017; Morellón et al., 2018; Zielhofer et al., 2019). Some records offer opposing perspectives on paleohumidity conditions, due to the study region's exposure to both Mediterranean and North Atlantic air masses (Moreno et al., 2017, 2021). However, many of these terrestrial records reflect larger climate trends and often lack the resolution needed to accurately capture the regional response to the 8.2 $\mathrm{kyr}$ event. It is therefore likely that these long-term changes are more influenced by summer insolation than by the meltwater discharge into the North Atlantic Ocean, as suggested by Kilhavn et al. (2022). Thus, other stalagmite records from the region (Baldini et al., 2019; Domínguez-Villar et al., 2009; Kilhavn et al., 2022) (Fig. ) and the combination of the carbonate isotopic composition (${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ and ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{c}}$; Bernal-Wormull et al., 2023) and the FI record from Mendukilo stalagmites offer a better understanding of the regional response during this colder-than-average Holocene period, which was characterized by increased humidity and changes in moisture source composition (Domínguez-Villar et al., 2009; Kilhavn et al., 2022; Bernal-Wormull et al., 2023).

6 Conclusions

The Ostolo and Mendukilo speleothems provide a replicated and precisely dated record of paleotemperature in NE Iberia for the past 16.5 $\mathrm{kyr}\mathrm{BP}$. The OM-FIT record contributes novel, non-biogenic evidence of rapid temperature transitions and emphasizes the significant connections between abrupt air temperature shifts in this southwestern European terrestrial record and changes in deep-water ocean circulation and meltwater input during the last deglaciation and the Greenlandian period (on millennial and/or centennial time scales). Our findings indicate temperatures for GS-2.1a up to 6.7 $\pm$ 2.8 $\mathrm{°}\mathrm{C}$ lower than those for GI-1 and present-day conditions and constrain the regional response of HE-1 between 16.2 and 15.8 $\mathrm{kyr}\mathrm{BP}$ with meltwater contributions from the North Atlantic. The sharp rise in temperatures during the GS-2.1a/GI-1 transition was quantitatively comparable to other records from SW Europe. Temperatures during GI-1 were equivalent to those of the Holocene, with a minimum observed at 14.10 $\pm$ 0.09 $\mathrm{kyr}\mathrm{BP}$ during GI-1d. All these shifts and noticeable cooling events observed in the OM-FIT were not only synchronous with climate variations associated with subarctic freshening documented by other records of Europe but also synchronous with glacier advances in the Pyrenees during GS-2.1a and GI-1d. The rapid temperature changes at early GS-1 and the onset of the Holocene recorded by OM-FIT are consistent with to those reported from other parts of Europe and also confirm quantitatively equally high temperatures for both GI-1 and the onset of the Holocene. Neither ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ nor OM-FIT reveal significant millennial-scale changes during the Holocene. The 8.2 $\mathrm{kyr}$ event is recorded between 8.29 and 8.10 $\pm$ 0.04 $\mathrm{kyr}$ in the ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}$ record, centered between 8.31 $\pm$ 0.06 and 8.28 $\pm$ 0.08 $\mathrm{kyr}\mathrm{BP}$ in the OM-FIT record, synchronous with Greenland ice-core data and well-dated records from central and SW Europe.

Future quantitative temperature studies using well-dated speleothems from SW Europe are necessary to corroborate the chronology and better delineate the drastic centennial-scale temperature changes during the last deglaciation where the OM-FIT only has a few data points (e.g., HE1 and GI-1d). On the other hand, the intricate interactions of the coupled ocean–atmosphere system during the Holocene, particularly in its later stages and on a millennial scale, underline the need for further research into the driving mechanisms and feedbacks that are crucial in influencing natural temperature variability during the recent interglacial in southern Europe.

Appendix A Table A1

${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$ measurements of Ostolo samples. The ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$ values were corrected for the ice-volume effect during the deglaciation period covered by the Ostolo speleothems. The age indicated for each sample corresponds with the time span covered by the respective calcite blocks sampled from the stalagmites used for fluid inclusion measurements (without taking into account the age model uncertainty). IV $=$ ice volume effect.

FI sample	Water	Water	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	Mean	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	Mean ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	Age ($\mathrm{kyr}\mathrm{BP}$)
	amount	content	($\mathrm{‰}\mathrm{VSMOW}$)	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	std. dev.	error	adjusted for IV
	($\mathrm{\mu }\mathrm{L}$)	($\mathrm{\mu }\mathrm{L}{\mathrm{g}}^{-\mathrm{1}}$)		($\mathrm{‰}\mathrm{VSMOW}$)			($\mathrm{‰}\mathrm{VSMOW}$)
OST1-16.1A	0.52	0.27	$-$50.9	$-$49.7	1.6	1.6	$-$57.1	16.06 $\pm$ 0.06
OST1-16.1B	0.69	0.45	$-$48.6
OST1-15.2A	0.04	0.06	$-$55.0	$-$51.4	3.4	3.4	$-$58.0	15.16 $\pm$ 0.05
OST1-15.2B	0.18	0.18	$-$51.1
OST1-15.2C	0.87	0.31	$-$48.3
OST1-14.6A	0.39	0.39	$-$26.0	$-$25.4	0.9	1.5	$-$31.4	14.57 $\pm$ 0.05
OST1-14.6B	0.86	0.79	$-$24.8
OST1-14.2A	0.11	0.44	$-$43.7	$-$43.4	0.4	1.5	$-$49.1	14.20 $\pm$ 0.02
OST1-14.2B	0.10	0.29	$-$43.1
OST1-13.0	0.57	0.37	$-$32.5	$-$32.5	n/a	2.7	$-$37.0	13.02 $\pm$ 0.04
OST1-10.9A	0.19	0.17	$-$29.3	$-$26.6	3.8	3.8	$-$28.9	10.95 $\pm$ 0.20
OST1-10.9B	0.12	0.40	$-$23.9
OST3-16.4	0.11	0.37	$-$46.1	$-$46.1	n/a	2.7	$-$53.6	16.40 $\pm$ 0.11
OST3-14.3	0.21	0.18	$-$26.0	$-$26.0	n/a	2.7	$-$31.8	14.30 $\pm$ 0.09
OST3-14.1A	0.09	0.07	$-$36.7	$-$35.3	1.6	1.6	$-$40.9	14.11 $\pm$ 0.09
OST3-14.1B	0.19	0.13	$-$33.5
OST3-14.1C	0.20	0.19	$-$35.7
OST3-13.5A	0.14	0.13	$-$34.5	$-$34.9	2.2	2.2	$-$39.9	13.50 $\pm$ 0.09
OST3-13.5B	0.18	0.14	$-$33.0
OST3-13.5C	0.16	0.11	$-$37.3
OST3-13.0	0.20	0.14	$-$31.0	$-$31.0	n/a	2.7	$-$35.4	13.00 $\pm$ 0.08

Table A1

Continued.

FI sample	Water	Water	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	Mean	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	Mean ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	Age ($\mathrm{kyr}\mathrm{BP}$)
	amount	content	($\mathrm{‰}\mathrm{VSMOW}$)	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	std. dev.	error	adjusted for IV
	($\mathrm{\mu }\mathrm{L}$)	($\mathrm{\mu }\mathrm{L}{\mathrm{g}}^{-\mathrm{1}}$)		($\mathrm{‰}\mathrm{VSMOW}$)			($\mathrm{‰}\mathrm{VSMOW}$)
OST3-12.9	0.12	0.12	$-$42.5	$-$42.5	n/a	2.7	$-$46.9	12.90 $\pm$ 0.06
OST3-12.8	0.17	0.14	$-$43.8	$-$43.8	n/a	2.7	$-$48.0	12.80 $\pm$ 0.08
OST3-11.7A	0.11	0.12	$-$26.2	$-$25.0	1.6	1.6	$-$27.9	11.67 $\pm$ 0.02
OST3-11.7B	0.24	0.27	$-$23.8
OST3-11.6A	0.09	0.08	$-$35.2	$-$39.7	6.4	6.4	$-$42.6	11.60 $\pm$ 0.02
OST3-11.6B	0.14	0.15	$-$44.2
OST3-11.5A	0.28	0.25	$-$31.9	$-$32.8	1.2	1.5	$-$35.6	11.49 $\pm$ 0.01
OST3-11.5B	0.32	0.24	$-$33.6
OST3-11.3A	0.12	0.14	$-$39.6	$-$40.0	0.6	1.5	$-$42.6	11.30 $\pm$ 0.02
OST3-11.3B	0.15	0.11	$-$40.5
OST2-16.7	0.11	0.37	$-$46.1	$-$46.1	n/a	2.7	$-$57.1	16.70 $\pm$ 0.07
OST2-16.4A	0.34	0.31	$-$52.4	$-$49.8	2.7	2.7	$-$57.3	16.40 $\pm$ 0.05
OST2-16.4B	0.27	0.23	$-$47.0
OST2-16.4C	0.39	0.90	$-$49.9
OST2-15.8A	0.10	0.12	$-$57.9	$-$59.7	2.6	2.6	$-$66.8	15.80 $\pm$ 0.07
OST2-15.8B	0.09	0.08	$-$61.6
OST2-15.3	0.33	0.17	$-$46.8	$-$46.8	n/a	2.7	$-$53.5	15.31 $\pm$ 0.08
OST2-14.7	0.18	0.14	$-$38.5	$-$38.5	n/a	2.7	$-$44.7	14.71 $\pm$ 0.18
OST2-14.0A	0.10	0.15	$-$55.0	$-$50.5	6.4	6.4	$-$56.1	14.10 $\pm$ 0.09
OST2-14.0B	0.10	0.09	$-$45.9
OST2-13.0A	0.38	0.51	$-$23.2	$-$26.0	3.3	3.3	$-$30.5	13.00 $\pm$ 0.08
OST2-13.0B	0.72	0.94	$-$23.3
OST2-13.0C	0.47	0.61	$-$28.0
OST2-13.0D	0.40	0.44	$-$29.7
OST2-12.9A	0.08	0.14	$-$40.2	$-$40.9	1.0	1.5	$-$45.2	12.89 $\pm$ 0.07
OST2-12.9B	0.09	0.15	$-$41.5
OST2-12.5A	0.09	0.11	$-$50.5	$-$50.9	0.5	1.5	$-$54.7	12.50 $\pm$ 0.10
OST2-12.5B	0.11	0.12	$-$51.2
OST2-12.3A	0.23	0.32	$-$45.3	$-$50.1	6.9	6.9	$-$53.7	12.29 $\pm$ 0.10
OST2-12.3B	0.19	0.23	$-$55.0
OST2-11.8	0.16	0.11	$-$44.5	$-$44.5	n/a	2.7	$-$47.6	11.80 $\pm$ 0.03
OST2-11.65A	0.25	0.17	$-$38.0	$-$36.7	1.8	1.8	$-$39.7	11.65 $\pm$ 0.02
OST2-11.65B	0.43	0.23	$-$35.4
OST2-11.5A	0.18	0.12	$-$28.6	$-$28.7	0.1	1.5	$-$31.5	11.50 $\pm$ 0.01
OST2-11.5B	0.22	0.41	$-$28.7
OST2-10.9A	0.09	0.24	$-$44.4	$-$39.6	6.8	6.8	$-$41.9	10.90 $\pm$ 0.08
OST2-10.9B	0.22	0.38	$-$34.8

n/a: not applicable.

Table A2

FI ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$ measurements of Mendukilo samples. The ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$ values were corrected for the ice-volume effect during the period covered by the Mendukilo speleothems. The age indicated for each sample corresponds with the time span covered by the respective calcite blocks sampled from the stalagmites used for fluid inclusion measurements (without taking into account the age model uncertainty).

FI sample	Water	Water	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	Mean	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	Mean ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	Age ($\mathrm{kyr}\mathrm{BP}$)
	amount	content	($\mathrm{‰}\mathrm{VSMOW}$)	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	std. dev.	error	adjusted for IV
	($\mathrm{\mu }\mathrm{L}$)	($\mathrm{\mu }\mathrm{L}{\mathrm{g}}^{-\mathrm{1}}$)		($\mathrm{‰}\mathrm{VSMOW}$)			($\mathrm{‰}\mathrm{VSMOW}$)
Men2-botA	0.61	0.29	$-$37.7	$-$40.4	2.9	2.9	$-$44.9	12.90 $\pm$ 0.10
Men2-botB	0.54	0.33	$-$37.6
Men2-botC	0.18	0.11	$-$44.4
Men2-botD	0.28	0.21	$-$40.4
Men2-botE	0.30	0.18	$-$41.9
Men2-0A	0.18	0.19	$-$48.9	$-$47.8	1.2	1.5	$-$52.0	12.78 $\pm$ 0.10
Men2-0B	0.32	0.21	$-$47.9
Men2-0C	0.30	0.29	$-$46.4
Men2-5A	0.16	0.27	$-$55.4	$-$57.3	2.7	2.7	$-$61.3	12.65 $\pm$ 0.10
Men2-5B	0.14	0.13	$-$59.3
Men2-10A	0.33	0.20	$-$47.0	$-$48.1	1.5	1.5	$-$52.0	12.51 $\pm$ 0.10
Men2-10B	0.12	0.11	$-$49.2
Men2-17A	0.07	0.08	$-$46.2	$-$47.7	3.8	3.8	$-$51.4	12.32 $\pm$ 0.10
Men2-17B	0.21	0.19	$-$52.9
Men2-17C	0.07	0.09	$-$48.0
Men2-17D	0.35	0.21	$-$43.9
Men2-22A	0.13	0.11	$-$44.5	$-$46.6	4.1	4.1	$-$50.1	12.21 $\pm$ 0.10
Men2-22B	0.15	0.11	$-$51.4
Men2-22C	0.23	0.12	$-$43.9
Men2-27A	0.22	0.22	$-$46.6	$-$48.1	1.9	1.9	$-$51.5	12.08 $\pm$ 0.10
Men2-27B	0.26	0.29	$-$47.7
Men2-27C	0.22	0.22	$-$50.9
Men2-27D	0.12	0.10	$-$47.2
Men2-35A	0.71	0.62	$-$43.1	$-$45.3	1.9	1.9	$-$48.5	11.83 $\pm$ 0.10
Men2-35B	0.40	0.29	$-$45.8
Men2-35C	0.62	0.55	$-$44.8
Men2-35D	0.16	0.14	$-$47.6
Men2-43A	0.15	0.15	$-$46.5	$-$44.3	3.8	3.8	$-$47.1	11.59 $\pm$ 0.08
Men2-43B	0.23	0.26	$-$38.6
Men2-43C	0.14	0.15	$-$46.3
Men2-43D	0.14	0.18	$-$45.8
Men2-47A	0.25	0.36	$-$48.0	$-$50.7	2.4	2.4	$-$53.5	11.50 $\pm$ 0.08
Men2-47B	0.26	0.36	$-$52.1
Men2-47C	0.13	0.18	$-$52.0
Men2-48A	0.15	0.17	$-$35.1	$-$33.5	2.3	2.3	$-$36.3	11.48 $\pm$ 0.08
Men2-48B	0.51	0.59	$-$31.9
Men2-52A	0.18	0.22	$-$39.0	$-$40.5	1.5	1.5	$-$43.2	11.37 $\pm$ 0.08
Men2-52B	0.13	0.18	$-$40.6
Men2-52C	0.67	0.64	$-$41.9
Men2-62A	0.46	0.74	$-$43.5	$-$39.6	3.2	3.2	$-$42.1	11.20 $\pm$ 0.08
Men2-62B	0.30	0.44	$-$36.9
Men2-62C	0.30	0.63	$-$36.9
Men2-62D	0.10	0.12	$-$40.9
Men2-73A	0.18	0.21	$-$35.4	$-$35.2	0.3	1.5	$-$37.5	11.01 $\pm$ 0.06
Men2-73B	0.68	0.76	$-$34.9
Men2-78A	0.23	0.25	$-$43.7	$-$43.6	2.6	2.6	$-$45.9	10.93 $\pm$ 0.06
Men2-78B	0.10	0.18	$-$40.9
Men2-78C	0.16	0.17	$-$46.2
Men2-85A	0.16	0.21	$-$46.3	$-$46.4	1.9	1.9	$-$48.6	10.81 $\pm$ 0.06
Men2-85B	0.27	0.26	$-$48.4
Men2-85C	0.18	0.19	$-$44.6
Men2-92A	0.32	0.3	$-$35.2	$-$37.8	2.2	2.2	$-$39.9	10.69 $\pm$ 0.06
Men2-92B	0.21	0.22	$-$38.6
Men2-92C	0.20	0.25	$-$39.5
Men2-97A	0.26	0.29	$-$43.2	$-$43.3	1.5	1.5	$-$45.3	10.60 $\pm$ 0.06
Men2-97B	0.18	0.21	$-$41.8
Men2-97C	0.14	0.12	$-$44.8
Men2-108A	0.13	0.13	$-$36.1	$-$37.0	1.4	1.5	$-$39.0	10.41 $\pm$ 0.06
Men2-108B	0.35	0.38	$-$38.0
Men2-116A	0.14	0.17	$-$46.2	$-$42.8	3.7	3.7	$-$44.7	10.28 $\pm$ 0.06
Men2-116B	0.33	0.35	$-$38.8
Men2-116C	0.17	0.14	$-$43.5

Table A2

Continued.

FI sample	Water	Water	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	Mean	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	Mean ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	Age ($\mathrm{kyr}\mathrm{BP}$)
	amount	content	($\mathrm{‰}\mathrm{VSMOW}$)	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	std. dev.	error	adjusted for IV
	($\mathrm{\mu }\mathrm{L}$)	($\mathrm{\mu }\mathrm{L}{\mathrm{g}}^{-\mathrm{1}}$)		($\mathrm{‰}\mathrm{VSMOW}$)			($\mathrm{‰}\mathrm{VSMOW}$)
Men2-122A	0.24	0.3	$-$49.3	$-$46.4	2.5	2.5	$-$48.2	10.07 $\pm$ 0.06
Men2-122B	0.12	0.14	$-$45.0
Men2-122C	0.55	0.48	$-$44.9
Men2-128A	0.15	0.17	$-$32.0	$-$32.9	4.1	4.1	$-$34.6	9.96 $\pm$ 0.08
Men2-128B	0.12	0.14	$-$29.3
Men2-128C	0.44	0.41	$-$37.4
Men2-134A	0.32	0.33	$-$41.3	$-$40.6	0.8	1.5	$-$42.2	9.84 $\pm$ 0.08
Men2-134B	0.39	0.46	$-$39.7
Men2-134C	0.25	0.18	$-$40.8
Men2-155A	0.66	0.61	$-$37.2	$-$38.2	6.4	6.4	$-$39.6	9.43 $\pm$ 0.08
Men2-155B	0.35	0.34	$-$32.4
Men2-155C	0.39	0.36	$-$45.1
Men2-162A	0.19	0.15	$-$53.9	$-$49.3	3.2	3.2	$-$50.6	9.29 $\pm$ 0.08
Men2-162B	0.21	0.16	$-$47.3
Men2-162C	0.11	0.13	$-$47.2
Men2-162D	0.21	0.13	$-$48.6
Men2-177A	0.33	0.32	$-$39.3	$-$39.2	0.2	1.5	$-$40.4	9.01 $\pm$ 0.08
Men2-177B	0.35	0.23	$-$39.1
Men2-185A	0.10	0.10	$-$44.2	$-$44.1	3.4	3.4	$-$45.3	8.85 $\pm$ 0.08
Men2-185B	0.19	0.10	$-$47.5
Men2-185C	0.21	0.12	$-$40.7
Men2-192A	0.14	0.16	$-$43.0	$-$38.3	4.3	4.3	$-$39.4	8.72 $\pm$ 0.08
Men2-192B	0.17	0.16	$-$34.6
Men2-192C	0.38	0.19	$-$37.4
Men2-200A	0.38	0.35	$-$40.4	$-$42.3	1.7	1.7	$-$43.4	8.57 $\pm$ 0.08
Men2-200B	0.37	0.24	$-$42.9
Men2-200C	0.22	0.17	$-$43.6
Men2-207A	0.11	0.09	$-$43.7	$-$43.2	0.6	1.5	$-$44.2	8.43 $\pm$ 0.08
Men2-207B	0.19	0.15	$-$42.8
Men2-215A	0.15	0.12	$-$48.7	$-$52.6	4.0	4.0	$-$53.5	8.28 $\pm$ 0.08
Men2-215B	0.18	0.14	$-$56.1
Men2-215C	0.25	0.15	$-$55.8
Men2-215D	0.25	0.18	$-$49.6
Men2-251A	0.23	0.16	$-$42.1	$-$39.4	2.4	2.4	$-$40.2	7.54 $\pm$ 0.08
Men2-251B	0.46	0.35	$-$38.2
Men2-251C	0.17	0.09	$-$37.8
Men2-267A	0.13	0.12	$-$43.7	$-$49.4	4.0	4.0	$-$50.1	7.20 $\pm$ 0.08
Men2-267B	0.16	0.10	$-$52.8
Men2-267C	0.16	0.10	$-$51.7
Men2-267D	0.20	0.13	$-$49.4
Men2-282A	0.15	0.11	$-$44.2	$-$45.6	1.7	1.7	$-$46.3	6.88 $\pm$ 0.08
Men2-282B	0.38	0.27	$-$45.1
Men2-282C	0.20	0.11	$-$47.5
Men2-295A	0.19	0.15	$-$47.5	$-$46.5	3.1	3.1	$-$47.2	6.59 $\pm$ 0.08
Men2-295B	0.22	0.21	$-$42.1
Men2-295C	0.16	0.12	$-$47.0
Men2-295D	0.17	0.13	$-$49.5
Men2-301A	0.11	0.08	$-$49.7	$-$49.9	1.3	1.5	$-$50.5	6.47 $\pm$ 0.08
Men2-301B	0.21	0.14	$-$48.8
Men2-301C	0.14	0.08	$-$51.3
Men2-307A	0.08	0.05	$-$43.5	$-$44.2	1.8	1.8	$-$44.8	6.34 $\pm$ 0.08
Men2-307B	0.16	0.11	$-$46.3
Men2-307C	0.16	0.11	$-$42.8
Men2-310A	0.20	0.18	$-$49.8	$-$46.1	5.1	5.1	$-$46.7	6.28 $\pm$ 0.08
Men2-310B	0.25	0.21	$-$40.3
Men2-310C	0.18	0.12	$-$48.2
Men2-312A	0.38	0.25	$-$45.1	$-$44.2	0.8	1.5	$-$44.8	6.20 $\pm$ 0.08
Men2-312B	0.52	0.35	$-$44.1
Men2-312C	0.33	0.34	$-$43.5
Men5-10A	0.82	0.94	$-$41.8	$-$38.7	2.9	2.9	$-$39.8	8.72 $\pm$ 0.06
Men5-10B	0.19	0.34	$-$36.0
Men5-10C	0.19	0.19	$-$38.4
Men5-20A	0.23	0.24	$-$40.1	$-$39.8	1.6	1.6	$-$40.8	8.58 $\pm$ 0.06
Men5-20B	0.28	0.26	$-$41.2
Men5-20C	0.21	0.19	$-$38.0

Table A2

Continued.

FI sample	Water	Water	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	Mean	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	Mean ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	Age ($\mathrm{kyr}\mathrm{BP}$)
	amount	content	($\mathrm{‰}\mathrm{VSMOW}$)	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	std. dev.	error	adjusted for IV
	($\mathrm{\mu }\mathrm{L}$)	($\mathrm{\mu }\mathrm{L}{\mathrm{g}}^{-\mathrm{1}}$)		($\mathrm{‰}\mathrm{VSMOW}$)			($\mathrm{‰}\mathrm{VSMOW}$)
Men5-30A	0.15	0.16	$-$44.4	$-$40.4	3.5	3.5	$-$41.4	8.45 $\pm$ 0.06
Men5-30B	0.44	0.34	$-$37.7
Men5-30C	0.24	0.18	$-$39.0
Men5-40A	0.39	0.35	$-$48.6	$-$55.3	6.5	6.5	$-$56.3	8.31 $\pm$ 0.06
Men5-40B	0.11	0.10	$-$61.7
Men5-40C	0.11	0.10	$-$55.6
Men5-50A	0.23	0.20	$-$37.8	$-$39.8	2.3	2.3	$-$40.7	8.17 $\pm$ 0.06
Men5-50B	0.50	0.50	$-$42.3
Men5-50C	0.34	0.36	$-$39.2
Men5-60A	0.30	0.25	$-$42.8	$-$44.0	2.8	2.8	$-$44.9	8.04 $\pm$ 0.06
Men5-60B	0.49	0.40	$-$47.1
Men5-60C	0.20	0.18	$-$42.0
Men5-70A	0.14	0.16	$-$39.5	$-$42.6	3.9	3.9	$-$43.5	7.90 $\pm$ 0.06
Men5-70B	0.29	0.29	$-$43.9
Men5-70C	0.09	0.11	$-$47.6
Men5-70D	0.20	0.16	$-$39.6
Men5-75A	0.16	0.16	$-$50.8	$-$49.5	1.2	1.5	$-$50.3	7.84 $\pm$ 0.06
Men5-75B	0.24	0.24	$-$48.5
Men5-75C	0.27	0.27	$-$49.1
Men5-80A	0.18	0.18	$-$33.7	$-$37.4	3.8	3.8	$-$38.2	7.77 $\pm$ 0.06
Men5-80B	0.33	0.30	$-$37.0
Men5-80C	0.23	0.21	$-$41.4
Men5-90A	0.22	0.24	$-$41.4	$-$44.0	4.5	4.5	$-$44.8	7.63 $\pm$ 0.06
Men5-90B	0.20	0.19	$-$48.9
Men5-90C	0.12	0.19	$-$46.4
Men5-90D	0.29	0.23	$-$39.2
Men5-100A	0.12	0.14	$-$37.3	$-$40.5	4.5	4.5	$-$41.3	7.50 $\pm$ 0.06
Men5-100B	0.10	0.14	$-$45.7
Men5-100C	0.21	0.16	$-$38.6
Men5-110A	0.25	0.22	$-$40.9	$-$36.4	4.0	4.0	$-$37.2	7.37 $\pm$ 0.06
Men5-110B	0.20	0.17	$-$33.1
Men5-110C	0.56	0.44	$-$35.1
Men5-120A	0.21	0.21	$-$42.9	$-$44.0	1.1	1.5	$-$44.7	7.23 $\pm$ 0.06
Men5-120B	0.63	0.73	$-$44.8
Men5-120C	0.73	0.87	$-$43.3
Men5-120D	0.32	0.55	$-$45.0
Men5-130A	0.19	0.16	$-$42.1	$-$43.8	3.0	3.0	$-$44.5	7.08 $\pm$ 0.06
Men5-130B	0.34	0.27	$-$47.2
Men5-130C	0.21	0.22	$-$42.1
Men5-140A	0.18	0.23	$-$49.9	$-$48.3	2.3	2.3	$-$49.0	6.93 $\pm$ 0.06
Men5-140B	0.20	0.23	$-$46.7
Men5-150A	0.26	0.26	$-$42.2	$-$43.9	2.7	2.7	$-$44.5	6.75 $\pm$ 0.06
Men5-150B	0.30	0.23	$-$47.0
Men5-150C	0.21	0.19	$-$42.4
Men5-160A	0.35	0.32	$-$40.5	$-$42.4	2.7	2.7	$-$43.0	6.58 $\pm$ 0.06
Men5-160B	0.47	0.39	$-$45.5
Men5-160C	0.30	0.21	$-$41.2
Men5-170A	0.22	0.22	$-$35.0	$-$38.4	2.9	2.9	$-$39.1	6.45 $\pm$ 0.06
Men5-170B	0.27	0.27	$-$37.4
Men5-170C	0.20	0.41	$-$39.8
Men5-170D	0.19	0.14	$-$41.5
Men5-180A	0.39	0.3	$-$44.0	$-$42.3	1.7	1.7	$-$42.9	6.31 $\pm$ 0.08
Men5-180B	0.3	0.3	$-$40.7
Men5-180C	0.45	0.34	$-$42.2
Men5-190A	0.21	0.19	$-$36.5	$-$38.0	2.1	2.1	$-$38.5	6.11 $\pm$ 0.08
Men5-190B	0.28	0.21	$-$39.5
Men5-200A	0.19	0.17	$-$47.2	$-$43.9	4.8	4.8	$-$44.4	5.92 $\pm$ 0.08
Men5-200B	0.17	0.14	$-$46.0
Men5-200C	0.30	0.22	$-$38.3
Men5-210A	0.49	0.45	$-$36.9	$-$38.9	3.1	3.1	$-$39.4	5.72 $\pm$ 0.08
Men5-210B	0.29	0.27	$-$42.4
Men5-210C	0.36	0.3	$-$37.4
Men5-220A	0.34	0.68	$-$39.8	$-$44.7	4.3	4.3	$-$45.2	5.51 $\pm$ 0.08
Men5-220B	0.32	0.45	$-$46.1
Men5-220C	0.24	0.32	$-$48.2

Table A2

Continued.

FI sample	Water	Water	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	Mean	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	Mean ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	Age ($\mathrm{kyr}\mathrm{BP}$)
	amount	content	($\mathrm{‰}\mathrm{VSMOW}$)	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	std. dev.	error	adjusted for IV
	($\mathrm{\mu }\mathrm{L}$)	($\mathrm{\mu }\mathrm{L}{\mathrm{g}}^{-\mathrm{1}}$)		($\mathrm{‰}\mathrm{VSMOW}$)			($\mathrm{‰}\mathrm{VSMOW}$)
Men5-230A	0.78	0.65	$-$36.0	$-$39.1	3.3	3.3	$-$39.5	5.31 $\pm$ 0.08
Men5-230B	0.39	0.39	$-$37.5
Men5-230C	0.26	0.42	$-$39.3
Men5-230D	0.18	0.20	$-$43.5
Men5-240A	0.28	0.23	$-$45.7	$-$46.1	2.8	2.8	$-$46.5	5.11 $\pm$ 0.08
Men5-240B	0.12	0.12	$-$45.8
Men5-240C	0.08	0.13	$-$43.0
Men5-240D	0.19	0.32	$-$49.9
Men5-250A	0.15	0.14	$-$42.3	$-$39.3	2.6	2.6	$-$39.6	4.90 $\pm$ 0.08
Men5-250B	0.17	0.15	$-$37.7
Men5-250C	0.36	0.51	$-$37.8
Men5-260A	0.32	0.29	$-$44.6	$-$43.7	1.2	1.5	$-$44.1	4.70 $\pm$ 0.08
Men5-260B	0.27	0.25	$-$42.9
Men5-280A	0.3	0.49	$-$48.2	$-$48.9	1.0	1.5	$-$49.3	4.29 $\pm$ 0.08
Men5-280B	0.39	0.66	$-$49.6
Men5-290A	0.26	0.21	$-$46.1	$-$45.3	3.2	3.2	$-$45.6	4.08 $\pm$ 0.08
Men5-290B	0.15	0.13	$-$48.0
Men5-290C	0.34	0.30	$-$41.8
Men5-300A	0.19	0.17	$-$40.5	$-$38.6	1.9	1.9	$-$38.9	3.88 $\pm$ 0.08
Men5-300B	0.20	0.16	$-$36.7
Men5-300C	0.46	0.46	$-$38.7
Men5-310A	0.32	0.27	$-$44.1	$-$43.8	1.1	1.5	$-$44.1	3.67 $\pm$ 0.08
Men5-310B	0.23	0.20	$-$42.6
Men5-310C	0.40	0.32	$-$44.7
Men5-330A	0.18	0.14	$-$39.6	$-$41.3	3.8	3.8	$-$41.6	3.26 $\pm$ 0.08
Men5-330B	0.75	0.51	$-$45.7
Men5-330C	0.31	0.52	$-$38.7
Men5-340A	0.21	0.17	$-$40.1	$-$44.0	3.4	3.4	$-$44.2	3.06 $\pm$ 0.08
Men5-340B	0.51	0.51	$-$46.3
Men5-340C	0.77	0.72	$-$45.6
Men5-350A	0.34	0.31	$-$40.7	$-$41.6	2.9	2.9	$-$41.9	2.85 $\pm$ 0.08
Men5-350B	0.14	0.17	$-$44.9
Men5-350C	0.26	0.24	$-$39.4
Men5-360A	0.24	0.19	$-$40.4	$-$37.6	3.1	3.1	$-$37.8	2.65 $\pm$ 0.08
Men5-360B	0.25	0.22	$-$34.3
Men5-360C	0.47	0.38	$-$38.3
Men5-380A	0.29	0.24	$-$37.4	$-$37.1	1.7	1.7	$-$37.3	2.37 $\pm$ 0.06
Men5-380B	0.18	0.16	$-$34.9
Men5-380C	0.42	0.50	$-$39.0
Men5-380D	0.28	0.24	$-$37.2
Men5-390A	0.45	0.40	$-$40.8	$-$47.7	4.9	4.9	$-$47.8	2.25 $\pm$ 0.04
Men5-390B	0.21	0.27	$-$48.5
Men5-390C	0.15	0.16	$-$52.4
Men5-390D	0.15	0.12	$-$49.0
Men5-430A	0.20	0.20	$-$37.1	$-$34.0	3.1	3.1	$-$34.2	1.84 $\pm$ 0.04
Men5-430B	0.25	0.19	$-$33.9
Men5-430C	0.43	0.43	$-$31.0
Men5-440A	0.30	0.25	$-$32.6	$-$36.0	3.4	3.4	$-$36.1	1.73 $\pm$ 0.04
Men5-440B	0.26	0.20	$-$36.1
Men5-440C	0.45	0.37	$-$39.3
Men5-450A	0.13	0.21	$-$46.9	$-$46.5	0.6	1.5	$-$46.6	1.63 $\pm$ 0.04
Men5-450B	0.21	0.24	$-$46.1
Men5-460A	0.39	0.34	$-$39.2	$-$38.0	2.3	2.3	$-$38.2	1.53 $\pm$ 0.04
Men5-460B	0.15	0.15	$-$35.4
Men5-460C	0.50	0.40	$-$39.4
Men5-470A	0.34	0.31	$-$35.3	$-$39.6	3.9	3.9	$-$39.7	1.43 $\pm$ 0.04
Men5-470B	0.47	0.56	$-$40.5
Men5-470C	0.16	0.17	$-$42.9
Men5-480A	0.26	0.29	$-$45.0	$-$41.2	3.8	3.8	$-$41.4	1.34 $\pm$ 0.04
Men5-480B	0.19	0.21	$-$37.4
Men5-480C	0.14	0.12	$-$41.2
Men5-490A	0.18	0.13	$-$40.5	$-$41.5	1.7	1.7	$-$41.6	1.24 $\pm$ 0.04
Men5-490B	0.21	0.17	$-$43.4
Men5-490C	0.38	0.32	$-$40.5
Men5-505A	0.33	0.29	$-$34.4	$-$35.5	1.1	1.5	$-$35.6	1.10 $\pm$ 0.04
Men5-505B	0.15	0.13	$-$35.4
Men5-505C	0.14	0.12	$-$36.6

Table A2

Continued.

FI sample	Water	Water	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	Mean	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	Mean ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	Age ($\mathrm{kyr}\mathrm{BP}$)
	amount	content	($\mathrm{‰}\mathrm{VSMOW}$)	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	std. dev.	error	adjusted for IV
	($\mathrm{\mu }\mathrm{L}$)	($\mathrm{\mu }\mathrm{L}{\mathrm{g}}^{-\mathrm{1}}$)		($\mathrm{‰}\mathrm{VSMOW}$)			($\mathrm{‰}\mathrm{VSMOW}$)
Men5-515A	0.18	0.15	$-$43.4	$-$46.3	3.0	3.0	$-$46.4	0.99 $\pm$ 0.04
Men5-515B	0.14	0.12	$-$49.4
Men5-515C	0.26	0.20	$-$46.0
Men5-525A	0.21	0.20	$-$36.4	$-$38.4	3.2	3.2	$-$38.6	0.90 $\pm$ 0.04
Men5-525B	0.21	0.19	$-$36.8
Men5-525C	0.48	0.39	$-$42.1
Men5-540A	0.22	0.19	$-$43.0	$-$45.0	1.9	1.9	$-$45.1	0.75 $\pm$ 0.04
Men5-540B	0.29	0.25	$-$46.7
Men5-540C	0.18	0.13	$-$45.1
Men5-550A	0.16	0.15	$-$46.7	$-$46.4	1.6	1.6	$-$46.5	0.65 $\pm$ 0.04
Men5-550B	0.41	0.39	$-$44.6
Men5-550C	0.18	0.20	$-$47.8
Men5-560A	0.36	0.37	$-$42.8	$-$45.4	3.3	3.3	$-$45.5	0.55 $\pm$ 0.04
Men5-560B	0.22	0.23	$-$44.3
Men5-560C	0.29	0.26	$-$49.1
Men5-570A	0.56	0.45	$-$44.4	$-$43.1	3.0	3.0	$-$43.3	0.45 $\pm$ 0.04
Men5-570B	0.57	0.43	$-$39.7
Men5-570C	0.43	0.37	$-$45.3
Men5-580A	0.53	0.46	$-$41.7	$-$43.0	2.1	2.1	$-$43.1	0.35 $\pm$ 0.04
Men5-580B	0.47	0.50	$-$45.5
Men5-580C	0.48	0.45	$-$42.0
Men5-590A	0.15	0.13	$-$50.4	$-$46.9	3.2	3.2	$-$47.0	0.25 $\pm$ 0.04
Men5-590B	0.17	0.15	$-$46.0
Men5-590C	0.56	0.59	$-$44.2
Men5-600A	0.26	0.27	$-$47.9	$-$43.7	3.7	3.7	$-$43.8	0.15 $\pm$ 0.04
Men5-600B	0.45	0.54	$-$41.1
Men5-600C	0.54	0.47	$-$42.1

Table A3

Paleotemperatures obtained from ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$ data using the OM-FIT transfer function.

Sample	Age	Mean ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$-	Temp.	Error $\pm$ ($\mathrm{°}\mathrm{C}$)
(stratigraphic	($\mathrm{kyr}\mathrm{BP}$)	adjusted for IV	corrected	OM-FIT
order)		($\mathrm{‰}\mathrm{VSMOW}$)	error $\pm$ (‰)	($\mathrm{°}\mathrm{C}$)
OST2-16.7	16.70 $\pm$ 0.07	$-$57.1	2.7	8.0	2.2
OST2-16.4	16.40 $\pm$ 0.05	$-$57.3	2.7	7.9	2.2
OST1-16.1	16.06 $\pm$ 0.06	$-$57.1	1.6	8.0	2.1
OST2-15.8	15.80 $\pm$ 0.07	$-$66.8	2.6	5.5	2.3
OST2-15.3	15.31 $\pm$ 0.08	$-$53.5	2.7	8.9	2.1
OST1-15.2	15.16 $\pm$ 0.05	$-$58.0	3.4	7.7	2.2
OST2-14.7	14.78 $\pm$ 0.18	$-$44.7	2.7	11.1	2.0
OST1-14.6	14.57 $\pm$ 0.05	$-$31.4	1.5	14.4	1.7
OST3-14.3	14.30 $\pm$ 0.09	$-$31.8	2.7	14.3	1.8
OST1-14.2	14.20 $\pm$ 0.02	$-$49.1	1.5	10.0	1.9
OST3-14.1	14.11 $\pm$ 0.09	$-$40.9	1.6	12.0	1.8
OST2-14.0	14.10 $\pm$ 0.09	$-$56.1	6.4	8.2	2.6
OST3-13.5	13.50 $\pm$ 0.09	$-$39.9	2.2	12.3	1.8
OST2-13.0	13.00 $\pm$ 0.08	$-$30.5	3.3	14.6	1.8
Men2-bot	12.90 $\pm$ 0.10	$-$44.9	2.9	11.8	2.0
OST2-12.9	12.89 $\pm$ 0.07	$-$45.2	1.5	10.9	1.9
OST3-12.8	12.80 $\pm$ 0.08	$-$48.0	2.7	10.3	2.0
Men2-0	12.78 $\pm$ 0.10	$-$52.0	1.5	10.0	2.0
Men2-5	12.65 $\pm$ 0.10	$-$61.3	2.7	7.7	2.3
Men2-10	12.51 $\pm$ 0.10	$-$52.0	1.5	10.0	2.0
OST2-12.5	12.50 $\pm$ 0.10	$-$54.7	1.5	8.6	2.0
Men2-17	12.32 $\pm$ 0.10	$-$51.4	3.8	10.2	2.2
OST2-12.3	12.29 $\pm$ 0.10	$-$53.7	6.9	8.8	2.6
Men2-22	12.21 $\pm$ 0.10	$-$50.1	4.1	10.5	2.2
Men2-27	12.08 $\pm$ 0.10	$-$51.5	1.9	10.2	2.0
Men2-35	11.83 $\pm$ 0.10	$-$48.5	1.9	10.9	2.0
OST2-11.8	11.80 $\pm$ 0.03	$-$47.6	2.7	10.3	2.0
OST3-11.7	11.67 $\pm$ 0.02	$-$27.9	1.6	15.3	1.6
OST2-11.65	11.65 $\pm$ 0.02	$-$39.7	1.8	12.3	1.8
OST3-11.6	11.60 $\pm$ 0.02	$-$42.6	6.4	11.6	2.4
Men2-43	11.59 $\pm$ 0.08	$-$47.1	3.8	11.2	2.1
Men2-47	11.50 $\pm$ 0.08	$-$53.5	2.4	9.6	2.1
OST2-11.5	11.50 $\pm$ 0.01	$-$31.5	1.5	14.4	1.7
Men2-48	11.48 $\pm$ 0.08	$-$36.3	2.3	13.9	1.8
Men2-52	11.37 $\pm$ 0.08	$-$43.2	1.5	12.2	1.9
OST3-11.3	11.30 $\pm$ 0.02	$-$42.6	1.5	11.6	1.8
Men2-62	11.20 $\pm$ 0.08	$-$42.1	3.2	12.5	2.0
Men2-73	11.01 $\pm$ 0.06	$-$37.5	1.5	13.6	1.8
OST1-10.9	10.95 $\pm$ 0.20	$-$28.9	3.8	15.0	1.8
Men2-78	10.93 $\pm$ 0.06	$-$45.9	2.6	11.6	2.0
OST2-10.9	10.90 $\pm$ 0.08	$-$41.9	6.8	11.8	2.5
Men2-85	10.81 $\pm$ 0.06	$-$48.6	1.9	10.9	2.0
Men2-92	10.69 $\pm$ 0.06	$-$39.9	2.2	13.1	1.8
Men2-97	10.60 $\pm$ 0.06	$-$45.3	1.5	11.7	1.9
Men2-108	10.41 $\pm$ 0.06	$-$39.0	1.5	13.3	1.8
Men2-116	10.28 $\pm$ 0.06	$-$44.7	3.7	11.9	2.1
Men2-122	10.07 $\pm$ 0.06	$-$48.2	2.5	11.0	2.0
Men2-128	9.96 $\pm$ 0.08	$-$34.6	4.1	14.4	2.0
Men2-134	9.84 $\pm$ 0.08	$-$42.2	1.5	12.5	1.8
Men2-155	9.43 $\pm$ 0.08	$-$39.6	6.4	13.1	2.4
Men2-162	9.29 $\pm$ 0.08	$-$50.6	3.2	10.4	2.1
Men2-177	9.01 $\pm$ 0.08	$-$40.4	1.5	12.9	1.8
Men2-185	8.85 $\pm$ 0.08	$-$45.3	3.4	11.7	2.0
Men5-10	8.72 $\pm$ 0.06	$-$39.8	2.9	13.1	1.9
Men5-20	8.58 $\pm$ 0.06	$-$40.8	1.6	12.8	1.8
Men2-200	8.57 $\pm$ 0.08	$-$43.4	1.7	12.2	1.9
Men5-30	8.45 $\pm$ 0.06	$-$41.4	3.5	12.7	2.0
Men2-207	8.43 $\pm$ 0.08	$-$44.2	1.5	12.0	1.9
Men5-40	8.31 $\pm$ 0.06	$-$56.3	6.5	9.0	2.6
Men2-215	8.28 $\pm$ 0.08	$-$53.5	4.0	9.6	2.2
Men5-50	8.17 $\pm$ 0.06	$-$40.7	2.3	12.9	1.9
Men5-60	8.04 $\pm$ 0.06	$-$44.9	2.8	11.8	2.0
Men5-70	7.90 $\pm$ 0.06	$-$43.5	3.9	12.2	2.1
Men5-75	7.84 $\pm$ 0.06	$-$50.3	1.5	10.4	2.0
Men5-80	7.77 $\pm$ 0.06	$-$38.2	3.8	13.5	2.0
Men5-90	7.63 $\pm$ 0.06	$-$44.8	4.5	11.8	2.2
Men2-251	7.54 $\pm$ 0.08	$-$40.2	2.4	13.0	1.9
Men5-100	7.50 $\pm$ 0.06	$-$41.3	4.5	12.7	2.1
Men5-110	7.37 $\pm$ 0.06	$-$37.2	4.0	13.7	2.0
Men5-120	7.23 $\pm$ 0.06	$-$44.7	1.5	11.8	1.9
Men2-267	7.20 $\pm$ 0.08	$-$50.1	4.0	10.5	2.2
Men5-130	7.08 $\pm$ 0.06	$-$44.5	3.0	11.9	2.0
Men5-140	6.93 $\pm$ 0.06	$-$49.0	2.3	10.8	2.0
Men2-282	6.88 $\pm$ 0.08	$-$46.3	1.7	11.5	1.9
Men5-150	6.75 $\pm$ 0.06	$-$44.5	2.7	11.9	2.0
Men2-295	6.59 $\pm$ 0.08	$-$47.2	3.1	11.2	2.0
Men5-160	6.58 $\pm$ 0.06	$-$43.0	2.7	12.3	1.9

Table A3

Continued.

Sample	Age	Mean ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$	${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}$-	Temp.	Error $\pm$ ($\mathrm{°}\mathrm{C}$)
(stratigraphic	($\mathrm{kyr}\mathrm{BP}$)	adjusted for IV	corrected	OM-FIT
order)		($\mathrm{‰}\mathrm{VSMOW}$)	error $\pm$ (‰)	($\mathrm{°}\mathrm{C}$)
Men2-301	6.47 $\pm$ 0.08	$-$50.5	1.5	10.4	2.0
Men5-170	6.45 $\pm$ 0.06	$-$39.1	2.9	13.3	1.9
Men2-307	6.34 $\pm$ 0.08	$-$44.8	1.8	11.8	1.9
Men5-180	6.31 $\pm$ 0.08	$-$42.9	1.7	12.3	1.9
Men2-310	6.28 $\pm$ 0.08	$-$46.7	5.1	11.4	2.3
Men2-312	6.20 $\pm$ 0.08	$-$44.8	1.5	11.8	1.9
Men5-190	6.11 $\pm$ 0.08	$-$38.5	2.1	13.4	1.8
Men5-200	5.92 $\pm$ 0.08	$-$44.4	4.8	11.9	2.2
Men5-210	5.72 $\pm$ 0.08	$-$39.4	3.1	13.2	1.9
Men5-220	5.51 $\pm$ 0.08	$-$45.2	4.3	11.7	2.1
Men5-230	5.31 $\pm$ 0.08	$-$39.5	3.3	13.2	1.9
Men5-240	5.11 $\pm$ 0.08	$-$46.5	2.8	11.4	2.0
Men5-250	4.90 $\pm$ 0.08	$-$39.6	2.6	13.1	1.9
Men5-260	4.70 $\pm$ 0.08	$-$44.1	1.5	12.0	1.9
Men5-280	4.29 $\pm$ 0.08	$-$49.3	1.5	10.7	2.0
Men5-290	4.08 $\pm$ 0.08	$-$45.6	3.2	11.6	2.0
Men5-300	3.88 $\pm$ 0.08	$-$38.9	1.9	13.3	1.8
Men5-310	3.67 $\pm$ 0.08	$-$44.1	1.5	12.0	1.9
Men5-330	3.26 $\pm$ 0.08	$-$41.6	3.8	12.6	2.0
Men5-340	3.06 $\pm$ 0.08	$-$44.2	3.4	12.0	2.0
Men5-350	2.85 $\pm$ 0.08	$-$41.9	2.9	12.6	1.9
Men5-360	2.65 $\pm$ 0.08	$-$37.8	3.1	13.6	1.9
Men5-380	2.37 $\pm$ 0.06	$-$37.3	1.7	13.7	1.8
Men5-390	2.25 $\pm$ 0.04	$-$47.8	4.9	11.1	2.3
Men5-430	1.84 $\pm$ 0.04	$-$34.2	3.1	14.5	1.8
Men5-440	1.73 $\pm$ 0.04	$-$36.1	3.4	13.9	1.9
Men5-450	1.63 $\pm$ 0.04	$-$46.6	1.5	11.4	1.9
Men5-460	1.53 $\pm$ 0.04	$-$38.2	2.3	13.5	1.8
Men5-470	1.43 $\pm$ 0.04	$-$39.7	3.9	13.1	2.0
Men5-480	1.34 $\pm$ 0.04	$-$41.4	3.8	12.7	2.0
Men5-490	1.24 $\pm$ 0.04	$-$41.6	1.7	12.6	1.8
Men5-505	1.10 $\pm$ 0.04	$-$35.6	1.5	14.1	1.7
Men5-515	0.99 $\pm$ 0.04	$-$46.4	3.0	11.4	2.0
Men5-525	0.90 $\pm$ 0.04	$-$38.6	3.2	13.4	1.9
Men5-540	0.75 $\pm$ 0.04	$-$45.1	1.9	11.8	1.9
Men5-550	0.65 $\pm$ 0.04	$-$46.5	1.6	11.4	1.9
Men5-560	0.55 $\pm$ 0.04	$-$45.5	3.3	11.6	2.0
Men5-570	0.45 $\pm$ 0.04	$-$43.3	3.0	12.2	2.0
Men5-580	0.35 $\pm$ 0.04	$-$43.1	2.1	12.2	1.9
Men5-590	0.25 $\pm$ 0.04	$-$47.0	3.2	11.3	2.0
Men5-600	0.15 $\pm$ 0.04	$-$43.8	3.7	12.1	2.0

The paleotemperature calculation using Eq. () requires more than one variable to solve; therefore it is necessary to propagate the errors to properly determine the uncertainty. To exemplify this error propagation, $a$ and $b$ are the measured variables (e.g., ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{FI}}$, ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{d}}$, ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}/\mathrm{T}$, and MAAT at “Las Güixas” tourist cave) for the calculation of the OM-FIT temperatures, and ${\mathit{\sigma }}_a$ and ${\mathit{\sigma }}_b$ are the standard deviations of those variables:

if $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\text{OM-FIT}=a+b\left(\text{in the case of addition or substraction}\right)\\ {\displaystyle }& {\displaystyle }\text{OM-FIT}={\displaystyle \frac{a}{b}}\left(\text{in the case of multiplication or division}\right).\end{array}$$

Equation () has subtractions, additions, and a division as mathematical operations. The standard deviation of the OM-FIT (${\mathit{\sigma }}_{\text{OM-FIT}}$) was calculated as follows: $$\begin{array}{c}{\displaystyle }{\mathit{\sigma }}_{\text{OM-FIT}}=\sqrt{{\mathit{\sigma }}_a+{\mathit{\sigma }}_b}\left(\text{in the case of addition or substraction}\right)\\ {\displaystyle }\begin{array}{rl}{\displaystyle \frac{{\mathit{\sigma }}_{\text{OM-FIT}}}{\text{OM-FIT}}}=\sqrt{{\left({\displaystyle \frac{{\mathit{\sigma }}_a}{a}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\mathit{\sigma }}_b}{b}}\right)}^{\mathrm{2}}}& (\text{in the case of multiplication}\\ & \text{or division}).\end{array}\end{array}$$

The errors involved in the interval means were calculated using $${\mathit{\sigma }}_{\text{interval}}=\sqrt{{\displaystyle \frac{{\mathit{\sigma }}_{\text{mean}\mathrm{1}}+{\mathit{\sigma }}_{\text{mean}\mathrm{2}}}{n}}}.$$ The uncertainty of a change is given as the sum of the ${\mathit{\sigma }}_{\text{interval}}$.

Data availability

The speleothem ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ data that support the findings of this study are available as an Excel file in the Supplement, and all the fluid inclusion data are available at Zenodo (10.5281/zenodo.15847530) and will later be integrated into the SISAL database (Wormull, 2025).

The supplement related to this article is available online at https://doi.org/10.5194/cp-21-1235-2025-supplement.

Author contributions

JBW, AM, MB, CPM and MA contributed to design this research project. JLBW, CS, YD, EI, and IC provided the isotopic data. JLBW, CPM, RLE, and HC provided the chronological data. JBW, AM, and EI provided the thin sections and/or contributed to the petrographic characterization. JBW, AM, MB, CPM, MA, RJ, and EI helped during fieldwork. All authors contributed to the writing of the manuscript.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Acknowledgements

We are grateful to the guides and workers of the Mendukilo cave, Miren Larburu and Amaia Govillar, for helping with the monitoring work in the cave and preserving the sampling points. We are also grateful to Manuela Wimmer and Marlene Steck for support during lab work at Innsbruck University and to all people who helped during fieldwork in the Ostolo cave (Iñaki Altzuri and Koldo Sanchez). We would like to acknowledge the use of Servicio de Apoyo a la Investigacion, Zaragoza, and the staff of the IsoTOPIK laboratory at University of Burgos. Isabel Cacho thanks the Catalan Institution for Research and Advanced Studies (ICREA) academia program from the Generalitat de Catalunya.

Financial support

This research has been supported by the Spanish Agencia Estatal de Investigación (AEI) (grant nos. PID2019–106050RB-I00 (PYCACHU) and PID2022-139101OB-I00 (TEMPURA)). We acknowledge support of the publication fee by the CSIC Open Access Publication Support Initiative through its Unit of Information Resources for Research (URICI).The article processing charges for this open-access publication were covered by the CSIC Open Access Publication Support Initiative through its Unit of Information Resources for Research (URICI).

Review statement

This paper was edited by Russell Drysdale and reviewed by two anonymous referees.