# Lateglacial paleoglacier and paleoclimate reconstructions in the north-western Italian Alps

Elena Serra, Fabio Magrani, Pierre Valla, Natacha Gribenski, Julien

Carcaillet, David Lundbek Egholm

► To cite this version:

Elena Serra, Fabio Magrani, Pierre Valla, Natacha Gribenski, Julien Carcaillet, et al. Lateglacial paleoglacier and paleoclimate reconstructions in the north-western Italian Alps. Quaternary Science Reviews, 2022, 298, pp.107822. 10.1016/j.quascirev.2022.107822. hal-03927106

HAL Id: hal-03927106

https://hal.science/hal-03927106

Submitted on 8 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Lateglacial paleoglacier and paleoclimate reconstructions in the north-western Italian Alps

<figure><img src="image_2.png" /></figure>

Elena Serra<sup>a, *</sup>, Fabio Magrani<sup>a</sup>, Pierre G. Valla<sup>b, a</sup>, Natacha Gribenski<sup>a</sup>, Julien Carcaillet<sup>b</sup>, David Lundbek Egholm<sup>c</sup>

## ARTICLE INFO

**Article history:**

Received 4 July 2022

Received in revised form

3 October 2022

Accepted 11 October 2022

Available online 9 November 2022

**Keywords:**

Lateglacial

¹⁰Be cosmogenic dating

Glacier modelling

Paleoclimate reconstruction

Western European Alps

## ABSTRACT

Lateglacial (19.0–11.7 ka ago) paleo-temperature records in the Western European Alps document several short-term cooling episodes within the general post-Last Glacial Maximum warming trend, in phase with North Hemisphere climatic oscillations. Alpine paleo-precipitation reconstructions are instead rare, and further constraints are needed in order to assess whether Lateglacial cold periods were associated with a modified atmospheric circulation pattern over the European Alps. The Alpine paleo-glacial record offers a quantitative framework to investigate Lateglacial paleoclimatic conditions and glacier sensitivity in response to both temperature and precipitation changes. Through the combination of ¹⁰Be surface-exposure dating of glacial landforms and deposits constraining ice front and surface, together with numerical glacier simulations (iSOSIA), our study aims to reconstruct in space and time the main Lateglacial ice stages and associated paleo-climatic conditions in three tributary valleys (Valpelline, Valsavarenche and Val di Cogne) located within the Dora Baltea catchment (north-western Italian Alps). Our dating-modelling approach reveals in all the three investigated sectors two distinct paleoglacier stages at ca. 13 and 11 ka, documenting respectively the ice configurations at the transition between the Oldest Dryas cold period and the Bølling-Allerød interstadial, and between the Younger Dryas cold period and the early Holocene warming. Numerical ice-simulation outcomes suggest a similar-to-today precipitation pattern during the two ice stages, with either same absolute or homogeneously decreased precipitation values over the Dora Baltea catchment, although a unique quantitative solution of paleo-precipitation magnitudes could not be constrained from our dataset. Using present-day precipitation pattern and magnitude, our results provide paleo-temperature offsets from present-day between –3.3 and –3.8 °C for the older glacial stage, and –2.7 to –3.2 °C for the younger glacial stage, coinciding with the upper range of paleo-temperature reconstructions from other paleoclimatic proxies. Finally, Alpine glaciers' sensitivity to climate fluctuations differ significantly between the investigated catchments, with a much higher glacier sensitivity to changes in the Equilibrium Line Altitude in the northern (Valpelline) compared to the southern (Valsavarenche and Val di Cogne) tributaries, reflecting different topographic and/or climatic conditions between the three valleys.

© 2022 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license

(http://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Paleoclimate reconstructions in the European Alps suggest major changes in atmospheric circulation during the Last Glacial Maximum (LGM, 26.5–19.0 ka ago; Clark et al., 2009), with dominant southwesterly moisture advection from the Mediterranean

(Florineth and Schlüchter, 2000; Kuhlemann et al., 2008; Luetscher et al., 2015; Becker et al., 2016; Monegato et al., 2017). Such enhanced precipitation in south-western Europe (Ludwig et al., 2016) was proposed to result from the topographic anomaly of the expanding European and Laurentide ice sheets in response to global temperature cooling, causing a southward migration of the North Atlantic storm track (Merz et al., 2015), and in turn triggering the last maximum expansion of Alpine glaciers (Wirsig et al., 2016; Monegato et al., 2017; Gribenski et al., 2021).

<footnote><sup>a</sup> Institute of Geological Sciences and Oeschger Centre for Climate Change Research, University of Bern, Switzerland</footnote>

<footnote><sup>b</sup> Institute of Earth Sciences (ISTerre), Université Grenoble Alpes, Université Savoie Mont Blanc, CNRS, IRD, IFSTAR, France</footnote>

<footnote><sup>c</sup> Department of Geoscience, Aarhus University, Denmark</footnote>

<footnote>* Corresponding author.
E-mail address: elena.serra@geo.unibe.ch (E. Serra).</footnote>
Post-LGM gradual retreat of the North Hemisphere ice sheets and associated northward migration of the North Atlantic storm track induced the re-establishment of north-westerly circulation over the Alps, with modern dominant moisture coming from the Atlantic (Merz et al., 2015) already in place at ca. 22 ka (Luetscher et al., 2015). Reduced moisture supply from the Mediterranean Sea, together with globally-increasing temperatures (Rasmussen et al., 2014), initiated general ice-retreat from the Alpine foreland at 24–19 ka and ice-thinning in the internal Alpine massifs since ca. 18 ka (Wirsig et al., 2016 and reference therein; Lehmann et al., 2020). However, overall glacier decay during the Lateglacial period (19.0–11.7 ka ago) was not continuous but interrupted by multiple stages of stillstand or re-advance (so-called Alpine Lateglacial stadials; Ivy-Ochs, 2015 and reference therein; Federici et al., 2016; Protin et al., 2021; Serra et al., 2022) associated to short-term periods of climatic deterioration (i.e. general cooling trend and enhanced climatic variability), recorded by paleoclimate proxies in the European Alps (see review in Heiri et al., 2014). These short-term cooling episodes punctuating the Lateglacial were recognized as analogous to the global stadials (GS) or cold interstadial sub-events (GI) of the Greenland ice core stratigraphy (Rasmussen et al., 2014). They include (1) a first post-LGM climatic deterioration associated to the Heinrich 1 ice-rafting event in the North Atlantic at 17.5–15.4 ka (GS-2.1a; Stanford et al., 2011), (2) several short cold events interrupting the Bølling-Allerød interstadial between 14.6 and 12.8 ka (~3–4 °C warming equivalent to GI-1; Heiri et al., 2014) with temperature drops of 0.5–1.5 °C (Older Dryas or Aegelsee Oscillation at 14.0–13.9 ka, equivalent to GI-1d; Gerzensee Oscillation at 13.3–13.0 ka, equivalent to GI-1b; Lotter et al., 2012), and (3) a last and abrupt cooling of 1.5–3 °C during the Younger Dryas (12.9–11.7 ka; GS-1), before the onset of the Holocene (11.7 ka-present; Heiri et al., 2014).

While detailed paleo-temperature records are available for the Lateglacial (Heiri et al., 2014 and reference therein), Alpine paleo-precipitation reconstructions are rare (Luetscher et al., 2015; Protin et al., 2019; Baroni et al., 2021) and further constraints are therefore needed in order to assess whether Lateglacial cold periods were associated with spatial variations in Alpine paleo-precipitation. Such a spatial re-organization in atmospheric circulation, similar to the LGM, has been suggested for the Younger Dryas globally (Bakke et al., 2009) and over Europe (Rea et al., 2020). The Alpine paleoglacial record offers a quantitative framework to further investigate Lateglacial paleoclimatic conditions and paleo-precipitation changes (e.g. Kerschner and Ivy-Ochs, 2008; Davis et al., 2009; Protin et al., 2019; Baroni et al., 2021), since glacier mass-balance is highly sensitive to both temperature and precipitation (Oerlemans, 2005) as reflected by the Equilibrium Line Altitude (ELA; i.e. the elevation to which annual ablation equals annual accumulation). However, uncertainties still remain due to the low spatial and temporal resolution of the Alpine Lateglacial records, and to the difficulty in decoupling the relative contribution of past temperature and precipitation changes from paleoglacier fluctuations (e.g. Protin et al., 2019; Baroni et al., 2021).

With the aim of assessing Alpine paleoglacier sensitivity to Lateglacial paleoclimate and potential changes in paleo-precipitation, our study focuses on the post-LGM glacial history of three tributary valleys within the Dora Baltea catchment (north-western Italian Alps; Fig. 1), namely Valpelline, Valsavarenche and Val di Cogne. We combined $^{10}$Be surface-exposure dating of glacial landforms and deposits together with numerical glacier simulations in order to constrain the paleoglacier extent and thickness during Lateglacial ice stages in the three studied tributary valleys. Using iSOSIA numerical simulations (e.g. Egholm et al., 2011) to model steady-state paleoglacier configurations, we derived paleoclimate estimates (i.e. past changes in precipitation and

temperature) and paleo-ELAs from our paleoglacier reconstructions in order to discuss the sensitivity of different Alpine paleoglaciers in the Dora Baltea catchment to climate forcing (Kerschner and Ivy-Ochs, 2008).

## 2. Study area

Our study focuses on three major tributary valleys of the Dora Baltea catchment: Valpelline, Valsavarenche and Val di Cogne (Fig. 1). The Dora Baltea catchment is a large drainage system (around 3400 km² and 170 km long) located in the north-western Italian Alps, with the three studied valleys being situated north (Valpelline) and south (Valsavarenche and Val di Cogne) of the catchment, and connecting to the main Dora Baltea stream ~50 km downstream of its source. The three catchments are around 20–30 km long, with total surface areas of ~150 km² (Valsavarenche) and 260–270 km² (Val di Cogne, Valpelline). They extend from the Dora Baltea valley floor downstream (~600 m a.s.l.) to major 4000-m Alpine peaks in the valley heads (i.e. Dent d'Hérens and Gran Paradiso massifs; Fig. 1), with an average altitude of 2350–2500 m a.s.l.

Present-day temperature conditions are similar in the three studied catchments, but with a significant mean annual temperature gradient across each catchment due to the large elevation range encompassed, varying from 9 °C (valley floor) to –6 °C (4000 m-high mountain peaks), and with ~8 °C of seasonal amplitude variation (Regione Autonoma Valle d'Aosta, 2019). Annual precipitation is comparatively low at the outlets of all the three tributary catchments (around 800 mm/yr), while moderately wetter conditions are observed in the uppermost sectors of Valpelline (1480 mm/yr) compared to Valsavarenche and Val di Cogne (1220 mm/yr; Isotta et al., 2014; Fig. S1A). Modern glaciers cover 10–13% of the studied catchment areas and are concentrated in the valley heads, above 2000 m a.s.l. (GlaRiskAlp Project, http://www.glariskalp.eu). Mass-balance data series for the Grand Etret glacier (Valsavarenche, Fig. 4A for location) for the 1999–2019 period document an ELA at around 2900 m a.s.l. (data from Parco Nazionale del Gran Paradiso, Fig. S2B), in line with other modern glaciers within the Dora Baltea catchment (Baroni et al., 2021).

Throughout the Quaternary, valley glaciers from Valpelline, Valsavarenche and Val di Cogne repeatedly fluctuated and connected with the extensive Dora Baltea glacial system during periods of major glaciations, when the Dora Baltea glacier occupied the entire catchment down to the Po Plain. The Ivrea morainic amphitheatre is a clear evidence for these periods of maximum ice extension (Fig. 1) and was lastly abandoned after the LGM (Gianotti et al., 2015). Following the LGM, the Valpelline, Valsavarenche and Val di Cogne tributaries were all disconnected from the retreating Dora Baltea glacier after the Alpine Gschnitz stadial (ca. 15 ka; Serra et al., 2022). During the Younger Dryas (ca. 12 ka) cooling event, paleoglaciers remained confined in the upper source catchments (Baroni et al., 2021). Subglacial deposits and ice-margin morainic landforms from the Little Ice Age (LIA; 1250–1860 CE), representing the maximum Holocene extent reached by local glaciers (Baroni et al., 2021), are well preserved in all the studied catchments and were mapped within the GlaRiskAlp Project (http://www.glariskalp.eu, Figs. 3–5).

## 3. Methods

### 3.1. $^{10}$Be surface-exposure dating

We performed remote-sensing based geomorphological mapping (high-resolution DEMs and orthophotos from Regione Autonoma Valle d'Aosta) cross-checked with fieldwork
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Study area and sampling locations within the Dora Baltea catchment (mosaic DEM from Regione Autonoma Valle d'Aosta, Regione Piemonte, swisstopo, and Institut Geographique National). Black line marks the extent of the Dora Baltea catchment within the mountain front. Solid green, red and blue lines delimit the Valpelline, Valsavarenche and Val di Cogne tributary catchments, respectively, which are investigated in the present study. Red and yellow dots indicate new sampling locations (glacially-polished bedrocks and morainic boulders, respectively) from this study. Green dots show morainic boulders from Baroni et al. (2021) used for discussion. Black boxes highlight the extent of Figs. 3–5. Present-day glaciers (GlaRiskAlp Project, http://www.glariskalp.eu), main topographic peaks and main rivers are indicated. Bottom-left inset shows location of the Dora Baltea catchment (yellow open box) within the European Alps, with the LGM ice extent (light blue; Ehlers and Gibbard, 2004). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

observations in the upper catchments of Valpelline, Valsavarenche and Val di Cogne (Fig. 1) to identify potential glacial landforms and deposits marking past glaciers geometries and suitable for ¹⁰Be surface exposure dating. In total, twenty-four glacially-polished bedrocks and five morainic boulders were sampled (Fig. 2 and Table 1). Our sampling strategy was driven both by field access restrictions and/or preservation of glacial landforms and deposits, and by our objective to acquire geomorphic constraints on both glacier extent and thickness. To constrain paleoglacier ice-front fluctuations where no frontal moraines have been preserved, such as in Valpelline and Valsavarenche valleys, we targeted glacially-polished surfaces on bedrock knobs along the valley floor (VALP01, 03, 11, 12; VSAV14, 15; Figs. 3 and 4). In the upper part of Valsavarenche, we also targeted erratic boulders on the top of a high-elevation latero-frontal moraine (VSAV06-11; Figs. 2C and 4), to complement the existing dating of morainic boulders (Valsavarenche and Val di Cogne, Baroni et al., 2021, Fig. 4A and B and 5A). In addition, glacially-polished bedrock surfaces were sampled along the valley flanks, following altitudinal transects from the valley bottom to just below the trimline (Fig. 2A and B) with the aim to provide ice-thickness constraints (e.g. Wirsig et al., 2016). In total, five altitudinal transects were performed in the three catchments (VALP04-06; VALP07-11; VSAV01-04; VSAV05, 08, 12, 13; and COGNE01, 03, 04; Fig. 2A and B, 3–5 and S3).

Samples were collected on sub-vertical polished bedrock surfaces distant from soil coverage or on top of >1-m height morainic

boulders. Coherent surfaces with evidence for minimal weathering were selected. Crushing and sieving were performed in order to separate the 250–400 µm grainsize fraction, from which pure quartz was isolated by following modified procedure based on Kohl and Nishiizumi (1992). Quartz purification was performed at the Institute of Geological Sciences (University of Bern, Switzerland). For ¹⁰Be extraction through ion-exchange chromatography, conventional chemical treatments from Brown et al. (1991) and Merchel and Herpers (1999) were adapted and performed at the GeoThermoChronology platform (ISTerre - University Grenoble Alpes, France). ¹⁰Be/⁹Be ratios were measured at the ASTER French National AMS facility (CEREGE, Aix-en-Provence, France; Arnold et al., 2010) and calibrated against the in-house Be standard (isotope ratio 1.191 × 10¹¹; Braucher et al., 2015). Full process blank ¹⁰Be/⁹Be ratios of 5.4 ± 0.6 × 10⁻¹⁵ and 5.5 ± 1.0 × 10⁻¹⁵ were used to correct ¹⁰Be concentrations of VALP and VSAV/COGNE samples, respectively.

The online CREp program (Martin et al., 2017; https://crep.otelo.univ-lorraine.fr/#/init) was used in order to calculate the ¹⁰Be surface-exposure ages of our samples (Table 1), and to re-calculate, for consistency, previously-published ¹⁰Be surface-exposure ages (Baroni et al., 2021; see Figs. 4 and 5 for locations and Table S2 for details). A ¹⁰Be production rate by neutron spallation at sea-level and high-latitude (SLHL) of 4.16 ± 0.10 at g⁻¹ yr⁻¹ (Claude et al., 2014) was used and scaled at the sampling location with the LSDn scaling scheme (Lifton et al., 2014). Corrections for
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Field photographs of sampling locations for ¹⁰Be surface-exposure dating. A) Glacially-polished bedrock altitudinal transect targeted in Valpelline tributary catchment (right valley side; Fig. 3B for location). Light blue dashed line represents mapped trimline. B) Glacially-polished bedrock altitudinal transect targeted in Valsavarenche tributary catchment (left valley side; Fig. 4B for location). Light blue dashed line represents mapped trimline. C) Latero-frontal moraine of Valsavarenche tributary catchment (Fig. 4B for location; orange dashed line shows the morainic ridge) with sampled morainic boulders. D) Highest-elevation polished bedrock from the altitudinal transect targeted in Val di Cogne tributary catchment (Fig. 5 for location). Sample altitudes are indicated, and sample colour code is same as in Fig. 1. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

atmospheric pressure according to the ERA-40 reanalysis data set (Uppala et al., 2005) and for geomagnetic field fluctuations according to the Lifton-VDM2016 geomagnetic database (Lifton, 2016) were integrated in the scaling scheme. Field measurements were used to calculate topographic shielding correction at each sampling location based on Dunne et al. (1999). ¹⁰Be surface-exposure ages were corrected for an erosion rate of 0.1 mm ka⁻¹ on the basis of field observations of surface striations and low intergranular relief, suggesting low erosion rates and consistent with previous works conducted in the Dora Baltea catchment (Wirsig et al., 2016; Serra et al., 2022). No snow correction was applied due to the uncertainties in estimating temporal and spatial variability of snow accumulation in mountainous settings. While snow-cover shielding for glacially-polished bedrock samples is likely negligible, as these are sub-vertical bedrock surfaces collected along steep valley sides, ¹⁰Be surface-exposure ages from morainic boulders and valley-bottom bedrock knobs might be influenced by long-term snow shielding, and need therefore to be considered as minimum exposure ages. Calculated ¹⁰Be surface-exposure ages will be reported in the following with external uncertainties, given the large spatial distribution of our samples between the three different tributary valleys.

In order to identify potential distinct glacial stages, we assessed correlation between ice-front and ice-surface constraints based on

clusters from normal kernel density estimates (KDEs; based on Lowell, 1995) of individual ¹⁰Be surface-exposure ages, and on sample elevation and geomorphic/stratigraphic information. For each identified cluster, individual ¹⁰Be age KDEs were summed (Lowell, 1995), with mode and standard deviation of the summed KDE defining the time period of the glacial stage. Samples at different locations along the valley and/or different elevations (especially glacially-polished bedrock surfaces), but grouped in one ¹⁰Be age cluster, would show the magnitude of glacier retreat/thinning shortly after the ice stage. We also considered that the upper age range of each cluster represents the onset of glacier retreat/thinning and thus a minimum time estimate for the glacial stage.

## 3.2. Paleoglacier modelling

### 3.2.1. Ice-flow model

Our paleoglacier simulations are based on the ice-flow model iSOSIA (depth-Integrated Second Order Shallow Ice Approximation), a finite volume solver with explicit time integration that allows the computation of ice flow and subglacial hydrology under climate forcing (Egholm et al., 2011, 2012). The numerical integration of “higher-order” effects ensures the preservation of good accuracy compared to Full-Stokes model predictions for steep and
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** $^{10}$Be surface-exposure ages from Valpelline tributary catchment (location in Fig. 1; modified 2 m-resolution DEM from Regione Autonoma Valle d'Aosta; LIA glacier extent from GlaRiskAlp Project, http://www.glariskalp.eu). **A)** Glacially-polished bedrock sample locations in the upstream part of Valpelline. The most downstream valley floor polished-bedrock knob is indicated (VALP01) with its respective $^{10}$Be surface-exposure age. Black box highlights the extent of the zoom presented in Fig. 3B, green line depicts the extent of the studied catchment. **B)** Locations and $^{10}$Be surface-exposure ages of samples collected along the valley floor and altitudinal transects. **C)** Individual (dashed lines) and summed (continuous lines) kernel density estimates (KDEs) of $^{10}$Be surface-exposure ages. The modes and uncertainties of the summed KDEs are reported in bold black font. Colour code (valid also for sample names and ages reported in panels A and B) highlights the two age clusters (red and blue labels, see main text for details). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

rugged topography (Egholm et al., 2011), while allowing large-scale simulations for paleoglacial questions. The ice dynamics is based on Glen's flow law (stress exponent $n = 3$) and depth-averaged horizontal flow velocities depend non-linearly on ice thickness, ice-surface gradients and horizontal stress gradients in the longitudinal and transverse directions of ice flow (Egholm et al., 2011). In our simulations, we followed the approach adopted by Ugelvig et al. (2018) for modelling ice sliding and subglacial hydrology, although we ignored short-term daily transients in water pressure and considered only steady-state solutions for subglacial hydrology and effective pressure (see details and parameters in Magrani et al., 2022).

In iSOSIA, the climatic input is based on a simple mass-balance approach using a Positive-Degree-Day model (PDD), which is a function of the mean-annual temperature, the amplitude of the annual temperature variation, and the mean annual precipitation rate. Snow avalanching in iSOSIA (Scherler and Egholm, 2020) transports new accumulated snow downwards based on a given threshold in ice-surface slope. Any change in ice thickness is therefore computed as a balance between ice flow, ice ablation and accumulation, while mass conservation assumes a constant ice density spatially and at depth. Ice fluxes are computed by vertically-integrating horizontal ice-flow velocities at cell boundaries.

### 3.2.2. Modelling set-up

In the present study, we performed numerical simulations to assess the climatic conditions (temperature and precipitation) associated with the glacial stages reconstructed from geochronology and paleoglacier spatial configuration (section 3.1). We first calibrated the PDD mass-balance model against modern mass-balance observations from two glaciers (Fig. S2): the Argentière glacier (north-western side of Mont Blanc massif; 1976–2019 data

series from Glacioclim network, https://glacioclim.osug.fr) and the Grand Etrèt glacier (Valsavarenche, Fig. 4A; 1999–2019 data series from Annual Glaciological Surveys of Parco Nazionale Gran Paradoiso). We imposed in the PDD model modern sea-level temperature of 14.8 °C (NOAA, 2016), annual precipitation rates of 2 m/yr and 1.2 m/yr for Argentière and Grand Etrèt, respectively (Isotta et al., 2014), annual temperature variation amplitudes (dTa) of 8 °C and atmospheric temperature lapse rate of 0.0065 °C/m (as in Protin et al., 2019). We then performed a trial-and-error process to best fit the observed modern mass-balance data (Fig. S2). By simultaneously fitting the two mass-balance datasets, we obtained a melting positive degree day factor (mPDD) of 5.6 mm w. e. $°C^{-1}$ d$^{-1}$ that was adopted in the subsequent iSOSIA simulations.

Simulations were run using the ice-corrected bedrock DEM of the Dora Baltea catchment with a resolution of 60 m (Viani et al., 2020), and a critical snow accumulation slope (i.e. minimum slope for avalanching) of 0.7, which was found to best reproduce ice accumulation in cirque glaciers during preliminary tests. We ran simulations with fixed climate forcing until steady-state conditions (i.e. maximum one grid cell oscillations in ice extent), generally attained after 1 ka simulation time. In a first series of test simulations, we used the present-day precipitation grid from Isotta et al. (2014) (downscaled to 60-m resolution; Fig. S1A) for computing locally maximum snow accumulation rate, and varied sea-level temperatures from 14.1 to 14.5 °C (0.1 °C increment), with the aim to model the LIA glacier extents (Fig. S4). These initial simulations allowed to set a reference sea-level temperature for the LIA (Fig. S4) for comparison with subsequent simulations of older paleoglacial stages.

Based on Lateglacial paleo-temperature estimates from independent Alpine proxies (Heiri et al., 2014 and references therein), we then tested different ice simulations by varying sea-level temperatures from 11.4 to 14.4 °C (0.2 °C increment). We first computed
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** $^{10}$Be surface-exposure ages from Valsavarenche tributary catchment (location in Fig. 1; modified 2 m-resolution DEM from Regione Autonoma Valle d'Aosta; LIA glacier extent from GlaRiskAlp Project, http://www.glariskalp.eu). A) Sample locations and (re-)calculated $^{10}$Be surface-exposure ages of morainic boulders (yellow: present study; green: after Baroni et al., 2021) and polished bedrocks (red) in the valley head. Orange lines represent Lateglacial morainic ridges (after Baroni et al., 2021). Black box highlights the extent of the zoom presented in Fig. 4B, red line depicts the extent of the studied catchment. B) Zoom on the polished-bedrock altitudinal transects and on the latero-frontal moraine sampled in this study. C) Individual (dashed lines) and summed (continuous lines) KDEs of $^{10}$Be surface-exposure ages. The modes and uncertainties of the summed KDEs are reported in bold black font. Colour code (valid also for sample names and ages in panels A and B) highlights the two age clusters (red and blue labels, see main text for details). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

paleo-temperature scenarios for which simulated ice configurations best fitted paleoglacier front and thickness constraints for the glacial stages identified in all the three catchments. In a second experiment, we conducted ice simulations for only one glacial stage (the most recent; i.e., the Younger Dryas, cf. section 4.1) using the same sea-level temperature range but different precipitation inputs: (1) uniformly decreased precipitation across the three catchments by 25% and 50% compared to present-day conditions, or (2) spatially-variable precipitation pattern across the Dora Baltea drainage system, estimated from large-scale climate models for the Younger Dryas (PaleoClim data downscaled to 60-m resolution; Brown et al., 2018; Fig. S1B). This latter scenario was only tested for Valpelline and Valsavarenche catchments.

# 4. Results

## 4.1. $^{10}$Be surface-exposure dating

$^{10}$Be surface-exposure ages newly-acquired (Table 1) or recalculated (Table S2) are shown in Figs. 3–5 and S5. Two distinct age populations clustering around ca. 13 and 11 ka are observed in all the three tributary valleys, based on KDE distributions and coherent past glacier geometries (cf. section 4.2).

In Valpelline, the summed KDE mode of the first age cluster is at $13.0^{-0.5}_{+1.1}$ ka (red, Fig. 3C) and includes the furthest downstream polished-bedrock knob in the valley floor (VALP01, ~16 km from the

LIA ice limit; Fig. 3A), together with the two lowest samples from the left valley side transect (VALP04-05, 2344–2447 m a.s.l.; Fig. 3B) and the highest sample from the right valley side and most upstream transect (VALP07, 2524 m a.s.l.; Fig. 3B). The latter is also the oldest age of the cluster ($13.9 \pm 0.6$ ka) and is located just below the glacial trimline (Fig. 2A). The second age population has a summed KDE mode of $10.6^{-0.4}_{+1.2}$ ka (blue, Fig. 3C) and includes all valley-bottom samples upstream of VALP01 (VALP03, 12, 11, ~7–3 km from the LIA ice limit), the highest sample from the left valley side transect (VALP06, 2521 m) and the samples along the right valley side transect below VALP07 (VALP08-10, 2334 to 2177 m a.s.l.; Figs. 2A–3B). The oldest and highest sample of the cluster is VALP06 ($12.0 \pm 0.6$ ka).

In Valsavarenche, the distribution of $^{10}$Be surface-exposure ages (including new and literature data from Baroni et al., 2021) also indicates two populations with summed KDE modes of $13.1^{-0.6}_{+1.2}$ and $11.2^{-1.3}_{+0.6}$ ka (Fig. 4C). The oldest population includes the two polished-bedrock knobs sampled along the main valley floor (VSAV14-15; Fig. 4A), at ~6 and 3.5 km from the LIA ice limit, all but one (lowest) samples from the left-side altitudinal transect (VSAV01-03, 2297–2659 m a.s.l.; Figs. 2B and 4B), and all bedrock samples from the right valley side (VSAV05, 08, 12–13, 1997–2663 m; Fig. 4B). In addition, this cluster also includes two morainic boulders deposited by tributary glaciers expanding on the right side of Valsavarenche (GP34.18 and 31.18, Baroni et al., 2021, Fig. 4A). The oldest age of the cluster is obtained for a valley-floor
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** $^{10}$Be surface-exposure ages from upper Val di Cogne tributary catchment (location in Fig. 1; modified 2 m-resolution DEM from Regione Autonoma Valle d'Aosta; LIA glacier extent from GlaRiskAlp Project, http://www.glariskalp.eu). A) Sample locations and (re-)calculated $^{10}$Be surface-exposure ages of morainic boulders (green, after Baroni et al., 2021) and polished bedrocks (red, this study) in the upper valley. Orange lines represent Lateglacial morainic ridges (after Baroni et al., 2021). B) Individual (dashed lines) and summed (continuous lines) KDEs of $^{10}$Be surface-exposure ages. The modes and uncertainties of the summed KDE are reported in bold black font. Colour code (valid also for sample names and ages in panel A) highlights the two age clusters (red and blue labels, see main text for details). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

**Table 1**

$^{10}$Be surface-exposure dating for new samples collected in Valpelline (VALP), Valsavarensche (VSAV) and Val di Cogne (CONGE) tributary catchments. Sample locations, topographic shielding, average sample thickness, $^{10}$Be/$^{9}$Be blank corrected ratios, $^{10}$Be concentrations and exposure ages are reported. Sample density is assumed to be 2.65 g cm$^{-3}$ for all samples. Mass of quartz dissolved, mass of Be carrier and non-blank corrected $^{10}$Be/$^{9}$Be ratios are reported in Table S1. 1 Topographic shielding correction according to Dunne et al. (1999). 2 $^{10}$Be/$^{9}$Be ratios of batch-specific analytical blanks used for the correction are 5.4 ± 0.6 × 10$^{-15}$ (VALP samples) and 5.5 ± 1.0 × 10$^{-15}$ (VSAV and COGNE samples). 3 Ages are reported with external uncertainties (i.e. including both analytical errors and production-rate uncertainties). Ages were calculated with a SLHL $^{10}$Be production rate of 4.16 ± 0.10 at g$^{-1}$ yr$^{-1}$ (Claude et al., 2014) and the LSDn scaling scheme (Lifton et al., 2014) and considering an estimated surface erosion rate of 0.1 mm ka$^{-1}$ (consistent with Wirsig et al., 2016 and Serra et al., 2022). No snow-cover correction was applied.

<table><thead><tr><th>Sample Name</th><th>Location WGS 84 (°N/°E)</th><th>Elevation (m a.s.l.)</th><th>Topographic shielding<sup>1</sup></th><th>Sample thickness (cm)</th><th><sup>10</sup>Be/<sup>9</sup>Be blank corrected (10<sup>-14</sup>)<sup>2</sup></th><th><sup>10</sup>Be/<sup>9</sup>Be uncertainty (%)</th><th><sup>10</sup>Be concentration (10<sup>5</sup> at g<sup>-1</sup>)</th><th><sup>10</sup>Be exposure age (ka)<sup>3</sup></th></tr></thead><tbody><tr><td>VALP01</td><td>45.8499/7.3863</td><td>1471</td><td>0.959</td><td>4</td><td>14.3</td><td>3.30</td><td>1.89 ± 0.06</td><td>13.8 ± 0.5</td></tr><tr><td>VALP03</td><td>45.8952/7.4752</td><td>1834</td><td>0.930</td><td>2</td><td>6.5</td><td>5.79</td><td>1.85 ± 0.11</td><td>10.9 ± 0.7</td></tr><tr><td>VALP04</td><td>45.8973/7.5069</td><td>2344</td><td>0.915</td><td>2.5</td><td>24.3</td><td>3.35</td><td>3.11 ± 0.10</td><td>12.7 ± 0.5</td></tr><tr><td>VALP05</td><td>45.8956/7.5071</td><td>2447</td><td>0.946</td><td>2.5</td><td>26.4</td><td>3.39</td><td>3.50 ± 0.12</td><td>12.8 ± 0.5</td></tr><tr><td>VALP06</td><td>45.8941/7.5082</td><td>2521</td><td>0.935</td><td>6</td><td>25.9</td><td>3.98</td><td>3.33 ± 0.13</td><td>12.0 ± 0.6</td></tr><tr><td>VALP07</td><td>45.9217/7.5142</td><td>2524</td><td>0.950</td><td>2.5</td><td>27.1</td><td>3.31</td><td>4.06 ± 0.13</td><td>13.9 ± 0.6</td></tr><tr><td>VALP08</td><td>45.9214/7.5169</td><td>2334</td><td>0.851</td><td>2.5</td><td>17.2</td><td>3.34</td><td>2.31 ± 0.08</td><td>10.2 ± 0.4</td></tr><tr><td>VALP09</td><td>45.9209/7.5176</td><td>2270</td><td>0.848</td><td>2.5</td><td>13.3</td><td>3.31</td><td>2.25 ± 0.08</td><td>10.5 ± 0.4</td></tr><tr><td>VALP10</td><td>45.9197/7.5192</td><td>2177</td><td>0.959</td><td>2.5</td><td>18.9</td><td>3.31</td><td>2.42 ± 0.08</td><td>10.7 ± 0.4</td></tr><tr><td>VALP11</td><td>45.9171/7.5243</td><td>1982</td><td>0.955</td><td>2</td><td>14.8</td><td>3.23</td><td>2.07 ± 0.07</td><td>10.6 ± 0.4</td></tr><tr><td>VALP12</td><td>45.9092/7.5100</td><td>1990</td><td>0.798</td><td>2</td><td>5.11</td><td>5.3</td><td>1.92 ± 0.10</td><td>11.6 ± 0.7</td></tr><tr><td>VSAV01</td><td>45.5027/7.1958</td><td>2659</td><td>0.961</td><td>2.5</td><td>26.4</td><td>3.19</td><td>4.16 ± 0.13</td><td>13.0 ± 0.5</td></tr><tr><td>VSAV02</td><td>45.5036/7.2027</td><td>2446</td><td>0.968</td><td>2.5</td><td>23.7</td><td>3.17</td><td>3.60 ± 0.12</td><td>13.0 ± 0.5</td></tr><tr><td>VSAV03</td><td>45.5098/7.2021</td><td>2297</td><td>0.933</td><td>2</td><td>21.1</td><td>4.13</td><td>3.40 ± 0.14</td><td>14.1 ± 0.7</td></tr><tr><td>VSAV04</td><td>45.5113/7.2026</td><td>2212</td><td>0.817</td><td>2.5</td><td>15.8</td><td>3.39</td><td>2.37 ± 0.08</td><td>12.0 ± 0.5</td></tr><tr><td>VSAV05</td><td>45.5117/7.2063</td><td>1997</td><td>0.822</td><td>2</td><td>15.9</td><td>3.17</td><td>2.28 ± 0.07</td><td>13.5 ± 0.5</td></tr><tr><td>VSAV06</td><td>45.5140/7.2231</td><td>2638</td><td>0.987</td><td>2.5</td><td>22.8</td><td>3.40</td><td>3.65 ± 0.13</td><td>11.3 ± 0.4</td></tr><tr><td>VSAV07</td><td>45.5133/7.2258</td><td>2675</td><td>0.990</td><td>3</td><td>21.5</td><td>3.26</td><td>3.70 ± 0.12</td><td>11.1 ± 0.4</td></tr><tr><td>VSAV08</td><td>45.5191/7.2238</td><td>2663</td><td>0.984</td><td>2.5</td><td>26.7</td><td>3.23</td><td>4.53 ± 0.15</td><td>13.7 ± 0.5</td></tr><tr><td>VSAV09</td><td>45.5160/7.2212</td><td>2580</td><td>0.988</td><td>2.5</td><td>22.2</td><td>3.16</td><td>3.53 ± 0.11</td><td>11.4 ± 0.4</td></tr><tr><td>VSAV10</td><td>45.5156/7.2211</td><td>2578</td><td>0.992</td><td>2</td><td>22.9</td><td>3.37</td><td>3.76 ± 0.13</td><td>12.0 ± 0.5</td></tr><tr><td>VSAV11</td><td>45.5155/7.2207</td><td>2576</td><td>0.992</td><td>2</td><td>20.1</td><td>3.74</td><td>3.32 ± 0.13</td><td>10.7 ± 0.5</td></tr><tr><td>VSAV12</td><td>45.5142/7.2160</td><td>2400</td><td>0.952</td><td>2.5</td><td>21.0</td><td>3.15</td><td>3.36 ± 0.11</td><td>12.8 ± 0.5</td></tr><tr><td>VSAV13</td><td>45.5145/7.2106</td><td>2204</td><td>0.966</td><td>2.5</td><td>21.7</td><td>3.18</td><td>3.30 ± 0.11</td><td>14.2 ± 0.6</td></tr><tr><td>VSAV14</td><td>45.5145/7.2106</td><td>1958</td><td>0.957</td><td>2</td><td>19.3</td><td>6.35</td><td>2.80 ± 0.18</td><td>14.6 ± 1.0</td></tr><tr><td>VSAV15</td><td>45.5453/7.2119</td><td>1841</td><td>0.945</td><td>2</td><td>13.9</td><td>3.28</td><td>2.30 ± 0.08</td><td>13.3 ± 0.5</td></tr><tr><td>COGNE01</td><td>45.5684/7.3432</td><td>2468</td><td>0.983</td><td>2</td><td>23.7</td><td>8.84</td><td>3.80 ± 0.34</td><td>13.3 ± 1.2</td></tr><tr><td>COGNE02</td><td>45.5705/7.3428</td><td>2434</td><td>0.984</td><td>2</td><td>27.9</td><td>3.56</td><td>4.19 ± 0.15</td><td>14.9 ± 0.6</td></tr><tr><td>COGNE04</td><td>7.3428/7.3398</td><td>1919</td><td>0.930</td><td>2</td><td>14.8</td><td>4.65</td><td>2.26 ± 0.11</td><td>12.6 ± 0.7</td></tr></tbody></table>
bedrock knob at 14.6 ± 1.0 ka (VSAV14), while the highest sample (VSAV01) has an age of 13.0 ± 0.5 ka. The second and younger age cluster comprises the lowest bedrock sample from the left-side transect (VSAV04, 2012 m a.s.l.; Fig. 4A), five morainic boulders sampled in this study (VSAV06-11; Figs. 2C-4B) and additional five morainic boulders from Baroni et al. (2021) (GP24.18, 36.18, 37.18, 05.17 and 11.17; Fig. 4A and B). All sampled boulders are located on morainic ridges deposited by tributary glaciers nested on the right slopes of the Valsavarenche upper catchment, and which at that time reached 200–800 m downstream the LIA extent (Fig. 4). The oldest age of this second cluster, at 12.0 ± 0.5 ka, corresponds to both the lowest bedrock sample in the Valsavarenche main valley (VSAV04) and a morainic boulder sample from a right-side tributary (VSAV10).

Lastly, in Val di Cogne, all bedrock samples from the right-side altitudinal transect (COGNE01, 02, 04, 1919–2468 m a.s.l.; Fig. 2D–Fig. 5A), and one morainic boulder from a left-side tributary (GP16.18, 2350 m a.s.l.; Baroni et al., 2021, Fig. 5A), can be grouped in one age cluster around $12.8^{-0.5}_{+2.1}$ ka (Fig. 5B). The oldest age is at 14.9 ± 0.6 ka (COGNE02, 2434 m), while the highest sample (COGNE01, 2468 m, Fig. 2D) has an age of 13.3 ± 1.2 ka. A younger age cluster ($10.2^{-0.8}_{+1.7}$; Fig. 5B) consists of two morainic boulders in the valley bottom, around 5 km from the LIA ice limit (GP01.17, 02.17; Fig. 5A), and another morainic boulder from a left-side tributary (GP19.18, 2659 m a.s.l.; Fig. 5A). For this cluster, the oldest age is at 11.8 ± 0.5 ka from the valley-floor morainic boulder (GP01.17).

## 4.2. Paleoglacier simulations

From our calibrated iSOSIA simulations, we identified multiple climate scenarios (pairs of temperature and precipitation; Table 2) allowing to best match paleoglacier extent and thickness for the two identified glacial stages at ca. 13 and 11 ka, both recognized in all the three tributary valleys (Figs. 6–8). We acknowledge that, when polished bedrock exposure ages were used as ice front and thickness constraints due to the absence of morainic ridges (i.e. for the older stage in all investigated catchments, as well as for the younger stage in Valpelline), the obtained paleoglacier reconstructions provide only a minimum extent/thickness reached during a phase of glacier retreat. This implies that the obtained paleo-temperature and precipitation couples represent the paleo-climatic conditions at the transition between fullstadial (cold) and interstadial (warm) periods (cf. section 5.2).

In Table 2, for each paleo-precipitation input (i.e. present-day pattern, homogeneously-decreased precipitation, or PaleoClim data, Figs. S1A–B), the corresponding paleo-temperatures are indicated as temperature offsets ($\Delta T_{LIA}$) from the reference LIA sea-level temperature (14.5°; Fig. S4). For the old glacial stage (ca. 13 ka), simulations were conducted only with present-day precipitation, therefore only one set of precipitation rate and temperature is shown in Table 2.

<caption>**Table 2**</caption>

Temperature and precipitation inputs for paleoglacier simulations in the three studied tributary catchments. Temperature offsets are compared to the model-calibrated LIA temperature (Fig. S2, $\Delta T_{LIA}$ in the text). We report $\Delta T_{LIA}$ values for the YD/EH (young) glacial stage, and in brackets in italics $\Delta T_{LIA}$ values for the OD/BA (old) glacial stage (obtained only with same precipitation as today). Spatially-variable precipitation pattern across the Dora Baltea drainage system is based on PaleoClim data (Brown et al., 2018; Fig. S1B) and output results are reported for Valpelline and Valsavarenche catchments. See text for details and discussion.

<table><thead><tr><th rowspan="2">Ratio of present-day precipitation</th><th colspan="3">Best-fitting temperature anomaly (°C)</th></tr><tr><th>Valpelline</th><th>Valsavarenche</th><th>Val di Cogne</th></tr></thead><tbody><tr><td>1</td><td>-1.7 (-2.5)</td><td>-1.7 (-2.3)</td><td>-1.7 (-2.3)</td></tr><tr><td>0.75</td><td>-2.3</td><td>-2.3</td><td>-2.3</td></tr><tr><td>0.5</td><td>-3.1</td><td>-2.9</td><td>-3.1</td></tr><tr><td>Spatially-variable (PaleoClim)</td><td>-2.9</td><td>-2.1</td><td>/</td></tr></tbody></table>

When using the present-day precipitation pattern, the best-matching scenarios for the reconstructed ice extent/thickness were obtained by applying a $\Delta T_{LIA}$ of 2.5 °C (Valpelline, Figs. 6A) and 2.3 °C (Valsavarenche and Val di Cogne; Figs. 7A–8A) for the older glacial stage, and a $\Delta T_{LIA}$ of 1.7 °C for the younger glacial stage (same in all three valleys; Figs 6B, 7B and 8B). As illustrated by the longitudinal and cross-section profiles along the Valpelline catchment (Fig. 9), different temperature inputs are primarily reflected by changes in glacier extent, while glacier thickness is less sensitive to colder conditions (i.e. -2.3 to -2.7 °C; Fig. 9B). For the older glacial stage, ice reconstructions based on the longitudinal profile (matching our ice extent markers; Fig. 9A) are thus more easily constrained than based on ice-thickness differences from the cross-section profiles (Fig. 9B). While for the younger glacial stage, both ice-thickness (altitudinal transect) and ice-extent constraints allow to distinguish the best-matching paleo-temperature scenario.

For the younger glacial stage, different temperature offsets (Table 2) allow the best-matching between modelled and reconstructed ice extent/thickness, depending on applying homogeneously lowered present-day precipitation (regional decrease in precipitation by 25 and 50%) or imposing a different precipitation pattern (PaleoClim data; Brown et al., 2018). Similar temperature offsets for the three catchments were able to reproduce the reconstructed ice configurations when using a similar to present-day precipitation pattern decreased by a factor of 25% (-2.3 °C) and 50% (-2.9 to -3.1 °C; Table 2). Significantly different temperature drops were instead obtained between Valpelline and Valsavarenche catchments (~0.8 °C difference, Table 2) when applying the PaleoClim precipitation pattern. The PaleoClim dataset suggests highly-variable precipitation difference compared to present-day conditions within the three investigated catchments (Fig. S1B). In Valpelline, a 40–50% decrease in precipitation is observed mostly in the accumulation zone, while the precipitation decrease is more homogeneous across the Valsavarenche and Val di Cogne catchments. Among the best-matching climatic scenarios, obtained with the different precipitation patterns (Table 2), our paleoglacier simulations do not allow to discriminate between the different scenarios as the ice configurations do not vary significantly (Fig. S6).

Paleo-ELA estimates were inferred based on our simulations (Table 3). The reconstructed LIA paleo-ELAs are 3112 m for all the three valleys, resulting in ELA depressions ($\Delta$ELA) of 385 m (Valpelline) and 354 m (Valsavarenche and Val di Cogne) for the older glacial stage, and of 262 m in all three tributary valleys for the younger glacial stage.

# 5. Discussion

## 5.1. Lateglacial paleoglacial stages

Based on our $^{10}$Be surface-exposure age datasets (Figs. 3–5, Table 1) together with numerical ice simulations (Figs. 6–8), we propose the spatial reconstructions of two main Lateglacial paleoglacier stages at ca. 13 and 11 ka (Fig. S5). These two stages, present in all the three investigated tributary valleys of the main Dora Baltea catchment, are consistent with Lateglacial stadials identified across the European Alps (Ivy-Ochs, 2015), and also correlated with North Hemisphere climate oscillation records (Rasmussen et al., 2014).

The first recorded stage (old glacial stage at ca. 13 ka) is mainly derived from polished bedrock surfaces, and therefore rather reflects a phase of ice retreat and thinning which started at least ~14 ka ago, as suggested by $^{10}$Be exposure ages at the highest elevations
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Valpelline paleoglacier simulation results obtained with present-day precipitation and temperature drop of -2.5 °C (A) and -1.7° (B) compared to the model-calibrated LIA temperature. Red and blue dots represent polished-bedrock samples used to fit the modelled ice extent and thickness. Colour code distinguishes samples from the older (red) and younger (blue) glacial stages (sample ages and clusters in Fig. 3). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Valsavarenche paleoglacier simulation results obtained with present-day precipitation and temperature drop of -2.3° (A) and -1.7° (B) compared to the model-calibrated LIA temperature. Red and blue dots represent morainic boulders and polished-bedrock samples used to fit the modelled ice extent and thickness. Colour code distinguishes samples from the older (red) and younger (blue) glacial stages (sample ages and clusters in Fig. 4). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

or furthest downstream along the three investigated valleys. We relate this constrained timing to the beginning of the Bølling-Allerød interstadial (BA, 14.6–12.8 ka; Heiri et al., 2014), which indicates a phase of general warming following the Oldest Dryas cold period (OD, 19–14.6 ka; Ivy-Ochs et al., 2008). It is interesting to note that a significantly earlier onset of Lateglacial ice thinning, at ca. 18 ka, has been recorded in other locations of the High European Alps (Wirsig et al., 2016 and references therein; Lehmann et al., 2020). We cannot however certify whether the time

difference for our study area indicates an actual delay in ice thinning compared to the other sectors of the High Alps, reflecting differences in glacier sensitivity to post-LGM climate change, or if the delayed onset in glacier thinning is related to the locations of our sampling sites (i.e., upper catchment area and close to the trimline). Our numerical ice experiments (Fig. 9) also show that glacier variations in response to temperature changes are primarily reflected in glacier extent (Fig. 9A), while glacier thickness in the upper catchment is only weakly sensitive to different temperature
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Val di Cogne paleoglacier simulation results obtained with present-day precipitation and temperature drop of -2.3° (A) and -1.7° (B) compared to the model-calibrated LIA temperature. Red and blue dots represent morainic boulders and polished-bedrock samples used to fit the modelled ice extent and thickness. Colour code distinguishes samples from the older (red) and younger (blue) glacial stages (sample ages and clusters in Fig. 5). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Modelled Valpelline paleoglacier longitudinal (A) and cross-section (B) profiles obtained with present-day precipitation and different temperature offsets from the model-calibrated LIA conditions. The longitudinal profile has been taken along the modern Valpelline valley floor, and the cross-section profile has been selected at the location of the altitudinal transect of VALP07-11 (Fig. 3B). Red and blue dots represent polished-bedrock samples used to fit the modelled ice extent and thickness. The samples which are not located along the longitudinal profile (see Figs. 3–6 for sample locations) were projected (Fig. 9A) for better visualization. Colour code distinguishes samples from the older (red) and younger (blue) glacial stages (sample ages and clusters in Fig. 3). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

**Table 3**

Summary of the paleo-ELAs and ΔELA values (compared to the LIA) obtained in the three studied tributary catchments, based on numerical simulations with present-day precipitation. Three different ice stages have been investigated: LIA, OD/BA and YD/EH. Similar ELA and ΔELA values were obtained when fitting the same ice front/thickness constraints with different temperature/precipitation couples.

<table><thead><tr><th rowspan="2"></th><th colspan="3">Paleo-ELA (m a.s.l.)</th><th colspan="2">ΔELA (m)</th></tr><tr><th>LIA</th><th>OD/BA stage</th><th>YD/EH stage</th><th>OD/BA stage</th><th>YD/EH stage</th></tr></thead><tbody><tr><td>Valpelline</td><td>3112</td><td>2727</td><td>2850</td><td>385</td><td>262</td></tr><tr><td>Valsavarenche</td><td>3112</td><td>2758</td><td>2850</td><td>354</td><td>262</td></tr><tr><td>Val di Cogne</td><td>3112</td><td>2758</td><td>2850</td><td>354</td><td>262</td></tr></tbody></table>

scenarios (Fig. 9B). We therefore cannot exclude that an earlier onset in ice lowering may have been recorded more downstream in our studied valleys (Magrani et al., 2022) but has been delayed at

our sampling locations. Such delayed timing for ice thinning is also consistent with ¹⁰Be exposure ages from bedrock surfaces close to the trimline in upstream Alpine catchments (Kelly et al., 2006; Böhlert et al., 2011; Hippe et al., 2014; Wirsig et al., 2016), associated with the persistence of thick valley glacier at high elevations despite the general post-LGM retreat of glaciers from the piedmont areas. Further investigation of the above hypotheses (i.e. actual delay in post-LGM ice thinning vs. location of sampling sites) would require detailed sampling around the trimline at different catchment locations, which is beyond the scope of this study. Additional investigation could exclude potential late-exhumation (Wirsig et al., 2016) of high-elevation samples and/or late ice persistence of small tributary cirque glaciers (e.g. potentially for VALP04, 05 and COGNE02; Figs. 3 and 5).

In Valsavarenche and Val di Cogne, three erratic boulders on
moraine ridges also present ¹⁰Be exposure ages around 13 ka (GP31.18, 34.18, 16.18, Baroni et al., 2021, Figs. 4 and 5). Based on their elevation and orientation, we interpret these morainic ridges to have been built during a short cooling interval within the Bølling-Allerød interstadial (Older Dryas or Aegelsee Oscillation at 14.0–13.9 ka, or Gerzensee Oscillation at 13.3–13.0 ka; Lotter et al., 2012). These morainic ridges could therefore correspond to the Alpine Daun glacial stage (Ivy-Ochs, 2015), also recognized in the Mont Blanc sector of the Dora Baltea catchment (Serra et al., 2022). Based on ice simulation results (Fig. 8A), we also propose that the latero-frontal moraine along the Val di Cogne was possibly built during the Daun ice re-advance/stillstand, despite its apparent younger age based on boulder ¹⁰Be exposure dating (GP01.17, 02.17; Fig. 5). Our ice simulations could only match this moraine location when applying temperature/precipitation conditions in agreement with other ¹⁰Be ages of 13–14 ka in Val di Cogne (Fig. 8A) as well as in Valsavarenche and Valpelline. Although such observation may require further dating constraints in Val di Cogne, we propose that the apparent young ¹⁰Be exposure ages from the latero-frontal moraine along the Val di Cogne may result from late exhumation following post-depositional moraine degradation.

For the second glacial stage, all three tributaries suggest a common timing at ca. 11 ka that we relate to the Younger Dryas cold phase (YD, 12.9–11.7 ka; Heiri et al., 2014) and its transition into the early Holocene warming (EH, 11.7–8.0 ka; Heiri et al., 2014). In Valsavarenche and Val di Cogne (Figs. 4 and 5), this paleoglacial stage is mainly recorded from morainic boulders with ¹⁰Be surface-exposure ages of 11–12 ka (new samples VSAV06-11, and recalculated ages of GP11.17, 05.17, 19.18 from Baroni et al., 2021). We interpret this age cluster as documenting ice stillstand/re-advance in response to the YD, being described in the Alps as the Egesen stadial (Ivy-Ochs, 2015). Ice simulation results also suggest that the Valsavarenche moraine with younger ¹⁰Be surface-exposure ages of 9–10 ka (GP36.18, 37.18; recalculated from Baroni et al., 2021) may belong to the YD paleoglacier configuration (Fig. 7B). Thus, apparent younger ages might result from small ice fluctuations of similar extent occurring at the YD-early Holocene transition (Baroni et al., 2021; Protin et al., 2021) or from sample late exhumation. A YD re-advance of a southern tributary cirque glacier is also suggested in Valpelline by our ice simulation results, in agreement with the ¹⁰Be surface-exposure age of sample VALP06 around 12 ka. All the other polished bedrock samples with ¹⁰Be surface-exposure ages around 12–11 ka in Valsavarenche and Valpelline (Figs. 3 and 4) record the phase of ice withdrawal induced by rapid warming following the end of the YD. Post-YD ice decay was most probably abrupt, as indicated in Valpelline (Figs. 3 and 9) where polished bedrock samples collected over a longitudinal distance of ~5 km and an altitude range of ~350 m present similar ¹⁰Be surface-exposure ages of ca. 10–11 ka (Figs. 3 and 9).

In summary, our combined results from ¹⁰Be exposure dating and numerical simulations support two main distinct Lateglacial paleoglacier configurations (1) at the transition between the Oldest Dryas (OD) cold period (based on the polished bedrock sites) and the Bølling-Allerød (BA) interstadial (stage OD/BA hereafter; Figs. 6–8A), with potential minor influence of the Older Dryas or Gerzensee Oscillation short cooling interval (from three morainic samples), and (2) at the transition between the Younger Dryas (YD) cold period and the early Holocene (EH) warming (stage YD/EH hereafter; Figs. 6B–8B). In the following section, we investigate the potential climatic conditions (precipitation/temperature input scenarios in ice simulations) for these two glacial stages, in order to discuss glacier sensitivity to climate forcing in the three investigated valleys.

## 5.2. Paleoclimatic interpretation and glacier sensitivity

We combined our paleoglacier reconstructions from geomorphological mapping and ¹⁰Be surface-exposure dating with the iSOSIA ice-flow model to numerically investigate the variability in sensitivity of Alpine glaciers to climate forcing. The modelling approach applied in this study was chosen as best trade-off that allows for catchment-scale computation of paleoglacier extents using high-resolution bedrock topography (60-m resolution in our simulations). However, our modelling approach also needed some simplifications, including steady-state subglacial hydrology (e.g. Ugelvig et al., 2018) but without transient meltwater fluctuations or other factors that may accelerate basal sliding in transition periods and differ between catchments. Also, the complex physics controlling ice melting (i.e. solar radiation or heat transfer at the ice surface and within the ice) that has been rather empirically calibrated on modern mass-balance observations is lumped into one single PDD factor in our mass-balance model (Fig. S2). Moreover, we used the ice-corrected bedrock DEM of the Dora Baltea catchment (Viani et al., 2020) as glacier bed input, which may present up to 30% uncertainties in bedrock elevations for currently glaciated zones (Linsbauer et al., 2012) and may differ from Lateglacial topographic conditions (e.g. different valley sediment or lake infills). However, we consider that small changes in topographic conditions would only have a minor influence on the ice simulation outputs and would not change the first-order paleoclimatic interpretation from our paleoglacier reconstructions and simulations. Lastly, we had to assume glacier steady-state conditions to model paleoglacier configurations, which rather represent post-stadial ice retreat stages, given our chronology on mostly polished bedrocks. We therefore acknowledge that the obtained paleo-temperature and precipitation couples represent the paleoclimatic conditions at the transition between full stadial (OD and YD) and interstadial (BA and EH) periods. Hence, potentially colder or warmer paleo-temperatures could have been in place during the full stadials and interstadials, respectively.

The reconstructed ΔELAs (compared to the LIA, Table 3) for both the OD/BA and the YD/EH stages, based on numerical simulations with present-day precipitation pattern, are in good agreement with values obtained from previous paleoglacial studies in nearby Alpine sectors (Protin et al., 2019; Serra et al., 2022) and in the Gran Paradiso massif (Baroni et al., 2021). The ΔELA depressions of 354 m (Valsavarenche and Val di Cogne) and 385 m (Valpelline) for the OD/BA stage, fit within the range of ΔELA reported for the Daun glacial stage in the Mont Blanc sector of the Dora Baltea catchment (365–478 m; Serra et al., 2022). For the YD/EH stage, our ΔELA estimate of 262 m obtained in all three valleys is in agreement with previous ΔELA estimates for the Egesen stadial from Valsavarenche and Val di Cogne (~200 m; Baroni et al., 2021) and from the Argentière glacier on the north-western side of the Mont Blanc massif (100–300 m; Protin et al., 2019).

Output paleo-temperatures for the OD/BA and YD/EH ice stages, assuming present-day precipitation (Figs. 6–8), are analogous in Valpelline, Valsavarenche and Val di Cogne (Table 2) and are hereafter discussed in comparison to Alpine paleo-temperature reconstructions from other paleoclimate records. For the OD/BA transition, our temperature anomaly compared to the model-calibrated LIA temperature (ΔT<sub>LIA</sub>) is -2.3 °C for Valsavarenche and Val di Cogne and -2.5 °C for Valpelline. According to data collected in the HISTALP project (http://www.zamg.ac.at/histalp/; Auer et al., 2007), LIA temperature depression in the Alps compared to modern value are between -1 and -1.5 °C. Therefore, based on our LIA-calibrated simulations, output paleo-temperature anomaly compared to present-day (ΔT<sub>Present</sub>) for the OD/BA stage would range between -3.3 and -3.8 °C.
Chironomid assemblages in other Alpine regions record distinctly the OD cold period from the BA warming, with reconstructed OD/present-day difference in mean summer temperatures between -3° and -8 °C (Heiri and Millet, 2005; Laroque and Finsiger, 2008; Samartin et al., 2012). We can estimate an increase in seasonal amplitude of up to 2 °C compared to present-day, due to larger drop in paleo-winter temperature compared to summer temperature around the YD period, as suggested by pollen reconstructions in SW and SE Europe (Davis et al., 2003). As a consequence, our ΔT<sub>Present</sub> estimates of -3.3 to -3.8 °C would be then on the upper limit of the summer temperature decrease indicated from chironomids during the OD period.

For the YD/EH transition, we obtained a ΔT<sub>LIA</sub> of -1.7 °C for all the three studied catchments, corresponding to ΔT<sub>Present</sub> around -2.7 to -3.2 °C. As for the OD/BA transition, this is within the upper range of paleo-temperature reconstructions from chironomid Alpine records (between -2° and -6 °C summer temperature anomaly compared to present-day; Samartin et al., 2012; Ilyashuk et al., 2009; assuming up to 2 °C increase in seasonal amplitude; Davis et al., 2003) and pollen regional temperature reconstruction (around -3 °C of area-average mean annual temperature anomaly compared to today; Davis et al., 2003). When compared with paleo-temperature estimates from the nearby Mont Blanc massif (Argentière glacier; Protin et al., 2019), our ΔT<sub>LIA</sub> (-1.7 °C) differs from the estimates of Protin et al. (2019) who proposed ΔT<sub>LIA</sub> of -2.8 °C for the Argentière glacier. Although we used similar combined dating-modelling approach as Protin et al. (2019), our ice model and calibration approaches differ, which may explain the output differences in estimated paleo-temperatures. Our ΔT<sub>Present</sub> estimates are also lower than estimates from treeline records (around -5 °C; Heiri et al., 2014). However, it cannot be excluded that differences between our ΔT estimates and from other records could be due to spatial variability in temperature anomalies across the European Alps (Bartlein et al., 2011).

At last, only a small difference in ΔT<sub>LIA/present</sub> (i.e. ~0.8 °C) be-
tween the OD/BA and YD/EH periods is observed from our paleo-
glacier reconstructions, despite very different ice configurations
(Figs. 6–8). Such temperature difference between the OD/BA and
YD/EH periods is in overall agreement with data trends from the
chironomid records (Heiri and Millet, 2005; Laroque and Finsiger,
2008; Samartin et al., 2012), indicating that the rapid warming of
the BA interstadial resulted in ~3 °C temperature increase
compared to the OD cold period (Heiri et al., 2014), subsequently
followed by cooling of similar magnitude (1.5–3°C; Heiri et al.,
2014) at the end of the BA. This would have led to the re-
establishment of cold stadial conditions during the YD, as also
shown by data from Greenland ice core stratigraphy (Rasmussen
et al., 2014), but resulted in significantly different paleoglacier ex-
tents between OD and YD (Figs. 6–8).

Simulation results from the YD/EH stage using different pre-
cipitation scenarios do not support a spatially-variable precipita-
tion pattern during the YD across the Dora Baltea catchment, as
suggested from PaleoClim reconstructions (Fig. S1B). Indeed, when
imposing spatially-variable precipitation rates with significantly
wetter YD conditions in Valsavarenche than in Valpelline (Paleo-
Clim data; Brown et al., 2018; Fig. S1B), different paleo-temperature
anomalies are needed in order to fit to the YD paleoglacier con-
figurations in the two valleys (i.e. 0.8 °C difference; Table 2), which
appears unrealistic considering the spatial proximity of the two
tributary valleys. We therefore favour, from our combined paleo-
glacier reconstructions and ice simulations, YD paleoclimate sce-
narios with a similar-to-today precipitation pattern (i.e. wetter
condition in Valpelline than in Valsavarenche/Val di Cogne;
Fig. S1A), as supported by the good agreement in paleo-

temperature anomaly estimates between the three tributary
catchments when assuming modern or homogeneously-decreased
precipitation (Table 2). It is however difficult to quantitatively
assess which precipitation scenario (similar or decreased compared
to present-day) was potentially in place during the YD, since our
range of different paleoclimatic inputs (temperature-precipitation
couples; Table 2) provide similar ice extents/thicknesses, which all
satisfactorily fit the geomorphological constraints in the three
valleys, with only minor changes in glacier extent and thickness
(Fig. S6). Similarly, all the ΔT<sub>LIA/present</sub> based on uniformly-decreased
YD precipitation (Table 2) are within the range of estimated paleo-
temperatures inferred from independent Alpine paleoclimatic
proxies (between -2° and -6 °C, see discussion above). We how-
ever note that our ΔT<sub>LIA/present</sub> when assuming present-day pre-
cipitation lie in the upper range limit of paleo-temperature
anomalies compared to chironomids archives or other re-
constructions from paleoglaciers or global climate model outputs
such as PaleoClim (averaged temperature depression of around 5 °C
between present-day and the YD, within our study area), which
hence may argue in favour of scenarios with reduced precipitation
(and larger ΔT) compared to present-day.

While our simulation results suggest similar precipitation pattern compared to today across the Dora Baltea (central) catchment, they however do not provide evidence neither in favour nor against the hypothesis of a larger regional spatial re-organization in atmospheric circulation with south-westerly moisture advection over the European Alps during the YD, as suggested from previous studies using paleoclimate modelling and glacier ELAs-atmospheric temperature proxies compilation at the scale of the entire southern Europe (Rea et al., 2020 and references therein). Our investigated Dora Baltea tributary catchments are most certainly too spatially close to capture any potential signal of regional-scale change in atmospheric circulation. On the other hand, the observed larger ΔT<sub>present</sub> estimates for the Argentière glacier (north-western side of the Mont Blanc massif, ~40–50 km west/north-west of our study area; Protin et al., 2019), would support the idea of YD atmospheric re-organization with precipitation decrease along an East-West gradient as proposed from other paleoclimatic studies (Brown et al., 2018; Rea et al., 2020).

Finally, we observe a clear difference in glaciers' sensitivity to Lateglacial and Holocene climate variations between the three investigated catchments. When assuming present-day precipitation, the spatially-uniform change in paleoclimate between OD/BA and YD/EH (+0.8 °C difference between the two paleoglacial stages), and between YD/EH and LIA (+2 °C difference between YD and LIA), resulted in a much larger ice-front retreat along the main valley in Valpelline (around 8 and 7.5 km retreat, respectively; Fig. 6) than in Valsavarenche (around 4 and 3 km retreat, respectively; Fig. 7) and Val di Cogne (around 2 and 3 km retreat, respectively; Fig. 8). Such discrepancy in glacier response to uniform regional climate change can highlight different topographic and/or local climatic conditions between the three investigated catchments. Catchment hypsometry analyses for Valpelline and Valsavarenche indicate significantly larger catchment area above the paleo-ELA estimates in Valsavarenche compared to Valpelline (Fig. S7). As such, for a similar change around these paleo-ELA estimates, the different catchment hypsometric distributions would result in larger paleoglacier fluctuations in Valpelline compared to Valsavarenche. Additionally, the present-day precipitation distribution across the Dora Baltea catchment (Fig. S1A), shows moderately-wetter conditions over the northern tributaries (Valpelline) compared to the southern tributaries (Valsavarenche and Val di Cogne). This spatial difference in precipitation, if stable over the Lateglacial-Holocene period as suggested from our simulation outputs, had likely played a significant role in governing glaciers'
sensitivity across the Dora Baltea catchment. Indeed, (paleo)glaciers from the wetter Valpelline catchment (Fig. S1A) are more sensitive to climate change (i.e. higher potential to accumulate snow depending on temperature conditions), therefore being more responsive to change in (paleo-)ELAs than in the drier Valsavarenche/Val di Cogne catchments.

## 6. Conclusions

In this study, we combined paleoglacial reconstructions from $^{10}$Be surface-exposure dating of glacial landforms and deposits together with numerical glacier simulations in three tributary valleys of the Dora Baltea catchment (north-western Italian Alps). Our results allowed to provide quantitative constraints on the timing, ice configuration and potential paleoclimatic conditions for two Lateglacial ice stages in this Alpine area, while illustrating the benefits of the integrated dating-modelling approach.

In all the three investigated valleys (Valpelline, Valsavarenche and Val di Cogne), we constrained the onset of post-LGM ice thinning at ca. 14 ka, in apparent delay compared to other locations of the High Alps. We propose that this timing difference could be due to actual delay in ice thinning associated with spatial differences in glaciers' sensitivity to post-LGM climate change across the Alps, or alternatively to our spatial sampling strategy with specific targets in the upper catchment parts, where numerical ice simulations showed a limited ice-thickness response to early Lateglacial climate warming.

Our combined geomorphology-dating-modelling approach revealed in all the three sectors two distinct paleoglacier stages at ca. 13 and 11 ka, consistent with reported Lateglacial stadials identified across the European Alps and in correlation with North Hemisphere climate oscillations. The old ice stage (OD/BA) documents the ice configuration at the transition between the Oldest Dryas cold period (19–14.6 ka) and the Bølling-Allerød interstadial (14.6–12.8 ka), with potential minor influence of the Older Dryas or Aegelsee Oscillation short cooling intervals (14.0–13.9 ka or 13.3–13.0 ka, respectively). The young ice stage (YD/EH) represents paleoglacier oscillations at the transition between the Younger Dryas cold period (12.9–11.7 ka) and the early Holocene warming (11.7–8.0 ka). For both paleoglacial stages, we estimated paleo-ELAs in agreement with literature values in nearby alpine sectors, with $\Delta$ELA (compared to LIA conditions) of around 355–385 m for the OD/BA stage and 260 m for the YD/EH stage.

Paleoclimatic conditions associated with the reconstructed glacial stages were assessed based on the temperature/precipitation inputs used in iSOSIA ice-flow simulations to best fit ice surface and thickness constraints from our geomorphological observations. Our simulation outcomes suggest a similar-to-today precipitation pattern (i.e. same absolute values or homogeneously decreased across the investigated tributary valleys) over the Dora Baltea catchment for the two glacial stages, leading to analogous temperature anomalies compared to the Little Ice Age ($\Delta$T$_{LIA}$) in all the three valleys. $\Delta$T$_{LIA}$ anomalies obtained when assuming present-day precipitation (–2.3/-2.5 °C and –1.7 °C, for the OD/BA and YD/EH stages respectively) lie in the upper range of paleo-temperature reconstructions from other paleoclimatic proxies during the same Lateglacial periods. Although this may result from our steady-state assumption for paleoglacier reconstructions and modelling, this outcome may also argue in favour of scenarios with reduced precipitation compared to today (and thus larger associated $\Delta$T$_{LIA}$), but a priori not for a spatial change in precipitation pattern over the Dora Baltea area. However, our results do not exclude changes in precipitation pattern at a larger regional scale during the YD period, as suggested by previous studies.

Finally, we observed a clear difference in glaciers' sensitivity to

Lateglacial and Holocene climate variations between the three investigated catchments, with much larger ice-front retreats in Valpelline than in Valsavarenche and Val di Cogne, in response to a spatially-uniform change in paleoclimatic conditions. We explain such discrepancy as related to different topographic (hypsometry) and possibly climatic conditions (precipitation distribution) between the three investigated catchments: glaciers from the relatively lower-elevation and wetter Valpelline catchment had been more responsive to change in (paleo-)ELAs than glaciers in the higher-elevation and drier catchments of Valsavarenche/Val di Cogne. Our study hence highlights the variability in alpine glaciers' response to common climate forcing due to local factors, which are critical to take into account in order to make reliable paleoclimate inferences from paleoglacier records.

## Author contributions

ES, FM, PGV and NG designed the study. ES, FM, PGV and NG performed field investigations and sample collection. ES performed $^{10}$Be cosmogenic analysis with supervision of JC. FM and ES performed numerical glacier simulations with supervision of PGV and DE. ES and FM wrote the manuscript with input from all co-authors.

## Declaration of competing interest

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

## Data availability

Data will be made available on request.

## Acknowledgements

The authors warmly thank A. J. Déjanaz for the help during fieldwork. J. Gajic is thanked for support for $^{10}$Be surface-exposure dating sample preparation. The authors warmly thank the ASTER team for the measurements performed at the ASTER AMS facility (CEREGE, Aix-en- Provence). A. Ribolini, one anonymous reviewer and the editor C. O'Cofaigh are warmly thanked for constructive reviews and suggestions which improved the quality of our manuscript. This study was supported by the Swiss National Science Foundation SNSF (Grant PP00P2_170559) and the French ANR-PIA programme (ANR-18-MPGA-0006). All data are available in tables in the main text and online supporting information.

## Appendix B. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.quascirev.2022.107822.

## References

* Arnold, M., Merchel, S., Bourlès, D.L., Braucher, R., Benedetti, L., Finkel, R.C., Aumâtre, G., Gottdang, A., Klein, M., 2010. The French accelerator mass spectrometry facility ASTER: improved performance and developments. Nucl. Ins. Trum. Methods Phys. Res. Sect. B Beam Interact. Mater. Atoms 268, 1954–1959. https://doi.org/10.1016/j.nimb.2010.02.107.

* Auer, I., Böhm, R., Jurkovic, A., Lipa, W., Orlik, A., Potzmann, R., Schöner, W., Ungersböck, M., Matulla, C., Briffa, K., Jones, P., 2007. HISTALP—historical instrumental climatological surface time series of the Greater Alpine Region. Int. J. Climatol. 27 (1), 17–46. https://doi.org/10.1002/joc.1377.

* Bakke, J., Lie, Ø., Heegaard, E., Dokken, T., Haug, G.H., Birks, H.H., Dulski, P., Nilsen, T., 2009. Rapid oceanic and atmospheric changes during the Younger Dryas cold period. Nat. Geosci. 2 (3), 202–205. https://doi.org/10.1038/NGEO439.

* Baroni, C., Gennaro, S., Cristina, M., Ivy-ochs, S., Christl, M., Cerrato, R., Orombelli, G., 2021. Last Lateglacial glacier advance in the Gran Paradiso Group reveals
relatively drier climatic conditions established in the Western Alps since at least the Younger Dryas. Quat. Sci. Rev. 255, 106815. https://doi.org/10.1016/j.quascirev.2021.106815.

Bartlein, P.J., Harrison, S.P., Brewer, S., Connor, S., Davis, B.A.S., Gajewski, K., Guiot, J., Harrison-Prentice, T.I., Henderson, A., Peyron, O., Prentice, I.C., 2011. Pollen-based continental climate reconstructions at 6 and 21 ka: a global synthesis. Clim. Dynam. 37 (3), 775-802. https://doi.org/10.1007/s00382-010-0904-1.

Becker, P., Seguinot, J., Jouvet, G., Funk, M., 2016. Last Glacial Maximum precipita-
tion pattern in the Alps inferred from glacier modelling. Geograph. Helv. 71 (3),
173-187. https://doi.org/10.5194/gh-71-173-2016.

Braucher, R., Guillou, V., Bourlès, D.L., Arnold, M., Aumaître, G., Keddadouche, K., Nottoli, E., 2015. Preparation of ASTER in-house 10Be/9Be standard solutions. Nucl. Instrum. Methods Phys. Res. Sect. B Beam Interact. Mater. Atoms 361, 335-340. https://doi.org/10.1016/j.nimb.2015.06.012.

Brown, E.T., Edmond, J.M., Raisbeck, G.M., Yiou, F., Kurz, M.D., Brook, E.J., 1991. Examination of surface exposure ages of Antarctic moraines using in situ produced 10Be and 26Al. Geochem. Cosmochim. Acta 55, 2269-2283. https://doi.org/10.1016/0016-7037(91)90103-C.

Brown, J.L., Hill, D.J., Dolan, A.M., Carnaval, A.C., Haywood, A.M., 2018. PaleoClim, high spatial resolution paleoclimate surfaces for global land areas. Nature - Scientific Data 5, 180254. https://doi.org/10.1038/sdata.2018.254.

Böhlert, R., Egli, M., Maisch, M., Brandova, D., Ivy-Ochs, S., Kubik, P.W., Haeberli, W., 2011. Application of a combination of dating techniques to reconstruct the Lateglacial and early Holocene landscape history of the Alba region (eastern Switzerland), Geomorphology 127, 1-13. https://doi.org/10.1016/j.geomorph.2010.10.034.

Clark, P.U., Dyke, A.S., Shakun, J.D., Carlson, A.E., Clark, J., Wohlfarth, B., Mitrovica, J.X., Hostetler, S.W., McCabe, A.M., 2009. The last glacial maximum. Science 325, 710-714. https://doi.org/10.1126/science.1172873.

Claude, A., Ivy-Ochs, S., Kober, F., Antognini, M., Salcher, B., Kubik, P.W., 2014. The Chronicon landslide (Valle Leventina, southern Swiss Alps): age and evolution. Swiss J. Geosci. 107, 273-291. https://doi.org/10.1007/s00015-014-0170-z.

Davis, B.A., Brewer, S., Stevenson, A.C., Guiot, J., 2003. The temperature of Europe during the Holocene reconstructed from pollen data. Quat. Sci. Rev. 22 (15-17), 1701-1716. https://doi.org/10.1016/S0277-3791(03)00173-2.

Davis, P.T., Menounos, B., Osborn, G., 2009. Holocene and latest Pleistocene alpine glacier fluctuations: a global perspective. Quat. Sci. Rev. 28, 2021-2033. https://doi.org/10.1016/j.quascirev.2009.05.020.

Dunne, J., Elmore, D., Muzikar, P., 1999. Scaling factors for the rates of production of cosmogenic nuclides for geometric shielding and attenuation at depth on sloped surfaces. Geomorphology 27, 3-11. https://doi.org/10.1016/S0169-555X(98)00086-5.

Egholm, D.L., Knudsen, M.F., Clark, C.D., Lesemann, J.E., 2011. Modeling the flow of glaciers in steep terrains: the integrated second-order shallow ice approximation (iOSIA). J. Geophys. Res. Earth Surf. 116, 1-16. https://doi.org/10.1029/2010JF001900.

Egholm, D.L., Pedersen, V.K., Knudsen, M.F., Larsen, N.K., 2012. Coupling the flow of ice, water, and sediment in a glacial landscape evolution model. Geomorphology 141-142, 47-66. https://doi.org/10.1016/j.quascirev.2011.12.019.

Federici, P.R., Ribolini, A., Spagnolo, M., 2016. Glacial history of the maritime Alps from the last glacial maximum to the Little ice age. Geol. Soc. London, Spec. Publ. 433, 137-159. https://doi.org/10.1144/sp433.9.

Florineth, D., Schlüchter, C., 2000. Alpine evidence for atmospheric circulation patterns in Europe during the last glacial maximum. Quat. Res. 54 (3), 295-308. https://doi.org/10.1006/qres.2000.2169.

Gianotti, F., Forno, M.G., Ivy-Ochs, S., Monegato, G., Pini, R., Ravazzi, C., 2015. Stratigraphy of the ivrea morainic amphitheatre (NW Italy): an updated synthesis. Alp. Mediterr. Quat. 28, 29-58.

Heiri, O., Koinig, K.A., Spötl, C., Barrett, S., Brauer, A., Drescher-Schneider, R., Gaar, D., Ivy-Ochs, S., Kerschner, H., Luetscher, M., Moran, A., Nicolussi, K., Preusser, F., Schmidt, R., Schoeneich, P., Schwörer, C., Spraße, T., Terhorst, B., Tinner, W., 2014. Palaeoclimate records 60-8 ka in the Austrian and Swiss Alps and their forelands. Quat. Sci. Rev. 106, 186-205. https://doi.org/10.1016/j.quascirev.2014.05.021.

Heiri, O., Millet, L., 2005. Reconstruction of late glacial summer temperatures from chironomid assemblages in lac lautrey (jura, France). J. Quat. Sci. 20 (1), 33-44. https://doi.org/10.1002/jqs.895.

Hippe, K., Ivy-Ochs, S., Kober, F., Zasadni, J., Wieler, R., Wacker, L., Kubik, P.W., Schlüchter, C., 2014. Chronology of Lateglacial ice flow reorganization and deglaciation in the Gotthard Pass area, Central Swiss Alps, based on cosmogenic 10Be and in situ 14C. Quat. Geochronol. 19, 14-26. https://doi.org/10.1016/j.quageo.2013.03.003.

Ilyashuk, B., Gobet, E., Heiri, O., Lotter, A.F., van Leeuwen, J.F., van der Knaap, W.O., Ilyashuk, E., Oberli, F., Ammann, B., 2009. Lateglacial environmental and climatic changes at the Maloja Pass, Central Swiss Alps, as recorded by chronomids and pollen. Quat. Sci. Rev. 28 (13), 1340-1353. https://doi.org/10.1016/j.quascirev.2009.01.007.

I, F.A., Frei, C., Weilguni, V., Perțec Tadić, M., Lassegues, P., Rudolf, B., Pavan, V., Cacciamani, C., Antolini, G., Ratto, S.M., Munari, M., 2014. The climate of daily precipitation in the Alps: development and analysis of a high-resolution grid dataset from pan-Alpine rain-gauge data. Int. J. Climatol. 34 (5), 1657-1675. https://doi.org/10.1002/joc.3794.

Ivy-Ochs, S., 2015. Glacial variations in the European Alps at the end of the last deglaciation. Cuadernos Invest. Geogr. 41, 295. https://doi.org/10.18172/cig.2750.

Ivy-Ochs, S., Kerschner, H., Reuther, A., Preusser, F., Heine, K., Maisch, M.,

Kubik, P.W., Schlüchter, C., 2008. Chronology of the last glacial cycle in the European Alps. J. Quat. Sci.: Published for the Quaternary Research Association 23 (6), 559-573. https://doi.org/10.1002/jqs.1202.

Kelly, M.A., Ivy-Ochs, S., Kubik, P., von Blanckenburg, F., Schlüchter, C., 2006. Chronology of deglaciation based on 10Be dates of glacial erosional features in the Grimsel Pass region, central Swiss Alps. Boreas 35. https://doi.org/10.1080/03009480600690829, 634e643.

Kerschner, H., Ivy-Ochs, S., 2008. Palaeocimate from glaciers: examples from the eastern Alps during the alpine lateglacial and early Holocene. Global Planet. Change 60, 58-71. https://doi.org/10.1016/j.gloplacha.2006.07.034.

Kohl, C., Nishiizumi, K., 1992. Chemical isolation of quartz for measurement of in-situ -produced cosmogenic nuclides. Geochem. Cosmochim. Acta 56, 3583-3587. https://doi.org/10.1016/0016-7037(92)90401-4.

Kuhlemann, J., Rohling, E.J., Krumrei, I., Kubik, P., Ivy-Ochs, S., Kucera, M., 2008. Regional synthesis of mediterranean atmospheric circulation during the last glacial maximum. Science 321 (5894), 1338-1340. https://doi:10.1126/science.1157638.

Larocque, I., Finsinger, W., 2008. Late-glacial chironomid-based temperature re-
constructions for Lago Piccolo di Avigliana in the southwestern Alps (Italy).
Palaeogeogr. Palaeoclimatol. Palaeoecol. 257, 207-223. https://doi.org/10.1016/
j.palaeo.2007.10.021.

Lehmann, B., Herman, F., Valla, P.G., King, G.E., Biswas, R.H., Ivy-Ochs, S., Steinemann, O., Christl, M., 2020. Postglacial erosion of bedrock surfaces and deglaciation timing: new insights from the Mont Blanc massif (western Alps). Geology 48 (2), 139-144.

Lifton, N., 2016. Implications of two Holocene time-dependent geomagnetic models for cosmogenic nuclide production rate scaling. Earth Planet Sci. Lett. 433, 257-268. https://doi.org/10.1016/j.EPSL.2015.11.006.

Lifton, N., Sato, T., Dunai, T.J., 2014. Scaling in situ cosmogenic nuclide production rates using analytical approximations to atmospheric cosmic-ray fluxes. Earth Planet Sci. Lett. 386, 149-160. https://doi.org/10.1016/j.EPSL.2013.10.052.

Linsbauer, A., Paul, F., Haeberli, W., 2012. Modeling glacier thickness distribution and bed topography over entire mountain ranges with GlabTop: application of a fast and robust approach. J. Geophys. Res. 117, F03007. https://doi:10.1029/2011JF002313.

Lotter, A.F., Heiri, O., Brooks, S., van Leeuwen, J.F., Eicher, U., Ammann, B., 2012. Rapid summer temperature changes during Termination 1a: high-resolution multi-proxy climate reconstructions from Gerzensee (Switzerland). Quat. Sci. Rev. 36, 103-113. https://doi.org/10.1016/j.quascirev.2010.06.022.

Lowell, T.V., 1995. The application of radiocarbon age estimates to the dating of glacial sequences: an example from the Miami sublobe, Ohio. U.S.A. Quat. Sci. Rev. 14, 85-99. https://doi.org/10.1016/0277-3791(94)00113-.

Ludwig, P., Schaffernicht, E.J., Shao, Y., Pinto, J.G., 2016. Regional atmospheric cir-
culation over Europe during the last glacial maximum and its links to precip-
itation. J. Geophys. Res. Atmos. 121 (5), 2130-2145. https://doi.org/10.1002/
2015JD024444.

Luetscher, M., Boch, S., Sodemann, H., Spötl, C., Cheng, H., Edwards, R.L., Frisia, S., Hof, F., Müller, W., 2015. North Atlantic storm track changes during the last glacial maximum recorded by alpine speleothems. Nat. Commun. 6 (1), 1-6. https://doi:10.1038/ncomms734.

Magrani, F.G., Valla, P.V., Egholm, D.L., 2022. Modelling Alpine glacier geometry and subglacial erosion patterns in response to contrasting climatic forcing. Earth Surf. Process. Landforms 47, 1054-1072. https://doi.org/10.1002/esp.5302.

Martin, L.C.P., Blard, P.-H., Balco, G., Lavé, J., Delunel, R., Lifton, N., Laurent, V., 2017. The CREp program and the ICE-D production rate calibration database: a fully parameterizable and updated online tool to compute cosmic-ray exposure ages. Quat. Geochronol. 38, 25-49. https://doi.org/10.1016/j.QUAGEO.2016.11.006.

Merchel, S., Herpers, U., 1999. An update on radiochemical separation techniques for the determination of long-lived radionuclides via accelerator mass spectrometry. Radiochim. Acta 84. https://doi.org/10.1524/ract.1999.84.4.215.

Merz, N., Raible, C.C., Woollings, P., 2015. North Atlantic eddy-driven jet in interglacial and glacial winter climates. J. Clim. 28 (10), 3977-3997. https://doi.org/10.1175/JCLI-D-14-00255.1.

Monegato, G., Scardia, G., Hajdas, I., Rizzini, F., Piccin, A., 2017. The Alpine LGM in the boreal ice-sheets game. Sci. Rep. 7, 1-8. https://doi.org/10.1038/s41598-017-02148-7.

NOAA National Centers for Environmental Information, 2016. Monthly Global Climate Report for Annual 2016 published online January 2017, retrieved on November 3, 2022 from. https://www.ncei.noaa.gov/access/monitoring/monthly-report/global/201613.

Oerlemans, J., 2005. Atmospheric science: extracting a climate signal from 169
glacier records. Science 308, 675-677. https://doi.org/10.1126/science.1107046.

Protin, M., Schimmelpfennig, I., Mugnier, J.L., Buonchristiani, J.F., Le Roy, M., Pohl, B., Moreau, L., , ASTER Team, Aumaître, G., Bourlès, D., Keddadouche, K., 2021. Millennial-scale deglaciation across the European Alps at the transition between the Younger Dryas and the Early Holocene–evidence from a new cosmogenic nuclide chronology. Boreas 50, 671–685. https://doi.org/10.1111/bor.12519. ISSN 0300-9483.

Protin, M., Schimmelpfennig, I., Mugnier, J.L., Ravanel, L., Le Roy, M., Deline, P., Favier, V., Buonchristiani, J.F., Aumaître, G., Bourlès, D.L., Keddadouche, K., 2019. Climatic reconstruction for the younger dryas/early Holocene transition and the Little ice age based on paleo-extents of Argentière glacier (French Alps). Quat. Sci. Rev. 221, 105863. https://doi.org/10.1016/j.quascirev.2019.105863.

Rasmussen, S.O., Bigler, M., Blockley, S.P., Blunier, T., Buchardt, S.L., Clausen, H.B., Cvijanovic, I., Dahl-Jensen, D., Johnsen, S.J., Fischer, H., Gkinis, V., 2014.
A stratigraphic framework for abrupt climatic changes during the Last Glacial period based on three synchronized Greenland ice-core records: refining and extending the INTIMATE event stratigraphy. Quat. Sci. Rev. 106, 14–28. https://doi.org/10.1016/j.quascirev.2014.09.007.

Rea, B.R., Pellitero, R., Spagnolo, M., Hughes, P., Ivy-Ochs, S., Renssen, H., Ribolini, A., Bakke, J., Lukas, S., Braithwaite, R.J., 2020. Atmospheric circulation over Europe during the younger Dryas. Sci. Adv. 6 (50), 4844. https://doi.org/10.1126/sciadv.aba4844.

Samartin, S., Heiri, O., Vescovi, E., Brooks, S.J., Tinner, W., 2012. Lateglacial and early Holocene summer temperatures in the southern Swiss Alps reconstructed using fossil chironomids. J. Quat. Sci. 27 (3), 279–289. https://doi.org/10.1002/jqs.1542.

Scherler, D., Egholm, D.L., 2020. Production and transport of supraglacial debris: insights from cosmogenic $^{10}$Be and numerical modeling, chhota shigri glacier, Indian himalaya. J. Geophys. Res. Earth Surf. 125, 1–26. https://doi.org/10.1029/2020JF005586.

Serra, E., Valla, P.G., Gribenski, N., Carcaillet, Deline, P., 2022. Post-LGM glacial and geomorphic evolution of the Dora Baltea valley (western Italian Alps). Quat. Sci. Rev. 282, 107446. https://doi.org/10.1016/j.quascirev.2022.107446.

Stanford, J.D., Rohling, E.J., Bacon, S., Roberts, A.P., Grousset, F.E., Bolshaw, M., 2011. A new concept for the paleoceanographic evolution of Heinrich event 1 in the

North Atlantic. Quat. Sci. Rev. 30, 1047–1066. https://doi.org/10.1016/j.quascirev.2011.02.003.

Ugelvig, S.V., Egholm, D.L., Anderson, R.S., Iverson, N.R., 2018. Glacial erosion driven by variations in meltwater drainage. J. Geophys. Res. Earth Surf. 123, 2863–2877. https://doi.org/10.1029/2018JF004.

Uppala, S.M., Källberg, P.W., Simmons, A.J., Andrae, U., da Costa Bechtold, V., Fiorino, M., Gibson, J.K., Haseler, J., Hernandez, A., Kelly, G.A., Li, X., Onogi, K., Saarinen, S., Sokka, N., Allan, R.P., Andersson, E., Arpe, K., Balmaseda, M.A., Beljaars, A.C.M., van de Berg, L., Bidlot, J., Bormann, N., Caires, S., Chevallier, F., Dethof, A., Dragosavac, M., Fisher, M., Fuentes, M., Hagemann, S., Hólm, E., Hoskins, B.J., Isaksen, L., Janssen, P.A.E.M., Jenne, R., McNally, A.P., Mahfouf, J.F., Morcrette, J.J., Rayner, N.A., Saunders, R.W., Simon, P., Sterl, A., Trenberth, K.E., Untch, A., Vasiljevic, D., Viterbo, P., Woollen, J., 2005. The ERA-40 re-analysis. Q. J. R. Meteorol. Soc. 131, 2961–3012. https://doi.org/10.1256/qj.04.

Viani, C., Machguth, H., Huggel, C., Godio, A., Franco, D., Perotti, L., Giardino, M., 2020. Potential future lakes from continued glacier shrinkage in the aosta valley region (western Alps, Italy). Geomorphology 355, 107068. https://doi.org/10.1016/j geomorph.2020.107068.

Wirsig, C., Zasadni, J., Christl, M., Akçar, N., Ivy-Ochs, S., 2016. Dating the onset of LGM ice surface lowering in the High Alps. Quat. Sci. Rev. 143, 37–50. https://doi.org/10.1016/j.quascirev.2016.05.00.