Rapid deglaciation during the Bølling-Allerød Interstadial in the Central Pyrenees and associated glacial and periglacial landforms

a Department of Geography, Universitat de Barcelona, Catalonia, Spain b Centre for Geographical Studies, IGOT, Universidade de Lisboa, Lisbon, Portugal c Department of Geography, Complutense University of Madrid, Madrid, Spain d Aix-Marseille Université, CNRS, IRD, INRAE, Coll. France, UM 34 CEREGE, Aix-en-Provence, France e Department of Geography & Centre for Northern Studies, Université Laval, Quebec, Canada


Introduction
Termination-1 (T-1), the period spanning from the end of the Last Glacial Maximum (LGM, 19-20 ka;Clark et al., 2009) to the onset of the Holocene (11.7 ka; Denton et al., 2014), saw a massive worldwide glacial retreat that favored a large-scale reorganisation of oceanic and atmospheric circulation patterns, global sea level rise, redefinition of coastlines, shifts in land cover and ecosystems, and changes in greenhouse gas concentrations. An accurate comprehension of the spatial and temporal patterns of environmental change that occurred during the last major deglacial period, particularly glacial oscillations, can provide insights into rapid landscape readjustment, which allows the significance of recent trends due to warming climates to be assessed in a longer-term context (Oliva et al., 2019;Oliva et al., 2021). A better understanding of the sensitivity of glaciers in mid-latitude mountains to rapidly changing past climates can thus help to assess the magnitude of future changes and the fate of mountain glaciers in these regions.
Whereas the long-term, global-scale glacial retreat that occurred during T-1 was favored by greenhouse gas increases, regional glacial advances and retreats followed forcings at smaller scales (Denton et al., 2014). In the Northern Hemisphere, ice core records from Greenland reveal a sequence of alternating colder and warmer periods during T-1: the cold period, the Oldest Dryas (OD; 17.5-14.6 ka), was followed by the much warmer Bølling-Allerød (B-A) Interstadial (14.6-12.9 ka) and a subsequent return to colder conditions during the Younger Dryas (YD; 12.9-11.7 ka) (Rasmussen et al., 2014). In mid-latitude regions, such as the Iberian Peninsula, glaciers shrank in response to the long-term warming that was recorded during T-1, although colder millennial-scale phases favored re-expansion and warmer periods triggered accelerated shrinking (Oliva et al., 2019). While local maximum ice extents during the last glacial cycle occurred asynchronously in different mountain ranges in Iberia (Oliva et al., 2019), glacial oscillations during T-1 followed very similar patterns in response to changing climate in the North Atlantic region (Buizert et al., 2018;Rea et al., 2020). For areas where chronological data are available, such as the Central Range (Carrasco et al., 2015;Palacios et al., 2012), Pyrenees Palacios et al., 2017bPalacios et al., , 2017aPalacios et al., , 2015b, Cantabrian Mountains (Rodríguez-Rodríguez et al., 2017), and Sierra Nevada (Gómez-Ortiz et al., 2012;Palacios et al., 2016), glaciers generally advanced during the OD and YD, and retreated during the B-A and after the YD during the Early Holocene.
In the Pyrenees, where the present work focuses, glacial records in several valleys (Gállego, Ésera, Noguera Ribargoçana, Ariège) have shown evidence of two periods with glacial expansion during T-1 (Crest et al., 2017;Jomelli et al., 2020;Palacios et al., 2017bPalacios et al., , 2015aPallàs et al., 2006). The first phase of glacial advance after the massive retreat of the LGM took place during the OD, with ice tongues up to 15 km long reaching the mountain fronts. The second phase, in the late Pleistocene, occurred during the YD with the presence of small glaciers up to 4 km long which developed from north-facing cirques of the highest massifs, above 2200 m (García-Ruiz et al., 2016). In between, the B-A saw a rapid glacial retreat, with the ice disappearing from most of the highest massifs by 15-14 ka and the formation of rock glaciers and debris-covered glaciers in many cirques in response to intense paraglacial readjustment, mainly in relatively low-altitude cirques below 2800 m Palacios et al., 2017b). By the end of the YD at 11-10.5 ka, almost all YD glaciers had disappeared and a new generation of rock glaciers formed in the recently deglaciated cirques .
Despite recent advances in our understanding of climatic and environmental consequences during T-1, the spatial and temporal patterns of glacial response in some European mountains is still poorly known, particularly in the Mediterranean region. In this sense, the Pyrenees, located in the transitional area between Atlantic and Mediterranean climatic influence, and between the southern and northern European mountains, constitute a mountain range of great glacial and climatic relevance. However, the chronology of glacial oscillations in several valleys of the Central Pyrenees is still uncertain, particularly for T-1. In addition, the origin of the great variety of glacial and periglacial landforms, especially in relatively low-altitude cirques, is unknown, as is the chronology of their formation. Data from this area and this period is therefore needed to shed light on the prevailing paleoclimatic conditions, as well as the atmospheric configuration driving glacial oscillations during T-1. To this end, we had the following specific objectives: -To provide new data on glacial oscillations during T-1 for the Upper Garonne Basin, where absolute deglaciation ages are still lacking. -To compare results from the Upper Garonne Basin with the timing of deglaciation as well as the age of formation and stabilization of the different glacial and periglacial landforms that exist in many of the cirques at lower altitudes in the Pyrenees.
-To compare glacial evolution in this mountain range with that which occurred in the Iberian mountains as well as other southern European ranges, and contrast their spatio-temporal patterns. This is needed to frame the glacial response within the paleoclimatic evolution of Europe during deglaciation, as inferred from natural archives and climate models.

Study area
The Pyrenees are the largest mountain range in the Iberian Peninsula, spanning a 400 km W-E transect. The central part of the range contains the highest massifs, with peaks exceeding 3000 m a.s.l., such as the Maladeta (3404 m), Monte Perdido (3355 m) and Posets (3371 m). This research focuses on the Bacivèr Cirque, located in the upper Garonne valley ( Fig. 1), at latitudes 42°40′ N-42°42′ N and longitudes 0°57′ E-1°00′ E.
The relief is structured by the U-shaped Garonne valley that drains towards the N-NW, receiving inflow from tributary rivers that drain adjacent glacial valleys with headwaters in peaks ranging from 2800 to 3000 m. Glacial cirques in the Upper Garonne valley are predominantly NE-exposed and the altitude of their floors ranges between 2200 and 2400 m . The floor of the Bacivèr Cirque, where this study focuses, is located at similar altitudes but faces W; it extends over 10 km 2 with the highest altitudes exceeding 2600 m at the Marimanya (2675 m; 42°42′33″N 1°00′44″E), La Llança (2658 m) and Bacivèr (2642 m) peaks and the lowest at 1850-1900 m at the Beret Plateau ( Fig. 1). This area represents the hydrological divide between rivers draining to the Atlantic via the Garonne River and those flowing to the Mediterranean Sea via the Noguera Pallaresa.
At present, the mean air annual temperature (MAAT) at the nearby Bonaigua station (2266 m) is 3°C whereas the annual precipitation totals 1227 mm, mostly in the form of snow that falls during the cold months of the year. Snow on the ground generally persists for 7-8 months of the year. The treeline lies between 2200 and 2300 m, where Pinus nigra is replaced at higher elevations by alpine meadows on the cirque floor and barren rocky terrain elsewhere. The lithology of the area is mainly composed of Carboniferous granodiorites and granites, but there are also Silurian marbles at the lower margins of the Bacivèr Cirque and limestones intercalated with Devonian slates on the Beret Plateau (Institut Cartogràfic i Geològic de Catalunya, 2017).
The landscape of the Bacivèr Cirque includes a wide range of glacial and periglacial landforms that are inherited from the last Pleistocene glacial cycle and the subsequent deglaciation (Fernandes et al., 2017). The W-NW aspect of this compound cirque and the high elevation of its floor (2200-2400 m) resulted in abundant snowfall accumulation and its subsequent transformation into ice during the cold Pleistocene phases. However, the chronology of the local maximum ice extent of the last glacial cycle is not yet known, although recent studies encompassing the entire Garonne paleoglacier show evidence that the entire cirque was largely covered by ice during that phase (Fernandes et al., 2017). This is confirmed by the existence of glacially polished sur-faces~200 m above the Bacivèr Cirque floor. Post-LGM warming favored the shrinking of that glacier, which receded and lost thickness during the last deglaciation. This is confirmed on the Beret Plateau, where glacial striae on polished surfaces show different directions which suggest that the glacier descending from the Bacivèr Cirque, among others, diverged either to the main Garonne Glacier (SW) or the Noguera Pallaresa Glacier (N). Within the Bacivèr Cirque, there is a variety of erosive and depositional landforms, including moraine complexes indicative of multiple readvances within the final stages of the long-term deglaciation. As observed in other neighboring valleys, the melting of the last glaciers favored the development of rock glaciers during the paraglacial stage (Knight, 2019;Knight et al., 2019), although they are inactive under present-day climate conditions . The southern fringe of the Bacivèr Cirque is included in the Baqueira-Beret ski resort domain, which has altered some slopes to expand winter sports facilities.

Methodology
In order to reconstruct past glacial oscillations since the onset of deglaciation in the Bacivèr Cirque, we used an integrated geomorphological and geochronological approach. Field work was conducted in June 2016, when the absence of snow cover enabled the identification of different geomorphic features and collection of samples for Cosmic-Ray Exposure (CRE) dating.

Geomorphological mapping
We produced a geomorphological map at 1:20000 scale by using an ArcGIS 10.7 (ESRI) database. The map was based on: (i) stereoscopic photo-interpretation though out the Iberpix 4 online anaglyph map viewer (https://www.ign.es/iberpix2/visor/); (ii) visual inspection of satellite imagery from Google Earth; and (iii) the combination of orthophotomaps (0.25 m cellsize) and high-resolution LIDAR digital elevation models (density points of 0.5-2 m −2 ) obtained from the 'Institut Cartogràfic i Geològic de Catalunya' (http://www.icc.cat/ appdownloads). Special attention was paid to glacial landforms and related features, which were outlined and symbolized according to Joly (1997). The generated geomorphological map was then validated in the field, as we surveyed the entire area with a focus on glacial and periglacial landforms.

Field strategy and sampling
We collected a total of 17 samples for CRE dating using a hammer and a chisel. We targeted boulders belonging to moraines, a debriscovered glacier and glacially polished outcrops of granites and granodiorites. We aimed to ensure the optimal exposure of the sampling sites to the cosmic-ray flux, and thus selected flat-topped and gentle surfaces (<20°) of rock outcrops rather than steep slopes or sharp crests. The thickness of the extracted samples ranged from 2 to 4.5 cm (Table 1). To account for any shielding due to the surrounding topography, the topographic shielding factor was calculated for all sampling sites by means of the ArcGIS toolbox devised by Li (2018) that implements well-known routines explained in Dunne et al. (1999), and only needs a point shapefile of the sampling sites, including the strike and the dip of the sampled surfaces, and a digital elevation model (DEM).

Laboratory procedures and exposure age calculation
Before the chemical processing of the samples, they were crushed and sieved to the 0.25-0.8 mm fraction at the Physical Geography Laboratory (Universidad Complutense de Madrid, Spain). Thereafter, we treated the rock samples at the Laboratoire National des Nucléides Cosmogéniques (LN 2 C) of the Centre Européen de Recherche et d'Enseignement des Géosciences de l'Environnement (CEREGE; Aixen-Provence, France). In accordance with the quartz-rich lithology of the samples, they were processed for the extraction of the in situ produced cosmogenic nuclide 10 Be.
In order to remove magnetic minerals, we conducted a magnetic separation through a "Frantz LB-1" separator. Once the non-magnetic fraction was isolated, it underwent several rounds of chemical attacks with a concentrated mixture of hydrochloric (1/3 HCl) and hexafluorosilicic (2/3 H 2 SiF 6 ) acids to dissolve and discard the nonquartz minerals. Subsequently, four successive partial dissolutions of the remaining minerals with concentrated hydrofluoric acid (HF) helped dissolve the remaining impurities (e.g. non-dissolved feldspar minerals) and removed atmospheric 10 Be. As a result, samples yielded from 8 to 22 g of purified quartz (Table 2). Just before the total dissolution of quartz, 150 μL of an in-house manufactured (from a phenakite crystal) 9 Be carrier solution (spike, concentration: 3025 ± 9 μg g −1 ; Merchel et al., 2008) were added to the samples. The purified quartz was subsequently dissolved by acid leaching with 48% concentrated HF (3.6 mL per g of quartz + 30 mL in excess). Following the total dissolution, we evaporated the resulting solutions until dryness, and recovered samples with HCl (7.1 M). The Be samples were then precipitated at PH = 9 to beryllium hydroxide (Be(OH) 2 ) by means of ammonia (NH 3 ), and separated from other elements in resin columns: an Dowex 1 × 8 anionic exchange column to remove elements such as Fe, Mn and Ti, and a Dowex 50WX8 cationic exchange column to discard B and recover Be (Merchel and Herpers, 1999). The final eluted Be was precipitated again, and the Be precipitate was dried and oxidized to BeO at 700°C. Finally, the targets for accelerator mass spectrometer (AMS) measurements were prepared by mixing the niobium powder with the BeO keeping an approximate 1:1 proportion and pressing the mixture into copper cathodes.
The targets were analysed at the Accelerator pour les Sciences de la Terre, Environnement et Risques (ASTER) national AMS facility at CEREGE in order to measure the 10 Be/ 9 Be ratio from which the 10 Be concentration was later inferred ( Table 2). The AMS measurements were calibrated against the in-house standard STD-11 with an assigned 10 Be/ 9 Be ratio of (1.191 ± 0.013) × 10 −11 (Braucher et al., 2015). The analytical 1σ uncertainties include uncertainties in the AMS counting statistics and an external 0.5% AMS error (Arnold et al., 2010) and the uncertainty related to the chemical blank correction. The 10 Be half-life considered was (1.387 ± 0.0012) × 10 6 years (Chmeleff et al., 2010;Korschinek et al., 2010).
We calculated 10 Be exposure ages by using the CREp online calculator (Martin et al., 2017; available online at: http://crep.crpg.cnrs-nancy. fr/#/), where we selected the following settings: LSD (Lifton-Sato-Dunai) elevation/latitude scaling scheme , ERA40 atmospheric model (Uppala et al., 2005) and geomagnetic database based on the LSD framework . These settings yielded a world-wide mean 10 Be production rate at sea level high latitude (SLHL) of 3.98 ± 0.22 atoms g −1 yr −1 . Exposure age results of the samples are shown in Table 2, together with their full 1σ uncertainties (derived from analytical and production rate uncertainties) and their analytical uncertainties only. The uncertainties discussed throughout the text are given with their analytical uncertainties only, for internal comparison. Due to the high analytical uncertainty, arising from low AMS currents and counting statistics, sample A-17 was discarded from the discussion. In order to evaluate the impact of potential erosion on the exposure ages, we assumed a steady erosion rate (1 mm ka −1 ) with a conservative maximum value according to André (2002). The impact of snow cover was estimated by extrapolating the current snow duration in the area, with an annual average of 7.5 months at 2200-2300 m and a mean thickness of 100 cm based on the data series of Bonaigua station ( Fig. 1) available since 1997-1998 (Servei Meteorològic de Catalunya; http://www.igc.cat/web/ca/allaus_gruix_neu_v2.php?e=bonaigua&t= totes) (Table 3). We applied eq. 3.76 in Gosse and Phillips (2001) to calculate the snow correction factor (Gosse and Phillips, 2001). Erosion and snow corrections as a whole resulted in ages older by~9%, with a minor contribution from the erosion correction (<2%). However, throughout the text we use the uncorrected ages in order to enable comparison with other areas, and considering that past snow cover thickness and duration are unknown. It must have changed significantly since the cirque's deglaciation, with an alternation of colder/warmer periods associated with changing moisture conditions.

Paleoglacier reconstruction and Equilibrium-Line Altitude (ELA) calculation
Paleoglaciers were reconstructed to model their spatial extent and paleotopography during the considered time periods. Threedimensional glacier reconstruction was produced using the 'GLaRe' Table 2 AMS analytical data and calculated exposure ages. 10 Be/ 9 Be ratios were inferred from measurements at the ASTER AMS facility. Individual ages are shown with their full uncertainties (including analytical AMS uncertainty and production rate uncertainty) and analytical uncertainty only within brackets. Arithmetic mean ages are given with their full uncertainties (including standard deviation and production rate uncertainty) and standard deviations only in brackets. Ages in grey italics correspond to potential outliers and thus are rejected and excluded from the interpretation and discussion. 10  ArcGIS toolbox devised by Pellitero et al. (2016). It estimates past ice thickness along a flowline by applying the perfect-plasticity physicalbased numerical model of Van der Veen (1999) following the calculation routines later proposed by Benn and Hulton (2010). The toolbox only requires a flowline, a tentative paleoglacier geometry (approached as a basin whose boundaries are defined by the position of lateral and/or frontal moraines) and a digital elevation model. Ice thickness was modelled by using an average shear stress of 100 kPa (Paterson, 1994;Benn and Hulton, 2010). Ice thicknesses were corrected based on shape factors (F-factor) obtained from a number of representative cross-sections in order to reduce the error of modelled values to <10% . From these procedures, we produced DEMs of the paleoglaciers during different stages. Later, from those DEMs, ELAs were calculated by using the automatic toolbox developed by Pellitero et al. (2015). The selected methods were the Accumulation Area Ratio (Porter, 1975) and the Area Altitude Balance Ratio (AABR; Osmaston, 2005). When applying the AABR, we considered two BR, namely: 1.75 ± 0.71 (global) and 1.9 ± 0.81 for mid-latitude maritime glaciers (Rea, 2009). We selected the latter, given the location of the study area in the mid-latitudes and the local influence of the Atlantic air masses. Conversely, for the AAR, the ratio 0.6 ± 0.05 (Porter, 1975) has been applied.

Results
The spatial distribution of glacial and periglacial landforms across the Bacivèr Cirque suggests the occurrence of a sequence of periods during the deglaciation of the cirque. The timing of these phases is constrained by 17 10 Be CRE ages inferred from glacial and periglacial records (Table 2).

Geomorphological evidence and CRE sampling strategy
Bacivèr is a large amphitheatre-shaped cirque (5 km long, 3.8 km wide) that forms the headwaters of the Malo River, which drains towards the Beret Plateau and flows into the Garonne River 6 km below. This complex cirque can be divided in three large geomorphological units ( Fig. 2): (i) the peaks and walls that form the head of the amphitheatre, (ii) the set of glacial, periglacial and paraglacial landforms distributed at the foot of the walls, and (iii) the large flat floor, composed of polished bedrock surfaces with small depressions and scattered erratic boulders. The cirque was heavily glaciated during the last glacial cycle, and the deglaciation following the LGM left widespread glacial and periglacial records (Fig. 3): (i) The cirque floor shows traces of intense glacial abrasion, with very scarce glacial deposits, mostly in the form of small, vegetated moraine ridges and some sparse erratic boulders. The short vertical distance between the cirque's walls and its floor implied low debris supply, which may explain the relatively small size of the moraine systems distributed across the cirque. Glacial erosion on the cirque floor generated several staggered, over-deepened basins that were occupied by lakes following deglaciation. From the central-lower part of the cirque, we collected three samples for CRE dating, two from polished surfaces (A-14, A-15) and one from an erratic boulder (A-10). The cirque gradually narrows downstream, becoming a steep U-shaped glacial valley until the confluence with the Beret Plateau at 1860 m. This lowest section mainly includes polished surfaces, with few remnants of moraine deposits due to the topographical setting, which favors intense erosion processes that have destroyed glacial accumulations. Here, we collected two samples from a polished bedrock surface (A-16, A-17) that is indicative of the onset of deglaciation of the bottom of the Bacivèr valley.
(ii) The Bacivèr Cirque constitutes a compound cirque (Barr and Spagnolo, 2015) composed of five smaller cirques with NW, N, E, SE and SW aspects (Fig. 3). Depending on the deepening of the individual cirques, distinct geomorphological sequences are present: (a) in steep and heavily carved out landforms (e.g. SE-and SW-exposed cirques), there are several generations of rock glaciers distributed in the upper part of the cirque basins that still show well-preserved furrows and ridges, with vegetation colonizing the lowest crests that reflects their current inactivity; in addition, some moraine ridges surround these landforms and enclose these small cirques distributed above the cirque floor; in the E cirque, there is also a small, collapsed, debris-covered glacier with longitudinal ridges and furrows; (b) elongated moraine ridges Highest moraines A-1 12.7 ± 0.8 (0.5) 12.9 ± 0.8 (0.5) 13.7 ± 0.9 (0.5) 13.9 ± 0.9 (0.5) A-2 12. fill the depressions in small glacial hollows at the contact between the cirque floor and the steep but short rock slopes (e.g. NW and E units).
The chronological study focused particularly on the eastern side of the main Bacivèr Cirque, at the foot of the summit ridge stretching N\ \S between Tuc de Marimanya and Cap del Muntanyó d'Àrreu (2602 m), where the main glacial and periglacial deposits are distributed. The northern side forms a 1 km-long, straight profile with peaks ranging from 2630 to 2530 m, and an elongate moraine, parallel to the wall at a distance of only 200 m, stretching over 800 m at altitudes between 2410 and 2460 m; here, we collected samples for CRE dating (A-11, A-12) from two moraine boulders. The southern hollow is more excavated by glacial and periglacial processes, including four small tongue-shaped moraines filling the hollows with polished surfaces among them. We collected two samples from the northernmost moraine ridge only 180 m from the wall (A-1, A-2), two samples from the central part of the polished bedrock 180 m from the wall (A-3, A-4), two samples from the external moraine of a debris-covered glacier 160 m from the wall (A-5, A-7), and one sample from a boulder located on a ridge inside this moraine 120 m from the wall (A-6). Moreover, we collected three samples from polished bedrock surfaces at distances of 100-200 m from the moraines and~300-400 from the wall, in both areas (A-8, A-9, A-13).
(iii) The peaks surrounding the cirque have relatively homogeneous altitudes ranging from 2500 to 2650 m, with small elevation differences of 150-200 m with respect to the base of the cirque floor. The cirque walls are steep and covered by a thick debris mantle generated by the intense frost shattering of the upper rock outcrops.

Geochronological data
The 17 samples collected for CRE dating yielded ages spanning from the B-A to the Mid-Late Holocene (Fig. 4). Glacial retreat during the last deglaciation in the Bacivèr valley started following the LGM according to the lowermost samples, which yielded ages of 18.6 ± 1.2 (A-17) and 14.2 ± 0.4 ka (A-16), respectively. Sample A-17 is considered an outlier because it is 3-4 ka older and its analytical uncertainty is significantly higher compared to the remaining samples of the dataset, and therefore is not further discussed.
Post-LGM glacial shrinking exposed the central-lowest part of the cirque floor, as revealed by the samples collected from the polished bedrock that yielded ages of 13.5 ± 0.4 (A-14) and 14.1 ± 0.4 ka (A-15) (13.8 ± 0.4 ka; n = 2) as well as from a scattered erratic boulder 15.4 ± 0.6 ka (A-10). The samples collected from polished bedrock in the upper part of the cirques show ages very similar to those of the central-lower sector. The sample from the northern side yielded an age of 15.2 ± 0.6 ka (A-13), whereas the samples from the southern sector of the Bacivèr Cirque were aged 15-14 ka (A-8 and A-9, respectively).
The two sampled boulders from the moraine of the northern side of the cirque returned ages of 10.0 ± 0.4 (A-12) and 13.3 ± 0.5 ka (A-11), although the former is not consistent with the chronostratigraphic sequence and thus was rejected. By contrast, the more robust chronology of the southern cirque unit gave slightly older ages of 12.7 ± 0.5 (A-1) and 12.8 ± 0.4 ka (A-2) (12.8 ± 0.5 ka; n = 2) for the boulders of the northernmost moraine. Samples from the polished bedrock dividing the moraine ridges within the cirque yielded an age range of~16-14 ka (A-3 and A-4, respectively). These are slightly older than the ages of boulders on the external moraine of a debris-covered glacier that were dated at 15-14 ka (A-5 and A-7, respectively). Finally, one sample from the external ridge inside this landform was dated at 7.2 ± 0.3 ka (A-6).

Discussion
The combination of geomorphological observations with 10 Be-dated glacial features distributed across the Bacivèr Cirque enabled us to infer space-time patterns of deglaciation in the cirque since the Late Glacial. In addition, this chronology also indicates the timing of formation of a debris-covered glacier located on the highest slopes of the cirque:

Chronology of the deglaciation and geomorphological significance
The age of polished glacial surfaces indicated that the deglaciation of the mouth of the cirque -where it gradually turns into a narrow Ushaped glacial valley -had occurred by 14.2 ± 0.8 ka, and consequently that the Bacivèr tributary glacier was disconnected from the main Garonne paleoglacier flowing from the upper Beret Plateau downslope through the Aran valley prior to this time (A-16; Fig. 5).
Glacial retreat was a very rapid process during the B-A as the entire cirque was deglaciated by~15-14 ka, according to polished surfaces from the central-lowest part of the cirque floor. These surfaces yielded a deglaciation age (13.8 ± 0.4 ka) that overlaps with that inferred from samples collected from the highest areas of the cirque, both from the eastern (~15-14 ka) and northern sides (15.0 ± 0.7 ka). There are minor inconsistencies between the CRE dataset (see Table 2) and the geomorphological stratigraphy, as some boulders are slightly older than the ages of the polished bedrock surfaces on which they rest. This may be due to nuclide inheritance: (i) the boulder may have fallen on the glacier during a period of glacial shrinking that exposed the highest section of the cirque walls, and it was subsequently transported supraglacially a few hundred meters prior to deposition on the cirque floor (Çiner et al., 2017); or (ii) the boulder was not reworked enough (inefficient erosion) given its closeness to the cirque headwall and the subsequent limited bearing and transport distance (see García-Ruiz, 1979). Also, variations in the local snow cover, which is much thicker on bedrock surfaces than on boulders that protrude considerably from the ground, might be responsible for the age differences. We therefore consider that within the geomorphological and analytical uncertainties of the CRE dating, the exposure of erratic boulders and related polished bedrock started roughly at the same time. Ice-moulded surfaces within the individual glacial cirques provided similar deglaciation ages ( Fig. 3; 16-14 ka). This phase, however, is compatible with the occurrence of periods of relative glacial stability, which may have generated the small moraine ridges distributed on the main cirque floor (Fig. 3), as well as the small debris-covered glacier located in the eastern hollow of the Bacivèr Cirque. Here, two boulders returned slightly younger ages (~15-14 ka; see Table 2) than the neighboring polished bedrock (~16-14 ka). These exposure ages from the highest surfaces suggest that the rapid glacial shrinking recorded during the B-A probably favored the disappearance of glaciers in the cirque bỹ 15-14 ka, with some stagnant ice masses in favorable topographical settings under the protection of the wall and with an intense debris supply. In short, the considerable homogeneity of the CRE results in the area seems to indicate a rapid retreat of 3 km of the glacier at the beginning of B-A, with only small glaciers (200-300 m long) sheltered by the wall surviving, and under an intense paraglacial activity. This abundant debris supply caused the residual glaciers to evolve towards a range of typologies, including debris-covered glaciers in some cirques, and rock glaciers in others. These landforms are widespread features in deglacierizing mountains that transition from glaciated into paraglacial landscapes . The coincidence of the ages of the boulders belonging to these formations and those of the rock platforms where they rested indicates the short time elapsed between the general deglaciation and the collapse of these residual small glaciers, which determined the stabilization of their deposits (Fig. 3).
Exposure ages of moraine boulders revealed that the highest moraine was abandoned and finally stabilized at 12.8 ± 0.5 ka, whereas one sample from the northernmost dated landform (Fig. 3)    slightly older age of 13.3 ± 0.5 ka; the average of the three moraine samples results in an age of 12.9 ± 0.3 ka. Remarkably, we highlight that these are minimum ages that do not include erosion and snow corrections (Table 3). The uncertainty range of CRE dating and the small number samples makes it difficult to differentiate whether the moraine stabilized at the end of the B-A or at the beginning of the YD. The absence of further glacial landforms at higher elevations suggests that the Bacivèr Cirque has not accumulated glacial ice since~1 2.9 ka. However, the degradation of these glaciers enhanced paraglacial dynamics and the persistence of small ice masses as small debris-covered glaciers. A sample collected from a ridge of one such landform yielded an exposure age of 7.2 ± 0.3 ka. We interpret this as the age of stabilization of the boulder once the inner ice melted away and mobility ceased, which might be representative of the collapse/ final melting of the inner sector of the debris-covered glacier (Fig. 6). The location of the sample on a ridge 100 m distant from the wall

Late Quaternary glacial dynamics in the Central Pyrenees and Iberian Peninsula
The Iberian ranges are among the mountains where knowledge about Pleistocene glacial evolution has most improved in the last decade (Oliva et al., 2019). New chronological information about their glacial history from recent years has been complemented with a better knowledge of the spatial and temporal dynamics of periglacial processes that Fig. 7. Normalized probability distribution functions (PDF) of exposure ages vs. temperature evolution since the LGM based on the δ 18 O record from the NGRIP ice core from Greenland (time periods are defined after Rasmussen et al. (2014)). The plots of the units result from the sum of the individual PDF of the samples belonging to them. reshaped the landscape fashioned by Quaternary glaciers, particularly during the Holocene .
In this Iberian context, together with the Sierra Nevada massif and the Cantabrian Mountains, the Pyrenees stand out as the mountain ranges where the numerous glacial studies have produced the best chronologies of glacial oscillations (Table 4). There is still an open debate and divergent information with regards to the timing of the maximum glacial advance during the last glacial cycle on both the northern and southern slope of the Pyrenees (Delmas et al., 2021). There is more consensus, and much more homogeneous information, about the timing of glacial advances and retreats during colder and warmer periods during T-1. However, most research has been conducted in the major glacial valleys and highest mountain cirques, whereas lower-altitude catchments remain poorly investigated. To shed light on the role of altitude and topographical conditions controlling glacial oscillations during T-1, we reconstructed the deglaciation process in the upper Garonne Basin focusing on one valley of its headwaters, the Bacivèr Cirque. The fact that its cirque summits reach only 2600-2650 m determined a pattern of deglaciation similar in both extent and timing to some of the tributaries of the Gállego Valley (Palacios et al., 2017b) or the Ariège Valley  -where summits are around 2600-2700 mrather than that reported for glaciers developed at the foot of summits around 3000 m (García-Ruiz et al., 2016;Palacios et al., 2017b;Pallàs et al., 2006).
In the Pyrenees, the highest peaks exceeding 2900 m remained covered by extensive, thick and strongly erosive ice caps during the Late Quaternary, while in lower ranges between 2400 and 2800 m glaciers were much thinner. Consequently, many of these ridges emerged above the glaciers as nunataks (Delmas et al., 2021). This is of great importance for understanding the local geomorphological evolution during T-1, especially on substrates of crystalline rocks. The cirques of the highest peaks were intensely eroded by glacial processes, which removed the entire weathering mantle. As such, when glaciers retreated, paraglacial dynamics were less intense in these areas. By contrast, in lower-altitude massifs, the weathered mantle has been partially preserved on the walls (less affected by glacial erosion) and when glaciers retreated, paraglacial processes were comparatively more intense. For this reason, a greater diversity of glacial and periglacial sedimentary landforms, which can provide a detailed picture of the environmental evolution during deglaciation, is found in the interior of low-altitude crystalline cirques. This geomorphological pattern has also been identified in other mid-altitude cirques of the Central Pyrenees, such as Catieras, Piniecho and Brazato (Palacios et al., 2017b), as well as in other relatively high cirques in the Eastern Pyrenees, including Malniu and Perafita .
The cold conditions that prevailed during the OD favored the reoccupation of the valley floors by glaciers at~17-16 ka (Fig. 7), as occurred in several valleys in the Central and Eastern Pyrenees (Palacios et al., 2017a) (Table 4). At this time, the upper Garonne Basin must still have been extensively glaciated, with moraine systems located in the main valley also fed by the glacier descending from the Bacivèr Cirque that contributed ice to the Garonne paleoglacier (Fernandes et al., 2017). However, this phase was reconstructed based on the identification of moraine systems alone, and no direct ages are yet available for the upper Garonne valley.
As temperatures rose by 3-5°C in western Europe during the B-A Interstadial (Clark et al., 2012), the Garonne paleoglacier receded rapidly and several small alpine or cirque glaciers remained, individualized within the headwaters of the highest valleys. CRE dates of the lowest sections of the Bacivèr valley suggest that the disconnection between the Beret Plateau (i.e. Garonne paleoglacier) and the Bacivèr paleoglaciers took place prior to 14.2 ± 0.8 ka. At that time, the glacier had a length of 3.8 km and had shrunk significantly since the maximum ice expansion of the last glacial cycle, when the Garonne paleoglacier reached a length of 89 km (Fernandes et al., 2017). In fact, the overlapping between the samples collected from the cirque floor and the polished bedrock surfaces of the higher glacial thresholds suggests that deglaciation was a very rapid process and that the entire cirque (~10 km 2 ) was ice-free by~15-14 ka. The presence of small (undated) moraine ridges distributed across the Bacivèr Cirque (Figs. 2c and 3) reveals the occurrence of short periods of glacial advance or standstills during the B-A. However, as small moraines, these ridges must have formed during short periods favoring shifts in glacier mass balance that resulted in stabilization or even very limited glacier growth within an overall trend of warming and accelerated retreat (see e.g. Chandler et al., 2016). The partial overlapping of the CRE uncertainty ranges impedes precise constraints of particular events within the cirque's deglaciation chronology, which has also occurred in other massifs in the Pyrenees where the transition between different T-1 cold/warm periods cannot be detected with high precision (e.g. Pallàs et al., 2006). In any case, glaciers in the Pyrenees during the B-A may only have persisted in the highest northern cirques of the highest massifs whereas low-altitude catchments must have been ice-free during this period (Oliva et al., 2019), as most of the glaciated massifs were during the OD in the Mediterranean region (Palacios et al., 2017b(Palacios et al., , 2017a. In the Alps, glacial extents were significantly reduced during the B-A and glaciers only persisted in the highest sectors of the highest valleys (Ivy-Ochs, 2015).
The cooling from the B-A to the YD in Western Europe has been quantified at 5-10°C (Clark et al., 2012), although it might have been lower in the Iberian Peninsula: terrestrial records suggest a temperature decrease of 2.5°C (Iriarte-Chiapusso et al., 2016) whereas marine records from the Alborán Sea (Cacho et al., 2002) and the Portuguese coast (Rodrigues et al., 2010) point to temperatures lower by 4-5°C. In parallel to the lower temperatures that generally prevailed during the YD, significant hydrological shifts were recorded (Bartolomé et al., 2015;Cheng et al., 2020;Rea et al., 2020) that must have affected the ELA in the Central Pyrenees. In the case of the Bacivèr cirque, we assume that these small glaciers readvanced or stabilized at the transition between the B-A and the YD, probably due to the protection of a thick debris cover on a glacier that did not reach the cirque floor. These cirque glaciers were small (Fig. 5), with lengths of 0.2-0.3 km and surface extents of 13-18 ha, and with their fronts at elevations of~2400 m. These glaciers represented only~0.3% of the length of the Garonne paleoglacier during the maximum ice extent of the last glacial cycle (Fernandes et al., 2017). This pattern is similar to that in in the Upper Gállego Valley (Catieras and Piniecho cirques) where the B-A frontal moraines are located 0.4-0.6 km from the headwall, at 2300-2350 m, in cirques where the maximum elevations are~2700 m, slightly higher than in the Bacivèr Cirque (Palacios et al., 2017a). Similar to the Bacivèr cirque, in the Alps, glacial expansion during the B-A/YD transition left two moraine systems that were dated betweeñ 13.5 and~12 ka (Ivy-Ochs, 2015); in the highest mountains in Greece on Mount Olympus moraines stabilized at~13.5-11.7 ka (Styllas et al., 2018) and on Mount Chelmos at~13.1-10.5 ka (Pope et al., 2017); and in Anatolian mountains two phases occurred at~13 and~11.5 ka, both in the Kaçkar Mountain range (Akçar et al., 2007) and Uludağ Mountain (Zahno et al., 2010). Based on the reconstruction of the two dated moraine ridges, during the transition between the B-A and the YD, the ELA in the Bacivèr Cirque must have been at 2485 and 2504 m, respectively (average 2495 m; Table 5). The current regional 0°C isotherm in the Central Pyrenees lies at~2950 m , which roughly coincides with the ELAs in some still glaciated massifs such as Monte Perdido and Maladeta (Chueca et al., 2005). Therefore, assuming an ELA depression of 455 m with respect to the current regional estimate of 2950 m, an average lapse rate of 0.65°C 100 m −1 and no change in precipitationas inferred for some areas in the Mediterranean region, such as the Maritime Alps (Spagnolo and Ribolini, 2019)summer temperatures must have beeñ 3.0°C lower than at present during the transition between the B-A and the YD to allow the formation or stabilization of such small glaciers. Mass balance models suggest slightly lower temperatures in the Ariège Valley, ranging from 3.9 and 5.1°C, based on the reconstruction of glaciers from the moraine systems of two cirques .
Overall, our results and those from other studies in the Pyrenees show that the time spanning from the early B-A to the YD was a major driver of landscape change in the high sectors of the Pyrenees as: (i) prevailing warm conditions promoted the definitive disappearance of glaciers in most cirques, particularly in low-to mid-altitude cirques, (ii) glacial shrinking favored the formation of debris-covered glaciers that extend over the cirque floors, and (iii) glacial retreat was followed by very intense paraglacial dynamics that favored the formation of permafrost-related landforms such as rock glaciers and protalus lobes  as well as abundant slope failures in formerly glaciated areas (Fernandes et al., 2020).
The temperature increase of~4°C in western Europe recorded at the onset of the Holocene (Clark et al., 2012) favored the disappearance of YD glaciers. However, as detected in some cirques at the foot of the highest peaks, at~2900-3000 m in the Central and Eastern Pyrenees (Table 4), intense paraglacial adjustment following glacial shrinking led to the formation of rock glaciers Oliva et al., 2016). Indeed, the deep glacial hollows excavated in steep slopes surrounding the Bacivèr Cirque included well-developed rock glaciers. Interestingly, in an east-facing hollow that was less excavated by glacial and periglacial erosion, a debris-covered glacier formed once climate conditions became unfavorable for glacial activity. According to the exposure age of the sampled boulder (sample A-6) suggesting its geomorphic stabilization, this landform must have contained glacial ice until at least~7.2 ka, when it finally melted away during the Holocene Thermal Maximum (Renssen et al., 2009). This timing is very similar to the stabilization of rock glaciers in the Central and Eastern Pyrenees García-Ruiz et al., 2016;Palacios et al., 2015a), as well as to the final collapse of the debris-covered glaciers that existed in the Sierra de la Demanda, Iberian Range, that persisted until 7-6 ka . As revealed by present-day analogs, these landforms can undergo significant morphodynamic changes in deglacierizing mountains and may have a long residence time in the landscape . Ice can thus persist protected beneath the debris cover in sheltered areas, little affected by atmospheric temperatures, as also detected in Iceland Fernández-Fernández et al., 2020;Tanarro et al., 2019).
However, recent research has also shown that debris-covered glaciers can be more sensitive to climate variability than initially thought, as some features may respond rapidly to changes in temperature and precipitation (Charton et al., 2020). Therefore, when interpreting past glacial oscillations both climate changes as well as paraglacial processes must be taken into account. This means that the great intensity of the paraglacial readjustment in these cirques, under relatively lowaltitude peaks in crystalline rocks and where the rock walls retain most of the weathered mantle, interferes with the effects of climate on glacier dynamics. In any case, it should be highlighted that the same chronological pattern for the deglacial period was also observed in cirques with similar characteristics in the Central (Palacios et al., 2017b) and Eastern Pyrenees . The new debris cover on the glacier surface, supplied by paraglacial processes, can promote glacier advance or long-term stability by drastically reducing the ablation and shifting mass balance (Herreid and Pellicciotti, 2020). As such, paraglacial processes can interrupt the deglaciation process to some extent even without climatic influence (Hambrey et al., 2019(Hambrey et al., , 2008Jones et al., 2019;Rowan et al., 2015). The factors controlling the advances/retreats of debris-covered glaciers, as well as their collapse, are not fully understood, and further research is needed to clarify their response to short-and long-term climate trends. This is particularly important in the current global warming scenario, where glacier shrinking is accelerating and paraglacial processes are delivering large amounts of sediment to the surface of glaciers, considering that 7.3% of the area of mountain glaciers is debris-covered and this percentage is expected to increase in response to increasing temperatures (Herreid and Pellicciotti, 2020).
Since the inferred disappearance of glacial ice of the cirque at~7.2 ka, the natural evolution of the landscape of the Bacivèr Cirque has been driven mainly by nival processes and periglacial dynamics (Fernandes et al., 2017). Current morphodynamics are associated with the occurrence of a seasonal frost ground thermal regime, as permafrost conditions are only found extensively in areas above 2900 m in the Central Pyrenees , much higher than the highest summits of the Bacivèr Cirque.
In summary, the reconstructed temporal pattern of deglaciation of the Bacivèr Cirque during T-1 is fully consistent with the timing reported in other sectors of the Pyrenees Palacios et al., 2017b) and other Iberian mountains (see summary in Oliva et al., 2019) with accelerated glacier retreat during the B-A, and the subsequent activation of paraglacial processes and a minor glacial re-expansion during the B-A/YD transition. No disparity in timing of maximum ice extent across the Mediterranean mountains is observed during the last deglaciation. The ages of glacial advances and retreats in Iberia during T-1 are similar to those that occurred in centralnorthern Europe, such as Iceland , the British Isles (Barth et al., 2016), and the Alps (Ivy-Ochs, 2015; Moran et al., 2016). The chronology is also similar to that reported in other temperate European mountains, such as the Anatolian Peninsula (Köse et al., 2019;Sarıkaya et al., 2017), the Carpathians (Gheorghiu et al., 2015;Makos et al., 2018), the southern Balkans (Styllas et al., 2018), the Dinaric Alps (Žebre et al., 2019), and Atlas Mountains . Indeed, the timing of oscillations resembles the past evolution of the Scandinavian Ice Sheet, which retreated considerably during the B-A and regrew and advanced tens of km as temperatures declined during the transition towards the YD (Greenwood et al., 2015;Hughes et al., 2016;Mangerud et al., 2016). All in all, the data seems to indicate that deglaciation followed very homogeneous climatic patterns throughout Europe, with minor differences imposed by local topoclimatic conditions.

Conclusions
One of the main attractions of the natural heritage of the Central Pyrenees is its glacial landscape. During Quaternary colder phases, glaciers have fashioned the valleys and cirques, which were subsequently reshaped by periglacial conditions during the warmer interglacial periods. Despite the fact that the highest massifs of the Central Pyrenees are some of the mid-latitude mountain areas where glacial evolution has been studied in the greatest detail, the timing of development of glacial phenomena remains poorly investigated for cirques at lower altitudes. This is the case of the Bacivèr Cirque, from where we introduced CRE ages of glacial landforms that are indicative of glacial oscillations during T-1 in the upper Garonne valley.
The Atlantic-influenced upper Garonne Basin favored the development of the longest glacier in the entire range (~90 km) despite having its headwaters at elevations of 2600-2650 m, 600-800 m lower than the highest peaks in the Pyrenees. The chronology of deglaciation of one of the cirques in the headwaters of the Garonne paleoglacier also contributes to a better understanding of glacial oscillations on the northern slope of the Pyrenees. The Bacivèr Cirque, with its wide cirque floor located at 2200-2400 m, was located above the regional ELA during the maximum ice extent of the last glacial cycle. As temperatures increased following the LGM, glaciers rapidly retreated to the headwaters of the highest valleys and the glacier flowing downvalley from the Bacivèr Cirque became disconnected from the Garonne Glacier prior to~15-14 ka, when the lowest part of the cirque became ice free. Glacial recession was enhanced during the B-A at~15-14 ka, with the cirque likely being fully deglaciated. Small cirque glaciers formed or remained at the foot of the wall during the transition between the B-A and the YD. These small glaciers were affected by paraglacial readjustment of the slopes, which triggered their transformation into rock glaciers and a debris-covered glacier. Their fronts collapsed almost immediately, but in some cases, their upper sections remained active during the Early Holocene, until at least~7 ka; since then, periglacial slope processes and nival activity have shaped the highest parts of the massif. The chronology of glacial advances and retreats during T-1 reconstructed from the Bacivèr Cirque is similar to that reported from other lower Pyrenean glaciers. Therefore, the period spanning from the early B-A to the transition towards the YD dramatically transformed the mountain landscape of the Pyrenees and favored the development of the great geomorphological diversity of glacial and periglacial landforms that exists today in many cirques of this mountain range.
The glacial and periglacial landscape of the Bacivèr Cirque inherited from past periods and slightly reshaped during the current interglacial is also being intensely transformed by human activities. The expansion of the neighboring ski resort, with the ski slopes and associated infrastructure, has already destroyed geomorphic evidence. The findings presented in this study are thus clear evidence of the richness of crucial information that this cirque contains, and of the need to preserve its landscape for future generations.

Declaration of competing interest
The authors of this manuscript declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.