Geochronological reconstruction of the glacial evolution in the Ésera valley (Central Pyrenees) during the Late Pleistocene and Holocene

Since the MIS (Marine Isotope Stage) 4, glaciers in the Pyrenees were subject to intricate dynamics, encompassing some advances and retreats not always aligned with the fluctuations observed in other European glaciers. The Ésera valley, located in the Central Pyrenees (northern Spain), provides a distinctive opportunity to reconstruct past climate and glacier fluctuations in high-mountain regions since the Pyrenean Last Glacial Maximum (PLGM). Previous studies of glacial evolution in this valley have employed a variety of methods, including the analysis of glacial lake sediments and detailed geomorphological studies of glacial landforms that are integrated here with new. Measurements of cosmogenic ¹⁰Be exposure ages, together with calculations of the equilibrium line altitude (ELA), with the objective of reconstructing the evolution of the Ésera glacier and the associated environmental implications. Following the PLGM, at approximately 75 ka in the Ésera valley, the Ésera glacier commenced a period of retreat during part of the Marine Isotopic Stage (MIS) 3, reaching a point of stabilisation at approximately 47 ka at the location of the Pllan d'Están proglacial lake. Subsequently, a new glacial advance resulted forming the Llanos del Hospital moraine (∼ 16 ka), a glacial deposit located a lower altitude in the valley than Pllan d'Están lake. During that time interval, and including the global Last Glacial Maximum (LGM), lacustrine sediments deposited at Pllan d'Están paleolake indicate that sedimentation took place in a subglacial environment. Following the conclusion of the Oldest Dryas period (OD; ∼ 16 ka) and continuing into the Early Holocene, the Ésera glacier underwent a rapid retreat. The Little Ice Age (LIA) represented the last cold period documented in the Ésera valley, after which the glacier has exhibited a persistent retreatment. The ELA analyses indicate that the temperature in the Ésera valley increased by 3.6 ± 0.45 °C over the past 16 ka, which resulted in the retreat of the glacier front from 1750 to 3000 metres above sea level (m a.s.l.).

1 Introduction

Glaciers are highly effective indicators of climate variability due to their sensitivity to climate change, particularly in relation to temperature and precipitation. Consequently, their thickness and surface variations, in spite they also depend on the catchment hypsometry, are regarded as one of the most informative proxies for climate change in mountainous regions (Beniston et al., 2018; Braithwaite and Hughes, 2020). A number of approaches have been developed for the reconstruction of past glacier fluctuations, which are typically based on the analysis of glacial geomorphologies, including till or moraine deposits, erratic boulders and polished bedrock. One approach, arguably among the most widely used methods over the past couple of decades, is to date glacial landforms with cosmogenic isotopes. The application of cosmogenic exposure dating has become a fundamental technique for the reconstruction of glacial history, including the advance and retreat of glaciers (Balco, 2020; Ivy-Ochs and Briner, 2014). Despite the discontinuous nature of glacial deposits, the dating of erratic boulders on moraine ridges, or polished bedrock via exposure dating, is a commonly used method for studying the dynamics of glaciers and the climate changes that occurred during the Late Pleistocene. When possible, that discontinuous information is complemented by the study of lacustrine sediments deposited on lakes whose dynamic is associated to the glaciers (Moreno et al., 2010).

The equilibrium line altitude (ELA) is also a highly effective method for inferring past climate variations at regional and local scales, identifying the altitudinal boundary between the accumulation zone (where processes favoring snow deposition are dominant) and the ablation zone (where the glacier loses ice mass). The ELA is highly sensitive to climatic changes, and also it depends on the shape and hypsometry of the glaciers (even though the evolution of the glaciers is also controlled by climate change; Quesada-Román et al., 2020; Sagredo et al., 2014; Zemp et al., 2007). This allows to deducing past climate variability through the variations in reconstructed paleoELAs (Dahl and Nesje, 1992; Sissons and Sutherland, 1976; Sutherland, 1984). In order to determine the paleoELA, it is necessary to have accurate knowledge of the glacier extension at a given time period. This is usually achieved through the geomorphological mapping and dating of till and moraines. There are several methods to determine the ELA for a glacier, but the most accurate method is the Area × Altitude Balance Ratio (AABR; Benn et al., 2005), although this method presents challenges associated with the reconstruction of past glacier extensions, particularly the identification and dating of geomorphological features that are not always well preserved to reconstruct the surface of glaciers, overthought this approach remains a valuable method for understanding past thermal changes in mountainous environments (Pellitero et al., 2019).

The Pyrenees currently host the southernmost glaciers in Europe (Grunewald and Scheithauer, 2010), which are classified as very small due to their limited extension (< 0.5 km²; Huss and Fischer, 2016). However, during the Pyrenean Last Glacial Maximum (PLGM) period, i.e. 60–70 ka BP, corresponding to Marine Isotopic Stage (MIS) 4, few glaciers covered a significant portion of the region above 800 m a.s.l. in the Western Pyrenees, constituting a prominent element of landscape modeling in the Mediterranean area (Oliva et al., 2019), as reported by many OSL dates (Lewis et al., 2009) but also by TCN dating (see Delmas et al., 2022 for review). Numerous previous studies have proposed that deglaciation in the Pyrenees would have commenced shortly after the PLGM thus occurring earlier than in other European mountains. Subsequent to the PLGM, a number of glacial phases have occurred. For example, during the MIS 4, glacial stages are preserved in the Ara valley (Bartolomé et al., 2021) or in the Gállego valley (Lewis et al., 2009). Other phases are less well preserved and are the subject of some controversy or inconsistency with other European regions. For example, the Last Glacial Maximum (LGM) was defined as an interval centered on 21 ka (Mix et al., 2001), characterised by a global increase in glacier extent (Hughes et al., 2013). However, in the Pyrenees and other Iberian mountains (Palacios et al., 2011, 2012, 2016), this period was also characterised by low temperatures but high aridity, which seemingly prevented a significant growth of glaciers, as occurred in many other Mediterranean mountains (Allard et al., 2021; Batbaatar et al., 2018). Accordingly, just few moraines were dated from this period in certain Eastern Pyrenean valleys (Delmas, 2015; Delmas et al., 2008, 2011, 2022; Pallàs et al., 2006, 2010), while no date is yet available for the Western Pyrenees.

The Ésera valley is currently the most glaciated valley in the Pyrenees (Rico et al., 2017; Vidaller et al., 2021) due to its geographic features, such as its high elevation, and its location, under the complex influence from both the Mediterranean and Atlantic climates. In the context of current climate change, measurements of the extent and thickness of the modern glaciers (namely, Aneto, Maladeta and Tempestades glaciers) were carried out to investigate the recent evolution of the glaciers in the valley (Campos et al., 2021; Jiménez-Vaquero, 2016; Mora et al., 2006; Pastor Argüello, 2013; Rico et al., 2017; Vidaller et al., 2021, 2023). In contrast, there is a paucity of studies that have focused on the past extent of the glaciers. Indeed, there are only two examples that have dated moraines and polished bedrock from the Maladeta massif via cosmogenic exposure dating (Crest et al., 2017; Reixach, 2022). In addition, the use of different geomorphological maps (1 : 50 000) of the Maladeta massif may prove beneficial in determining the evolution of the glaciers in the Ésera valley (Bordonau, 1992, 1993; Chueca-Cía and Julián-Andrés, 2008; García-Ruiz et al., 1992; Martínez de Pisón, 1989, 1990). In a recent study, (Vidaller et al., 2024a) produced a detailed geomorphological map of the Maladeta massif, with the objective of facilitating the identification of glacier deposits and landforms at a higher resolution (1 : 15 000). The sediments from the two paleolakes in the valley have also been the subject of study. The glaciolacustrine rhythmites of Barrancs lake (2360 m a.s.l.) have been observed to cover the last 3 ka BP (Copons et al., 1997; Copons and Bordonau, 1997). On the other hand, the sedimentary sequence of the Pllan d'Están paleolake (1840 m a.s.l.) evidenced main climatic variations that have occurred over the past 47 ka (Vidaller et al., 2024b).

In this paper, we present the measurements of cosmogenic ¹⁰Be from granitic erratic blocks, boulders from moraine ridges and a polished quartzite bedrock from the Ésera valley, together with ELA calculations, with the aim of reconstructing the chronological sequence of the Late Pleistocene and Holocene Ésera glacier fluctuations deglaciation of the Ésera glacier. Subsequently, we examine the paleoenvironmental conditions, integrating our cosmogenic exposure ages with the preceding sedimentary study of the Pllan d'Están paleolake. Additionally, paleotemperature estimation from different glacial stages in the Ésera valley will facilitate a more quantitative determination of the paleoenvironmental condition.

2 Study area

The Ésera valley is situated in the Central Pyrenees in northern Spain (Fig. 1a). At the headwaters, the southernmost slope of the valley corresponds to the northern slope of the Maladeta massif, which is home to the largest concentration of the last remaining glaciers in the Pyrenees. The massif comprises over 40 peaks exceeding 3000 m a.s.l., including some of the highest in the Pyrenees, such as the Aneto peak, which reaches 3404 m a.s.l. A distinctive feature of this valley is its circuitous route around the Maladeta massif, which results in a notable shift in orientation from north-northwest to south-southeast to north-northwest in a span of less than 10 km.

Figure 1 Location of the Ésera valley and the Maladeta massif. (a) Location of the study area in the European context. Pink pin indicates the location of the Maladeta massif in the Pyrenees. Hillshade obtained from the European Environment Agency. (b) Major geological characteristics of the Maladeta massif and the Ésera valley (modified from Vidaller et al., 2024a). Names of the main locations considered in this work are indicated. Viewpoint of panel (c) is also marked. Hillshade obtained from Spanish National Geographic Institute. (c) Overview of the headwaters of Ésera valley and northern slope of the Maladeta massif with the key locations mentioned in this paper. The picture is taken by looking to the south (note the position of the North as indicated by the arrow).

In geological terms, the study area is located at the northern boundary of the tardihercinic granitic batholith of the Maladeta massif, which is one of the batholiths that constitute the Axial Zone of the Pyrenees, and it is situated above the Paleozoic sedimentary and metamorphic units. The lower areas of the valley are characterised by the presence of Devonian limestones, which are intercalated with Carboniferous shales and quartzites. The northern slope of the valley has been subject to influence from both the Hercynian and Alpine orogenies, as well as the outcropping of Ordovician shales, quartzites and sandstones. The northern slope does not attain the same elevations as the southern slope. The upper part of the Maladeta massif is composed of granite and granodiorite of the batholith (García-Sansegundo et al., 2013; Ríos-Aragüés et al., 2002) (Fig. 1b).

Attending to the geomorphological characteristics, the Maladeta and Ésera landscape has been mostly modelled by the action of glaciers along the Quaternary (Copons and Bordonau, 1997; García-Ruiz et al., 1992; Martínez de Pisón, 1989, 1990). The glacial and periglacial modelling is dominating in the massif and especially visible in granites, while limestones often show karstic landforms (Vidaller et al., 2024a). The variation in lithologies from granite at the top and limestone at the bottom of the valley allows us to identify the source area of the glacial deposits (Fig. 1c), which in turn enables the reconstruction of paleoglacial activities in the region. The prevalence of limestone in extensive areas has resulted in the formation of diverse karst features, including sinkholes, poljes and ponors, which have a significant impact on the hydrological system. Numerous depressions within the valley have been filled with glacial or lacustrine sediments (and may also have undergone erosion due to glacial advances), which, in conjunction with glacially-transported deposits (e.g. till, moraines), serve as a valuable archive of past climate records.

3 Methods

3.1 Terrestrial cosmogenic nuclide exposure dating

To determine the timing and position of glaciers during the last glacial cycle in the Ésera valley, a total of 14 samples were collected for cosmogenic exposure dating. The samples included 7 granitic till boulders (PDE-1, -2, -4, -5, -6, -7, -12), 6 boulders on moraine ridges (PDE-8, -9, -10, -11, -13, -14) and one sample (PDE-3) from a polished quartzite bedrock (Fig. 2, Table 1; see Table S1 in the Supplement for field data). In order to select the most appropriate samples, the detailed geomorphological map of the Maladeta massif (Vidaller et al., 2024a) was consulted. Granite boulders large enough to avoid post-depositional movements due to glacier melting and those associated to moraine ridges were considered for dating. In addition, special attention was paid to ensure that the boulders had a flat surface to guarantee the highest solar incidence. The samples were collected from a range of elevations, with ten originating from the middle and lower reaches of the Ésera valley at elevations between 1744 and 1995 m a.s.l., and the remaining four from the upper part of the valley in the Maladeta massif at elevations between 2039 and 2634 m a.s.l. (Fig. 2). Two samples (PDE-1 and -2) were collected in the vicinity of Pllan d'Están paleolake and are directly comparable to the chronologies reconstructed from the lacustrine sequence in our previous study (Vidaller et al., 2024b).

Figure 2 Location of the cosmogenic samples and photographs of the boulders and the landscape for the geomorphological context. In the map, the dots represent the location of each sample, the colour refers to each type of sample (yellow for moraine ridge boulder, orange for till boulders, green for erratic boulders and purple for polished surface). The pink squares of the map show the area of the landscapes showed in each photo. Four samples (PDE-1, -2, -5 and -9) are showed in detail. The simplified geomorphological map was obtained and modified from Vidaller et al. (2024a). Hillshade obtained from Spanish National Geographic Institute.

Table 1 Results of cosmogenic ¹⁰Be measurements and exposure ages from Ésera valley. Uncertainties are in 1σ: Internal (Int.) errors include only analytical errors and external (Ext.) errors include also systematic errors (such as errors associated with production rate and half-life). When compared to other geochronological data, external errors must be considered.

^a All sample and field details are summarized in Table S1. ^b Details of laboratory data are described in text and summarized in Table S2. ^c Calculated assuming no erosion and no snow shielding. ^d Calculated assuming erosion rates at 3 mm ka⁻¹ and snow shielding considering respective boulder heights, where for PDE1 and PDE2, which are located in the valley bottom, no boulder heights were considered (see text). ^e Distance between the sample location and the front of Aneto glacier in 2023 (except for PDE7, its distance relative to the Maladeta glacier in 2023).

The samples were processed at the Cosmogenic Nuclide Dating Laboratory in CENIEH, Burgos, Spain, and nuclide measurements were carried out at ASTER, CEREGE, Aix-en-Provence, France (see Table S2 for details of sample preparation and AMS data). The sample preparation protocol was modified from previously established procedures (Fujioka et al., 2022; Kohl and Nishiizumi, 1992). In summary, the rock samples were crushed and sieved into fractions of grain size between 0.25 and 0.50 mm, followed by density separations at 2.62 and 2.68 g cm⁻³, leaching by aqua regia and finally etching by diluted HF to obtain pure quartz. A quantity of between 12 and 30 g of the pure quartz was spiked with approximately 0.25 to 0.27 mg of ⁹Be and dissolved in concentrated HF. Subsequent to dissolution of the quartz, small aliquots were obtained for analysis of aluminium by ICP-OES at CENIEH. Beryllium (and aluminium for a bedrock sample, PDE-3) was isolated from unwanted cations via sequential anion and cation exchange chromatography. Subsequently, hydroxides were precipitated with Suprapur (NH₃ aq.). The Be- (and Al-) hydroxides were then transported to ASTER, where they underwent further processing to become BeO (and Al₂O₃). Finally, they were compressed into Cu cathodes with Nb powder (and Ag powder for Al₂O₃) for AMS measurements.

AMS ¹⁰Be $/$ ⁹Be (and ²⁶Al $/$ ²⁷Al) measurements were conducted using the ASTER 5MV accelerator at the Centre de Recherche en Écologie et Gestion de l'Environnement (CEREGE) (Arnold et al., 2013). The ¹⁰Be $/$ ⁹Be measurements were normalised against the STD-11 standard, with a nominal ratio of 1.191 × 10⁻¹¹, while the ²⁶Al $/$ ²⁷Al measurements were normalised against SM-Al-1. One standard was used, with a nominal ratio of 7.401 × 10⁻¹² (Arnold et al., 2010; Braucher et al., 2015; Merchel and Bremser, 2004). The measurements of chemistry procedural blanks indicate a ratio of 2.6–4.0 × 10⁻¹⁵ for ¹⁰Be $/$ ⁹Be and 1.6–5.5 × 10⁻¹⁵ for ²⁶Al $/$ ²⁷Al. These values are typically ∼ 2 % and ∼ 1 % of the measured ratios for the samples, respectively. The ¹⁰Be (and ²⁶Al) concentrations were calculated from the measured ratios after standard normalisation and blank correction. The final errors were calculated from AMS uncertainties, including counting statistics, standard reproducibility, the error on the standard nominal ratio, and blank correction. Additionally, a 1 % error in the ⁹Be carrier concentration was considered, or in the case of ²⁶Al, a 3 % error in the intrinsic ²⁷Al measurement by ICP-OES. These errors were calculated in quadrature. The surface exposure ages were calculated using version 3 of the CRONUS-Earth online calculator (http://hess.ess.washington.edu/, last access: February 2024; Balco et al., 2008), employing the LSDn scaling scheme (Lifton et al., 2014) (see Table S3). The snow shielding factor was calculated in accordance with the methodology outlined by (Gosse and Phillips, 2001), with an assumed average snow density of 0.4 g cm⁻³ and the utilisation of monthly snow thickness data from the Ésera valley region between 2007 and 2023 (AEMET database). The calculated age estimates are considered to be the most accurate when the erosion rate of 3 mm ka⁻¹ and the snow shielding correction are applied. It should be noted that PDE-1 and PDE-2 are situated in the lowest part of the valley (Fig. 2), where snow accumulation is likely to have been greater than the average. Therefore, for these two samples, the calculated age estimates are the most precise when the snow shielding factors are calculated without considering the boulder heights.

It is noted that the mass attenuation length for high-energy neutrons responsible for spallation reactions producing cosmogenic ¹⁰Be may be shorter in snow (109 g cm⁻²) than the conventional value assumed here (160 g cm⁻² ), due to hydrogen moderation (Zweck et al., 2013). Calculating exposure ages with this lower value shows minimal differences for most samples (within 1σ errors) and only slight differences for Samples PDE-1 and PDE-2 (within 2σ errors; see Table S4 in the Supplement).

3.2 PaleoELAs

The reconstruction of paleoELAs of different glacial stages during the Late Pleistocene and Holocene represents a valuable tool for calculating temperature variations (Pellitero et al., 2019). Although there are numerous methods for determining the ELA, this study employs the area-altitude balance ratio (AABR) developed by Osmaston (2005), which is widely regarded as the most accurate (Pellitero et al., 2015; Serrano et al., 2004). To determine the paleoELAs using this method, the ELA Calculation Tool for ArcGIS (Pellitero et al., 2015) was employed, with a balance ratio (BR, defined as the ratio between the accumulation area and the ablation area of a glacier) of 1.29 (Barr et al., 2022). For each calculated ELA, the glacier surface was reconstructed, considering the frontal moraines, till and glacial thresholds previously mapped and dated with cosmogenic isotopes. To determine the lateral extension of the Ésera glacier for each extent lateral moraines, till and in some cases very clear thresholds that avoid the expansion of very thin layer of ice previously mapped were considered (Vidaller et al., 2024a). The topographic reconstruction was conducted using a digital elevation model (DEM) with a resolution of 2 m per pixel, obtained from the Centro Nacional de Información Geográfica (CNIG). Given the favorable state of preservation of the moraines from the Little Ice Age (LIA), the LIA ELAs obtained using the AABR method were compared with the LIA ELAs obtained with the MELA method (Mean Equilibrium-Line Altitude; Serrano et al., 2012). This method considers that elevation of the upper limit of the frontolateral moraines coincides with the ELA of that moment, as was at this elevation where the glacier began to melt. Unfortunately, drawings or photographs of enough quality for the LIA period to better constrain the ELA are not available, so the best option was to use the glacial geomorphology features to reconstruct the glacier surface. The 2023 ELA was obtained from orthophotos captured during a drone survey conducted in September 2023. To determine the ELA in 2023, the line that separated the accumulation area (snow area) from the ablation area (ice area) was drawn manually from UAV photos, and we calculated as the mean elevation of this line for the three glaciers.

In order to ascertain the variation in temperature, it was necessary to consider the temperature lapse rate (TLR). The TLR was calculated considering the annual average of temperature (for yeas 2020 and 2022) from weather stations situated at Aneto peak (3404 m a.s.l.), Renclusa hut (2140 m a.s.l.), Besurta hut (1929 m a.s.l.) and Llanos del Hospital (1750 m a.s.l.), obtained from Posets-Maladeta Natural Park, AEMET and Clima y Nieve Pirineos database. Accordingly, the TLR in this valley is 0.525 °C (100 m)⁻¹ (see Fig. S2), which is consistent with the TLR considered in the Pyrenees (García-Ruiz et al., 2024). Using this gradient and the varying paleoELA elevations of each phase with the 2023 ELA, the temperature variation can be calculated using the following equation (Eq. 1):

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }T=\mathrm{TLR}\times ({\mathrm{ELA}}_{\mathrm{2023}}-{\mathrm{ELA}}_{\mathrm{phase}})\end{array}$$

The reconstruction of paleotemperature has an estimated error of ±0.45 °C (Barr et al., 2022), originated by the imprecision in the BR value. Also, the TLR present some inaccuracy, because we cannot assure that this gradient has been constant trough time as it is very dependent on atmospheric humidity. Even thought, it is the only available approximation to obtain a proxy of temperature variation between periods, as shown by Serrano-Cañadas and González-Trueba (2004) and Vidaller (2018).

4 Results

4.1 ¹⁰Be Exposure ages

The ¹⁰Be concentrations in the samples are presented in Table 1 for reference. The measured ¹⁰Be concentrations range from 0.13 to 2.98 × 10⁵ atoms g⁻¹ (Qz), with typical errors of approximately 4 %. However, some samples exhibit higher errors, reaching 8 %–13 %, which can be attributed to two main factors: firstly, the limited availability of quartz mass (in PDE-3), and secondly, the relatively low ¹⁰Be concentrations observed in PDE-8 and PDE-9. The ²⁶Al concentration was determined for one bedrock sample, PDE-3, and was 1.52 × 10⁶ atoms g⁻¹ (Qz), with an associated 20 % error (see Table S2). The ²⁶Al $/$ ¹⁰Be ratio for PDE-3 was calculated to be 8.0 ± 1.7, which was indistinguishable from the nominal production rate ratio of 6.75 within the stated uncertainty (Balco et al., 2008). It can therefore be concluded that this bedrock surface has most likely experienced a simple exposure history, without any prolonged burial period (> 2–300 ka).

The model minimum ¹⁰Be exposure ages, assuming no erosion and no snow shielding, range from 420 ± 36 years (PDE-8) and 390 ± 50 years (PDE-9) at the near Aneto peak to approximately 15 ka at the furthest point in the Ésera valley (PDE-13 and -14). It can be reasonably assumed that erosion rates of 3 mm ka⁻¹ would result in an increase in exposure ages of between 2 % and 4 % (see Table S3). Furthermore, the application of snow shielding would also result in an increase in the calculated ages of between 0 % and 6 % for the majority of samples, with the exception of PDE-2 and -3, for which the offsets were larger (8 % and 13 %, respectively).

Figure 3 Location of the cosmogenic samples (yellow dots) and ¹⁰Be exposure ages (pink squares) along the Ésera valley (geomorphological map modified from Vidaller et al., 2024a). Best ages (Table 1) are shown in ka. Hillshade obtained from a 5 m digital elevation model obtained from the Spanish National Geographic Institute.

In general, the obtained ages were grouped into four clusters (Fig. 4a), which corresponded to well-known periods during the last deglaciation (Clark et al., 2012; Oliva et al., 2019). The first pair of samples corresponded to the Oldest Dryas (OD) (PDE-13 and -14), the second group were eight samples dated for the Bølling–Allerød period (B–A) (PDE-1 to -7 and -12), the third group belonged to the Early Holocene (PDE-10 and -11), and the last pair matched the LIA (PDE-8 and -9). The samples from the furthest downstream location (PDE-13 and -14), dated at 16.26 ± 0.67 and 16.64 ± 0.62 ka and situated at Llanos del Hospital (Fig. 3), were retrieved from a moraine that was formed during the OD. Eight samples located between Pllan d'Están and Aiguallut paleolakes (PDE1-7 and -12) exhibit highly similar ages (between 12.54 ± 0.45 and 14.34 ± 0.52 ka), within the confines of their respective error margins (Fig. 4a). These samples, located all of them in a disperse till, represent the rapid deglaciation of the Ésera valley during the B–A period (14.6–12.9 ka). In the case of samples PDE-10 and -11, obtained from a frontal moraine, these samples were associated to an equilibrium glacial stage in the Early Holocene. Finally, the samples collected from a moraine at the highest elevation in this study (PDE-8 and -9) indicated an age of approximately 0.4 ka, which suggests the extent of the LIA for the Aneto glacier (Figs. 3 and 4 and Table 1).

Figure 4 (a) Distribution of dates according to the elevation of the samples. Clusters of ages are marked in yellow dash circles indicating the period associated to these samples. Right graph represents an enlargement of the left graph excluding LIA samples. (b) Distribution of dates according to the distance from the Aneto glacier front in 2023. Right graph represents an enlargement of the left graph excluding LIA samples, in order to focus on the retreatment and hold on of the deglaciation of the Ésera glacier.

The transition between the OD and the Allerød was marked by a retreat of the Ésera glacier of 2.4 km in ∼ 2.7 ka, equal to a deglaciation rate of 0.9 km ka⁻¹, indicating the onset of a progressively warmer period. Most of the samples from Pllan d'Están and Aiguallut paleolakes were dated the final stage of the Allerød period, suggesting the fastest deglaciation of this valley, ∼ 2.5 km in ∼ 1 ka (Fig. 4b). Between the end of the Allerød period (∼ 12.3 ka; PDE-12) and the Early Holocene (∼ 10.6 ka; PDE-10 and -11), the glacier retreated slower, at a rate of ∼ 1.0 km ka⁻¹ (retreated for 1.9 km in 1.7 ka), although the rate could have been even slower during the cold reversal corresponding to the Younger Dryas (YD, 12.9–11.7 ka; Palacios et al., 2016). Between the Early Holocene and the LIA, the glacier retreated more slowly at 0.2 km ka⁻¹ (from positions of the samples PDE-10 and -11, to PDE-8 and -9; Fig. 4b).

4.2 PaleoELAs

The reconstruction of the glacial phases during the deglaciation of the Ésera valley based on geomorphological evidence has enabled the definition of the paleoELAs and, to a reasonable approximation, the calculation of the temperature variation between those phases and the present day (2023). This has been achieved considering the ELA obtained with the AABR method for each moment and the TLR calculated with Eq. (1). The results are presented in Table 2 for reference. In order to evaluate the methodology employed, a comparison was made between the ELA determined from drone survey images in 2023 and the theoretical ELA obtained with the AABR method. This comparison revealed a high degree of similarity between the two calculations, with an average difference of only 9 m. Furthermore, a similar exercise can be conducted using the LIA moraines, comparing the paleoELA obtained with the MELA method and the paleoELA obtained with the AABR method. In this instance, when solely considering the largest glaciers (Maladeta, Aneto and Tempestades), which exhibit the most well-preserved moraines and thus provide the most robust calculations, the discrepancy between the theoretical paleoELA of the LIA and the elevation of the moraine did not exceed 50 m, thereby validating our calculations.

Table 2 Weighted average paleoELA (with the AABR method and BR 1.29) of each phase and variation of temperature (increase in all the cases) from each phase respect to year 2023. The uncertainty associate to the ELA corresponds with the standard deviation of the ELA of the glaciers of each period. In the same way, the uncertainty of the variation of temperature is measured as the variation of temperature associated to the standard deviation of the ELA for the different massifs and glaciers in the valley. These data do not include the ±0.45 °C error associated to the imprecision in the BR value (Barr et al., 2022) which is constant along the whole period. The phases mentioned in this table correspond with Fig. 5 phases. Each dated age corresponds to a specific moment during a climatic period, which does not imply that there were other situations during that period that were not recorded. The TLR was calculated considering the mean annual temperature of 2020 and 2022.

The lacustrine sedimentary sequence retrieved in Pllan d'Están paleolake (Fig. 1) established the position of the glacier front at 47 ka, coinciding with the initial deposition of lacustrine sediments in a proglacial lake (Vidaller et al., 2024a). A comparison of the paleoELA for that time and the current ELA allows the determination of a change of more than 3 °C on average (Table 2). Similarly, during the post-LGM deglaciation process of the Ésera valley, other glacial phases and the associated paleoELAs and temperature changes were calculated. At 16 ka, the Ésera glacier reached Llanos del Hospital, as indicated by PDE-13 and -14 isotopic cosmogenic dates, representing a glacier advance of 2 km down valley from the Pllan d'Están paleolake (Fig. 1). In accordance with the aforementioned evidence, the recovered sediments from Pllan d'Están indicated a deposition in a subglacial environment (Vidaller et al., 2024a). From the LIA to the present (year 2023), temperatures have risen by over 1 °C (based on paleoELAs), particularly in recent decades, leading to the rapid melting of the remaining glaciers.

5 Discussion

The combination of the new cosmogenic ¹⁰Be exposure ages and the reconstruction of paleotemperatures presented in this study, with the previous studies on the geomorphology of the valley (Vidaller et al., 2024a) and the Pllan d'Están sedimentary sequence (Fig. 5) (Vidaller et al., 2024b) allows us to reconstruct the evolutionary history of the Ésera glacier with a great detail (Fig. 6) and contribute to assess the regional deglaciation by comparing with previously published records (e.g. El Portalet lacustrine sequence; González-Sampériz et al., 2006).

Figure 5 Summary of the results obtained in the study of the sedimentary record of Pllan d'Están (modified from Vidaller et al, 2024b). Yellow lines in the age column refers to a hiatus in the sedimentation. Purple graphics show two geochemical analyses of the sedimentary record. Green graphs show the percentage of mesophytes and steppe taxa analysed from pollen grains founded in the sediments.

Figure 6 Evolution of the Ésera glacier since the PLGM. Yellow dots mark the cosmogenic and OSL dates used for the reconstruction of the glacier extensions. Panel (a) illustrates the cold phase of maximum extent of the Ésera valley 75 ka ago. Panel (b) shows a warmer phase during the MIS3 when Pllan d'Están paleolake was uncovered by ice for a time. Panel (c) represents a new cold period and an advance of the Ésera glacier, whose front was located at Llanos del Hospital during the OD. (d) At the onset of the Allerød a very fast deglaciation started. (e) At the end of the Allerød, the Ésera valley had split in several cirques' glaciers. Panel (f) shows the area covered by ice at the Early Holocene. Panel (g) illustrates the extension of the glaciers of the north face of the Maladeta massif during the LIA. Panel (h) represents the current situation where the glaciers are very small and withdrawn to the cirque's walls. See Fig. S3 in the Supplement to see the evolution of Ésera glacier from a bird's-eye view and the position of the moraines and till used for the delimitation of the glacier extension for each phase.

5.1 Pre-16 ka phase

The oldest evidence of glacial evolution at the headwaters of the Ésera valley is located at the base of the Pllan d'Están paleolake sedimentary sequence characterized as till deposit (Figs. 5, 6 and 7). This till has been dated by OSL to 74.9 ka ± 7.3 (Vidaller et al., 2024b), as illustrated in Figs. 6a and 7a. This sediment is likely to represent the Pyrenean Last Glacial Maximum (PLGM) phase, as described in previous studies for other Pyrenean valleys in spite is not yet clear the extension of the Ésera glacier at that time. Examples of the PLGM include OSL dating from glaciofluvial deposits dated in the Aragón valley (68 ± 7 ka; ∼ 900 m a.s.l.; García-Ruiz et al., 2013), the Ara valley (55 ± 4.5 ka; Sancho et al., 2018) and the Cinca valley (weighted average age from 3 OSL dates of 64 ± 11 ka; 790 m a.s.l.); terraces associated with this glacial phase in the Gállego valley (66 ± 4 and 74 ± 10 ka; Peña et al., 2004) or a clastic deposit in the Granito cave with an age of 71.8 ± 5.6 ka (Bartolomé et al., 2021). In Pyrenees, there are other evidences of glacial advances at the beginning of the Late Pleistocene (Calvet et al., 2011; Pallàs et al., 2010). Therefore, the last glacial maximum extent in the Pyrenees does not correspond in time with the global Last Glacial Maximum (LGM) considering a European context, which is estimated to have occurred between 23 and 18 ka implying the greatest advance in European glaciers, erasing the glaciological record of previous periods (Cutler et al., 2003; Toucanne et al., 2022).

Figure 7 Changes in sedimentation mechanism in Pllan d'Están overdeepening basin. Panel (a) represents the first phase, dated at 75 ka when the glacier was excavating the basin and some basal till was deposited. Panel (b) shows the sedimentation in Pllan d'Están paleolake when it was a proglacial lake during the relatively warm period of the MIS 3 (∼ 47 ka). The mm-thick lamination suggests the sedimentation under proglacial lake condition. Panel (c) represents the sedimentation in Pllan d'Están as a subglacial lake under the Ésera glacier. Banded sediments are associated to subglacial lakes. Both images of the sediment are obtained from the Pllan d'Están core, from different sections (Vidaller et al., 2024b).

Following its maximum extent, the Ésera glacier is thought to have undergone advances and retreats associated with the abrupt temperature oscillations that characterised the MIS 3 (60–27 ka; Dansgaard et al., 1993). However, there is no information regarding glacial activity in the valley until 46.7 ka, when lacustrine sedimentation commenced in Pllan d'Están (Figs. 5 and 6b) (Vidaller et al., 2024b). This age correlates with the onset of one of the most significant warm interstadials during the MIS 3, the Greenland Interstadial (GI) 12 (Rasmussen et al., 2014). This suggests that the inception of the major deglaciation phase in the Ésera valley may coincide with the GI-12. At that time, Pllan d'Están was a proglacial bedrock-dammed lake (Figs. 5, 6b, 7b), a type of proglacial lake common in large valleys (Otto, 2019). The glacier front was situated at an elevation of 1840 m a.s.l. In general, sedimentation in proglacial lakes is primarily influenced by the proximity to the glacier, resulting in the formation of varves and rhythmites (Carrivick and Tweed, 2013). In Pllan d'Están, the initial lacustrine sediments were laminated and contained some carbonate eroded from the limestones in close proximity to the lake (Vidaller et al., 2024b). A similar pattern has been identified in other high-altitude lakes, including Marboré lake in the Central Pyrenees (Leunda et al., 2017) and Enol lake in the Cantabrian mountain range (Moreno et al., 2010). In these cases, the carbonate in the sediments has been preserved during relatively warm periods, when the dissolution of carbonates was prevented. It is noteworthy that the paleoELA at this age was estimated to be 2517 m a.s.l. (Table 2), representing a temperature offset of 3.1 °C colder than the present. This temperature variation is consistent with marine sediment cores from the eastern margin of the Iberian Peninsula and the Western Mediterranean, which also indicate that GI-12 was the warmest interstadial (although colder than the present) within MIS 3 (Cacho et al., 1999; Martrat et al., 2007).

The carbonate content in Pllan d'Están sediments exhibited a precipitous decline following 34.8 ka BP marking a colder period (Fig. 5). The concomitant shifts in sedimentary facies and the low sedimentation rate are indicative of a subglacial lake environment (Vidaller et al., 2024b). This interpretation of very cold conditions is corroborated by palynological findings in other lacustrine sequences, including El Portalet (González-Sampériz et al., 2006) and El Cañizar de Villarquemado (González-Sampériz et al., 2020) for the same time period. Furthermore, the existence of moraines in the Gállego valley at an elevation of 830 m a.s.l., dated using OSL to 36 ka (Lewis et al., 2009), also indicates a period of low temperatures and glacial expansion.

The LGM in a European context is linked to the most significant expansion of glaciers across northwest Europe, where glaciers reached their maximum extent during the Late Würmian period (MIS 2: 35–14 ka; Florineth and Schlüchter, 2000; Ivy-Ochs et al., 2008; Buoncristiani and Campy, 2011). While not ubiquitous across Europe, the advance was consistent in its magnitude. For instance, in the Maritime Alps, the estimated ELAs are as much as 450 m higher than those observed in the northern Apennines or the Corsican mountains. This suggests that the Italian side of the Maritime Alps experienced a relatively arid climate during the LGM (Rettig et al., 2024). In spite the general importance of the LGM period in many European glaciers, this time period presents a more complex picture in the Pyrenees. In the headwaters of the Ésera valley, for instance, there are no published dates that correspond to this period (23–18 ka), but surely Pllan d'Están was covered by ice, still being a subglacial lake (Fig. 7c). We then assume that the remains of a possible advance of the glacier during the LGM were eroded by following glacier phases (e.g. during the OD). In the Eastern Pyrenees, however, some LGM cosmogenic exposure ages from moraines at an elevation of 2183 m a.s.l. and evidence of the cold stage from lacustrine deposits have been reported. The existence of such sequences has been documented in the literature (Delmas, 2015; Delmas et al., 2008, 2011, 2022; Pallàs et al., 2006, 2010; Rodes et al., 2008). The most supported explanation for this difference in the LGM glacier expression from west to east in the Pyrenees is related to the available humidity. The distribution of temperature and precipitation is mainly controlled by the elevation and the proximity to the Atlantic Ocean or the Mediterranean Sea (Cuadrat et al., 2007; García-Ruiz et al., 2015). The transition between the two climates is approximately in the Central Pyrenees, around the Maladeta massif.

5.2 Phase Llanos del Hospital (16 ka)

A significant indicator of the evolution of the Ésera glacier is the moraine dated at 16 ka in Llanos del Hospital at an elevation of 1750 m a.s.l. (PDE-13 and -14; Table 1 and Fig. 6). These ages coincide with a well-documented glacial advance in Europe, the OD (approximately 17.5–15 ka), when Alpine glacier fronts reached 1400 m a.s.l. (Darnault et al., 2012; Ivy-Ochs et al., 2006; Kerschner and Ivy-Ochs, 2008). The OD in the Pyrenees has been dated using moraines from the Eastern Pyrenean valleys (Andrés et al., 2018; Pallàs et al., 2010; Reixach et al., 2021), in particular the Carlit massif (15 ka; Delmas et al., 2008), the Màniga-La Feixa platform (15.5 ka; Pallàs et al., 2010), the Arriège valley (Delmas, 2015; Delmas et al., 2011, 2012), the Noguera-Ribagorzana valley (Pallàs et al., 2006) and the Caldarés valley (dates from 17.6 to 15.9 ka at an elevation of 2200–1400 m a.s.l.; Palacios et al., 2017b) have also yielded similar results. Furthermore, the OD is well characterised in Pyrenean lakes and speleothems, indicating a sharp decrease in temperature (Bernal-Wormull et al., 2021) and an increase in aridity (González-Sampériz et al., 2006; Morellón et al., 2009).

Reixach (2022) measured ¹⁰Be exposure ages from nearby sites at Llanos del Hospital. Their ages (18–10 ka) are scattered and have relatively large errors (6 %–14 %), and therefore the timing of glacial advance and retreat was not well constrained. Our new exposure ages (16.6 ± 0.6 and 16.3 ± 0.7 ka; PDE-13 and -14 in Table 1), which are concordant between the two samples and more precise, largely agree with the previous data when their ages are recalculated using the same parameters we used in this study (e.g. ¹⁰Be production rate, geographical scaling scheme, snow density, post-depositional denudation rate).

The continuous sedimentary lacustrine sequence of Pllan d'Están and the location of the 16 ka moraine at a lower elevation in the valley suggest a subglacial lake environment between 34–13 ka (Figs. 5, 6c and 7c) where sediments are delivered into the lake through a complex subglacial hydrological system (Smith et al., 2018). Although it is very difficult to distinguish between proglacial and subglacial sediments, the sediments of Pllan d'Están have different facies (Fig. 7), which help to distinguish the two subenvironments. Proglacial sediments are usually characterised by the alternation of clay (summer months) and silt (winter months) laminations forming varves (Palmer et al., 2008; Ringberg and Erlström, 1999; Smith and Ashley, 1985). The rhythmic nature of the lamination is also marked by different grain sizes, indicating strong seasonality in a glacial environment (Leonard and Reasoner, 1999; Ohlendorf et al., 2003). The coarse, angular silt laminae are deposited during the melting season, whereas the fine silt laminae are deposited during the ice-covered season when fine particles settle by suspension (Carrivick and Tweed, 2013). In subglacial lakes, sediments are usually homogenised towards the top of the sequence, with sand bands and laminae (Livingstone et al., 2015), as occurred in the Pllan d'Están sedimentary sequence after 34 ka BP (Vidaller et al., 2024b) (Fig. 7b and c). Furthermore, the sedimentation rate is lower than in the other lacustrine sequences of the Pyrenees, such as El Portalet (González-Sampériz et al., 2006) or the Tramacastilla lakes (García-Ruiz et al., 2003; Jalut et al., 1982; Montserrat-Martí, 1992), which supports the subglacial lake condition at Pllan d'Están during ∼ 34–13 ka BP. Therefore, sedimentary and chronological evidences indicate that Pllan d'Están was covered by ice at the time of the global LGM and also the OD (Figs. 5 and 6c and d).

The paleoELA at 16 ka was at 2410 m a.s.l., which implies that the average temperature was 3.6 °C lower than the present. Also, the paleoELA calculated is consistent with similar studied in the western Alps, where the paleoELA was at 2428 m a.s.l. (Serra et al., 2022). Similarly, studies of marine sediments from the Iberian margin and the Western Mediterranean indicate a comparable variation in temperature since the Heinrich Event 1 (Cacho et al., 2001; Martrat et al., 2007). Additionally, other terrestrial records from the Iberian Peninsula indicate the occurrence of a particularly cold period at 16 ka (Bernal-Wormull et al., 2021; Pérez-Mejías et al., 2021).

5.3 Phase Pllan d'Están (13.7–12.8 ka)

The B–A period was characterized at the scale of the northern hemisphere by an increase in temperature and a change towards more humid conditions (Clark et al., 2012). In Iberia, the B–A is observed in many terrestrial (González-Sampériz et al., 2006; Bernal-Wormull et al., 2021; Vidaller et al., 2024b) and marine records (Fletcher et al., 2010) where a trend to higher temperatures is evidenced. In the Ésera valley. The majority of the cosmogenic samples included in this study correspond to the B–A period (14.6–12.9 ka) and were situated between Pllan d'Están basin and Aiguallut ponor (PDE1-7 and -12; Fig. 6). The cosmogenic dates indicate that the Ésera glacier commenced its retreat following the OD, approximately 16 ka. Therefore, during the Allerød period, the glacier front retreated from 1750 m a.s.l. (Llanos del Hospital) to 1995 m a.s.l. (Pllan d'Están), representing a rapid retreat of almost 4 km (at a rate of 2.3 km ka⁻¹). The absence of moraines in this area, with only sporadic till present, lends further credence to the hypothesis that the glacier retreated rapidly during this period, never attaining a state of equilibrium. This improvement in climatic conditions during the B–A period, was reflected by an increase in the carbon content of the sediment and by a change from steppe species to mesophytes as show the pollen record (Fig. 5).

Similarly, other valleys in the Pyrenees exhibited comparable rapid retreat during the B–A, including the Bacivèr cirque (Oliva et al., 2021) and the Gállego valley. In the latter, the glacial front retreated from 1500 to 2200 m a.s.l. between 14.6–11.7 ka (Palacios et al., 2015a). Similarly, glaciers situated in the Eastern Pyrenees (Delmas et al., 2008, 2011, 2023; Oliva et al., 2021) and the Ariège valley (Jomelli et al., 2020; Reixach et al., 2021) retreated back to cirques during the B–A period. The shrinkage of the glaciers resulted in the formation of rock glaciers and debris covered glaciers, below an elevation of 2800 m a.s.l. in numerous areas within the Pyrenees (Andrés et al., 2018; Oliva et al., 2021; Palacios et al., 2017a). Therefore, the B–A period was characterised by a rapid process of deglaciation, whereby glaciers in valley bottoms undergo a transformation into cirque glaciers. In the Ésera valley, this phenomenon was exemplified by the Maladeta-Alba, Aneto-Tempestades and Escaletes glaciers (Fig. 6d and e).

In terms of temperature, at the end of the Allerød period, the paleoELA in the headwaters of the Ésera valley was recorded at 2645 m a.s.l., exhibiting a 1.3 °C increase compared to the previous OD. This value is considerably lower than the estimated increase of 3–5 °C derived from pollen records across Europe (Clark et al., 2012) or in marine cores surrounding the Iberian Peninsula (an increase of 5 °C from Heinrich Event 1 to the Allerød period; Cacho et al., 2001; Martrat et al., 2007; Bernal-Wormull et al., 2025). PaleoELAs have been demonstrated to serve as effective proxies in the determination of temperature variations during periods for which instrumental records are unavailable. However, it should be noted that the method is not without its limitations. In certain instances, such as the present case, the temperature variation values obtained through the utilisation of paleoELAs do not exactly correspond with those obtained through the application of alternative proxies. Among other limitations in the method, we are aware of the likely changes in the temperature lapse rate with elevation and with the season throughout the different considered periods.

5.4 Phase Salterillo-Barrancs (11 ka)

The ¹⁰Be ages from PDE-10 and -11 (11.1 ± 0.6 and 10.8 ± 0.4 ka, respectively; Table 1) indicate that the Ésera glacier front was situated at 2046 m a.s.l. during the Early Holocene (Fig. 6f). These ages are somewhat older than the few preceding cosmogenic dates from analogous locations (∼ 4–8 ka with one outlier at ∼ 16 ka; Reixach, 2022). Similar ages have been reported from other regions within the Pyrenees. For example, in the La Cerdanya mountains (Southeastern Pyrenees), a rock glacier boulder was dated to 10.0 ± 0.4 ka (³⁶Cl exposure ages) at an elevation of 2490 m a.s.l. (Palacios et al., 2015b). Similarly, in the Noguera-Ribagorzana valley at an elevation of 1721 m a.s.l., some cosmogenic ages have been reported, with the most recent being 10.4 ka (Pallàs et al., 2006). In the Carlit massif (Delmas et al., 2008), dates between 11–10.7 ka were obtained at an elevation of 2170–2180 m a.s.l. Thus, some polished bedrocks at an elevation of 2549–2719 m a.s.l. in the Gállego valley (Central Pyrenees), were dated with ¹⁰Be and ³⁶Cl exposure ages, resulting in a date of 10.6 ± 1.3 ka (Palacios et al., 2017a). The discrepancy in altitude between these samples may be attributed to their location in a western position, which is subject to a greater influence from an Atlantic climate.

In terms of temperature, the paleoELA of this phase indicates a difference of 1.7 °C relative to the present (Table 2), and 0.7 °C higher than the end of the B–A period. At a global scale, an increase of 1.3 °C was recorded since the Early Holocene to current conditions (Shakun et al., 2012; Snyder, 2016). In the Eastern Swiss Alps, paleoELAs calculation showed an increase of 1.8 °C (Joerin et al., 2008). This value contrasts again with the temperature change defined between the B–A and the onset of the Holocene in the Alboran sea, a change in sea surface temperature (SST) of ∼ 3 °C (Cacho et al., 2001). The one study based on chironomids in the Pyrenees reveals a temperature difference of 1 °C between the Early Holocene and the present (Tarrats et al., 2018), more coherent with our estimations. Other studies based on chironomids too, but in the Australian and Swiss Alps, calculate an increase of 4 °C since the Early Holocene to nowadays (mean July temperature; Ilyashuk et al., 2011) in the Austrian case and 1.7 °C in the swiss case (Samartin et al., 2012).

5.5 Phase Little Ice Age (0.4 ka): the last advance of the glaciers

The LIA represents the final cold period documented in numerous mountainous regions across the globe, occurring during the 14th and 19th centuries (García-Ruiz et al., 2020; Grove, 2004; Oliva et al., 2018a; Solomina et al., 2016). This cold pulse was characterised by low summer temperatures and snowy winters (Matthews and Briffa, 2005), implying a global cooling of 1–2 °C lower than the present (Dyurgerov and Meier, 2000; Grove, 2004), with a likely larger temperature change in mountainous regions such as the Pyrenees. Some studies based on paleoELAs have determined that temperatures have increased by 2.5 °C in the Eriste massif (Central Pyrenees; Vidaller, 2018), 2 °C in the Tendeñera and Sabocos mountain range (López-Moreno, 2000) and 0.9 °C for the Posets massif (Serrano and Martín-Moreno, 2018), which is very similar to our own calculation of +1.1 °C obtained from the Maladeta massif (Table 2). As indicated by data from the Spanish National Meteorology Agency (AEMET) and the Catalonia Meteorological Service (SMC), the temperature increase since the conclusion of the LIA in the Pyrenees region has been 1.1 °C (Pérez-Zanón et al., 2017).

Since the LIA maximum, dated in this study with approx. 400 years (Fig. 6g), temperatures have increased almost continuously, and the glaciers have shrunk rapidly globally (Oliva et al., 2018b; Zemp et al., 2015). During the LIA there were 52 glaciers in the Pyrenees, that covered 2060 ha of glacier surface (Rico et al., 2017). Considering the case of Aneto glacier, the largest glacier in the Pyrenees, it has lost 64.7 % of its area in the last 41 years (period 1981–2022), and its ice thickness has decreased by an average of 30.5 m (Vidaller et al., 2023). In addition, the occurrence of extremely hot and dry years, such as those observed in 2022 and 2023 summers, has accelerated the melting processes leading to a drastic degradation of the glacier and posing a high risk for its survival (Vidaller et al., 2023). In fact, at the current ELA, located at 3098.6 m a.s.l., very close to the walls of the cirque, snow/ice is absent in some areas indicating the demise of current glaciers and their accelerated melting.

6 Conclusions

The ice fluctuation of the headwaters of the Ésera valley since the PLGM was a complex process, comprising both advances and rapid retreats. The chronological sequence of deglaciation was determined through the analysis of dates of cosmogenic ¹⁰Be from granitic boulders and polished bedrock along the headwaters of the Ésera valley. The excellent reproducibility of our data at each site, for example at Pllan d'Están (PDE-1, -2), in the vicinity of the Besurta hut (PDE-5, -6), at the LIA moraine (PDE-8, -9), the upper slope of Aiguallut (PDE-10, -11), or Llanos de Hospital (PDE-13, -14), makes a robust chronology of glacial phases in the Ésera valley. Furthermore, paleoELA calculations were conducted for various phases in order to estimate associated temperature changes.

Following the maximum extent of the Ésera glacier during the PLGM at approximately 75 ka, advances and retreats associated with climate oscillations during MIS 3 are likely to have occurred (60–27 ka). At the beginning of GI-12, a significant deglaciation occurred in the Ésera valley, marked by the formation of a proglacial lake in Pllan d'Están at approximately 47 ka. A subsequent glacial advance resulted in the lake being covered by the glacier once more time, and Pllan d'Están became a subglacial lake for an extended time period while sedimentation continued through the input of water and sediment via glacier crevasses. This subglacial lake condition was identified for the first time in the Pyrenees by combining the new ¹⁰Be exposure ages presented in this study with those from previous sedimentological studies of the paleolake.

An important phase of the evolution of the Ésera glacier is marked by a moraine dated at 16 ka in Llanos del Hospital, which represents the OD advance. This occurred when the temperature was 3.6 °C colder than the present day. During the B–A interstadial, the climate exhibited a large warming trend, particularly rapid during the Allerød. Our ELA analysis indicates an increase of approximately 0.8 °C in temperature (from the onset to the end of the Allerød) over a period of approximately 1000 years, accompanied by a rapid retreat of the Ésera glacier from Pllan d'Están to near Aiguallut ponor (2.5 km). This study offers a contrasting perspective to the majority of European paleoclimate records, which indicate an abrupt transition at the onset of the Bølling period. The Early Holocene was characterised by an increase in temperature and the transformation of the Ésera glacier from a valley glacier into smaller cirque glaciers. The LIA (0.4 ka) was the final cold period recorded in these mountains and resulted in the formation of the largest and most well-preserved moraines in the valley.