David Palacios
email: davidp@ucm.es

Antonio Gómez-Ortiz

Jesús Alcalá-Reygosa

Nuria Andrés

Marc Oliva

Luis M Tanarro

Ferran Salvador-Franch

Irene Schimmelpfennig

José M Fernández-Fernández

Laëtitia Leanni

José M Fernández- Fernández

Laëtitia Léanni

The challenging application of cosmogenic dating methods in residual glacial landforms: The case

Keywords: Sierra Nevada, cosmic-ray exposure dating, beryllium 10, Little Ice Age, Penultimate Glacial Cycle

come   L'archive ouverte pluridisciplinaire

Introduction

Scientific knowledge on the glacial evolution of the Iberian mountains has greatly advanced over the last decade based on the application of Cosmic-Ray Exposure (CRE) dating on moraines, erratic boulders and polished surfaces. This method allowed inferring spatiotemporal patterns of glaciation since the maximum extent of glaciers during the last Pleistocene glacial cycle, with a focus on the Last Glacial Maximum (LGM) [START_REF] García-Ruiz | La cronología del último ciclo glaciar en las montañas del sur de Europa[END_REF], Oldest Dryas (Palacios et al., 2017a) and Younger Dryas [START_REF] García-Ruiz | Mountain glacier evolution in the Iberian Peninsula during the Younger Dryas[END_REF] in the main Iberian ranges such as the Cantabrian Mountains [START_REF] Rodríguez-Rodríguez | Chronology of glaciations in the Cantabrian Mountains (NW Iberia) during the Last Glacial Cycle based on in situ-produced 10 Be[END_REF], Pyrenees [START_REF] Delmas | The last maximum ice extent and subsequent deglaciation of the Pyrenees: an overview of recent research[END_REF], Central System [START_REF] Domínguez-Villar | Early maximum extent of paleoglaciers from Mediterranean mountains during the LPG[END_REF], Iberian System [START_REF] Fernández-Fernández | Chronological and geomorphological investigation of fossil debris-covered glaciers in relation to deglaciation processes: A case study in the Sierra de La Demanda, northern Spain[END_REF] and Sierra Nevada [START_REF] Palacios | A Timing and new geomorphologic evidence of the last deglaciation stages in Sierra Nevada (southern Spain) Quat[END_REF].

The existence and distribution of landforms and deposits, their delimitation and mapping, generated by glacial advances prior to the last Pleistocene glaciation are known since the first geomorphological descriptions carried out in most Iberian mountains in the late 19 th and early 20 th centuries, such as in the Pyrenees [START_REF] Penck | La période glaciaire dans les Pyrénées[END_REF][START_REF] Nussbaum | Sur les traces des glaciers quaternaires dans la région de l Aragón[END_REF][START_REF] Barrère | Equilibrie glacier actuel et quaternaire dans l´Ouest des Pyrénèes Centrales[END_REF], Cantabrian Mountains [START_REF] Penck | Die Picos de Europa und das kantabrische Gebirge[END_REF][START_REF] Obermaier | Estudio de los glaciares de los Picos de Europa[END_REF], Central System [START_REF] Penck | Die Picos de Europa und das kantabrische Gebirge[END_REF][START_REF] Obermaier | Los glaciares cuaternarios de la Sierra de Guadarrama[END_REF] and Sierra Nevada [START_REF] Quelle | Beiträge zur Kenntnis der spanischen Sierra Nevada[END_REF][START_REF] Obermaier | Los glaciares cuaternarios de Sierra Nevada[END_REF].

However, the present knowledge on calendar of development of these landforms is scarce, and it is based on CRE dating only in very few cases. In fact, the only CRE data related to the penultimate glacial cycle were obtained with 21 Ne in the NW of the Iberian Peninsula [START_REF] Vidal-Romaní | The glaciation of Serra de Queixa-Invernadoiro and Serra do Gerês-Xurés, NW Iberia. A critical review and a cosmogenic nuclide ( 10 Be and 21 Ne) chronology[END_REF]. These authors found glacial evidence (moraine boulders and polished surfaces) of a maximum advance that occurred at 155 ± 30 ka 21 Ne age in the Serra de Queixa, whereas polished bedrock surfaces were dated at 231 ± 48 ka and 131 ± 31 ka 21 Ne ages in the Serra da Geres. In addition, erratic blocks were dated at 113.9 ± 7.1 ka 10 Be age in the Porma Valley, Cantabrian Mountains [START_REF] Rodríguez-Rodríguez | Chronology of glaciations in the Cantabrian Mountains (NW Iberia) during the Last Glacial Cycle based on in situ-produced 10 Be[END_REF]. Another erratic boulder, with a minimum 10 Be age of 122.2 ± 4.9 ka was found in Ariège valley (northern Pyrenees) [START_REF] Delmas | Palaeogeography and 10Be exposure-age chronology of Middle and Late Pleistocene glacier systems in the northern Pyrenees: implications for reconstructing regional palaeoclimates[END_REF]. This scarcity of CRE data from old moraines may be due to the problem of moraine degradation after being affected by intensive erosion, as it has been observed even in younger moraines of other regions [START_REF] Putkonen | Accuracy of cosmogenic ages for moraines[END_REF][START_REF] Briner | Cosmogenic exposure dating of late Pleistocene moraine stabilization in Alaska[END_REF][START_REF] Putkonen | Landscape evolution degrades the geologic signature of past glaciations[END_REF][START_REF] Balco | Contributions and unrealized potential contributions of cosmogenic-nuclide exposure dating to glacier chronology, 1990-2010[END_REF]Heyman et al., 2011).

The rest of dates corresponding to glacial landforms formed during the penultimate glaciation in the Iberian Peninsula were carried out by other dating methods. [START_REF] Villa | Age and significance of the Quaternary cemented deposits of the Duje Valley (Picos de Europa, Northern Spain)[END_REF] were able to date with U / Th the limestone of cemented calcareous breccia in Duje valley, Picos de Europa, also in the Cantabrian Mountains. This breccia deposit lies on glacially abraded surfaces and is covered by moraines dated between 394 ± 50 ka and 276 ± 23 ka / 394 ± 50 ka U / Th age. In the southern slope of the Central Pyrenees, Optically Stimulated Luminescence (OSL) dating was applied to sandy layers in fluvioglacial terraces in two nearby valleys that show the following ages: 151 ± 11 ka in the Gállego valley [START_REF] Lewis | Climatic implications of correlated Upper Pleistocene glacial and fluvial deposits on the Cinca and Gállego rivers (NE Spain) based on OSL dating and soil stratigraphy[END_REF] and 263 ± 21 ka and 171 ± 22 ka in the Aragón valley [START_REF] García-Ruiz | Glacial and fluvial deposits in the Aragón valley, central-western Pyrenees: chronology of the Pyrenean late pleistocene glaciers[END_REF].

In other Mediterranean mountains, landforms originated from glaciations prior to the Late Pleistocene Glaciation (LPG) have been dated by U / Th, 230 Th and OSL methods in the Balkan region (Hughes et al., 2006a(Hughes et al., , 2006b[START_REF] Hughes | Middle Pleistocene cold stage climates in the Mediterranean: new evidence from the glacial record[END_REF], Dinaric Alps [START_REF] Hughes | The glacial history of the Dinaric Alps, Montenegro[END_REF] and

Apennines [START_REF] Kotarba | On the age of Campo Imperatore glaciations, Gran Sasso Massif, Central Italy[END_REF], but not through CRE dating. In comparison to the large dataset of CRE ages existing for reconstructing the LPG in the Alps, there are only a few erratic boulders related to the penultimate glacial cycle dated by this method [START_REF] Graf | Multiple advances of Alpine glaciers into the Jura Mountains in the Northwestern Switzerland[END_REF], as it is also the case of the Himalayas [START_REF] Schaefer | Cosmogenic beryllium-10 and neon-21 dating of late Pleistocene glaciations in Nyalam, monsoonal Himalayas[END_REF]. On the contrary, the application of 230 Th-dates on the stalagmite record in the Alps allowed determining with high resolution the last maximum ice advance of the penultimate glacial cycle at 133.1 ± 0.7 and 131.9 ± 0.6 ka. These dates coincide with Heinrich stadial 11 and the onset of the penultimate deglaciation or Termination II (TII) which started at 131.8 ± 0.6 ka [START_REF] Häuselmann | Timing and nature of the penultimate deglaciation in a high alpine stalagmite from Switzerland[END_REF] and reached its maximum intensity at 130.9 ± 0.9 ka [START_REF] Moseley | Termination-II interstadial/stadial climate change recorded in two stalagmites from the north European Alps[END_REF].

However, CRE methods have never been applied to date the Little Ice Age (LIA) in the Iberian Peninsula, although this cold event has been extensively documented in this region [START_REF] González-Trueba | Little Ice Age glaciation and current glaciers in the Iberian Peninsula[END_REF]Oliva et al., 2018a). CRE dating showed that the LIA was the most extensive Holocene glacial advance in most Iberian mountains (Palacios et al., 2017b), with the exception of some cirques in the Central Pyrenees [START_REF] García-Ruiz | Holocene and 'Little Ice Age' glacial activity in the Marboré Cirque, Monte Perdido Massif, Central Spanish Pyrenees[END_REF]. In fact, after pioneer studies in New Zealand [START_REF] Schaefer | High frequency Holocene glacier fluctuations in New Zealand differ from the northern signature[END_REF] and the Andes [START_REF] Licciardi | Holocene glacier fluctuations in the Peruvian Andes indicate northern climate linkages[END_REF], the use of CRE dating methods to date landforms of Neoglacial and LIA age has just been recently applied in the Alps [START_REF] Schimmelpfennig | Holocene glacier culminations in the Western Alps and their hemispheric relevance[END_REF][START_REF] Schimmelpfennig | A chronology of Holocene and Little Ice Age glacier culminations of the Steingletscher, Central Alps, Switzerland, based on high-sensitivity beryllium-10 moraine dating[END_REF][START_REF] Roy | 10 Be exposure dating of the timing of Neoglacial glacier advances in the Ecrins-Pelvoux massif, southern French Alps[END_REF], Arctic [START_REF] Jomelli | Paradoxical cold conditions during the medieval climate anomaly in the Western Arctic[END_REF][START_REF] Young | Glacier maxima in Baffin Bay during the Medieval Warm Period coeval with Norse settlement[END_REF] and Central Asia mountains [START_REF] Dong | Cosmogenic 10 Be surface exposure dating of Little Ice Age glacial events in the Mount Jaggang area, central Tibet[END_REF][START_REF] Li | Cosmogenic 10 Be constraints on Little Ice Age glacial advances in the eastern Tian Shan, China[END_REF].

Radiocarbon dating on glacial, fluvioglacial and lacustrine sediments, as well as speleothems, tree rings and historical documents established the onset of the LIA in the Iberian Peninsula at 1300 common era (CE), with the coldest conditions of ca. 2 ºC below present-day values at 1620-1715 CE and the last cold episodes at 1760-1800 and 1815-1835 CE (Oliva et al., 2018a). In Sierra Nevada there were also small glaciers during the LIA that have been reconstructed based on radiocarbon dating of lake sediments [START_REF] Oliva | Late Holocene environmental dynamics and climate variability in a Mediterranean high mountain environment (Sierra Nevada, Spain) inferred from lake sediments and historical sources[END_REF] and through historical documents [START_REF] Gómez-Ortiz | Evidences from historical documents of landscape evolution after Little Ice Age of a Mediterranean high mountain area, Sierra Nevada, Spain (eighteenth to twentieth centuries)[END_REF], 2018).

The objective of this work is to apply CRE methods to glacial landforms originated during advances prior to LPG and to very young moraines as well as to review and discuss those found in this research in comparison to other similar studies. Sierra Nevada constitutes a unique environment to examine the use of CRE dating as a potential tool to reconstruct the chronology of old and young glaciations, due to the well-constrained glacial chronological record of the LPG and the existence of glacial landforms originated before and after of this glaciation in this mountains.

Study area

Sierra Nevada is located in the SE of the Iberian Peninsula and runs SW-NE parallel to the Mediterranean Sea (Fig. 1). This massif includes the highest peaks of the Iberian Peninsula:

Mulhacén (3478 m a.s.l.; 37°03'12"N, 3°18'41"W) and Veleta (3398 m; 37°03'02"N 3°20'54"W). The area is composed of Paleozoic metamorphic rocks, mainly micaschists, profoundly tectonized during the Alpine orogeny, with the presence of intensive foliation and joints that makes them susceptible to intense weathering [START_REF] Messerli | Beiträge zur Geomorphologie der Sierra Nevada (Andalusien)[END_REF][START_REF] Sanz De Galdeano | Nature and impact of the Neotectonic deformation in the western Sierra Nevada (Spain)[END_REF]. Climate conditions in Sierra Nevada are characteristic of a semiarid mid-latitude mountain environment, with a Mean Annual Air Temperature (MAAT) at 2016a). The MAAT at the top of the Veleta peak between 2002 and 2013 was 0.08 °C, with a slight increase of 0.12 °C along this period (Oliva et al., 2016a).

The Quaternary landscape evolution of Sierra Nevada has been largely studied over the last decades (Fig. 1, Table 1). From the first basic geomorphological observations of the late 19 th and early 20 th centuries to the recent multi-dating approaches, our understanding of the glacial processes shaping the landscape of the highest lands has significantly improved (i.e. [START_REF] Gómez-Ortiz | The deglaciation of Sierra Nevada (Spain): synthesis of current knowledge and new contributions[END_REF].

During the 1960s and 1970s several authors suggested the existence of older glaciations occurred prior to the LPG, formerly known as the Riss glaciation in Alpine terminology.

Their observations focused on the existence of highly eroded moraine remnants and fluvioglacial deposits at very low altitudes, mainly in the southern slope of the massif [START_REF] Hempel | Límites altitudinales geomorfológicos en Sierra Nevada[END_REF][START_REF] Messerli | Beiträge zur Geomorphologie der Sierra Nevada (Andalusien)[END_REF][START_REF] Lhenaff | Recherches géomorphologiques sur les Cordillères Bétiques centrooccidentales (Espagne)[END_REF][START_REF] Sánchez-Gómez | Aplicación del estudio de suelos a la dinámica de la cuenca del río Lanjarón. Relación suelos-geomorfología[END_REF]). However, this hypothesis has not been validated until now though direct dating methods. Subsequently, temperature increase recorded worldwide [START_REF] Clark | The Last Glacial Maximum[END_REF] favored a massive deglaciation of the massif, although two stages of glacial expansion occurred during the transition towards the Holocene, namely during the Oldest Dryas and Younger Dryas. During these phases, glaciers expanded significantly and flowed down-valleys several kilometers, both in southern and northern slopes [START_REF] Gómez-Ortiz | The deglaciation of the Sierra Nevada (southern Spain)[END_REF][START_REF] Palacios | A Timing and new geomorphologic evidence of the last deglaciation stages in Sierra Nevada (southern Spain) Quat[END_REF].

The Holocene was characterized by warmer temperatures that promoted the gradual melting of these glaciers during the Early Holocene [START_REF] Palacios | A Timing and new geomorphologic evidence of the last deglaciation stages in Sierra Nevada (southern Spain) Quat[END_REF]. The formerly glaciated environments were then subjected to paraglacial processes, which favoured the development of rock glaciers associated to permafrost conditions (Oliva et al., 2016b;2018b). These permafrost-related landforms finally stabilized at 6-7 ka during the Holocene Climate Optimum, when warmer conditions must have conditioned permafrost thawing (Palacios et al., 2016).

Climate variability intensified in the Northern Hemisphere during the Late Holocene [START_REF] Mayewski | Holocene climate variability[END_REF]. In Sierra Nevada, this was reflected in alternating cold and warm phases [START_REF] Oliva | Environmental evolution in Sierra Nevada (South Spain) since the LPG based on multi-proxy records[END_REF]. During the coldest stages, some glacier spots developed in the highest northern cirques, such as in the Mulhacén cirque, where lake sediments from La Mosca Lake revealed evidence of the existence of a glacier within the cirque at ca. 2.8-2.7 ka cal BP, 1.4-1.2 ka cal BP, and during the LIA between 1440 and 1710 CE [START_REF] Oliva | Late Holocene environmental dynamics and climate variability in a Mediterranean high mountain environment (Sierra Nevada, Spain) inferred from lake sediments and historical sources[END_REF]. For the last cold stage, historical documents also showed evidence of the prevailing colder conditions in Sierra Nevada, with descriptions, sketches and pictures of glaciers and permanent snow fields at the foot of the highest peaks of the westernmost mountains of Sierra Nevada between the late 17 th and the early 20 th century. The last glacier was located in the Veleta cirque and constituted the southernmost glacier in Europe during the LIA, although it finally melted away during the mid-20 th century [START_REF] Gómez-Ortiz | Location of permafrost in marginal regions: Corral del Veleta, Sierra Nevada, Spain[END_REF], 2009, 2018;[START_REF] Oliva | Late Holocene environmental dynamics and climate variability in a Mediterranean high mountain environment (Sierra Nevada, Spain) inferred from lake sediments and historical sources[END_REF]Oliva et al., 2018a). However, no CRE dating is yet available from the different LIA glacial stages in Sierra Nevada.

Methods

Sampling strategy

With the aim of dating glacial landforms formed potentially before the LPG, this research focuses on the Naute valley, a steep mountain ravine located on the southern slope of the Mulhacén peak (Fig. 1, 2, 3 and 4). The valley results from the confluence of two glacial headwaters -Río Seco and Mulhacén valleys-converging at 2300 m of altitude. In previous studies these two high mountain valleys were selected to reconstruct glacial stages of the LPG using CRE methods [START_REF] Gómez-Ortiz | The deglaciation of the Sierra Nevada (southern Spain)[END_REF][START_REF] Palacios | A Timing and new geomorphologic evidence of the last deglaciation stages in Sierra Nevada (southern Spain) Quat[END_REF], but landforms distributed in the lower Naute valley affected by intense erosion processes were not examined. In Río Seco and Mulhacén valleys landforms from the LPG have been already dated using the in situ cosmogenic nuclide 36 Cl , such as the front of a fossil rock glacier (9.0 ± 0.8 ka) resting on polished rock surfaces (11.9 ± 1.1 ka) in the Río Seco cirque, as well as of a fossil rock glacier (6.3 ± 0.5 ka and 13.1 ± 1.3 ka) resting on a polished rock surface (11.2 ± 1.1 ka and 12.7 ± 1.2 ka) at the head of the Mulhacén valley [START_REF] Palacios | A Timing and new geomorphologic evidence of the last deglaciation stages in Sierra Nevada (southern Spain) Quat[END_REF]. Just before their confluence, both valleys include a sequence of moraine ridges, where 5 blocks were dated using 36 Cl (29.9 ± 2.6 ka, 16.9 ± 1.8 ka, 14.9 ± 1.5 ka, 13.7 ± 1.3 ka, 11.7 ± 1.6 ka; [START_REF] Palacios | A Timing and new geomorphologic evidence of the last deglaciation stages in Sierra Nevada (southern Spain) Quat[END_REF]. These authors interpreted these ridges as a polygenic moraine formed by successive advances generated from just before the LGM to the Oldest Dryas. In the Naute valley, immediately below the confluence of these two tributaries, a lateral moraine is preserved in the left sector of the valley at elevations between 2150 and 2300 m. This moraine is placed below one of the ridges where a boulder was dated at 29.9 ± 2.6 ka [START_REF] Palacios | A Timing and new geomorphologic evidence of the last deglaciation stages in Sierra Nevada (southern Spain) Quat[END_REF], and therefore must be older in age, possibly from a previous glaciation. For this reason, this lower moraine was considered as one of the objectives of the work, in order to apply the CRE dating method.

To date glacial oscillations during the LIA by CRE dating methods, we focused our research on the Veleta cirque, located under the northern wall of Veleta peak (Fig. [START_REF] Schimmelpfennig | A chronology of Holocene and Little Ice Age glacier culminations of the Steingletscher, Central Alps, Switzerland, based on high-sensitivity beryllium-10 moraine dating[END_REF]5,6). Historical documents provided abundant evidence of a glacier existing at the floor of this cirque during the LIA until the first half of the 20 th century [START_REF] Gómez-Ortiz | Location of permafrost in marginal regions: Corral del Veleta, Sierra Nevada, Spain[END_REF], 2009, 2018). The cirque is closed by a large frontal moraine whose age is unknown, although due to its large size it has been considered to have formed during the Last Termination or during Holocene advances [START_REF] Oliva | Environmental evolution in Sierra Nevada (South Spain) since the LPG based on multi-proxy records[END_REF](Oliva et al., , 2018a)).

Sampling and analytical procedures

Seven samples were collected from moraine boulders using a hammer and a chisel: four located in Naute valley (NAUT-1, -2, -3, -4) between 2164 and 2263 m, and three in the Veleta cirque (SN-11-1, -2, -3) between 3061 and 3095 m (Table 1). All samples were taken from the flat-topped surfaces of boulders > 1 m high, located on the crests of the moraines.

Field data and sample characteristics of the seven samples are listed in Table 2.

Physical and chemical sample preparation and beryllium measurements were carried out at the Centre Européen de Recherche et d Enseignement des Géosciences de l Environnement (CEREGE, France). Lichen, moss and other organic fragments were removed from the samples with a brush. Samples were crushed with a roller grinder and sieved to retrieve the grain size fraction 0.25-1 mm. All samples are quartz-bearing micaschists, and therefore we selected in situ-produced cosmogenic 10 Be dating to determine the ages of the moraines. To isolate the quartz from the bulk rock, the samples went repeatedly through a magnetic separator (Frantz LB-1) until all magnetic minerals were discarded. Subsequently, the nonmagnetic fraction experienced successive chemical attacks with a mixture of concentrated hydrochloric (HCl) and hexafluorosilic (H 2 SiF 6 ) acids to dissolve the remaining non-quartz minerals. The residual impurities were dissolved when the sample grains were decontaminated from meteoric 10 Be by three successive partial dissolutions with concentrated hydrofluoric acid (HF).

The subsequent beryllium extraction protocol is adapted from [START_REF] Brown | Examination of surface exposure ages of Antarctic moraines using in-situ produced 10 Be and 26 Al[END_REF] and [START_REF] Merchel | An update on radiochemical separation techniques for the determination of longlived radionuclides via Accelerator Mass Spectrometry[END_REF] chemical procedures. The samples yielded between 18 and 32 g of purified quartz (Table 3). About 100 µl of a 9 Be carrier solution with a 9 Be concentration of dissolved in HF. A chemistry blank was prepared along with the seven samples. Following evaporation of the resulting solution, the samples were recovered in a hydrochloric acid solution and beryllium was precipitated to Be(OH) 2 with ammonia before and after elution through an anionic exchange column (Dowex 1X8) to remove iron. Following the methods described by [START_REF] Merchel | An update on radiochemical separation techniques for the determination of longlived radionuclides via Accelerator Mass Spectrometry[END_REF], a cationic exchange column (Dowex 50WX8) was used to remove boron and to separate the Be from other elements. Beryllium was precipitated with ammonia to Be(OH) 2 and the resulting precipitate was oxidized to BeO at 700 °C. Then, this final BeO was mixed with niobium powder and loaded on cathodes for analysis of the 10 Be/ 9 Be ratios at the French national Accelerator Mass Spectrometry (AMS) facility ASTER (CEREGE) [START_REF] Arnold | The French accelerator mass spectrometry facility ASTER: improved performance and developments[END_REF]. The measurements were calibrated against in-house standard STD-11 using an assigned 10 Be/ 9 Be ratio of 1.191 (±0.013) x 10 -11 [START_REF] Braucher | Preparation of ASTER in-house 10 Be/ 9 Be standard solutions[END_REF]. Sample 10 Be/ 9 Be ratios were corrected for the chemical blank background by subtracting the measured chemistry blank 10 Be/ 9 Be ratio (Table 3). Analytical 1 sigma uncertainties include uncertainties in AMS counting statistics, the uncertainty in the standard 10 Be/ 9 Be ratio, an external AMS error of 0.5% [START_REF] Arnold | The French accelerator mass spectrometry facility ASTER: improved performance and developments[END_REF], and the uncertainty in the chemical blank measurement. A 10 Be half-life of 1.387 (± 0.01) x 10 6 years was used [START_REF] Chmeleff | Determination of the 10 Be half-life by multicollector ICP-MS and liquid scintillation counting[END_REF][START_REF] Korschinek | A new value for the half-life of 10Be by heavy-ion elastic recoil detection and liquid scintillation counting[END_REF].

Age calculation ( 10 Be age computation)

We calculated the 10 Be CRE ages with the online CREP exposure age calculator [START_REF] Martin | The CREp program and the ICE-D production rate calibration database: a fully parameterizable and updated online tool to compute cosmic-ray exposure ages[END_REF]; http://crep.crpg.cnrs-nancy.fr), using the LSD scaling model [START_REF] Lifton | Scaling in situ cosmogenic nuclide production rates using analytical approximations to atmospheric cosmic-ray fluxes[END_REF] that is similar to other previous empirical models [START_REF] Borchers | Geological calibration of spallation production rates in the CRONUS-Earth project[END_REF]. Age calculations have also considered the ERA40 atmospheric model [START_REF] Uppala | The ERA-40 reanalysis[END_REF] and LSD Framework geomagnetic database [START_REF] Lifton | Scaling in situ cosmogenic nuclide production rates using analytical approximations to atmospheric cosmic-ray fluxes[END_REF]. As there is no regional production rate available, we used the worldwide mean 10 Be spallation production rate of 3.99 ± 0.22 atoms g - 1 yr -1 , as calibrated in the ICE-D production rate database linked to CREp. For all samples, a rock density of 2.7 g cm -3 was considered and the topographic shielding factor of each sample was calculated (Table 2).

In the moraine of the Veleta cirque, from where samples SN-11-1, SN-11-2 and SN-11-3

were extracted, snow cover remains for a large part of the year [START_REF] Gómez-Ortiz | Evidences from historical documents of landscape evolution after Little Ice Age of a Mediterranean high mountain area, Sierra Nevada, Spain (eighteenth to twentieth centuries)[END_REF], and there is information on the properties and persistence of the snow cover for the last decades [START_REF] Herrero | Evaposublimation from the snow in the Mediterranean mountains of Sierra Nevada (Spain)[END_REF], which allows considering snow shielding effect for these samples. The correction for snow cover has been calculated applying the equation by [START_REF] Gosse | Terrestrial in situ cosmogenic nuclides: theory and application[END_REF], including local parameters [START_REF] Herrero | Evaposublimation from the snow in the Mediterranean mountains of Sierra Nevada (Spain)[END_REF], and considering a snow thickness on moraine boulders ranging from 30 to 130 cm during 8 months per year. We used an average value of snow density of 0.3 g cm -3 , considering that the attenuation length for fast neutrons in snow is 109 g cm -2 [START_REF] Zweck | Snow shielding factors for cosmogenic nuclide dating inferred from Monte Carlo neutron transport simulations[END_REF][START_REF] Delunel | Snow shielding factors for cosmogenic nuclide dating inferred from longterm neutron detector monitoring[END_REF]. Taking into account these characteristics, we applied a snow shielding factor of 0.865 to the affected samples. We have not applied a snow shielding factor to the old samples of Naute moraine, since the coverage of snow in this area is much shorter and data are not available for a period encompassing several tens of thousands of years.

Table 3 also includes the ages determined with version 3 of the online exposure age calculator formerly known as the CRONUS-Earth online exposure age calculator [START_REF] Balco | A complete and easily accessible means of calculating surface exposure ages or erosion rates from 10 Be and 26 Al measurements[END_REF][START_REF] Balco | CRONUS-Earth online exposure age calculator[END_REF]. These ages have been calculated using the default calibration data set based on the ICE-D calibration database and the time-dependent "LSDn" scaling method [START_REF] Lifton | Scaling in situ cosmogenic nuclide production rates using analytical approximations to atmospheric cosmic-ray fluxes[END_REF]. The average difference between the ages calculated with the CREP exposure age calculator and those calculated with the CRONUS-Earth (v.3) is 2.4% for NAUT samples and 6.9% for SN-11 samples.

The ages of the SN-11 samples calculated with the snow shielding factor turn out to be between 8.8 and 16.9% higher in the online CREP exposure age calculator [START_REF] Martin | The CREp program and the ICE-D production rate calibration database: a fully parameterizable and updated online tool to compute cosmic-ray exposure ages[END_REF] and between the 7.9 and 16.9% higher in the online exposure age calculator formerly known as the CRONUS-Earth online exposure age calculator v. 3 [START_REF] Balco | A complete and easily accessible means of calculating surface exposure ages or erosion rates from 10 Be and 26 Al measurements[END_REF][START_REF] Balco | CRONUS-Earth online exposure age calculator[END_REF]. The CREP exposure calculator ages and the analytical uncertainties are used in the text and in the figures of this manuscript and, only in the case of the Veleta cirque, with the application of the snow shielding factor (Table 3).

Results

Geomorphological analysis of the Naute moraine and CRE dating results

To select the most suitable boulders for CRE dating, an exhaustive survey was carried out of the Naute moraine system, formed potentially during a glacial advance before the LPG. The results of this survey are expressed in a detailed geomorphological map (Fig. 2). The moraine is located on the left bank of the Naute river (west) and the Peñón Grande (east) between 2150 and 2300 m. It has been intensely reshaped by postglacial environmental dynamics: (i) fluvial processes driven by steep streams have deeply eroded the moraine ridge, and (ii)

gullying processes have washed away a large part of its abundant fine sediments. As a result, a large number of boulders -highly weathered and fragmented due to their abundant foliation planes and joints-are found on the surface of the moraine (Fig. 3).

To sample the boulders that best represent the original surface of the moraine, we focused on blocks that: (i) present aligned upper edges indicating the maximum possible height of the original moraine crest, and (ii) are inserted into the moraine, which ensures that they have not been mobilized since the moraine stabilization. Using these criteria, only four boulders were found with these characteristics (NAUT-1, 2, 3 and 4) (Fig. 4). The rest of the boulders showed traces of having been affected by erosion and exhumed.

Samples from two boulders showed a similar age: NAUT-1 yielded 134.8 ± 6.2 ka and NAUT-4 129.2 ± 5.3 ka. The other two boulders showed different and substantially younger ages: NAUT-2 obtained 53.7 ± 2.9 ka and NAUT-4 43.8 ± 2.3 ka (Figs 2, 3 and Table 3).

Geomorphological analysis of the Veleta cirque moraine and CRE dating results

The Veleta cirque moraine is composed of a large proportion of fine sediments with a few large blocks on the surface, possibly deposited by different glacial advances. Within the area enclosed by this moraine, there are large rock avalanche deposits fallen from the Veleta northern wall and a small rock-glacier where buried glacial ice is still preserved under the debris cover, though it presents evidence of accelerated degradation. On the other hand, the fine sediments of the moraine, resulting from heavily weathered micaschists, are being washed away by runoff activity. This process is particularly intense in late spring and early summer with snowmelt runoff, and in autumn with torrential rain events. The outer slope of the moraine is being reshaped by frequent debris flows, which also facilitate the outcropping and exhumation of large blocks buried in the moraine. In addition, there are some metric-size boulders fallen from the Veleta north wall and deposited on the surface of the moraine, normally during the winter season when snow completely fills the floor of the cirque. These rocks slide on the snow and ice surface and stabilize on the moraine ridge. On the inner slope of the moraine, accumulations are found of superimposed till blocks of different sizes. In contrast to the main moraine, these till deposits include a smaller proportion of fine particles but conserve a recent glacial imprint as revealed by the presence of flutes. This deposit exceeds the limits of the main moraine on the eastern side of the cirque, where the moraine height is lower (Fig. 5 and6). Evidence from historical sources confirms that the distribution of this recent till deposit coincides with the surface covered by the glacier during the 19 th century (Fig. 7; [START_REF] Bide | Deuxième excursion dans la Sierra Nevada[END_REF][START_REF] Gómez-Ortiz | The geographical interest of historical documents to interpret the scientific evolution of the glacier existing in the Veleta cirque (Sierra Nevada, Spain) during the Little Ice Age[END_REF].

Before selecting boulders for CRE dating, a very detailed geomorphological survey of the large polygenic moraine system enclosing the Veleta cirque was conducted. A detailed geomorphological map shows the main landforms observed in the cirque (Fig. 5). Schematic transects show the sectors of the moraine affected by current or recent post-glacial processes and derived landforms (Fig. 6). Considering these geomorphological characteristics, we avoided collecting samples from blocks corresponding to very recently deglaciated surfaces because: (i) we already know that this area was occupied by the glacier during the 19 th century and (ii) they can result from the frequent rockfalls from the headwall. As a result, only three CRE datable boulders were found anchored in the moraine and protruding sufficiently not to have been covered by other sediments since their deposition (SN11-1, 2 and 3) (Fig. 7). Thus, they cannot have fallen from the headwall, or have been covered by the glacier in a subsequent readvance. and Table 3).

Discussion

The existence of a glaciation prior to the last Pleistocene glaciation in Sierra Nevada

The results of 10 Be CRE dating yielded ages of ca. 130-135 ka for the two samples collected from boulders in the moraine downslope of the maximum ice advance in Sierra Nevada [START_REF] Gómez-Ortiz | The deglaciation of the Sierra Nevada (southern Spain)[END_REF]. These ages suggest that the moraine may have been formed

during the last maximum advance of the penultimate glaciation. The other two boulders of the same moraine yield much younger ages, around 45-50 ka (Table 3), and thus can be considered outliers. On the other hand, none of the studies in Sierra Nevada reported boulder samples affected by nuclide inheritance [START_REF] Gómez-Ortiz | The deglaciation of the Sierra Nevada (southern Spain)[END_REF][START_REF] Delmas | The last maximum ice extent and subsequent deglaciation of the Pyrenees: an overview of recent research[END_REF][START_REF] Palacios | A Timing and new geomorphologic evidence of the last deglaciation stages in Sierra Nevada (southern Spain) Quat[END_REF], and therefore this problem can be rejected in the older samples.

In Sierra Nevada, no surfaces polished and striated bedrock surfaces where found in this study between the moraine ridges dated ca. 130-135 ka and those dated in previous studies between ca. 30 and ca. 15 ka [START_REF] Palacios | A Timing and new geomorphologic evidence of the last deglaciation stages in Sierra Nevada (southern Spain) Quat[END_REF] , 2016, 2017), suggest using a minimum of 5 boulders from the same moraine to ensure obtaining the formation age of the landform. Thus, where most of the results show dispersed data, the ages are rejected, whereas an average value is obtained with the rest, interpreted as indicative of the development age of the moraine. In this study we attempted to apply this five-boulder criterion, but this proved impossible due to the deterioration of the moraine since the penultimate glaciation, as the landform has been intensely eroded and few potentially datable boulders are available. The headwaters of several torrents have eroded the ridge summit, leaving the moraine covered with loose blocks. When trying to reconstruct the original surface of the moraine we found only four embedded boulders (Fig. 4).

The CRE ages of ca. 130-135 ka obtained in Sierra Nevada coincide with the last cold period of the penultimate glaciation in the Alps as suggested by 230 Th dating on speleothems that obtained 131.8 ± 0.6 ka [START_REF] Häuselmann | Timing and nature of the penultimate deglaciation in a high alpine stalagmite from Switzerland[END_REF][START_REF] Moseley | Termination-II interstadial/stadial climate change recorded in two stalagmites from the north European Alps[END_REF]. This stage occurred during Heinrich stadial 11 [START_REF] Oppo | Evolution and demise of the last interglacial warmth in the subpolar North Atlantic[END_REF], parallel to the end of the last cold phase in the Alboran sea surface temperatures (near Sierra Nevada, in SE Iberia) before the beginning of TII [START_REF] Martrat | Similarities and dissimilarities between the last two deglaciations and interglaciations in the North Atlantic region[END_REF][START_REF] Jiménez-Amat | Offset timing of climate oscillations during the last two glacial-interglacial transitions connected with large-scale freshwater perturbation[END_REF]. By that time, the onset of warmer temperatures favored the end of the last cold period of the penultimate glaciation in Greenland [START_REF] Grant | Rapid coupling between ice volume and polar temperature over the past 150,000 years[END_REF] as well as the initial melting of large ice sheets at high latitudes in the Northern Hemisphere [START_REF] Govin | Sequence of events from the onset to the demise of the Last Interglacial: Evaluating strengths and limitations of chronologies used in climatic archives[END_REF], leading to sea level rise [START_REF] Rohling | Differences between the last two glacial maxima and implications for icesheet, 18 O, and sea-level reconstructions[END_REF].

The formation and stabilization of the Naute moraine may therefore be contemporary with what is called the Penultimate Glacial Maximum (PGM), a parallel concept to the LGM for the last glaciation, centered on 140 ka [START_REF] Rohling | Differences between the last two glacial maxima and implications for icesheet, 18 O, and sea-level reconstructions[END_REF]. The deposition age of the Naute moraine is similar to the stabilization age of the maximum glacial advance moraine in the Sierra de Queixa dated 155 ± 30 ka 21 Ne [START_REF] Vidal-Romaní | The glaciation of Serra de Queixa-Invernadoiro and Serra do Gerês-Xurés, NW Iberia. A critical review and a cosmogenic nuclide ( 10 Be and 21 Ne) chronology[END_REF]. The erratic boulder dated 113.9 ± 7.1 ka 10 Be in the Porma Valley, Cantabrian Mountains, could also be related to the deglaciation process following the PGM [START_REF] Rodríguez-Rodríguez | Chronology of glaciations in the Cantabrian Mountains (NW Iberia) during the Last Glacial Cycle based on in situ-produced 10 Be[END_REF]. The same is true for the erratic boulder in the Ariège valley (northern Pyrenees), with a minimum 10 Be age of 122.2 ± 4.9 ka [START_REF] Delmas | Palaeogeography and 10Be exposure-age chronology of Middle and Late Pleistocene glacier systems in the northern Pyrenees: implications for reconstructing regional palaeoclimates[END_REF], slightly after the formation of the Naute moraine from the penultimate glaciation. In addition, fluvioglacial terraces deposited simultaneously on the outermost moraine ridges in the Gállego and Aragón valleys, Central Pyrenees, were dated 151 ± 11ka and 171 ± 22 ka OSL ages by other authors [START_REF] Lewis | Climatic implications of correlated Upper Pleistocene glacial and fluvial deposits on the Cinca and Gállego rivers (NE Spain) based on OSL dating and soil stratigraphy[END_REF][START_REF] García-Ruiz | Glacial and fluvial deposits in the Aragón valley, central-western Pyrenees: chronology of the Pyrenean late pleistocene glaciers[END_REF], and therefore may represent the same glacial stage that led to the formation of the Naute moraine. However, despite some evidence suggesting glacial activity between ca. 170 and 130 ka in the Iberian Peninsula, the small number of reliably dated terrestial records makes it difficult to infer spatiotemporal glaciation patterns occurring before the LPG. This is a common issue in other Mediterranean mountains where glacial advances concurrent with the development of the Naute moraine have only been dated through the application of U series to secondary carbonates in cemented moraines in mountains in Greece (Hughes et al., 2006a(Hughes et al., , b, 2007)), Dinaric Alps [START_REF] Hughes | Middle Pleistocene cold stage climates in the Mediterranean: new evidence from the glacial record[END_REF] and the Gran Sasso Massif, Central Italy [START_REF] Kotarba | On the age of Campo Imperatore glaciations, Gran Sasso Massif, Central Italy[END_REF].

A case study similar to the Naute moraine was undertaken in the Jura mountains, where several erratic boulders were found scattered beyond the Alpine glaciers LGM extent. In this case the oldest ages were very similar to those of the Naute moraine (129.7 ± 4.7 ka, 143.2 ± 8.2 ka and 144.0 ± 5.3 ka 10 Be ages), although just as in Sierra Nevada, other boulders in the same deposit yielded younger ages [START_REF] Graf | Multiple advances of Alpine glaciers into the Jura Mountains in the Northwestern Switzerland[END_REF]. After their experience applying CRE dating to several moraine boulders in the Himalayas, [START_REF] Schaefer | Cosmogenic beryllium-10 and neon-21 dating of late Pleistocene glaciations in Nyalam, monsoonal Himalayas[END_REF] highlighted the role of erosion processes in the exhumation of blocks in old moraines. As a result, these authors claim that the oldest, not the average age of the moraine boulders reflects the best minimum age estimation for the moraine deposition, in agreement with previous authors [START_REF] Hallet | Surface Dating of dynamic landforms: Young boulders on aging moraines[END_REF]. [START_REF] Briner | Cosmogenic exposure dating of late Pleistocene moraine stabilization in Alaska[END_REF] proposed the same criterion after analyzing a large number of moraine boulders, presumably from the penultimate glaciation in Alaska, and found many problems of boulder exhumations. These authors also defend an oldest-age method to provide the closest approximation to the moraine stabilization age. Similarly, and from previous experience, [START_REF] Balco | Contributions and unrealized potential contributions of cosmogenic-nuclide exposure dating to glacier chronology, 1990-2010[END_REF] defined the most appropriate strategies when applying CRE methods to degraded moraines, which have been considered in this study.

From our results we conclude that the difficulty of finding suitable boulders to apply CRE methods in moraines from old glaciations should not discourage their use. Dating an adequate number of samples is very difficult for moraines prior to the LPG, where erosive activity over time may have exhumed most of the boulders. In the case of Sierra Nevada, the profoundly tectonized micaschists have been intensely affected by erosion, which makes it more difficult to find suitable boulders for dating [START_REF] Gómez-Ortiz | The deglaciation of Sierra Nevada (Spain): synthesis of current knowledge and new contributions[END_REF]. Nevertheless, Sierra Nevada is an arid mountain, where ice margin deposits preserved on gentle slopes make it possible to obtain a well-dated glacial sequence using the CRE method [START_REF] Pallàs | Small, isolated glacial catchments as priority target for cosmogenic surface dating of Pleistocene climate fluctuations[END_REF] (Fig. 9).

Chronology of LIA moraine formation in Sierra Nevada

It was also difficult to find suitable blocks for CRE dating in young moraines in Sierra Nevada, as they are found in very dynamic geomorphological environments on steep slopes of 200 m high (up to 300 m from the peaks to the cirque floors) with thick and long-lasting snow cover. On the one hand, these factors accelerate erosion processes in the moraine, often reducing its size, vertical development and the slope gradient; on the other hand, they favor new sediment deposits on the moraine surface. These difficulties are also observed in most

Mediterranean mountains, where LIA glaciers did not exceed the limits of the cirques and developed at the foot of vertical walls [START_REF] Hughes | Little Ice Age glaciers in the Mediterranean mountains[END_REF][START_REF] Hughes | Little Ice Age glaciers and climate in the Mediterranean mountains: a new analysis[END_REF]. Previous studies adverted to the polygenic character of the Veleta cirque moraine, related to the effects of an Early Holocene glacial advance that pushed downslope, accumulating rockfall deposits from the end of TI [START_REF] Gómez-Ortiz | The deglaciation of the Sierra Nevada (southern Spain)[END_REF][START_REF] Palacios | A Timing and new geomorphologic evidence of the last deglaciation stages in Sierra Nevada (southern Spain) Quat[END_REF]. According to these authors, this moraine system functioned as a barrier to subsequent advances, which accumulated sediments on the inner slope, generating a single, large, polygenic moraine.

Using the criteria defined in Methods above, we obtained three coherent results to infer the age of the moraine, although dates were obtained for only three boulders. The high uncertainty levels in our results prevent us from relating the age of each sample to LIA periods (Oliva et al., 2018a). These ages could be related to advances from 1300 CE to 1610/1680 CE, contemporary with the LIA maximum glacial advance in the Alps [START_REF] Holzhauser | Glacier and lake-level variations in west-central Europe over the last 3500 years[END_REF] and the lowest temperatures in the Iberian Peninsula during the Maunder Minimum (Oliva et al., 2018a).

The inner slope of the moraine therefore developed during the LIA, although the large moraine system may also include boulders from other Neoglacial glacial stages, such as those detected in the Mulhacén cirque, where lake sediments show evidence of a glacier within the cirque ca. 2.8-2.7 ka, 1.4-1.2 ka cal BP and LIA [START_REF] Oliva | Late Holocene environmental dynamics and climate variability in a Mediterranean high mountain environment (Sierra Nevada, Spain) inferred from lake sediments and historical sources[END_REF].

New results obtained in Mount Olympus (Greece) found similar results using 36 Cl dating method, with one sample of 0.64 ± 0.08 ka [START_REF] Styllas | Lateglacial and Holocene history of the northeast Mediterranean mountains-New insights from in situ-produced 36Cl-based cosmic ray exposure dating of paleo-glacier deposits on Mount Olympus, Greece[END_REF]. However, most studies focusing on the use of CRE dating methods for dating LIA glacial advances have been carried

out in the Alps. In the Western Alps, LIA moraines are very close to Neoglacial (Late Holocene) moraines and may even overlap them [START_REF] Schimmelpfennig | Holocene glacier culminations in the Western Alps and their hemispheric relevance[END_REF][START_REF] Hughes | Little Ice Age glaciers in the Mediterranean mountains[END_REF]; as in the Veleta cirque, LIA glacial advances have left a single polygenic ridge with dated boulders showing a wide range of dates: 1430 ± 32, 1534 ± 28, and 1829 ± 11 CE 10 Be (ages not considering snow shielding) [START_REF] Schimmelpfennig | Holocene glacier culminations in the Western Alps and their hemispheric relevance[END_REF]. These dates could correspond to the three most important glacial advances in the Alps during the LIA: 1300-1380, 1600-1670 and 1800-1860 CE [START_REF] Holzhauser | Glacier and lake-level variations in west-central Europe over the last 3500 years[END_REF], similar to those detected in the Pyrenees (Oliva et al., 2018a). Recent studies in the Central Alps showed a similar pattern when dating LIA moraine boulders -also very close to or overlapping Neoglacial moraines-forming single polygenic crests [START_REF] Schimmelpfennig | A chronology of Holocene and Little Ice Age glacier culminations of the Steingletscher, Central Alps, Switzerland, based on high-sensitivity beryllium-10 moraine dating[END_REF]. These authors dated 14 boulders in a single LIA moraine ridge using 10 Be, assuming negligible inheritance for the boulders, and reported dates ranging from 1430 to 1870 CE [START_REF] Schimmelpfennig | A chronology of Holocene and Little Ice Age glacier culminations of the Steingletscher, Central Alps, Switzerland, based on high-sensitivity beryllium-10 moraine dating[END_REF]. Le [START_REF] Roy | 10 Be exposure dating of the timing of Neoglacial glacier advances in the Ecrins-Pelvoux massif, southern French Alps[END_REF] dated Neoglacial moraines near several present-day glacier fronts in the French Alps and obtained 10 Be ages that are all consistent with the Late Holocene period (~4-1 ka), but do not follow the logical chronostratigraphic moraine sequence. These authors highlighted the problems related to surface exhumation and erosion of many moraine ridges.

Some attempts to date the moraines of the LIA have been also made in other high mountain ranges such as the Cordillera Vilcabamba, Central Andes where [START_REF] Licciardi | Holocene glacier fluctuations in the Peruvian Andes indicate northern climate linkages[END_REF] applied

Be CRE dating for the first time. Samples collected from a single massive ridge showed average dates ranging from 1740 ± 30 to 1810 ± 20 CE, though samples of old dates were rejected as being outliers. They suggest that the moraines found in several valleys of this massif may be the result of the successive accumulation of multiple Late Holocene glacial expansions. Due to the short exposure duration of our samples, their 10 Be ages have 25-30% analytical uncertainties. These uncertainties are higher than the age difference of 10-16% obtained when a snow shielding correction factor is applied, as was proposed in early studies [START_REF] Benson | The probable importance of snow and sediment shielding on cosmogenic ages of north-central Colorado Pinedale and pre-Pinedale moraines[END_REF][START_REF] Schildgen | Simulation of snow shielding corrections for cosmogenic nuclide surface exposure studies[END_REF]. In any case, this study defends the importance of considering the snow cover of the moraines studied, both when applying CRE methods, and to discover the degree of degradation to which they have been subjected. The location of the glacial fronts during the LIA suggests MAAT ca. 1 ºC below current values in mid-latitude mountain regions (Oliva et al., 2018a), determining longer snow cover duration in the highlands of the glaciated massifs and strengthening shielding from cosmogenic radiation (Fig. 7). Longlasting snow cover also intensifies nivation-related erosion processes and can intensely affect glacial deposits [START_REF] Christiansen | Nivation forms and processes in unconsolidated sediments, NE Greenland. Earth Surface Processes and Landforms[END_REF][START_REF] Palacios | Distribution and effectiveness of nivation in Mediterranean mountains: Peñalara (Spain)[END_REF].

On the other hand, as we observed in the Veleta cirque, and as described in several papers cited above, the intensity of the paraglacial processes [START_REF] Ballantyne | Paraglacial geomorphology[END_REF][START_REF] Oliva | Coupling patterns between paraglacial and permafrost degradation responses in Antarctica[END_REF] in the current deglaciation phase from the LIA advance is very high, with a complex superposition of processes that destroy the glacial landforms and accumulate new deposits on them. Therefore, as we saw, it was difficult to find suitable boulders for CRE dating in the Veleta cirque, despite its young age. In this context, the long-term study and monitoring of the dynamics of the LIA moraines in full paraglacial phase, should serve to evaluate the time needed for moraine stabilization, once it has been abandoned by glacial retreat (Zreda andPhillips, 1994,1995;[START_REF] Putkonen | Degradation of unconsolidated Quaternary landforms in the western North America[END_REF][START_REF] Putkonen | Landscape evolution degrades the geologic signature of past glaciations[END_REF]Heyman et al., 2011) (Fig. 9).

In areas where LIA moraines are constrained by steep slopes, as is very common in Alpine systems, at least two, or possibly three, major LIA advances overlap, in some cases with other Neoglacial moraines, to form a single polygenic ridge, as occurred in the Veleta cirque. This was highlighted in the mid-1980s with the concepts of obliterative overlap [START_REF] Gibbons | Probability of moraine survival in a succession of glacial advances[END_REF] and distal-flank accretion [START_REF] Osborn | Lateral-moraine stratigraphy and neoglacial history of Bugaboo Glacier, British Columbia[END_REF] applied to Neoglacial advances.

Conclusions

The results of this research supply evidence of the very few CRE dates available on very old glacial landforms formed before the last Pleistocene glaciation, and also the few very young dates available, such as those derived from LIA glaciation, compared with the large number of dates obtained for Termination I. For both periods, the selection of boulders suitable for sampling was the critical issue in Sierra Nevada, due to intense postglacial environmental dynamics such as occurred during the last major deglaciation, which impedes the collection of a statistically significant number of samples. In the case of the LIA, moraines are still undergoing intense paraglacial readjustment to the new geomorphological setting, with very intense erosion and sediment redistribution of the unconsolidated moraine particles. In addition, moraine accretion processes during successive glacial advances are detected from geomorphological observations. In moraines originating during the penultimate glaciation, the long time lapse since their formation has facilitated intense erosion and reshaping of these landforms, with very few stable blocks remaining on the surface since their original stabilization. Periods of certain geomorphic stability may have occurred between the very old and the younger glaciations (Fig. 9).

The results of this research have enabled us to verify the existence of moraines deposited in glaciations prior to LPG in Sierra Nevada, as many of the first geomorphologists who studied the glacial landforms of this mountains suggested, and confirm similar CRE and OSL ages established in other mountains in the northern Iberian Peninsula.

In addition, CRE dates are obtained for the first time in the Iberian Peninsula in LIA moraines, with results suggesting that maximum advances during this period occurred between the 14 th and 17 th centuries, as is the case in many other mountains where the same method has been applied. Be sample exposure ages.

Table 1. Current knowledge of Quaternary glacial stages in Sierra Nevada.

Stage Chronology Processes and landforms

Pre-Last Glaciation Unknown

Possible existence of eroded moraines from glaciations occurred before the Last Glaciation as well as glacio-fluvial sediments distributed at lower elevations than moraines formed during the Last Glaciation [START_REF] Hempel | Límites altitudinales geomorfológicos en Sierra Nevada[END_REF][START_REF] Messerli | Beiträge zur Geomorphologie der Sierra Nevada (Andalusien)[END_REF][START_REF] Lhenaff | Recherches géomorphologiques sur les Cordillères Bétiques centrooccidentales (Espagne)[END_REF][START_REF] Sánchez-Gómez | Aplicación del estudio de suelos a la dinámica de la cuenca del río Lanjarón. Relación suelos-geomorfología[END_REF].

Last Glaciation (Maximum Ice Extent; MIE)

Two glacial advances at ca. 32-30 ka and 20-19 ka Development of alpine valleys at elevations between 2000 and 2500 m in northern and southern slopes, respectively [START_REF] Gómez-Ortiz | The deglaciation of the Sierra Nevada (southern Spain)[END_REF], 2013[START_REF] Delmas | The last maximum ice extent and subsequent deglaciation of the Pyrenees: an overview of recent research[END_REF][START_REF] Oliva | Environmental evolution in Sierra Nevada (South Spain) since the LPG based on multi-proxy records[END_REF][START_REF] Palacios | A Timing and new geomorphologic evidence of the last deglaciation stages in Sierra Nevada (southern Spain) Quat[END_REF].

The glaciated environment during the MIE encompassed 105 km 2 [START_REF] Palma | Spatial characterization of glacial and periglacial landforms in the highlands of Sierra Nevada (Spain)[END_REF].

Deglaciation

Two glacial advances during the Oldest Dryas (OD) at ca. 17 ka and Younger Dryas (YD) at ca. 13-12 ka Glacial retreat at ca. 19 ka followed by two phases of glacial development until the Holocene, mainly in the highest valleys from the northern slope of the massif [START_REF] Gómez-Ortiz | The deglaciation of the Sierra Nevada (southern Spain)[END_REF], 2013;[START_REF] Oliva | Late Holocene environmental dynamics and climate variability in a Mediterranean high mountain environment (Sierra Nevada, Spain) inferred from lake sediments and historical sources[END_REF][START_REF] Oliva | Environmental evolution in Sierra Nevada (South Spain) since the LPG based on multi-proxy records[END_REF][START_REF] Palacios | A Timing and new geomorphologic evidence of the last deglaciation stages in Sierra Nevada (southern Spain) Quat[END_REF].

Holocene

Glacial retreat until the Early Holocene at ca. 10-9 ka Three stages with formation of glaciers occurred during the Late Holocene at ca. 2.8-2.7, 1.4-1.2 ka cal BP and LIA

The glaciers formed during the YD finally melted and paraglacial activity favoured the development of rock glaciers in these areas [START_REF] Gómez-Ortiz | The deglaciation of the Sierra Nevada (southern Spain)[END_REF], 2013;Oliva et al., 2016b;[START_REF] Palacios | A Timing and new geomorphologic evidence of the last deglaciation stages in Sierra Nevada (southern Spain) Quat[END_REF].

Lake sediments suggest the development of a small glacier in the Mulhacén cirque, which is reflected in several moraine ridges distributed across the cirque floor [START_REF] Oliva | Late Holocene environmental dynamics and climate variability in a Mediterranean high mountain environment (Sierra Nevada, Spain) inferred from lake sediments and historical sources[END_REF]Oliva et al., 2016a).

Little Ice Age (LIA)

From 1440 to 1710 AD a glacier existed in the Mulhacén cirque, and another glacier developed in the Corral de Veleta cirque until the mid-20 th century

In the highest northern cirques between Mulhacén and Veleta peaks, historical sources and geomorphic evidence (i.e. moraines) shows evidence of the existence of several small glaciers that gradually disappeared during the 19 th and 20 th centuries [START_REF] Gómez-Ortiz | Location of permafrost in marginal regions: Corral del Veleta, Sierra Nevada, Spain[END_REF], 2009, 2018;[START_REF] Oliva | Late Holocene environmental dynamics and climate variability in a Mediterranean high mountain environment (Sierra Nevada, Spain) inferred from lake sediments and historical sources[END_REF]Oliva et al., 2018a). a Carrier has a concentration of 3025 µg 9 Be/g. b Sample 10 Be/ 9 Be ratios corrected for batch-specific analytical blank ratio of (4.5 ± 1.3) x 10 -15 . c 10 Be ages calculated with the online CREp exposure age calculator [START_REF] Martin | The CREp program and the ICE-D production rate calibration database: a fully parameterizable and updated online tool to compute cosmic-ray exposure ages[END_REF]. 1sigma errors include analytical and production rate uncertainties; analytical errors shown in brackets.

d 10
Be ages calculated with the online exposure age calculator formerly known as the CRONUS-Earth online exposure age calculator version 3 [START_REF] Balco | A complete and easily accessible means of calculating surface exposure ages or erosion rates from 10 Be and 26 Al measurements[END_REF]. 1sigma errors include analytical and production rate uncertainties; analytical errors shown in brackets.

  By contrast, the LPG has been recently well-constrained with CRE techniques. Gómez-Ortiz et al. (2012) and Palacios et al. (2016) reported two major glacial advances occurred at ca. 30-32 ka and 19-20 ka. The second stage of glacial advance occurred in phase to the global Last Glacial Maximum and reached almost the same extent than during the previous glacial stage.

  Sample SN-11-1, collected on the top of the moraine, very close to 20 th century till deposits, obtained an age of 340 ± 120 years, corresponding to the year 1675 ± 120 CE. Sample SN-11-2, collected on the crest of the lobate ridge of the moraine, yielded an age of 720 ± 270 years (1295 ± 270 CE). Finally, sample SN-11-3, collected on the top of the lateral moraine ridge in the western part of the cirque, obtained an age of 400 ± 110 years (1615 ± 110 CE) (Figs 5, 6

Figure captions Figure 1 .

 captions1 Figure captions

Figure 2 .

 2 Figure 2. Geomorphological map of the Naute valley and location of samples. We show CRE

Figure 3 .

 3 Figure 3. Results of the Naute moraine samples of this work, together with pictures of some

Figure 4 .

 4 Figure 4. Geomorphological sketch of the two moraines existing in the lower part of the

Figure 5 .

 5 Figure 5. Geomorphological map of the Veleta cirque and sample location. We show CRE

Figure 6 .

 6 Figure 6. (A) Schematic transect of the Veleta cirque including the different landforms

Figure 7 .

 7 Figure 7. (A) Photograph of Corral of the Veleta cirque from the West in September, 2016.

Figure 8 .

 8 Figure 8. (A) Veleta cirque with its 19 th century moraine (source: Google Earth), (B) Vertical

Figure 9 .

 9 Figure 9. Sketch representing the erosion processes affecting a moraine: (A) Intense mass

Figure

  Figure (Color) Click here to download high resolution image

  , i.e. rock surfaces which may have been icecovered prior to the LPG. All the polished bedrock surfaces dated in Sierra Nevada -in many

cases with abundant and well-preserved fresh striations-correspond to the deglaciation process during Termination I

[START_REF] Gómez-Ortiz | The deglaciation of the Sierra Nevada (southern Spain)[END_REF][START_REF] Palacios | A Timing and new geomorphologic evidence of the last deglaciation stages in Sierra Nevada (southern Spain) Quat[END_REF]

. The fact that moraines from a previous glaciation are found in close proximity to those of the last glacial cycle without intermediate polished surfaces led us to apply CRE dating only on boulder surfaces in the old moraines. Recent research in mountain areas, such as in the Alps (e.g.

[START_REF] Chenet | 10 Be cosmic-ray exposure dating of moraines and rock avalanches in the Upper Romanche valley (French Alps): Evidence of two glacial advances during the Late Glacial/Holocene transition[END_REF]

, the Andes (e.g.

[START_REF] Martini | Late Pleistocene glacial fluctuations in Cordillera Oriental, subtropical Andes[END_REF]

, Central Asia (e.g.

[START_REF] Batbaatar | Asynchronous glaciations in arid continental climate[END_REF]

and the Iberian Peninsula (Rodríguez-Rodríguez et al.

  [START_REF] Schaefer | High frequency Holocene glacier fluctuations in New Zealand differ from the northern signature[END_REF] also applied10 Be dating in the Mueller valley in the Southern Alps, New Zealand, to date many LIA moraine boulders. The geomorphological setting is similar to that described in the European Alps, with a single LIA moraine system overlapping with various Late Holocene moraine ridges. Boulders taken from different sites in the LIA moraine reported dates of 1350 ± 60, 1600 ± 50, 1780 ± 40 and 1820 ± 20 CE 10 Be, which coincide with cold periods inferred from tree-ring data. Li et al. (2016) applied 10 Be dating in LIA moraines in Tian Shan mountains, Central Asia; interestingly, they found that boulders in moraines from glaciers smaller than 1.0 km 2 show very old ages because of nuclide inheritance. The same problems were found in the small Veleta cirque, where the wall is very close to the moraine, but we discarded this possibility following the results obtained in this present study and those from other previous work in Sierra Nevada, because of the intense, continuous rock fall activity on its wall. For larger glaciers, Li et al. (2016) found evidence of proposed a recent glacial maximum during the Medieval Climate Anomaly, when glaciers in Europe receded. Jomelli et al. (2016) applied in situ cosmogenic 36 Cl to date three moraines in Lyngmarksbraeen glacier (West Greenland), which were deposited during the last millennium in a relatively flat area not constrained by topography, the opposite that happens in Alpine environments. These authors obtained advances with average ages of 1200 ± 130 CE, 1450 ± 90 CE and 1720 ± 60 CE for the most external and internal ridges. The most recent advances coincide with glacial expansion in European high mountains, but the older moraine may have developed during the Medieval Climate Anomaly (Jomelli et al., 2016). In the sub-Antarctic Kerguelen islands, Jomelli et al. (2017) recently dated a moraine with 36 Cl that is probably related to the LIA.

	a major advance ca. 1600 ± 100 CE 10 at ca. 1480 ± 55 CE 10 Be. Dong et al. (2017) dated two different LIA ridges in Tibetan Be, and also a remarkable early LIA glacial expansion mountains with 10 Be, collecting four samples from each ridge. Their results provide a range of 1480±139 to 1975±31 CE 10 Be dates for the entire moraine system, very similar to the Veleta cirque results. Young et al. (2015) dated several stable moraine ridges close to the current alpine glacier snouts in Baffin Island and western Greenland and found a logical geomorphological order from 1040±40 CE to 1700±40 CE 10 Be ages with a retreat from 1750 we concluded that, in general, a first glacial advance was recorded during the 14 th century, with the major ice expansion during the 17 th century and a last minor readvance at the CE, and From a review of the available literature on LIA moraine dating in different mountain ranges, beginning of the 19 th century.

Table captions Table 1 .

 captions1 Current knowledge of Quaternary glacial stages in Sierra Nevada.

	Table 2. Field data and characteristics of sampled boulders dated with 10	Be in Sierra Nevada
	(*boulders are covered by snow thickness of 30 to 130 cm during 8 months per year). See text
	for details.

Table 3 .

 3 Analytical data and10

Table Table 2

 Table . Field data and characteristics of sampled boulders dated with10 Be in Sierra Nevada (*boulders are covered by snow thickness of 30 to 130 cm during 8 months per year). See text for details.

	Table
	Sample	Latitude (ºN)	Longitude (ºW)	Elevation (m a.s.l.)	Thickness (cm)	Topographic shielding factor	Snow shielding factor
	NAUT-1	37.0196111	3.33047222	2164	1.5	0.97832
	NAUT-2	37.0201111	3.33036111	2179	1.5	0.98414
	NAUT-3	37.0212500	3.32977778	2211	1.5	0.98221
	NAUT-4	37.0225556	3.32866667	2263	3.0	0.97817
	*SN-11-1	37.0595833	3.36569444	3076	2.5	0.80979	0.70067
	*SN-11-2	37.0605833	3.36747222	3061	2.0	0.80979	0.70067
	*SN-11-3	37.0606944	3.36897222	3095	2.0	0.80979	0.70067

Table 3 .

 3 Analytical data and 10 Be sample exposure ages

			9 Be
	Sample	Quartz (g)	carrier solution
			(mg)

a 10 Be/ 9 Be (x10 13 ) b [ 10 Be] (x10 3 at.g -1 ) 10 Be age (CREp c ) (ka) 10 Be age (C-E d ) (ka) 10 Be age (CREp c ) (ka) With snow shielding factor 10 Be age (C-E d ) (ka)

									With snow shielding
									factor
	NAUT-1	18.32	98.9	23.3 ± 1.0	2530 ± 110	134.8 ± 10 (6.2)	130.4 ± 10 (6.1)
	NAUT-2	21.98	100.2	11.2 ± 0.58	1030 ± 53	53.7 ± 4.3 (2.9)	54.5 ± 4.3 (2.9)
	NAUT-3	22.53	100.5	9.62 ± 0.55	868 ± 50	43.8 ± 3.2 (2.3)	43.0 ± 3.6 (2.5)
	NAUT-4	22.62	100.7	28.6 ± 1.2	2580 ± 110	129.2 ± 8.9 (5.4)	125.4 ± 9.4 (5.5)
	SN-11-1	30.91	100.9	0.127 ±0.052	8.4 ± 3.5	0.31 ± 0.10 (0.10)	0.28 ± 0.11 (0.11)	0.34 ± 0.12 (0.12)	0.32 ± 0.13 (0.13)
	SN-11-2	32.79	100.7	0.285 ± 0.086	17.7 ± 5.3	0.60 ± 0.22 (0.21)	0.59 ± 0.20 (0.18)	0.72 ± 0.28 (0.27)	0.71 ± 0.22 (0.22)
	SN-11-3	31.39	101.0	0.159 ± 0.050	10.3 ± 3.3	0.36 ± 0.09 (0.09)	0.35 ± 0.11 (0.10)	0.40 ± 0.11 (0.11)	0.38 ± 0.12 (0.12)
	Blank		99.4	0.0450 ± 0.013

Table Figure (Color) Click here to download high resolution image

 Figure

Acknowledgments

The research was carried out within the MOUNTAIN WARMING project (CGL2015-65813-R). Marc Oliva is grateful for the support of the Ramón y Cajal research program (RYC-2015-17597) and the ANTALP research group (Antarctic, Arctic, Alpine Environments, 2017-SGR-1102). The 10 Be measurements were performed at the ASTER AMS National facility (CEREGE, Aix en Provence) which is supported by the INSU/CNRS, the ANR h gh he "P je h ma i e d e cellence" g am f he "E i emen d e cellence" ASTER-CEREGE action and IRD. The authors express their deep gratitude to Dr. Magali Delmas and two anonymous reviewers whose detailed and interesting suggestions have helped to improve our manuscript.