Marcelo Fernandes
email: marcelo.fernandes@live.com

M Oliva

G Vieira

D Palacios

J M Fernández-Fernández

Irene Schimmelpfennig

D Antoniades

J Garcia-Oteyza

Schimmelpfennig

Aster Team

Glacial oscillations during the Bølling-Allerød Interstadial-Younger Dryas transition in the Ruda Valley, Central Pyrenees

Keywords: Central Pyrenees, Bølling-Allerød, Younger Dryas, Cosmic-Ray Exposure dating, deglaciation 58 59

come

L'archive ouverte

Introduction 60

The present-day landscape of the highest mid-latitude mountain environments is dominated by the effects of 61 glacial dynamics that prevailed during the last glacial cycle and of post-glacial processes that reshaped glacial 62 landforms during the subsequent deglaciation. Following the end of the Last Glacial Maximum (LGM) at 19-63 20 ka (Clark et al., 2009), most glaciated mountain regions underwent a rapid deglaciation only interrupted by 64 brief periods of glacial growth. Glacial disappearance favoured paraglacial readjustment and intense 65 periglacial and slope processes, which reshaped the already deglaciated valleys and cirques [START_REF] Ballantyne | After the ice: Holocene geomorphic activity in the Scottish Highlands, Scottish Barr, I[END_REF]. 66

The current spatial distribution of glacial and periglacial features in the high elevations of most mid-latitude 67 ice-free mountain ranges -such as the Pyrenees, where this research focuses on -results from the climatic 68 shifts during Termination-1 (hereafter T-1), the period spanning from 20-19 to 11.7 ka (Denton et al., 2014).

69

At present, mid-latitude mountains such as those of the Mediterranean region are only sparsely glaciated, with 70 small glaciers located in sheltered cirques (Hughes, 2014). However, there are several other mountain ranges 71 located between 30 and 50º latitude North, such as the Alps or the Himalaya, which are still heavily glaciated 72 and record recent dramatic glacial recession trends (Huss et al., 2021). An accurate comprehension of the 73 climatic sensitivity of current glaciated mountain landscapes can help to understand the timing and magnitude 74 of glacial change during the rapid climate shifts of T-1. The alternation of millennial-scale cold/warm periods 75

during T-1 promoted a temporal pattern of glacial advances and retreats across the Northern Hemisphere 76 [START_REF] Allard | Heinrich Stadial aridity forced Mediterranean-wide glacier 588 retreat in the last cold stage[END_REF]Rea et al., 2020). Ice core records from the Greenland interior show the following 77 succession of climatic phases in the North Atlantic region during T-1: a cold climate regime during the Oldest 78 Dryas (OD; GS-2.1a: 17.5-14.6 ka), followed by the general warm Bølling-Allerød (B-A) Interstadial (GI-1: 79 14.6-12.9 ka), with a subsequent intensification of cold conditions during the Younger Dryas (YD; GS-1: 12.9-80 11.7 ka) (Rasmussen et al., 2014).

81

In the Iberian mountains, recent advances in constraining the chronology of glacial activity during T-1 using 82

Cosmic-Ray Exposure (CRE) dating have revealed contrasting spatio-temporal patterns (García- Ruiz et al., 83 2016b;Oliva et al., 2019;Palacios et al., 2017a). Whereas the maximum ice extent occurred asynchronously, 84 glacial oscillations during T-1 show a more synchronous response: colder phases favoured glacial growth, 85 whereas warm periods promoted glacial recession (Oliva et al., 2019) and the evolution from debris-free 86 glaciers to rock glaciers such as those in the Iberian Range (García-Ruiz et al., 2020), Cantabrian Mountains 87 (Rodríguez-Rodríguez et al., 2016) and Central (Palacios et al, 2017a, b) and Eastern Pyrenees (Andrés et al., 88 2018). To date, there is geomorphic evidence of glacial advances during the OD and YD, as well as of glacial 89 retreat during the B-A and following the YD in the Cantabrian Mountains (Rodríguez-Rodríguez et al., 2017), 90

Central Range (Carrasco et al., 2015;Palacios et al., 2012b), and Sierra Nevada (Gómez-Ortiz et al., 2012;91 Palacios et al., 2016). However, CRE ages require a careful interpretation and support of geomorphological 92 criteria, as their uncertainty ranges may be greater than the time spanning the T-1 phases (Oliva et al., 2021a).

93

The Pyrenees were the most extensively glaciated mountain system during the maximum ice extent of the last 94 glaciation in Iberia, hosting ~50% of the peninsula's glacial surface (Oliva et al., 2021b). The chronology of 95 glacier fluctuations during T-1 seems to have followed a similar temporal pattern, although available data is 96 principally derived from the range's southern slope [START_REF] Andrés | Timing of deglaciation and rock glacier origin in the southeastern Pyrenees: a 594 review and new data[END_REF]Palacios et al., 2015b[START_REF] Oliva | Following ice retreat, the slopes were affected by intense 420 readjustments (large, catastrophic rock slope failures, rock slope deformations, rockfalls, and rockslides), 421 particularly in areas where the local lithostructure was prone to instability[END_REF]Palacios et al., , 2017a)). 97

Pyrenean glaciers readvanced to the valley bottoms during the OD, with ice tongues > 15 km long. During the 98 B-A, by 15-14 ka, high temperatures favoured rapid deglaciation and glaciers disappeared from relatively low 99 cirques shaped under peaks of ~2600 m a.s.l. (Oliva et al., 2021a). A subsequent readvance occurred during 100 the YD, with smaller glaciers of only ~4 km length confined within the cirques of the highest massifs above 101 2200 m (Crest et al., 2017;García-Ruiz et al., 2016a;Jomelli et al., 2020;[START_REF] Oliva | Following ice retreat, the slopes were affected by intense 420 readjustments (large, catastrophic rock slope failures, rock slope deformations, rockfalls, and rockslides), 421 particularly in areas where the local lithostructure was prone to instability[END_REF]Palacios et al., , 2015a;;Pallàs 102 et al., 2006). Following both cold periods, glacial retreat triggered intense paraglacial dynamics at ~15-14 ka 103 and ~11-10.5 ka, respectively, which led to the formation of rock glaciers and debris-covered glaciers, 104

particularly in cirques below 2800 m [START_REF] Andrés | Timing of deglaciation and rock glacier origin in the southeastern Pyrenees: a 594 review and new data[END_REF][START_REF] Oliva | Following ice retreat, the slopes were affected by intense 420 readjustments (large, catastrophic rock slope failures, rock slope deformations, rockfalls, and rockslides), 421 particularly in areas where the local lithostructure was prone to instability[END_REF].

105

The timing of glacial oscillations is still poorly constrained on the northern slope of the Pyrenees, such as in 106

the Upper Garonne Basin where this work focuses on. Recently, Reixach et al. (2021) have summarized the 107 glacial evolution from 19 to 12 ka in the northern Pyrenees, reporting a spatial and temporal pattern for T-1 108

glacial oscillations similar to that on the southern slopes of the range. Jomelli et al. (2020) confirmed glacial 109 growth during the YD in the Ariège Valley, where glaciers persisted in the highest cirques until the early 110

Holocene. Another recent study that focused on the Bacivèr Cirque, near our study area, revealed a rapid 111 deglaciation during the B-A, followed by glacial readvance and the formation of the highest moraines at the 112 B-A/YD transition at ~12.9 ka; subsequently, intense paraglacial readjustment favoured the formation of 113 debris-covered glaciers and rock glaciers inside these moraine systems (Oliva et al., 2021a).

114

In order to complement our understanding of the glacial history of the Pyrenees, this study uses new data to 115 improve the deglaciation chronology of the Upper Garonne Basin, the largest glacial system in this mountain 116 range, during the local last glacial maximum (lLGM; Fernandes et al., 2017). Glacial retreat in the Ruda Valley,117 where this work focuses on, left widespread geomorphic evidence of glacial and paraglacial origin that can be 118 used to reconstruct the last stages of long-term deglaciation in the Central Pyrenees.

119

Our study, therefore, has the following objectives: The Pyrenees stretch E-W over 400 km across north Iberia, with the highest peaks (> 3000 m) in the central 130 part of the range. This study focuses on the Ruda Valley, located between latitudes 42º40' N -42º42' N and 131 longitudes 0º57' E -1º00' E (Fig. 1), that constitutes the headwaters of the Upper Garonne Basin. The area is 132 included in the peripheral protection area of the National Park of Aigüestortes and Sant Maurici Lake.

133 Figure 1 134

The Upper Garonne Valley shows the typical U-shaped cross-section of glacial valleys, and the Garonne River 135 flows to the N-NW through the central axis of the Pyrenees. Several tributary rivers drain the glacial valleys 136 descending from glacial cirques carved below peaks from 2800 to 3000 m. Glacial cirques are prevailingly 137 NE-oriented, with floors generally between 2200 and 2400 m (Lopes et al., 2018).

138

The U-shaped Ruda Valley, drained by the Ruda River that runs 12 km northwards, constitutes the headwaters 139 of the Upper Garonne Basin. Its catchment encompasses 34 km 2 , with an elevation range of ~1300 m between 140 Saboredo Peak (2827 m) and the village of Baqueira (~1500 m). The Saboredo Cirque is located in the 141 headwaters of the Ruda Valley; the cirque is N-exposed, extending over 5.7 km 2 , with its floor at ~2250-2400 142 m. Saboredo (2827 m), Sendrosa (2702 m) and Locampo (2658 m) peaks are the highest elevations of the 143 cirque, while its toe is located at ~2250 m.

144

The present-day mean annual air temperature (MAAT) at Port de la Bonaigua station (2266 m) is 3 ºC, with 145 the regional 0 ºC isotherm in the Central Pyrenees estimated to be at 2950 m (López- Moreno et al., 2019). 146

Annual precipitation reaches 1227 mm and falls mostly as snow (average depth of 1 m) from October to May, 147 with a snow duration of 7-8 months in the cirque floor. The combination of low temperatures and abundant 148 snowfall determines the scarce vegetation cover above the treeline (2200-2300 m), where barren rocky terrain 149 is widespread, together with small alpine meadows.

150

The landscape of the Ruda Valley resulted from the interaction of tectonics and vertical incision, both during 151

Quaternary glaciations and interglacial periods, the latter dominated by alluvial and mass-wasting processes 152 (Fernandes et al., 2017). We produced a geomorphological map focusing on glacial and periglacial landforms based on recent general 170 geomorphological maps from the area (Fernandes et al., 2017;Fernandes et al., 2021). Here we used the same 171 approach based on the combination of aerial orthophotographs (0.25 m cell size) and high-resolution LIDAR 172 digital elevation models (point density of 0.5-2 m -2 ) from the "Institut Cartogràfic i Geològic de Catalunya" 173 (http://www.icc.cat/appdownloads). The symbology of geomorphological landforms followed Joly (1997).

174

Field strategy and sampling 175

We collected 12 granite and granodiorite samples for CRE dating, using a hammer and chisel, from elevations 176 spanning from 1862 to 2470 m (Table 1). To ensure that samples were intact since glacial retreat, we targeted 177 at well-anchored moraine boulders and glacially polished bedrock that stood out within their surroundings, 178 which minimizes the risk of having been buried under sediments (Fig. 2). To ensure the best cosmic-ray flux 179 reception of the sampling sites, we restricted sample extraction to flat-topped and gentle surfaces, and sampled 180 a thickness of up to 3-5 cm. Corrections for partial shielding due to the surrounding topography were 181 implemented using the ArcGIS toolbox by Li (2018), which applies equations describing the effects of 182 geometric shielding on the sampling sites and a model to predict the nuclide production rate at depth in sloped 183 surfaces (Dunne et al. 1999); it required only a point shapefile of the sampling sites including strike/dip of the 184 sampling surfaces and a digital elevation model of 5 m spatial resolution.

185 subsequently attacked with 4-6 successive partial dissolutions with HF to remove the atmospheric 10 Be and 199 discard the remaining impurities. In order to be sure that the quartz was pure, samples were examined under a 200 binocular microscope. As a result, ~20 g of pure quartz per sample was used to extract the 10 Be (Table 2). After 201 adding ~150 µL of an in-house manufactured (from a phenakite crystal) 9 Be carrier solution (spike, 202 concentration: 3025 ± 9 µg g 1 ; Merchel et al., 2008), the purified quartz was dissolved with 48% concentrated 203 HF (3.5 mL per g of quartz + 30 mL in excess). Following total dissolution, the resulting solution was 204 evaporated until dryness, and the residue was recovered with HCl (7.1 molar). Be was then precipitated at 205 pH=8 to beryllium hydroxide (Be(OH) 2 ) by means of ammonia (NH 3 ), and separated from other elements in 206 resin columns: a Dowex 1x8 anionic exchange column to remove elements such as Fe, Mn and Ti, and a 207 Dowex 50W x8 cationic exchange column to discard B and separate Be from Al (Merchel and Herpers, 1999).

208

The final eluted Be was again precipitated, dried and oxidized to BeO at 700 ºC. Finally, the targets for the 209 accelerator mass spectrometer (AMS) measurements were prepared by mixing with niobium powder at an 210 approximate 1:1 proportion and pressing the mixture into copper cathodes.

211

The targets were analysed at the 'Accelerator pour les Sciences de la Terre, Environnement et Risques'

212

(ASTER) national AMS facility at CEREGE in order to measure the 10 Be/ 9 Be ratio from which the 10 Be 213 concentration was later inferred (Table 2). AMS measurements were calibrated against the in-house standard 214 STD-11 with an assigned 10 Be/ 9 Be ratio of (1.191 ± 0.013) × 10 11 (Braucher et al., 2015). Based on calibration 215 with the previous ASTER AMS standard NIST-27900 (SRM4325), the new STD-11 is fully compatible with 216 07KNSTD. The analytical 1σ uncertainties include uncertainties in the AMS counting statistics, an external 217 0.5% AMS error [START_REF] Arnold | The French accelerator mass spectrometry facility ASTER: Improved performance 597 and developments[END_REF] and the uncertainty related to the chemical blank correction. The 10 Be 218 half-life considered was (1.3870 ± 0.0012) × 10 6 years (Chmeleff et al., 2010;Korschinek et al., 2010).

219

Table 2 220

We calculated 10 Be exposure ages by using the CREp online calculator (Martin et al., 2017; available online 221 at: http://crep.crpg.cnrs-nancy.fr/#/), where we selected the following settings: LSD (Lifton-Sato-Dunai) 222 elevation/latitude scaling scheme (Lifton et al., 2014), ERA40 atmospheric model (Uppala et al., 2005) and 223 geomagnetic database based on the LSD framework (Lifton et al., 2014). These settings yielded a world-wide 224 mean sea level high latitude (SLHL) 10 Be production rate of 3.98 ± 0.22 atoms g 1 yr 1 . Exposure ages of the 225 samples with their 1σ full and analytical uncertainties are shown in Table 2. The uncertainties discussed 226

throughout the text include analytical and production rate error unless otherwise stated.

227

In order to evaluate the potential impact of erosion on the exposure ages, we assumed a steady erosion rate (1 228 mm ka -1 ) following a conservative maximum value in such a lithological setting according to [START_REF] André | Rates of Postglacial rock weathering on glacially scoured outcrops (Abisko-Riksgränsen 590 area, 68°N)[END_REF].

229

The impact of snow cover was explored assuming that snow currently remains in the area an annual average 230 of 7.5 months at 2200-2300 m with a mean thickness of 100 cm -and assuming a snow density of 0.2 g cm -3 231 (Styllas et al., 2018) -(Table 3), based on the available data from the Bonaigua weather station (Fig. of cold/warm periods and their associated moisture supplies. Finally, the chi-squared test according to Ward 238 and Wilson (1978) was applied to the samples belonging to each geomorphological unit established in the 239 geomorphological map (Table 1) in order to detect potential exposure age outliers (i.e. "too old" or "too 240 young"; Heyman et al., 2011). No outliers were detected (Table 3).

241

Table 3 242 3.4 Glacier reconstruction and Equilibrium-Line Altitude (ELA) calculation

243

A three-dimensional paleoglacier reconstruction was performed for the different glacial phases using the 244 'GLaRe' ArcGIS toolbox devised by Pellitero et al. (2016). Former ice thickness was estimated by applying a 245 perfect-plasticity physical-based rheological model along flowlines from the termini to the headwall (Van der 246 Veen, 1999;Benn and Hulton, 2010). Such a toolbox requires only a flowline, a tentative paleoglacier 247 geometry -approximated from lateral or frontal moraines -and a digital elevation model. Ice thickness was 248 modelled by using a constant shear stress of 100 kPa (Paterson, 1994;Benn and Hulton;2010). The effect of 249 lateral side drag was corrected by integrating shape factors (F-factor) obtained from cross-sections (Schilling 250 and Hollin, 1981).

251

ELAs were calculated by using the automatic toolbox 'ELA calculation' developed by Pellitero et al. (2015).

252

The selected methods were the Accumulation Area Ratio (Porter, 1975; AARs: 0.6 ± 0.05) and the Area 253

Altitude Balance Ratio (AABR; Osmaston, 2005). We considered two balance ratios for the AABR method: 254

1.75 ± 0.71 (global) and 1.9 ± 0.81 (mid-latitude marine glaciers) (Rea, 2009). The average ELAs were used 255 to discuss the paleoclimatic implications of the identified glacier advance/stillstand phases.

256

Results

257

The The 12 CRE samples collected from the Ruda Valley returned a chronology spanning from ~15 ± 0.9 to 11.8 299 ± 0.7 ka, with all ages being consistent with the geomorphological distribution of the dated surfaces (Fig. 4).

300

Polished surfaces from the highest parts of the Ruda Valley located ~200-250 m below the cirques provided 301 CRE ages of 15.0 ± 0.9 (ARAN-40) and 13.8 ± 0.8 ka (ARAN-39, Fig. 4A) (average: 14.4 ± 1.2 ka). A 302 similar age was obtained from a sample collected from a truncated spur dividing the Sendrosa and Saboredo 303 cirques, which yielded 14.8 ± 0.9 ka (ARAN-31). A sample collected from a polished surface at the edge of 304 the Saboredo cirque floor (Fig. 3) reported an age of 14.0 ± 0.8 ka (ARAN-34, Fig. 4D).

305

Moraine boulders from the external moraine located in the lower overdeepened basin yielded ages of 13.3 ± 306 0.8 (ARAN-38) and 13.7 ± 1.0 ka (ARAN-37) (average: 13.5 ± 0.9 ka). A boulder from the internal moraine 307 ridge located 110 m above the former moraine provided an age of 13.1 ± 0.8 ka (ARAN-27); this age was 308 compared with a sample from a polished surface which correlates with the internal moraine (Fig. 4C) that 309 returned an age of 13.0 ± 0.8 ka (ARAN-28), which agrees with the latter.

315

The ELAs of the three reconstructed paleoglacier extents were calculated for the Ruda catchment. The 316 paleoglacier existing at 13.5 ± 0.9 ka deposited a moraine system at 2080 m, with an average ELA located at 317 2461 m (Table 4). The paleoglacier at 13.0 ± 0.8 ka formed a moraine at 2190 m, with an average ELA at 2505 318 m. Finally, the Sendrosa Cirque moraine located at 2470 m was dated at 12.6 ± 1.3 ka, with an ELA of 2571 319 m (Table 4).

320

Table 4 321 5. Discussion 322

Interpretation of the geomorphological observations and CRE results

323

The geomorphological analysis supported by 10 Be dating of glacial landforms show the spatial and temporal 324 patterns of the deglaciation between the B-A and the early YD in the Ruda Valley (Fig. 5). Although the CRE 325 dataset is consistent with sequence of glacial landforms, uncertainty remains in assigning ages either to the 326 late B-A or to the early YD.

327

The ages of samples from polished surfaces at the valley at 1860-1900 m (ARAN-40 and ARAN-39; 14.4 ± 328 1.2 ka) suggest that the beginning of glacial retreat in the Ruda Valley took place during the B-A. As glaciers 329 retreated, they became confined within their respective cirques, as revealed by the age of 14.8 ± 0.9 ka (ARAN-330 31) of a polished surface on the ridge dividing the Sendrosa and Saboredo cirques (Fig. 3). The sample from 331 an outcrop on the lower edge of the cirque (ARAN-34; 14.0 ± 0.8 ka) indicates that ice covered the entire 332

Saboredo cirque floor until ~14 ka.

333

The sampled boulders of the external moraine system (ARAN-38 and ARAN-37) showed statistical coherence 334 based on the chi-squared test following Ward and Wilson (1978) and suggested a glacier advance or stillstand 335 at 13.5 ± 0.9 ka. The age of 13.0 ± 0.8 ka from a moraine boulder (ARAN-27) and an adjacent polished surface 336

(ARAN-28) of the internal moraine system indicate another short advance or stillstand, which is also consistent 337 with the geochronological sequence. These results from staggered moraine systems at ~2080 and ~2190 m 338 provided chronologically consistent mean ages within the B-A. But is based on their distinct relative position 339 that they have been ascribed to two different stages. Nevertheless, during this period, glaciers flowed from the 340 Saboredo and Sendrosa cirques and generated moraines on the lower overdeepened basins.

341

According to the ages of samples from the highest surfaces of the Saboredo Cirque floor at 2320-2360 m 342 (ARAN-35 and ARAN-36; 12.7 ± 0.8 ka), the entire cirque floor was deglaciated at the onset of the YD. These 343 ages are consistent with those of the samples from the highest dated moraine cirques in the Sendrosa Cirque 344 at 2460-2470 m (ARAN-32 and ARAN-33), which was also abandoned by the ice by 12.6 ± 1.3 ka. No 345 statistical inconsistence was detected in the exposure ages from this moraine.

346

Figure 5 347

The chronology of deglaciation 348

During the lLGM, the Garonne glacier flowed ~80 km reaching the Loures-Barousse-Barbazan basin, where 349 it formed a terminal moraine system (Fernandes et al., 2017;Stange et al., 2014) indicating that by the early B-A the Garonne glacier represented only ~5% of the length that it had during the 358 lLGM (Fernandes et al., 2017).

359

As glaciers retreated and thinned during the B-A, ice masses became individualized within their respective 360 cirques, and the valley glaciers flowing from Sendrosa and Saboredo cirques merged only in the lowest 361 overdeepened basins (Fig. 3). As conditions were less favourable, glaciers receded to higher elevations and the 362 Sendrosa and Saboredo paleoglaciers became isolated in their respective basins by 14.8 ± 0.9 ka (Fig. 3). The 363 ice covered the entire Saboredo cirque floor until ~14 ka, as revealed by the sample collected from an outcrop 364 on the edge of the cirque (14.0 ± 0.8 ka), indicating that glacial recession in the Ruda Valley was underway at 365 ~15-14 ka. Similarly, a rapid deglaciation was also recorded at ~15-14 ka in the nearby Bacivèr Cirque, which 366 favoured its complete deglaciation (Fig. 1; Oliva et al., 2021a).

367

The strong clustering of ages at between ~13.5 and 12.7 ka confirms the dynamic environmental response of 368 the Ruda palaeoglaciers to the abrupt climate transition between the B-A and the early YD. However, our CRE 369 dataset cannot assess the magnitude and patterns of glacial recession during the early B-A in the Ruda Valley, 370

or whether the entire cirque was deglaciated and glaciers reformed and subsequently advanced. Our results 371 suggest, however, that the millennial-scale trend towards glacier recession during the B-A was interrupted by 372 two glacial advances or stillstands with moraine development at 13.5 ± 0.9 and 13.0 ± 0.8 ka (Fig. 5). Samples 373 from moraines located at average elevations of ~2080 and ~2190 m, respectively, indicated the existence of 374 valley glaciers 3.5 and 2.5 km long during the middle of the B-A and during the transition from the B-A to the 375 YD. Indeed, the reconstructed glacial surface for both bases was of 4.4 and 3.8 km 2 , respectively, with a 376 maximum ice thickness of 300 and 200 m for each of these phases (Fig. 6a,b).

377

Figure 6 378 Evidence of glacial recession during the middle of the B-A was also recorded on polished surfaces from the 379

neighbouring Bacivèr Cirque at 13.5 ± 0.8 ka (Oliva et al., 2021a). Local topographical factors, such as the 380 shape of the cirque floor (convex versus concave areas) determined the pattern of glacial flow of thin ice 381 masses during the glacial advances/stillstands that occurred during the B-A. Glaciers flowed downwards from 382 the cirque through concave areas across the western slope, favouring the development of the moraine arcs and 383 crests in lower overdeepened basins (Fig. 3). It is likely that the topography played a major role in determining 384 the process and timing of the deglaciation in both cirques: the accumulation area above 2400 m is larger in the 385 Saboredo Cirque (552 ha) with respect to the Bacivèr Cirque (244 ha) and its prevailing aspect is also more 386 favourable for ice persistence (N vs W).

387

The Saboredo Cirque floor became ice free at the beginning of the YD, as revealed by the samples from 388 polished surfaces of the cirque floor at ~2350 m (average age = 12.7 ± 0.8 ka). The evidence of glacial striations 389 on the rock surface make unlikely the possibility of nuclide inheritance at the sampling sites. Glaciers might 390 thus also have persisted during the early YD in the highest hollows of the Saboredo Cirque, as was the case in 391 the northern cirque of Sendrosa Peak. Here, the cirque moraine at 2470 m, formed by a 0.4 km-long glacier 392

(with a surface of 0.7 km 2 and a thickness of 90 m; Fig. 6c), was abandoned at 12.6 ± 1.3 ka. Inside the highest 393 glacial cirques of the Ruda Valley, there are numerous relict rock glaciers and protalus lobes indicative of the 394 past occurrence of permafrost during their formation (Fernandes et al., 2018). As in most of the Pyrenees, these 395 are probably glacier-derived permafrost features associated with the paraglacial phase, when glacial shrinking 396 favoured wall debuttressing, a high debris supply and the subsequent burial of the residual ice masses (Andrés 397 et al., 2018;Oliva et al., 2016).

398

In short, by the early B-A the Garonne Paleoglacier had receded ~95% of its ~80 km length of the lLGM 399 (Fernandes et al., 2017;Stange et al., 2014). The Ruda Paleoglacier receded 4 km in only ~1.7 ka (~14.4 to 400 ~12.7 ka) from an elevation of ~1800 to >2400-2500 m. This trend towards glacial shrinking before the final 401 ice disappearance was interrupted by three phases of glacier advance or stillstands at 13.5 ± 0.9, 13.0 ± 0.8 and 402

12.6 ± 1.3 ka, which generated moraines at ~2080, 2190 and 2470 m, respectively. Iberian mountain ranges (García- Ruiz et al., 2010;Palacios and García-Ruiz, 2015). There is increasing 408 evidence that the timing of glacial oscillations during T-1 across Europe followed a similar sequence (Palacios 409 et al., 2021): a significant glacier advance took place during the OD followed by a rapid retreat during the B-410

A and a subsequent more limited glacier advance during the YD (Oliva et al., 2019).

411

Dated evidence of glacial responses in the Pyrenees to Northern Hemisphere cold events during T-1 have 412 mostly focused on the southern slope (Oliva et al., 2019), although there is an increasing body of new data 413 from the northern slope of the range that confirms glacial advances and stillstands there during the OD and 414 YD, and glacial shrinking during the B-A (Reixach et al., 2021). The elevation of the peaks as well as the 415 topography of the cirques played a major role in the duration and intensity of the glaciation during T-1 in the 416 high massifs of the Pyrenees. Cirques surrounded by peaks above 2800-2900 m likely included glaciers during 417 most of T-1 until the onset of the Holocene (see e.g. Jomelli et al., 2020;Pallàs et al., 2006;Reixach et al., 418 2021), while those situated at the foot of lower peaks (i.e. 2500-2600 m) were completely deglaciated during 419 the B-A (Oliva et al., 2021a;[START_REF] Oliva | Following ice retreat, the slopes were affected by intense 420 readjustments (large, catastrophic rock slope failures, rock slope deformations, rockfalls, and rockslides), 421 particularly in areas where the local lithostructure was prone to instability[END_REF]. Following ice retreat, the slopes were affected by intense 420 readjustments (large, catastrophic rock slope failures, rock slope deformations, rockfalls, and rockslides), 421

particularly in areas where the local lithostructure was prone to instability (Fernandes et al., 2020).

422

In the Central Pyrenees, a period of glacial growth occurred during the OD after the massive deglaciation that 423 took place following the LGM. Moraines dated at 17-15 ka have been reported in the highest valleys of the 424 main massifs, such as in the Ariège (Delmas et al., 2011;Jomelli et al., 2020;Reixach et al., 2021), Ésera 425 (Crest et al., 2017), Gállego (Palacios et al., 2015a) and Ésera valleys in the Central Pyrenees (Crest et al., 2017;[START_REF] Oliva | Following ice retreat, the slopes were affected by intense 420 readjustments (large, catastrophic rock slope failures, rock slope deformations, rockfalls, and rockslides), 421 particularly in areas where the local lithostructure was prone to instability[END_REF].

440

No moraine remnants have been found in the lowest and central sections of the Ruda Valley, and the prevailing 441 lithology of slates and marbles in this area is not favourable for the preservation of glacially abraded surfaces. 442

Indeed, the neighbouring Bacivèr Cirque was completely deglaciated also by ~15-14 ka (Oliva et al., 2021a). 443 Thus, our data cannot confirm whether small glaciers persisted in the highest parts of the Ruda Valley during 444 the B-A, although CRE ages support climatic conditions during the early B-A unfavourable for the 445 maintenance of the OD ice masses (Fig. 7). In other valleys of the Pyrenees, glaciers also receded rapidly at 446 ~14.6-14 ka as shown by CRE ages from polished surfaces and scattered erratic boulders in the Ariège (14.4 447 ± 1.1, 14.1 ± 0.7, 14.0 ± 0.8 ka; Delmas et al., 2011), Gállego (14.6 ± 2.3 ka;Palacios et al., 2015a) and Ésera 448 valleys (14.3 ± 0.5 ka; Crest et al., 2017). Locally, some massifs of the Pyrenees were fully deglaciated during 449 the B-A, with small glaciers persisting only in the highest cirques (Delmas, 2015). A similar pattern was 450 detected in other Iberian ranges, such as in the Sierra Nevada (Palacios et al., 2016) -Wormull et al., 2021). Marine records from the western Iberian margin also pointed to a 5-6ºC of 458 warming at the onset of the B-A (Martrat et al., 2007). These warmer conditions also favoured the recession 459 of glaciers across most European mountain ranges, such as in the Alps at 15.9-14.3 cal ka BP (Ivy-Ochs, 2015), 460

Tatra Mountains at 14.8-14.2 ka (Zasadni et al., 2020), Taurus Mountains at 14.6 ± 2.8 ka (Ciner and Sarikaya, 461 2015), as well as the Scandinavian Ice Sheet, which recorded a shrinking trend starting at 14.6 cal ka BP 462 (Mangerud et al., 2013).

463

However, the B-A was a climatically unstable period characterised by an alternation of warm and cold 464 temperatures (Rasmussen et al., 2014). The gradual decrease of temperatures from the early to the late B-A 465 was depicted by periods of enhanced cold conditions in the high latitudes of the Northern Hemisphere, which 466 also favoured periods of glacial growth in the Pyrenees (Table 5). Indeed, CRE results indicate the formation 467 of moraines in the Ruda Valley at 13.5 ± 0.9 and 13.0 ± 0.8 ka (Fig. 7). These CRE ages are consistent with 468 the chronostratigraphic sequence but must be taken as approximate due to their uncertainty ranges and the low 469 number of dated samples at each moraine. At that time, small ice tongues ~3 km long flowed downslope from 470 the main Saboredo Cirque as well as from other neighbouring glacial hollows to overdeepened basins, forming 471 moraines at 2080 and 2190 m, respectively. Whereas evidence of phases of glacier advance or stillstands during 472 the transition from the late B-A to the early YD were already found in other cirques in the Central Pyrenees, 473 such as in the Bacivèr Cirque (Oliva et al., 2021a), and in the Cuerpo de Hombre Valley (Sierra de Gredos, 474

Central Range), with moraine stabilization at 13.1 ka (Carrasco et al., 2015), this is the first study that reports 475 geochronological evidence of a phase with glacial development in the second half of the B-A after the long-476

term cooling initiated at ~14.5 ka (Fig. 7; Table 5). A few studies in Europe have also found evidence of glacier 477 advances or stillstands during the middle B-A, such as in the Tatra Mountains where moraines stabilized at 478

13.4 ± 0.5 ka (Engel et al., 2017). By contrast, geomorphic evidence of a period of glacial expansion during 479 the transition between the B-A to the early YD was detected in several regions, such as in the 480

Hardangerfjorden-Herdla area, western Norway, where a till deposit suggested a glacier advance of the 481

Scandinavian Ice sheet at 13.5-13.0 ka (Mangerud et al., 2016), and in the West Highland and Mull ice fields 482 of the Scottish Highlands, where 18 14 C dates from organic fragments found in basal tills revealed a glacial 483 readvance from the end of the B-A or earliest YD (Bromley et al., 2018), although these results are considered 484 controversial (Small and Fabel, 2016).

485

Table 5 486

Based on the position of these two moraine ridges, the ELA in the Ruda Valley glacier must have been located 487 at 2461 and 2505 m during the glacial advances of the second part of the B-A and the late B-A or early YD, 488 respectively. Considering that the regional ELA in the Pyrenees currently lies at approximately 3100 m 489 (Jomelli et al., 2020;René, 2011), the ELA must have been located 639-595 m lower than present-day 490 conditions. Assuming an average lapse rate of 0.65 °C 100 m -1 and no change in precipitation, we infer that 491 summer temperatures must have been 4.2 and 3.9 ºC lower than present-day during these two phases (Table 492 4).

493

Figure 7 494

The temperature decrease during the early YD favoured the persistence of glaciers inside the highest hollows 495 in the Ruda Valley, depositing moraines at ~2400-2500 m, which were later abandoned by the ice at 12.6 ± 496 1.3 ka (Fig 7). The existence of glaciers at high altitudes of the Saboredo Cirque is also confirmed by polished 497 surfaces at elevations of 2200-2300 m that were dated at 12.7 ± 0.8 ka. In the neighbouring Bacivèr Cirque, 498 moraines at 2400 m were abandoned by the ice at 12.8 ka (Oliva et al., 2021a). The elevation difference 499 between the ELA during the YD (2571 m) in the Ruda Valley and present-day values indicates that summer 500 temperatures were ~3.4 ºC lower during that cold phase (Table 4), which is very similar to the ~3.0 ºC inferred 501 for the Bacivèr Cirque (Oliva et al., 2021a). Indeed, topoclimatic factors must have played a key role 502 maintaining small glaciers during the YD at the foot of steep cirque walls, where shading and snow 503 accumulation allowed the longer persistence of reduced ice masses (Boston and Lukas, 2019). In the Pyrenees, 504 YD glacial advances have been detected in cirques of the Central (Delmas, 2009;Delmas et al., 2008;García-505 Ruiz et al., 2016a) and Eastern sectors (Crest et al., 2017;Jomelli et al., 2020;Pallàs et al., 2010;Reixach et 506 al., 2021), where 2-3.5 km-long glaciers formed at 12.9-11.2 ka. Based on chronostratigraphic data from a 507 variety of sources, evidence of YD glacial advances was also recorded in the highest cirques of other Iberian 508 massifs, although no direct ages are yet available (García- Ruiz et al., 2016a). This is true for the NW ranges 509 (Cowton et al., 2009) and the Cantabrian Mountains (Pellitero et al., 2019;Serrano et al., 2015), where 510 moraines above 1800 and 2000 m, respectively, were associated with this cold phase. In the Sierra Nevada,

511

YD moraines dated at 13-12 ka are located in the highest cirques facing east from 2350 (north face) to 2800 m 512 (south face) (Oliva et al., 2014;Palacios et al., 2016). These glacial advances were driven by a temperature 513 cooling quantified at 2-4 ºC in western Europe that was triggered by a slowdown of the overturning circulation 514 in the North Atlantic (Renssen et al., 2018). In the Pyrenees, the YD was the coldest period of the last 515 deglaciation recorded by a decrease of 18 O in the Ostolo and Seso speleothems that might reflect the 516 temperature decrease of ~5ºC (Bartolomé et al., 2015;Bernal-Wormull et al., 2021). Colder conditions also 517 favoured glacial growth in other mountainous regions of Europe, with widespread advances across the Alps 518 between 13.5 and 12 ka (Ivy-Ochs, 2015;Ivy-Ochs et al., 2009), in the Tatra Mountains where moraines 519 stabilized at 11.9 ± 0.5 ka (Engel et al., 2017) and in Turkey where they formed at 12.6 ± 2.3 ka (Ciner and 520 Sarikaya, 2015).

521

The final deglaciation of the cirques favoured paraglacial dynamics and the formation of the numerous rock 522 glaciers existing today in the Saboredo and Sendrosa cirques, as well as in the lateral cirques distributed along 523 the Ruda Valley (Fig. 3). These features are currently inactive under the present-day climate regime (Fernandes 524 et al., 2018), with abundant lichen coverage and forest colonization, particularly on the lowest ridges. In the 525

Pyrenees, as well as in other Iberian ranges, cirque wall readjustment following deglaciation triggered different 526 landforms, such as debris-covered glaciers or rock glaciers, that remained active until the early Holocene 527 (Oliva et al., 2019). This is the case of the Bacivèr cirque, where the debris-covered glacier that formed inside 528 the highest YD moraines stabilized at 7.2 ± 0.5 ka (Oliva et al., 2021a), or of several other massifs in the 529

Central and Eastern Pyrenees where rock glaciers were active until well into the Holocene [START_REF] Andrés | Timing of deglaciation and rock glacier origin in the southeastern Pyrenees: a 594 review and new data[END_REF]530 García-Ruiz et al., 2016a). A similar pattern was also detected in other Iberian massifs, such as in the Iberian 531

Range [START_REF] Aumaître | Chronological and geomorphological investigation of fossil debris-602 covered glaciers in relation to deglaciation processes: A case study in the Sierra de La Demanda, 603 northern Spain[END_REF] or Sierra Nevada (Palacios et al., 2016), where these paraglacial landforms 532 formed immediately after the deglaciation of the cirques, and remained active until the complete melting of 533 their relict ice during the Holocene Thermal Maximum. 534

Conclusions

535

The Central Pyrenees were extensively glaciated during Quaternary glacial phases and were only slightly 536 transformed by periglacial, nival and slope processes during subsequent interglacial periods such as the 537

Holocene. As a result, the landscape in the Central Pyrenees includes a wide range of geomorphological 538 features that constitutes a prominent natural heritage that contributed to the designation of most protected areas 539 in this mountain range. This is the case of the Ruda Valley, located on the periphery of the National Park of 540

Aigüestortes and Sant Maurici Lake. The valley is one of the most scenic glacial valleys of the Pyrenees, and 541 includes a wide variety of glacial and periglacial landforms, particularly above 1800 m. However, in contrast 542 to the southern slope of the Central Pyrenees where the timing of the glacial evolution has been widely studied, 543 on the northern side there are still temporal and spatial gaps in the chronology of glacial activity. In this work, 544

we determined the magnitude and timing of T-1 glacial oscillations in the Upper Garonne Basin by producing 545

12 new CRE ages of glacial landforms (moraine boulders and glacially polished surfaces) that provide evidence 546 of the spatio-temporal pattern of glacial oscillations from the B-A to the YD.

547

The Upper Garonne Basin hosted the largest glacial system of the Pyrenees during the last glacial cycle thanks 548 to its Atlantic-influenced climate and the high altitude of the peaks in the central part of this range. The long-549 term post-LGM warming favoured glacial recession, particularly during the early B-A, and by 14.4 ± 1.2 ka 550 the Ruda Paleoglacier was already disconnected from the main Garonne Paleoglacier. At that time, the glacier 551 front was located at 1860-1900 m, and some areas of the Saboredo Cirque were already ice-free as a result of 552 ice thinning. The B-A was also a period of great climatic variability, and the small glaciers in the Ruda Valley 553

showed rapid response times with formation of moraines during the coldest B-A phases in north Iberia. Phases 554 of glacier advance or stillstands promoted the development of moraines at 13.5 ± 0.9 ka, as well as during the 555 transition between the late B-A and the onset of the YD at 13.0 ± 0.8 ka. The YD probably favoured the 556 persistence of small glaciers during an early stage of this cold phase; there is evidence of glacial retreat at the 557 uppermost part of the cirque floor at 12.7 ± 0. ARAN-40 15.0 ± 0.9 (0.5) 15.2 ± 0.9 (0.5) 16.1 ± 1.0 (0.5) 16.4 ± 1.0 (0.6) ARAN-39 13.8 ± 0.8 (0.4) 14.0 ± 0.8 (0.4) 14.8 ± 0.9 (0.4) 15.0 ± 0.9 (0.5)

Valley -external moraines 13.5 ± 0.9 13.7 ± 0.9 14.5 ± 1.0 14.7 ± 1.0 ARAN-38 13.3 ± 0.8 (0.4) 13.5 ± 0.8 (0.4) 14.4 ± 0.9 (0.4) 14.6 ± 0.9 (0.4) ARAN-37 13.7 ± 1.0 (0.6) 13.8 ± 1.0 (0.7) 14.7 ± 1.0 (0.7) 14.9 ± 1.0 (0.7)

Valley -internal moraines 13.0 ± 0.8 13.2 ± 0.8 14.1 ± 0.9 14.2 ± 0.9 ARAN-28 13.0 ± 0.8 (0.4) 13.1 ± 0.8 (0.4) 14.0 ± 0.9 (0.5) 14.2 ± 0.9 (0.5) ARAN-27 13.1 ± 0.8 (0.4) 13.3 ± 0.8 (0.4) 14.1 ± 0.8 (0.4) 14.3 ± 0.9 (0.4)

Cirque -polished surfaces 12.7 ± 0.7 12.9 ± 0.8 13.7 ± 0.8 13.9 ± 0.8 ARAN-36 12.7 ± 0.8 (0.4) 12.9 ± 0.8 (0.4) 13.7 ± 0.8 (0.4) 13.9 ± 0.8 (0.4) ARAN-35 12.8 ± 0.8 (0.4) 12.9 ± 0.8 (0.4) 13.8 ± 0.8 (0.4) 14.0 ± 0.8 (0.4) ARAN-34 14.0 ± 0.8 (0.4) 14.2 ± 0.9 (0.4) 15.1 ± 0.9 (0.5) 15.3 ± 0.9 (0.5) ARAN-31 14.8 ± 0.9 (0.4) 15.1 ± 0.9 (0.4) 16.0 ± 0.9 (0.4) 16.2 ± 0.9 (0.4)

Cirque -moraines 12.6 ± 1.3 12.7 ± 1.4 13.5 ± 1.4 13.7 ± 1.5

ARAN-33 11.8 ± 0.7 (0.4) 12.0 ± 0.8 (0.4) 12.7 ± 0.8 (0.4) 12.9 ± 0.8 (0.4) ARAN-32 13.3 ± 1.0 (0.7) 13.5 ± 1.0 (0.7) 14.3 ± 1.0 (0.7) 14.5 ± 1.0 (0.7)

Phase (ages in ka)

Climate trend in Greenland based on ice core records (Rasmussen et al., 2014)

Glacier response in the Ruda Valley

Climate in the Iberian Peninsula based on other independent proxies Abrupt warming (3-5ºC) GI-1d (14.1-14.0) Short cold reversal (1-2ºC) Rapid recession at ~15-14 ka

Gradual warming in the northern Iberia at 15.4-13.4 cal ka BP (Moreno et al., 2010). Warming conditions in the Pyrenees at 14.7-14.0 cal ka BP (González-Sampériz et al., 2006). Significant warming at the beginning (2-4ºC) followed by a gradual cooling at (~2ºC). Intense cold conditions at 13.7-13.6 ka Glacier advance or stillstand with moraine formation at 13.5 ka

Colder conditions starting at 13.5 cal ka BP in northern Iberia (Moreno et al., 2010). Cold climate in the Pyrenees at 14-13.4 cal ka BP (González-Sampériz et al., 2006).

GI-1b (13.3-13.1) Short cold reversal (1-2ºC) -Short periods of colder and drier conditions (Moreno et al., 2010). Expansion of coldclimate forest species in the Pyrenees at 13.4-13.1 cal ka BP (González-Sampériz et al., 2006).

GI-1a (13.1-12.9) The Pyrenees were the most extensively glaciated mountain system during the maximum ice extent of the last 94 glaciation in Iberia, hosting ~50% of the peninsula's glacial surface (Oliva et al., 2021b). The chronology of 95 glacier fluctuations during T-1 seems to have followed a similar temporal pattern, although available data is 96 principally derived from the range's southern slope [START_REF] Andrés | Timing of deglaciation and rock glacier origin in the southeastern Pyrenees: a 594 review and new data[END_REF]Palacios et al., 2015b[START_REF] Oliva | Following ice retreat, the slopes were affected by intense 420 readjustments (large, catastrophic rock slope failures, rock slope deformations, rockfalls, and rockslides), 421 particularly in areas where the local lithostructure was prone to instability[END_REF]Palacios et al., , 2017a)).

97

Pyrenean glaciers readvanced to the valley bottoms during the OD, with ice tongues > 15 km long. During the 98 B-A, by 15-14 ka, high temperatures favoured rapid deglaciation and glaciers disappeared from relatively low 99 cirques shaped under peaks of ~2600 m a.s.l. (Oliva et al., 2021a). A subsequent readvance occurred during 100 the YD, with smaller glaciers of only ~4 km length confined within the cirques of the highest massifs above 101 2200 m (Crest et al., 2017;García-Ruiz et al., 2016a;Jomelli et al., 2020;[START_REF] Oliva | Following ice retreat, the slopes were affected by intense 420 readjustments (large, catastrophic rock slope failures, rock slope deformations, rockfalls, and rockslides), 421 particularly in areas where the local lithostructure was prone to instability[END_REF]Palacios et al., , 2015a;;Pallàs 102 et al., 2006). Following both cold periods, glacial retreat triggered intense paraglacial dynamics at ~15-14 ka 103 and ~11-10.5 ka, respectively, which led to the formation of rock glaciers and debris-covered glaciers, 104

particularly in cirques below 2800 m [START_REF] Andrés | Timing of deglaciation and rock glacier origin in the southeastern Pyrenees: a 594 review and new data[END_REF][START_REF] Oliva | Following ice retreat, the slopes were affected by intense 420 readjustments (large, catastrophic rock slope failures, rock slope deformations, rockfalls, and rockslides), 421 particularly in areas where the local lithostructure was prone to instability[END_REF].

105

The Holocene. Another recent study that focused on the Bacivèr Cirque, near our study area, revealed a rapid 111 deglaciation during the B-A, followed by glacial readvance and the formation of the highest moraines at the 112 B-A/YD transition at ~12.9 ka; subsequently, intense paraglacial readjustment favoured the formation of 113 debris-covered glaciers and rock glaciers inside these moraine systems (Oliva et al., 2021a).

114

In order to complement our understanding of the glacial history of the Pyrenees, this study uses new data to 115 improve the deglaciation chronology of the Upper Garonne Basin, the largest glacial system in this mountain 116 range, during the local last glacial maximum (lLGM; Fernandes et al., 2017). Glacial retreat in the Ruda Valley, 117

where this work focuses on, left widespread geomorphic evidence of glacial and paraglacial origin that can be 118 used to reconstruct the last stages of long-term deglaciation in the Central Pyrenees. The Pyrenees stretch E-W over 400 km across north Iberia, with the highest peaks (> 3000 m) in the central 130 part of the range. This study focuses on the Ruda Valley, located between latitudes 42º40' N -42º42' N and 131 longitudes 0º57' E -1º00' E (Fig. 1), that constitutes the headwaters of the Upper Garonne Basin. The area is 132 included in the peripheral protection area of the National Park of Aigüestortes and Sant Maurici Lake.

133 Figure 1 134

The Upper Garonne Valley shows the typical U-shaped cross-section of glacial valleys, and the Garonne River 135 flows to the N-NW through the central axis of the Pyrenees. Several tributary rivers drain the glacial valleys 136 descending from glacial cirques carved below peaks from 2800 to 3000 m. Glacial cirques are prevailingly 137 NE-oriented, with floors generally between 2200 and 2400 m (Lopes et al., 2018).

138

The 2019, when snow cover was limited to the concave and sheltered areas at the highest elevations (>2500 m), 167 which enabled better identification of the glacial landforms, as well as sampling for CRE dating.

168

Geomorphological mapping 169

We produced a geomorphological map focusing on glacial and periglacial landforms based on recent general 170 geomorphological maps from the area (Fernandes et al., 2017;Fernandes et al., 2021). Here we used the same 171 approach based on the combination of aerial orthophotographs (0.25 m cell size) and high-resolution LIDAR 172 digital elevation models (point density of 0.5-2 m -2 ) from the "Institut Cartogràfic i Geològic de Catalunya" 173 (http://www.icc.cat/appdownloads). The symbology of geomorphological landforms followed Joly (1997).

174

Field strategy and sampling 175

We collected 12 granite and granodiorite samples for CRE dating, using a hammer and chisel, from elevations 176 spanning from 1862 to 2470 m (Table 1). To ensure that samples were intact since glacial retreat, we targeted 177 at well-anchored moraine boulders and glacially polished bedrock that stood out within their surroundings, 178 which minimizes the risk of having been buried under sediments (Fig. 2). To ensure the best cosmic-ray flux 179 reception of the sampling sites, we restricted sample extraction to flat-topped and gentle surfaces, and sampled 180 a thickness of up to 3-5 cm. Corrections for partial shielding due to the surrounding topography were 181 implemented using the ArcGIS toolbox by Li (2018), which applies equations describing the effects of 182 geometric shielding on the sampling sites and a model to predict the nuclide production rate at depth in sloped 183 surfaces (Dunne et al. 1999); it required only a point shapefile of the sampling sites including strike/dip of the 184 sampling surfaces and a digital elevation model of 5 m spatial resolution.

185 Table 1 186 Figure 2 3. with the quartz-rich lithology of the samples, they were treated for the extraction of the in-situ-produced 193 cosmogenic nuclide 10 Be.

194

First, magnetic minerals were discarded through magnetic separation conducted in a "Frantz LB-1" magnetic 195 separator. Once the non-magnetic fraction was isolated, it underwent several rounds of chemical attacks with 196 a 1:2 mixture of concentrated hydrochloric (HCl) and hexafluorosilicic (H 2 SiF 6 , 2/3) acids to dissolve and 197 discard the non-quartz minerals. The remaining minerals (i.e. non-dissolved feldspar minerals) were 198 subsequently attacked with 4-6 successive partial dissolutions with HF to remove the atmospheric 10 Be and 199 discard the remaining impurities. In order to be sure that the quartz was pure, samples were examined under a 200 binocular microscope. As a result, ~20 g of pure quartz per sample was used to extract the 10 Be (Table 2). After 201 adding ~150 µL of an in-house manufactured (from a phenakite crystal) 9 Be carrier solution (spike, 202 concentration: 3025 ± 9 µg g 1 ; Merchel et al., 2008), the purified quartz was dissolved with 48% concentrated 203 HF (3.5 mL per g of quartz + 30 mL in excess). Following total dissolution, the resulting solution was 204 evaporated until dryness, and the residue was recovered with HCl (7.1 molar). Be was then precipitated at 205 pH=8 to beryllium hydroxide (Be(OH) 2 ) by means of ammonia (NH 3 ), and separated from other elements in 206 resin columns: a Dowex 1x8 anionic exchange column to remove elements such as Fe, Mn and Ti, and a 207 Dowex 50W x8 cationic exchange column to discard B and separate Be from Al (Merchel and Herpers, 1999).

208

The final eluted Be was again precipitated, dried and oxidized to BeO at 700 ºC. Finally, the targets for the 209 accelerator mass spectrometer (AMS) measurements were prepared by mixing with niobium powder at an 210 approximate 1:1 proportion and pressing the mixture into copper cathodes.

211

The targets were analysed at the 'Accelerator pour les Sciences de la Terre, Environnement et Risques' 212

(ASTER) national AMS facility at CEREGE in order to measure the 10 Be/ 9 Be ratio from which the 10 Be 213 concentration was later inferred (Table 2). AMS measurements were calibrated against the in-house standard 214 STD-11 with an assigned 10 Be/ 9 Be ratio of (1.191 ± 0.013) × 10 11 (Braucher et al., 2015). Based on calibration 215 with the previous ASTER AMS standard NIST-27900 (SRM4325), the new STD-11 is fully compatible with 216 07KNSTD. The analytical 1σ uncertainties include uncertainties in the AMS counting statistics, an external 217 0.5% AMS error [START_REF] Arnold | The French accelerator mass spectrometry facility ASTER: Improved performance 597 and developments[END_REF] and the uncertainty related to the chemical blank correction. The 10 Be 218 half-life considered was (1.3870 ± 0.0012) × 10 6 years (Chmeleff et al., 2010;Korschinek et al., 2010).

219

Table 2 220

We calculated 10 Be exposure ages by using the CREp online calculator (Martin et al., 2017; available online 221 at: http://crep.crpg.cnrs-nancy.fr/#/), where we selected the following settings: LSD (Lifton-Sato-Dunai) 222 elevation/latitude scaling scheme (Lifton et al., 2014), ERA40 atmospheric model (Uppala et al., 2005) and 223 geomagnetic database based on the LSD framework (Lifton et al., 2014). These settings yielded a world-wide 224 mean sea level high latitude (SLHL) 10 Be production rate of 3.98 ± 0.22 atoms g 1 yr 1 . Exposure ages of the 225 samples with their 1σ full and analytical uncertainties are shown in Table 2. The uncertainties discussed 226

throughout the text include analytical and production rate error unless otherwise stated.

227

In order to evaluate the potential impact of erosion on the exposure ages, we assumed a steady erosion rate (1 228 mm ka -1 ) following a conservative maximum value in such a lithological setting according to [START_REF] André | Rates of Postglacial rock weathering on glacially scoured outcrops (Abisko-Riksgränsen 590 area, 68°N)[END_REF].

229

The impact of snow cover was explored assuming that snow currently remains in the area an annual average 230 of 7.5 months at 2200-2300 m with a mean thickness of 100 cm -and assuming a snow density of 0.2 g cm -3 231 (Styllas et al., 2018) -(Table 3), based on the available data from the Bonaigua weather station (Fig. (Sendrosa and Saboredo peaks) to 2525 m in the mountain divide that may have functioned as a glacial 281 difluence pass between the Saboredo (Garonne basin) and the Ratera cirques (Noguera Pallaresa basin) 282 (Fig. 3). Saboredo is a complex cirque [START_REF] Barr | Glacial cirques as palaeoenvironmental indicators: Their potential and 607 limitations[END_REF] that includes several minor hollows 283 carved out on the marginal flanks of the peaks. They are separated by truncated spurs with polished surfaces, 284

where one sample (ARAN-31) was collected to examine the timing of glacial thinning at the cirque level.

285

The highest moraines in the Ruda Valley are located at ~2400-2500 m, at the foot of the highest peaks. A 286 prominent cirque moraine at an altitude of 2470 m is located on the eastern slope of Sendrosa Peak; here, 287 we collected two samples (ARAN-33 and ARAN-32). Above the upper moraines, fed by abundant debris 288 supply from the cirque rock walls, rock glaciers and protalus lobes are found. The cirque floors extend 289 between 2200 and 2400 m and preserve features indicative of intense glacial erosion (e.g. overdeepened 290

basins or glacial thresholds) and glacial deposits (e.g. moraines or erratic boulders). In the Saboredo Cirque, 291 glacial abrasion and periglacial processes shaped a large amphitheatre with several overdeepened basins.

292

They are surrounded by polished surfaces and occupied by peatlands and lakes. Three samples were 293 collected from polished surfaces from 2360 to 2270 m (Fig. 3) (ARAN-36; ARAN-35; ARAN-34; Fig. 2; 294 Table 1). Glacial retreat left scattered boulders and till across the area.

295

Figure 3 296 Figure 4 297

Geochronological data 298

The 12 CRE samples collected from the Ruda Valley returned a chronology spanning from ~15 ± 0.9 to 11.8 299 ± 0.7 ka, with all ages being consistent with the geomorphological distribution of the dated surfaces (Fig. 4).

300

Polished surfaces from the highest parts of the Ruda Valley located ~200-250 m below the cirques provided 301 CRE ages of 15.0 ± 0.9 (ARAN-40) and 13.8 ± 0.8 ka (ARAN-39, Fig. 4A) (average: 14.4 ± 1.2 ka). A 302 similar age was obtained from a sample collected from a truncated spur dividing the Sendrosa and Saboredo 303 cirques, which yielded 14.8 ± 0.9 ka (ARAN-31). A sample collected from a polished surface at the edge of 304 the Saboredo cirque floor (Fig. 3) reported an age of 14.0 ± 0.8 ka (ARAN-34, Fig. 4D).

305

Moraine boulders from the external moraine located in the lower overdeepened basin yielded ages of 13.3 ± 306 0.8 (ARAN-38) and 13.7 ± 1.0 ka (ARAN-37) (average: 13.5 ± 0.9 ka). A boulder from the internal moraine 307 ridge located 110 m above the former moraine provided an age of 13.1 ± 0.8 ka (ARAN-27); this age was 308 compared with a sample from a polished surface which correlates with the internal moraine (Fig. 4C) that 309 returned an age of 13.0 ± 0.8 ka (ARAN-28), which agrees with the latter.

315

The ELAs of the three reconstructed paleoglacier extents were calculated for the Ruda catchment. The 316 paleoglacier existing at 13.5 ± 0.9 ka deposited a moraine system at 2080 m, with an average ELA located at 317 2461 m (Table 4). The paleoglacier at 13.0 ± 0.8 ka formed a moraine at 2190 m, with an average ELA at 2505 318 m. Finally, the Sendrosa Cirque moraine located at 2470 m was dated at 12.6 ± 1.3 ka, with an ELA of 2571 319 m (Table 4).

320

Table 4 321 5. Discussion 322

Interpretation of the geomorphological observations and CRE results

323

The geomorphological analysis supported by 10 Be dating of glacial landforms show the spatial and temporal 324 patterns of the deglaciation between the B-A and the early YD in the Ruda Valley (Fig. 5). Although the CRE 325 dataset is consistent with sequence of glacial landforms, uncertainty remains in assigning ages either to the 326 late B-A or to the early YD.

327

The ages of samples from polished surfaces at the valley at 1860-1900 m (ARAN-40 and ARAN-39; 14.4 ± 328 1.2 ka) suggest that the beginning of glacial retreat in the Ruda Valley took place during the B-A. As glaciers 329 retreated, they became confined within their respective cirques, as revealed by the age of 14.8 ± 0.9 ka (ARAN-330 31) of a polished surface on the ridge dividing the Sendrosa and Saboredo cirques (Fig. 3). The sample from 331 an outcrop on the lower edge of the cirque (ARAN-34; 14.0 ± 0.8 ka) indicates that ice covered the entire 332

Saboredo cirque floor until ~14 ka.

333

The sampled boulders of the external moraine system (ARAN-38 and ARAN-37) showed statistical coherence 334 based on the chi-squared test following Ward and Wilson (1978) and suggested a glacier advance or stillstand 335 at 13.5 ± 0.9 ka. The age of 13.0 ± 0.8 ka from a moraine boulder (ARAN-27) and an adjacent polished surface 336

(ARAN-28) of the internal moraine system indicate another short advance or stillstand, which is also consistent 337 with the geochronological sequence. These results from staggered moraine systems at ~2080 and ~2190 m 338 provided chronologically consistent mean ages within the B-A. But is based on their distinct relative position 339 that they have been ascribed to two different stages. Nevertheless, during this period, glaciers flowed from the 340 Saboredo and Sendrosa cirques and generated moraines on the lower overdeepened basins.

341

According to the ages of samples from the highest surfaces of the Saboredo Cirque floor at 2320-2360 m 342 (ARAN-35 and ARAN-36; 12.7 ± 0.8 ka), the entire cirque floor was deglaciated at the onset of the YD. These 343 ages are consistent with those of the samples from the highest dated moraine cirques in the Sendrosa Cirque 344 at 2460-2470 m (ARAN-32 and ARAN-33), which was also abandoned by the ice by 12.6 ± 1.3 ka. No 345 statistical inconsistence was detected in the exposure ages from this moraine.

346

Figure 5 347

The chronology of deglaciation 348

During the lLGM, the Garonne glacier flowed ~80 km reaching the Loures-Barousse-Barbazan basin, where 349 it formed a terminal moraine system (Fernandes et al., 2017;Stange et al., 2014) indicating that by the early B-A the Garonne glacier represented only ~5% of the length that it had during the 358 lLGM (Fernandes et al., 2017).

359

As glaciers retreated and thinned during the B-A, ice masses became individualized within their respective 360 cirques, and the valley glaciers flowing from Sendrosa and Saboredo cirques merged only in the lowest 361 overdeepened basins (Fig. 3). As conditions were less favourable, glaciers receded to higher elevations and the 362 Sendrosa and Saboredo paleoglaciers became isolated in their respective basins by 14.8 ± 0.9 ka (Fig. 3). The 363 ice covered the entire Saboredo cirque floor until ~14 ka, as revealed by the sample collected from an outcrop 364 on the edge of the cirque (14.0 ± 0.8 ka), indicating that glacial recession in the Ruda Valley was underway at 365 ~15-14 ka. Similarly, a rapid deglaciation was also recorded at ~15-14 ka in the nearby Bacivèr Cirque, which 366 favoured its complete deglaciation (Fig. 1; Oliva et al., 2021a).

367

The strong clustering of ages at between ~13.5 and 12.7 ka confirms the dynamic environmental response of 368 the Ruda palaeoglaciers to the abrupt climate transition between the B-A and the early YD. However, our CRE 369 dataset cannot assess the magnitude and patterns of glacial recession during the early B-A in the Ruda Valley, 370 or whether the entire cirque was deglaciated and glaciers reformed and subsequently advanced. Our results 371 suggest, however, that the millennial-scale trend towards glacier recession during the B-A was interrupted by 372 two glacial advances or stillstands with moraine development at 13.5 ± 0.9 and 13.0 ± 0.8 ka (Fig. 5). Samples 373 from moraines located at average elevations of ~2080 and ~2190 m, respectively, indicated the existence of 374 valley glaciers 3.5 and 2.5 km long during the middle of the B-A and during the transition from the B-A to the 375 YD. Indeed, the reconstructed glacial surface for both bases was of 4.4 and 3.8 km 2 , respectively, with a 376 maximum ice thickness of 300 and 200 m for each of these phases (Fig. 6a,b).

377

Figure 6 378 Evidence of glacial recession during the middle of the B-A was also recorded on polished surfaces from the 379

neighbouring Bacivèr Cirque at 13.5 ± 0.8 ka (Oliva et al., 2021a). Local topographical factors, such as the 380 shape of the cirque floor (convex versus concave areas) determined the pattern of glacial flow of thin ice 381 masses during the glacial advances/stillstands that occurred during the B-A. Glaciers flowed downwards from 382 the cirque through concave areas across the western slope, favouring the development of the moraine arcs and 383 crests in lower overdeepened basins (Fig. 3). It is likely that the topography played a major role in determining 384 the process and timing of the deglaciation in both cirques: the accumulation area above 2400 m is larger in the 385 Saboredo Cirque (552 ha) with respect to the Bacivèr Cirque (244 ha) and its prevailing aspect is also more 386 favourable for ice persistence (N vs W).

387

The Saboredo Cirque floor became ice free at the beginning of the YD, as revealed by the samples from 388 polished surfaces of the cirque floor at ~2350 m (average age = 12.7 ± 0.8 ka). The evidence of glacial striations 389 on the rock surface make unlikely the possibility of nuclide inheritance at the sampling sites. Glaciers might 390 thus also have persisted during the early YD in the highest hollows of the Saboredo Cirque, as was the case in 391 the northern cirque of Sendrosa Peak. Here, the cirque moraine at 2470 m, formed by a 0.4 km-long glacier 392

(with a surface of 0.7 km 2 and a thickness of 90 m; Fig. 6c), was abandoned at 12.6 ± 1.3 ka. Inside the highest 393 glacial cirques of the Ruda Valley, there are numerous relict rock glaciers and protalus lobes indicative of the 394 past occurrence of permafrost during their formation (Fernandes et al., 2018). As in most of the Pyrenees, these 395 are probably glacier-derived permafrost features associated with the paraglacial phase, when glacial shrinking 396 favoured wall debuttressing, a high debris supply and the subsequent burial of the residual ice masses (Andrés 397 et al., 2018;Oliva et al., 2016).

398

In short, by the early B-A the Garonne Paleoglacier had receded ~95% of its ~80 km length of the lLGM 399 (Fernandes et al., 2017;Stange et al., 2014). The Ruda Paleoglacier receded 4 km in only ~1.7 ka (~14.4 to 400 ~12.7 ka) from an elevation of ~1800 to >2400-2500 m. This trend towards glacial shrinking before the final 401 ice disappearance was interrupted by three phases of glacier advance or stillstands at 13.5 ± 0.9, 13.0 ± 0.8 and 402

12.6 ± 1.3 ka, which generated moraines at ~2080, 2190 and 2470 m, respectively. Iberian mountain ranges (García- Ruiz et al., 2010;Palacios and García-Ruiz, 2015). There is increasing 408 evidence that the timing of glacial oscillations during T-1 across Europe followed a similar sequence (Palacios 409 et al., 2021): a significant glacier advance took place during the OD followed by a rapid retreat during the B-410

A and a subsequent more limited glacier advance during the YD (Oliva et al., 2019).

411

Dated evidence of glacial responses in the Pyrenees to Northern Hemisphere cold events during T-1 have 412 mostly focused on the southern slope (Oliva et al., 2019), although there is an increasing body of new data 413 from the northern slope of the range that confirms glacial advances and stillstands there during the OD and 414 YD, and glacial shrinking during the B-A (Reixach et al., 2021). The elevation of the peaks as well as the 415 topography of the cirques played a major role in the duration and intensity of the glaciation during T-1 in the 463 However, the B-A was a climatically unstable period characterised by an alternation of warm and cold 464 temperatures (Rasmussen et al., 2014). The gradual decrease of temperatures from the early to the late B-A 465 was depicted by periods of enhanced cold conditions in the high latitudes of the Northern Hemisphere, which 466 also favoured periods of glacial growth in the Pyrenees (Table 5). Indeed, CRE results indicate the formation 467 of moraines in the Ruda Valley at 13.5 ± 0.9 and 13.0 ± 0.8 ka (Fig. 7). These CRE ages are consistent with 468 the chronostratigraphic sequence but must be taken as approximate due to their uncertainty ranges and the low 469 number of dated samples at each moraine. At that time, small ice tongues ~3 km long flowed downslope from 470 the main Saboredo Cirque as well as from other neighbouring glacial hollows to overdeepened basins, forming 471 moraines at 2080 and 2190 m, respectively. Whereas evidence of phases of glacier advance or stillstands during 472 the transition from the late B-A to the early YD were already found in other cirques in the Central Pyrenees, 473 such as in the Bacivèr Cirque (Oliva et al., 2021a), and in the Cuerpo de Hombre Valley (Sierra de Gredos, 474

Central Range), with moraine stabilization at 13.1 ka (Carrasco et al., 2015), this is the first study that reports 475 geochronological evidence of a phase with glacial development in the second half of the B-A after the long-476

term cooling initiated at ~14.5 ka (Fig. 7; Table 5). A few studies in Europe have also found evidence of glacier 477 advances or stillstands during the middle B-A, such as in the Tatra Mountains where moraines stabilized at 478

13.4 ± 0.5 ka (Engel et al., 2017). By contrast, geomorphic evidence of a period of glacial expansion during 479 the transition between the B-A to the early YD was detected in several regions, such as in the 480

Hardangerfjorden-Herdla area, western Norway, where a till deposit suggested a glacier advance of the 481

Scandinavian Ice sheet at 13.5-13.0 ka (Mangerud et al., 2016), and in the West Highland and Mull ice fields 482 of the Scottish Highlands, where 18 14 C dates from organic fragments found in basal tills revealed a glacial 483 readvance from the end of the B-A or earliest YD (Bromley et al., 2018), although these results are considered 484 controversial (Small and Fabel, 2016).

485

Table 5 486

Based on the position of these two moraine ridges, the ELA in the Ruda Valley glacier must have been located 487 at 2461 and 2505 m during the glacial advances of the second part of the B-A and the late B-A or early YD, 488 respectively. Considering that the regional ELA in the Pyrenees currently lies at approximately 3100 m 489 (Jomelli et al., 2020;René, 2011), the ELA must have been located 639-595 m lower than present-day 490 conditions. Assuming an average lapse rate of 0.65 °C 100 m -1 and no change in precipitation, we infer that 491 summer temperatures must have been 4.2 and 3.9 ºC lower than present-day during these two phases (Table 492 4).

493

Figure 7 494

The temperature decrease during the early YD favoured the persistence of glaciers inside the highest hollows 495 in the Ruda Valley, depositing moraines at ~2400-2500 m, which were later abandoned by the ice at 12.6 ± 496 1.3 ka (Fig 7). The existence of glaciers at high altitudes of the Saboredo Cirque is also confirmed by polished 497 surfaces at elevations of 2200-2300 m that were dated at 12.7 ± 0.8 ka. In the neighbouring Bacivèr Cirque, 498 moraines at 2400 m were abandoned by the ice at 12.8 ka (Oliva et al., 2021a). The elevation difference 499 between the ELA during the YD (2571 m) in the Ruda Valley and present-day values indicates that summer 500 temperatures were ~3.4 ºC lower during that cold phase (Table 4), which is very similar to the ~3.0 ºC inferred 501 for the Bacivèr Cirque (Oliva et al., 2021a). Indeed, topoclimatic factors must have played a key role 502 maintaining small glaciers during the YD at the foot of steep cirque walls, where shading and snow 503 accumulation allowed the longer persistence of reduced ice masses (Boston and Lukas, 2019). In the Pyrenees, 504 YD glacial advances have been detected in cirques of the Central (Delmas, 2009;Delmas et al., 2008;García-505 Ruiz et al., 2016a) and Eastern sectors (Crest et al., 2017;Jomelli et al., 2020;Pallàs et al., 2010;Reixach et 506 al., 2021), where 2-3.5 km-long glaciers formed at 12.9-11.2 ka. Based on chronostratigraphic data from a 507 variety of sources, evidence of YD glacial advances was also recorded in the highest cirques of other Iberian 508 massifs, although no direct ages are yet available (García- Ruiz et al., 2016a). This is true for the NW ranges 509 (Cowton et al., 2009) and the Cantabrian Mountains (Pellitero et al., 2019;Serrano et al., 2015), where 510 moraines above 1800 and 2000 m, respectively, were associated with this cold phase. In the Sierra Nevada,

511

YD moraines dated at 13-12 ka are located in the highest cirques facing east from 2350 (north face) to 2800 m 512 (south face) (Oliva et al., 2014;Palacios et al., 2016). These glacial advances were driven by a temperature 513 cooling quantified at 2-4 ºC in western Europe that was triggered by a slowdown of the overturning circulation 514 in the North Atlantic (Renssen et al., 2018). In the Pyrenees, the YD was the coldest period of the last 515 deglaciation recorded by a decrease of 18 O in the Ostolo and Seso speleothems that might reflect the 516 temperature decrease of ~5ºC (Bartolomé et al., 2015;Bernal-Wormull et al., 2021). Colder conditions also 517 favoured glacial growth in other mountainous regions of Europe, with widespread advances across the Alps 518 between 13.5 and 12 ka (Ivy-Ochs, 2015;Ivy-Ochs et al., 2009), in the Tatra Mountains where moraines 519 stabilized at 11.9 ± 0.5 ka (Engel et al., 2017) and in Turkey where they formed at 12.6 ± 2.3 ka (Ciner and 520 Sarikaya, 2015).

521

The final deglaciation of the cirques favoured paraglacial dynamics and the formation of the numerous rock 522 glaciers existing today in the Saboredo and Sendrosa cirques, as well as in the lateral cirques distributed along 523 the Ruda Valley (Fig. 3). These features are currently inactive under the present-day climate regime (Fernandes 524 et al., 2018), with abundant lichen coverage and forest colonization, particularly on the lowest ridges. In the 525

Pyrenees, as well as in other Iberian ranges, cirque wall readjustment following deglaciation triggered different 526 landforms, such as debris-covered glaciers or rock glaciers, that remained active until the early Holocene 527 (Oliva et al., 2019). This is the case of the Bacivèr cirque, where the debris-covered glacier that formed inside 528 the highest YD moraines stabilized at 7.2 ± 0.5 ka (Oliva et al., 2021a), or of several other massifs in the 529

Central and Eastern Pyrenees where rock glaciers were active until well into the Holocene [START_REF] Andrés | Timing of deglaciation and rock glacier origin in the southeastern Pyrenees: a 594 review and new data[END_REF]530 García-Ruiz et al., 2016a). A similar pattern was also detected in other Iberian massifs, such as in the Iberian 531

Range [START_REF] Aumaître | Chronological and geomorphological investigation of fossil debris-602 covered glaciers in relation to deglaciation processes: A case study in the Sierra de La Demanda, 603 northern Spain[END_REF] or Sierra Nevada (Palacios et al., 2016), where these paraglacial landforms 532 formed immediately after the deglaciation of the cirques, and remained active until the complete melting of 533 their relict ice during the Holocene Thermal Maximum. 534

Conclusions

535

The Central Pyrenees were extensively glaciated during Quaternary glacial phases and were only slightly 536 transformed by periglacial, nival and slope processes during subsequent interglacial periods such as the 537

Holocene. As a result, the landscape in the Central Pyrenees includes a wide range of geomorphological 538 features that constitutes a prominent natural heritage that contributed to the designation of most protected areas 539 in this mountain range. This is the case of the Ruda Valley, located on the periphery of the National Park of 540

Aigüestortes and Sant Maurici Lake. The valley is one of the most scenic glacial valleys of the Pyrenees, and 541 includes a wide variety of glacial and periglacial landforms, particularly above 1800 m. However, in contrast 542 to the southern slope of the Central Pyrenees where the timing of the glacial evolution has been widely studied, 543 on the northern side there are still temporal and spatial gaps in the chronology of glacial activity. In this work, 544 we determined the magnitude and timing of T-1 glacial oscillations in the Upper Garonne Basin by producing 545 12 new CRE ages of glacial landforms (moraine boulders and glacially polished surfaces) that provide evidence 546 of the spatio-temporal pattern of glacial oscillations from the B-A to the YD.

547

The Upper Garonne Basin hosted the largest glacial system of the Pyrenees during the last glacial cycle thanks 548 to its Atlantic-influenced climate and the high altitude of the peaks in the central part of this range. The long-549 term post-LGM warming favoured glacial recession, particularly during the early B-A, and by 14.4 ± 1.2 ka 550 the Ruda Paleoglacier was already disconnected from the main Garonne Paleoglacier. At that time, the glacier 551 front was located at 1860-1900 m, and some areas of the Saboredo Cirque were already ice-free as a result of 552 ice thinning. The B-A was also a period of great climatic variability, and the small glaciers in the Ruda Valley 553

showed rapid response times with formation of moraines during the coldest B-A phases in north Iberia. Phases 554 of glacier advance or stillstands promoted the development of moraines at 13.5 ± 0.9 ka, as well as during the 555 transition between the late B-A and the onset of the YD at 13.0 ± 0.8 ka. The YD probably favoured the 556 persistence of small glaciers during an early stage of this cold phase; there is evidence of glacial retreat at the 557 uppermost part of the cirque floor at 12.7 ± 0.7 ka (~2300-2350 m) and moraine abandonment of the highest

120-

  To date the glacial landforms of the Upper Garonne Basin. 121 -To discuss the T-1 deglaciation in the Upper Garonne Basin and its geomorphological implications in 122 the context of the known chronologies from the Pyrenees southern slope and from other Iberian and 123 European mountain ranges. 124 -To examine the paleoclimatic and paleoenvironmental significance of the glacial response of the 125 northern slope of the Pyrenees and compare it with the paleoclimatic evolution inferred from other 126 records in the North Atlantic region.

  on the chronology of the last deglaciation in the Pyrenees has experienced substantial 406 advances over the last years in parallel to progress on glacial chronological constraints in other European and 407

Figure 1 . 888 Figure 2 . 902 Figure 6 .Figure 7 .

 1888290267 Geographical setting of the Ruda Valley. (A) Location in Europe. (B) Location in the Pyrenees with 886 the distribution of the studies mentioned in the discussion (green stars). (C) Limits and topography of the valley 887 with the distribution and view direction of the photos in Figure 2. Examples of the main geomorphological landforms of glacial origin in the Ruda Valley, together 889 with the location and typology of the dated samples: (A) terminal moraine damming a paleolake at 2080 m; 890 (B) N-S perspective of the Saboredo Cirque floor showing several lakes distributed in staggered overdeepened 891 basins; (C) polished bedrock surface including striations that reveal the former glacial flow; and (D) roches 892 moutonnées located in the forested valley bottom of the valley at 1860 m. 893 Figure 3. (A) Geomorphological map of the study area, with the location of the CRE samples; and (B) 894 enlargement of the highest parts of the cirque including the location of the 10 Be ages. 895 Figure 4. Examples of sampled moraines and polished bedrock surfaces, together with CRE ages (ka). 896 Figure 5. Reconstruction of the glacial extent during the different time stages: (A) Extensive glaciers occupying 897 the entire Upper Garonne Basin during the lLGM based on the reconstruction of Fernandes et al. (2017); (B) 898 Glacial extent during the mid B-A at ~13.5 ka in the Ruda Valley; (C) Glacier flowing from the Saboredo 899 Cirque during the transition between the late B-A and the early YD at ~13 ka; and (D) Small glaciers confined 900 within the highest hollows next to the Sendrosa (D) and Saboredo (E) peaks during the early YD at ~12.7/12.6 901 ka. Schematic evolution of the glacial phases recorded in the Ruda Valley during the last glacial cycle: 903 (A) extensive and thick glaciers inundated the valley and only few nunataks protruded the ice during the lLGM; 904 (B) warm temperatures in the early B-A favoured glacial recession, and subsequently the lowest areas of the 905 valley as well as some rock outcrops of the cirque became exposed; (C) glacial advances/stillstands formed 906 moraines during the B-A and at the transition to the YD; (D) small glaciers persisted during the early YD in 907 the highest north-facing cirques; and (E) subsequent temperature increases favoured complete deglaciation and 908 triggered paraglacial dynamics, which promoted the formation of rock glaciers in the recently deglaciated ice-909 free cirques. 910 Normalized probability distribution functions (PDF) of surface exposure ages vs. temperature 911 evolution from the OD to the early Holocene based on the δ 18 O record from the NGRIP ice core from 912 Greenland and the δ 18 O record from the speleothem of the Ostolo cave in north Iberia (time periods are defined 913 after Rasmussen et al. (2014)). The plots of the units result from the sum of the individual PDF of the samples 914 belonging to them; cold (warm) phases are represented by the blue (red) bands.

119

  Our study, therefore, has the following objectives: 120 -To date the glacial landforms of the Upper Garonne Basin. 121 -To discuss the T-1 deglaciation in the Upper Garonne Basin and its geomorphological implications in 122 the context of the known chronologies from the Pyrenees southern slope and from other Iberian and 123 European mountain ranges. 124 -To examine the paleoclimatic and paleoenvironmental significance of the glacial response of the 125 northern slope of the Pyrenees and compare it with the paleoclimatic evolution inferred from other 126 records in the North Atlantic region.

  on the chronology of the last deglaciation in the Pyrenees has experienced substantial 406 advances over the last years in parallel to progress on glacial chronological constraints in other European and 407

Figure 1 . 888 Figure 2 . 902 Figure 6 .Figure 7 .

 1888290267 Geographical setting of the Ruda Valley. (A) Location in Europe. (B) Location in the Pyrenees with 886 the distribution of the studies mentioned in the discussion (green stars). (C) Limits and topography of the valley 887 with the distribution and view direction of the photos in Figure 2. Examples of the main geomorphological landforms of glacial origin in the Ruda Valley, together 889 with the location and typology of the dated samples: (A) terminal moraine damming a paleolake at 2080 m; 890 (B) N-S perspective of the Saboredo Cirque floor showing several lakes distributed in staggered overdeepened 891 basins; (C) polished bedrock surface including striations that reveal the former glacial flow; and (D) roches 892 moutonnées located in the forested valley bottom of the valley at 1860 m. 893 Figure 3. (A) Geomorphological map of the study area, with the location of the CRE samples; and (B) 894 enlargement of the highest parts of the cirque including the location of the 10 Be ages. 895 Figure 4. Examples of sampled moraines and polished bedrock surfaces, together with CRE ages (ka). 896 Figure 5. Reconstruction of the glacial extent during the different time stages: (A) Extensive glaciers occupying 897 the entire Upper Garonne Basin during the lLGM based on the reconstruction of Fernandes et al. (2017); (B) 898 Glacial extent during the mid B-A at ~13.5 ka in the Ruda Valley; (C) Glacier flowing from the Saboredo 899 Cirque during the transition between the late B-A and the early YD at ~13 ka; and (D) Small glaciers confined 900 within the highest hollows next to the Sendrosa (D) and Saboredo (E) peaks during the early YD at ~12.7/12.6 901 ka. Schematic evolution of the glacial phases recorded in the Ruda Valley during the last glacial cycle: 903 (A) extensive and thick glaciers inundated the valley and only few nunataks protruded the ice during the lLGM; 904 (B) warm temperatures in the early B-A favoured glacial recession, and subsequently the lowest areas of the 905 valley as well as some rock outcrops of the cirque became exposed; (C) glacial advances/stillstands formed 906 moraines during the B-A and at the transition to the YD; (D) small glaciers persisted during the early YD in 907 the highest north-facing cirques; and (E) subsequent temperature increases favoured complete deglaciation and 908 triggered paraglacial dynamics, which promoted the formation of rock glaciers in the recently deglaciated ice-909 free cirques. 910 Normalized probability distribution functions (PDF) of surface exposure ages vs. temperature 911 evolution from the OD to the early Holocene based on the δ 18 O record from the NGRIP ice core from 912 Greenland and the δ 18 O record from the speleothem of the Ostolo cave in north Iberia (time periods are defined 913 after Rasmussen et al. (2014)). The plots of the units result from the sum of the individual PDF of the samples 914 belonging to them; cold (warm) phases are represented by the blue (red) bands.

Figure 1 .

 1 Figure 1. Geographical setting of the Ruda Valley. (A) Location in Europe. (B) Location in the Pyrenees with the distribution of the studies mentioned in the discussion (green stars). (C) Limits and topography of the valley with the distribution and view direction of the photos in Figure 2. 2476x1757mm (39 x 39 DPI)

  Glacial and periglacial landforms are preserved on different lithologies across the 153 valley, composed of Carboniferous granodiorites and granites in the highest sections, with Cambrian-Devonian 154 metamorphic (hornfels, marbles, slates, quartzites, schists) and sedimentary rocks (limestones, sandstones, 155 conglomerates and lutites) in the lowest parts of the valley (ICGC, 2017; Quesada and Oliveira, 2019). 156 However, despite abundant geomorphic evidence revealing that the Garonne Valley presented the longest 157 glacier in the Pyrenees during the lLGM (ca. 80 km; Fernandes et al., 2017), the chronological framework of 158 glacial activity in the entire Upper Garonne Basin is still unknown. The headwaters of the Ruda Valley include 159 a rich sequence of moraine systems and polished bedrock surfaces -as well as relict rock glaciers and protalus 160 lobes (Fernandes et al., 2018) -that can be used to improve the temporal resolution for the chronosequence of 161 the final deglaciation of the upper Garonne Valley.

162 163 3. Methodology 164 We reconstructed the last deglaciation chronology of the Ruda Valley in the Upper Garonne Basin by 165 integrating geomorphological and chronological techniques. Field work was conducted during the summer of 166 2019, when snow cover was limited to the concave and sheltered areas at the highest elevations (>2500 m), 167 which enabled better identification of the glacial landforms, as well as sampling for CRE dating.

168 3.1 Geomorphological mapping 169

Table 1

 1

	186	Figure 2

Laboratory procedures and exposure age calculation 188

  Samples were crushed and sieved (>200 g) at the Physical Geography Laboratory of the Universidad 189 Complutense de Madrid (Spain) to the 125-500 µm fraction, which was subsequently processed at the 190 'Laboratoire National des Nucléides Cosmogéniques' (LN 2 C) of the 'Centre Européen de Recherche et 191 d'Enseignement des Géosciences de l'Environnement' (CEREGE; Aix-en-Provence, France). In accordance 192 with the quartz-rich lithology of the samples, they were treated for the extraction of the in-situ-produced 193 cosmogenic nuclide 10 Be.

	194	First, magnetic minerals were discarded through magnetic separation conducted in a "Frantz LB-1" magnetic
	195	separator. Once the non-magnetic fraction was isolated, it underwent several rounds of chemical attacks with
	196	a 1:2 mixture of concentrated hydrochloric (HCl) and hexafluorosilicic (H 2 SiF 6 , 2/3) acids to dissolve and
	197	discard the non-quartz minerals. The remaining minerals (i.e. non-dissolved feldspar minerals) were
	198

  resulted in ages that were ~9% older (~1% erosion; ~8% snow cover). However, to enable 235 comparisons with other areas, we use the non-corrected ages throughout the text. Moreover, both the thickness 236 and duration of snow cover must have varied significantly since the cirque's deglaciation, given the alternation 237

								1)
	232	available	since	1997-1998	(Servei	Meteorològic	de	Catalunya;
	233

http://www.igc.cat/web/ca/allaus_gruix_neu_v2.php?e=bonaigua&t=totes). Erosion and snow corrections, as a whole,

2 Geochronological data 298

		(ii) The highest peaks of the Ruda Valley surrounding the Saboredo Cirque range from 2700-2800 m
	280	(Sendrosa and Saboredo peaks) to 2525 m in the mountain divide that may have functioned as a glacial
	281	difluence pass between the Saboredo (Garonne basin) and the Ratera cirques (Noguera Pallaresa basin)
	282	(Fig. 3). Saboredo is a complex cirque (Barr and Spagnolo, 2015) that includes several minor hollows
	283	carved out on the marginal flanks of the peaks. They are separated by truncated spurs with polished surfaces,
	284	where one sample (ARAN-31) was collected to examine the timing of glacial thinning at the cirque level.
	285	The highest moraines in the Ruda Valley are located at ~2400-2500 m, at the foot of the highest peaks. A
	286	prominent cirque moraine at an altitude of 2470 m is located on the eastern slope of Sendrosa Peak; here,
	287	we collected two samples (ARAN-33 and ARAN-32). Above the upper moraines, fed by abundant debris
	288	supply from the cirque rock walls, rock glaciers and protalus lobes are found. The cirque floors extend
	289	between 2200 and 2400 m and preserve features indicative of intense glacial erosion (e.g. overdeepened
	290	basins or glacial thresholds) and glacial deposits (e.g. moraines or erratic boulders). In the Saboredo Cirque,
	291	glacial abrasion and periglacial processes shaped a large amphitheatre with several overdeepened basins.
	292	They are surrounded by polished surfaces and occupied by peatlands and lakes. Three samples were
	293	collected from polished surfaces from 2360 to 2270 m (Fig. 3) (ARAN-36; ARAN-35; ARAN-34; Fig. 2;
	294	Table 1). Glacial retreat left scattered boulders and till across the area.
	295	Figure 3
	296	Figure 4
	297	4.
		Ruda Valley preserves glacial and periglacial landforms indicative of glacial evolution during the last
	258	deglaciation in the Upper Garonne Basin. The chronological framework of this evolution is unveiled by a
	259	dataset of 12 10 Be CRE ages obtained from glacial landforms.
	260	4.1 Geomorphological evidence
	261	The wide range of glacial features from the highest sectors to the Ruda Valley floor shows evidence that the
	262	entire valley was heavily glaciated during the lLGM and recorded several advances or stillstands within the
	263	long-term deglaciation trend. The Ruda Valley is divided into two main geomorphological units (Fig. 3):
	264	(i) The valley floor and mountain slopes follow the typical U-shaped morphology of a glacial valley. The
	265	steep slopes, descending from ~2200 to 1800 m, connect the lateral cirques with the valley floor. The
	266	hillslopes of the valley start at the upper edge of hanging glacial cirques that define the glacial trimline,
	267	marking the uppermost extent of recent glaciations. Hillslopes are concave to rectilinear 300-400 m-long
	268	sections steeper than 30º that are mostly colonized by forest (i.e. Pinus nigra). The upper part of the valley
	269	is defined by two staggered, overdeepened basins filled with peatlands at average elevations of 2150 and
	270	2070 m. At their margins, two moraine systems occur at 2190 and 2080 m. Two samples were collected
	271	from boulders from the external frontal moraine (ARAN-38 and ARAN-37) and the other two were
	272	obtained from the internal moraine unit: one from a boulder on the moraine ridge (ARAN-27), and the other
	273	from a polished surface with glacial striae (ARAN-28) and moraine boulders (Fig. 4). The lowest part of
	274	the valley floor, with a relatively flat slope characteristic of U-shaped glacial valleys and descending from
	275	~1800 to 1500 m, is colonized by alpine meadows and coniferous forests. Isolated roches moutonnées are
	276	found in the highest section of the valley bottom where granites prevail (> ~1800 m), while they are absent
	277	in lower areas, where slates and marbles are the dominant lithology. Two samples were collected from
	278	polished surfaces at 1900 and 1860 m (ARAN-39 and ARAN-40).

3 Reconstructed paleoglaciers and ELAs

	310	Two samples from polished surfaces of the Saboredo Cirque above 2300 m yielded ages of 12.8 ± 0.8 (ARAN-
	311	35) and 12.7 ± 0.8 ka (ARAN-36) (12.7 ± 0.8 ka). The highest moraines are located inside the highest cirque
	312	hollows, such as the Sendrosa hollow. Here, two moraine boulders yielded ages of 11.8 ± 0.7 (ARAN-33) and
	313	13.3 ± 1.0 ka (ARAN-32) (12.6 ± 1.3 ka).
	314	4.

  , Malniu[START_REF] Andrés | Timing of deglaciation and rock glacier origin in the southeastern Pyrenees: a 594 review and new data[END_REF] Palacios et al., 2015b; Pallàs 426 et al., 2010), Têt(Delmas et al., 2008; Tomkins et al., 2018), Arànser[START_REF] Andrés | Timing of deglaciation and rock glacier origin in the southeastern Pyrenees: a 594 review and new data[END_REF] Palacios et al., 427 2015b), and Bassiès(Crest et al., 2017). During that time, glaciers flowed downslope from the headwaters of 428 the highest valleys, with valley glaciers often exceeding 10-15 km(Oliva et al., 2019). Consequently, it is 429 likely that, by that time, most of the Ruda Valley was still glaciated, with the ice tongue probably merging 430 with the main Garonne Paleoglacier that was also fed by neighbouring glacier tributaries, such as those 431 descending from the Beret area(Oliva et al., 2021a). Some of the lateral moraine remnants distributed on the 432 left side of the Garonne Valley, ca. 20 km from the highest parts of the Saboredo Cirque, formerly attributed 433 to the OD(Fernandes et al., 2021; Fernandes et al., 2017) may be also correlated with some of the highest 434 moraines and erratic boulders distributed in the Ruda Valley (Fig.3). Other Iberian mountain ranges include 435 geomorphic evidence of glacial advances during the OD (Palacios et al.

, 2017a), such as in the Monasterio 436

Valley in the Cantabrian Mountains

(Serrano et al., 2016)

, Peña Negra Valley in the Iberian Range (García-437

[START_REF] Ruiz | Peñalara massif in the Central Range (Palacios et al[END_REF]

, Peñalara massif in the Central Range

(Palacios et al., 2012a)

, San Juan Valley in the Sierra 438

Nevada

(Palacios et al., 2016)

, En Garcia Valley in the Eastern Pyrenees

(Reixach et al., 2021)

, and Gállego 439

  , Cantabrian Mountains 451 (Rodríguez-Rodríguez et al., 2017) and the Iberian Range (García-Ruiz et al., 2020), where ice masses shrank 452 significantly by ~15-14 ka. Glacial retreat in the Iberian Mountains was triggered by the sudden temperature 453 increase of 3-5 ºC at the onset of the B-A in western Europe, as inferred from climate models (Renssen and 454 Isarin, 2001). In parallel, significant shifts occurred in the Atlantic Meridional Overturning Circulation (Obase 455 and Abe-Ouchi, 2019). In north Iberia, a maximum increase of the mean annual temperature of 7.5ºC at the 456 onset of GI-1 has been estimated in northern Iberia, based on the Ostolo cave speleothem isotopic record 457 (Bernal

  7 ka (~2300-2350 m) and moraine abandonment of the highest 558 cirque moraines (2470 m) at 12.6 ± 1.3 ka. Subsequently, glacial disappearance must have promoted 559 paraglacial dynamics and the formation of the (currently relict) rock glaciers existing within these moraines. 560 Whereas the impact of the B-A deglaciation on the current landscape of the Pyrenean cirques was already 561 known, the sequence of glacial advances and retreats during the B-A and YD presented in this study provides 562 one of the most accurate deglacial histories in southern Europe. This opens new perspectives on the 563 paleoclimatic evolution in the Iberian Peninsula, in particular, and the Mediterranean region in general. The 564 period spanning from the B-A to the YD has been shown to be a major driver of transformation of the mountain 565 landscapes of the Pyrenees and a key phase to better frame Holocene environmental dynamics in the highest 566 lands of this mountain range. Future studies with new chronological data from other Pyrenean valleys and 567 Iberian mountains will help to better understand T-1 glacier advances and retreats as well as the climate 568 mechanisms behind millennial-scale glacial oscillations in southern Europe. This work was funded by the Research Group ANTALP (Antarctic, Arctic, Alpine Environments; 2017-SGR-571 1102), the Government of Catalonia and the Centro de Estudos Geográficos/IGOT -University of Lisbon (FCT 572 I.P. UIDB/00295/2020 and UIDP/00295/2020). The research topics complement those of the project 573 PALEOGREEN (CTM2017-87976-P) funded by the Spanish Ministry of Economy and Competitiveness and 574 the project NUNANTAR funded by the Fundação para a Ciência e Tecnologia of Portugal (02/SAICT/2017 -575 32002). Marcelo Fernandes holds a PhD fellowship of the Fundação para a Ciência e Tecnologia of Portugal 576 (FCT -UIDB/00295/2020); Marc Oliva is supported by the Ramón y Cajal Program (RYC-2015-17597) and 577 José M. Fernández-Fernández is supported by a postdoctoral grant within the NUNANTAR project. 10 Be 578 measurements were performed at the ASTER AMS national facility (CEREGE, Aix-en-Provence), which is 579 supported by the INSU/CNRS and the ANR through the "Projets thématiques d'excellence" program for the 580 "Equipements d'excellence" ASTER-CEREGE action and IRD. This research is also framed within the College 581 on Polar and Extreme Environments (Polar2E) of the University of Lisbon. We also thank the National Park 582of Aigüestortes and Sant Maurici Lake for providing field access to the study sites and the reviewers for their 583 constructive comments that helped to improve the quality of an earlier version of the manuscript.

	569	Acknowledgements
	570
	584
	585	7. References
	586
	587	Allard, J.L., Hughes, P., Woodward, J.C., 2021. Heinrich Stadial aridity forced Mediterranean-wide glacier
	588	retreat in the last cold stage. Nat. Geosci. https://doi.org/10.1038/s41561-021-00703-6
	589	André, M.F., 2002. Rates of Postglacial rock weathering on glacially scoured outcrops (Abisko-Riksgränsen
	590	area, 68°N). Geogr. Ann. Ser. A Phys. Geogr. 84, 139-150. https://doi.org/10.1111/j.0435-
	591	3676.2002.00168.x
	592	Andrés, N., Gómez-Ortiz, A., Fernández-Fernández, J.M., Tanarro García, L.M., Salvador-Franch, F., Oliva,
	593	M., Palacios, D., 2018. Timing of deglaciation and rock glacier origin in the southeastern Pyrenees: a
	594	review and new data. Boreas 47, 1050-1071. https://doi.org/10.1111/bor.12324
	595	Arnold, M., Merchel, S., Bourlès, D.L., Braucher, R., Benedetti, L., Finkel, R.C., Aumaître, G., Gottdang,
	596	A., Klein, M., 2010. The French accelerator mass spectrometry facility ASTER: Improved performance
	597	and developments. Nucl. Instruments Methods Phys. Res. Sect. B Beam Interact. with Mater. Atoms
	598	268, 1954-1959. https://doi.org/10.1016/j.nimb.2010.02.107
	599	Aumaître, G., Bourlès, D., Keddadouche, K., Schimmelpfennig, I., Léanni, L., Fernández-Fernández, J.M.,
	600	Palacios, D.

, Andrés, N., Úbeda, J., García-Ruiz, J.M., Gómez-Villar, A., Santos-González, J., Álvarez-601 Martínez, J., Arnáez, J., 2017. Chronological and geomorphological investigation of fossil debris-602 covered glaciers in relation to deglaciation processes: A case study in the Sierra de La Demanda, 603 northern Spain. Quat. Sci. Rev. 170, 232-249. https://doi.org/10.1016/j.quascirev.2017.06.034 604 Ballantyne, C.K., 2008. After the ice: Holocene geomorphic activity in the Scottish Highlands, Scottish

Table 1 . Geographic location of samples, topographic shielding factor and sample thickness.

 1

	Sample name	Landform	Latitude (DD)	Longitude (DD)	Elevation (m a.s.l.) a Topographic shielding	Thickness (cm)	Lithology	Size of the boulders
						factor			largest axis (m)
	Valley polished surfaces
	ARAN-40	Polished surface	42.6416	0.9713	1862	0.9630	3.7	Granodiorite
	ARAN-39	Polished surface	42.6393	0.9701	1900	0.9640	5.0	Granodiorite
	Valley external moraines
	ARAN-38	Moraine boulder	42.6296	0.9624	2079	0.9759	4.2	Granite	3.5
	ARAN-37	Moraine boulder	42.6299	0.9620	2075	0.9752	3.0	Granite	3
	Valley internal moraines
	ARAN-28	Polished surface	42.6263	0.9593	2188	1.0000	3.5	Granite
	ARAN-27	Moraine boulder	42.6257	0.9599	2184	0.9852	5.0	Granite	4
	Cirque polished surfaces
	ARAN-36	Polished surface	42.6111	0.9681	2360	0.9806	3.5	Granite
	ARAN-35	Polished surface	42.6152	0.9630	2316	0.9875	3.2	Granite
	ARAN-34	Polished surface	42.6224	0.9596	2266	0.9901	5.0	Granite
	ARAN-31	Polished surface	42.6274	0.9536	2378	1.0000	3.0	Granite
	Cirque moraines
	ARAN-33	Moraine boulder	42.6275	0.9487	2470	0.9841	3.0	Granite	5
	ARAN-32	Moraine boulder	42.6277	0.9480	2467	0.9789	4.0	Granite	6

a Elevations are derived from the 5 m Digital Elevation Model of the Spanish "Instituto Geográfico Nacional" and are subjected to a vertical accuracy of ±5 m.

Table 2 . AMS analytical data and calculated exposure ages. 10 Be/ 9 Be ratios were inferred from measurements at the ASTER AMS facility. Chemistry blank correction was done by subtracting the number of 10 Be atoms in the blank from that in the sample. 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. 10 Be samples analytical AMS data

 2

	Sample name	Quartz weight (g) Mass of carrier	ASTER AMS	10 Be/ 9 Be (10 -14 )	Blank correction	[ 10 Be] (10 4 atoms g -1 )	Age (ka)	Mean age (ka)
			( 9 Be mg)	cathode number		(%)
	Valley polished surfaces
	ARAN-40 ARAN-39	20.44 21.32	0.46 0.45	CHBL CHBK	16.34 ± 0.56 15.92 ± 0.49	2.11 2.18	23.85 ± 0.83 22.15 ± 0.70	15.0 ± 0.9 (0.5) 13.8 ± 0.8 (0.4)	14.4 ± 1.2 (0.8)
	Valley external moraines
	ARAN-38 ARAN-37	21.11 21.38	0.46 0.46	CHBJ CHBI	17.56 ± 0.54 18.08 ± 0.87	1.95 1.09	24.99 ± 0.79 25.74 ± 1.25	13.3 ± 0.8 (0.4) 13.7 ± 1.0 (0.6)	13.5 ± 0.9 (0.2)
	Valley internal moraines
	ARAN-28 ARAN-27	20.97 20.89	0.46 0.45	CHBB CHBA	18.79 ± 0.62 18.64 ± 0.58	1.05 1.85	27.13 ± 0.91 26.61 ± 0.84	13.0 ± 0.8 (0.4) 13.1 ± 0.8 (0.4)	13.0 ± 0.8 (0.1)
	Cirque polished surfaces
	ARAN-36 ARAN-35	21.51 21.43	0.46 0.47	CHBH CHBG	21.24 ± 0.65 20.18 ± 0.62	1.62 1.65	29.53 ± 0.92 29.05 ± 0.91	12.7 ± 0.8 (0.4) 12.8 ± 0.8 (0.4)	12.7 ± 0.8 (0.0)
	ARAN-34	20.85	0.46	CHBF	21.13 ± 0.66	1.63	30.35 ± 0.96	14.0 ± 0.8 (0.4)
	ARAN-31	21.40	0.44	CHBC	26.54 ± 0.72	1.35	35.74 ± 0.98	14.8 ± 0.9 (0.4)
	Cirque moraines
	ARAN-33 ARAN-32	20.26 21.63	0.45 0.45	CHBE CHBD	20.30 ± 0.68 24.19 ± 1.26	1.70 0.85	29.92 ± 1.02 33.15 ± 1.75	11.8 ± 0.7 (0.4) 13.3 ± 1.0 (0.7)	12.6 ± 1.3 (1.0)
	Chemistry blank details
	Blank name	Processed with	mass of carrier	ASTER AMS	10 Be/ 9 Be (10 -14 )	[ 10 Be] (10 4 atoms)
			( 9 Be mg)	cathode number
		ARAN-27, 31, 33,
	BK-1	34, 35, 36, 38, 39,	0.46	CHAT	0.34 ± 0.04	10.30 ± 1.24	-	-
		40
	BK-2	ARAN-32	0.46	CHBM	0.20 ± 0.03	6.13 ± 0.80	-	-
	BK-3	ARAN-28, 37	0.46	IGHI	0.20 ± 0.03	5.98 ± 0.93	-	-

Table 3 . Exposure ages according to erosion and snow cover corrections. Exposure ages (arithmetic mean, in ka) Sample name No correction Erosion correction

 3

	Snow correction	Erosion + snow
		correction

Table 4 . Reconstructed ELAs (m a.s.l.) for the middle, late B-A and YD moraines using the AAR and AABR methods.

 4

	Phase	Glacier	AAR (0.6±0.2)	AABR (1.9±0.81)	AABR (1.75±0.71)	Average ELA
	Valley external moraines (middle B-A; GI-1c)	Ruda	2481 ± 20	2449 -30/+20	2454 -30/+15	2461
	Valley internal moraines (transition B-A/YD; GI-1b)	Ruda	2515 -15/+10	2498 -25/+15	2503 -25/+20	2505
	Cirque moraines (early YD)	Sendrosa	2568 -10/+5	2571 -15/+5	2576 ± 10	2571

Table 5 . Glacial response in the Ruda Valley during the B-A Interstidal according to the reconstructed climate in Greenland and in Iberian Peninsula.

 5

  timing of glacial oscillations is still poorly constrained on the northern slope of the Pyrenees, such as in 106 the Upper Garonne Basin where this work focuses on. Recently, Reixach et al. (2021) have summarized the 107 glacial evolution from 19 to 12 ka in the northern Pyrenees, reporting a spatial and temporal pattern for T-1 108 glacial oscillations similar to that on the southern slopes of the range. Jomelli et al. (2020) confirmed glacial 109 growth during the YD in the Ariège Valley, where glaciers persisted in the highest cirques until the early 110

  Saboredo (2827 m), Sendrosa (2702 m) and Locampo (2658 m) peaks are the highest elevations of the 143 cirque, while its toe is located at ~2250 m.144The present-day mean annual air temperature (MAAT) at Port de la Bonaigua station (2266 m) is 3 ºC, with 145 the regional 0 ºC isotherm in the Central Pyrenees estimated to be at 2950 m (López-Moreno et al., 2019).146Annual precipitation reaches 1227 mm and falls mostly as snow (average depth of 1 m) from October to May, 147 with a snow duration of 7-8 months in the cirque floor. The combination of low temperatures and abundant 148 snowfall determines the scarce vegetation cover above the treeline (2200-2300 m), where barren rocky terrain 149 is widespread, together with small alpine meadows.150The landscape of the Ruda Valley resulted from the interaction of tectonics and vertical incision, both during 151Quaternary glaciations and interglacial periods, the latter dominated by alluvial and mass-wasting processes 152(Fernandes et al., 2017). Glacial and periglacial landforms are preserved on different lithologies across the 153 valley, composed of Carboniferous granodiorites and granites in the highest sections, with Cambrian-Devonian 154 metamorphic (hornfels, marbles, slates, quartzites, schists) and sedimentary rocks (limestones, sandstones, 155 conglomerates and lutites) in the lowest parts of the valley (ICGC, 2017; Quesada and Oliveira, 2019). 156 However, despite abundant geomorphic evidence revealing that the Garonne Valley presented the longest 157 glacier in the Pyrenees during the lLGM (ca. 80 km; Fernandes et al., 2017), the chronological framework of 158 glacial activity in the entire Upper Garonne Basin is still unknown. The headwaters of the Ruda Valley include 159 a rich sequence of moraine systems and polished bedrock surfaces -as well as relict rock glaciers and protalus 160 lobes (Fernandes et al., 2018) -that can be used to improve the temporal resolution for the chronosequence of 161 the final deglaciation of the upper Garonne Valley. We reconstructed the last deglaciation chronology of the Ruda Valley in the Upper Garonne Basin by 165 integrating geomorphological and chronological techniques. Field work was conducted during the summer of 166

	162
	163	3. Methodology
	164

U-shaped Ruda Valley, drained by the Ruda River that runs 12 km northwards, constitutes the headwaters 139 of the Upper Garonne Basin. Its catchment encompasses 34 km 2 , with an elevation range of ~1300 m between 140 Saboredo Peak (2827 m) and the village of Baqueira (~1500 m). The Saboredo Cirque is located in the 141 headwaters of the Ruda Valley; the cirque is N-exposed, extending over 5.7 km 2 , with its floor at ~2250-2400 142 m.

3 Laboratory procedures and exposure age calculation 188

  Samples were crushed and sieved (>200 g) at the Physical Geography Laboratory of the Universidad 189 Complutense de Madrid (Spain) to the 125-500 µm fraction, which was subsequently processed at the 190 'Laboratoire National des Nucléides Cosmogéniques' (LN 2 C) of the 'Centre Européen de Recherche et 191 d'Enseignement des Géosciences de l'Environnement' (CEREGE; Aix-en-Provence, France). In accordance 192

  The highest peaks of the Ruda Valley surrounding the Saboredo Cirque range from 2700-2800 m 280

								1)
	232	available	since	1997-1998	(Servei	Meteorològic	de	Catalunya;
	233

http://www.igc.cat/web/ca/allaus_gruix_neu_v2.php?e=bonaigua&t=totes). Erosion and snow corrections, as

(ii)

  . Until the onset of 350 deglaciation, 300-400 m thick ice sat in the Ruda Valley, as shown by the highest lateral moraines between 351 the Ruda and Aiguamòg valleys. This ice thickness must have favoured the occurrence of warm-based glaciers, 352 thus implying basal sliding evidenced by striations on the bedrock surface. Hence, no nuclide inheritance from 353 previous exposure periods is expected at the sampled surfaces. 354 Samples from polished surfaces from the bottom of the Ruda Valley at 1860-1900 m indicate that by 14.4 ± 355 1.2 ka, i.e. during the B-A, the glacier was retreating towards the glacial cirque, and was already disconnected 356 from the main Ruda glacier. Consequently, one third of the Ruda Valley was already deglaciated at that time, 357

Table 2 . AMS analytical data and calculated exposure ages. 10 Be/ 9 Be ratios were inferred from measurements at the ASTER AMS facility. Chemistry blank correction was done by subtracting the number of 10 Be atoms in the blank from that in the sample. 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. 10 Be samples analytical AMS data

 2

	Sample name	Quartz weight (g) Mass of carrier	ASTER AMS	10 Be/ 9 Be (10 -14 )	Blank correction	[ 10 Be] (10 4 atoms g -1 )	Age (ka)	Mean age (ka)
			( 9 Be mg)	cathode number		(%)
	Valley polished surfaces
	ARAN-40 ARAN-39	20.44 21.32	0.46 0.45	CHBL CHBK	16.34 ± 0.56 15.92 ± 0.49	2.11 2.18	23.85 ± 0.83 22.15 ± 0.70	15.0 ± 0.9 (0.5) 13.8 ± 0.8 (0.4)	14.4 ± 1.2 (0.8)
	Valley external moraines
	ARAN-38 ARAN-37	21.11 21.38	0.46 0.46	CHBJ CHBI	17.56 ± 0.54 18.08 ± 0.87	1.95 1.09	24.99 ± 0.79 25.74 ± 1.25	13.3 ± 0.8 (0.4) 13.7 ± 1.0 (0.6)	13.5 ± 0.9 (0.2)
	Valley internal moraines
	ARAN-28 ARAN-27	20.97 20.89	0.46 0.45	CHBB CHBA	18.79 ± 0.62 18.64 ± 0.58	1.05 1.85	27.13 ± 0.91 26.61 ± 0.84	13.0 ± 0.8 (0.4) 13.1 ± 0.8 (0.4)	13.0 ± 0.8 (0.1)
	Cirque polished surfaces
	ARAN-36 ARAN-35	21.51 21.43	0.46 0.47	CHBH CHBG	21.24 ± 0.65 20.18 ± 0.62	1.62 1.65	29.53 ± 0.92 29.05 ± 0.91	12.7 ± 0.8 (0.4) 12.8 ± 0.8 (0.4)	12.7 ± 0.8 (0.0)
	ARAN-34	20.85	0.46	CHBF	21.13 ± 0.66	1.63	30.35 ± 0.96	14.0 ± 0.8 (0.4)
	ARAN-31	21.40	0.44	CHBC	26.54 ± 0.72	1.35	35.74 ± 0.98	14.8 ± 0.9 (0.4)
	Cirque moraines
	ARAN-33 ARAN-32	20.26 21.63	0.45 0.45	CHBE CHBD	20.30 ± 0.68 24.19 ± 1.26	1.70 0.85	29.92 ± 1.02 33.15 ± 1.75	11.8 ± 0.7 (0.4) 13.3 ± 1.0 (0.7)	12.6 ± 1.3 (1.0)
	Chemistry blank details
	Blank name	Processed with	mass of carrier	ASTER AMS	10 Be/ 9 Be (10 -14 )	[ 10 Be] (10 4 atoms)
			( 9 Be mg)	cathode number
		ARAN-27, 31, 33,
	BK-1	34, 35, 36, 38, 39,	0.46	CHAT	0.34 ± 0.04	10.30 ± 1.24	-	-
		40
	BK-2	ARAN-32	0.46	CHBM	0.20 ± 0.03	6.13 ± 0.80	-	-
	BK-3	ARAN-28, 37	0.46	IGHI	0.20 ± 0.03	5.98 ± 0.93	-	-

Table 3 . Exposure ages according to erosion and snow cover corrections. Exposure ages (arithmetic mean, in ka) Sample name No correction Erosion correction

 3

	Snow correction	Erosion + snow
		correction

http://mc.manuscriptcentral.com/jqs Journal of Quaternary Science

Abstract

43

The Upper Garonne Basin included the largest glacial system in the Pyrenees during the last glacial cycle.

44

Within the long-term glacial retreat during Termination-1 (T-1), glacier fluctuations left geomorphic evidence 45 in the area. However, the chronology of T-1 glacial oscillations on the northern slope of the Central Pyrenees 46 is still poorly constrained. Here, we introduce new geomorphological observations and a 12-sample dataset of 47 10 Be Cosmic-Ray Exposure (CRE) ages from the Ruda Valley. This U-shaped valley, surrounded by peaks 48 exceeding 2800 m a.s.l., includes a sequence of moraines and polished surfaces that enabled to reconstruct the 49 chronology of the last deglaciation. Following the maximum ice extent, warmer conditions prevailing at ~15-50 14 ka, during the Bølling-Allerød (B-A) Interstadial, favoured glacial retreat in the Ruda Valley. Within the 51 B-A, glaciers experienced two phases of advance/stillstand with moraine formation at 13.5 and 13.0 ka. (Clark et al., 2009), most glaciated mountain regions underwent a rapid deglaciation only interrupted by 64 brief periods of glacial growth. Glacial disappearance favoured paraglacial readjustment and intense 65 periglacial and slope processes, which reshaped the already deglaciated valleys and cirques [START_REF] Ballantyne | After the ice: Holocene geomorphic activity in the Scottish Highlands, Scottish Barr, I[END_REF]. 66

The current spatial distribution of glacial and periglacial features in the high elevations of most mid-latitude 67 ice-free mountain ranges -such as the Pyrenees, where this research focuses on -results from the climatic 68 shifts during Termination-1 (hereafter T-1), the period spanning from 20-19 to 11.7 ka (Denton et al., 2014).

69

At present, mid-latitude mountains such as those of the Mediterranean region are only sparsely glaciated, with 70 small glaciers located in sheltered cirques (Hughes, 2014). However, there are several other mountain ranges 71 located between 30 and 50º latitude North, such as the Alps or the Himalaya, which are still heavily glaciated 72 and record recent dramatic glacial recession trends (Huss et al., 2021). An accurate comprehension of the 73 climatic sensitivity of current glaciated mountain landscapes can help to understand the timing and magnitude 74 of glacial change during the rapid climate shifts of T-1. The alternation of millennial-scale cold/warm periods 75

during T-1 promoted a temporal pattern of glacial advances and retreats across the Northern Hemisphere 76 [START_REF] Allard | Heinrich Stadial aridity forced Mediterranean-wide glacier 588 retreat in the last cold stage[END_REF]Rea et al., 2020). Ice core records from the Greenland interior show the following 77 succession of climatic phases in the North Atlantic region during T-1: a cold climate regime during the Oldest 78 Dryas (OD; GS-2.1a: 17.5-14.6 ka), followed by the general warm Bølling-Allerød (B-A) Interstadial (GI-1: Gradual warming in the northern Iberia at 15.4-13.4 cal ka BP (Moreno et al., 2010). Warming conditions in the Pyrenees at 14.7-14.0 cal ka BP (González-Sampériz et al., 2006). Significant warming at the beginning (2-4ºC) followed by a gradual cooling at (~2ºC). Intense cold conditions at 13.7-13.6 ka Glacier advance or stillstand with moraine formation at 13.5 ka

List of tables

Colder conditions starting at 13.5 cal ka BP in northern Iberia (Moreno et al., 2010). Cold climate in the Pyrenees at 14-13.4 cal ka BP (González-Sampériz et al., 2006).

GI-1b (13.3-13.1) Short cold reversal (1-2ºC) -Short periods of colder and drier conditions (Moreno et al., 2010). Expansion of coldclimate forest species in the Pyrenees at 13.4-13.1 cal ka BP (González-Sampériz et al., 2006).

GI-1a (13.1-12.9) Early warming followed by a general cooling

Glacier advance or stillstand with moraine formation at 13.0 ka Cold and dry conditions in the Pyrenees starting at 13.1 cal ka BP (González-Sampériz et al., 2006).

http://mc.manuscriptcentral.com/jqs Journal of Quaternary Science