Alpine rock glacier activity over Holocene to modern timescales (western French Alps)

Active rock glaciers are some of the most frequent cryospheric landforms in midlatitude high-elevation mountain ranges. Their activity strongly influences the hydrology and geomorphology of alpine environments over short (years to decades) and long (centuries to millennia) timescales. Being conspicuous expressions of mountain permafrost and important water reserves in the form of ground ice, rock glaciers are seen as increasingly important actors in the geomorphological and hydrological evolution of mountain systems, especially in the context of current climate change. Over geological timescales, rock glaciers both reflect paleoclimate conditions and transport rock boulders produced by headwall erosion, and they therefore participate in shaping high mountain slopes. However, the dynamics of rock glaciers and their evolution over different timescales remain under-constrained.

In this study, we adopt a multi-method approach, including field observations, remote sensing, and geochronology, to investigate the rock glacier system of the Vallon de la Route (Combeynot Massif, western French Alps). Remotely sensed images and correlation techniques are used to document the displacement field of the rock glacier over timescales ranging from days to decades. Additionally, to estimate displacement over periods from centuries to millennia, we employ terrestrial cosmogenic nuclide (quartz ¹⁰Be) surface-exposure dating on rock boulder surfaces located along the central flow line of the rock glacier, targeting different longitudinal positions from the headwall to the rock glacier terminus.

The remote sensing analysis demonstrates that between 1960 and 2018 the two lower units of the rock glacier were motionless, the transitional unit presented an integrated surface velocity of 0.03±0.02 m a⁻¹, and the two upper active units above 2600 m a.s.l. showed a velocity between 0.14±0.08 and 0.15±0.05 m a⁻¹. Our results show ¹⁰Be surface-exposure ages ranging from 13.10±0.51 to 1.88±0.14 ka. The spatial distribution of dated rock glacier boulders reveals a first-order inverse correlation between ¹⁰Be surface-exposure age and elevation and a positive correlation with horizontal distance to the headwall. These observations support the hypothesis of rock boulders falling from the headwall and remaining on the glacier surface as they are transported down valley, and they may therefore be used to estimate rock glacier surface velocity over geological timescales. Our results also suggest that the rock glacier is characterized by two major phases of activity. The first phase, starting around 12 ka, displays a ¹⁰Be age gradient with a rock glacier surface velocity of about 0.45 m a⁻¹, following a quiescent period between ca. 6.2 and 3.4 ka before the emplacement of the present-day upper two active units. Climatic conditions have favored an integrated rock glacier motion of around 0.18 m a⁻¹ between 3.4 ka and present day. These results allow us to quantify back-wearing rates of the headwall of between 1.0 and 2.5 mm a⁻¹, higher than catchment-integrated denudation rates estimated over millennial timescales. This suggests that the rock glacier system promotes the maintenance of high rock wall erosion by acting as debris conveyor and allowing freshly exposed bedrock surfaces to be affected by erosion processes.

1 Introduction and motivations

Rock glaciers are important geomorphic structures influencing the evolution of high-elevation mountain environments. They are lobate or tongue-shaped assemblages of angular rock debris and ice. Rock glaciers move by the deformation of internal ice, conveying large caliber sediments from high-elevation steep slopes and cirque headwalls to their terminus at lower elevations (Barsch, 1977; Giardino and Vitek, 1988). In the context of current climate change, rock glaciers are considered one of the most resilient cryospheric bodies in alpine environments thanks to the insulating effect of their rocky carapace (Jones et al., 2019). Indeed, they represent an important water storage when mountain glaciers have disappeared (e.g., Jones et al., 2019). However, in situ measurements in the last few decades have shown rock glacier acceleration and destabilization associated with increasing air temperature in the European Alps (e.g., Marcer et al., 2021). Over geological timescales, rock glaciers participate actively in the development of asymmetrical mountain crests by conveying rock from lee side headwalls (where rockfall is the primary source of debris) to lower elevations in the valley (Gilbert, 1904; Johnson, 1980). Although rock glaciers have received considerable attention in the last couple decades, being catalogued in several geographic areas (see Jones et al., 2019, for the latest review), the formation and evolution of rock glaciers are still a subject of debate, with two main holistic views (see Haeberli et al., 2006; Berthling, 2011). On the one hand, rock glaciers are seen as periglacial features in which ice forms and ice content increases within debris interstices and thereby forms an ice–rock mixture that creeps by the influence of gravity and sufficient slope (e.g., Wahrhaftig and Cox, 1959; Ikeda et al., 2008). On the other hand, rock glaciers are also thought to be the remnants of glaciers with a deforming ice core protected by a continuous debris cover (e.g., Whalley, 1974; Monnier and Kinnard, 2015; Anderson et al., 2018).

The development of rock glaciers is a long process taking decades to thousands of years (e.g., Berthling, 2011). Their morphology, activity, and dynamics reflect present and past climates (e.g., temperature and precipitation fluctuations) and geomorphological forcing (rock and snow avalanching, bedrock structural patterns; Ikeda and Matsuoka, 2006; Kellerer-Pirklbauer and Rieckh, 2016; Jones et al., 2019). Rock glacier activity is categorized into active, transitional, and relict modes and has been recently updated based on geomorphological indicators (Delaloye and Echelard, 2020). An active rock glacier presents movement over most of its surface, whereas a transitional rock glacier presents low magnitude of movement only detectable by in situ or remote sensing measurement and/or is restricted to non-dominant areas. Finally, a relict rock glacier has no detectable movement and no morphological evidence of recent movement and/or ice content (Delaloye and Echelard, 2020). Rock glaciers have been documented in the European Alps that have accelerated with increasing temperature (Delaloye et al., 2010; Kellerer-Pirklbauer, 2017; Wirz et al., 2016; Eriksen et al., 2018; Kenner et al., 2018; Marcer et al., 2021), but when the ice content falls below a critical saturation threshold, rock glaciers stop creeping, turning from active into transitional and eventually relict mode (Sandeman and Ballantyne, 1996). Their activity is also controlled by the geomorphology of the surrounding topography. For instance, it has been suggested that when the rock boulder delivery rate and debris and ice incorporation become insufficient to sustain the insulation of the ice-rich part, the activity of the rock glacier will decrease and stop regardless of the rock glacier thermal state (Amschwand et al., 2021).

The relationships between external forcings and the activity of rock glaciers need to be better understood. Consequently, their past activity and in particular past vs. modern rock glacier surface velocity estimates must be quantified from annual to millennial timescales. This will allow us to assess how ongoing climate change has and will affect rock glaciers. Analytical advances over the past decades have allowed significant progress based on remote sensing tools for monitoring changes on high-mountain landforms (e.g., Necsoiu et al., 2016; Vivero and Lambiel, 2019; Blöthe et al., 2021; Robson et al., 2022). Indeed, methods such as lidar (Micheletti et al., 2017), InSAR (e.g., Liu et al., 2013; Barboux et al., 2014; Strozzi et al., 2020), aerial photogrammetry (e.g., Kaab et al., 1997), and unpiloted aerial vehicle systems (e.g., Dall'Asta et al., 2017; Vivero and Lambiel, 2019) have made a remarkable improvement to the temporal and spatial resolution of surface velocity surveys for rock glaciers. Recent studies have shown the feasibility of using high-resolution digital elevation models (DEMs) and ortho-rectified images produced from historical aerial and satellite images to reconstruct the surface velocity of rock glaciers over the last 7 decades (Fleischer et al., 2021; Vivero et al., 2021; Kääb et al., 2021; Cusicanqui et al., 2021). Extrapolations from short-term surface velocities have been used to estimate the rock glacier formation time and to reconstruct their activity over longer timescales (Kaab et al., 1997; Frauenfelder and Kááb, 2000; Bodin, 2013). However, it remains difficult to assess such extrapolations and to accurately constrain the long-term dynamics and morphological changes of rock glaciers without reliable estimates over centennial to millennial timescales.

To improve our understanding of rock glacier long-term dynamics and potential forcing mechanisms, relative and absolute dating methods have been applied on both active and relict rock glaciers (e.g., Haeberli, 2013; Amschwand et al., 2021). In rare cases, radiocarbon dating has been used on lacustrine sediments or trees buried by rock glaciers (Paasche et al., 2007) or on vegetal macrofossils found in old permafrost cores from a rock glacier (Krainer et al., 2015). Schmidt hammer methods have been employed to estimate the surface-exposure age of boulders from numerous rock glaciers (European Alps, Pyrenees, Norway, Iceland, and New Zealand) but such an approach requires local calibration surfaces and often only provides relative dating (e.g., Böhlert et al., 2011; Scapozza et al., 2014; Matthews and Wilson, 2015; Winkler and Lambiel, 2018). Similarly, lichenometry has been applied successfully on rock glaciers with stable rock boulders at the surface, although absolute dating requires calibration of this technique (e.g., Konrad et al., 1999). Optically stimulated luminescence has been used to quantify the travel time of buried fine sediments in rock glaciers (Swiss Alps), but large uncertainties potentially coming from pre-burial bleaching of fine sediments make this approach challenging to apply at a larger scale (Fuchs et al., 2013).

Terrestrial cosmogenic nuclide (TCN) dating has been successfully applied to constrain the exposure time of rock boulders at the surface of relict rock glaciers and their stabilization in the European Alps (Hippolyte et al., 2009; Steinemann et al., 2020), the Iberian Peninsula (Rodríguez-Rodríguez et al., 2017; Andrés et al., 2018; Palacios et al., 2020; García-Ruiz et al., 2000), Scotland (Sandeman and Ballantyne, 1996), and Iceland (Fernández-Fernández et al., 2020), demonstrating their potential as independent paleoclimate archives to reconstruct past permafrost development and to identify activity phases of rock glaciers (Andrés et al., 2018). Cossart et al. (2010) combined ¹⁰Be surface-exposure dating and weathering rind thickness to document three main generations of a rock glacier in the southern French Alps. Recently, two studies (in Iceland and Switzerland) have applied TCN dating on rock glacier systems composed of both active and relict units (Fernández-Fernández et al., 2020; Amschwand et al., 2021), showing deactivation and stabilization of the rock glacier at lower elevations and greater distances from the headwall.

The onset of rock glacier development in the high-elevation parts of the European Alps is thought to have started after the onset of glacier retreat following the Last Glacial Maximum (around 19–18 ka in the European Alps, e.g., Ivy-Ochs, 2015; Monegato et al., 2017; Lehmann et al., 2020). Chronologies of rock glacier development in the Alps of Austria, central Switzerland, and France have shown different rock glacier generations: during the Lateglacial (ca. 16 ka; Steinemann et al., 2020), during or shortly after the Younger Dryas (ca. 12 ka; Cossart et al., 2010; Steinemann et al., 2020; Charton et al., 2021), and during the late Holocene (Amschwand et al., 2021), probably at the end of Subboreal period ($\mathrm{5.2}/\mathrm{5.0}$–$\mathrm{4.3}/\mathrm{4.2}$ ka; Cossart et al., 2010) when the high-elevation cirques became ice-free.

The goal of this study is to reconstruct the activity and the surface velocities at different timescales of the rock glacier system of the Vallon de la Route (Combeynot Massif, western French Alps). A remote sensing approach utilizing image correlation over photogrammetric products allows us to reconstruct the surface displacement field of the rock glacier over the last 6 decades. Over longer periods (10³ to 10⁴ years), we apply TCN dating (quartz ¹⁰Be) to rock boulder surfaces at different positions along the central flow line of the rock glacier, from its terminus to its highest part, allowing the conversion of the ¹⁰Be surface-exposure ages into long-term surface displacement estimates. By discussing our estimates of rock glacier surface kinematics at different timescales, we show that it is possible to reconstruct the history of past activity of the rock glacier and to use rock glaciers as independent paleoclimate and paleo-geomorphological proxies to reveal the evolution of alpine environments.

Figure 1 (a) Regional map of the Combeynot Massif showing the Vallon de la Route rock glacier and surrounding topography (satellite image from © Google Earth 2020). Outline of the rock glacier (red), hydrography (cyan), and ridge line (white) from Bodin (2007). The inset in (a) shows location of the Vallon de la Route rock glacier within western Europe. (b) Map of the Vallon de la Route rock glacier, with the outline of the rock glacier, the main furrow structures, the outline of the bedrock outcrops (from Bodin, 2007), and the location and orientation of pictures presented in Fig. 2. The satellite image from Bing Aerial © Microsoft.

2 Study site

The Combeynot Massif (45^∘0^′ N–6^∘2^′ E) represents the northeastern part of the Ecrins Pelvoux Massif, located in the western French Alps (Fig. 1a). The Ecrins Pelvoux Massif has a high alpine topography, with its highest summit being the only one above 4000 m a.s.l. (Barre des Ecrins, 4101 m a.s.l.) and valley bottoms around 1000–1500 m a.s.l. Widespread U-shaped valley profiles, hanging valleys and glacial trimlines illustrate the imprint of Quaternary glaciations on the massif (Delunel et al., 2010; Valla et al., 2010; Le Roy et al., 2017). Paleo-glacier reconstructions since the Last Glacial Maximum have been previously constrained using mapping, interpolation of glacial features, and TCN dating applied on moraine deposits and glacially polished bedrock surfaces (Delunel, 2010). Evidence for Egesen stadial readvances (the alpine counterpart of the Younger Dryas) have been mapped and dated in several catchments of the massif (e.g., Coûteaux and Edouard, 1987; Charton et al., 2021; Chenet et al., 2016; Hofmann et al., 2019). Based on surface-exposure dating of moraine deposits, glacial advances in the Ecrins Pelvoux Massif have been identified from the Lateglacial to the early Holocene (around 11 ka, Hofmann et al., 2019). Glacial advances during the Neoglacial (from ca. 4.3 ka) have also been reconstructed (Le Roy et al., 2017). Modern glacierized terrain is characterized by small cirque and slope ice bodies, which covered 68.6 km² in 2009 (Gardent et al., 2014). The two largest valley glaciers remaining today are the Girose Glacier (5.1 km²) and the Glacier Blanc (4.8 km²); most of the other glaciers are cirques or debris-covered glaciers (Gardent et al., 2014; Fig. 1b in Le Roy et al., 2017).

Ranging from 1670 to 3155 m a.s.l., the Combeynot Massif hosts 33 active and transitional rock glacier systems and 38 relict rock glaciers (Bodin, 2013). These landforms range from 2000 to 2850 m a.s.l.; the mean elevation of the rooting zones of active rock glaciers is about 2700 m a.s.l., whereas the mean elevation of their frontal positions is 2620 m a.s.l. (Bodin, 2013). The Laurichard rock glacier, on the northern side of the Combeynot Massif (Fig. 1a), is the site of one of the longest geodetic surveys for surface velocity in the European Alps (since the late 1970s; Francou and Reynaud, 1992; Thibert et al., 2018).

The Combeynot Massif is comprised of a slice of granitic intrusion confined in volcanic–sedimentary gneiss; a flysch layer of the ultra-dauphinoise zone locally covers its eastern side (Barbier et al., 1973). The characteristic macro-crystalline scale fragility of the crystalline bedrock can be related to pre-Hercynian hydrothermal activity. A network of NNW–SSE faults and a high density of diaclases cut the Combeynot Massif, producing meter- and sub-meter-scale jointing of the bedrock (Francou and Reynaud, 1992). Thick superficial deposits (mostly coarse material) are the consequence of gravitational and nivo-periglacial processes (cryoclastic and avalanche activity; Francou, 1982). From a present-day climatic point of view, the Combeynot Massif is located in the transition zone between areas influenced by the Mediterranean climate and areas with a more Atlantic climate. Consequently, the local climatic setting is characterized by western frontal incursions and rainfall coming from the Italian side of the Alps and summer periods with low rainfall (Bodin, 2013). Cusicanqui et al. (2021) have performed weather reanalysis of data provided by the S2M (SAFRAN data from Durand et al., 2009a, b; Vernay et al., 2020). Their results show annual average temperature of 1.3±0.76 ^∘C and snowfall of 791±169 mm yr⁻¹ for the period 1958–2018 at the Laurichard site (elevation 2400 m a.s.l., northern exposure, 20^∘ slope; Cusicanqui et al., 2021). The average air temperature for the 1958–1990 period is 1.0 against 1.7 ^∘C for the period 1991–2018, revealing a change trend of +0.23 ^∘C per decade (Cusicanqui et al., 2021). Inversely, the mean number of days of snow cover during winter went from 221 for the 1958–1990 period to 200 for the 1990–2018 period (Cusicanqui et al., 2021).

Figure 2 (a) Photograph taken looking toward the NE and a geomorphological interpretation (b) of the rock glacier units III, IV, and V, with bedrock outcrops, talus, ridge, furrow and fronts. (c) Photograph taken looking toward the NEE of the relict unit I and sampled boulder (VR16). (d) Photograph taken looking toward the NW of the ridge I where two different boulders were sampled (VR5 and VR6). Picture locations and orientations are indicated in Fig. 1b.

The Vallon de la Route catchment is occupied by neither debris-free nor debris-covered glaciers, but is sheltering a rock glacier system of about 0.674 km². The total catchment area is about 4.41 km², is southwest facing, and ranges from 1960 to 3155 m a.s.l. It is bounded upstream by the highest peaks of the massif, the Tete de Pradieu (2879 m a.s.l.), the Roc Noir de Combeynot (3112 m a.s.l.), and the West Pic of Combeynot (3155 m a.s.l.) and downstream by the torrent of the Rif de la Planche (Fig. 1a and b). The rock glacier system is overhung by a debris source with area of about 5.351×10⁵ m² composed of leucogranitic bedrock (Fig. 1b). Geoelectrical measurements (performed in 2006) on the rock glacier have shown that the active layer reaches a maximum thickness of 9 m at 2630 m a.s.l. and that the ice-rich layer does not exceed 15 m and may consist of ice-cemented debris, with occasional layers of higher ice content (Bodin, 2013). The rock glacier system of the Vallon de la Route (Fig. 1) was chosen because it presents the ideal attributes of a rock glacier for this study, namely (i) it has an active snout and sharp edges, (ii) its total length is on the order of 1 km long (with both (i) and (ii) allowing remote sensing analysis), (iii) it is located in the middle of a valley (to be able to identify from where debris are coming from), and (iv) it has one single bedrock type in the headwall source (i.e., leucogranite allowing the use of TCN dating).

3 Methods

A multi-method approach combining geomorphological mapping and identification, remote sensing, and geochronology was used to reconstruct the history of activity of the rock glacier system of the Vallon de la Route (Combeynot Massif, western French Alps).

Figure 3 Surface displacement map (1960–2018) of the Vallon de la Route rock glacier. The red line is the longitudinal transect used to extract surface displacement (Fig. 5a). The solid black lines outline the different units of the rock glacier system. The dashed black line delimits the stable terrain control areas, and the red dots show the locations of the manual control points used to quantify uncertainties in the image correlation. Non-mapped areas are the result of the filtering process as described in Sect. 3.1.3.

3.1 Geomorphological mapping and identification and image correlation

3.1.1 Geomorphological mapping and identification

Geomorphological recognition was performed using the protocol described in the “Towards standard guidelines for inventorying rock glaciers, Baseline concepts” document of the International Permafrost Association Action group for rock glacier inventories and kinematics (Delanoye and Echelard, 2020). In situ visual inspection, such as observation of the steepness of the front, description of ridge and furrow topography, and size and shape distribution of the debris cover, was performed alongside a high-resolution lidar digital surface model (DSM, 0.5 m resolution) survey realized by SINTEGRA and landscape image analysis to geomorphologically classify the different landforms and their connection with each other. Five units (I, II, III, IV, and V) were identified from bottom to top according to in situ geomorphological observations such as their elevation, slope, vegetation cover, continuity, and the apparent activity of their landforms (elevation, average slope, and covered area were determined using the 0.5 m lidar high-resolution DSM; Figs. 1–4). Units are separated by either ridge and furrow topography or front, which are the expression of the gravity-driven buckle folding of rock glacier morphology (Frehner et al., 2015). In the present study, we focus our sampling strategy on 13 ridges annotated using the letters A to M (Table 1 and Figs. 2 and 6).

Figure 4 (a–f) Histograms of the surface displacements and median (±1σ) values extracted for each unit area and the entire rock glacier. (g) Histograms of the surface displacements and median (±1σ) values extracted from the entire control areas as shown in Figs. 3 and A1. Independent histograms of surface displacement values for each control area are presented in Fig. A3. (h) Distance between the two orthomosaics (1960 and 2018) manually estimated using a control point (stable features) as shown in Figs. 3 and A1 and Table 2.

Figure 5 (a) Surface displacement (1960–2018) extracted following the longitudinal transect (red line in Fig. 3) with identification of units I to V. The blue line represents the 10 m interpolation of the raw data. (b) Median surface displacement in each rock glacier unit with its ±1σ variability (see Fig. 4 for histograms). The cross-hatched pattern represents the detection limit defined by median value of the control areas (Fig. 4) used as a threshold value for detection of rock glacier surface displacement.

Figure 6 Map of the Vallon de la Route rock glacier, with units and sample locations. The hillshade DEM in the inset has been produced from a high-resolution lidar DSM (0.5 m resolution) survey realized by SINTEGRA (17 August 2012); the white lines show elevation isolines derived from the same data. The red line is the longitudinal transect used to extract surface displacement (Fig. 5a). Individual ¹⁰Be surface-exposure ages are shown with 1 standard deviation (Table 2). The lower-right inset shows the location of the study area within the Combeynot Massif (red box). The upper-left inset shows the samples (black dots) and ridges (red lines) distributed over the main furrow structures (black lines). The background is the orthomosaic of 2018 computed using a tri-stereo Pléaides 0.7 m acquisition.

Table 1 Geographic details within the rock glacier system of the Vallon de la Route (Combeynot Massif, France) of the samples collected for ¹⁰Be surface-exposure dating. Units represent areas of the rock glacier that have been defined geomorphologically. Ridges are also indicated (Fig. 3); for six of the ridges, two replicate boulders have been sampled. Ridge codes go from A at the lowest elevation to M at the highest elevation. The distance to the headwall has been measured following the center line of the rock glacier starting at the foot of the headwall (red line in Figs. 3 and 6).

3.1.2 Orthomosaic production

The reconstruction of the rock glacier surface displacement over decadal timescales is done using image correlation protocol between different orthomosaics. Here, we compared time series of two different orthomosaics (1960 and 2018). The oldest one is built from three historical aerial photographs acquired on 22 August 1960 and distributed by the French National Institute of Geographic Information and Forestry (IGN, https://remonterletemps.ign.fr/, last access: 15 April 2022) as black and white images of the original scanned films, with an average ground resolution of 0.59 m. The processing of those three images is based on the structure-from-motion and multi-view stereo photogrammetry approach implemented into the Agisoft Metashape software. The orthomosaic of 2018 was computed using tri-stereo Pléaides 0.7 m resolution acquisition (acquired on 12 August 2018) and an Ames Stereo Pipeline (ASP) (Shean et al., 2016) without ground control points (GCPs). ASP uses rational polynomial coefficients (RPCs) provided with the Pléaides images, eliminating the requirement of many high-accuracy GCPs. The 1960 orthomosaics were computed with Agisoft Metashape (version 1.6) software using 14 GCPs with coordinates collected from the IGN map service (https://www.geoportail.gouv.fr/, last access: 15 April 2022) and elevations using a high-resolution lidar DSM (0.5 m resolution) survey realized by SINTEGRA (17 August 2012). Details about the coordinates of the GCPs are given in Fig. A1 and Table A1. As the 2018 orthomosaic presented slight shifts due to processing without GCPs, the corresponding DSMs were then co-registered using a previous lidar high-resolution DSM and following the Nuth and Kääb (2011) methodology. Consequently, orthomosaics were then also shifted (translation only) with co-registration values (x and y displacements). All original orthomosaics were resampled at a 0.5 m resolution and set in a common 3465 pixel by 3541 pixel grid system.

3.1.3 Image correlation for surface displacement measurement

The 2D displacements of the rock glacier between orthomosaic pairs (1960 and 2018) were computed using the IMCORR module within the SAGA toolbox in QGIS (Scambos et al., 1992). The feature-tracking algorithm retrieves pixel patterns between two georeferenced images and attempts to match small subscenes (called “chips”) and produce shapefiles (points and lines) containing the 2D surface displacements. The program uses a normalized cross-covariance method based on a fast Fourier transform (Scambos et al., 1992). In the present study, several parameters of the algorithm were tested before settling on the following parameters: a search chip size of 128 pixels, a reference chip size of 64 pixels, and a grid space of 10 m.

The obtained surface displacements were first filtered with a threshold of 100 pixels for error on x- and y-direction estimates (IMCORR xerr and yerr values) removing about 1 % of the initial values. We then manually filtered according to (i) different local spatial coverage and artifacts (e.g., related to random local similarity of the coarse blocky surface; Bodin et al., 2018), (ii) lack of consistency of the displacement between neighboring vectors (difference > 30^∘), and (iii) outlier displacement values (Cusicanqui et al., 2021). Finally, 5.4 % of the points were removed. The ±1σ variability is calculated using all the pixels for each unit together with the median displacement. The quality of the results mentioned above is assessed through two analyses. First, the displacements obtained on the rock glacier system are compared to the measured displacements of six control areas where no displacement should be observed. Details of those control areas are given in Figs. 3, A1, and A2 and Table A3. These control areas were chosen to be outside of and around the rock glacier system and out of the scree field. The absence of movement (solifluction, creeping, landsliding) was determined by visually inspecting historical aerial photographs collected from the IGN-France and the two orthomosaics. Next, the potential mismatches between 11 control points in the two orthomosaics were manually measured. These manual control points were chosen based on topographic features of the surface of our studied area (blocks, cliff structures), and their stability was estimated using the historical aerial photographs collected from the IGN-France.

3.2 ¹⁰Be surface-exposure dating

TCN surface-exposure dating is based on the observation that when cosmic rays reach the Earth's surface, they produce cosmogenic isotopes in specific targets, such as the production of Beryllium-10 (¹⁰Be) in quartz minerals (e.g., Gosse and Phillips, 2001; Lifton et al., 2014). The in situ production of ¹⁰Be in quartz occurs predominantly within a few meters of Earth's surface and decreases exponentially with depth, such that by knowing a measured concentration of ¹⁰Be in the first centimeter of a rock surface and the local production rate of ¹⁰Be in such a rock, it is possible to calculate an apparent surface-exposure age (Portenga and Bierman, 2011). In this section, we detail how samples were collected, prepared, and analyzed, and we explain how external processes affecting the dataset can be estimated.

3.2.1 Sampling

The samples were collected with approximately 0.5 kg of rock material from the 13 ridges of interest (Tables 1 and 2; Fig. 6). For six of these ridges (A, D, G, I, K and L, Table 1 and Fig. 6), two different boulders were sampled to evaluate the reproducibility of our dating approach (Fig. 2c). The boulders were chosen following the central flow line, which was defined to be both at the center of the rock glacier width and perpendicular to the main ridge and furrow topography (red lines in Figs. 3 and 6) from the terminus to the high-elevation active lobes of the rock glacier system. Sampled boulders were chosen on the top of the ridges to minimize topographic shielding, snow-cover effect, and complex exposure histories (i.e., sediment or clast cover and late exhumation; Fig. 2c and d). Suitable boulders are large (>1.5 m) and in a stable position (Fig. 2c and d). Appropriate rock surfaces do not show signs of intense weathering or recent chipping. Sampling was done using a hammer, chisel, and small electric circular saw over 4 d (28 September–1 October 2020). The sampling details of each rock boulder, including coordinates, elevation, distance to the headwall, height of the sample from the ground, size of the boulder, and topographical shielding are summarized in Tables 1 and 2.

Table 2 Sample details, analytical data related to ¹⁰Be measurements and surface-exposure ages for the rock glacier system of the Vallon de la Route (Combeynot Massif, France), and inputs for the CREp (Cosmic Ray Exposure program) calculator (https://crep.otelo.univ-lorraine.fr/, last access: 3 March 2022; Martin et al., 2017). We used the production rate (4.16±0.10 at g⁻¹ a⁻¹) derived by Claude et al. (2014) at the Chironico landslide site. The ¹⁰Be surface-exposure ages are presented with ±1σ external error and ±1σ internal error (in brackets). Shielding correction includes the topographic shielding due to the surrounding landscape, and the dip of the sampled surface is calculated with the CRONUS-Earth online calculators (Balco et al., 2008, http://hess.ess.washington.edu/math, last access: 3 March 2022). The density of the rock boulder samples is assumed to be 2.75 g cm⁻³. The ¹⁰Be concentrations were corrected with a blank ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{9}}\mathrm{Be}$ ratio of $\mathrm{6.28}\pm \mathrm{0.534}\times {\mathrm{10}}^{-\mathrm{15}}$. Snow cover correction was calculated using the Gosse and Phillips (2001) equation with a snow density of 0.3 g cm⁻³, an attenuation length for fast neutrons in snow of 150 g cm⁻², and a cover of 50 cm of snow for 6 months of the year.

3.2.2 TCN preparation and ¹⁰Be measurement

Samples were crushed and sieved to retain the 200–500 µm grain size fraction. Beryllium extraction was performed at the GTC platform (ISTerre, France) using a chemical protocol adapted from Brown et al. (1991) and Merchel and Herpers (1999). Magnetic separation was used to isolate the quartz fraction, followed by successive leaching in an ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SiF}}_{\mathrm{6}}/\mathrm{HCl}$ mixture. In order to speed up the purification of quartz and save leaching cycles, magnetic separation with fine magnetite powder was performed between leaching cycles to remove partially altered minerals. Meteoric Be purification was achieved with three sequential dissolutions using diluted hydrofluoric acid (HF; Kohl and Nishiizumi, 1992). The purified quartz samples (13–26 g for each individual sample) were completely dissolved in concentrated HF after being spiked with ∼510 mg of a 998 mg L⁻¹ Be carrier solution (Scharlab ICP Standard, batch 16107901) to fix the ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{9}}\mathrm{Be}$ ratio (Table 2). After HF evaporation, perchloric and nitric acids were added and evaporated to remove organic compounds and fluorides. Beryllium was extracted by successive alkaline precipitations of Be(OH)₂ alternated with separation on anion and cation columns. Samples were then oxidized at 700 ^∘C for 1 h, and the final BeO mixture was mixed with Nb powder and loaded into nickel cathodes. The ¹⁰Be concentrations were measured at ASTER national facility (Cerege, France) against standard BeO_STD-11 ($\mathrm{1.191}\pm \mathrm{0.013}\times {\mathrm{10}}^{\mathrm{11}}$; Braucher et al., 2013) and were corrected for the full process blank with a ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{9}}\mathrm{Be}$ ratio of $\mathrm{6.278}\pm \mathrm{0.534}\times {\mathrm{10}}^{-\mathrm{15}}$.

3.2.3 Surface-exposure age calculation

Surface-exposure ages were computed with the CREp online calculator (Martin et al., 2017) using the Lifton–Sato–Dunai (LSD) scaling scheme (Lifton et al., 2014), the ERA40 atmospheric model (Uppala et al., 2005) and the Lifton VDM 2016 geomagnetic database (Pavón-Carrasco et al., 2014). We used the Chironico landslide production rate (with high-latitude sea level value of 4.16±0.10 atoms per gram per year; Claude et al., 2014), scaled according to the sample longitude, latitude, and elevation.

The production rate was corrected for sample thickness (Table 2) and density (2.75 g cm⁻³). Shielding correction includes the topographic shielding due to surrounding landscape and the dip of the sampled surface calculated with the CRONUS-Earth online calculator (Balco et al., 2008; http://hess.ess.washington.edu/math, last access: 3 March 2022). In addition, we explore the influence of snow cover attenuation using the Gosse and Phillips (2001) equation with a snow density of 0.3 g cm⁻³ and an attenuation length for fast neutrons in snow of 150 g cm⁻² (Delunel et al., 2014). According to a previous study that dated a rock avalanche less than 1 km north of our site (Chenet et al., 2016), we use an estimate of 50 cm cover of snow for 6 months of the year, the values that are most often cited in the literature for the Alpine regions for these altitudes (Hormes et al., 2008; Ivy-Ochs et al., 2006; Kelly et al., 2004; Schindelwig et al., 2012; Wirsig et al., 2016; Chenet et al., 2016).

3.2.4 Inheritance and pre-exposure estimation

The measured ¹⁰Be concentrations of rock glacier boulder surfaces should always be interpreted with caution as multiple external processes can affect them. Surface erosion can cause a depletion of ¹⁰Be concentration at the rock boulder surface, as can complex exposure histories (discontinuous exposure, snow and sediment cover), both of which would lead to an underestimation of the “accurate” ¹⁰Be surface-exposure age. Note that in this study we do not consider the effect of boulder surface erosion; consequently, our reported ¹⁰Be surface-exposure ages must be seen as minimum estimates. On the other hand, inheritance (i.e., headwall pre-exposure before rock collapse on the rock glacier) will lead to overestimation of the ¹⁰Be surface-exposure age. We choose two approaches to quantify the potential inheritance and pre-exposure bias. The first one is to use linear regression between ¹⁰Be surface-exposure age and distance to the headwall. Without inheritance, boulders located at the contact between the headwall and the talus slope (distance to the headwall equals to zero) should have a negligible ¹⁰Be concentration. Any ¹⁰Be concentration given by linear regression at the headwall would be interpreted as time spent (inheritance) in the cliff before the rockfall event. The second approach is to compare the ¹⁰Be concentration of samples collected on the same ridge. Any difference in ¹⁰Be concentration could be interpreted as a difference in time spent on the cliff face before the rockfall event. In this case, the inheritance time is calculated using the ¹⁰Be concentration difference between replicates using the ¹⁰Be production rate at the elevation of 2997 m a.s.l. (corresponding to the middle elevation of the cliff source).

3.3 Surface velocity estimation

Armed with both geochronological and remote sensing datasets, the time-averaged surface velocity of the rock glacier system can be estimated. When considering ¹⁰Be surface-exposure dating, the mean surface speeds are calculated by dividing the distance to the headwall of either the individual sample or the unit (mean distance to the headwall evaluated at all pixels inside within the unit limits) by their corresponding individual or median ¹⁰Be surface-exposure ages. Here we assumed that the ¹⁰Be surface-exposure age represents exposure between the time of the rockfall event that delivered the block to the surface of the headwall and talus transition and its arrival at the sample site. In contrast, the modern velocities derived from remote sensing analysis are calculated by dividing the measured displacement or the median values for each unit by the time between the two orthomosaics, in this case 58 years (1960–2018 period).

4 Results

According to the geomorphological mapping and identification approach described in Sect. 3.1.1, the rock glacier system was divided into five different units in which we apply the approach presented above to reconstruct the chronology of the rock glacier displacement since the onset of the Holocene.

4.1 Field observations

At the lower part of the rock glacier system, i.e., units I and II, with elevation ranges of 2585–2521 and 2718–2557 m a.s.l., respectively (and covering about 0.214 and 0.177 km², respectively), clearly display relict morphologies with no geomorphological evidence of recent surface movement. Their debris layers (on scales from meters to several meters) are highly covered by superficial material and vegetation (Fig. 2c). Unit III, spanning from 2685 to 2556 m a.s.l. and covering 0.155 km², was evaluated to be transitional. Indeed, the evidence of downslope movement is less visible (but still detectable) than for units IV and V, as the ridge-and-furrow topography is less prominent. Unit IV (ranging from 2735 and 2626 m a.s.l.) and unit V (ranging from 2867 and 2685 m a.s.l.) were classified as active and cover areas of about 0.542 and 0.729 km², respectively. They display evidence of downslope creep movement such as steep fronts, strongly marked ridge-and-furrow topography, absence of vegetation cover, and active layers composed of angular debris that are tens of centimeters to meters in size. These two units are talus connected, meaning that they are part of a downslope sequence including headwall, talus slope, and rock glacier (Fig. 2a and b). The delivery of debris is likely accomplished by rockfall activity, surface runoff, debris flow, and/or avalanche events from the headwall bedrock. The horizontal limit between the talus slope and rock glacier units can be estimated with about 50 m of uncertainty. This top to bottom organization from active to transitional and relict units is common in alpine settings (Amschwand et al., 2021). The surfaces of the boulders evolve along the rock glacier longitudinal transect. The boulders of units I and II are rounded and display quartz phenocrystals, a rugged surface, millimeter-scale weathered crust, and about 80 % lichen cover (Fig. 2c). On the other hand, boulders of units IV and V are more angular, with only about 10 % to 30 % of lichen cover and less obvious surface weathering features (Fig. 2d).

Table 3 The ¹⁰Be surface-exposure ages for each ridge of the rock glacier system of the Vallon de la Route (Combeynot Massif, France), where n represents the number of samples per ridge. For the ridges with replicates, the median values are reported with the standard variation ±1σ. For the inheritance estimate, the difference in ¹⁰Be concentration of each pair of replicates has been used, and ¹⁰Be surface-exposure ages have been recalculated assuming an origin in the headwall (at an elevation of 2997 m a.s.l.).

4.2 Image correlation

As described in Sect. 3.2, we used a pair of orthomosaics to reconstruct surface displacement of the rock glacier system. The 1960–2018 pair is the most suited to reconstruct the activity of the rock glacier system over decadal timescales. As the aim of this study is to compare modern to Holocene rock glacier activity, we focus on the integrated displacement over the longest and best-quality time series available, instead of focusing on shorter-scale variations in surface displacement over modern timescales (e.g., Cusicanqui et al., 2021). Figures 3–5 present the results obtained using the IMCORR feature-tracking module. The surface displacement of the control areas (dashed outlined area in Figs. 3 and A1) within stable terrain shows a median displacement of 0.79±0.43 m (Figs. 4g, A2 and Table A3). This value represents the accumulation of error from the orthomosaic production and the image correlation procedure. The quality of the orthomosaic production can be assessed using the statistics on the GCPs showing a median absolute error of 0.57±0.34 m (Table A1) and the manual control points presenting a median mismatch distance between the two orthomosaics of about 1.04±0.45 m (Fig. 4h and Table A2). This last value, being the highest of the three error estimations, is used hereafter as a threshold value to control the confidence level of our remote sensing analysis and should be considered the detection limit. Consequently, all rock glacier areas showing surface displacement lower than 1.04 m are consequently considered below the detection level (dashed area in Fig. 5).

Surface displacements calculated over the entire rock glacier system for the last 58 years (i.e., between 1960 and 2018) show a maximum value of 20.4 m, with a median displacement of 1.3 m and a standard deviation of 2.0 m over the entire rock glacier area (Figs. 3, 4f and Table A3). Note that those estimates integrate displacements over 58 years and do not allow us to assess whether the displacements have been steady or not. The spatial distribution (Fig. 3) and the longitudinal transect (Fig. 5a, location shown by the red line in Fig. 3) show that significant surface displacements are concentrated in the upper part of the rock glacier system. Indeed, units IV and V show median surface displacements over the 1960–2018 period of 8.5±2.9 and 8.3±4.9 m, respectively (values calculated for all the values inside each unit outline; Figs. 3, 4d–e and Tables 4, A3). This agrees with our classification as “active” from geomorphological observations (Fig. 2 and Sect. 4.1). Unit III presents surface displacements of 1.9±1.4 m (Figs. 3, 4c and Tables 4, A3), which is above the detection limit. Our classification of transitional activity (Sect. 4.1) thus seems appropriate. Finally, units I and II, with median displacement of 1.2±0.4 and 0.9±0.6 m, respectively, are similar within ±1σ and below the detection limit (Figs. 3, 4a–b and Tables 4, A3). These parts of the rock glacier system can thus be considered to be immobile over this period and correspond well with the relict classification determined from geomorphological observations (Sect. 4.1). The rock glacier becomes inactive around 945 m from the headwall, corresponding to an elevation of about 2600 m a.s.l. In addition, two displacement peaks can be observed in the upper part of the rock glacier system (Figs. 3 and 5a), one in unit IV and another in unit V, potentially indicating different debris sources for the two units.

Table 4 Median value of the ¹⁰Be surface-exposure ages for each unit of the rock glacier system of the Vallon de la Route (Combeynot Massif, France), where n represents the number of samples per unit. Surface velocity calculated from ¹⁰Be surface-exposure dating (distance from the headwall divided by ¹⁰Be surface-exposure age) and from image correlation (IMCORR, SAGA package on QGIS) of orthoimages (1960–2018 interval).

Figure 7 The ¹⁰Be surface-exposure ages of individual samples (a, b) and units (c, d, median values) plotted against elevation (a, c) and horizontal distance to the headwall (b, d). The dashed red and blue lines represent the linear regressions for cluster 1 (units I, II, and III) and 2 (units IV and V), respectively.

4.3 Surface-exposure dating results

Figures 6 and 7, together with Tables 2–4, present the analytical results of ¹⁰Be surface-exposure dating for each individual boulder sample, as well as for statistics within ridges and units. Our ¹⁰Be age results range from 13.10±0.51 to 1.88±0.14 ka for the entire dataset. The correction for snow cover shielding ranges between 7 % and 9 % between samples. In view of the controversy over whether wind could remove snow from moraine and rock glacier ridges during the Holocene period (Federici et al., 2008; Schimmelpfennig et al., 2014; Moran et al., 2016; Chenet et al., 2016), we refrain from correcting the output ¹⁰Be ages for snow cover in the following discussion. This implies that our reported ¹⁰Be ages, with neither snow cover nor surface erosion correction, should be considered as minimum estimates.

Our results clearly reveal a first-order inverse correlation between ¹⁰Be surface-exposure age and elevation and a positive correlation between ¹⁰Be surface-exposure age and horizontal distance from the headwall (Fig. 7). These correlations remain valid when we consider the median values of the ¹⁰Be surface-exposure ages for each unit (Fig. 7c and d, Table 3) and ridge (Fig. A3, Table 4).

Figure 8 Probability distributions of the individual ¹⁰Be surface-exposure age, sum, and median for each ridge (similar results are shown in Fig. A4 for each unit).

Visual inspection of the ¹⁰Be age dataset allows the identification of two clusters: cluster 1 includes the lowermost units (I, II, and II), whereas cluster 2 combines the highest units (IV and V). Cluster 1 presents ¹⁰Be ages between 9.25±0.40 and 13.1±0.51 ka, and cluster 2 shows ¹⁰Be ages between 1.88±0.14 and 4.88±0.29 ka. These results agree with the geomorphological classification we proposed for these rock glacier units, in which units I and II are viewed as relict, unit III is viewed as transitional, and units IV and V are viewed as active.

To assess the reproducibility of our dating approach, we sampled two different boulders on six of the rock glacier ridges (A, D, G, I, K and L; Figs. 2d, 6, and 8). The minimum and maximum horizontal distances between two replicates are about 8 and 82 m for ridges D (samples VR12 and VR13) and L (samples VR18 and VR19), respectively, a minimum elevation difference of about 2 m for ridges G (VR8 and VR9) and A (VR16 and VR17), and a maximum elevation difference of 14 m for ridge D (Fig. 6 and Tables 1 and 2). Significant variability in ¹⁰Be surface-exposure age occurs at the ridge scale, although it does not affect the correlations discussed above (Table 3 and Fig. 7). Ridge K presents the higher age variability (99 %, ¹⁰Be surface-exposure ages of 1.32±0.21 and 4.88±0.29 ka for samples VR2 and VR3, respectively; see Table 3 and Fig. 8d). The age variability for the other ridges is correlated with elevation and anticorrelated with distance to the headwall (26 %, 13 %, 10 %, 8 %, and 2 % of age variability for two samples per ridge for ridges L, I, G, D, and A, respectively; see Table 3 and Fig. 8). Only ridge A displays variability that is smaller than the absolute uncertainty on individual ¹⁰Be surface-exposure ages and may therefore be considered non-significant. The same pattern is observable for variability at the scale of the units (Table 4 and Fig. A4). Finally, samples from cluster 2 show much higher variability than samples from cluster 1.

Figure 9 (a) Rock glacier surface velocity from ¹⁰Be surface-exposure dating (distance from the headwall divided by the ¹⁰Be surface-exposure age, shown as red squares) and from orthoimage correlation (IMCORR, SAGA package in QGIS, 1960–2018 interval, shown as a blue line). (b) Median surface velocities are presented for each independent method and each individual unit alongside their standard deviation (±1σ). The dashed pattern represents the detection limit (0.026 m a⁻¹) defined by median value of the control areas in Fig. 3 (used as a threshold value to detect rock glacier surface movement).

4.4 Surface velocity

Figure 9 and Table 4 compile and illustrate the rock glacier surface velocities calculated from the ¹⁰Be surface-exposure dating and from the correlation of aerial and satellite orthomosaics. The surface velocities based upon ¹⁰Be surface-exposure dating range from 0.08±0.004 to 0.33±0.05 m a⁻¹ with a median value of 0.13 m a⁻¹ (Fig. 9). When we calculate the median value for the different units, the surface velocities range from 0.09±0.01 to 0.18±0.11 m a⁻¹ (Table 4).

For the remote sensing analysis, we define a detection limit of 0.02 m a⁻¹ corresponding to the median mismatch distance between the manual control point integrated over the 1960–2018 period (Figs. 4h, 5 and Table A2). As the displacements of units I and II show surface velocities identical to the detection limit, we consider them immobile over the last 6 decades. Measurable motion occurs above 2600 m a.s.l., with velocities of 0.03±0.03 m a⁻¹ in unit III. The upper units display higher velocities of the same order, about 0.15±0.05 and 0.14±0.08 m a⁻¹ for units IV and V, respectively (Fig. 9).

5 Discussion

The surface displacement reconstructions of the rock glacier system of the Vallon de la Route from both image correlation and ¹⁰Be surface-exposure dating provide interesting and original insights into the applicability of such a methodology on rock glacier landforms. It also suggests potential feedback between rock glacier activity, past climate, and geomorphological processes such as headwall erosion. Here we discuss the implications of the results obtained at the Vallon de la Route rock glacier system.

5.1 Inheritance or pre-exposure and loss or incomplete exposure

The evaluation of the source elevation of debris and the estimation of the inheritance and pre-exposure of investigated boulders was performed using linear regression of the ¹⁰Be dataset presented in Fig. 7, as explained in Sect. 3.4 (Amschwand et al., 2021). For instance, when we calculate the intercept of cluster 2 regression (i.e., elevation at which the ¹⁰Be surface-exposure age is null), we obtain an elevation of 2737 m a.s.l. Whereas if we include all the samples together (dotted red and black lines in Fig. 7a, respectively), we obtain 2748 m a.s.l. From our geomorphological observations, the elevation at which the talus slope connects to the headwall is close to 2880 m a.s.l. (mean elevation of the foot of upper headwall), which may safely be considered the elevation at which debris is delivered to the rock glacier. The difference between these elevations could lead to the interpretation that the ¹⁰Be surface-exposure ages are underestimated. We argue that this discrepancy can instead be explained by the fact that our sampling strategy targeted the biggest boulders at the surface of the rock glacier so that the likelihood of any burial event was minimized. The sampled meter- to meters-scale boulders might have rolled farther from the cliff and might therefore be incorporated onto the rock glacier surface at a higher distance and lower elevation than the present-day limit between the talus and the headwall. The relationship between ¹⁰Be surface-exposure age and elevation is also dependent on the relation between elevation and distance along the rock glacier (i.e., hypsometric distribution of the rock glacier surface) and the potential inheritance and pre-exposure effects on the measured ¹⁰Be concentrations.

Linear regression between horizontal distance from the headwall and ¹⁰Be surface-exposure age (Fig. 7b) allows us to quantify potential inheritance and pre-exposure bias (e.g., Amschwand et al., 2021b). By considering cluster 2, the samples have experienced inheritance and pre-exposure of about 2.16 ka (intercept of dashed red line in Fig. 7b). This assumes (i) that blocks fall at the talus–headwall contact and (ii) that the displacement rate is continuous over the temporal range of the considered ¹⁰Be surface-exposure ages. Once again, this inheritance and pre-exposure estimate must be put in proper geomorphic perspective, as units IV and V do not share the same headwall source (with potentially different slope aspects and thus different erosion rates). When regressed individually for units IV and V, the estimated inheritances are about 2.79 and 1.59 ka for units IV and V, respectively. For the above calculations, we only used samples from cluster 2, being the youngest and closest to the headwall. Samples from cluster 1, with the oldest ¹⁰Be surface-exposure ages and the greatest distances from the headwall, could also involve other biases that include non-continuous displacement rate over this timescale and loss or incomplete exposure due to surface erosion or tilting and burial of the sampled surface.

The variabilities of our ¹⁰Be surface-exposure ages of the ridges and the units are lower at low elevation and far from the headwall (Tables 3–4 and Figs. 8–A4). A first interpretation is that of the samples whose variability is smoothed between samples during transport and exposure on the rock glacier (cluster 1). We interpret this as highlighting that the events of tilting, burial, and/or erosion of the sampled boulders do not strongly influence the reported ¹⁰Be surface-exposure ages and that the variability likely arises instead from differing exposure times on the headwall prior to rockfall delivery onto the rock glacier surface. Secondly, the importance of inheritance and pre-exposure events would have less importance for the oldest ¹⁰Be surface-exposure ages than for the youngest ¹⁰Be surface-exposure ages. The high variability in cluster 2 of ¹⁰Be surface-exposure ages could be explained by variation in ¹⁰Be inheritance due to pre-exposure in the headwall. Using the difference in ¹⁰Be concentration between replicates, we estimated an inheritance. The calculated age differences are about 0.17±0.01, 0.67±0.03, 1.06±0.04, 0.39±0.03, 2.88±0.47, and 0.88±0.04 ka for ridges A, D, G, I, K, and L, respectively (Table 3). The median value of those results is 0.78±0.97 ka and can now be compared with the inheritance estimate of 2.16 ka derived from using linear regression between ¹⁰Be surface-exposure age and distance to the headwall of cluster 2.

We interpret the majority of the observed variability in ¹⁰Be surface-exposure ages as representing the stochastic nature of rockfall events. This leads to both different residence times of boulders in the headwall before rockfall and different sites of incorporation of boulders in the talus–rock glacier system. Interestingly, all sample ¹⁰Be surface-exposure ages suggest low inheritance compared to other settings in the European Alps. In the Mont Blanc Massif, for example, the more competent granitic spurs result in potential inheritance of >10 ka (Gallach et al., 2018), with a commensurately lower rate of debris supply and lower frequency of rockfall events (see Sect. 5.3 for a discussion about headwall erosion rates).

5.2 Surface velocity comparison and reconstruction

The 1960–2018 rock glacier surface velocities calculated from the correlation of aerial and satellite orthomosaics validate our proposed geomorphological classification for the activity of the different units (units I and II are relict, unit III is transitional, and units IV and V are active). This activity is occurring above 2600 m a.s.l. While no correlation between the distance from the headwall and the surface velocity is clearly visible, the variability in surface velocity is significantly higher for units IV and V. The difference between units IV and V likely reflects their different debris and snow avalanche sources; they may therefore have independent age and surface velocity profiles. This could therefore lead to an overestimation of the distance to the headwall for unit IV (as the central line is defined with respect to unit V, Fig. 6) and as a consequence to an overestimation of its surface velocity.

The velocity of the two upper units above 2600 m a.s.l., which show surface velocity of about ∼0.15 m a⁻¹ between 1960 and 2018, are slower than the reconstructed surface velocity of the Laurichard rock glacier, located in the adjacent cirque 1 km north of our study site and facing north (Cusicanqui et al., 2021). In this study, we quantified an acceleration of the average surface velocity changing from 0.5±0.09 m a⁻¹ for the 1952–1994 period to 1±0.09 m a⁻¹ for the period 2013–2017 for this landform ranging from 2430 to 2630 m a.s.l. The difference of activity between the two sites could be explained by the control of insolation and mean annual temperature on the permafrost conditions, which are more favorable to the north-facing slopes (Laurichard) than the southeast-facing slopes (Vallon de la Route). This has been highlighted by the Permafrost Favorability Index distribution of the area (Marcer et al., 2017). In addition, Marcer et al. (2021) have estimated the rock glacier kinematics over the past 7 decades for the entire French Alps using aerial orthoimagery. Mean displacement rates increased from 0.3 m a⁻¹ (for the period from 1948–1952 to 2001–2004) to 0.97 m a⁻¹ (for the period between 2001–2004 and 2008–2009) to 1.25 m a⁻¹ (from 2008–2009 to 2015–2017). Note that the values obtained in our study site are below the detection limit of this regional reconstruction (0.52 m a⁻¹ for the period from 2008–2009 to 2015–2017; Marcer et al., 2021).

Figure 10 Schematic interpretation of the evolution of the ¹⁰Be surface-exposure age patterns as a function of distance from the headwall, considering two phases of rock glacier activity and stochastic rockfall delivery of boulders. See Sect. 5.3 for details.

The velocity obtained by integrating ¹⁰Be surface-exposure age over the distance to the headwall, ranging from 0.08±0.004 to 0.33±0.05 m a⁻¹, are about the same order of magnitude as that obtained by Amschwand et al. (2021) using the same approach (∼0.3 m a⁻¹). Comparing the surface velocity obtained with our two datasets (orthomosaics correlation and ¹⁰Be surface-exposure dating), we see that active units (IV and V) share similar surface velocities between long-term and short-term approaches (blue and red datasets in Fig. 9). Integration of the short-term surface velocities over the late Holocene appears to predict the ¹⁰Be surface-exposure ages of investigated rock glacier boulders well. This suggests that the climatic and geomorphological conditions controlling the activity of the rock glacier have been stable above 2600 m a.s.l. over the last ca. 5 ka.

These observations should be put in a spatial perspective. The remote sensing analysis results in an estimate of the mean surface velocity over the entire area of the unit. In contrast, the velocity estimated from the ¹⁰Be surface-exposure dating was calculated from samples collected at the center of the rock glacier system, where the surface velocity is likely to be the fastest regarding a transversal cross section. The median velocity of the unit will be lower than the maximum centerline speed, as lower thicknesses at the margins should slow the surface speeds. It is therefore expected that the ¹⁰Be-based method would likely yield faster speeds than the remote sensing method. Consequently, the agreement and relationship between the two datasets and the two timescales should be interpreted with caution.

5.3 History of rock glacier activity

The Vallon de la Route cirque is occupied by a rock glacier system with well-defined rock glacier geomorphological attributes such as steep fronts, margins, ridges, and furrows topography (Figs. 1b and 2). No evidence of former occupation of the cirque by a clean ice or debris-covered glacier is visible (no moraine nor polished bedrock surface). Consequently, we interpret the correlations presented in Sect. 4.3 between the ¹⁰Be surface-exposure age and distance to the headwall to support the hypothesis that rock boulders originate from the headwall and are then transported downward on the surface of the rock glacier: the further from the headwall (and the lower the elevation) the boulder is, the older its ¹⁰Be surface-exposure age will be.

Following this reasoning and acknowledging the ¹⁰Be surface-exposure age distribution along the rock glacier, we propose a possible history for the rock glacier activity that includes two pulses of constant surface velocity. Figure 10 schematically presents our interpretation of the two clusters of ¹⁰Be surface-exposure ages according to their distances from the headwall (Figs. 7b and 11b). During the first phase of activity, boulders fall from the headwall onto the surface of the rock glacier (Fig. 10a). The random distance from the headwall at which the boulder is incorporated into the rock glacier represents the stochasticity of rockfall travel. Therefore, the ¹⁰Be inheritance, which corresponds to the residence time of the rock in the headwall, is also stochastic. As rockfall-derived boulders are transported down valley, both their ¹⁰Be surface-exposure ages and distances from headwall increase from these initial values (red lines in Fig. 10a). When the motion halts (presumably because the rock glacier thins beyond some threshold thickness), this first phase of activity ends and boulders remain stationary while their ¹⁰Be surface-exposure ages continue to increase (Fig. 10b). During this phase of inactivity, we consider neither snow nor rock avalanches to be active. Finally, a new phase of activity begins at the base of the talus (red points in Fig. 10c), and the new rock glacier overrides the up-valley boulders on the relict rock glacier (shadow points in Fig. 10c). In this conceptual model, we assume that the first phase of activity transported the boulders further downstream than the second phase of activity.

Figure 11 (a) Activity phases of the rock glacier, defined as the evolution of the rock glacier surface velocity in time (explored histories and the best-fitting history are shown as grey and green lines, respectively). (b) Relationship between ¹⁰Be surface-exposure age and distance to the headwall. Observed values of the clusters 1 and 2 are shown as red and blue squares, respectively. Modeled and best-fit values are shown in grey and green, respectively. (c) Likelihood distribution of the inversion exploring the time at which the first phase of activity ends (t2). See Sect. 5.3 for details.

We therefore argue that cluster 1, corresponding to ¹⁰Be surface-exposure ages of units I, II, and III, represents a first phase of activity of the rock glacier and that cluster 2, with units IV and V, represents a second phase of activity. To constrain both the timing and the surface velocities of these phases of activity, we numerically simulate the evolution of ¹⁰Be surface-exposure ages of boulders during their movement at the surface of the rock glacier (Fig. 11b). To represent its stochasticity, we prescribed the inheritance (Fig. 10b) with random values between 0 and 2.16 ka (as determined in Sect. 5.1). In the same way, the distance of incorporation of boulders on the rock glacier surface (Fig. 10a) is randomly sampled between 0 and a maximum of 100 m. In this model, three different times must be prescribed. The initiation of the first phase of activity is set at 12.1 ka (t1 in Fig. 11a), which is the ¹⁰Be surface-exposure median age of unit I. The second phase of activity is set to start at 3.4 ka (t3 in Fig. 11a) as this is the ¹⁰Be surface-exposure median age of unit III and is still active now. Only the time at which the first phase of activity ends cannot be directly extracted from the data (t2 in Fig. 11a). Consequently, we simulate the ¹⁰Be surface-exposure age structure of the rock glacier complex for 100 values of t2 ranging from t3 (3.4 ka) to the youngest age of cluster 1 (i.e., 9.25 ka for sample VR8). The velocity of phase 1 is calculated using the maximum distance a block traveled at the surface of the rock glacier (1720 m) and the activity duration of phase 1 (t1 and t2). The velocity of phase 2 is fixed at 0.22 m a⁻¹, a value calculated using the maximum distance a block traveled at the surface of the rock glacier during this phase (740 m) and the duration of activity (t3 in Fig. 11a).

The 100 simulations are evaluated against the measured ¹⁰Be concentrations using chi-square per degree of freedom, ${\mathit{\chi }}_{\mathit{\upsilon }}^{\mathrm{2}}=\frac{{\mathit{\chi }}^{\mathrm{2}}}{\mathit{\upsilon }}$. The chi-square is a weighted sum of squared deviations: ${\mathit{\chi }}^{\mathrm{2}}=\sum _i\frac{{\left(O_i-C_i\right)}^{\mathrm{2}}}{{\mathit{\sigma }}_i^{\mathrm{2}}}$, where σ is the variance on our ¹⁰Be dataset, O is the observations, and C is the modeled data. The degree of freedom, $\mathit{\upsilon }=n-m$, is equal to the number of observations (n) minus the number of fitted parameters (m; here there are four parameters: maximum inheritance, maximum distance of incorporation of a boulder on the rock glacier, the initiation of phase 1, i.e., t1, and the initiation of phase 2, i.e., t3). The likelihood probability function is then calculated as $L=\mathrm{1}/\exp \left(\frac{{\mathit{\chi }}_{\mathit{\upsilon }}^{\mathrm{2}}}{\mathrm{2}}\right)$ and normalized with its maximum to extract the median value and the standard variation (±1σ) of t2 (Fig. 11c). The inversion results suggest that the first phase of activity lasted from 12.1 to 6.26±1.96 (±1σ) ka, with a surface velocity of 0.29±0.15 (±1σ) m a⁻¹. The second phase of activity starts at 3.4 ka and has a surface velocity of 0.22 m a⁻¹ (Fig. 11a). The most recent phase of activity overrides the ¹⁰Be surface-exposure ages of the two upper units. We now discuss how these two phases of rock glacier activity can be connected to what is known about the paleo-environmental conditions in the western European Alps.

5.4 Reconstruction of paleo-environmental conditions

In the European Alps, the end of the Lateglacial (i.e., Egesen) led to readvancing of the mountain glaciers, reaching a maximum extent around 12 ka for both the eastern and western Alps (e.g., Ivy-Ochs et al., 2008; Protin et al., 2019; Hofmann et al., 2019). Directly downstream of the Vallon de la Route catchment, ¹⁰Be surface-exposure ages of moraines show ages of 13.0±1.1 and 12.4±1.5 ka, providing evidence for two stages of glacial advance or standstill at the end of the Lateglacial (Chenet et al., 2016). In a southern valley of the Ecrins Pelvoux Massif, morainic deposits at Pré de la Chaumette (down valley from Rougnoux Valley) have been dated at 12.5±0.6 ka (Hofmann et al., 2019). Immediately after the onset of the Alpine glacier retreat (12.2±1.5 ka in the same valley, Chenet et al., 2016), several advance episodes lasting ∼1 ka were identified in the Ecrins Massif. Dating in the southern part of this massif has shown glacial activity during the Lateglacial that may have lasted until the early Holocene before final glacial retreat (around 11 ka, Hofmann et al., 2019). Cossart et al. (2010) reported histories of glacier retreat and rock glacier generation in the Clarée Valley (about 10 km to the east of our study site). They identified three generations of rock glacier development during the second half of the Holocene, ranging in elevation from ∼2400 to 2800 m a.s.l. Recent dating of Charton et al., (2021) on two rock glaciers located ∼3 km to the north of our site and at an elevation of about 2050 m a.s.l. reveals ¹⁰Be surface-exposure ages of ca. 11 ka. They interpreted the ¹⁰Be surface-exposure ages as marking the end of activity of the rock glacier.

In the present study, we interpret the ¹⁰Be surface-exposure ages as being the sum of its residence time on the headwall cliff, the time spent traveling on the surface of the rock glacier, and the time since deactivation of the relict portion of the rock glacier for the relict units. We argue that rock boulders remain at the surface of the rock glacier while being transported down valley. This is supported by the small variability in ¹⁰Be surface-exposure ages obtained from the ridge replicates far from the headwall, which implies little occurrence of tilting and burial events (cf. Sect. 5.1). This is also supported by the rock boulder weathering evolution along the rock glacier, which displays more weathered surfaces far from the headwall (cf. Sect. 4.1). Consequently, in the Vallon de la Route catchment, the first phase of rock glacier activity appears to start around 12.1 ka from our oldest ¹⁰Be surface-exposure age (median value of the unit I). We suggest that this coincides with the final glacier retreat at the end of the Lateglacial at the onset of the warm period marking the Younger Dryas–Holocene transition (e.g., Liu et al., 2014). The upper mountain catchments and cirques then became free of glaciers, allowing the headwall and scree field to feed a rock glacier with debris thickness sufficient to insulate the ice. In this case, the rock glacier development would be geomorphically controlled by contrast to a climatic control (Cossart et al., 2010). As presented above, this first phase of activity would have ended around 6.26±1.96 (±1σ) ka. According to our reconstruction, the second generation of rock glacier development occurred starting at about 3.4 ka. This is earlier than the estimate proposed by Bodin (2013) from the relationship between slope and velocity based upon an estimated time of 1.7 ka for debris to reach the front of the unit V lying at around 2740 m a.s.l.

Whereas lateral glacier moraines dated using ¹⁰Be surface-exposure approach suggest several minor glacier re-advances between ca. 4.25 and 0.92 ka in the main glacierized valleys of the Ecrins-Pelvoux Massif (Le Roy et al., 2017), there is no evidence for glacial re-occupation during the Neoglacial (Little Ice Age) periods in the Vallon de la Route Massif. This specific cirque does not share an upstream connection with any of the main glacierized valleys of the massif. The headwall and scree taluses were ice free and could therefore feed the rock glacier system with debris and snow avalanches, maintaining the rock glacier during the last 3.4 ka, when the Neoglacial (Little Ice Age) climate was favorable for glacier and rock glacier activity.

Consequently, the inspection of the age structure of our rock glacier suggests two episodes of motion (Fig. 11a). The first phase, starting around 12.1 ka, displays a gradient in age with rock glacier surface velocity of about 0.45 m a⁻¹. The rock glacier activity then declines and stops at 6.26±1.96 ka. By around 3.4 ka, the climate again becomes conducive to rock glacier motion at elevations above 2600 m a.s.l. where the presently active upper two units have been emplaced. Even if the climate fluctuated during the late Holocene (e.g., Liu et al., 2014), the integrated velocities calculated with the ¹⁰Be surface-exposure ages reveal that the surface velocity of 0.18 m a⁻¹ (Fig. 11a) agrees with modern estimates.

5.5 Headwall erosion and implications

The reconstruction of the rock glacier activity provides a way to quantify the erosion rate of the surrounding headwalls over Holocene timescales (e.g., Humlum, 2000; Amschwand et al., 2021). Bodin (2013) has performed geophysical measurement of the rock glacier area and could determine a maximum thickness of the active layer of about 9 m and a maximum thickness of the ice-rich layer of about 15 m at 2630 m a.s.l. The entire area of the rock glacier is 6.745×10⁵ m². An approximation of the total volume of debris, considering a total thickness of between 9.5 m (active layer thickness of 5 m and ice-rich layer of 4.5 m thickness at 2720 m a.s.l.; Bodin, 2013) and 24 m (combining maximum of active layer thickness of 9 m and ice-rich layer of 15 m thickness at 2630 m a.s.l.; Bodin, 2013), gives 3.37×10⁶ and 13.49×10⁶ m³, respectively. Regarding the low surface velocity estimated between 1960 and 2018, between 0.14 and 0.03 m a⁻¹ over 42 % of the total area (units III, IV, and V), and the lack of movement for the other 58 % of the total rock glacier area (units I and II), we can assume a negligible ice concentration over the full volume of the rock glacier system. Taking this volume of debris over 12.1 ka (the median ¹⁰Be surface-exposure age of unit I) and considering that all boulders are derived from bedrock exposed above the rock glacier system (about 5.351×10⁵ m²), we can calculate a mean rate of erosion of the headwalls of between 1.0 and 2.5 mm a⁻¹. These results agree with estimates of erosion rate (∼1.2–4.1 mm a⁻¹) from the granodioritic headwall of Bleis Marscha rock glacier in the eastern part of the Swiss Alps (Amschwand et al., 2021). The catchment-wide denudation rate of the Ecrins-Pelvoux Massif has been estimated to range from around 0.3 to 1.1 mm a⁻¹ on millennial timescales using ¹⁰Be concentrations in stream sediment (Delunel et al., 2010), suggesting that frost-cracking processes strongly control the post-glacial topographic evolution of midlatitude mountain belts. The high erosion rates estimated in our study highlight that the steep rock walls that serve as the sources for debris on the rock glacier are retreating rapidly. This may be aided by the downstream conveyance of boulders and debris by the rock glaciers that prevent the headwalls from burying themselves in their own debris. This system therefore promotes the maintenance of high rockwall erosion rates, the development of cirques, and the possibility of distinctly asymmetric mountain ridges where the local climate is more conducive to rock glacier development on one side of a ridge than the other (Gilbert, 1904).

6 Conclusion

In this study, we quantitatively constrain the surface displacement field of an alpine rock glacier system over Holocene and modern timescales by using both remote sensing and geochronological datasets. The ¹⁰Be surface-exposure dating of individual boulders along the main center line of the rock glacier reveals ages from 13.1 to 1.8 ka, corresponding to elevations of 2535 and 2751 m a.s.l., respectively. Our first-order observation shows an inverse correlation between ¹⁰Be surface-exposure age and elevation and a positive correlation between ¹⁰Be surface-exposure age and distance from the headwall. This confirms the simple conceptual model in which rock debris falls from the headwall and remains at the surface as they are transported down valley by the rock glacier. Comparison of replicates from the transverse ridges along the rock glacier shows that loss and incomplete exposure due to surface erosion, burial, or tilting of the boulders is negligible. These replicates also show that ¹⁰Be concentrations of boulders close to the headwall can vary, which in turn provides constraint on the inheritance and pre-exposure of rock boulders. We estimate the possible maximum inheritance of 2.16 ka in our study area, corresponding to the residence time of boulders in the headwall.

Comparison of orthoimages from both aerial (1960) and satellite (2018) surveys shows that the rock glacier system is composed of the two uppermost active units with a surface velocity of about 0.14 m a⁻¹ at elevations from 2867 and 2626 m a.s.l. and a transitional unit with a surface velocity of about 0.03 m a⁻¹ at elevations between 2685 and 2556 m a.s.l. Analysis of a stable area outside of the rock glacier system constrains the detection limit to be 0.02 m a⁻¹. Reported values of less than this detection threshold imply that the downstream part of the rock glacier below 2600 m a.s.l. is presently immobile, confirming our geomorphic analysis of the feature as relict. The comparison of the surface velocities estimated using the ¹⁰Be surface-exposure dating relative to distance to the headwall and from the surface displacement integrated over the 1960–2018 period between the orthoimage surveys shows that late Holocene and modern velocities are comparable for the active units of the rock glacier system.

Comparison of these results for the entire rock glacier allows us to propose an activity history for the Vallon de la Route rock glacier that consists of two main phases of surface displacement. The first episode lasted between about 12.1 and 6.26±1.96 ka, with onset around the end of the Younger Dryas cooling event, when the cirques became ice free, allowing the headwall and scree field to feed the rock glacier with debris with insulation of the ice beneath. After a period of quiescence, the second phase of activity started around 3.4 ka and continues towards the present, possibly attributed to the more favorable climate of the Neoglacial (Little Ice Age) period. Finally, we use the surface velocities obtained using ¹⁰Be surface-exposure dating to reconstruct the erosion rate of the headwalls. The suggested erosion rates are between 1.0 and 2.5 mm a⁻¹. These are higher than catchment-wide denudation rates estimated over millennial timescales over the entire Ecrins-Pelvoux Massif, suggesting that the rock glacier system promotes the maintenance of high rock wall erosion (back-wearing) rates and the development of cirques. To go further in reconstructing the paleo-environmental conditions of this specific region, physically based numerical modeling of rock glacier evolution (e.g., Anderson et al., 2018) should be applied using the existing topography, the spatial patterns of ¹⁰Be surface-exposure ages, and the modern surface velocities as modeling targets.

Appendix A

Figure A1 Map showing the location of the ground control points, the control areas, and the manual control points used in the remote sensing analysis. Hillshade DEM has been produced from high-resolution lidar DSM (0.5 m resolution) survey realized by SINTEGRA (17 August 2012).

Table A1 Statistics of the ground control points (GCPs) used for the co-registrations of the two orthomosaics. Coordinates are mentioned in the EPSG:2154–RGF93 v1/Lambert-93 systems.

Figure A2 Histograms and median values of the surface displacement for each control area, as shown in Figs. 3 and A1.

Figure A3 The ¹⁰Be surface-exposure ages of individual ridges (a, b) plotted against elevation (a) and horizontal distance to the headwall (b). The dashed red and blue lines represent the linear regressions for cluster 1 (units I, II, and III) and 2 (units IV and V), respectively.

Figure A4 Probability distributions of the ¹⁰Be surface-exposure ages for individual samples with sum, mean, median, and standard deviation given for each unit.

Table A2 Results of the manual control point assessment between the orthomosaics of 1960 and 2018. Locations of the points are shown in Figs. 3 and A1. Coordinates are mentioned in the EPSG:2154 – RGF93 v1/Lambert-93 systems.

Table A3 Statistics of the rock glacier, units, and control areas from the IMCORR analysis. The total area of both orthomosaics is 3 067 391 m².