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

In the Northern Hemisphere, the last 16.5 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 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 ± 2.8 °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 ± 2.8 °C. The end of this cold phase and the onset of the Holocene are marked by a rapid warming of about 5 °C, reaching a maximum at 11.66 ± 0.03 kyr BP. The OM-FIT record also exhibits abrupt events during the Holocene (e.g., the 8.2 kyr event), which are also reflected in the δ¹⁸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 (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 δ²H_FI composition and uses the δ²H_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 δ²H_FI–temperature relationship is established using monitoring data, and the approach is most effective in settings where δ²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 (δ²H_FI, δ¹⁸O_FI) into temperature may bias temperature estimates; (iii) the relationship between δ²H_FI and δ¹⁸O_FI may have changed over time; and (iv) the δ¹⁸O_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 δ²H_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).

Here, we assess the air temperature evolution in the northern Iberian Peninsula over the last 16.5 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. 1). Although only about 28 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 ± 0.8 °C; > 2000 mm yr⁻¹) compared to Mendukilo (12.2 ± 0.4 °C; ∼ 1365 mm yr⁻¹). This temperature difference is even more pronounced inside the caves: the average annual cave air temperature in Ostolo is 13 °C, while in Mendukilo, it is 8.8 °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 δ¹⁸O and δ²H values of drip waters feeding the stalagmites studied here remain relatively constant, with mean values of −7.7 ± 0.4 ‰ and −45.3 ± 2.9 ‰ Vienna Standard Mean Ocean Water (VSMOW), respectively (1σ 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 δ¹⁸O and δ²H values of drip water are also similarly stable, with mean values of −6.3 ± 0.2 ‰ and −37.8 ± 1.6 ‰ 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) δ²H and δ¹⁸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) δ¹⁸O values of precipitation events and their respective temperature at Villanúa (Giménez et al., 2021). (c) δ²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.

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 km east of the Ostolo and Mendukilo caves. This show cave, located in the Central South Pyrenees (Fig. 1), experiences a transitional Mediterranean–Oceanic climate (Cfb of the Köppen–Geiger climate classification) with a MAAT of 11 °C and around 1100 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 δ¹⁸O and δ²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 δ²H = 7.56 ⋅ δ¹⁸O + 4.33 (n = 210; R² = 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. 2a). In general, the isotopic composition of rainfall correlates with air temperature for the 2-year period (n = 210; R² = 0.44, Fig. 2b) and shows moderate correlation with relative humidity and a weaker correlation with rainfall amount at the event scale when performing Spearman's correlation (r_s; n = 180; between rainfall amount and δ¹⁸O [δ²H]: r_s = −0.27 [−0.25]; between temperature and δ¹⁸O [δ²H]: r_s = 0.70 [0.69]; between relative humidity at the rainfall site and δ¹⁸O [δ²H]: r_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 µ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 δ²H_FI (287 subsamples of Mendukilo stalagmites and 69 of Ostolo samples). Between 0.3 and 2.5 g of calcite was used to ensure a sufficiently high water yield (0.1–1 µ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 δ²H_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). δ²H_FI values are reported in per mil relative to VSMOW. The average long-term precision (1σ uncertainty) of replicate measurements of an in-house calcite standard is ± 1.5 ‰ for δ²H_FI for water amounts between 0.1 and 1 µL.

δ²H_FI is regarded as a more robust proxy of paleotemperature than δ¹⁸O_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, δ¹⁸O_FI values obtained with the Innsbruck FI setup can be inaccurate for samples of low water content (< 0.1 µL; Dublyansky and Spötl, 2009). Therefore, we only used δ²H_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.

Primary FIs were observed in all stalagmite samples (Fig. 3). The Mendukilo samples contain considerably more FIs compared to those from Ostolo, mainly concentrated along growth layers (Fig. 3a). In the Mendukilo stalagmites, primary inter-crystalline (10–30 µm; Fig. 3b) and intra-crystalline (10 to > 100 µm; Fig. 3c) 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. 3c; 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. 3d and e). The intra-crystalline FIs in Ostolo samples are, on average, smaller than those in the Mendukilo stalagmites (10–40 µm) and predominantly exhibit pyriform or rounded shapes (Fig. 3f). Petrographic observations confirm that the FIs in these samples are primary, well preserved, and suitable for their stable isotopic analysis.

Figure 4 (a) δ²H of FI water (δ²H_FI) and δ¹⁸O of calcite (δ¹⁸O_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). δ²H_FI values are corrected for the ice-volume effect (Bintanja et al., 2005), with vertical error bars representing isotope measurements errors and 1σ from repeated measurements. The yellow dashed line in the upper graphs of each panel indicates the annual mean δ²H value in drip water for each cave. Modeled $\mathrm{U}/\mathrm{Th}$ ages with 2σ 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.

4.2 Last deglaciation and the Holocene δ¹⁸O speleothem record

The chronology of the Ostolo stalagmites spans the last deglaciation between 16.5 and 11.7 kyr BP with high precision due to their very high ²³⁸U concentrations (10–80 ppm). The carbonate δ¹⁸O (δ¹⁸O_c) profiles show consistency among the three stalagmites (Fig. 4). 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 δ¹⁸O_c values recorded by the other two stalagmites (Fig. 4). On the other hand, the MEN stalagmites, despite having lower ²³⁸U concentrations (100–350 ppb), also have lower detrital ²³²Th contents, enabling robust age models for both stalagmites. These models cover various intervals of the Holocene and GS-1 with good overlap (Fig. 4), specifically as follows: (i) MEN-2 grew between 12.8 and 6.3 kyr BP, with δ¹⁸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 kyr BP). (ii) MEN-5 spans the last 8.8 kyr and presents prominent negative values during certain short events (e.g., 8.2 kyr 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 ± 2.7 ‰ for δ²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. 3d and e), which restricted subsampling of some individual growth layers. All MEN measurements were duplicated, triplicated, or even quadruplicated. The δ²H values of sub-samples of MEN-2 and MEN-5 (ranging between −34 ‰ and −61 ‰) with water contents of 0.1 to 1 µL replicated within 2.7 ‰.

δ²H_FI values for the Holocene and GI-1 (Tables A1 and A2) are comparable to cave drip waters at Mendukilo and Ostolo caves (Fig. 4). In contrast, values are more negative during GS-1 and GS-2.1a (Fig. 4). GS-2.1a is represented by 8 OST subsamples with a mean δ²H_FI value of −50 ‰. One of these values, dated to 15.80 ± 0.05 kyr BP, is even more depleted (−59.7 ± 2.6 ‰). Values become less negative rapidly at 14.57 ± 0.05 kyr BP (Fig. 4; mean during GI-1: −37 ‰). This trend is interrupted in the three OST stalagmites at 14.13 ± 0.09 kyr BP, leading to more negative values (between −43 ‰ and −51 ‰). During GS-1, the δ²H_FI values decrease again (Fig. 4), 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 kyr BP, MEN-2 and MEN-5 δ²H_FI values show excellent correlation (Fig. 4). There is no significant variation between the Greenlandian (−42 ‰), Northgrippian (−43 ‰), and Meghalayan (−42 ‰). Despite these relatively stable δ²H_FI values throughout the Holocene substages, a short negative shift is identified in two samples dated at 8.28 ± 0.08 (−53.5 ± 4.0 ‰) and 8.31 ± 0.06 (−56.3 ± 6.5 ‰) kyr BP.

5 Discussion

5.1 Interpretation of the δ¹⁸O signal

At isotopic equilibrium, the δ¹⁸O_c value is related to the δ¹⁸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 δ¹⁸O_c records may reflect changes in the δ¹⁸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 δ¹⁸O_c signal is coherent with air temperature changes throughout the deglaciation period (Bernal-Wormull et al., 2021). The overall δ¹⁸O_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 δ¹⁸O_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. 5.4.1), which may account for the negative excursion in the Ostolo δ¹⁸O_c record during Heinrich event 1 (HE1) at 16.2–16.0 kyr BP, with values reaching as low as −8.9 ‰ (Bernal-Wormull et al., 2021; Fig. 4).

Monitoring of the Mendukilo cave suggests that the MEN δ¹⁸O_c record captured an annual signal, which is primarily influenced by rainfall amount (Bernal-Wormull et al., 2023). Conversely, the MEN δ¹⁸O_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 δ¹⁸O_c signal (Bernal-Wormull et al., 2023). Additionally, during the earlier part of the record (13–8 kyr 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 δ¹⁸O_c records is a −0.7 ‰ anomaly (relative to the Holocene mean of −5.4 ‰) observed at 8.11 and 7.00 kyr BP (Fig. 4). These two events of anomalously low δ¹⁸O_c values likely reflect rapid, short-lived decreases in temperature and in the δ¹⁸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 δ²H (δ²H_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 (δ¹⁸O_r and δ²H_r) observed from July 2017 to June 2019 (n = 210). The observed correlation between δ¹⁸O_r, δ²H_r, and air temperature is verified at biannual scale above Las Güixas cave, with significant correlation between MAAT and the values of δ¹⁸O_r and δ²H_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 δ²H_r and air temperature at the time of precipitation (same data series) results in a r_s = 0.69 (p ≪ 0.01) by applying this multiple regression model (Giménez et al., 2021).

δ¹⁸O and δ²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 δ¹⁸O and δ²H values in the ocean water as more fresh water is stored as ice on continents (Lachniet, 2009). δ²H_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 δ¹⁸O (Bintanja et al., 2005) converted to δ²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 (δ²H_d; Ostolo δ²H_d = −37.8 ‰; Mendukilo δ²H_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 δ²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 m for δ¹⁸O, i.e., 1.6 ‰ per 100 m for δ²H_p (Poage, 2001). Uncertainties of the OM-FIT record reflect isotope measurement errors and the selected paleotemperature transfer function. The uncertainties associated with δ²H_FI, δ²H_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 A3 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) δ¹⁸O_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.

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

Our δ²H_FI values (Tables A1 and A2; Fig. 4) 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 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 δ¹⁸O_r and MAAT (δ¹⁸O/T) is 0.55 ± 0.05 ‰ °C⁻¹ (in the case of δ²H_r, the ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}/\mathrm{T}$ is 4.00 ± 0.31 ‰ °C⁻¹) for the “Las Güixas” tourist cave in Villanúa (Giménez et al., 2021), which is consistent with the average European δ¹⁸O/T gradient of 0.59 ± 0.08 ‰ °C⁻¹ (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 kyr BP (Fig. 5). The OM-FIT is calculated using the corrected δ²H_FI values, δ²H_d, MAAT (T_modern), and the modern-day ${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}/\mathrm{T}$ gradient derived from the LMWL of rainfall isotopes:

$$\begin{array}{}\text{(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}}}}\end{array}$$

Equation (1) expresses the fact that the paleotemperatures were lower than present-day temperatures; hence the calculated temperature difference should be subtracted from T_modern. The temperature reconstruction with Eq. (1) is based on the mean relationship of 4.00 ‰ °C⁻¹ (${\mathit{\delta }}^{\mathrm{2}}\mathrm{H}/T_{\text{gradient}}$; Fig. 2c) 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 δ²H_d and T_modern (see Sect. 2) that correspond to the δ²H_FI values of the speleothems from the same cave when using Eq. (1). The final calculated uncertainty in the paleotemperature ranges from 1.6 to 2.6 °C (Table A3).

Figure 6 Paleotemperature reconstructions over the last 16.5 kyr BP, spanning Greenland to SW Europe, along with speleothem δ¹⁸O records from the Iberian Peninsula. (a) NGRIP δ¹⁸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 δ¹⁸O record, compared to the benthic δ¹⁸O record from the Iberian margin (δ¹⁸O_NISA-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 δ¹⁸O from Globigerina bulloides of marine core MD99-2334K (green dashed line; Skinner and Shackleton, 2006). (g) δ¹⁸Oc records from Mendukilo and Ostolo (Bernal-Wormull et al., 2021, 2023). (h) LV5 δ¹⁸O record from Kaite cave (northern Iberia; Domínguez-Villar et al., 2017). (i) GAR-01 δ¹⁸O record from La Garma cave (northern Iberia; (Baldini et al., 2019). (j) SIR-14 δ¹⁸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 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.

5.4 Temperature regime of northern Spain based on OM-FIT

5.4.1 Last deglaciation

The Ostolo cave δ¹⁸O_c and δ²H_FI records (Fig. 4) and the OM-FIT (Fig. 5) 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 δ¹⁸O (Rasmussen et al., 2014) and temperature reconstructions (Kindler et al., 2014) (Fig. 6) 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 ± 0.07 kyr BP and a temperature decrease of approximately 2.5 °C relative to the GS-2.1a average (Fig. 5). This negative temperature excursion is supported by a single point in the δ²H_FI record from OST2 (Figs. 4 and 5). The same δ²H_FI values were obtained in further analyses of the same stalagmite and OST1 but unfortunately with water contents < 0.1 µL. Nevertheless, this anomaly coincides with the most negative δ¹⁸O_c values in stalagmites OST1 and OST2 (Figs. 4 and 5), 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 ± 0.3 kyr 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 δ¹⁸O_c anomaly between 16.2 and 15.9 kyr BP (Pérez-Mejías et al., 2021). Furthermore, this change during HE1 has also been observed in δ¹⁸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. 6). The negative excursion in the Ostolo δ¹⁸O_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. 6). 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 ± 2.8 °C occurred at the onset of GI-1 (Fig. 5). This sudden temperature shift, occurring between 14.78 ± 0.18 and 14.56 ± 0.05 kyr BP in the OM-FIT record (Fig. 5), aligns with data from Greenland ice cores marking the beginning of GI-1 (14.64 ± 0.19 kyr 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 °C) from the West-Central Pyrenees, Millet et al., 2012), and an increase in MAAT (from ca. 12.2 to ca. 18.6 °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 δ²H_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 δ¹⁸O_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 °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 ‰ °C⁻¹, Affolter et al., 2019), estimated a warming of about 5.5 °C (4.1–8.4 °C) (Li et al., 2020) for this transition.

During GI-1, the δ²H_FI record is marked by higher δ²H values and similar temperatures in the OM-FIT record compared to the onset of the Holocene (Fig. 5). As observed in the OST δ¹⁸O_c record, δ²H_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. 6) and OST δ²H_FI values (Fig. 4) 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. 6), 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 δ²H_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 ± 0.09 kyr BP in the OM-FIT record. This suggests that the OM-FIT minimum during GI-1, also registered at 14.10 ± 0.03 kyr BP in the OST δ¹⁸O_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 ± 2.5 °C in the OM-FIT record), occurring just after the GI-1e warm phase (Fig. 5). 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 kyr B.P., the δ²H_FI decrease (Fig. 4) records a cooling of 6.1 ± 2.8 °C in the OM-FIT record (Fig. 5), 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 δ¹⁸O_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 °C; Muñoz Sobrino et al., 2013) and the Central Pyrenees (1.5–2 °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 °C; Li et al., 2020) and the Jura Mountains (4.3 ± 0.8 °C; Affolter et al., 2019) and (ii) a drop in sea-surface temperatures of around 4 °C off the Iberian coast at 12.9 kyr 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 ∼5 °C (Fig. 5), peaking at 11.67 ± 0.02 kyr BP. This warming aligns (within uncertainties) with the record from Greenland ice cores (14.65 ± 0.10 kyr 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 δ¹⁸O_c data during the GI-1/GS-1 and GS-1/Holocene onset transitions is less pronounced compared to OST δ¹⁸O_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 δ¹⁸O_c, as already reported in other speleothem records from this region (e.g., Baldini et al., 2019). In contrast, the δ¹⁸O_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 δ¹⁸O_c compared to the OST δ¹⁸O_c record during GS-1, a cold and dry period (Fletcher et al., 2010). Nevertheless, the MEN δ²H_FI records capture important changes during the GI-1/GS-1 and GS-1/Holocene transitions and correlate quite well with the δ²H_FI data from OST (Fig. 4).

5.4.2 Holocene

As mentioned above, the Holocene section of the OM-FIT record (Fig. 5) is based on δ²H_FI values of the MEN stalagmites (Fig. 4). This record not only captures variability in δ¹⁸O_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 °C in certain intervals. However, these variations are close to the uncertainty range of the OM-FIT record (± 1.6 to ± 2.6 °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 ∼10 kyr BP), albeit with high variability. Between 11.7 and 10.9 kyr BP, OM-FIT records a mean temperature of 13.1 ± 1.9 °C, just 0.6 °C higher than the mean temperature value for the GI-1 of 12.5 ± 1.9 °C (Fig. 5). 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 δ¹⁸O_c record.

The OM-FIT does not capture a clear cooling trend after the Holocene Thermal Maximum (HTM) compared to the δ¹⁸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. 6), 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 δ¹⁸O_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. 6). 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 δ¹⁸O_c values of MEN-2 and MEN-5 (Fig. 6). For example, one of the lowest Holocene OM-FIT temperatures (9.6 °C) occurred at 11.50 ± 0.08 kyr BP (mean temperature at the onset of the Holocene, 12.3 ± 1.8 °C), corresponding within age uncertainties to the Preboreal Oscillation (11.4 kyr) recorded in Greenland ice cores (11.27 ± 0.03 kyr BP, based on the new ice-core chronology – Seierstad et al., 2014) and by MC-FIT in Switzerland (11.37 ± 0.15 kyr BP – Affolter et al., 2019) (Fig. 6). Another example is the 9.2 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 δ²H_FI value of −50.6 ‰ in MEN-2 and an OM-FIT temperature of 10.4 ± 2.1 °C at 9.29 ± 0.08 kyr BP (Fig. 6). However, it is absent from the δ¹⁸O_c record of MEN-2, and previous research suggests that the climate in northern Spain was likely considerably warmer and wetter ∼ 9 ka BP (Baldini et al., 2019; Morellón et al., 2018; Tarrats et al., 2018). This observation supports that the less variable δ¹⁸O_c signal in Mendukilo cave is influenced by short-lived decreases in δ¹⁸O_sw and changes in humidity, while the OM-FIT shows the quantitative change in temperature during the 9.2 kyr event (Fig. 6).

Catastrophic meltwater discharge during the “8.2 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 kyr event overlapped a multi-centennial cool period from 8.29 to 8.10 ± 0.04 kyr BP recorded by MEN δ¹⁸O_c, characterized by an abrupt drop in temperature of about ∼ 3.0 °C between 8.31 ± 0.06 and 8.28 ± 0.08 kyr BP in the OM-FIT record (Fig. 6). 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 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 ∼ 1–1.7 °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 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. 6) and the combination of the carbonate isotopic composition (δ¹⁸O_c and δ¹³C_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 kyr 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 ± 2.8 °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 kyr 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 ± 0.09 kyr 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 δ¹⁸O_c nor OM-FIT reveal significant millennial-scale changes during the Holocene. The 8.2 kyr event is recorded between 8.29 and 8.10 ± 0.04 kyr in the δ¹⁸O_c record, centered between 8.31 ± 0.06 and 8.28 ± 0.08 kyr 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 δ²H measurements of Ostolo samples. The δ²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.

n/a: not applicable.

Table A2 FI δ²H measurements of Mendukilo samples. The δ²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).

Table A3 Paleotemperatures obtained from δ²H_FI data using the OM-FIT transfer function.

The paleotemperature calculation using Eq. (1) 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., δ²H_FI, δ²H_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 σ_a and σ_b are the standard deviations of those variables:

if

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

Equation (1) has subtractions, additions, and a division as mathematical operations. The standard deviation of the OM-FIT (σ_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\hspace{0.17em}\hspace{0.17em}the\hspace{0.17em}\hspace{0.17em}case\hspace{0.17em}\hspace{0.17em}of\hspace{0.17em}\hspace{0.17em}addition\hspace{0.17em}\hspace{0.17em}or\hspace{0.17em}\hspace{0.17em}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\hspace{0.17em}\hspace{0.17em}the\hspace{0.17em}\hspace{0.17em}case\hspace{0.17em}\hspace{0.17em}of\hspace{0.17em}\hspace{0.17em}multiplication}\\ & \text{or\hspace{0.17em}\hspace{0.17em}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 σ_interval.