Reconstructing hydroclimate changes over the past 2500 years using speleothems from Pyrenean caves (NE Spain)

Reconstructing of past hydroclimates at regional scales during the Common Era (CE) is necessary to place the current warming in the context of natural climate variability. Here we present a composite record of oxygen isotope variations during last 2500 years based on eight stalagmites from four caves in the central Pyrenees (NE Spain) dominated by temperature variations, with the amount of precipitation playing a minor role. The dataset is compared with other Iberian reconstructions that show a high degree of internal coherence with respect to variability at the centennial scale. The Roman Period (RP) (especially 0–200 CE), the Medieval Climate Anomaly (MCA), and part of the Little Ice Age (LIA) represent the warmest periods, while the coldest decades occurred during the Dark Ages (DA) and most of the LIA intervals (e.g., 520–550 CE and 1800–1850 CE). Importantly, the LIA cooling or the MCA warming were not continuous or uniform and exhibited high decadal variability. The Industrial Era (IE) shows an overall warming trend although with marked cycles and partial stabilization during the last 2 decades (1990–2010). The strong coherence between the speleothem data, European temperature reconstructions and global tree-ring data informs about the regional representativeness of this new record as Pyrenean past climate variations. Solar variability, likely through its impact on the North Atlantic Oscillation, and major volcanic eruptions appear to be the two main drivers of climate in southwestern Europe during the past 2.5 millennia.

1 Introduction

Global surface temperatures in the first 2 decades of the 21st century (2001–2020) were 0.84 to 1.10 °C warmer than 1850–1900 CE (IPCC, 2021). There is strong evidence that anthropogenic global warming is unprecedented in terms of absolute temperatures and spatial consistency over the past 2000 years (Ahmed et al., 2013; Konecky et al., 2020). On the contrary, pre-industrial temperatures were less spatially coherent, and further work is needed to explain the regional expression of climate change (Mann, 2021; Neukom et al., 2019). Obtaining new and high-quality records in terms of resolution, dating, and regional representativeness is thus critical for characterizing natural climate variability on decadal to centennial scales (PAGES2k Consortium et al., 2017).

High mountains are particularly sensitive regions to climate change, and among them the Pyrenees occupy a crucial frontier position in southern Europe, influenced by both Mediterranean and Atlantic climates. In the Pyrenees, the temperature has increased by more than 1.5 °C since 1882, as shown by the longest time series from the Pic du Midi observatory (Bücher and Dessens, 1991; Dessens and Bücher, 1995). Recent studies confirm this warming trend, showing an increase of 0.1 °C per decade during the last century in the Central Pyrenees (Pérez-Zanón et al., 2017) or even 0.28 °C per decade if only the 1959–2015 period is considered (Observatorio Pirenaico de Cambio Global, 2018). Long-term snow depth observations (starting in 1955) show a statistically significant decline, especially at elevations above 2000 m a.s.l. (López-Moreno et al., 2020). This fact, together with the increase in temperature, has caused the glaciated area in the Pyrenees to decrease by 21.9 % in the last decade (Vidaller et al., 2021), changing from 2060 ha during the Little Ice Age (LIA) to 242 ha in 2016 (Rico et al., 2017). Recent studies on one of the emblematic glaciers in the Pyrenees, the Monte Perdido glacier, show that the current ice retreat is unprecedent in the last 2000 years, as this glacier survived previous warm periods such as the Medieval Climate Anomaly (MCA) and the Roman Period (RP) (Moreno et al., 2021a).

The study of sediment records from lakes in the Pyrenees, where considerable variations in water level, water chemistry, and biological processes have occurred due to changes in effective moisture and temperature, is an excellent approach to reconstruct past climate variability (González-Sampériz et al., 2017). Recently, a comprehensive study in six high-altitude Pyrenean lakes indicates unprecedented changes in the lithogenic and organic carbon fluxes since 1950 CE, suggesting an increase in algal productivity likely favored by warmer temperatures and higher nutrient deposition associated with the Great Acceleration (Vicente de Vera García et al., 2023), a period when human-driven global, social, technological, and environmental changes intensifying dramatically (Steffen et al., 2015). Marine records off the Iberian coast show a clear long-term cooling trend, from 0 CE to the beginning of the 20th century, probably reflecting the decline in Northern Hemisphere summer insolation that began after the Holocene optimum (Abrantes et al., 2017). Unfortunately, it is not possible to record decadal temperature changes from the studied proxies of these lake or marine records, so other archives allowing higher chronological robustness and larger resolution are required.

The Central Pyrenees are largely composed of limestones and host numerous caves, some of which are rich in speleothems, thus making it possible to reconstruct the past climate by studying stalagmites from different caves. Unfortunately, despite the high potential of stalagmite with annual to sub-annual resolution in the CE, it is extremely difficult to obtain high-resolution and well-replicated records. In most cases, the CE period spans only a few centimeters, limiting the number of samples drilled for high-precision U–Th dating (PAGES Hydro2k Consortium, 2017). In addition to this chronological challenge, the interpretation of oxygen isotopes of speleothems ($\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\mathrm{c}}$) from southern Europe is also complex (Moreno et al., 2021b). Recent studies of Pyrenean stalagmites covering the last deglaciation indicate the important role of changes in annual temperature in the variability of $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\mathrm{c}}$ (Bartolomé et al., 2015a; Bernal-Wormull et al., 2021). However, correct interpretation of $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\mathrm{c}}$ proxies requires a sound understanding of the influence of climate variables on carbonate deposition in caves through monitoring (e.g., Pérez-Mejías et al., 2018) and calibration to the instrumental period (Mangini et al., 2005; Tadros et al., 2022).

In this study, we provide high-resolution $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\mathrm{c}}$ data for eight stalagmites from four different caves in the Central Pyrenees, allowing us to construct a stacked curve of climate variability for the last 2500 years with potential regional representativeness. These eight stalagmites allow climate changes during the CE to be studied in reasonably robust chronological framework. Monitoring and calibration of $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\mathrm{c}}$ with instrumental data for the two youngest stalagmites suggests that the $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\mathrm{c}}$ variability primarily reflects annual temperatures, while precipitation (e.g., amount of precipitation, seasonality, source) played a role during certain periods. This new record represents an excellent opportunity to characterize natural temperature changes in this region on decadal to centennial scales for the last 2500 years and compare them with other approaches to examine their regional representativeness.

2 Study sites

2.1 Geological setting, climate, and vegetation

This study of speleothems is located in the central sector of the Pyrenees, in northeastern Iberia (Fig. 1a and b). All caves are located in the Sobrarbe Geopark, close to or at the borders of the Ordesa and Monte Perdido National Park, formed in Mesozoic and Cenozoic limestones and at different altitudes (Fig. 1c). This area has a steep topography due to the high altitudinal gradient and constitutes the largest limestone massif in Europe (with 22 peaks above 3000 m a.s.l.).

Figure 1 (a) Location of regional speleothem records covering last 2500 years to be compared with the samples studied in the Pyrenees (red rectangle, enlarged in b). (b) Location of caves (orange circles) and other nearby records from northern Spain. See legend for the different types of available paleoclimate archives. (c) Location of the four studied caves in the Central Pyrenees of NE Spain in the vicinity of the Ordesa and Monte Perdido National Park. The source base map includes a digital elevation model and hillshade derived from Mapzen Global Terrain and coastline, boundaries, and geographic lines from NaturalEarthData.com.

The climate is Mediterranean according to the Köppen classification. However, the high relief influences the climate of this high-altitude area, which is accurately described as humid sub-Mediterranean because of higher rainfall than the typically Mediterranean climate, particularly for the caves above 1000 m a.s.l. where annual precipitation is above 1000–1200 mm and falls mostly as snow. In lower-altitude caves (e.g., Seso Cave) mean annual precipitation is 900 mm, concentrated in spring and fall. Mean air temperatures range from 0.5 to 15 °C, depending on the altitude.

Around the caves, in the valleys, there are mid-mountain forests dominated by Pinus sylvestris and Quercus ilex, as well as shrublands, whereas the highlands are characterized by exposed rock with sparse vegetation such as meadows.

2.2 Cave locations

Seso cave ($\mathrm{42}\mathit{°}{\mathrm{27}}^{\prime }{\mathrm{23.08}}^{\prime \prime }$ N, $\mathrm{0}\mathit{°}{\mathrm{02}}^{\prime }{\mathrm{23.18}}^{\prime \prime }$ E, 794 m a.s.l.) is formed in the eastern flank of the Boltaña Anticline, close to Boltaña village. The cave developed in unsoluble marly strata between limestone beds of Eocene age. The cave system consists of two longitudinal shallow galleries (2–3 m of limestone thickness over the cave) controlled by the bedding and the main set of joints. Formation of this shallow cave involved the mechanical removal of large amounts of marl under vadose conditions which took place about 60–40 ka BP (Bartolomé et al., 2015b). Subsequently, calcite speleothems formed, which became more abundant during the Holocene. The average annual temperature inside Seso cave is ∼11.8 °C.

Las Gloces cave ($\mathrm{42}\mathit{°}{\mathrm{35}}^{\prime }{\mathrm{40}}^{\prime \prime }$ N, $\mathrm{0}\mathit{°}{\mathrm{1}}^{\prime }{\mathrm{41}}^{\prime \prime }$ W, 1243 m a.s.l.) is located on the border of the Ordesa National Park, next to Fanlo village. The cave formed in limestones of Early Eocene age. The limestone's thickness above the cave is ∼20–30 m. Two galleries form the cave. The upper one preserves phreatic features and hosts the majority of speleothems located in a small room, while vadose morphologies characterize the lower gallery. The average annual temperature where the stalagmites were taken is ∼9.8 °C.

B-1 cave ($\mathrm{42}\mathit{°}{\mathrm{36}}^{\prime }{\mathrm{0.2}}^{\prime \prime }$ N, $\mathrm{0}\mathit{°}{\mathrm{7}}^{\prime }{\mathrm{46}}^{\prime \prime }$ E; 1090 m a.s.l.) is the lower entrance of the Las Fuentes de Escuaín karstic system and acts as the collector of all water drained by the system. This system comprises more than 40 km of galleries and shows a vertical extension of −1150 m. It drains an area of ∼15 km² and developed mostly in Eocene limestones. Since a river runs through the cave, several detrital sequences appear, as well as speleothems, affected by floods. The cave is then well ventilated and shows annual temperature variations in response to the seasonal ventilation changes and seasonal flooding. The studied sample was obtained in a fossil gallery that is not currently influenced by flooding and with an average annual temperature of ∼9.5 °C.

Pot au Feu cave (42°31.48^′ N, 0°14.26^′ W; 996 m a.s.l.) is located in the Irués river valley in the Cotiella massif. The host rock is an Upper Cretaceous limestone. Hydrogeologically, the cave belongs to the high mountain unconfined karst Cotiella-Turbón aquifer but is located in a non-active level. The cave comprises horizontal galleries and small rooms connected by shafts formed by phreatic circulation. Some rooms are well decorated by large speleothems. The limestone thickness over the gallery where the stalagmite was collected is approximately 800 m.

2.3 Cave climate

Understanding the modern microclimatic and hydrological conditions of caves is important for a sound interpretation of speleothem proxy data (Genty et al., 2014; Lachniet, 2009; Moreno et al., 2014). Specifically, the transfer of the stable isotopic signal from the rainfall to the dripwater and eventually to the studied stalagmite is influenced by different processes in the atmosphere, soil, and epikarst. Our preliminary results for the Pyrenees show a seasonal pattern of precipitation isotopes consistent with the annual temperature cycle (Moreno et al., 2021a). These data also suggest an interannual temperature–δ¹⁸O relationship of 0.47 ‰ °C⁻¹ (Giménez et al., 2021) that is only partially compensated by the −0.18 ‰ °C⁻¹ due to the water–calcite isotope fractionation (Tremaine et al., 2011), thus allowing us to use δ¹⁸O in speleothems as a temperature indicator in this region (see also Bartolomé et al., 2015a; Bernal-Wormull et al., 2021).

From the four studied caves, the best monitored one is Seso cave, where a detailed monitoring survey was conducted including analyses of δ¹⁸O variability in rainfall, soil water, dripwater, and farmed calcite (Bartolomé, 2016). Seso cave developed under just few meters of rock, while the other caves are much deeper, allowing a faster response to rainfall variability in Seso dripwaters and speleothems. Monitoring carried out in Seso cave indicates a relationship between temperature and δ¹⁸O of rainfall observed at seasonal scale, while rainfall isotopic composition is slightly modulated by the amount of precipitation (Bartolomé et al., 2015a).

3 Methods

3.1 Speleothem samples

This study is based on eight stalagmites from four different caves in the Central Pyrenees (Fig. 1c, Table 1). The specimens were cut parallel to the growth axis and the central segment was sampled for U–Th dating, stable isotopes (δ¹⁸O and δ¹³C), and $\mathrm{Mg}/\mathrm{Ca}$. Furthermore, the ¹⁴C activity of multiple samples from the top of stalagmites MIC and XEV (both from Seso cave and underneath active drips) was determined in order to detect the atmospheric bomb peak induced by the nuclear tests in 1945–1963.

Table 1 Sample characteristics.

Four small stalagmites were obtained from Seso cave, all showing fine laminations consisting of pairs of dark-compact and light-porous laminae, but they were difficult to count due to their irregular pattern. The four Seso stalagmites show medium to high porosity in some intervals, usually more frequent towards the top. MIC (8.5 cm long) and XEV (26 cm long, composed of two stacked stalagmites – Fig. A1a) were sampled from base to top. In stalagmites CHA (8.5 cm long) and CLA (10.5 cm long), the uppermost interval was discarded due to the poor chronological control and associated with a possible hiatus above a macroscopic discontinuity (Fig. A1a).

Stalagmites ISA (13.5 cm long, with a visual hiatus at 7 cm above the base) and LUC (23.3 cm long, also with a hiatus at 12.5 cm above the base) were sampled in Las Gloces cave (Fig. A1b). Both are candle-shaped examples with a slight tilt in the growth axis above their respective hiatus. One stalagmite, TAR, was obtained from B1 cave and is an overgrowth over an older stalagmite composed of 7.5 cm of white carbonate that is slightly laminated towards the top (Fig. A1c). Finally, a 80 cm long stalagmite (JAR) was obtained from Pot au Feu cave. It is a candle-shaped example that is laminated and lacks macroscopic hiatuses (Fig. A1d).

3.2 Stable isotope and $\mathrm{Mg}/\mathrm{Ca}$ analyses

Samples for stable isotopic (δ¹⁸O and δ¹³C) analyses were microdrilled at 1 mm resolution along the growth axis of seven of the eight speleothems (JAR from Pot au Feu was sampled every 5 mm) using a 0.5 mm tungsten carbide dental bur. The first batch of the isotopic analyses was analyzed at the University of Barcelona (Scientific-Technical Services), Spain, using a Finnigan-MAT 252 mass spectrometer, linked to a Kiel Carbonate Device III, with a reproducibility of 0.02 ‰ for δ¹³C and 0.06 ‰ for δ¹⁸O. Calibration to Vienna Pee Dee Belemnite (VPDB) was carried out by means of the NBS-19 standard. A second batch was analyzed at the University of Innsbruck using a ThermoFisher Delta V Plus isotope ratio mass spectrometer coupled to a ThermoFisher GasBench II. Calibration of the instrument was accomplished using international reference materials, and the results are also reported relative to VPDB. Long-term precision on the 1-sigma level is 0.06 ‰ and 0.08 ‰ for δ¹³C and δ¹⁸O, respectively (Spötl, 2011).

The elemental chemical composition was analyzed in the eight stalagmites (every 1 mm in Las Gloces, Seso, and B1 stalagmites and every 5 mm in JAR from Pot au Feu cave) using matrix-matched standards on an inductively coupled plasma atomic emission spectrometer (Thermo ICAP DUO 6300 at the Pyrenean Institute of Ecology) following the procedure described in Moreno et al. (2010). Reported ratios are from measurement of Ca (315.8 nm) and Mg (279.5 nm), all of which are in radial mode.

3.3 U–Th dating and ¹⁴C bomb peak

A total of 55 samples were prepared for U–Th dating, according to the U and Th chemical procedures described in Edwards et al. (1987). Sample portions characterized by high porosity and voids were avoided to minimize the effect of open system behavior and possible age inversions. From those 55 samples, 45 were measured at the University of Minnesota (USA) and at the Xian' Jiaotong University (China), while 10 samples were analyzed at the University of Melbourne (Australia) (samples of JAR) using the methodology described in Hellstrom (2006). In the three laboratories, measurements were performed using a MC-ICP-MS (Thermo-Finnigan Neptune or Nu Instruments) following previously described methods (Cheng et al., 2013).

Due to the low U content (Table 2), the U–Th ages are not precise enough to obtain an accurate chronology for the recent speleothem growth (see large errors in top samples in Fig. A1). Therefore, the ¹⁴C “bomb peak” method was applied to the MIC and XEV stalagmites that were actively growing in Seso cave at the time of collection (2010 and 2013, respectively), confirmed by $\mathrm{U}/\mathrm{Th}$ ages, albeit of low precision. We drilled 10 and 8 subsamples for MIC and XEV, respectively (Fig. 2a and b), and ¹⁴C activities were measured using a novel online sampling and analysis method combining laser ablation with accelerator mass spectrometry (LA-AMS) at ETH Zurich (Welte et al., 2016). LA-AMS allows us to produce spatially resolved ¹⁴C profiles of carbonate minerals with a precision of 1 % for modern samples. The background measured on ¹⁴C-free marble (F¹⁴C=0.011 ± 0.002) is low, and reference carbonate material is well reproduced. This method relies on the exploitation of the global anthropogenic increase in atmospheric ¹⁴C resulting from nuclear testing predominately in the 1950s and 1960s as a chronological marker in the mid to late 20th century (e.g., Genty et al., 1998; Hua et al., 2012). Atmospheric ¹⁴C concentrations began to rise in 1955 CE, peaking in the Northern Hemisphere (NH) in 1963 CE (Reimer, 2004). Because 80 % to 90 % of the carbon found in most speleothems comes from soil CO₂, this being linked to atmosphere CO₂, it is likely that speleothem ¹⁴C activity is close to the atmospheric ¹⁴C activity or at least to the soil activity (Markowska et al., 2019). Thus, the point where the ¹⁴C concentration begins to rise, the highest concentration point, and the date when the speleothem was removed from the cave (if actively dripping) were used as chronological anchor points (Fig. 2a and b).

Table 2 The ²³⁰Th dating results of the eight stalagmites examined in this study (data from the University of Minnesota, University of Xi'an, and University of Melbourne). Analytical errors are 2σ of the mean. The sample marked in italic font (Isa-6) was discarded due to the high error.

^a Activity ratios determined after Hellstrom (2003) using the decay constants of Cheng et al. (2000). ^b Age in ka before present corrected for initial ²³⁰Th using Eq. (1) of (Hellstrom, 2006) and initial [${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$] of 0.9 ± 0.4. ^c Initial [${}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}$] calculated using corrected age. U decay constants are ${\mathit{\lambda }}_{\mathrm{238}}=\mathrm{1.55125}\times {\mathrm{10}}^{-\mathrm{10}}$ (Jaffey et al., 1971) and ${\mathit{\lambda }}_{\mathrm{234}}=\mathrm{2.82206}\times {\mathrm{10}}^{-\mathrm{6}}$ (Cheng et al., 2013). Th decay constant is ${\mathit{\lambda }}_{\mathrm{230}}=\mathrm{9.1705}\times {\mathrm{10}}^{-\mathrm{6}}$ (Cheng et al., 2013). ^d $\mathit{\delta }{}^{\mathrm{234}}\mathrm{U}=\left(\right[{}^{\mathrm{234}}\mathrm{U}{/}^{\mathrm{238}}\mathrm{U}{]}_{\text{activity}}-\mathrm{1})\times \mathrm{1000}$. ^e $\mathit{\delta }{{}^{\mathrm{234}}\mathrm{U}}_{\text{initial}}$ was calculated based on ²³⁰Th age (T), i.e., $\mathit{\delta }{{}^{\mathrm{234}}\mathrm{U}}_{\text{initial}}=\mathit{\delta }{{}^{\mathrm{234}}\mathrm{U}}_{\text{measured}}\times e^{{\mathit{\lambda }}_{\mathrm{234}}\times T}$. Corrected ²³⁰Th ages assume the initial ${}^{\mathrm{230}}\mathrm{Th}{/}^{\mathrm{232}}\mathrm{Th}$ atomic ratio of 4.4 ± $\mathrm{2.2}\times {\mathrm{10}}^{-\mathrm{6}}$. Those are the values for a material at secular equilibrium, with the bulk earth ${}^{\mathrm{232}}\mathrm{Th}{/}^{\mathrm{238}}\mathrm{U}$ value of 3.8. The errors are arbitrarily assumed to be 50 %. ^f BP stands for “Before Present” where the “Present” is defined as the year 1950 CE.

Figure 2 The ¹⁴C activity (expressed as F¹⁴C, following recommendations made in Reimer, 2004) of the top parts of stalagmites MIC (a) and XEV (b) from Seso Cave. The start of the increase in F¹⁴C and its maximum are recorded at 1955 and 1963 CE, respectively, in both stalagmites. (c) Composite δ¹⁸O record using Iscam with data from MIC and XEV stalagmites.

3.4 Age model

Age models were produced using StalAge software (Scholz and Hoffmann, 2011) for the eight speleothems (Fig. A1) using the U–Th dates presented in Table 2. In the ISA stalagmite, one date was discarded due to the large error (indicated in red in Table 2). During several intervals, two or more stalagmites grew contemporaneously, allowing us to test the reproducibility of the proxy records. We made the a priori assumption that the δ¹⁸O data of the selected stalagmites record a common rainfall and temperature signal, given that these caves were only 20 km apart (Fig. 1c). Following this, the records are combined with Iscam (Fohlmeister, 2012), a method that correlates dated proxy signals from several stalagmites, determines the most probable age–depth model, and calculates the age uncertainty for the combined record.

In order to minimize the effect of different absolute isotopic values and ranges of individual stalagmite data series, we detrended and normalized the δ¹⁸O series using Iscam. By doing so, the interpretation of absolute values will be precluded. Regarding the other parameters that can be changed in Iscam, we used point-wise linear interpolation and 1000 Monte Carlo simulations, and the smoothing window was fixed at 10 years. The stalagmites were included in the Iscam composite record from the oldest to the youngest one as this was the order that provided the highest correlation coefficients: JAR, LUC, ISA, TAR, CHA, CLA, XEV, and MIC. The ISA sample was treated as two parts (ISA top and ISA base) to account for the hiatus, while LUC was regarded as only one as StalAge does not suggest a hiatus in this stalagmite (Fig. A1b). For the two stalagmites that were active when collected, MIC and XEV, we also produced a composite record for the last 200 years using Iscam (Fig. 2c). The use of Iscam software minimized the age uncertainty being lower than the error in the U–Th dates. As an example, for last 600 years, the uncertainty is below 20 years. However, it may reach 100 years for some particular intervals (e.g., the century 1250–1350 CE).

In order to explore correlations among stalagmites from the same caves, we repeated the procedure to obtain a composite record for the four stalagmites from Seso cave (CHA, CLA, XEV, and MIC) (Fig. A2) and the two from Las Gloces cave (ISA and LUC) (Fig. A3). In those two cases, we did not detrend or normalize the individual records since they belong to the same cave and show the same range of δ¹⁸O values. These four records (composite records from Las Gloces and Seso caves and individual stalagmites from Pot au Feu and B1 caves) are shown in Fig. 3 and compared to the final composite record. The composite δ¹⁸O record is used in this paper as a proxy record for the Central Pyrenees climate of the last 2500 years. We have used approximate the onset and end of five time subperiods following previous literature (e.g., Sánchez-López et al., 2016): the end of the RP at 450 CE, the Dark Ages (DA) (450–850 CE), the MCA (850–1250 CE), the LIA (1250–1950 CE), and the Industrial Era (IE) (since 1850 CE).

Figure 3 Comparison of individual δ¹⁸O records from four Pyrenean caves, i.e., (a) Seso, (b) Pot au Feu, (c) B1, and (d) Las Gloces caves, and (e) the composite δ¹⁸O record produced using Iscam (black curve) for the last 2500 years. Generating Seso and Las Gloces curves required Iscam age modeling while Pot au Feu and B1 curves represent only one stalagmite, the age model of which was produced by StalAge modeling. Black double arrows indicate intervals with patterns present in all records. MCA stands for Medieval Climate Anomaly, and IE stands for Industrial Era.

3.5 Statistical analyses

Statistical analyses were carried out using PAST software (Hammer et al., 2001). The δ¹⁸O series and the instrumental climatic series were first resampled (linear interpolation) to obtain the same regular spacing (annual). Following this, correlation was computed using Spearman's rank correlation analysis, a nonparametric measure as an alternative to Pearson correlation analysis. This analysis was preferred to account for nonlinear relationships, with r indicating the correlation coefficient and p value, the probability value of that correlation. The Bonferroni test was applied to prevent data from spuriously appearing as statistically significant by making an adjustment during comparison testing (PAST software; Hammer et al., 2001).

4 Results

4.1 Age models and composite record

4.1.1 Detection of the bomb peak and composite record of the last 200 years

Stalagmites MIC and XEV from Seso cave were actively dripping when removed from the cave (in 2010 and 2013, respectively). Calcite deposited on glass plates placed below the two dripping points and collected seasonally until 2021 demonstrates that the drip water is supersaturated with respect to calcite and suggests that the top layer of both stalagmites was formed during the respective collection year (Fig. 2). Therefore, these two stalagmites were analyzed for their ¹⁴C activity to identify the “bomb peak” and improve the age model.

A strong increase in the ¹⁴C activity is registered in the MIC and XEV stalagmites at 16 mm and 40 mm depth from top (dft), respectively (Fig. 2a and b), with a rise in the fraction modern F¹⁴C, interpreted as the start of the mid-20th century atmospheric bomb peak. This allows defining the year 1955 CE, within ±2-year uncertainties, at 16 mm dft in MIC and 40 mm dft in XEV (Fig. 2). All radiocarbon bomb peaks published from speleothems show that the response of speleothem ¹⁴C activity to the increase in atmospheric radiocarbon activity occurred nearly simultaneously. However, whether the ¹⁴C activity peak in a speleothem can be assigned to the year 1963 CE depends on the soil properties and the thickness of the rock above the cave, as well as the delay in the transfer of the atmospheric ¹⁴C signal to the speleothem (Fohlmeister et al., 2011; Hua et al., 2017). In the case of Seso cave, which is just 2–3 m below the surface and in an area where soils are patchy and thin (Bartolomé, 2016), the transfer of the ¹⁴C signal was likely fast. We therefore place the year 1963 CE, within ±2 years uncertainties, at 11 mm dft in MIC and at 25 mm dft in XEV (Fig. 2a and b).

Since the two stalagmites MIC and XEV are the only ones in this study whose records extend to modern times, we compare them with the instrumental record in order to improve the interpretation of the stable isotope data. Thus, MIC and XEV δ¹⁸O data were first combined using Iscam (Fig. 2c). Using the parameters indicated in Sect. 3.3 but without normalizing the records (both stalagmites belong to the same cave and show the same range of δ¹⁸O values), the correlation of stalagmites MIC and XEV provided by Iscam software (r) is 0.81 (95 % significance). This composite δ¹⁸O record covers the last 200 years and has an amplitude of 0.9 ‰. The main feature (Fig. 2c) is a trend towards less negative values (indicated by a polynomial line in Fig. 2c).

4.1.2 StalAge models and Iscam stack

Age models obtained by StalAge for individual stalagmites indicate that the growth rate was quite stable, except of ISA and LUC, both from Las Gloces cave, where the growth rate changed after hiatuses (Fig. A1b). The temporal resolution of the stable isotope data allows us to explore changes occurring on a decadal scale (Table 1).

Using the parameters for constructing a composite record using Iscam (see Sect. 3), correlation (r) value (95 % significance) of stalagmite JAR and LUC is 0.48, 0.67 between ISA_base and the combined stack of JAR-LUC, 0.65 between ISA_top and the previous stack, 0.74 between TAR and the previous stack, 0.79 between CHA and the previous stack, 0.95 between CLA and the previous stack, 0.71 between XEV and the previous stack, and finally 0.53 between MIC and the previous stack. These values demonstrate a statistically significant correlation among the individual stalagmites and a higher correlation than between the original time series. The composite δ¹⁸O record was compared to the composite records from Seso (Fig. A3) and Las Gloces (Fig. A4) caves and the two individual stalagmites from the other two caves (Fig. 3). This comparison shows that many of the main features of the original records are also well recorded in the composite (Fig. 3). One example is the interval 530–550 CE during the Dark Ages characterized by relatively low δ¹⁸O values in Las Gloces and Pot au Feu cave records (black arrows in Fig. 3) or the interval at the end of the LIA (1675–1750 CE) with less negative δ¹⁸O values in Seso, B1, and Las Gloces cave records (this interval is recorded in five stalagmites: CHA, XEV, TAR, LUC, and ISA; see Fig. A1).

4.2 Individual isotopic and $\mathrm{Mg}/\mathrm{Ca}$ profiles and composite δ¹⁸O record

The isotopic (δ¹⁸O and δ¹³C) and $\mathrm{Mg}/\mathrm{Ca}$ profiles are shown for the eight stalagmites, using their StalAge models (Fig. A1) for the four caves studied (Seso, Las Gloces, B1 and Pot au Feu). In general, δ¹⁸O and δ¹³C are not well correlated ($r=\sim \mathrm{0.3}$–0.4; p values indicating no significant correlation) with the exception of TAR (r>0.8) and CHA (r=0.5). Generally, δ¹³C is better correlated with $\mathrm{Mg}/\mathrm{Ca}$, pointing to a hydrological link of these proxies via changes in prior calcite precipitation (PCP) associated with the longer residence time of the water in the soil and epikarst during dry periods (Genty et al., 2006; Moreno et al., 2010). A similar interpretation was suggested for other Holocene records from northeastern Spanish caves, such as speleothems from Molinos and Ejulve caves in the Iberian Range (Moreno et al., 2017) and records covering the last deglaciation in the Pyrenees (Bartolomé et al., 2015a). However, δ¹³C and $\mathrm{Mg}/\mathrm{Ca}$ are highly variable in absolute values and patterns among caves, and further studies are required to better constrain the climate-proxy transfer functions for two parameters. Therefore, we base our paleoclimate interpretations on the oxygen isotopes which are known to show a more robust response to regional climate change.

The composite δ¹⁸O record for the Central Pyrenees of the last 2500 years is shown in Fig. 3. The highest δ¹⁸O values of last 2500 years were reached during the RP (50 BCE–250 CE). The MCA is characterized by two intervals of relatively high values (900–950 CE and 1150–1250 CE), and the also LIA shows a one such interval (1675–1750 CE). In contrast, the Dark Ages are characterized by consistently low values. In fact, the most negative interval of last 2500 years is reached at ∼520 CE, a well-known cold episode related to volcanic eruptions (see Sect. 5.2). A long interval with low δ¹⁸O values corresponds to the onset of the LIA (1250–1500 CE, with two very negative excursions) as well as the end of the LIA (1750–1850 CE). The most remarkable feature of the MCA and LIA is the large centennial-scale variability. In fact, the LIA has a clear tripartite pattern, with two intervals of low values at the onset and end and less negative values in between. In contrast, the MCA pattern, although also tripartite, is characterized by two intervals of less negative values at the onset and end and a short period of low values in between. An interval with high δ¹⁸O values is observed since 1950 CE (Fig. 3).

5 Discussion

5.1 Interpretation of δ¹⁸O data

Under equilibrium conditions, the δ¹⁸O value of speleothem carbonate is related to just two variables: the δ¹⁸O value of the drip water, and the cave temperature through its control on equilibrium isotope fractionation between water and calcite (Lachniet, 2009). Over the CE, air temperature in a given cave likely changed very little (<1 °C corresponding to ∼0.18 ‰ in stalagmite δ¹⁸O, following Tremaine et al., 2011) (PAGES Hydro2k Consortium, 2017) such that the observed δ¹⁸O variations in these Pyrenean speleothems of more than 1 ‰ are governed primarily by the δ¹⁸O variability of the drip water.

For a constant sea-surface $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}$ value, as it is expected for this time period, event-scale monitoring of the isotopic composition of oxygen in the rainwater ($\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\mathrm{r}}$) in different areas of the Iberian Peninsula constrains some of the drivers of rainfall isotopic fractionation (Moreno et al., 2021a). Recent rainfall monitoring surveys in the Central Pyrenees indicate that the values of $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\mathrm{r}}$ show an interannual dependence on temperature equivalent to 0.47 ‰–0.52 ‰ °C⁻¹, depending on the site (Giménez et al., 2021; Moreno et al., 2021b). This dependence is only partially offset by the empirical value of isotope fractionation during calcite precipitation (−0.18 ‰ °C⁻¹; Tremaine et al., 2011), thus allowing us to consider temperature as one important factor driving δ¹⁸O variability. Apparently, the rainfall amount does not strongly control the isotopic values at event scale, but when analyzing the $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\mathrm{r}}$ variation through time, added to the strong dependence on air temperature, it is clearly observed how the most intense rainfall events together with the longest-lasting rain events (several days) resulted in an isotopic lightening (Giménez et al., 2021). Thus, we consider dripwater $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{dw}}$ to be driving the $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\mathrm{c}}$ signal in the stalagmites, and very likely air temperature and precipitation amount have been modulating its variability along last 2500 years.

The δ¹⁸O composite record, based on the combination of MIC and XEV δ¹⁸O data, provides the opportunity to correlate with instrumental temperature and precipitation amount data (Figs. A4 and A5). It is worth noting that the chronological control of δ¹⁸O data is robust at decadal scale, thus limiting an annual accurate correlation. While temperature records in the region of the studied caves unfortunately are scarce and short (e.g., the Goriz hut station covers only the last 50 years, Fig. A4b), there are two exceptions. First, the homogenized MAAT dataset since 1882 from the Pic du Midi de Bigorre meteorological station (2860 m a.s.l. in the French Pyrenees) (Bücher and Dessens, 1991; Dessens and Bücher, 1995), which started in 1882 CE, is the currently longest one from the Pyrenees (Fig. A4c). Second, the temperature and precipitation reconstruction by Pérez-Zanón et al. (2017) based on 155 stations from the Central Pyrenees starting in 1910 CE (Fig. A4d). Comparing the MIC and XEV δ¹⁸O data with those temperature datasets a significant correlation is found with Pic du Midi de Bigorre mean annual minima temperature (σ_s=0.32; p-value<0.005). It is likely that the other temperature records were too short to generate a significant correlation.

Additionally, when comparing our δ¹⁸O stack with the HadCRU5 reconstruction for the mean Northern Hemisphere temperatures (Morice et al., 2021) (Fig. A4e), the correlation is higher and significant (σ_s=0.49; p-value<0.005). We suspect that the length of this last series (150 years), together with a large spatial scale, leads to a better correlation with the speleothem composite. However, a large part of the variance remains to be explained by other factors (i.e., precipitation changes in source, seasonality, or amount). Using these relationships as a guide and considering all the isotopic change related to temperature change, the observed variation of 0.30 ‰–0.32 ‰ in δ¹⁸O of our composite would represent a change of 1 °C (Fig. A4) that appears quite plausible for the studied period.

The influence of precipitation amount variability on the isotopic composition of speleothem composite is evident from 1970 to 1980 CE, a relatively cool interval in the Pyrenees characterized by a sustained decrease in the amount of precipitation (Pérez-Zanón et al., 2017) (Fig. A5, note the reversed axis for precipitation). For this interval, the relationship between the δ¹⁸O composite and temperature series is reversed, as the low precipitation leads to higher δ¹⁸O values (as if they represented higher air temperatures). On the contrary, a rapid increase in precipitation at ca. 1960 without any important change in temperature, results in a negative peak on the δ¹⁸O speleothem composite (Fig. A5). This shows that in spite of air temperature being an important factor influencing δ¹⁸O variability in speleothems from the Pyrenees, other processes such as the amount of precipitation, its seasonality distribution, or even its source(s) may be also a significant controlling factor (Priestley et al., 2023; Treble et al., 2022), especially when extreme values are reached (very dry or very wet time intervals), as was indicated by rainfall studies in the Pyrenees (Giménez et al., 2021; Moreno et al., 2021b). In any case, MIC and XEV δ¹⁸O data are not significantly correlated with any of the precipitation data from Fig. A5.

Finally, it is important to note that the δ¹⁸O values in the different caves varied at distinct ranges (Fig. 3). Thus, when producing the composite record, the δ¹⁸O profiles of the eight stalagmites were normalized and detrended with the aim of combining different caves. With such a procedure, comparing relative temperature changes coming from different time periods is not possible. Thus, for example, comparing the warming magnitude of the RP with the MCA or with the IE is not feasible since data were obtained from different caves and were previously normalized and detrended. Unfortunately, none of our stalagmites cover continuously from a warm period, i.e., the MCA, to current conditions to compare values. Therefore, the ability of current data to accurately quantify changes in temperature for last 2500 years in the Central Pyrenees is limited. Normalized δ¹⁸O composite record is evaluated in the context of previous local, regional, and global information.

5.2 Climate reconstruction for the last 2500 years

The Pyrenees is a region threatened by global warming, where the impact on biodiversity, elements of the mountain cryosphere such as glaciers or ice caves, and water resources has been increasing in recent decades (https://www.opcc-ctp.org, last access: 1 March 2024). The δ¹⁸O composite constructed using eight speleothems represents the first climate reconstruction based on speleothems for this region covering the last 2500 years and provides an excellent opportunity to reconstruct natural variability and disentangle the main driving mechanisms. We compare it first with other climate series from the Pyrenees and northern Iberia (Sect. 5.2.1) and then with available speleothems from Europe and western Mediterranean to obtain a regional overview (Sect. 5.2.2). Finally, a short discussion about the potential drivers of main observed changes is provided (Sect. 5.2.3).

5.2.1 The last 2500 years in the context of the Iberian Peninsula

Previous climate reconstructions for the CE from the Pyrenees were mostly based on lake records (e.g., González-Sampériz et al., 2017), tree-ring data (e.g., Büntgen et al., 2017), and a few data points from glaciers or ice caves (Moreno et al., 2021a; Oliva et al., 2018; Sancho et al., 2018; Leunda et al., 2019). Despite large variability, these records reveal a clear distinction between relatively cold (DA, LIA) and warm (RP, MCA) periods, which were generally characterized by high and low lake levels, respectively. The differences and similarities among Pyrenean records merit a more detailed evaluation, organized by chronological periods.

The Iberian–Roman period in the Pyrenees

Considering the last 2500 years, the RP stands out as a clear warm period from the speleothem composite record (Fig. 4a). In the Eastern Pyrenees, Redon Lake records low winter–spring temperatures with a warming trend at the end (Pla and Catalan, 2005, 2011), whereas the summer–autumn temperatures show a transition from cold to warm (Pla and Catalan, 2005, 2011). Not many high-resolution Pyrenean lake records exist for this period (e.g., Corella et al., 2016; Vegas-Vilarrúbia et al., 2022) and dendrochronological studies in this mountain range do not cover this time period. Thus, an interesting record to compare with is the A294 ice cave in the Cotiella massif (Sancho et al., 2018). This 9 m thick ice is divided into intervals of low and high snow accumulation, requiring moist and cold conditions to form. The fourth (and last) stage of this ice deposit indicates a high accumulation rate (Fig. 4d), thus a relatively humid and cold period, from 500 BCE to 62 CE. Afterwards, the record stopped reflecting the onset of a warmer and drier climate (Sancho et al., 2018) associated with the RP thermal maximum (Fig. 4a). Recently, not yet published observations indicate the ice deposit grew during the cold and wet years associated with the DA (Miguel Bartolomé, personal communication, 2023). In our speleothem composite, the RP is represented by Las Gloces and Pot au Feu stalagmites that show less negative values (Fig. 3), which suggest rather warm and probably dry conditions in the Central Pyrenees during the RP, particularly between 1 and 200 CE (Fig. 4). This is supported by data showing retreating glaciers in the Pyrenees at that time (Moreno et al., 2021a).

Figure 4 (a) Central Pyrenean δ¹⁸O composite record for the last 2500 years based on eight stalagmites from four caves. Blue bars mark intervals of δ¹⁸O values below −0.75, while yellow bars mark those with δ¹⁸O values above +0.75 (note this composite record was obtained from normalized records, so it varies among −3 and 3 without possibility of direct translation to absolute δ¹⁸O values). (b) Rainfall reconstructed from calcite layers from Montcortés lake in the Pre-Pyrenees (Corella et al., 2016). (c) Salinity reconstructed from geochemical data from Estanya lake in the Pre-Pyrenees (González-Sampériz et al., 2017; Morellón et al., 2012, 2011). (d) Snow and ice accumulation in ice cave A294 in the Cotiella massif of the Central Pyrenees (Sancho et al., 2018). (e) Pyrenean temperature reconstruction based on tree-ring data (Büntgen et al., 2017). MCA stands for Medieval Climate Anomaly, and IE stands for Industrial Era.

The Dark Ages in the Pyrenees

This period is characterized in our speleothem composite by cold and wet climates starting ca. 300 CE, with two particular cold events at 500–650 CE and 750–850 CE and a warmer and drier interval in between (650–750 CE) (Fig. 4a). Pyrenean lake records also point to cold and wet conditions but with a high heterogeneity and low resolution, thus preventing a detailed characterization of this time period (González-Sampériz et al., 2017). For example, Estanya Lake recorded a dominant dry climate between 500 and 750 CE (Fig. 4c), changing to higher lake levels afterwards (Morellón et al., 2009), a pattern that is quite coherent with the speleothem composite. Proxy data from Redon Lake suggest cold winter–spring temperatures in the Eastern Pyrenees during the DA (Pla and Catalan, 2005, 2011).

The Medieval Climate Anomaly in the Pyrenees

The large centennial-scale temperature variability recorded by the speleothem composite is particularly well expressed for the MCA and the LIA, with three distinct intervals of temperature changes (yellow and blue bands in Fig. 4a), thus revealing a more complex pattern as previously inferred by lower-resolution records (e.g., Moreno et al., 2012; Sánchez-López et al., 2016). The MCA has been interpreted as a “warm and dry” climate regime in the southern Pyrenees (Morellón et al., 2012) (Fig. 4c), characterized by low lake levels and more abundant xerophytic vegetation. Our new data show, however, that a colder (and maybe wetter) interval between 950 and 1050 CE separated two clear warm periods before (900–950 CE) and after (1150–1250 CE; Fig. 3) this point. This cold interval was also identified in the Redon Lake record as a sudden cooling about 1000 years ago (Pla and Catalan, 2005). Interestingly, this cold century was not observed by an increase in heavy precipitation in the Montcortés lake record (Fig. 4b).

The Little Ice Age in the Pyrenees

The LIA climate variability is well characterized in the Pyrenees thanks to records from glaciers, such as moraines associated with glacier advances, but also due to historical documents such as pictures or old photographs (Oliva et al., 2018). The available information indicates that the LIA glaciers in the Pyrenees occupied 3366 ha in 1876 and just 810 ha in 1984 and that these glaciers have lost 23.2 % of their volume considering only data from 2011 to 2020 (Hughes, 2018; Vidaller et al., 2021). In many Pyrenean valleys, more than one moraine belt was assigned to the LIA (García-Ruiz et al., 2014), but unfortunately the discontinuous character of these landforms and difficulties in dating them does not allow us to resolve the internal pattern of the LIA in the Pyrenees. A recent compilation of records across the Iberian mountains proposed several climate phases during the LIA (Oliva et al., 2018) that are well correlated with our speleothem composite (Fig. 4a). A first cooling phase lasted from the onset of the LIA (ca. 1200 CE) until 1480 CE, followed by relatively warm conditions from 1480 to 1570 CE. A second phase of gradual cooling occurred until 1600 CE followed by very cool conditions lasting until 1715 CE and coinciding with the Maunder Minimum (1645–1715 CE). In our speleothem composite, this interval is well defined as a cold period, but it was not the one with minimum δ¹⁸O values of the LIA (Fig. 4a). The first half of the 18th century was characterized by warm conditions, supported by many records compiled by Oliva et al. (2018). After 1760 and until the end of the LIA (ca. 1850 CE), a climate deterioration and more frequent extreme climate events were described. This last cold phase is also captured by the speleothem composite and may correspond to the Dalton Minimum (1790–1830 CE). It is characterized by large climate variability and lasted until about 1850 CE.

The Industrial Era in the Pyrenees

The Industrial Era (IE), defined as the last 150 years, is characterized in the Pyrenean speleothem composite by low temperatures that started to increase at about 1950 CE (Fig. 4a) in response to the Great Acceleration (Steffen et al., 2015) (yellow band in Fig. 4). This increase in temperature is well recorded in other Pyrenean climate archives, such as glaciers or lake records. Thus, the last 150 years were marked by a gradual glacier retreat since 1850 CE that accelerated specially after 1980 CE, considered a “tipping point” in glacier retreat not only on a Pyrenean scale (López-Moreno et al., 2016) but also on a global scale (Beniston et al., 2018). A decrease in heavy rainfall (Fig. 4b) and an increase in salinity (Fig. 4c) are well defined in Montcortés and Estanya lake records, respectively, during the IE, indicating a decrease in the amount of precipitation in a likely drier scenario. Besides these two lake records, high-altitude lakes show a significant increase in primary productivity during the last decades (Vicente de Vera García et al., 2023). These recent results demonstrate the combined impacts of climate change and increased human pressure in the Pyrenees. Coherently, the last 50 years are characterized by generally enriched δ¹⁸O values in our speleothem record (yellow bands in Fig. 4). However, the last two decades (our record ends in 2013, the year XEV sample was collected) are not the ones with the highest δ¹⁸O values (Fig. 4a) as also observed in tree-ring data from the Spanish Central Pyrenees (Büntgen et al., 2017) (Fig. 4e). One potential explanation for the lack of exceptionally high δ¹⁸O values would be a slight increase in precipitation amount. Thus, precipitation reconstruction for the Pyrenees during the last 2 decades indicate slightly higher values than those of previous decades (Pérez-Zanón et al., 2017, Fig. A5). Other factors, such as changes in the precipitation source or type (e.g., dominance of Atlantic frontal rainfall versus Mediterranean convective episodes) may be also behind the recorded δ¹⁸O values of last decades.

Figure 5 Comparison of European and W Mediterranean speleothem records covering the last 2500 years. (a) Winter NAO reconstruction based on growth rate of Irish speleothems (Baker et al., 2015). (b) Precipitation variability reconstructed for western Portugal (Thatcher et al., 2022), (c) the Balearic Islands (Cisneros et al., 2021), and (d) Morocco (Ait Brahim et al., 2019). Temperature variation reconstructed from (e) Sweden (Sundqvist et al., 2010), (f) central Europe (Affolter et al., 2019; Fohlmeister et al., 2012), (g) the European Alps (Mangini et al., 2005), and (h) northern Iberia (Martín-Chivelet et al., 2011). (i) Central Pyrenean δ¹⁸O composite record (this study). Black arrows indicate intervals of well-reproduced patterns during the Dark Ages and the Little Ice Age cold intervals. MCA stands for Medieval Climate Anomaly, and IE stands for Industrial Era.

5.2.2 Temperature variability in western Europe and the western Mediterranean during last 2500 years

The PAGES2k European temperature record is the most recent compilation of the last two millennia at European scale (PAGES 2k Consortium, 2013) and it is coherent with our speleothem composite for the Central Pyrenees (Fig. 6). This comparison shows a synchronicity for several of the warmest intervals of the CE, such as the first centuries CE in the RP, the 1150–1250 CE period within the MCA, and the last decades (marked as orange bars in Fig. 6). There are very few high-resolution speleothem records in Europe covering the CE (Comas-Bru et al., 2020); we selected nine speleothems records in Europe and northern Africa that cover the last 2500 years with robust chronology and decadal resolution (Fig. 5). One of these records is interpreted as North Atlantic oscillation (NAO) variability (Baker et al., 2015), three are paleo-precipitation reconstructions (Ait Brahim et al., 2019; Cisneros et al., 2021; Thatcher et al., 2022) and the other five are reflecting paleo-temperature variations (Affolter et al., 2019; Fohlmeister et al., 2012; Mangini et al., 2005; Martín-Chivelet et al., 2011; Sundqvist et al., 2010). Considering these differences in the interpretation and the fact these records are from different regions with different climates (from Sweden to Morocco), dissimilar profiles of paleoclimate variability can be expected. Still, some features are comparable and can be discussed to obtain a super-regional picture.

Figure 6 Global records and forcing mechanisms. (a) Volcanic forcing represented by the (nss)S (ppb) in the NEEM ice core (blue line). (b) Changes in the irradiance as a consequence of Northern Hemisphere volcanic eruptions (Sigl et al., 2015) (brown bars). (c) Sunspot numbers (Usoskin et al., 2014). (d) Probability of paleofloods in European temperate regions (Benito et al., 2015). (e) Worldwide tree-ring compilation (green line, running average width window = 15) (Sigl et al., 2015). (f) Temperature reconstruction from Europe, compiled by the PAGES2k group (red line, running average width window = 15) (PAGES 2k Consortium, 2013). (g) The NAO reconstruction for the central Iberian Peninsula (purple line) and the 95 % (light gray band) uncertainty intervals. (h) The Central Pyrenean δ¹⁸O composite record (this study). Light brown bars indicate warming periods during the Roman Period, at the end of the MCA, and in recent decades. Light blue bands mark cooling events during the DA, while dark blue bands mark solar minima (Oort, Wolf, Spörer, Maunder and Dalton). MCA stands for Medieval Climate Anomaly, and IE stands for Industrial Era.

The Roman period in Europe – western Mediterranean

In Europe, and particularly in the Mediterranean region, the RP is well known as a warm period (e.g., McCormick et al., 2012). The average sea surface temperature in the western Mediterranean Sea was 2 °C higher than the average temperature of later centuries (Margaritelli et al., 2020). Our composite, with high values of normalized δ¹⁸O values during the whole RP, and particularly from 0–200 CE, agrees with the scenario of warm temperatures (Fig. 5i). Speleothem data from the Balearic Islands (Cisneros et al., 2021) indicate a transition from humid to dry conditions along the Iberian–Roman Period (Fig. 5c). The dry period at the end of the RP in the Balearic record appears to be in agreement with a new speleothem record from northern Italy (Hu et al., 2022), suggesting that the observed drying trend was a possible contribution to the collapse of the Roman Empire in 476 CE. Records from Morocco (Ait Brahim et al., 2019) in contrast mark a humid trend at the end of the RP (Fig. 5d). Similarly, an increase in humidity was observed in southern Iberia during the RP (Jiménez-Moreno et al., 2013; Martín-Puertas et al., 2009), thus reflecting a large spatial heterogeneity in precipitation amount when comparing records from the north and south of the Mediterranean basin.

The Dark Ages in Europe – western Mediterranean

After the RP, the cold Dark Ages started (450–850 CE). Part of this period is known as the “Late Antique Little Ice Age” (LALIA), lasting from 536 to 670 CE, characterized by specially cold conditions in Europe (Büntgen et al., 2016). Our speleothem composite shows in general cold and wet conditions but with centennial-scale variability during the DA (Fig. 5). Three clear intervals can be defined, following the δ¹⁸O pattern of our composite, as well as speleothem records from the Alps (Mangini et al., 2005) and Central Europe (Affolter et al., 2019; Fohlmeister et al., 2012): an initial cooling phase corresponding to the LALIA (ca. 500–650 CE), a warming phase (ca. 650–750 CE), and a final cooling phase right before the onset of the warming associated with the MCA (ca. 750–850 CE). A δ¹³C speleothem record from three northern Iberian caves (Martín-Chivelet et al., 2011) shows a warming trend in the DA period but with internal variability that, within dating uncertainties, can be related to the three phases defined above (Fig. 5i). It is worth noting that the period with the most negative δ¹⁸O values recorded in the speleothem composite from the Pyrenees corresponds to the LALIA decades, a cooling period that provoked widespread social disruption in Europe, famine, and episodes of epidemic diseases (Peregrine, 2020).

The Medieval Climate Anomaly in Europe – western Mediterranean

The MCA was one of the warmest periods in continental Europe (and the western Mediterranean, Lüning et al., 2019) of the CE, usually dated to 900 to 1300 CE and characterized by warm (Goosse et al., 2012) and relatively dry conditions (Helama et al., 2009). The MCA was also characterized by a general glacier retreat, mainly associated with a decline in precipitation amount in the Alps (Holzhauser et al., 2016) and the Pyrenees (Moreno et al., 2021a). This scenario is supported by speleothem records from Europe and the western Mediterranean (Fig. 5), which all point to generally warm (Affolter et al., 2019; Fohlmeister et al., 2012; Mangini et al., 2005; Martín-Chivelet et al., 2011; Sundqvist et al., 2010) and/or dry conditions (Ait Brahim et al., 2019; Baker et al., 2015; Thatcher et al., 2022), even leading to speleothem growth stops as seen, for example, in the Balearic record (Cisneros et al., 2021). Previous studies have emphasized the complexity of the spatial and seasonal structure of the MCA in Europe (Goosse et al., 2012). The selected speleothem records underscore this complexity, particularly considering that in our Pyrenean composite one of the periods marked as cold-wet occurred during the MCA, ca. 950–1050 CE (Fig. 5). We propose that this cold interval represents the climate response to the Oort solar minimum in the Pyrenees, a time period characterized by low number of sunspots spanning 1010 to 1050 CE (Bard et al., 2000).

The Little Ice Age in Europe – western Mediterranean

The LIA is well known in Europe and the western Mediterranean region, characterized by cold temperatures and relatively humid conditions as recorded, for example, in chironomid-inferred summer temperatures (Ilyashuk et al., 2019), Mediterranean sea surface temperatures (Cisneros et al., 2016), the advance of alpine glaciers (Holzhauser et al., 2016), and the rise in lake levels (Magny, 2013). The LIA cooling, however, was not continuous and uniform in space and time. Regarding temperatures, many of the available reconstructions from the Alps (Trachsel et al., 2012), Scandinavia (Zawiska et al., 2017), and other regions of Europe (Luterbacher et al., 2016) provide evidence for a main LIA cooling phase that was divided into three parts: two cold intervals with a slightly warmer episode in between, with the most severe cooling during the 18th century (Ilyashuk et al., 2019). This pattern is also found in the two temperature records Iberian speleothems (this study and the one from Martín-Chivelet et al., 2011) and a temperature record from the Alps (Mangini et al., 2005) (Fig. 5, marked by arrows). The other European speleothem records show only two phases during the LIA: a longer and intense cooling period followed by a warming (Fig. 5, Affolter et al., 2019; Fohlmeister et al., 2012; Sundqvist et al., 2010). A tripartite pattern is recorded by humidity-sensitive speleothems from Portugal, with wet–dry–wet conditions in excellent agreement with the cold–warm–cold pattern in the Pyrenean record (this study), supporting the concept that this pattern is controlled by changes in intensity and N–S migration of the Azores High (Thatcher et al., 2022).

The Industrial Era in Europe – western Mediterranean

Between about 1870 CE and today, an increase in temperature is detected by European speleothem records (Fig. 5), as previously shown by the retreat of European glaciers (Beniston et al., 2018) and tree ring summer temperature records (Büntgen et al., 2011) as well as drought reconstructions (Büntgen et al., 2021). The impacts in Europe and the western Mediterranean of the current global warming trend, accelerated during last 50 years, are becoming more and more evident (Jacob et al., 2018; Naumann et al., 2021).

5.2.3 Drivers of past temperature variability in the Pyrenees

The good correlation and synchronicity between the PAGES2k European record and the Pyrenean composite (marked as orange bars in Fig. 6) supports the interpretation of temperature being the dominant factor in controlling the speleothem record. This centennial-scale correlation can be extended to a worldwide tree-ring compilation (Sigl et al., 2015) pointing to the presence of common warm periods in the Central Pyrenees. Interestingly, if precipitation amount was the dominant factor controlling the δ¹⁸O speleothem composite, it would be difficult to find a common signal at regional or even continental scale, as indicated by the overall good correlation shown in Fig. 6.

It is worth to mention the good correlation with several especially cold periods at decadal scale (blue bars in Fig. 6), such as the event at 540–550 CE (registered at 520 CE in the speleothem record) or two cold spikes at 800–850 CE at the end of the DA. We proposed that the cold event at ca. 540 CE (the coldest of the speleothem record) is related to a cataclysmic volcanic eruption that took place in Iceland in 536 CE and spewed ash across the Northern Hemisphere, together with the effect of two other massive eruptions in 540 and 547 CE (Fig. 6b, Sigl et al., 2015). An unprecedented, long-lasting, and spatially synchronized cooling was observed in European tree ring records associated with these large volcanic eruptions, corresponding to the LALIA period (Büntgen et al., 2016). Therefore, volcanic events, at least the large ones such as that in 536 CE in Iceland, have an effect driving temperature variations in the Pyrenean region.

There is also an evident synchrony between the European record and the Pyrenean speleothems in several of the more recent coldest intervals of the MCA and the LIA (dark blue bars in Fig. 6), probably a regional response to minima in solar irradiance as these events correspond to minima in sunspot numbers (Fig. 6c; Usoskin et al., 2014, 2016): 1010–1050 CE (Oort minimum), 1280–1350 CE (Wolf minimum), 1450–1550 CE (Spörer minimum), 1645–1715 CE (Maunder minimum), and 1790–1820 CE (Dalton minimum). Because variations in total solar irradiance are relatively small (on the order of a few tenths of W m⁻²), the mechanism that could result in a detectable cooling remains uncertain (Gray et al., 2010). While some studies discarded the idea that there has been a strong direct radiative influence of solar forcing on Northern Hemisphere temperatures in the past millennium (Schurer et al., 2014), other authors demonstrated a connection among solar variability and climate throughout changes in the large-scale atmospheric circulation of the Northern Hemisphere, such as the North Atlantic Oscillation (NAO) (Martin-Puertas et al., 2012). The NAO was proposed as a plausible mechanism to explain climate changes in Europe during the MCA vs. LIA periods through the study in combination of proxy records and model simulations (Trouet et al., 2009; Mann et al., 2009). Thus, it was postulated that the MCA–LIA transition included a weakening of the Atlantic Meridional Overturning Circulation (AMOC) and a transition to more negative NAO conditions, resulting in a strong cooling of the North Atlantic region and an increase in the storm intensity (Trouet et al., 2012).

Such a connection among solar irradiance and temperature over Europe is then manifested through a change in the pressure gradient in the Atlantic that resembles a negative phase of the NAO and results in lower temperatures over Europe but also in a southward shift in the storm tracks enhancing precipitation amount over central and southern Europe (Swingedouw et al., 2011). As solar irradiance decreases, colder temperatures over the Northern Hemisphere continents are observed, especially in winter (1 to 2 °C), in agreement with historical records and proxy data for surface temperatures (Shindell et al., 2001). Coherently, most episodes of flooding in northwestern and northern Europe region match with multi-decadal periods of grand solar minima and are thus also associated with the negative phase in the NAO index (Benito et al., 2015) (Fig. 6d).

In Iberia, the NAO index was embraced to explain the dryness during the MCA as observed in low-resolution records (Moreno et al., 2012). Further studies based on proxy reconstructions in Iberia explained those MCA–LIA differences by using interactions between the NAO and the East Atlantic (EA) phases (Sánchez-López et al., 2016). In that line, the persistence of NAO phases, for example, the dominance of positive index during medieval times has been questioned (Ortega et al., 2015), and the interactions with other atmospheric modes, together with the non-stationary character of these atmospheric patterns, are nowadays important issues to contemplate when providing a NAO reconstruction (Comas-Bru and Hernández, 2018). In Fig. 6g, the NAO reconstruction provided using a lake record in NW Iberia (Hernández et al., 2020) is compared with the speleothem Pyrenean record demonstrating a good connection. Not surprisingly, the lack of correlation for some periods could be associated with (i) chronological uncertainties in both records, (ii) different seasons recorded by the analyzed proxies, and (iii) distinct influence of NAO in the west and east of the IP.

6 Conclusions

The eight stalagmites presented in this study document significant climate changes on the decadal scale in the Central Pyrenees during the last 2500 years for the first time. The δ¹⁸O composite record is dominated by regional temperature changes, as suggested by monitoring data and by the correlation with observational temperature data from the Pyrenees and at a hemispheric scale. The precipitation amount may also play a role as shown by the comparison with Pyrenean lake records.

On a regional scale, there is a good agreement with other Pyrenean and Iberian records (lake levels, tree rings, and glacier advances) indicating a regional representativity of this new record. The RP stands out as a clear warm period, while the DA, MCA, and LIA exhibit a high centennial-scale variability with cold (e.g., 520–540 CE and 1750–1850 CE) and warm intervals (e.g., 900–950 CE and 1150–1250 CE) modulated by increases and decreases in the precipitation amount, respectively. In spite of temperature increases since 1950 CE, known as the Great Acceleration within the IE, the last 2 decades are not the ones with the highest δ¹⁸O values in the composite record, likely pointing to the secondary role played by precipitation amount.

On a European scale, the Pyrenean composite is in robust agreement with the PAGES2k temperature reconstructions, particularly during warm events. It shows some similarities with other speleothem reconstructions from the Alps and central and northern Europe, pointing to coherent patterns all over the continent for cold and wet and warm and dry periods of last 2500 years. This coherence is supported by synchronous changes with the sunspot number (low temperatures during solar minima), the North Atlantic Oscillation index (low NAO correlates with cold and wet decades), and major volcanic eruptions (e.g., several eruptions during LALIA).

Appendix A

Figure A1 Polished slabs, age–depth model using StalAge (left) and proxy profiles versus age (right) for the stalagmites used in this study arranged by cave (a Seso, b Las Gloces, c B1, and d Pot au Feu caves). Correlation coefficients among the three proxies are indicated based on Pearson correlation. Horizontal lines represent the age error for every data point, following StalAge uncertainty.

Figure A2 Construction of the composite δ¹⁸O record for Seso cave. In the upper graph, the individual δ¹⁸O profiles of the four Seso stalagmites are presented, using their StalAge models (XEV in red, CHA in blue, CLA in orange and MIC in green). Some records overlap (mostly between XEV and CHA and XEV and MIC). The composite δ¹⁸O record for Seso cave is shown in purple on the same y axis as the individual curves. The Central Pyrenees δ¹⁸O composite record is shown at the bottom of the graph.

Figure A3 Construction of the composite δ¹⁸O record for Las Gloces cave. In the upper graph, the δ¹⁸O profiles of the two Las Gloces stalagmites are presented using their StalAge models (ISA in red and LUC in blue). The composite δ¹⁸O record for this cave is shown in purple curve on the same y axis as the individual curves. The Central Pyrenees δ¹⁸O composite record is shown at the bottom of the graph. MCA stands for Medieval Climate Anomaly, and IE stands for Industrial Era.

Figure A4 Correlation of (a) composite δ¹⁸O record from MIC and XEV stalagmites with instrumental temperature records at local, regional, and global levels. (b) Mean annual average temperature (MAAT) from Goriz hub (AEMET data), (c) MAAT from Pic du Midi de Bigorre (Bücher and Dessens, 1991; Dessens and Bücher, 1995), (d) minimum annual temperature from the Pyrenees from the AEMET series (Pérez-Zanón et al., 2017), and (e) MAAT anomalies (with respect to 1961–1990) using the HadCRUT 5.0.1.0. dataset (Morice et al., 2021). Panel (f) shows δ¹⁸O values of the Pyrenees composite record (in a) compared to North Hemisphere mean annual temperatures (in e), showing a significant correlation.

Figure A5 Correlation of (a) composite δ¹⁸O records from MIC and XEV stalagmites with instrumental precipitation records at regional levels. (b) Annual precipitation from Goriz hub (AEMET data) and (c) precipitation anomalies from the Pyrenees from AEMET series (respect to 1961–1990 years) (Bücher and Dessens, 1991; Dessens and Bücher, 1995). No significant correlation is observed.