Cosmogenic ³He paleothermometry on post-LGM glacial bedrock within the central European Alps

Diffusion properties of cosmogenic ³He in quartz at Earth surface temperatures offer the potential to directly reconstruct the evolution of past in situ temperatures from formerly glaciated areas, which is important information for improving our understanding of glacier–climate interactions. In this study, we apply cosmogenic ³He paleothermometry to rock surfaces gradually exposed from the Last Glacial Maximum (LGM) to the Holocene period along two deglaciation profiles in the European Alps (Mont Blanc and Aar massifs). Laboratory experiments conducted on one representative sample per site indicate significant differences in ³He diffusion kinetics between the two sites, with quasi-linear Arrhenius behavior observed in quartz from the Mont Blanc site and complex Arrhenius behavior observed in quartz from the Aar site, which we interpret to indicate the presence of multiple diffusion domains (MDD). Assuming the same diffusion kinetics apply to all quartz samples along each profile, forward model simulations indicate that the cosmogenic ³He abundance in all the investigated samples should be at equilibrium with present-day temperature conditions. However, measured cosmogenic ³He concentrations in samples exposed since before the Holocene indicate an apparent ³He thermal signal significantly colder than today. This observed ³He thermal signal cannot be explained with a realistic post-LGM mean annual temperature evolution in the European Alps at the study sites. One hypothesis is that the diffusion kinetics and MDD model applied may not provide sufficiently accurate, quantitative paleo-temperature estimates in these samples; thus, while a pre-Holocene ³He thermal signal is indeed preserved in the quartz, the helium diffusivity would be lower at Alpine surface temperatures than our diffusion models predict. Alternatively, if the modeled helium diffusion kinetics is accurate, the observed ³He abundances may reflect a complex geomorphic and/or paleoclimatic evolution, with much more recent ground temperature changes associated with the degradation of alpine permafrost.

1 Introduction

This study applies cosmogenic noble gas paleothermometry (Tremblay et al., 2014a) to attempt to reconstruct temperature changes associated with gradual ice lowering following the Last Glacial Maximum (LGM, ca. 27–19 ka; Clark et al., 2009) in two sites of the high European Alps. Because glaciers are sensitive to both temperature and precipitation, obtaining information about in situ temperature conditions from an independent proxy is critical to disentangling the role of either variable in recorded glacier fluctuations and to adequately use these records for paleoclimate reconstructions. In particular, paleoglacier records can be used as direct, site-specific paleo-precipitation indicators (e.g., Kerschner et al., 2000; Kerschner and Ivy-Ochs, 2008; Martin et al., 2020) to trace changes in regional atmospheric circulation systems (Kuhlemann et al., 2008; Becker et al., 2016; Gribenski et al., 2021). More detailed information about paleoclimate conditions would moreover improve both our understanding of glacier responses to current climate change as well as our ability to anticipate glacier evolutions for proposed future climate scenarios (Zemp et al., 2006; Haeberli et al., 2022). Furthermore, direct temperature constraints associated with paleoglacier variations are also critical to our understanding of glacier erosion processes (Hallet, 1979), which have profoundly shaped high-latitude and mountain landscapes over 10³ to 10⁶-year timescales (Herman et al., 2021) and which seem to relate to climatic conditions, among other factors (Koppes et al., 2015; Cook et al., 2020).

Available data on the relationship between glacier geometry and climate, as well as between glacial erosion and climate, are largely biased toward present-day and historical time periods, therefore obliging us to rely on the assumption that modern to centennial records are representative of the range of variation and mechanistic trends between climate–glacier variation and erosion operating on geological time scales (Jaeger and Koppes, 2016). While combined records of paleoglacier geometry and erosion rates on Late Pleistocene timescales are growing due to the recent development of analytical and numerical techniques (e.g., Kapannusch et al., 2020; Mariotti et al., 2021), obtaining direct quantitative paleoclimate constraints from formerly glaciated areas remains challenging, even for regions with relatively well-known paleoglacial histories. In the European Alps, the most detailed paleoglacier record goes back to the Late Pleistocene ice maximum advance, dated around ∼ 26–24 ka in the northern and central Alps (Monegato et al., 2017), in line with the global LGM. During the LGM, ice spread to within several tens of kilometers of the piedmonts and reached more than 1000–1500 m thickness in the main valleys (Ivy-Ochs, 2015; Wirsig et al., 2016a; Serra et al., 2022). More restricted stages (i.e., Gschnitz, Daun, Egesen stadials; Ivy-Ochs, 2015) marking the gradual retreat (and thinning) of the ice into the upper catchments followed between the LGM and the Younger Dryas cooling event (YD, 12.8–11.7 ka; Heiri et al., 2014a). During the early Holocene (i.e., the last 11 ka; Heiri et al., 2014a), glaciers retreated quickly behind the position where the Little Ice Age moraines are located today and remained within these limits for the rest of the Holocene period (Heiri et al., 2014a).

The timing and pattern of paleoglacier variations in the European Alps are consistent with polar ice oxygen isotope (δ¹⁸O) records from the Northern Hemisphere (North Greenland Ice Core Project members, 2004), which indicate that maximum Greenland temperature anomalies of around −20 ^∘C were reached at 25–20 ka, followed by a gradual warming until ca. 10 ka, with the last pronounced isotopic excursion marked by a −15 ^∘C temperature anomaly occurring at ca. 12 ka in association with the YD event (Buizert et al., 2018). After the YD, temperatures stabilized around values similar to today, with only minor fluctuations (less than 2 ^∘C) throughout the remaining Holocene period (Buizert et al., 2018). High-resolution δ¹⁸O in Alpine speleothems similarly support a coupling between the northern European Alps and Greenland records (Moseley et al., 2020; Li et al., 2021).

While there is evidence for a temporal coupling, a direct, scaled translation of polar ice records over the Alps to obtain quantitative temperature and precipitation constraints is not valid. Indeed, major climate forcing components, such as ice sheet extent, atmospheric greenhouse gas concentrations, and changes in ocean circulation, also underwent large-scale changes between the LGM and the Holocene transition (Clark et al., 2012). This resulted in variable atmospheric circulation patterns (Eynaud et al., 2009) and variable latitudinal temperature gradients (Heiri et al., 2014b) in the Northern Hemisphere during this period. Existing past climate information from Alpine paleoenvironmental proxies is mainly qualitative, with only a few scarce and fragmented quantitative temperature and precipitation records available for the pre-Holocene period (Heiri et al., 2014a). These are mostly from pollen and chironomid proxy records located on the outer rim of the Alpine range from lake and peat archives (Heiri et al., 2014a), noble gas proxy records from groundwater and speleothems (Beyerle et al., 1998; Ghadiri et al., 2018), and tentative inverse glacial modeling (Kerschner and Ivy-Ochs, 2008; Becker et al., 2016; Seguinot et al., 2018), with some noticeable variability in derived paleoclimate information between and within proxy records. Proposed reconstructed mean temperature anomalies during the LGM hence vary from −11 to −14 ^∘C based on pollen reconstructions (Wu et al., 2007; Bartlein et al., 2011), −5 to −9 ^∘C based on noble gas groundwater records (Beyerle et al., 1998, Seltzer et al., 2021), and −8 to −15 ^∘C using glacial modeling studies calibrated on reconstructed ice limits and estimates of paleo-equilibrium line altitude (ELA; i.e., the elevation at which the annual net ice budget in a glacier equals zero; Allen et al., 2008; Becker et al., 2016; Seguinot et al., 2018; Višnjević et al., 2020). LGM precipitation conditions are even more uncertain, with estimates for precipitation anomalies varying widely between −20 % and −60 % (e.g., Peyron et al., 1998; Luetscher et al., 2015; Becker et al., 2016) and for which a differential north–south distribution pattern (Florineth and Schlüchter, 2000; Luetscher et al., 2015; Becker at al., 2016) is still debated (Seguinot et al., 2018; Višnjević et al., 2020). Similarly, little is known regarding climatic conditions during the Late Glacial period between the LGM and the YD; besides that, significantly lower summer temperatures (> 6 ^∘C negative anomalies) were still persisting before ca. 15 ka based on chironomid and treeline proxies (Heiri et al., 2014a). During the short-lived (∼ 1 kyr) YD cooling event, temperatures dropped, with mean annual anomalies varying between 2–3 and 5–9 ^∘C below present-day values, depending on the considered proxy between paleoglacial reconstructions (e.g., Protin et al., 2019; Baroni et al., 2021), lacustrine pollen assemblages (Magny et al., 2001), and noble gas speleothem records (Ghadiri et al., 2018; Affolter et al., 2019). On the other hand, for the Holocene period, all the available records are in general agreement, indicating that temperature conditions relatively similar to today prevailed, with only minor (less than 2 ^∘C) deviations (e.g., Davis et al., 2003; Heiri et al., 2014a; Ghadiri et al., 2018; Affolter et al., 2019).

There is a crucial lack of direct and quantitative in situ temperature constraints from within the Alpine massifs during the different reconstructed glacial stages since the LGM. In this study, we attempt to reconstruct paleotemperatures in the high Alps during the Late Glacial period by applying cosmogenic noble gas paleothermometry (Tremblay et al., 2014a). This method exploits the diffusive behavior of cosmogenic ³He in quartz minerals at Earth surface temperatures (Brook et al., 1993; Shuster and Farley, 2005). Using forward models of cosmogenic ³He production and thermally activated diffusive loss, quantitative constraints on the thermal history of an exposed rock surface can thus be inferred from the difference between surface-exposure ages derived from the diffusive ³He system and from a cosmogenic nuclide that does not experience diffusive loss (Tremblay et al., 2014a, b, 2018). Cosmogenic ³He paleothermometry provides a unique opportunity to obtain quantitative information about past temperatures from in situ rock surfaces located in the formerly glaciated Alps. Here, we explore the applicability of cosmogenic ³He paleothermometry along two deglaciation profiles in the northern and western Alps. The advantages of such sampling targets are (1) relatively simple exposure history of rock surfaces revealed between the LGM and YD (Wirsig et al., 2016b; Lehmann et al., 2020), with limited shadowing effect (e.g., steep surface, limited vegetation, or postglacial sediment cover), and (2) the access to sequences of gradually exposed but lithologically similar samples, enabling a semi-continuous record of a temperature-change history. Based on ³He analytical measurements and forward model simulations, we aim to investigate the sensitivity of the in situ quartz–³He system in two different high Alpine areas and its suitability for the preservation of a ³He thermal signal on Late Pleistocene timescales. We also compare our results to previous studies applying cosmogenic ³He paleothermometry elsewhere in the Alps to gain a further understanding of ³He diffusion behavior in quartz at Earth surface temperatures.

2 Study sites

Two sites located in major Alpine massifs were selected for this study and have been previously investigated for their deglaciation history: the Mont Blanc Trelaporte (MBTP) profile (Mont Blanc massif, France; Lehmann et al., 2020), located in the western Alps along the western flank of the Mer de Glace valley (NNE exposure); and the SW-exposed Gelmersee (GELM) ridge (Aar massif, Switzerland; Wirsig et al., 2016b), formed by a hanging valley on the east wall of the Haslital valley in the northern central Alps (Fig. 1, inset). Both sites have steep valley sides that are several hundred meters high with ∼ 30–35^∘ slopes and are characterized by smoothly abraded rock surfaces and “roche moutonnée”-like features molded by flowing glaciers. Homogeneous lithologies are exposed along the valley walls, with phenocrystalline granite of the Mont Blanc (Dobmeier, 1998) at the MBTP site and Aare granite (Labhart, 1977; Abrecht, 1994) as part of the Helvetic crystalline basement at the GELM site. At both sites, the upper parts of the valley sides are characterized by jagged rock surfaces resulting from active periglacial processes. The trimline, which is the transition between smooth and rough bedrock surfaces, is located at ∼ 2600 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ at the MBTP site and ∼ 2450 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ at the GELM site. This trimline either marks the upper limit of the LGM ice surface or was a subglacial boundary marking the limit between warm-based eroding ice and cold-based ice (Wirsig et al., 2016a). The Mer de Glace valley is still occupied by ice today, with the ice limit at ∼ 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ near our profile site, while the Haslital valley is fully deglaciated. Continuous permafrost is expected above ∼ 3000 m a.s.l in the north faces of the Mont Blanc massif (permafrost index ≥ 0.9; Magnin et al., 2015a) but can be found more discontinuously down to 2300 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (permafrost index ≥ 0.5) and as low as 1900 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ in especially favorable conditions (permafrost index ≥ 0.1). Along the Gelmersee ridge on the western side of the Haslital Valley, continuous permafrost is expected above ∼ 2700 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, with sporadic patches down to ∼ 2150 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (Boeckli et al., 2012b).

Figure 1 Mont Blanc Trelaporte (MBTP, a) and Gelmersee (GELM, b) deglaciation profiles since the Last Glacial Maximum (LGM), with the spatial distribution of samples collected for quartz ¹⁰Be (Lehmann et al., 2020; Wirsig et al., 2016b) and ³He (this study) analyses. Samples with an asterisk have been exposed for ∼ 10–11 kyr (i.e., the entire Holocene period). The inset map indicates the location of the two study sites within the European Alps and the extent of ice cover during the LGM (in blue; Ehlers et al., 2011).

Ice-surface lowering of around 400 (MBTP) to > 500 (GELM) meters between the LGM and the YD has been recorded using in situ ¹⁰Be cosmogenic exposure dating on bedrock surfaces collected at regular intervals along each profile, starting from just below the trimline (Figs. 1 and 2, Table 1; Lehmann et al., 2020; Wirsig et al., 2016b). In this study, new samples were collected for ³He experiments from the same rock surfaces and locations as the sampling sites previously collected for ¹⁰Be dating by Lehmann et al. (2020; MBTP profile, samples MBTP18-1, -2, -11, and -9, n=4) and Wirsig et al. (2016b; GELM profile, samples GELM18-12, -1, -5, -11, -6, and -9, n=6; Fig. 1, Table 1). All samples are from glacially scoured bedrock surfaces, except GELM18-11, which comes from the top of a ∼ 5 m high boulder of similar lithology deposited during the post-LGM ice-surface lowering (Wirsig et al., 2016b).

Figure 2 Apparent quartz ³He (this study) and ¹⁰Be (Lehmann et al., 2020; Wirsig et al., 2016b) exposure ages (a, c), and ${}^{\mathrm{3}}\mathrm{He}{/}^{\mathrm{10}}\mathrm{Be}$ exposure age ratios or retention (b, d) as a function of elevation along the two deglaciation profiles (MBTP: a, b; GELM: c, d).

Table 1 MBTP and GELM sample information.

^a Samples used for ³He diffusion experiments. ^b Re-calculated ¹⁰Be exposure ages (after Wirsig et al., 2016b and Lehmann et al., 2020) and calculated ³He apparent exposure ages using the non-time dependent scaling scheme of Stone (2000) and Balco et al. (2008), using SLHL production rates of 4.01 $\mathrm{at}{\mathrm{g}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ (¹⁰Be; Borchers et al., 2016) and of 116 $\mathrm{at}{\mathrm{g}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ (³He; Vermeesch et al., 2009) and assuming a rock density of 2.65 g cm⁻³. See the Supplement for the details of ¹⁰Be and ³He concentrations (Table S1 in the Supplement). ^c Estimated modern mean annual rock surface temperature (MARST) at ∼ 3 cm depth. ^d Modern effective diffusion temperature (EDT) calculated using E_a of 93.5 (MBTP) or 98.5 (GELM) kJ mol⁻¹ and using 10 ^∘C annual and 5 ^∘C diurnal amplitudes. See Sect. 3.5 for the detailed explanation regarding modern MARST and EDT estimates.

3 Methods

In addition to measurements of cosmogenic ¹⁰Be, to determine the exposure time of a rock surface, cosmogenic ³He paleothermometry requires at least two additional types of measurement. First, we need measurements of the cosmogenic ³He concentration in the quartz, which permits us to estimate the amount of cosmogenic ³He loss by diffusion during the exposure of the rock surface. Second, we need to measure sample-specific ³He diffusivities for varying temperatures in our quartz samples by conducting stepwise-heating experiments.

We also need several different types of models to obtain paleotemperature estimates from these datasets. First, we need models to obtain diffusion kinetics parameters (i.e., activation energy, E_a, and the length-scale-normalized diffusivity at infinite temperature, also known as the pre-exponential factor, $D_{\mathrm{0}}/a^{\mathrm{2}}$; Tremblay et al., 2014a, b) from our stepwise-heating experiments. Some quartz samples, including the samples we study here, show complex ³He diffusion behavior that prevents us from using a simple linear model to extract diffusion kinetics parameters. In these cases, we model the diffusion kinetics parameters using the multiple diffusion domain (MDD) model framework (Lovera et al., 1989, 1991). Second, we forward model the production and diffusion of cosmogenic ³He over the duration of each sample's surface exposure for different thermal histories to compare the predicted ³He concentrations with the cosmogenic ³He concentrations observed in our samples. In order to account for diurnal and seasonal temperature oscillations, effective diffusion temperatures (EDTs, Tremblay et al., 2014a) are used as temperature inputs in the forward models.

We describe the measurement and modeling methods we used to apply cosmogenic noble gas paleothermometry to our samples as well as how EDTs were derived in greater detail below.

3.1 Sample preparation

Rock samples were disaggregated using a high-voltage, pulse-based system (SELFRAG equipment, Institute of Geological Sciences, University of Bern) to optimize the breaking of the rock along crystal grain boundaries. After rinsing, quartz mineral grains were separated from other minerals (heavy minerals and feldspar) by magnetic separation and froth flotation (e.g., Nichols and Goehring, 2019). The quartz separates were etched in 1 % HF for 3 weeks at room temperature to ensure the removal of any adhering micro-mineral particles that have the potential to contaminate the ³He measurements. For both sites, the grain size distribution is centered around a diameter of 850 µm after removal of the fraction finer than 200 µm (Fig. S1 in the Supplement).

For cosmogenic ³He measurements, the 800–1000 µm grain size fraction (i.e., 400–500 µm radii) from each quartz separate was selected. We chose this large grain size fraction, because we anticipated it would have best preservation potential of a measurable ³He signal, given the expected range of thermal histories experienced by the MBTP and GELM samples (Brook and Kurz, 1993; Tremblay et al., 2014a). Three replicates per sample consisting of ∼ 100 mg of quartz were prepared for analysis of natural ³He concentrations.

One representative sample per profile was selected for stepwise-heating experiments to determine diffusion kinetics parameters: MBTP18-9 and GELM18-1. For these samples, 200 to 300 mg of quartz grains were visually selected under a binocular microscope to avoid obvious mineral inclusions and fractures. The selected grains were sent to the Francis H. Burr Proton Therapy Center at the Massachusetts General Hospital for proton beam irradiation (Shuster et al., 2004; Shuster and Farley, 2005) in February 2019. After several months of rest to lower the level of radioactivity, one individual coarse quartz grain with no obvious fractures, mineral inclusions, or fluid inclusions was selected from each irradiated sample to conduct stepwise-heating experiments. The grain selected from MBTP18-9 had a ∼ 700 µm diameter, while the grain selected from GELM18-1 had a ∼ 900 µm. Grain diameters were estimated using calibrated petrographic microscope measurements.

3.2 Helium measurements

Both bulk degassing measurements to determine the natural cosmogenic ³He abundances and the stepwise-heating experiments on proton-irradiated quartz grains to characterize ³He diffusion kinetics were carried out at the BGC Noble Gas Thermochronometry Lab (Berkeley, USA). The measurements were conducted with an MAP 215-50 sector field mass spectrometer following a similar procedure to Tremblay et al. (2014b). For bulk degassing measurements, the samples were packaged into tantalum packets and heated in two, 15 min-long heating steps at 800 and 1100 ^∘C with a diode laser, with temperature of the tantalum packet measured by pyrometry. The amounts of ³He and ⁴He released from each heating step were measured (Tables S1–S3). Hot blanks on empty tantalum packets were also analyzed, from which an averaged ³He blank correction of 7.7 × 10³ atoms was obtained. No ³He above blank levels was observed in any of the 1100 ^∘C heating steps. For stepwise-heating experiments on proton-irradiated grains, the selected grains were placed in contact with the tip of a bare wire K-type thermocouple inside small platinum–iridium (PtIr) packets. The PtIr packets were heated with a diode laser in a feedback control loop with the thermocouple. Each experiment included over 30 to 40 heating steps of varying durations, with heating step temperatures between 70 and 550 ^∘C and including at least one retrograde heating cycle (Tables S4 and S5). Blank measurements at room temperature were regularly conducted throughout the experiments for background subtraction from the measured raw signals, with averaged ³He blank corrections of 2.1 × 10⁴ atoms (MBTP18-9) and 4.9 × 10⁴ atoms (GELM18-1).

3.3 Diffusion kinetics determination

To obtain diffusion kinetics parameters (E_a and $\ln (D_{\mathrm{0}}/a^{\mathrm{2}})$) from the stepwise-heating experiments on proton-irradiated grains, we followed a multi-step procedure. For each proton-irradiated sample (one per site), we first produced an Arrhenius plot displaying the natural log of diffusivity D (scaled to the diffusion length scale a) as a function of inverse temperature (Fig. 3), calculated from the observed ³He release fractions using the equation of Fechtig and Kalbitzer (1966; in Tremblay al., 2014b). The Arrhenius plots for both samples are shown in Fig. 3. In the case of simple diffusion behavior that follows an Arrhenius relationship, all data points would form a single linear array in this plotting space, and single values of E_a and $\ln (D_{\mathrm{0}}/a^{\mathrm{2}}$) could be obtained from the slope and intercept of a linear fit, respectively. However, the Arrhenius plots for both samples show deviations from linearity (Fig. 3).

Figure 3 Arrhenius plots of ³He step-degassing experiments conducted on one representative sample per study site (a: MBTP18-9, b: GELM18-1). Gray circles show $\ln (D/a^{\mathrm{2}})$ values calculated from the laboratory experiments after Fechtig and Kalbitzer (1966). Diffusion kinetics (DK) parameters were calculated using a multi-step procedure. First, we fit the laboratory data assuming a range of E_a values (79 a and 86 b to 100 kJ mol⁻¹; gray lines) for a three-domain MDD model. The red lines show three-domain DK parameters that minimize the misfit between the observed and predicted diffusivities (E_a is proportional to the slope of each line, and $\ln (D_{\mathrm{0}}/a^{\mathrm{2}})$ is given by the intercept). Second, we assessed whether these three-domain DK parameters could predict the cosmogenic ³He concentrations recorded in the Holocene calibration samples (c: MBTP18-9, d: GELM18-6, and e: GELM18-9). The three-domain DK parameters that could reproduce the Holocene calibration data are represented with light blue lines, with the dark blue line indicating the final calibrated DK parameters producing the best match with the natural Holocene ³He concentration(s). The small black circles in panels (a) and (b) represent the $\ln (D/a^{\mathrm{2}})$ values modeled along the heating experiment schedule for this best-fit set of three-domain DK parameters. The vertical green line in panels (a) and (b) indicates the modern effective diffusion temperature (EDT) associated with the Holocene calibration sample(s). The gas fraction assigned to each domain for both laboratory-only (gray) and Holocene-calibrated (blue) DK parameters is also indicated along the model lines in panels (a) and (b).

Given the observation of nonlinear Arrhenius behavior, we follow the approach of Tremblay et al. (2014b) and use a multi-diffusion domain (MDD; Lovera et al., 1989, 1991) model framework to determine quartz ³He diffusion kinetics parameters for each study site. In the MDD framework, ³He diffusion is modeled as occurring in several non-interacting domains, with different effective diffusion length scales within the quartz grain (e.g., sub-grain fragments). Such a model can reproduce an observed nonlinear Arrhenius behavior and has been shown to be consistent with cosmogenic ³He observations in a geologic case study (Tremblay et al., 2014b). We use the code from Tremblay et al. (2021) to perform the MDD analysis on our stepwise-heating experiment data.

Preliminary MDD calculations were carried out to determine the E_a that best fitted the Arrhenius data points in the lower temperature range (∼ 70 to 100 ^∘C), assuming a single diffusion domain, as well as the (minimum) number of diffusion domains to explain the entire dataset (i.e., all heating steps). Additional iterative runs using the MDD model with the minimum number of domains inferred from the preliminary tests were then conducted for a range of increasing E_a up to 100 kJ mol⁻¹ (with a 0.5 kJ mol⁻¹ increment; Fig. 3). This range of E_a values is based on existing E_a estimates reported for quartz in the literature (Shuster and Farley, 2005; Tremblay et al., 2014b, 2018; Domingos et al., 2020). In each run, E_a was kept common to each domain (Lovera et al., 1991; Baxter, 2010), while $\ln (D_{\mathrm{0}}/a^{\mathrm{2}})$ and the gas fraction for each different domain were allowed to vary until the misfit coefficient was minimized between the simulated and observed $\ln (D/a^{\mathrm{2}})$ values for all the heating steps.

Stepwise-heating experiments conducted in a laboratory do not permit us to observe ³He diffusion behavior at Earth surface temperature range (i.e., from around −30 to 30 ^∘C). Considering this fact, we introduced an extra calibration step that uses the measured natural ³He concentration from the samples with Holocene-only exposure (from both GELM and MBTP sites) to constrain the diffusion kinetics previously inferred from MDD models. The rationale for this is as follows: based on independent global and regional paleoclimate proxy records, samples exposed during the Holocene have experienced relatively stable average temperature conditions, with only minor variations (i.e., less than 2 ^∘C; see Sect. 1). Assuming no complex exposure history, the ³He signal recorded in these samples should therefore be representative of ³He diffusion occurring at a constant temperature equivalent to the modern effective diffusion temperature (EDT) at each sample site (see Sect. 3.4 for full definition of EDT). We thus tested whether the sets of diffusion kinetic parameters from the MDD models could explain the natural ³He concentrations recorded in the Holocene samples from our two study sites (MBTP18-9, ¹⁰Be exposure age of ca. 11.6 ka; GELM18-9 and -6, ¹⁰Be exposure ages of ca. 10.2 and 11.2 ka). For each set of MDD diffusion kinetics parameters that minimized the misfit with the stepwise-heating experiment data, we used a forward model of ³He production and diffusion (Tremblay et al., 2021) to predict the concentration of ³He that would be expected for an exposure duration equivalent to the recalculated ¹⁰Be exposure age and for a constant temperature equivalent to the modern EDT of the Holocene sample(s) (insets Fig. 3). Diffusion kinetics parameters for which the modeled ³He concentration matched the observed cosmogenic ³He concentration in the Holocene calibration sample within error were retained. We considered the diffusion kinetics parameters that yielded the best match as the final calibrated diffusion kinetics parameters (Fig. 3). We assume that the Holocene-calibrated diffusion kinetics parameters apply to all the samples collected along each profile and use these diffusion kinetics parameters in subsequent calculations. We justify the assumption of common diffusion kinetics for samples from a same valley profile by the homogenous lithology observed between those samples.

3.4 Forward models of ³He production and diffusion

To explore what thermal histories could explain the observed cosmogenic ³He abundances in our quartz samples, we forward model the production and diffusion of cosmogenic ³He over the duration of each sample's surface exposure for different time–EDT scenarios using the approach of Tremblay et al. (2021). We model ³He production by cosmic ray incidence using a ³He production rate in quartz at sea level and a high latitude (SLHL) of 116 $\mathrm{at}{\mathrm{g}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ (Vermeesch et al., 2009), scaled to the sample geographic location and elevation according to the non-time-dependent scaling scheme of Stone (2000) and Balco et al. (2008). The duration of ³He production and diffusion in the simulations is constrained by each sample's ¹⁰Be exposure age. For consistency, the SLHL production rate of ¹⁰Be (4.01 $\mathrm{at}{\mathrm{g}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ by neutron spallation; Borchers et al., 2016) is also scaled with the Stone (2000) scaling scheme. Apparent ³He and ¹⁰Be exposure ages along the deglaciation profiles investigated in this study are (re-) calculated using the measured ³He (this study) as well as ¹⁰Be concentrations reported in the literature (previous studies, Wirsig et al., 2016b; Lehmann et al., 2020), assuming negligible erosion (Fig. 2; Tables 1 and S1).

3.5 Effective diffusion temperature estimates

Rock surfaces experience temperature fluctuations at the diurnal, seasonal, and longer (1 to 10⁵ years) timescales, which will all activate thermal diffusion of ³He in quartz (Tremblay et al., 2014a). Because ³He diffusivity increases exponentially with temperature, a constant model temperature required to explain a total ³He loss (i.e., corresponding to the mean diffusivity through time) from a geological sample will equal or exceed the actual mean temperature experienced at the rock surface. This temperature is called effective diffusion temperature (EDT; Christodoulides et al., 1971; Tremblay et al., 2014a) and is a function of the ³He diffusion activation energy E_a, the long-term mean (rock surface) temperature, and the different frequency temperature amplitudes.

We estimated the modern EDT at the different sampling sites as follows. Mean annual air temperatures (MAATs) at each sampling site along the MBTP and the GELM profiles were calculated by linear interpolation, assuming a lapse rate of 5 ${}^{\circ }\mathrm{C}{\mathrm{km}}^{-\mathrm{1}}$ (Grämiger et al., 2018) based on mean annual temperatures recorded by nearby reference weather stations at Chamonix (1042 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, ∼ 5 km west; period 1993–2012; Magnin et al., 2015a) and Grimsel-Hospiz (1980 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$; ∼ 5 km south; period 2010–2020, data MeteoSwiss), respectively. Mean annual rock surface temperatures (MARSTs) are typically higher than MAATs, with the difference amplified between south- and north-exposed slopes (Gruber et al., 2003). Boeckli et al. (2012a), based on 57 sensor measurements on snow-free rock slopes > 55^∘, showed that the measured difference between MAAT and MARST increased linearly from < 1 to up to 10 ^∘C, depending on potential incoming solar radiation (PISR), which is largely controlled by rock surface aspect and angle in addition to elevation. For moderately inclined surfaces, the difference between MARST and MAAT is expected to be reduced by ∼ 1–3 ^∘C due to micro-topography and snow-insolating effects (Hasler et al., 2011). To estimate MARSTs, we calculated the PISR at each sampling site using the area solar radiation tool (ArcGIS software, version 10.3.1) applied to a 30 m resolution digital elevation model (SRTM 1 Arc-Second data) at the study sites. The calculation was performed at hourly resolution using data from 1 year (2000), assuming no nebulosity and using a sky size of 512 cells (Magnin et al., 2015a; Mair et al., 2020). Based on the linear relationship between MAAT–MARST and PISR from Boeckli et al. (2012a), we estimated the average MARST–MAAT difference, assuming snow-free conditions at each site, from which we then subtracted 2 ^∘C to take into account snow-insulating and micro-topographic effects in moderately steep terrain (Hasler et al., 2011). Final differences between MAAT and MARST of +1 and +2.5 ^∘C were thus obtained for the north-exposed MBTP and the southwest-exposed GELM sites, respectively. These estimates are consistent with in situ MAAT and MARST measurements available in nearby areas with similar orientations, elevations, and slope inclinations (e.g., Gruber et al., 2004; Magnin et al., 2015a, b; Haberkorn et al., 2017; Grämiger et al., 2018; Guralnik et al., 2018) and were thus used to estimate the MARSTs at each sampling site.

Table 2 Diffusion kinetics parameters for MBTP and GELM sites.

^a Diffusion kinetics measurements made on one representative sample per profile: MBTP18-9 (350 µm spherical equivalent radius) and GELM18-1 (450 µm spherical equivalent radius). ^b Range of MDD diffusion kinetics parameters obtained by fitting laboratory experimental data and matching ³He concentrations (within 1σ error) from Holocene calibration samples. ^c Best-fitting MDD diffusion kinetics parameters obtained by fitting laboratory experimental data matching ³He concentrations from Holocene calibration samples. ^d MDD diffusion kinetics parameters based only on laboratory experimental data and providing the best match in the lower temperature range of the heating schedule (∼ 70–100 ^∘C).

A mean annual temperature amplitude of 10 ^∘C and diurnal amplitude of 5 ^∘C were adopted for the two sites, based on long-term (i.e., several years) temperature records from the Chamonix and Grimsel-Hospiz weather stations and from direct in situ rock surface measurements available in the Alps (Gruber et al., 2004; Magnin et al., 2015b; Grämiger et al., 2018; Guralnik et al., 2018; Mair et al., 2020). These estimates are consistent with the annual and diurnal amplitudes obtained from the spatially interpolated land surface climate dataset, WorldClim 2.0, based on gridded time series of meteorological data from available weather stations (target temporal range 1970–2000; 1 km resolution; Fick and Hijmans, 2017).

Past colder EDT inputs in forward ³He modeling for varying thermal histories are based on temperature anomalies from modern EDTs.

4 Results

First, we examine the characteristics of the ³He diffusion kinetics parameters we modeled for our quartz samples and explore the sensitivity of the ³He signal in those samples to Earth surface EDTs. We then present forward model results for the evolution of the cosmogenic ³He concentrations recorded along each deglaciation profile for two different sets of thermal histories. The first set of thermal histories we investigate assumes a constant EDT since the exposure of the sampled rock surfaces following ice retreat. We then investigate a set of more climatologically interesting thermal histories, wherein a change in EDT occurs at some point during the exposure time of each sample.

4.1 Diffusion kinetics and sensitivity tests

Figure 3 shows the range of diffusion kinetics parameters (E_a and ln($D_{\mathrm{0}}/a^{\mathrm{2}}\left)\right)$ that fit the laboratory stepwise-heating experiments (one representative sample for each site; Fig. 3a and b) and which permit us to reproduce the observed natural ³He concentrations from the Holocene calibration samples for a constant EDT equivalent to the modern EDT (Fig. 3c–e). Stepwise-heating experiment data indicate relatively first-order Arrhenius behavior for quartz ³He diffusion of MBTP18-9, with one dominant linear array accounting for ∼ 85 % of ³He release (Fig. 3a, Table 2). The remaining ∼ 15 % gas fraction is distributed within two additional minor diffusion domains, one of higher retentivity and one of lower retentivity (Tremblay et al., 2014b). GELM18-1 exhibits more complex quartz ³He diffusion behavior, with gas release distributed more equally between three linear arrays (Fig. 3b, Table 2) that can be interpreted as three or more distinct diffusion domains, with each domain contributing significantly to ³He retention over geological times.

In order to explore the theoretical sensitivity (and potential variability) of the MBTP and GELM quartz, we numerically evaluated the time required for the concentration of ³He in each sample to reach steady-state conditions (i.e., thermal loss balanced with cosmic-ray induced production gain) as a function of constant EDT. Forward simulations of ³He production and diffusion (Tremblay et al., 2014a, b, 2021) using the final Holocene-calibrated diffusion kinetic parameters were thus run for a range of isotherms representative of Earth surface EDTs (hereafter referred to as isoEDTs; tested range from −30 to 30 ^∘C), assuming a 450 µm radius and no initial ³He concentration. Equilibrium conditions were assumed to be reached once no significant change in ³He concentration was recorded (< 1 % kyr⁻¹). While we observe some variability in ³He diffusion behavior and derived diffusion kinetics parameters between MBTP and GELM quartz (Fig. 3, Table 2), results from sensitivity tests in terms of steady-state times are relatively similar. For isoEDTs between −10 and 10 ^∘C, bracketing approximately potential EDT values experienced along both deglaciation profiles between the LGM and today, the time predicted for ³He diffusion to reach equilibrium varies between ∼ 10 kyr (isoEDT of 10 ^∘C) and ∼ 20 kyr (isoEDT of −10 ^∘C; Fig. 4). Interestingly, while steady-state time estimates remain relatively constant for quartz from both sites at ca. 20 kyr for colder isoEDTs (−10 to −30 ^∘C), we observe a pronounced non-linear dependence for EDTs above 0 ^∘C, resulting in much shorter equilibrium times in the high EDTs range (less than 5 kyr for EDT above 20 ^∘C, Fig. 4).

Figure 4 Theoretical ³He steady-state time estimates for isoEDTs varying between −30 and 30 ^∘C using the final Holocene-calibrated diffusion kinetics parameters determined for each study site and assuming 450 µm grain radius.

4.2 ³He exposure ages and PaleoIsoEDTs

For each site, apparent ³He exposure ages are systematically lower (from 20 % to 75 %) than apparent ¹⁰Be exposure ages (Table 1, Fig. 2). The ¹⁰Be ages show a general decrease with decreasing elevation, in agreement with progressive ice thinning along a deglaciation profile in the high Alps during the Late Glacial. This trend is less evident for the apparent ³He ages, which overlap significantly within uncertainties above ∼ 2200 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (Fig. 2). The ³He retention (${}^{\mathrm{3}}\mathrm{He}{/}^{\mathrm{10}}\mathrm{Be}$ exposure age ratio) shows a clear decrease with decreasing elevation along the GELM profile (∼ 1500 to 2500 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$), which is not visible along the MBTP profile (similar retention between MBTP samples), which is also more restricted in elevation range (∼ 2100 to 2600 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$; Fig. 2).

To determine the apparent constant EDT (that we refer to as paleoIsoEDT) from the natural ³He signal recorded in each sample, forward models of ³He production and diffusion (implemented with the final Holocene-calibrated diffusion kinetics) were run for a time period equal to the sample's ¹⁰Be exposure age and for a range of isoEDTs (isothermal holding between −10 and 15 ^∘C for MBTP; −25 to 20 ^∘C for GELM, 1 ^∘C increment). The isoEDT leading to best-matching synthetic ³He concentration with the observed natural ³He concentration was retained as the paleoIsoEDT (Fig. 5, Table 1). As Holocene samples were used to calibrate the diffusion kinetics (see Sect. 4.1), it is expected that ³He-derived paleoIsoEDTs from these samples are equivalent to their respective modern EDTs. On the other hand, all pre-Holocene samples at both sites have paleoIsoEDTs that are lower than their corresponding modern EDTs (Fig. 5a and d; Table 1). For the MBTP profile, the difference between modern EDTs and paleoIsoEDTs varies from around 3 to 9 ^∘C. This difference is even greater for the GELM profile, where paleoIsoEDTs are around 12 to 20 ^∘C lower than their associated modern EDTs. Pre-Holocene samples are located well above (200 to 500 m) Holocene samples and all above 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$. While paleoIsoEDTs derived from the high-elevation, pre-Holocene samples agree within error for each site, they clearly depart from paleoIsoEDTs obtained from the low-elevation, Holocene sample(s) by ∼ 6^∘ (MBTP site) and ∼ 18 ^∘C (GELM site), based on the peak values from the obtained bimodal probability distributions (Fig. 5b and e). After correcting for temperature decreases with elevation (assuming a lapse rate of 5 ${}^{\circ }\mathrm{C}{\mathrm{km}}^{-\mathrm{1}}$), the difference in the paleoIsoEDTs between pre-Holocene and Holocene samples is still significant for GELM (∼ 10 to 20 ^∘C, Fig. 5f). For MBTP, although elevation-corrected paleoIsoEDTs from two high-elevation, pre-Holocene samples (MBTP18-2 and -11) are still clearly distinguishable from the low-elevation, Holocene sample (MBTP18-9; Fig. 5c), the probability distribution appears closer to unimodal, since the paleoIsoEDT from the highest sample (MBTP18-1) partially overlaps that from MBTP18-9 within uncertainty.

Figure 5 Distribution of ³He-derived paleoIsoEDTs along the MBTP (a–c) and GELM (d–f) deglaciation profiles. Holocene samples used for calibration are marked by asterisks. (a, d) PaleoIsoEDTs relative to modern EDTs (black dashed line is 1:1); (b, e) relationship between paleoIsoEDT and elevation; and (c, f) relationship between paleoIsoEDT and ¹⁰Be exposure age after correction for lapse rate (marked by red arrows, relative to the Holocene sample marked with a black asterisk). Green solid and dashed lines are modern EDTs for Holocene samples, with the solid line indicating the modern EDT taken as reference for the lapse rate correction. The thin lines represent the sums (yellow and orange) of the individual (gray) probability distributions of paleoIsoEDTs.

4.3 Forward simulations with time-varying EDT

Based on global and regional paleoenvironmental records, we can expect that pre-Holocene samples collected at the MBTP and GELM sites have experienced at least one main significant temperature change, marking the transition from (colder) Late Glacial to warmer and more stable Holocene conditions (cf. Sect. 1 for details).

Following this observation, we first investigate the theoretical time needed for the MBTP and GELM ³He quartz systems to re-adjust to a change in temperature for a warming scenario, assuming that these systems were already at steady-state conditions. Forward model simulations (450 µm radii assumed) were run for different time–EDT scenarios involving an initial EDT (ranging from −30 to 30 ^∘C; initial ³He concentration at steady state with initial EDT) followed by a step-warming event (+2, +5, +10 +15, +20, +25, +30 ^∘C; over 0.1 to 1 kyr, depending on the sensitivity of the quartz system for the considered EDT scenario), after which the resulting warmer EDT was maintained until full re-adjustment of the ³He-quartz system. We considered full re-adjustment to have occurred when modeled ³He concentrations following the step-warming event matched within 10 % the ³He concentrations expected for an isoEDT equivalent to the final (warmer) EDT (Fig. 6a). We present simulation results in Fig. 6b and c. For past EDTs < 0 ^∘C followed by a step warming up to 20 ^∘C, readjustment times are all longer than 10 kyr, considering either MBTP or GELM diffusion kinetics. Estimates of LGM-temperature anomalies suggested for the European Alps are equivalent to an apparent warming of 5 to 15 ^∘C (see Sect. 1). When considering EDT scenarios with a similar warming range applied to our study sites, with modern EDTs around 5 to 10 ^∘C (at pre-Holocene sampling sites; i.e., equivalent to initial past EDT of 0 to 5 and −10 to −5 ^∘C for 5 and 15 ^∘C warming step, respectively), our simulation outcomes show relatively long re-adjustment times from around 20 to 45 kyr (Fig. 6). We should note, however, that these times are maximum estimates, since we considered ³He quartz systems at steady-state conditions with initial cold EDTs before the warming event.

Figure 6 Conceptual approach (a) and output results for MBTP (b) and GELM (c) of ³He re-adjustment time (t_re-adjust) for one-step EDT change scenarios (temperature warming from 2 to 30 ^∘C), using the final diffusion kinetics parameters determined for each study site and assuming 450 µm grain radius. Calculations assume that ³He concentrations were already at steady state for past EDT conditions (i.e., as would be expected for infinite exposure time) prior to imposing the temperature change.

In a second set of forward model runs, we explore the thermal memory of ³He in quartz for a step-warming EDT scenario fixed in time that is more representative of the post-LGM paleoclimate history in the Alps, including: (1) an initial cold period starting at 24 ka, with an imposed EDT set 15 ^∘C lower compared to modern EDT (maximum LGM temperature anomalies, see Sect. 1); (2) a warming step to modern EDT that is either progressive from 24 to 10 ka or abrupt between 11 and 10 ka (i.e., consistent with a Younger Dryas–Holocene transition); and (3) stable conditions at the modern EDT throughout the Holocene (lasty10 kyr, Fig. 7a). Forward simulations of ³He diffusion and concentration evolution were conducted for each pre-Holocene sample following this scenario, with the time period and the time-dependent EDT variable set accordingly to each sample ¹⁰Be exposure age and modern EDT, respectively. For all GELM and MBTP samples, these forward simulations result in synthetic ³He concentrations significantly lower than their respective measured ³He concentrations. In Fig. 7b, we present the results for sample MBTP18-1, for which we observed the smallest difference between modern EDT and paleoIsoEDT (Fig. 5a; Table 1).

Figure 7 (a) Simplified warming EDT scenario since the LGM (∼ 24 ka), with progressive and abrupt EDT changes in light and dark blue lines, respectively; (b) synthetic evolution of ³He concentration (blue lines) compared with the natural ³He concentration recorded in MBTP18-1 (red circle). The ³He concentration evolution is also indicated for a constant-temperature scenario at the modern EDT and paleoIsoEDT (set in Fig. 5; black solid and dashed lines, respectively). We were unable to reproduce the observed natural ³He concentration for any samples with pre-Holocene exposure under this simplified LGM EDT scenario.

To further investigate potential effects of a larger EDT differences between modern and past conditions and/or a more recent EDT change, we performed an additional set of numerical simulations using step-warming EDT scenarios with more free parameters. Scenarios with an EDT change occurring from 10⁴ to 10² years ago and with a difference between past and modern EDTs (ΔEDT) up to 40 ^∘C were tested iteratively on each pre-Holocene sample, assuming no initial ³He concentration and with total exposure time and EDT variables adjusted accordingly, as described above. Scenarios for which we could reproduce the observed natural ³He concentration (within uncertainties) were accepted, resulting in a range of different possible scenarios with varying ΔEDT and time of EDT change for each pre-Holocene sample (Fig. 8). For both sites, we observed a similar pattern between ΔEDT and time of EDT change: the further back in time the EDT change occurs, the greater the ΔEDT that is needed to reproduce observed natural ³He concentrations. In addition, for any given time of EDT change, ΔEDTs tend to be inversely correlated with sample elevation and ¹⁰Be exposure age. Within these similarities, the two sites differ by the magnitude of the ΔEDT required to reproduce observed natural ³He concentrations. For example, along the MBTP site (Fig. 8a), ΔEDTs of 5 ^∘C occurring a few kyr ago are required to explain ³He concentrations measured in the highest and/or oldest sample (MBTP18-1), while ΔEDTs of 35 ^∘C occurring a few kyr ago are required to explain ³He concentrations measured in the lowest and youngest pre-Holocene sample (MBTP18-11). For the same sites, ΔEDTs of 3 and 15 ^∘C are required if the ΔEDT occurred within the last centuries. On the other hand, for the GELM site, our simulations found no ΔEDT solution if the EDT change is applied prior to 1 ka (within our ΔEDT limit of 40 ^∘C; except for GELM18-12; Fig. 8b). In the case of EDT change occurring within the last centuries, ΔEDTs for the GELM samples are significantly larger than for MBTP samples, with ΔEDTs between 15 and > 30 ^∘C required for the highest and/or oldest samples (GELM18-12 and -1). For the intermediate samples (GELM18-11 and -5), which also exhibit the greatest ³He-¹⁰Be age differences, numerical solutions could only be recovered for very recent EDT changes (≤ 200 yr) and with ΔEDT > 35 ^∘C.

Figure 8 One-step EDT change scenarios that reproduce the observed natural ³He concentration for each pre-Holocene MBTP (a) and GELM (b) sample, with ΔEDT solution as function of the time of EDT change. The error bars indicate all the possible ΔEDT solutions and the color circles indicate the best-matching scenario.

5 Discussion

5.1 Paleoclimatic interpretation of ³He signals

All studied samples indicate the preservation of a ³He concentration consistent with temperatures that are colder than present-day EDT conditions at both the MBTP and the GELM sites (paleoIsoEDTs ∼ 3–9 and ∼ 12–20 ^∘C lower than modern EDTs, respectively, Fig. 5). However, for both sites, the recorded ³He concentrations are apparently not concordant with simple time–EDT scenarios describing a plausible post-LGM mean temperature evolution in the European Alps (i.e., LGM mean temperature anomaly up to 15 ^∘C, Fig. 7). Even when allowing for a larger EDT difference between LGM and the present day (up to 40 ^∘C), modeled ³He concentrations remain significantly below the observed values at both sites. Such large EDT differences would not be supported by any set of mean temperature reconstructions for the European Alps since the LGM (e.g., Heiri et al., 2014a). Likewise, potential variation in seasonal temperature cannot contribute significantly to a larger pre-Holocene EDT anomaly. Indeed, global and regional paleoclimate studies rather suggest that a larger seasonal temperature amplitude occurred before the Holocene (e.g., Davis et al., 2003; Buizert et al., 2018), which would have the effect of increasing the paleoEDTs instead (Fig. 9).

Figure 9 Difference between EDT and mean annual temperature as a function of increasing seasonal temperature amplitude for different diurnal temperature amplitudes and assuming an E_a of 93.5 kJ mol⁻¹ (MBTP site, Holocene-calibrated diffusion kinetics). The yellow star indicates the conditions used to estimate the modern EDT at the sampling sites. Decreasing the annual and/or diurnal amplitude can yield to up to a ∼ 5 ^∘C decrease in the modern EDT. Similar results are obtained when using E_a = 98.5 kJ mol⁻¹ (GELM site, Holocene-calibrated diffusion kinetics).

We attribute the result of modeled ³He concentrations that are significantly lower than the observed ones to the damping effect of modeled exposure during the Holocene period, which is characterized by relatively stable mean temperature conditions similar to the present day. By damping effect, we mean that modeling ∼ 10 kyr of exposure at temperatures similar to today results in a partial to total readjustment of the cosmogenic ³He thermal signal, with little to no inherited signal memory of the prior exposure to colder Late Glacial conditions. This hypothesis appears to contradict our theoretical tests, which indicate that the ³He thermal signal inherited from past EDTs 10 to 15 ^∘C colder than today should a priori be (partly) preserved for 30–45 kyr under modern EDT conditions (Fig. 6). However, this time range relies on the assumption that bedrock surfaces were exposed for long enough to past colder conditions before the temperature change occurred in order to reach ³He steady-state concentrations. For both sample sites, the time required to reach steady state is around 20 kyr (Fig. 4). On the other hand, along the MBTP and GELM profiles, bedrock surfaces have not been exposed for more than 5–8 and 4–13 kyr before the Late Glacial–Holocene transition, respectively. This results in ³He accumulation up to 35 %–55 % (MBTP) and 30 %–85 % (GELM) of ³He steady-state concentrations when considering paleoIsoEDTs 10–15 ^∘C lower than present-day EDTs. In such a case, ³He re-adjustment time estimates to modern EDTs are predicted to be reduced by ∼90 % to 80 % for the MBTP site and by >90 % to 60 % for the GELM site, implying we should recover the dominance of Holocene temperature conditions in the ³He signal from the sampled bedrock surfaces.

Our observed ³He concentrations can be reproduced by forward simulations with an EDT change occurring on much more recent timescales (Fig. 8). For the MBTP site, a ΔEDT of 7 to 5 ^∘C within the last few thousand years to centuries predicts the observed natural ³He concentrations for two pre-Holocene samples: MBTP18-1 and -2. A ΔEDT of 12 to 8^∘ is required for MBTP18-11. Such a ΔEDT estimate, considering mean temperature fluctuations up to 2 ^∘C for the Holocene period (Davis et al., 2003), would also require variations in diurnal and/or annual temperature amplitudes to account for an additional 5 ^∘C ΔEDT. However, this would imply the lowering of both diurnal and annual temperature amplitudes to null before modern conditions (Fig. 9), which contradicts global and regional records that indicate an increased seasonality in the early Holocene compared to the present day (Davis et al., 2003; Buizert et al., 2018) and which would result in a larger ΔEDT. Furthermore, the forward simulations discussed here used diffusion kinetics calibrated on Holocene samples (Fig. 8). Therefore, allowing a significant EDT change over the last 10²–10³ years is in contradiction with our calibration approach (see Sect. 3). If instead we use diffusion kinetics solely derived from laboratory experiments without Holocene calibration (Fig. 3, Table 2), a ΔEDT of 15 ^∘C or greater is required to explain observed MBTP ³He concentrations for changes within the last 10³ years (Fig. S2b). Such large ΔEDTs are significantly greater than expected EDT variations from changes in mean annual temperatures and/or in annual or diurnal temperature amplitudes during the Holocene. Even greater ΔEDTs are needed to explain the observed GELM ³He concentrations using either diffusion kinetics approach (Fig. S3b). Both cases are clearly incompatible with plausible Holocene paleoclimatic histories.

5.2 Sources of uncertainty

Cosmogenic ³He paleothermometry is still in an early stage of development for applications to Quaternary geology (Tremblay et al., 2014a, b, 2018), and there are several aspects of our approach that are under-constrained, which could contribute to our estimates of unrealistically cold Late Glacial temperatures. Below, we explore how uncertainties related to (1) estimating modern EDTs, (2) interpreting cosmogenic nuclide measurements, and (3) determining helium diffusion kinetics in quartz could each have affected our production–diffusion modeling and therefore our paleotemperature estimates.

5.2.1 Modern EDT estimates

Potential uncertainties in modern EDT estimates, used to define the EDT of the recent and stable period in the step-warming EDT scenarios (Sect. 4.3), cannot be ruled out. In particular, it is not known to what extent present-day conditions (based on decadal direct air and ground temperature measurements; see Sect. 3.1) are representative over centennial to millennial time scales. Correcting for overestimated diurnal and/or annual temperature amplitudes and/or mean annual temperatures would result in lower modern EDTs (Fig. 9). Assuming an overestimate of 50 % in modern diurnal and annual temperature amplitudes and up to 2 ^∘C overestimate in MARST based on recorded mean temperature fluctuations (Davis et al., 2003; Ghadiri et al., 2018, 2020) and applied corrections to MAAT (see Sect. 3.1) would lead to ∼ 3.5 ^∘C lowering of modern EDTs for MBTP and GELM sites. Applying such an estimated correction to the recent-period EDT potentially permits us to resolve observed ³He concentrations for two of the MBTP samples (MBTP18-1 and -2), with ΔEDT of 5 to 10 ^∘C for a change occurring at ca. 10 ka (i.e., LGM scenario; Fig. S2c). It is also worth noting that natural ³He MBTP concentrations for those samples can be reproduced with minor ΔEDTs (≤ 1.5 ^∘C) over recent timescales (10²–10³ years). When using laboratory-derived diffusion kinetics without Holocene calibration, a priori more appropriate to explore recent EDT changes, only scenarios with more than −10 ^∘C ΔEDT within the last thousand years are accepted (Fig. S2d), inconsistent with paleoclimatic records over this recent time period. For GELM samples, correcting modern or recent EDT is not sufficient to reproduce the observed ³He concentrations with plausible ΔEDTs for an EDT change occurring at the Late Glacial–Holocene transition (ca. 10 ka; no solution) nor on more recent timescales (Fig. S3c and d).

5.2.2 Interpretation of cosmogenic nuclide measurements

Several sources of geological uncertainties may affect the results obtained in this study. First, our approach relies on the assumption that bedrock surfaces have experienced a simple exposure history along the time period recorded by ¹⁰Be concentrations, without pre-exposure or episodic coverage (i.e., non-erosive cold-based ice). Depth profiles of ¹⁰Be measurements on glacially polished bedrocks in the western Alps, with apparent exposure ages of 10–20 ka, indicate that an inherited ¹⁰Be concentration due to insufficient glacial erosion may persist and could lead to up to 9 % age overestimates (Prud'homme et al., 2020). Similarly, Wirsig et al. (2016b) suggested potential but limited pre-LGM (less than a few ka overestimate) inheritance for some GELM samples. While previous bedrock surface exposure would also imply an inherited ³He concentration, the latter would be subject to diffusion (partial or total) during glacier coverage, even at subzero temperatures and EDTs (Fig. S4). On the contrary, ¹⁰Be would experience only minor radioactive decay over 10–100 kyr timescales. This scenario of inheritance and/or complex exposure history would result in lower ³He concentrations recorded by bedrock surfaces regardless of the temperature history experienced by the rock surface during the total ¹⁰Be exposure period (i.e., lower ${}^{\mathrm{3}}\mathrm{He}{/}^{\mathrm{10}}\mathrm{Be}$ concentration ratio; Balco et al., 2016). This scenario is also valid for post-LGM episodic coverage. Such effects are, however, expected to be minor, considering the limited potential ¹⁰Be inheritance (< 10 %) from pre-LGM exposure as well as the unlikelihood of prolonged coverage of the relatively steep (i.e., no loose sediments or thick snow accumulation) and high (i.e., above tree-line) sampled bedrock surfaces. Moreover, attempting to correct for these processes would result in opposite effects than what we observed for MBTP and GELM samples, with even lower paleoIsoEDT estimates and greater ΔEDTs required for warming EDT scenarios to recover observed natural ³He concentrations.

An additional source of uncertainty is postglacial erosion of sampled bedrock surfaces, assumed to be negligible in this study. Based on a combined approach exploiting cosmogenic ¹⁰Be and optically stimulated luminescence (OSL) systems, Lehmann et al. (2020) suggested potential high postglacial erosion rates (above 3.5 mm kyr⁻¹) for low-elevation MBTP samples. Other regional estimates for crystalline bedrock commonly indicate Alpine postglacial erosion rates of 0.1 to 1 mm kyr⁻¹ (Kelly et al., 2006; Dielforder and Hetzel, 2014; Wirsig et al., 2016b), which is in line with estimates from other studies (André, 2002; Balco, 2011). Relatively low postglacial erosion rates are further supported along our study sites by the presence of still-visible glacial striations (Wirsig et al., 2016b). Applying an erosion correction (0.1 to 1 mm kyr⁻¹) will only moderately affect apparent ¹⁰Be exposure ages (< 1 ka change) and would result in lower predicted ³He concentrations compared to our observed ones.

In summary, geological uncertainties related to exposure history and postglacial surface erosion are generally small and overall do not resolve the significant discrepancy between the natural ³He signal recorded in pre-Holocene MBTP and GELM samples and modeled ³He concentrations from expected EDT histories.

On the other hand, some of the observed differences may relate to uncertainties regarding the ³He production rate ($P_{{}^{\mathrm{3}}\mathrm{He}}$) in quartz. Directly estimating $P_{{}^{\mathrm{3}}\mathrm{He}}$ in quartz from geological calibration sites is challenging, as ³He diffuses from quartz at Earth surface temperatures over a 10²–10⁴ year timescale. Alternative approaches using artificial targets (e.g., Vermeesch et al., 2009) or scaling $P_{{}^{\mathrm{3}}\mathrm{He}}$ measured in retentive minerals (i.e., olivine; e.g., Cerling and Craig, 1994; Goehring et al., 2010) have hence been used. While we adopted the Stone (2000)-scaled $P_{{}^{\mathrm{3}}\mathrm{He}}$ from Vermeesch et al. (2009; i.e., 116 $\mathrm{at}{\mathrm{g}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) in this study, a ∼ 10 % higher ³He production rate has also been proposed from olivine ³He measurements scaled to quartz (e.g., Masarik and Reedy, 1995; Ackert et al., 2011). Applying an increased $P_{{}^{\mathrm{3}}\mathrm{He}}$ (Stone-scaled $P_{{}^{\mathrm{3}}\mathrm{He}}$ = 128 $\mathrm{at}{\mathrm{g}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) in general leads to smaller ΔEDTs in order to match the measured ³He concentrations, as well as to an older range of possible times for the EDT change. For the MBTP site, however, we could not reproduce ³He concentrations for an EDT change at 10 ka (except for MBTP18-1; Fig. S2e, Holocene-calibrated diffusion kinetics). Likewise, for more recent changes (10²–10³ years; laboratory-derived diffusion kinetics without Holocene calibration; Fig. S2f), the resulting ΔEDTs (10 to 25 ^∘C) are still not compatible with plausible Holocene temperature conditions. Similar results were obtained for the GELM Late Glacial samples when adopting a 10 % increase in $P_{{}^{\mathrm{3}}\mathrm{He}}$ (Fig. S3e and f).

In addition to a higher cosmogenic ³He production rate, another possibility that we have not accounted for is non-cosmogenic sources of ³He, specifically nucleogenic ³He produced by (n,α) reactions with ⁶Li. Not accounting for nucleogenic ³He would result in lower true cosmogenic ³He concentrations, which would have the effect of reducing the ΔEDTs at our sample sites toward more realistic values. However, we think it is unlikely that there is significant nucleogenic ³He in our samples for several reasons. First, the ${}^{\mathrm{3}}\mathrm{He}{/}^{\mathrm{4}}\mathrm{He}$ ratios we measured during the 800 ^∘C heating step are on the order of 10⁻⁶ to 10⁻⁷ (Table S3). This is more than an order of magnitude above the ${}^{\mathrm{3}}\mathrm{He}{/}^{\mathrm{4}}\mathrm{He}$ ratio of ∼ 10⁻⁸ expected from $\mathrm{U}/\mathrm{Th}$ decay and ⁶Li neutron capture (Niedermann, 2002). Second, no ³He above the detection limit was measured in the 1100 ^∘C heating step despite nontrivial amounts of ⁴He being released in this step. This indicates that retentive mineral and fluid inclusions, if present in the samples, are not contributing a significant amount of non-cosmogenic ³He to the measured ³He amounts. Third, based on the diffusion kinetics of ³He in quartz, we anticipate that any nucleogenic ³He produced in the quartz itself over geologic timescales will be diffusively lost before the sampled rock surface is exhumed at near-surface temperatures. Furthermore, the production rate of nucleogenic ³He is low compared to the cosmogenic production rate of ³He. We do not have direct data of major and trace element for the MBTP and GELM samples in order to calculate the nucleogenic ³He production rate directly. However, we can say that a rough maximum estimate for the production rate of nucleogenic ³He in the GELM samples is ∼ 1 $\mathrm{at}{\mathrm{g}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$, which is based on a maximum Li concentration of 70 ppm for the Aare granite (Schaltegger and Krähenbühl, 1990) and the production rate estimate of Farley et al. (2006) for an “average” granite. This is 0.3 % of the local, scaled production rate of cosmogenic ³He for sample GELM18-9, which has the lowest cosmogenic ³He production rate of all of our samples. Given this maximum production rate estimate for nucleogenic ³He and using the Holocene-calibrated diffusion kinetics for our samples, we estimate that the maximum steady-state concentration of nucleogenic ³He is 2.8 × 10⁴ atoms g⁻¹, which is 2 orders of magnitude smaller than the measured ³He concentrations in our samples and well within the uncertainties of those measurements. It is therefore unlikely that not correcting for nucleogenic ³He affected our modeled ΔEDTs in any significant way.

5.2.3 ³He diffusion kinetics characterization

Noble gas diffusion in minerals is generally assumed to have an Arrhenius-type dependence on temperature, where diffusivity increases exponentially with temperature and inversely with the diffusion domain size (e.g., Baxter, 2010 and references therein). Interestingly, theoretical studies investigating the fundamentals of ³He diffusion in quartz predict considerably lower E_a (and much higher diffusivity) than expected when considering a perfect quartz crystal (∼ 20 to 50 kJ mol⁻¹; Kalashnikov et al., 2003; Lin et al., 2016; Domingos et al., 2020; Liu et al., 2021), the latter suggesting that no ³He should be retained over geological timescales at Earth surface temperatures. These results are, however, in contradiction with common observations of ³He retention in natural rock surfaces (e.g., Brook et al., 1993; Brook and Kurz, 1993; Tremblay et al., 2018) and with typical E_a values empirically determined from laboratory stepwise-heating experiments (between 70 and 100 kJ mol⁻¹; Shuster and Farley, 2005; Tremblay et al., 2014b). Furthermore, previous ³He stepwise-heating experiments conducted on quartz from various origins indicate a large variability in diffusion kinetics (i.e., E_a and D₀) and diffusion behavior, wherein some quartz samples exhibit complex ³He diffusion behavior, while others exhibit a simple, linear Arrhenius dependence (Tremblay et al., 2014b). Both the observed variability and the discrepancy with theoretical predictions suggest that ³He diffusion in natural quartz is largely governed by sample-specific crystal defects (e.g., structural defects, radiation damages; Domingos et al., 2020), advocating for the use of sample-specific diffusion kinetics (Tremblay et al., 2014b). Complex, non-linear diffusion behavior has been previously observed for argon diffusion in feldspar (e.g., Berger and York, 1981; Harrison and McDougall, 1982), which is analogous to the complex ³He diffusion behavior observed in some quartz samples. Lovera et al. (1989, 1991) proposed a multi-diffusion domain (MDD) model to account for complex argon diffusion behavior, which describes the simultaneous diffusion of discrete, non-interacting intracrystalline sub-domains (e.g., sub-grain fragments) characterized by different effective diffusion length scales. Tremblay et al. (2014b) applied the MDD model framework to ³He diffusion in quartz for samples that exhibited complex Arrhenius behavior, and we have adopted the same approach here.

However, it remains an open question as to whether MDD-type models are applicable to quartz ³He paleothermometry. In Antarctica (Pensacola Mountains), both a single-diffusion domain model using diffusion kinetics from Shuster and Farley (2005) and a two-domain model using kinetics from four local erratics could successfully explain the ³He signal observed in a series of Holocene samples (Tremblay et al., 2014a; Balco et al., 2016), with a similar predicted ³He concentration evolution between the two approaches over this timescale (Balco et al., 2016). However, each approach could only partially explain the ³He signal recorded in samples with older ¹⁰Be exposure ages, with complex exposure history and/or significant inter-sample variability in diffusion kinetics (e.g., different quartz sources for the sandstone lithology) likely acting as compounding factors (Balco et al., 2016). Additional quartz ³He analyses using a MDD model and sample-specific diffusion kinetics were recently conducted on moraine boulders from the Gesso Valley in the Italian Alps with LGM to Late Glacial chronologies (Tremblay et al., 2018). PaleoIsoEDTs within the range of their respective modern EDTs were obtained for two out of five samples, with no clear trend between paleoIsoEDTs and boulder (¹⁰Be) exposure ages or relative moraine age in addition to significant intra-moraine variability. Tremblay et al. (2018) highlighted multiple sources of potential uncertainties related to local shading effects (i.e., vegetation, snow cover, topography), grain-size scaling, and complex boulder exposure histories, which could have contributed to the observed ³He signal inconsistencies.

In this study, bedrock-surface samples were purposefully collected along high-elevation valley profiles progressively deglaciated between the LGM and Holocene, with the aim of limiting the potential for complex exposure (see Sect. 5.2). Diffusion kinetics parameters were measured on one representative sample per profile (MBTP18-9 and GELM18-1). Although inter-sample diffusion kinetics variability cannot be excluded, the apparent homogeneous igneous lithology along each profile supports the representativeness of our chosen sample per profile for diffusion kinetics experiments. Based on this first-order assumption, we noted different ³He diffusion trends between MBTP and GELM representative samples. MBTP quartz exhibits a nearly simple (i.e., linear; Fig. 3a) Arrhenius diffusion behavior, and measured ³He concentrations recorded along the MBTP profile can potentially be interpreted as being at quasi-equilibrium with respect to modern EDTs (despite a slight trend towards a colder signal) when considering the potential sources of uncertainty (e.g., Holocene EDT, $P_{{}^{\mathrm{3}}\mathrm{He}}$, Sect. 5.1 and 5.2; based on the Holocene-calibrated diffusion kinetics). On the contrary, GELM quartz is characterized by complex ³He diffusion behavior (Fig. 3b), and bedrock surfaces record a ³He thermal signal that is apparently well colder than their modern EDTs when using diffusion kinetics derived from a MDD framework (calibrated on Holocene samples). This apparent divergence cannot be resolved within the multiple sources of geologic uncertainties, nor can it be explained by plausible fluctuations in thermal variables (i.e., mean annual temperature and diurnal and/or annual amplitudes) during the Late Glacial and Holocene time periods. One possible interpretation of these results would be that the MDD model we applied to the GELM samples does not accurately represent the ³He diffusion in quartz that occurred during exposure time. This could be because the MDD model does not adequately represent the physical process of ³He diffusion in quartz. From a mineralogical perspective, it is indeed unclear if potential processes involved in the formation of sub(-grain) domains (e.g., cooling, alteration, deformation) are consistent with the assumed conditions of the MDD model, i.e., disconnected sub-domains with fixed volumes, Fickian and isotropic diffusion, and zero concentration boundary conditions (e.g., Lovera et al., 1991; Baxter, 2010). While MDD models have been successfully applied in a number of thermochronology applications (Reiners et al., 2005 and references therein), deformation processes may also lead to interconnected sub-grain microstructures (e.g., Reddy et al., 1999), in which case the MDD model may be inappropriate for obtaining accurate thermal constraints, as already acknowledged in the literature (e.g., Lovera et al., 2002; Harrison and Lovera, 2013). On the other hand, alternative diffusion models involving multi-path diffusion (e.g., Lee, 1995) also suffer from substantial theoretical and experimental gaps (Baxter, 2010; Harrison and Lovera, 2013).

Alternatively, we cannot rule out that an MDD model for quartz ³He paleothermometry (Tremblay et al., 2014b) is applicable on both MBTP and GELM quartz but that the diffusion kinetics is inaccurately constrained. The MDD models we implemented do not provide unique solutions to our laboratory-measured diffusion kinetics, which we then extrapolate down to Earth surface temperatures (< 30 ^∘C). This is illustrated by the significant difference between modern EDTs and estimated paleoIsoEDTs observed for Holocene samples (both MBTP and GELM sites) when using laboratory-derived diffusion kinetics without Holocene calibration (Fig. S5), which therefore supports the additional Holocene calibration step applied in this study.

5.3 Potential role of permafrost processes

At last, we must consider alternative environmental factors besides paleoclimate air temperature variation that may influence ground surface conditions and explain the apparently ³He colder signals recorded in our samples. We compiled all quartz ³He paleoIsoEDTs available in the European Alps (Tremblay et al., 2018; Guralnik et al., 2018; this study; Fig. 10). This compilation, which includes samples with exposure ages ranging from the LGM to the Holocene, reveals no apparent relationship between ³He paleoIsoEDT and ¹⁰Be exposure age (Fig. 10a). However, we do observe a negative correlation between sample paleoIsoEDT and elevation (Fig. 10b). Furthermore, while samples at low to moderate elevations have paleoIsoEDTs that are relatively consistent with their estimated modern EDTs along an apparent linear lapse rate (around −0.5 ^∘C per 100 m lapse rate), paleoIsoEDTs recorded in rock surfaces above ∼ 2200 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ clearly depart from modern EDTs and lapse rate trends, with significantly “colder” ³He signals. Although the compiled Alpine dataset is still limited, such an observed distribution raises the question of the influence of rock-surface elevation on ³He signal records. One hypothesis is that the recorded “colder” ³He signals in high-elevation samples may reflect recent changes in Alpine permafrost ground conditions. Indeed, bedrock-surface samples around or above ∼ 2200 m are located close to or in the lower range of sporadic to discontinuous permafrost distribution in the present-day Alps (Boeckli et al., 2012b; Magnin et al., 2015a). Recent warming after the Little Ice Age is expected to have led to permafrost degradation and restriction of its spatial distribution towards higher elevations (Magnin et al., 2015a, 2017). We hence cannot exclude that those high-elevation bedrock surfaces may have experienced permanent permafrost conditions until recently (i.e., last tens to hundreds of years), where the past MARSTs were thus lower (sub-zero range) than modern MARST estimates scaled on mean annual air temperature (Table 1; Sect. 3.1). In that case, the recent change in climate conditions over the last decades to centuries would have resulted in both mean annual temperature increases and the amplification of annual and diurnal temperature oscillations at the sampling sites greater than those constrained from air temperature records (Etzelmüller et al., 2020) due to the transition from a permafrost to a non-permafrost zone. This scenario would also be consistent with the apparent positive trend observed between ${}^{\mathrm{3}}\mathrm{He}{/}^{\mathrm{10}}\mathrm{Be}$ exposure age ratios and elevation, especially for the GELM site (Fig. 2). This apparent positive trend contrasts with the inverse relationship expected from rock surfaces with a temperature history dominantly controlled by post-LGM ice lowering and general atmospheric warming (i.e., ${}^{\mathrm{3}}\mathrm{He}{/}^{\mathrm{10}}\mathrm{Be}$ exposure age ratio decreasing with elevation). To test the hypothesis about recent permafrost degradation effects would, however, require further quartz ³He measurements at high-elevations and in other alpine/cold regions.

Figure 10 (a) Relationship between ³He paleoIsoEDT anomaly and ¹⁰Be exposure age from available data from moraine boulders and glacially scoured bedrock surfaces in the Alps (results from this study; Tremblay et al., 2018, Gesso, and Guralnik et al., 2018, Gotthard). (b) Relationship between ³He paleoIsoEDT and elevation for same dataset as in panel (a).

6 Conclusion

Paleoglacier fluctuations in alpine settings lack direct constraints of associated past temperature and/or precipitation conditions, which are essential to improving our understanding of the response of glaciers and (para)glacial processes to past and future climate changes. In this study, we applied quartz ³He cosmogenic paleothermometry to derive in situ paleo-temperature (EDT) estimates along two deglaciation sequences gradually exposed from the Last Glacial Maximum to the Holocene in the western and northern European Alps (Mont Blanc and Aar massifs, MBTP and GELM, respectively). Investigation of quartz ³He diffusion kinetics indicates a clear difference between the two study sites, with quasi-linear vs. complex diffusion behaviors for MBTP and GELM sites, respectively. Based on the assumption that the same diffusion kinetics parameters apply to all samples at each site, forward numerical simulations of ³He production and diffusion suggest that no thermal signal from the Late Glacial period should be preserved in investigated rock surfaces with brief exposure durations (several kyr) before the transition to relatively stable Holocene climatic conditions like the present day. However, all our rock-surface samples exposed prior to Holocene indicate an apparent ³He thermal signal significantly colder than present-day conditions. Our recorded ³He signals cannot be explained by realistic post-LGM mean annual temperature evolution in the European Alps (as recorded by other paleoclimatic proxies), nor by changes in annual and/or diurnal temperature oscillations at the study sites.

When accounting for potential uncertainties related to Holocene thermal conditions and the quartz ³He production rate, the ³He signals (ΔEDT) recorded along the MBTP site can potentially be interpreted to be close to equilibrium with present-day and Holocene conditions, with minor change in mean annual temperature or diurnal and annual temperature oscillations. However, ³He-derived paleo-EDTs along the GELM site remain distinctively colder than present-day conditions. One hypothesis is that the multi-diffusion domain models applied to characterize the observed complex diffusion behavior in the GELM quartz do not accurately quantify quartz ³He diffusion for the samples of this site throughout their exposure histories. Alternatively, if the generally colder trend recorded along both profiles is possible, the assumed quartz ³He diffusion kinetics may inaccurately extrapolate to Earth surface temperatures, precluding quantitative EDT constraints from the observed ³He abundances in these samples. Finally, considering the high elevations of the investigated rock-surface samples (> 2000 m), it is also possible that our ³He thermal signals result from much more recent changes in Alpine permafrost ground conditions during the past decades and centuries. While the data presented in this study demonstrate the promising use of ³He cosmogenic paleothermometry to quantify past environmental changes, additional ³He analyses in high-alpine or cold settings would be necessary to clarify to which phenomena the ³He thermal signal is most responsive, i.e., between Late Pleistocene ambient temperature variations and recent changes in permafrost distribution.