Glacier response to Holocene warmth inferred from in situ ¹⁰Be and ¹⁴C bedrock analyses in Steingletscher's forefield (central Swiss Alps)

Mid-latitude mountain glaciers are sensitive to local summer temperature changes. Chronologies of past glacier fluctuations based on the investigation of glacial landforms therefore allow for a better understanding of natural climate variability at local scale, which is relevant for the assessment of the ongoing anthropogenic climate warming. In this study, we focus on the Holocene, the current interglacial of the last 11 700 years, which remains a matter of dispute regarding its temperature evolution and underlying driving mechanisms. In particular, the nature and significance of the transition from the early to mid-Holocene and of the Holocene Thermal Maximum (HTM) are still debated. Here, we apply an emerging approach by combining in situ cosmogenic ¹⁰Be moraine and ¹⁰Be–¹⁴C bedrock dating from the same site, the forefield of Steingletscher (European Alps), and reconstruct the glacier's millennial recession and advance periods. The results suggest that, subsequent to the final deglaciation at ∼10 ka, the glacier was similar to or smaller than its 2000 CE extent for ∼7 kyr. At ∼3 ka, Steingletscher advanced to an extent slightly outside the maximum Little Ice Age (LIA) position and until the 19th century experienced sizes that were mainly confined between the LIA and 2000 CE extents. These findings agree with existing Holocene glacier chronologies and proxy records of summer temperatures in the Alps, suggesting that glaciers throughout the region were similar to or even smaller than their 2000 CE extent for most of the early and mid-Holocene. Although glaciers in the Alps are currently far from equilibrium with the accelerating anthropogenic warming, thus hindering a simple comparison of summer temperatures associated with modern and paleo-glacier sizes, our findings imply that the summer temperatures during most of the Holocene, including the HTM, were similar to those at the end of the 20th century. Further investigations are necessary to refine the magnitude of warming and the potential HTM seasonality.

1 Introduction

Mountain glaciers in most glacierized regions of the world, such as the European Alps, are currently rapidly retreating in response to accelerating global warming, driven by human-induced greenhouse gas emissions into the atmosphere (IPCC 2007, 2013, 2021). Small mountain glaciers are reliable indicators of regional climate variations on decadal to multi-millennial timescales, because their mass balance is sensitive to variations of meteorological parameters, in particular summer temperature and precipitation (Oerlemans, 2005). Investigating past glacier behavior and the underlying regional climate variability provides the opportunity to better understand the natural driving mechanisms within Earth's climate system and to help quantify the anthropogenic contribution to the ongoing climate evolution (e.g., Roe et al., 2021).

The current interglacial Holocene followed the end of the last glacial period ∼11 700 years ago and is characterized by moderate climate variations, including both colder-than-today and warmer phases (Mayewski et al., 2004; Wanner et al., 2008). In the northern mid- and high-latitudes, many studies provide evidence of several millennia of warm conditions during the early and mid-Holocene, generally referred to as the Holocene Thermal Maximum (HTM) (e.g., Renssen et al., 2009; Axford et al., 2013; Heiri et al., 2015; Kobashi et al., 2017). However, the occurrence of extended periods that were significantly warmer than recent decades are still debated (e.g., Marcott et al., 2013; Marsicek et al., 2018; Affolter et al., 2019; Kaufman et al., 2020; Bova et al., 2021). The response of mountain glaciers to these Holocene warm periods remains unclear, because records of when and how long mountain glaciers have receded to modern extents or beyond are still scarce and challenging, because much of the potential evidence is buried beneath ice.

The European Alps are one of the regions that are best documented in terms of Holocene glacier behavior (Ivy-Ochs et al., 2009; Solomina et al., 2015), but existing Holocene glacial chronologies are dominated by studies of moraines and thus large glacier extents that occurred during cold episodes (e.g., Schimmelpfennig et al., 2012, 2014; Moran et al., 2017; Le Roy et al., 2017; Protin et al., 2019, 2021; Braumann et al., 2020, 2021). Most of the existing constraints on the timing and amplitudes of glacier recessions come from discrete radiocarbon dates of sub-fossil wood and peat (e.g., Porter and Orombelli, 1985; Baroni and Orombelli, 1996; Nicolussi and Patzelt, 2000; Hormes et al., 2001, 2006; Deline and Orombelli, 2005; Joerin et al., 2008; Nicolussi and Schlüchter, 2012; Le Roy et al., 2015), and few studies provide records that characterize glacier extents during the majority of the Holocene, including periods of glacier retreat (Joerin et al., 2006; Luetscher et al., 2011; Badino et al., 2018).

A more recently developed and powerful approach to addressing the chronological reconstruction of millennial-scale Holocene glacier retreat relies on the measurement of in situ cosmogenic ¹⁴C exposure dating in deglaciated bedrock, combined with in situ cosmogenic ¹⁰Be or other dating techniques, so far applied in only a few studies around the globe (Goehring et al., 2011, and Wirsig et al., 2016, in the Alps; Schweinsberg et al., 2018, Pendleton et al., 2019, and Young et al., 2021, in the Greenland/Baffin region; Rand and Goehring, 2019, in Norway; Johnson et al., 2019, in Antarctica). This method provides the possibility to quantitatively derive the total duration that deglaciated bedrock has been exposed, i.e., was ice-free, and how long it was buried beneath ice throughout the Holocene. The pioneering studies applying in situ ¹⁴C–¹⁰Be exposure–burial dating showed that Rhône Glacier, located in the central Swiss Alps, was smaller than its ∼2005 CE extent for 6.4±0.5 kyr, i.e., for the majority of the Holocene (Goehring et al., 2011, 2013). The challenge of this approach arises from the need for additional chronological constraints to determine the specific number and timing of recession periods, if they were interrupted by successive glacier advances.

Figure 1 Maps of the study site based on shaded relief ALTI3D models by ©swisstopo. (a) Switzerland and the location of Steingletscher as a red dot (central Swiss Alps, 47 ^∘C). (b) Overview of the whole glacier catchment with the extents of Steingletscher and Steinlimigletscher in the year 2016. The red rectangular corresponds to panel (c). (c) Steingletscher's forefield with mapped Holocene moraines, their ¹⁰Be exposure ages and 1σ analytical uncertainties (white boxes; recalculated from Schimmelpfennig et al., 2014; one outlier in italic), and the new bedrock sample locations with their apparent ¹⁰Be and ¹⁴C exposure durations and 1σ analytical uncertainties (green boxes). Mean landform ages are shown in darker boxes with 1σ uncertainties including analytical and ¹⁰Be production rate uncertainties. Pink boxes give years of historically recorded moraine deposits.

In this study, we focus on Steingletscher, an often-visited, small mountain glacier in the central Swiss Alps (Fig. 1). Its evolution has been well-monitored since the end of the 19th century, showing that it has been constantly retreating, losing >1000 m of its length since 1985 CE as a response to the modern climate warming (GLAMOS, 2020). Its forefield has been subject to various scientific studies, including investigations of the glacier's responses to Holocene cold episodes (King, 1974; Schimmelpfennig et al., 2014). Schimmelpfennig et al. (2014) mapped and ¹⁰Be-dated the Holocene moraines in the forefield of Steingletscher (Fig. 1), providing evidence of several large glacier extents between the early Holocene and the end of the Little Ice Age (LIA, ∼14th to 19th century). This existing knowledge on the glacier's advances and its high sensitivity to warming make it an ideal target to investigate the more difficult question: how did this glacier respond to extended warm periods during the Holocene? We therefore present here new measurements of in situ cosmogenic ¹⁴C and ¹⁰Be in bedrock recently deglaciated in front of Steingletscher, generally following the approach applied at nearby Rhône Glacier (Goehring et al., 2011). We combine the new data on glacier recession with the previously published local Holocene moraine chronology, as well as with earlier published bracketing radiocarbon ages (King, 1974; Hormes et al., 2006) and historical and instrumental documentation of the recent glacier evolution. Our principal objectives are to (1) temporally constrain the Holocene intervals during which Steingletscher was at least as retracted as in modern times and (2) evaluate whether the Steingletscher and Rhône Glacier retreat histories are individual records or rather represent regional glacier responses to warming climate phases. We then put the result into the context of Holocene climate and glacier evolution in the Alps to test the significance of the HTM at regional scale.

2 Study site, previous chronological work and sampling strategy

Steingletscher is located in the eastern part of the Bernese Alps (∼47^∘ N, 8^∘ E) at an altitude of 2220 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ close to Susten Pass and is part of a larger glacier catchment that hosts another glacier, Steinlimigletscher (Fig. 1b). It had a length of 3.35 km in 2019 CE and an area of 7.6 km² in 2013 (GLAMOS, 2019, 2020). Metagranitoids, gneisses and amphibolites constitute the geology of the glacier's surroundings. During the reference period 1991–2020 CE, the climate at the nearest weather station (Meiringen, 588 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) was characterized by an annual mean temperature of 8.7 ^∘C, a monthly mean temperature range between −1.1 and 17.9 ^∘C, and an annual mean precipitation of 1341 mm (https://www.meteoswiss.admin.ch, last access: 10 January 2022). During the reference period 1981–2010 CE, snow cover with a thickness of >50 cm occurred on 2.9 d (https://www.meteoswiss.admin.ch, last access: 18 July 2021; climsheet 2.1.6/5 January 2021). In the forefield of Steingletscher, the mean July temperature was 9.5 ^∘C between 1991 and 2020 CE (information provided by costumer service MeteoSwiss).

Figure 2 Photographs taken during field work. (a, b) Bedrock surfaces sampled close to the summit of Hublen Plateau. (c) Bockberg riegel and the glacier terminus in the year 2016. (d) Location of sample STEI-16-10 on Chüebergli riegel with the view on the eastern side of the Steingletscher forefield, highlighting the LIA trimline and composite moraine.

Steingletscher's forefield, stretching almost linearly towards the north, features glacially smoothed hills, moraines, trimlines and the proglacial Steinsee, a lake located of the half-bowl-shaped distal part of the forefield. The catchment's outlet is located in the northwestern corner of this bowl and drains westward into Gadmen valley. On the left-lateral catchment flank, the up to 2090 m high “Plateau Hublen” and the lower “Plateau In Miseren” are characterized by glacially polished and lichen-covered bedrock knobs (roches moutonnées; Fig. 2a and b), interspersed with vegetated depressions, and small peat bogs and lakes (Fig. 1c). This landscape is overprinted by relicts of several moraine belts, i.e. the outer and the inner Hublen moraines and the In Miseren moraine, which were ¹⁰Be-dated at ∼11.0, ∼10.6 and ∼10.0 ka, from outer to inner (Schimmelpfennig et al., 2014). King (1974) obtained minimum radiocarbon ages for moraine formations from basal parts of peat bogs on Plateaus Hublen and In Miseren, which are in agreement with the ¹⁰Be Holocene moraine data (Fig. 3a). They indicate that cold conditions still persisted during the deglaciation after the Younger Dryas (YD; 12.8–11.7 ka; Rasmussen et al., 2006). However, further evidence of YD-related extents has not yet been identified. Here, we targeted two roches moutonnées located close to the highest point of Plateau Hublen, a few meters outboard of the outer Hublen moraine (Figs. 1 and 2a, b) with the objective to date the timing of initial deglaciation of this plateau.

Figure 3 Map of Holocene extents of Steingletscher based on ¹⁰Be moraine dating (see Fig. 1c) for the period between ∼11 ka and LIA, and on historical topographic maps and data from the Swiss Glacier Inventory (CH-INVGLAZ), adapted from Wirz (2007) and King (1974), for the period between 1850 and 1999 CE (a), on shaded relief ALTI3D model by ©swisstopo. Black stars represent locations of radiocarbon-dated organic material and corresponding calibrated ages from King (1974) (in blue boxes; black numbers are maximum ages and white numbers are minimum ages for moraine deposits or glacier retreats/advances; yellow ages correspond to certain pollen assemblages) and from Hormes et al. (2006) (in white box). Bedrock sample locations on Chüebergli and Bockberg riegels are shown as green dots, and samples STEI-16-10 and -12 are highlighted with blue rims. (b) Length measurements of Steingletscher since 1893 (from GLAMOS, 2020). Red and blue bars correspond to trends of glacier retreat and advance, respectively; intervals of glacier stagnation are in white.

The absence of moraines after ∼10 ka and throughout the mid-Holocene indicates a warmer climate at that time (Schimmelpfennig et al., 2014). This is supported by radiocarbon ages of ∼9 ka cal BP and from the mid-Holocene obtained from peat bogs on Plateau In Miseren (King, 1974; Table 2 in Schimmelpfennig et al., 2014; Fig. 3a). In addition, evidence of a substantially retracted mid-Holocene extent of Steingletscher comes from two wood fragments melted out from the glacier front between 1995 and 2000 CE and radiocarbon-dated at 5.3–4.8 and 4.8–4.6 ka cal BP (Hormes et al., 2006; Fig. 3a). In the same study, two organic silt samples from the forefield of the neighboring Steinlimigletscher, also collected between 1995 and 2000 CE, were radiocarbon-dated at 5.9–5.3 and 2.3–1.8 ka cal BP. Further information on the amplitude and duration of glacier recession during the early and mid-Holocene is missing, as geomorphic markers of glacier extents during that time were destroyed by the glacier re-advances during late Holocene cooling.

Evidence of late Holocene moraine formation comes from glacial boulders in the vicinity of the catchment outlet dated at ∼2.9 ka (Schimmelpfennig et al., 2014). This maximum glacier extent during the late Holocene is corroborated by radiocarbon dates of organic material from soil and peat on the right-lateral side of the catchment outlet, which provide bracketing ages for the glacier advance and retreat around 3 ka ago (King, 1974; Table 2 in Schimmelpfennig et al., 2014; Fig. 3a). One large moraine boulder located on a ridge inboard of the ∼3 ka moraine was ¹⁰Be-dated at ∼1.9 ka and might represent a glacier extent at that time (Schimmelpfennig et al., 2014).

Figure 4 Aerial photograph of Steingletscher terminus position in 1982 CE with indications of the timing of earlier glacier extents, from http://www.swisseduc.ch (last access: 10 January 2022). The location of the Chüebergli sample profile is indicated by the green box. Ice thickness in 1982 CE is estimated at ∼50 to ∼20 m above the analyzed samples.

The most evident geomorphic markers of glacier expansion in Steingletscher's forefield are those from the Little Ice Age, including a sharp composite moraine on the eastern side of the forefield (Fig. 2d), multiple moraine ridges, and a clearly visible trimline (Figs. 2d and 4). Boulders from the preserved moraines yield ¹⁰Be ages between ∼570 and 140 years, consistent with the period of the LIA (Schimmelpfennig et al., 2014).

Historical and instrumental records provide constraints on glacial extents during the general retreat between 1850 CE and the beginning of the 21st century, which are here mainly based on previous compilations by King (1974) and Wirz (2007) and on glacier length measurements (GLAMOS, 2020). Figure 3a shows the glacier outlines in the years 1850, 1920, 1933, 1973, 1988, 1999 and 2007. The length measurements of Steingletscher between the years 1893 and 2019 (Fig. 3b) indicate that the glacier had never retreated as much as in 2007 during that time. The most pronounced retreat of the 20th century occurred between 1960 and 1970, leading very briefly to a minimum extent that was comparable to that of the very beginning of the 21st century, but still slightly bigger than that in 2007 (Fig. 3b). During the subsequent advance of the 1970–1980s the glacier reached as far as into Steinsee. Since 1988, the glacier has been rapidly and continuously retreating and has lost >1000 m of its length in response to the modern global warming. Steingletscher is thus currently not in equilibrium with climate but lags behind the accelerating warming by up to several decades.

To investigate how Steingletscher responded to naturally driven Holocene warmth, specifically how long it was similar in size to its modern configuration during the Holocene, we targeted two bedrock riegels for in situ ¹⁴C–¹⁰Be exposure–burial dating. One riegel is located east of Chüebergli (“Chüebergli riegel” hereafter) and is ∼400 m long. The other riegel is located north of Bockberg (“Bockberg riegel” hereafter) and is ∼200 m long. Both riegels are characterized by glacially polished, recently deglaciated roches moutonnées (Fig. 2c and d) that form steep cliffs towards the north. Chüebergli riegel was completely covered under at least ∼140 m of ice during the LIA maximum, inferred from the altitude of the LIA composite moraine and trimline on the eastern catchment flank, and continued to be buried throughout most of the 20th century. Eight samples were collected on a transect on Chüebergli riegel from the highest and outmost bedrock surface down to the lowest bedrock surface in the glacial trough, following the sampling approach in Goehring et al. (2011). Note that the three lowest sample surfaces were still covered by ice in 2007 (Fig. 1c) but were ice-free by the sampling year 2010. From Bockberg riegel, two samples were analyzed (Fig. 1c).

After the glacier retreat in the 1960s, the outmost parts of the two riegels were temporarily ice-free, but ice-covered again in the 1970s–1980s by ∼20 m of ice (Figs. 3a and 4). Chüebergli riegel might have been completely ice-free briefly around 1970. Both riegels were deglaciated in the early 21st century, and Steingletscher has continued to retreat since.

3 Methodology

3.1 Fieldwork

The 14 bedrock samples were collected during field campaigns in 2010 and 2016, targeting glacially polished and striated surfaces that were free of sediment cover. To minimize the risk of significant snow and sediment cover in the past, slightly sloping surfaces were preferred (Fig. 2d). Rock surface pieces with average thicknesses of ∼2 to ∼4 cm (Table 1) were sampled using either hammer and chisel alone or in combination with a cordless angle grinder and a diamond blade. Primarily the quartz-rich parts of the gneissic lithologies were targeted. Latitude, longitude and elevation at the sample locations were recorded with a Trimble GeoTX GPS; the data reduction was conducted using the WGS coordinate system relative to EGM896 Geoid (Table 1). Correction factors for the shielding by the surrounding topography and the strike and dip of the sampled surfaces were determined from measurements using a handheld inclinometer. These correction factors range from 0.87 to 1.0 (Table 1).

Table 1 Sample field information and analytical ¹⁰Be data. Carrier solution “CEREGE” used at LN2C has a ⁹Be concentration of 3025 µg g⁻¹; ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{9}}\mathrm{Be}$ ratios measured at ASTER (Arnold et al., 2010) were normalized to in-house standard STD-11 with the ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{9}}\mathrm{Be}$ ratio of $\mathrm{1.191}(\pm \mathrm{0.013})\times {\mathrm{10}}^{-\mathrm{11}}$ (Braucher et al., 2015). LDEO carrier solutions “5.1” and “5.2” have ⁹Be concentrations of 1024 and 1031.78 µg g⁻¹, respectively; ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{9}}\mathrm{Be}$ ratios measured at CAMS were normalized to standard 07KNSTD with the ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{9}}\mathrm{Be}$ ratio $\mathrm{2.85}\times {\mathrm{10}}^{-\mathrm{12}}$ (Nishiizumi et al., 2007). A ¹⁰Be half-life of $\mathrm{1.387}(\pm \mathrm{0.01})\times {\mathrm{10}}^{\mathrm{6}}$ years was used (Chmeleff et al., 2010; Korschinek et al., 2010). Samples were corrected for chemistry blanks by subtracting the number of atoms ¹⁰Be in the corresponding blank from that of the sample.

3.2 Analytical methods

We separated and decontaminated quartz from 14 samples and extracted ¹⁰Be from the clean quartz after spiking it with a pure ⁹Be carrier (Table 1). Chemical processing was carried out at the Lamont-Doherty Earth Observatory (LDEO) Cosmogenic Nuclide Laboratory (New York, USA) according to the standard procedures described in Schaefer et al. (2009) and at the Laboratoire National des Nucléides Cosmogéniques (LN2C) at the Centre Européen de recherche et d'Enseignement des Géosciences de l'Environnement (CEREGE, Aix en Provence, France) following routine methods described for example in Protin et al. (2019). ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{9}}\mathrm{Be}$ ratios were measured at the Lawrence Livermore National Laboratory – Center for Accelerator Mass Spectrometry (LLNL-CAMS) and at Accélérateur pour les Sciences de la Terre, Environnement, Risques (ASTER) at CEREGE. All data related to the ¹⁰Be analyses are presented in Table 1.

In situ ¹⁴C was extracted from the quartz of the two bedrock samples from Chüebergli riegel with the highest ¹⁰Be concentrations (STEI-16-10 and -12). These extractions were performed at LDEO following the procedure described in Goehring et al. (2014) and Lamp et al. (2019), with two updates: the purified CO₂ sample and blank gas were diluted with only small amounts of ¹⁴C-free gas corresponding to ∼20 µg of C (Table 2); the sample and dilution gas mixtures were not converted into graphite, but sealed into pyrex break seals for ${}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}$ ratio measurements at the AixMICADAS facility at CEREGE using the ion source dedicated for gaseous samples (Bard et al., 2015; Tuna et al., 2018), thus avoiding the graphitization step. ¹⁴C concentrations were calculated following the method of Hippe and Lifton (2014) (Table 2).

Table 2 Analytical in situ ¹⁴C data. “Mass C (µg) sample” and “Mass C (µg) sample + diluted” are the masses of carbon released during extraction and after addition of ¹⁴C-free dilution gas, respectively. Oxalic acid OX-II is used as AMS standard, and all measurements are corrected for the AMS machine background blank with a ${}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}$ ratio of 0.0023. The percent of modern carbon value (pmC) includes the d¹³C correction, which is undone following the calculations of Hippe and Lifton (2014). Samples were corrected for the procedural blank by subtracting the number of atoms ¹⁴C in the blank from that of the samples. The apparent ¹⁴C ages do not account for subglacial erosion effects and are therefore minimum ages, referenced to the sampling year 2016. See text for parameters used for the age calculations.

3.3 Principles of the exposure–burial dating approach

Three types of glacial surfaces with different exposure histories are investigated in this study.

• i.

  Moraine boulders. They can in most cases be assumed to have a simple exposure history; i.e., they were free from cosmogenic nuclides at the moment of their stabilization after the glacier had retreated from its advance and have since been continuously exposed. In this case, the analysis of ¹⁰Be alone usually provides the exposure age of the surface. Based on the consistency of the ¹⁰Be moraine boulder ages of Schimmelpfennig et al. (2014), they seem to fulfill this condition.

• ii.

  Bedrock surfaces that remained continuously ice-free during the Holocene. They are located outboard of the maximum Holocene glacier extent. We assume that they were covered long enough and subglacially eroded deep enough during the ∼100 ka lasting last glacial period that ended with the YD that these surfaces were free from cosmogenic nuclides at the moment of their last deglaciation. Subsequently, they experienced a simple exposure history, if cover by sediment, soil or vegetation is negligible. The two dated bedrock samples from Hublen (STEI-16-1 and -2), located outboard of the outer Hublen moraine, are assumed to fulfill these conditions.

• iii.

  Bedrock surfaces that were alternately ice-free and ice-covered during the Holocene. They are located inboard of all Holocene moraines. Like the bedrock type ii, their cosmogenic nuclide inventory is assumed to have been set to zero during the last glacial period. If during the Holocene phases of no glacier cover were followed by glacier cover during 10²–10³-year periods and with moderate subglacial erosion, cosmogenic nuclide concentrations accumulated from several exposure periods but were reduced through the glacial erosion during ice cover. Consequently, the analysis of ¹⁰Be alone in this type of sample provides a minimum duration of the cumulative exposure period. The 10 samples from Chüebergli and Bockberg riegels correspond to this type of bedrock.

In the case of the bedrock type iii, the combined analysis of ¹⁰Be and ¹⁴C allows for solving for the unknown exposure duration and erosion depth, if the moment of initial deglaciation can be constrained (see Sect. 3.5). This is because ¹⁴C (half-life 5.7 kyr) decays much faster than ¹⁰Be (half-life 1.39 Myr, Chmeleff et al., 2010; Korschinek et al., 2010), and therefore their concentrations evolve differently as a function of exposure and burial: during exposure both nuclides accumulate at a nearly constant rate with the ${}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{10}}\mathrm{Be}$ concentration ratio close to ∼3; during burial under ice, production stops and only the decay of the ¹⁴C concentration is notable on these relatively short timescales, leading to lower concentration ratios (Fig. 5; Hippe, 2017). In addition, subglacial erosion of the initially exposed surface removes the superficial bedrock layers, thereby advecting less ¹⁴C and ¹⁰Be concentrated rock from depth to the surface (Fig. 5). If this subglacial erosion is negligible or largely dominated by abrasion, the ${}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{10}}\mathrm{Be}$ concentration ratio is less than 3 in the collected samples (Fig. 5), while higher ratios indicate surface quarrying by the glacier (Rand and Goehring, 2019). The latter is because production at depth of ¹⁴C is higher relative to that of ¹⁰Be due to a higher ¹⁴C contribution from muons, which penetrate deeper under the subsurface than neutrons (Hippe, 2017). Thus, in the case of moderate erosion rates dominated by abrasion, the apparent (i.e., non-burial- and erosion-corrected) ¹⁰Be and ¹⁴C ages provide both minimum exposure durations, with apparent ¹⁴C ages being younger than their ¹⁰Be counter parts, and the ${}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{10}}\mathrm{Be}$ concentration ratio less than ∼3.

Figure 5 General concept of the in situ ¹⁴C–¹⁰Be exposure–burial approach. Left panels show a hypothetical scenario of alternating large and retracted glacier extents, covering and uncovering a fictive sample location (red spot). Middle panels depict the position of the sample material (red rectangle) with regard to the bedrock surface. Right panels of stages (2) and (3) illustrate the trajectories of the in situ ¹⁴C and ¹⁰Be concentrations in the sample during its exposure and advection (exhumation) towards the surface (Lagrangian perspective; see Sect. 3.5). Stage (1): subsequent to continuous ice cover and deep subglacial erosion during the last glaciation, the bedrock is free from cosmogenic nuclides. Deposition of the early Holocene (EH) moraines and subsequent glacier retreat constrains the beginning of the Holocene glacier retreat, i.e., 10 ka in this example. Stage (2): the sample location is exposed for an unknown duration t_exp, with the sample material remaining at the same depth (no surface erosion), where ¹⁴C and ¹⁰Be accumulate constantly. Stage (3): the sample location is buried by ice, which abrades the bedrock surface, thus progressively exhuming the sample material. No production of ¹⁴C and ¹⁰Be due to shielding by ice cover. Radioactive decay leads to steady loss of ¹⁴C (short half-life of 5.7 kyr) but does not affect the ¹⁰Be concentration due to the long half-life of this nuclide of 1.4 Myr. The burial duration t_b is equal to the timing of initial EH glacier retreat minus t_exp. The dashed line represents the initially exposed surface elevation. Stage (4): exposed since the recent glacier retreat, the sample material is now at the surface where it is collected. E represents the unknown erosion depth, i.e., the thickness of rock layer that was removed by abrasion during stage (3). The lowest right panel shows a two-isotope plot with the thick black line representing continuous exposure (${}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{10}}\mathrm{Be}$ ratio = ∼3) and the grey area burial and subglacial abrasion (${}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{10}}\mathrm{Be}$ ratio <3). The ${}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{10}}\mathrm{Be}$ ratio and ¹⁰Be concentration of the fictive sample is plotted as a red dot.

We make the assumption that during the periods of burial, the ice cover was always thick enough at our sample locations to hinder significant ¹⁴C accumulation via muogenic production. Shielding by >70 m of ice is required to reduce ¹⁴C production to 1 % of its surface production, and under a thin ice cover of ∼13 m ¹⁴C is produced at 10 % compared to an ice-free surface, while ¹⁰Be is produced at only 1 % (Hippe, 2017). The photograph in Fig. 4 allowed us to estimate that during the glacier extent in 1982, ice cover was on the order of ∼20 to 50 m above the sample locations, suggesting that ¹⁴C production in the subglacial rock surfaces during episodes of ice cover should be small enough to not significantly affect the interpretation of our data.

This exposure–burial bedrock dating approach thus allows us to determine the cumulative duration that the glacier retreated beyond the sample locations during the Holocene. As reference for this minimum amplitude of retreat, we refer to the extent in modern times when the glacier uncovered the sample locations for the last time, i.e., ∼1999 CE for sample STEI-16-12 and ∼2007 CE for sample STEI-16-10 (Fig. 3a). As the time elapsed between these two years is insignificant compared to the centennial to millennial timescales investigated here, we simplify and refer to 2000 CE for both sample locations.

3.4 Calculations of simple exposure ages

The ¹⁰Be moraine and apparent bedrock ages discussed below were calculated with CREp (http://crep.crpg.cnrs-nancy.fr/#/, last access: 11 January 2022; Martin et al., 2017), choosing the local “Alpine” ¹⁰Be spallation production rate (Claude et al., 2014) (4.11±0.10 $\mathrm{atoms}{\mathrm{g}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ as calculated in CREp via the link to the ICE-D calibration database), the ERA40 atmosphere model (Uppala et al., 2005) and Lal–Stone time-corrected scaling (Lal, 1991; Stone, 2000; Balco et al., 2008) with atmospheric ¹⁰Be-based VDM (Muscheler et al., 2005; Valet et al., 2005). The apparent ¹⁴C bedrock ages were calculated accounting for spallogenic and muogenic production and radioactive decay, using the global ¹⁴C spallation production rate of 12.22±0.89 $\mathrm{atoms}{\mathrm{g}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$, the muon parameters from Balco (2017) and the half-life of 5730 years. The same atmosphere model and scaling methods as for ¹⁰Be were applied. A density of 2.7 g cm⁻³ is assumed for all samples.

As most of the recent studies in the Alps used the former CRONUS-Earth calculator by Balco et al. (2008) to calculate cosmogenic nuclide ages, we also show the ¹⁰Be moraine and apparent ¹⁰Be and ¹⁴C exposure durations of bedrock calculated with version 3 of this tool, choosing the same parameters as above regarding the ¹⁰Be production rate (calculated as 4.04±0.38 $\mathrm{atoms}{\mathrm{g}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ in the former CRONUS-Earth calculator using the calibration data from the ICE-D calibration database), atmosphere model and scaling scheme. Note that the calculator applies default parameters for the geomagnetic field model and the ¹⁴C production rate.

All ¹⁰Be and ¹⁴C exposure ages and durations are shown in Tables 2 and 3. ¹⁰Be exposure ages and durations calculated with both calculators differ by 1 %–2 % for early Holocene ages, and by at most 8 % for younger ages and exposure durations. ¹⁴C exposure durations differ by 3 %–5 %. Higher differences in the results are due to high late Holocene inter- and intra-variability in the geomagnetic field records used in both calculators.

Table 3 ¹⁰Be moraine ages and apparent ¹⁰Be bedrock ages at Steingletscher, with their full 1σ uncertainties. The uncertainties in parentheses are the analytical 1σ uncertainties in the case of individual ages and standard deviations in the case of mean ages. Moraine ages are recalculated from Schimmelpfennig et al. (2014), using the same methods as for the bedrock samples (see Sect. 3.4). One outlier is in italics. The apparent ¹⁰Be exposure ages do not account for subglacial erosion effects and are therefore minimum ages (with reference to the year of sampling). All moraine ages and bedrock ages from Hublen are referenced to 1950 CE.

Unless otherwise stated, in situ cosmogenic ages corresponding to glacial surface types i and ii are reported in the text, tables and figures with the unit ka (“thousand years ago”) with reference to 1950 CE (before BP), i.e., 60–66 years were deduced from their exposure ages. Exposure and burial durations of type iii bedrock are given with the unit kyr and refer to their year of sampling (2010 or 2016 CE).

3.5 Modeling of complex bedrock exposure history and erosion depth

Regarding the complex exposure history of bedrock type iii, three variables are unknown: the cumulative exposure duration t_exp, the cumulative burial duration t_b and the erosion depth E, while the two-nuclide system allows for two of these unknowns to be solved (Goehring et al., 2011). In our study, the age of the youngest early Holocene moraine (In Miseren moraine, ∼10 ka), provides the timing of initial deglaciation of Steingletscher's forefield, and therefore the cumulative burial duration can be constrained by $t_{\mathrm{b}}=\mathrm{10}\mathrm{ka}-t_{\text{exp}}$ (Fig. 5; Goehring et al., 2011). In the pioneer Rhône Glacier study by Goehring et al. (2011) the remaining two unknowns, t_exp and E, were determined for each sample through the classical Eulerian computation, i.e., using the two equations for ¹⁰Be and ¹⁴C production and loss at the rock surface. In a follow-up study, Goehring et al. (2013) proposed in addition an isochron Bayesian approach that allows several samples to be considered simultaneously with the purpose of reducing the uncertainties for the whole data set.

Here, we use the Lagrangian instead of the Eulerian calculation method (Knudsen et al., 2019), combined with a Monte Carlo–Markov chain (MCMC) inversion to constrain the values of t_exp and E. We define a forward model to predict ¹⁰Be and ¹⁴C concentrations as a function of t_exp and E in each of the two ¹⁰Be- and ¹⁴C-analyzed rock samples (STEI-16-10 and -12) during their advection toward the surface, as schematically illustrated in Fig. 5. For this, we apply the same production parameters and scaling method as for the simple exposure age calculations. We assume that the ¹⁰Be and ¹⁴C concentrations were set to zero due to deep and sustained glacial erosion during the previous glacial phase (Fig. 5). The onset of the exposure history is set at 10 ka, as constrained by In Miseren moraine. The subsequent history is divided into two steps. First, the surface is exposed for a duration t_exp, leading to steady accumulation of ¹⁰Be and ¹⁴C. Second, the surface is covered by ice until the present, leading to zero nuclide production, but bedrock erosion of a thickness E. The corresponding hypothetical ¹⁰Be and ¹⁴C concentration–time trajectories are shown in Fig. 5. We then use a standard MCMC approach to sample the parameter plane defined by t_exp and E with the Metropolis–Hastings algorithm and obtain the posterior distributions of these parameters for both samples. For each sample, we run six MCMC chains of length 10⁵, including a 10³ burn-in phase. For each parameter we report the average and standard deviation of the chain. The average erosion rate ε is subsequently determined by $\mathit{\epsilon }=E/t_{\mathrm{b}}$.

The advantage of using the Lagrangian instead of the Eulerian approach is that it is more flexible when calculating the changing nuclide concentrations in a bedrock sample during its advection to the surface (Knudsen et al., 2019). In our specific case it conveniently allows the depth-dependent production and decay of ¹⁴C to be considered as a function of time, instead of summing up nuclide production and loss in a virtual surface sample. However, we stress that in our application, no significant differences result from the two approaches. Figure 7l indeed depicts the concentration–time trajectories at the surface, i.e., using the Eulerian approach, both for the above-described single exposure–burial scenario and a more complex scenario with several exposure–burial alternations.

3.6 Recalibration of previously published radiocarbon dates

The radiocarbon ages previously published in King (1974) and Hormes et al. (2006) and discussed in this study were calibrated with the online program OxCal 4.4 (Bronk Ramsey, 2009), using its standard options and the IntCal20 calibration curve (Reimer et al., 2020), and are reported relative to the year 1950 CE. This changes the 2σ age intervals by at most 6 % compared to the calibration with the earlier OxCal 3.9 version presented in Hormes et al. (2006). The radiocarbon ages published in King (1974) were uncalibrated and can therefore not be compared to the ages calibrated in our study. Note that with regard to the original studies, the general interpretations of all radiocarbon ages discussed here remain unaffected from the (re)calibration.

4 Results

Tables 2 and 3 and Figs. 1c and 6a show the exposure ages and apparent exposure durations directly calculated from the measured ¹⁰Be and ¹⁴C concentrations. Also listed in Table 3 are all ¹⁰Be moraine boulder ages and moraine mean ages previously published in Schimmelpfennig et al. (2014) and recalculated with the methods presented in Sect. 3.4. In the tables, ages are shown with their full uncertainties (i.e., including analytical and production rate uncertainties) and analytical uncertainties only. In the text and in the figures, the individual ages are presented with their analytical uncertainties only, while moraine mean ages are presented in the text and in the figures with their full uncertainties, i.e., standard deviation and production rate error combined through simple error propagation (square root of the sum of their values in quadrature).

Figure 6 Exposure durations of bedrock samples inferred from measured in situ ¹⁰Be and ¹⁴C concentrations. (a) Apparent exposure durations calculated from ¹⁰Be concentrations of eight samples and ¹⁴C concentrations of two samples as a function of position in the glacial trough profile. The subglacial erosion effect during ice burial is not accounted for; calculations therefore yield minimum exposure durations. (b) Exposure durations and erosion depths modeled from combined in situ ¹⁰Be and ¹⁴C concentrations of samples STEI-16-10 and -12 (orange and green signatures). Solid grey curves represent possible exposure duration and erosion–depth combinations for samples with ¹⁰Be concentration only. Yellow lines indicate selected theoretical bedrock erosion rates.

The two ¹⁰Be bedrock samples from Plateau Hublen yield indistinguishable ages of 11.39±0.58 ka (STEI-16-1) and 11.52±0.56 ka (STEI-16-2) with an arithmetic mean age of 11.5±0.3 ka.

The two oldest apparent ¹⁰Be bedrock exposure durations from the Chüebergli riegel profile are 5.90±0.10 kyr (STEI-16-12) and 4.10±0.07 kyr (STEI-16-10). The other samples on this profile range between 2.67±0.07 kyr (STEI-5) and 0.05±0.01 kyr (STEI-2). All eight apparent ¹⁰Be ages from the profile present a general trend from high values at the outer margin towards low values at the inner, lowest sample locations (Fig. 6a). This trend is consistent with the higher ice flow velocity in the center of a glacial trough that leads to deeper subglacial bedrock erosion during periods of ice cover, thus reducing the ¹⁰Be surface concentrations at a higher rate (Goehring et al., 2011).

Only the two samples with the oldest apparent ¹⁰Be exposure durations (STEI-16-12 and STEI-16-10) were chosen for ¹⁴C analyses to ensure ¹⁴C measurements precise enough for a meaningful interpretation. The apparent ¹⁴C exposure durations of the two analyzed samples are 3.74±0.10 and 2.20±0.05 kyr, respectively; i.e., both are apparently younger than their ¹⁰Be counterparts. This trend and the low ${}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{10}}\mathrm{Be}$ concentration ratios of 1.67±0.05 and 1.62±0.04, respectively, are consistent with temporary burial of the surfaces and indicate that subglacial erosion was moderate and dominated by abrasion (see Sect. 3.3; Table 2).

The two apparent ¹⁰Be bedrock exposure durations from the Bockberg riegel are 0.05±0.02 kyr (STEI-16-7) and 0.62±0.03 kyr (STEI-16-5). We decided to not use them for ¹⁴C analyses due to the low cosmogenic nuclide inventory. These samples are therefore only briefly discussed in the following text.

We highlight the exceptionally young apparent ages of samples STEI-2 and STEI-16-7, prepared/measured at LDEO/CAMS and LN2C/ASTER, respectively (Table 1), showing that at these facilities we are able to detect ¹⁰Be signals that correspond to ages as young as ∼50 years with 20 %–40 % analytical uncertainty.

Modeling of the ¹⁰Be and ¹⁴C data yields cumulative exposure durations t_exp for STEI-16-12 and STEI-16-10 of 8.0±0.3 and 6.8±0.3 kyr (1σ analytical uncertainties), burial durations t_b of 2.1±0.1 and 3.2±0.1 kyr, and subglacial erosion depths E of 21±3 cm and 36±3 cm. This corresponds to erosion rates ε of 0.10 ± 0.02 mm yr⁻¹ and 0.11 ± 0.01 mm yr⁻¹, respectively. Following the assumption that samples on a transect parallel to the ice flow experienced the same exposure history on millennial timescales (Goehring et al., 2011), an average ice-free duration and standard deviation of 7.4±0.8 kyr can be deduced from the results of STEI-16-12 and STEI-16-10. Given that STEI-16-12 lies ∼160 m further outboard and 34 m higher in elevation than STEI-16-10, the longer exposure duration derived for STEI-16-12 could also either be due to shorter ice cover or be an artifact of a thinner ice cover above this sample, the latter potentially leading to small muogenic ¹⁴C production in the subglacial surface of the sample (see Sect. 3.3). Distinct exposure histories of the two sample locations are explored in Sect. 5.1.

Figure 6b shows the t_exp and E distributions that can explain the observed concentration of one nuclide alone, ¹⁰Be or ¹⁴C, resulting from the forward model (see Sect. 3.5) and represented by plain and dashed curves, respectively, for each sample (colored curves). The intersection of the two curves coincides with the average t_exp and E values, determined with the MCMC inversion, thus providing a visual check of the optimal model solution. Following this strategy, we added the same type of curves obtained from the ¹⁰Be concentrations alone of four other samples from Chüebergli riegel (grey curves; not shown are the results from the two innermost samples due to their low ¹⁰Be concentrations). This approach allows evaluation of the range of possible erosion depths and rates, illustrated by the yellow lines, even without combination with ¹⁴C analyses. Assuming that these four samples experienced a similar exposure–burial history to samples STEI-16-12 and STEI-16-10, their ¹⁰Be concentrations indicate subglacial erosion rates of >0.2 mm yr⁻¹ (Fig. 6b). This erosion rate estimate and the values inferred from the combined ¹⁰Be and ¹⁴C data (∼0.10–0.11 mm yr⁻¹) agree with those from Rhône Glacier (∼0.02–0.33 mm yr⁻¹) and the steep Trift riegel (0–>2 mm yr⁻¹) obtained with the same analytical method (Goehring et al., 2011; Steinemann et al., 2021).

5 Discussion

5.1 Holocene timing and duration of Steingletscher retreat

Given the similarity in ages of the Hublen bedrock (highest and stratigraphically outmost position; 11.5±0.3 ka) and the outer Hublen moraine (11.0±0.3 ka), we infer that these bedrock surfaces experienced a simple exposure history and that their age represents the timing of the first general deglaciation of the summit of Plateau Hublen at ∼11.5 ka, during the transition from the YD to the Holocene. This supports our assumption that the Hublen bedrock was free from significant ¹⁰Be inheritance from earlier interglacials (see Sect. 3.3), thus providing evidence that this is the case for any bedrock surface further inboard in Steingletscher's forefield.

This date of deglaciation of the plateau summit also implies that the glacier had a greater extent prior to ∼11.5 ka. However, chronological constraints on positions of Steingletscher prior to or during the YD do not currently exist.

Subsequent to the ∼11.5 ka deglaciation of the Hublen Plateau summit, Steingletscher experienced a slow oscillatory retreat for about 1.5 kyr, similar to other glaciers in the Alps (Protin et al., 2021), and deposited the three preserved early Holocene moraines on Plateaus Hublen and In Miseren. Steingletscher's forefield was thus most likely completely covered by ice until the deposition of the In Miseren moraine (10.0±0.9 ka). Chüebergli riegel was covered by at least 140 m of ice until then. Subsequently, the glacier retreated considerably, most probably for the first time in the Holocene. Although the amplitude of Steingletscher's recession at that time is unknown, we assume that Chüebergli and Bockberg riegels already became ice-free at ∼10 ka. This is supported by radiocarbon-dated subfossil wood relicts that melted out of several glaciers across the Swiss Alps between 1990 and 2006 CE, indicating that the tree line had been well above the wood sample localities and that these glaciers were equally small or smaller between ∼10 and ∼8.2 cal BP than between 1990 and 2006 CE (Hormes et al., 2001, 2006; Joerin et al., 2006, 2008; Nicolussi and Schluechter, 2012).

Figure 7 Holocene advance and retreat history of Steingletscher (l) in comparison with independent records. (a) Greenland summit temperatures (Kobashi et al., 2017). (b) Global temperature anomalies relative to 1961–1990 CE (Marcott et al., 2013). The black dashed line is the global annual mean temperature anomaly in 1991–2020 CE (https://data.giss.nasa.gov/gistemp/, last access: 10 January 2022). (c) Chironomid-based stacked summer temperature at 1000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (Heiri et al., 2015). The dashed lines are the modern July temperatures at the weather stations closest to Hinterburgsee and Stazersee, i.e. Meiringen (588 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) and Samedan (1708 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) (https://www.meteoswiss.admin.ch, last access: 10 January 2022), extrapolated to 1000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, using a lapse rate of 6 ${}^{\circ }\mathrm{C}/\mathrm{1}\mathrm{km}$ as in Heiri et al. (2015). (d) Pollen-inferred summer temperatures from the Rutor Glacier forefield (Italian Alps), modern July temperature and uncertainty as horizontal lines. Fluctuations of Rutor glacier in comparison to 1968 CE (Badino et al., 2018). (e) Speleothem-based elevation changes of Upper Grindelwald Glacier (central Swiss Alps) (Luetscher et al., 2011). (f) Recession periods of six Swiss glaciers based on radiocarbon-dated subfossil wood and peat (Joerin et al., 2006). (g, h) Glacier fluctuations of Great Aletsch (central Swiss Alps; length changes) and Mer de Glace (northern French Alps; altitude changes) based on radiocarbon-dating and dendrochronology (Holzhauser et al., 2005; Le Roy et al., 2015, respectively). The grey extension added to the curve of Great Aletsch is the loss in length until 2020 CE (http://www.glamos.ch/en/factsheet#/B36-26, last access: 10 January 2022). (i) Holocene advance and retreat scenario of Rhône Glacier, central Swiss Alps (Goehring et al., 2011). (j) Bracketing radiocarbon ages from Steingletscher forefield (King, 1974) with arrow lengths corresponding to the 2σ age intervals. Plain arrows stand for moraine deposits or glacier advances. Dashed arrows correspond to pollen-inferred climate cooling (blue) or warming (pink) trends. (k) Radiocarbon dates of organic samples melted out from Steingletscher and Steinlimigletscher, with line lengths corresponding to the 2σ age intervals (Hormes et al., 2006). (l) Steingletscher advance and retreat scenario (this study). Summed probability curves are mean ¹⁰Be moraine ages (Schimmelpfennig et al., 2014). The modeled in situ ¹⁴C and ¹⁰Be concentrations are trajectories at the bedrock surface (Eulerian perspective). Blue bands are periods dominated by glacier positions larger than in ∼2000 CE; pink bands are periods with dominantly smaller glacier extents.

Apart from two wood fragments collected in Steingletscher's forefield and radiocarbon-dated at 5.3–4.8 and 4.8–4.6 cal ka testifying to a glacier recession upstream of its ∼2000 CE extent, no further direct constraints exist for the timing and amplitudes of Steingletscher's fluctuations during the mid-Holocene (Hormes et al., 2006; Fig. 7k). Evidence of recession of the neighboring Steinlimigletscher at 5.9–5.3 cal ka suggests that Steingletscher was in a retreated position, too, at that time (Hormes et al., 2006). During the late Holocene, the glacier readvanced considerably at ∼3 ka ago and covered Chüebergli riegel again with ∼140 m of ice or more. This happened again from at the latest 0.6 ka onward during the period of the LIA until the general retreat trend from 1850 CE on, which eventually uncovered Chüebergli riegel completely around 2000 CE. Between ∼3 ka and the LIA, the exact timing of further Steingletscher fluctuations is unknown. The bracketing radiocarbon age of a wood fragment of ∼2.5–1.9 cal ka at the base of a ∼1.10 m deep peat profile on fluvioglacial deposits near the LIA glacier limit points to glacier retreat (King, 1974). Significant glacier recession to modern extents at that time has also been documented at other Alpine glaciers, in particular in the detailed late Holocene records at Great Aletsch (Holzhauser et al., 2005) and Mer de Glace (Le Roy et al., 2015) (Fig. 7g and h). At Steinlimigletscher, one organic silt sample radiocarbon-dated at ∼2.3–1.8 cal ka, attesting to glacier retreat beyond the ∼2000 CE extent, supports a local impact of this probably regional-scale warming event (Hormes et al., 2006; Fig. 7k). A subsequent readvance of Steingletscher is tentatively suggested by one ¹⁰Be boulder age of ∼1.9 ka (Schimmelpfennig et al., 2014). Another period of considerable retreat can be assumed for the Medieval Warm Period (∼1.3–0.7 ka), when several glaciers in the Alps are known to have retreated again to modern extents (Holzhauser et al., 2005; Le Roy et al., 2015; Fig. 7g and h).

Based on the ¹⁰Be and ¹⁴C data from Chüebergli riegel, we now explore the duration that the glacier had uncovered this area during the Holocene and put it into the context of the above-described glacier fluctuations. The apparent ¹⁰Be exposure durations from Chüebergli riegel alone already indicate that at least the outmost part of the profile has been exposed for a cumulative period of more than ∼5.9 kyr during the Holocene, i.e., that the glacier was smaller for that duration than its extent in 2000 CE. The general trend of decreasing apparent ¹⁰Be durations on the profile (by 2 orders of magnitude) indicates that the glacier temporarily covered and subglacially eroded the riegel again, otherwise similar apparent exposure durations would be expected for all samples on the profile. The erosion-corrected exposure durations modeled from the ¹⁰Be–¹⁴C data of the two samples indicate that the glacier was smaller than its 2000 CE extent for a total of ∼7.4 kyr during the Holocene. Based on the assumptions and knowledge of relatively rare and short periods of glacier retreat after 3 ka ago, it is most likely that almost all of this glacial retreat occurred between ∼10 and ∼3 ka. In Fig. 7l we propose two general scenarios of ice-free and ice-burial periods at Chüebergli riegel (i.e., glacier retreat and advance periods relative to 2000 CE) represented by theoretically reconstructed evolutions of ¹⁰Be and ¹⁴C concentrations at the surface of the locations where samples STEI-16-12 and -10 were collected. The simplest scenario, consistent with the measured concentrations of STEI-16-10, consists of a single ice-free period between 10 and 3 ka, and a continuously ice-covered period between ∼3 ka and modern times (orange curves). In the more complex scenario, which yields final theoretical concentrations similar to the measured concentrations of STEI-16-12, the glacier retreat beyond the 2000 CE glacier extent occurred not only between 10 and 3 ka, but also during the late Holocene warm periods at 2.5–2 and 1.3–0.7 ka (green curves). The sum of the exposure durations in this scenario amounts to 8.1 kyr and is thus indistinguishable from the exposure duration of 8.0±0.3 kyr analytically inferred from the ¹⁰Be and ¹⁴C measurements in sample STEI-16-12. For comparison, the scenario with a single exposure period between 10 and 2 ka is also shown for sample STEI-16-12 (also leading to consistent ¹⁰Be and ¹⁴C concentrations), which however is not consistent with the recorded glacier advance at ∼3 ka and therefore appears to be less realistic.

Taken together, it is most likely that frequent glacier oscillations led to temporary burial of Chüebergli riegel, and that during the ∼10–3 ka period, burial beneath ice occurred very rarely or not at all, while the ∼3 ka to 2000 CE period was dominated by Steingletscher advances and ice cover of Chüebergli riegel.

5.2 Comparison with Holocene glacier reconstructions in the European Alps

At a growing number of glacier sites, millennia-spanning records support the observed trend by providing evidence of significant and long-lasting glacier retreat during the early and mid-Holocene and progressive glacier re-advance during the late Holocene. Rutor Glacier in the western Italian Alps was smaller than its LIA maximum extent between 8.8 and 0.85 ka, and it was at least as contracted as in the 1960s between 8.8 and 3.7 ka (Badino et al., 2018; Porter and Orombelli, 1985; Fig. 7d). Upper Grindelwald Glacier experienced high-frequency elevation changes that were dominated by ice loss between 9.2 and 3.8 ka, followed by predominant advances close to the glacier's Holocene maximum (Luetscher et al., 2011; Fig. 7e). Radiocarbon-dated subfossil wood and peat fragments that had melted out of several retreating glaciers across Switzerland since 1990 CE revealed frequent and prolonged periods of recession between ∼10 and 3.5 ka followed by rare and short recessions during the late Holocene (Hormes et al., 2001; Joerin et al., 2006, 2008; Fig. 7f). Detailed reconstructions of late Holocene glacier fluctuations of Great Aletsch (central Swiss Alps) and Mer de Glace (northern French Alps), based on radiocarbon-dated founds, historical data and dendrochronology, reveal frequent large advances from ∼4–3 ka onward (Holzhauser et al., 2005; Le Roy et al., 2015) (Fig. 7g and h). Finally, Fig. 7i shows the Holocene retreat–advance scenario for Rhône Glacier (central Swiss Alps) inferred from combining the retreat duration of 6.4±0.5 kyr with chronological data on glacier advances from the Alps (Goehring et al., 2011, 2013). The general consistency of the findings from Steingletscher with the Rhône Glacier scenario and with the findings from the other locations in the Alps validate the approach applied here and confirms that the millennial retreat behavior of Steingletscher and Rhône Glacier represent regional glacial responses to Holocene climate warming.

Quantitative differences in the cumulative retreat durations inferred from the different methods seem to be consistent with the amplitude of glacier recession that is associated with the investigated sample material or location. While the in situ ¹⁴C–¹⁰Be dating approach suggests that glaciers were as retracted as the ∼2000 CE glacier extents for ∼7 kyr during the Holocene (Goehring et al., 2011; this study), the currently known recession periods inferred from the radiocarbon record by Joerin et al. (2006) add up to less, ∼5 ka. The reason could be that the subfossil organic material implies glacier retreat significantly upstream of the sampling locations and the ∼2000 CE glacier extents, as tree growth and peat development can only occur in deglaciated basins. This would roughly imply that since the final deglaciation (∼10 ka), glaciers in the Alps were similar in size to ∼2000 CE during a total of ∼2 kyr, while they were smaller than this extent for ∼5 kyr. However, we acknowledge that this interpretation is tentative and will need to be verified, as the observed differences in the cumulative retreat durations might also be inherent to uncertainties in the dating approaches. In particular, it is likely that some periods of retracted glaciers are still unknown because the associated radiocarbon-dateable material has not yet been discovered. Another source for differing results from the two methods could also derive from unaccounted-for in situ ¹⁴C production through thin ice (see Sect. 3.3). We also note that the existing data on Holocene glacier retreat do not allow verification of whether or not the glaciers completely vanished at some point during the Holocene.

5.3 Holocene glacier evolution in the Alps in the context of regional and global temperatures

Alpine summer temperature reconstructions based on chironomid and pollen assemblages are consistent with the Holocene glacier behavior in the Alps, showing a prolonged period of high temperatures during the early and mid-Holocene that might have been periodically and locally up to ∼1–3 ^∘C warmer than in 1981–2010 CE (e.g., Heiri et al., 2015; Badino et al., 2018; Fig. 7c and d). Also, the proxy-based global mean temperature reconstructions by Marcott et al. (2013) and Kaufman et al. (2020) reveal a very similar trend suggesting an early to mid-HTM followed by long-term cooling until the LIA (Fig. 7b). By contrast, model simulations of mean annual temperatures indicate that steady warming prevailed throughout the Holocene and that the recent decades are the warmest of the whole Holocene (Liu et al., 2014; Marsicek et al., 2018). Amongst the possible causes that have been proposed to explain these discrepancies, recent studies pay particular attention to the effect of seasonal biases (Marsicek et al., 2018; Affolter et al., 2019; Bova et al., 2021). According to this hypothesis, proxy-based global temperature reconstructions reflect warm-season rather than annual temperatures, because the growth of biogenic proxies is controlled by summer temperatures.

As mid-latitude glacier records at regional scale are mainly driven by summer temperature evolutions (Oerlemans, 2005; Solomina et al., 2015), our results corroborate the existence of an extended warm-season HTM. The fact that alpine glaciers are currently out of equilibrium with the accelerating anthropogenic warming, lagging behind by up to several decades, complicates a direct comparison of summer temperatures associated with glacier positions of the Holocene and the Anthropocene. Glaciers in similar settings and of similar size to Steingletscher have response times on the order of a few to ∼50 years (e.g., Oerlemans, 2012; Zekollari and Huybrechts, 2015), indicating that the summer temperatures responsible for Steingletscher's 2000 CE extent may have occurred in the middle or end of the 20th century, thus being 0.5–1 ^∘C less than in 2000 CE, according to the instrumental temperature record in the Alps (http://www.zamg.ac.at/histalp/, last access: 10 January 2022; Auer et al., 2007). Our data therefore imply that summer climate during the HTM was similarly warm to or warmer than during the second part of the 20th century. No further inferences can be drawn on the amplitude of warming.

The occurrence of an extended warm-season HTM seems to support the hypothesis of seasonality during the early and mid-Holocene. However, the amplitude of this seasonality cannot be determined from glacier chronologies. Further investigations are therefore necessary to resolve the controversial annual temperature evolution of the Holocene. A recent reconstruction of seasonally unbiased temperatures in Greenland, based on argon and nitrogen isotopes, provides evidence of several early and mid-Holocene episodes (amounting to 27 % of the Holocene) with temperatures above the average of the 1988–2015 CE period (Kobashi et al., 2017; Fig. 7a). This reconstruction is incompatible with the model-based steady annual warming and rather points to a hemispheric teleconnection with the trend of glacier fluctuations in the Alps.

In addition, we note that none of the continuous proxy records reveals a significant cooling event at ∼3 ka that could explain the deposition of late Holocene moraines outboard of the LIA extent of Steingletscher. Hints of a LIA-like cooling at 4–3 ka are only noticeable in the pollen-based summer temperature record at Rutor Glacier (Fig. 7d). This inconsistency will also need to be further investigated, because several other moraines of similar age are preserved across the Alps (Schimmelpfennig et al., 2012; Le Roy et al., 2017; Moran et al., 2017), while other Alpine records indicate glacier extents at ∼3 ka that are as short as in ∼ 2000 CE (Holzhauser et al., 2005; Le Roy et al., 2015; Fig. 7g and h).

Various drivers are relevant for Holocene climate change, i.e., external forcings at low (orbital summer insolation) and high frequency (volcanism and solar irradiance), feedback of the carbon cycle (greenhouse gases), and different climate boundary conditions linked to residual Northern Hemisphere ice sheets (Mayewski et al., 2004; Wanner et al., 2008). Our findings of Alpine glacier retreats and advances are in line with the current understanding of Holocene climate change. Orbital summer insolation modulates the long-term summer temperature evolution, thus driving millennial-scale glacier evolution in the northern middle and high latitudes (e.g., Solomina et al., 2015). Insolation is strongest in the early Holocene followed by progressive decrease, consistent with glacier retreat during the early and mid-Holocene and glacier re-expansion in the late Holocene. Volcanic eruptions and changes in solar irradiance superimpose centennial to decadal glacier fluctuations on the long-term trend during the late Holocene (e.g., Büntgen et al., 2016; Jomelli et al., 2016). Finally, while greenhouse gas concentrations were relatively stable over the Holocene, the accelerating anthropogenic greenhouse gas forcing has caused glaciers in the Alps and worldwide to retreat over the last century, with drastically increasing speed over the past few decades (Figs. 3b and 7g, h; e.g., Maurer et al., 2020; Roe et al., 2021; IPCC, 2007, 2013, 2021). The high sensitivity of Steingletscher to the moderate summer temperature amplitudes during the Holocene implies that the glacier will continue to melt and shrink dramatically and will most likely disappear if the human-induced warming is not reversed.

6 Conclusions

We find that Steingletscher responded highly sensitively to natural climate changes throughout the Holocene. It was as small as or smaller than its 2000 CE extent for a total of ∼7.4 kyr throughout the Holocene. No later than ∼10 ka, it shrank to its 2000 CE extent (or beyond) and advanced again to a LIA-like size at ∼3 ka, followed by expanded extents throughout much of the past 3000 years until the rapid general retreat that started in the 19th century. This Steingletscher record is consistent with the regional Holocene glacier evolution in the Alps suggesting that glaciers across the Alps were as small as or smaller than their extents around ∼2000 CE for most of the Holocene. The correlation between reconstructed summer temperature variability and the established glacier pattern demonstrates that the Alpine warm season temperatures between ∼10 and ∼3 ka, i.e., throughout the total of the HTM, were similar to or warmer than in recent decades. However, additional investigations are needed to fully understand whether or not the early and mid-Holocene was characterized by significant seasonality.

Uncertainty also remains with regard to the amplitude of glacier recessions and thus the magnitude of warming. The exact amplitude of glacier retreat can indeed not be inferred from the paired-nuclide approach applied here, because the dated bedrock does not delineate a past glacier extent. As a perspective, applying this approach at various distances from the current glacier front could be a valuable strategy to add further spatial constraints on the amplitudes of glacier recessions, provided subglacial bedrock erosion is low enough. In addition, the combination with complementary dating methods and glacier reconstruction approaches will help refine the long-term records in terms of both chronological and spatial constraints and thus add important knowledge to the Holocene glacier and climate picture in the region.