Introduction

Different types of rock slope instabilities, including both slow rock slope deformations and catastrophic failures, are commonly observed in formerly glaciated, alpine environments. Field evidence, monitoring and modelling studies suggest that the stability of large rock slopes in alpine environments is strongly controlled by their glacial history, including the extent and timing of glaciation, deglaciation and paraglacial landscape adjustment (Church and Ryder 1972; Ballantyne 2002; Ambrosi and Crosta 2006; Ballantyne et al. 2014). Moreover, several authors suggested that the catastrophic collapse of large slopes is often the outcome of long-term progressive failure processes, promoted and modulated by glacial erosion, deglaciation and associated stress perturbations, rock mass damage and hydrologic changes (Agliardi et al. 2001; Eberhardt et al. 2004; Bigot-cormier et al. 2005; Cossart et al. 2008; Grämiger et al. 2017; Riva et al. 2018; Fig. 1a).

Fig. 1
figure 1

Different processes occurring during deglaciation and promoting paraglacial slope instability. a Conceptual scheme of conditions, timescales and controlling factors of paraglacial rock slope evolution: stable vs. subcritically stressed conditions, sudden vs. progressive failure. “Early paraglacial” rock slope failures characterized by progressive evolution (b, Moosfluh, Swiss Alps, modified from GeoPraevent 2020)) or catastrophic failure (c, Eiger, Swiss Alps, modified from Oppikofer et al. 2008). Paraglacial rock slope failure characterized by a long-term stage of nucleation and slow continuous movements (d, Sechilienne, French Alps, modified from Cnrs-Omiv 2020), progressive rockslide differentiation (e, La Clapière, Franche Alps, modified from Durville et al. 2020) or catastrophic failure (f, Randa, Swiss Alps, modified from Gischig et al. 2011). See text for explanation

Glacial erosion leads to valley steepening, deepening and undercutting (Penck 1905; Herman et al. 2011; Sternai et al. 2013). These can have critical effects on rock slope stability by causing stress-redistribution within the slope (Augustinus 1995) and exposing persistent, potentially weak stratigraphic boundaries (Agliardi et al. 2019) or inherited persistent fractures, shear zones or pre-damaged slope sectors, which may promote slope failure by structural controls. During deglaciation, unloading due to the progressive removal of ice masses causes topographic stress tensor rotation (Harbour 1995). This may lead to tensile and/or shear fracture damage accumulation (Grämiger et al. 2017), increasing fracture intensity and, thus, enhanced rock mass permeability with downstream effects on the groundwater regime, also subject to changes of the meltwater supply (Riva et al. 2018). Rock damage accumulation during and after deglaciation results in progressive slope failure, mirrored by time-dependent displacements (slope creep; Emery 1978) and leading to the differentiation of rockslides increasingly sensitive to external factors. These include ice mass distribution and varying hydrogeological conditions, which greatly change across glacial, paraglacial and postglacial conditions, as well as rainfall and permafrost degradation (Church and Ryder 1972; Ballantyne 2002; Grämiger et al. 2017; Riva et al. 2018).

Slopes that experienced glacial erosion and deglaciation exhibit a full spectrum of failure modes and timing during the subsequent period of paraglacial adjustment (Bovis 1982; Ballantyne 2002). Based on direct observations or displacement monitoring datasets, some authors reported fast slope response during or soon after deglaciation. Kos et al. (2016) and Glueer et al. (2019) reported that the Moosfluh rock slope failure (Swiss Alps; Fig. 1b) started moving at low rates during the post-LIA (Little Ice Age) deglaciation of the Greater Aletsch glacier (Switzerland) until a critical acceleration stage in 2016 without massive collapse. Differently, deformation along the eastern Eiger flank (Swiss Alps, Fig. 1c) accelerated since 1935 (Messerli et al. 2018) and collapsed on July 2006, less than 150 years after the 1860 Lower Grindelwald glacier maximum (Oppikofer et al. 2008). Opposed to these “early paraglacial” rock slope instabilities, geomorphological analysis and absolute dating techniques (e.g. radiocarbon, OSL–optically-stimulated luminescence, Cosmogenic Radionuclide Dating; Pánek and Klimeš 2016) showed that slow rock slope deformations in alpine mountain ranges worldwide commonly activate thousands of years after the end of deglaciation and undergo very different evolution (Bovis 1990; Agliardi et al. 2009; Audemard et al. 2010; Crosta et al. 2013; Ballantyne et al. 2014; Riva et al. 2018). The Sechilienne landslide (Western Alps, France; Fig. 1d) nucleated at least 6 kyr after the valley was completely ice-free and kept on deforming without undergoing a progressive evolution (Le Roux et al. 2009). In the other hand, the La Clapière deep-seated gravitational slope deformation (Western Alps, France, Fig. 1e) started more than 3 kyr after the last deglaciation of the Tinée valley and evolved through successive stages of slope destabilization during the Holocene, until the differentiation of a large, active rockslide since the 1970s (Bigot-Cormier et al. 2005). Finally, the Randa rockslide in the Matter valley (Swiss Alps, Fig. 1f) failed suddenly in 1991, about 15 kyr after the post-LGM (Last Glacial Maximum) valley deglaciation, without significant precursor activity or geomorphic evidence of long-term deformation (Sartori et al. 2003; Willenberg et al. 2008).

Ongoing global warming caused the retreat of numerous glaciers (Clague 2009; Pachauri et al. 2014), compromising the stability of alpine slopes with major implications for landscape evolution and hazard assessments. Thus, a comprehensive modelling framework embracing the different styles and timing of paraglacial rock slope failure and their preconditioning factors is required to develop a predictive capability on the probable location, mode and timing of future large rock slope instabilities in presently deglaciating alpine areas. Such a modelling framework must be able to account for progressive slope failure processes, including progressive damage accumulation and time-dependent deformation; (Lacroix and Amitrano 2013; Riva et al. 2018) and their interactions with space- and time-variable ice loads and groundwater conditions, related to realistic, physically based scenarios of glaciation and deglaciation (e.g. volume, thickness, timing) that are usually poorly constrained in the field. In fact, while the maximum extent of glaciers is reconstructed by geomorphological evidence (Jäckli 1957; Florineth and Schlüchter 1998), the timing and rate of deglaciation since the LIA and later stages was inferred by modelling (Seguinot et al. 2018) and validated by absolute dating (Lehmann et al. 2018). However, regarding the post-LGM deglaciation, one can perform this operation only at few key locations (Wirsig et al. 2016). As a result, numerical modelling of the effects of the deglaciation on the slope stability is usually performed under the strong assumptions of constant ice melting rates (Agliardi et al. 2001; Grämiger et al. 2017; Riva et al. 2018) and without accounting for mass balance and ice dynamics.

Here, we use a novel numerical modelling approach that couples ice dynamics and mass balance (Braun et al. 1999; Sternai et al. 2013) to 2D time-dependent, damage-based finite-element progressive failure models (Amitrano and Helmstetter 2006; Riva et al. 2018). We systematically investigate the relationships among valley morphology, rock properties and deglaciation processes and their impact on the mode and timing of paraglacial large rock slope instability.

Methods

We simulated the long-term paraglacial evolution of simplified large rock slopes by modifying the time-dependent, damage-based 2D Finite-Element model DaDyn-RS (Riva et al. 2018), which couples continuum damage mechanics and time-to-failure laws. The model is forced by deglaciation time histories based on Stokes solutions of the ice flow, (Herman et al. 2011; Sternai et al. 2013). We performed a parametric study to evaluate the effects of initial ice thickness, deglaciation timing, slope strength and initial valley shape on the timing and mechanisms of rock slope failure nucleation and evolution under deglaciation forcing.

Deglaciation model

The ice dynamics during deglaciation is addressed using a numerical method, which was extensively described and applied in previous studies (Braun et al. 1999; Tomkin and Braun 2002; Herman and Braun 2008; Herman et al. 2011; Sternai et al. 2013). We provide here a brief overview of the equations governing the ice mass balance and flow, while more detail about the algorithm can be found in the above references.

Ice dynamics

The ice thickness, h, is calculated by solving the equation of ice mass conservation,

$$ \partial h/\partial t=\nabla \cdot Q+M $$

where t is time, Q is the ice flux and M is the ice net mass balance. M is calculated as a linear function of the surface temperature, Ts, which scales linearly with elevation:

$$ M=-\gamma Ts\kern0.5em \mathrm{with}\kern0.5em Ts={T}_0-\lambda \left(\mathrm{h}+\mathrm{z}\right) $$

where γ and λ are arbitrary constants, is the elevation and T0 is the temperature at sea level that varies through time to simulate glacial/interglacial cycles. Minimum (Mmin) and maximum (Mmax) rates of ice accumulation and ablation are also set to reproduce a realistic ice net mass balance. To describe ice sliding velocities, us, we apply a commonly used empirical law (Bindschadler et al. 1978; Van der Veen 1987; Paterson 1994)

$$ {u}_s={B}_s/{P}_e\kern0.5em \tau \kern0.5em {b}_3 $$

where the \( {\tau}_b \) is the ice basal shear stress, Bs is a sliding constant equal to 5x1013Pa3yr1m2 (Sternai et al. 2011) and pe is the effective pressure (i.e., pipw, where pi is the ice overburden pressure and pw is the water pressure). Assuming that glaciers are fed predominantly by snowfall, with a clean ablation area, we set Mmax, Mmin, γ and λ equal to 1 m/yr, −20 m/yr, 2 m/°C and 0.006 °C/m (Benn and Lehmkuhl 2000), respectively. Maximum T0 is set equal to 15°C and minimum T0 is set equal to 13, 11 and 9 °C to allow the ice thickness to cover the 20%, 50% and 80% of the valley slope for each assumed morphology.

Mechanical model: DaDyn-RS

DaDyn-RS is a damage-based, time-dependent 2D finite-element model conceived to simulate the progressive failure of large rock slopes subjected to gravitational loading and undergoing glaciation and deglaciation, accounting for the mechanisms and timing of damage accumulation, strain localization and hydro-mechanical differentiation of rockslides undergoing acceleration and potential catastrophic evolution (Riva et al. 2018). The model, implemented starting from the approach of Amitrano and Helmstetter (2006), combines a continuum damage mechanics law (Kachanov 1993) and a time-to-failure law (Atkinson 1984). The model is able to simulate both instantaneous and progressive failure of each finite element depending on the local stress state (i.e. critical vs. subcritical with respect to a Mohr-Coulomb damage threshold upscaled to fractured rock mass conditions). Damage is described in terms of the degradation of the effective rock deformation modulus at each damage event, simulating the effect of increasing crack density. In DaDyn-RS, damage-dependent permeability is generated based on dilatancy and connectivity criteria, and varies in space and time within the model (Riva et al. 2018). This allows simulating hydro-mechanical coupling within the model in a simplified form without a priori assumptions on slope hydrology and water table location. Thanks to this formulation, the model is able to simulate time-varying pore pressure distribution, strain localization within dilatant and compacting shear zones until the development of perched aquifers and the differentiation of rockslide masses accelerating until tertiary creep. Explicit model time is scaled to real time using absolute or relative constrains, e.g. absolute dating, climatic data and geomorphological evidence of deglaciation and slope failure timings (Lacroix and Amitrano 2013; Riva et al. 2018).

Deglaciation in DaDyn-RS

DaDyn-RS includes a specific routine simulating ice loading and removal to account for deglaciation in an explicitly time-dependent framework. The glacier applies transient hydrostatic loads to the FE mesh nodes belonging to the topography that vary as the maximum ice thickness is reduced. In previous works, the ice loads reduced following a linear time-dependent law (Lacroix and Amitrano 2013; Riva et al. 2018). Here, we implemented DaDyn-RS to account for realistic deglaciation histories based on the results from the deglaciation numerical model, respecting the ice mass balance and flow controls (Fig. 2a). The glacier thickness curves calculated at each model time step are imported into DaDyn-RS using a linear interpolation method (Fig. 2d). At each predicted time of failure, DaDyn-RS interrogates the interpolated deglaciation curve to derive the glacier thickness and the subsequent surface load to be applied. Furthermore, the deglaciation process is discretized in time intervals based on the initial glacier thickness and the deglaciation duration. This time interval, independent on the damage model time, and the predicted time of failure are compared and the shorter is used as next model time step, to avoid model instability.

Fig. 2
figure 2

Deglaciation and mechanical model setup. (a,b,c) Examples of valley deglaciation simulation respecting the ice mass balance and flow dynamics, extracted from the deglaciation model at different time steps. (d) Differences between constant-rate and simulated deglaciation curves; (e) DaDyn-RS model setup for the V-shaped valleys. U- and UV- shaped valleys geometry is represented respectively with yellow and white lines. Black dots on white lines represent the three monitoring points used in the analyses for each valley geometry. GSI values for each mesh elements and the initial glacier heights used in the analyses

Model setup

Model geometry was set up assuming a slope height of 1200 m, which can be considered as a representative valley scale relief of typical alpine valleys, with the valley bottom at 1500 m a.s.l. and the top of the slope at 2700 m a.s.l. (Fig. 2b). The cross-section of formerly glaciated alpine valleys can have a V-shape geometry (respectful of their fluvial origin) or a U-shaped geometry with or without suspended shoulders, depending on the amount of glacial erosional overprint resulting from the advance and retreat of glaciers during Quaternary glacial cycles (Penck 1905). Glacial reshaping of a valley is strongly conditioned by climate (e.g. frozen-based vs. temperate glaciers, Thomson et al. 2010), the ice flow velocity and local lithological and structural conditioning (Harbor 1995). Different authors show that the 2D morphology of glaciated valleys is well described by a parabolic topography that follows the eq. Y = aXb (Graf et al. 1970; Hirano and Aniya 1990; Harbor and Wheeler 1992; Harbor 1992, 1995; Svensson 2001), where b = 2 in U-shaped valleys (Hirano and Aniya 1988). The FR parameter described by Graf et al. (1970) can be used to define the parabola concavity and width.

Thus, three different half-space valley shapes were considered in the simulations, namely, V, U and UV-shaped valley. The V-shaped profile has a constant slope of 30° representing an average threshold hillslope value in alpine settings (Agliardi et al. 2013). The U and UV-shaped topography were thus derived using b and FR values for a typical slope carved in gneiss (Brook et al. 2004). The model domain was discretized using a three-noded triangular finite elements mesh of about 16,000 elements with maximum edge size of 60 m (i.e. 5% of the slope height). The optimal mesh size was selected to balance model accuracy and run time (Riva et al. 2018).

Mesh nodes located on the lateral boundaries were prevented to move in the horizontal direction, while those on the bottom boundary were fixed in both directions. Initial stress was initialized to the model by applying gravity. The damage increment parameter, D, and the time-to-failure law parameter, b, were selected from values calibrated by Riva et al. (2018). Simulations were stopped once the vertical displacement of the slope summit reaches 100 m (t100), in order to obtain full-developed slope deformation and failure mechanisms comparable among different simulations.

Simulations were performed by a computer equipped with Intel® Xeon® CPU 3.30 GHz, 64 GB RAM and running MATLAB R2019a 64bit. A single simulation required approximately between 48 and 56 h.

Parametric approach

We carried out parametric modelling considering different combinations of (a) material properties (strength and deformability); (b) initial glacier thickness; (c) deglaciation duration; and (d) valley geometry.

We considered three different sets of mechanical parameters, upscaled to fractured rock mass conditions using the Hoek-Brown approach (Hoek et al. 2002), typical of schist, gneiss and granitoid rocks (Table 1). These are the rock types more typically prone to large rock slope deformation in alpine contexts based on inventories statistics for the European Alps (Crosta et al. 2013; Agliardi et al. 2013).

Table 1 Numerical model parameters for the different slope materials

Values of the Geological Strength Index (GSI; Hoek et al. 2002) drawn from a normal distribution with a mean value of 80 and a standard deviation of 2 were selected to describe initial model conditions. In DaDyn-RS, a normal distribution of GSI values is used in initial model setup to account for material heterogeneity (Fig. 2e). The Poisson ratio was derived from sensitivity analyses performed with DaDyn RS by Riva (2017).

Three glacier thickness scenarios were simulated, namely, 20, 50 and 80% of the valley height to represent realistic alpine conditions. Deglaciation durations were chosen to vary between 1000 and 7000 years, considered reasonable duration for post-LGM deglaciation in the European Alps (Ivy-ochs et al. 2008; Wirsig et al. 2016). Dadyn-RS applies the time-to-failure law at every mobilized stress, thus the numerical model never reaches a true equilibrium state.

Combining different scenarios of valley shape (V, UV and U), glacier thickness (20, 50, 80% of slope relief) and deglaciation duration (1000, 3000, 7000 years) yield 27 different deglaciation histories. These were simulated for each considered rock type (schist, gneiss, granitoid), resulting in 81 DaDyn-RS simulation runs.

Results

Model results were extracted and analysed in terms of time and morphometry indicators (i.e. time t100, at which the 100-m vertical displacement threshold of the slope crest is reached, landslide area and depth, landslide toe elevation) that were then compared to trends recognized in natural instabilities in the European Alps, in order to evaluate the contribution by each factor to deformation timing and mechanisms. Rockslide differentiation was then evaluated by examining the presence of perched aquifer above the failure surface, displacements acceleration, damage end creep evolution.

Slope failure mechanisms

Different deformation and failure mechanisms arise from our simulation, depending on the considered valley shape. To compare the influence of valley morphology on instability evolution, water pressure distribution, shear strain and number of ruptures were extracted from each simulation results (Fig. 3). Landslide volume is a key parameter for risk management. Elevation differences between the valley bottom and the landslide toe directly influences the potential energy for slope movements and thus the probability of catastrophic events. Thus, for each simulation, we systematically quantified the unstable area (i.e. 2D proxy of volume) and the elevation difference between the elevations of shear zone daylighting and the valley bottom, and evaluated their dependence on initial glacier thickness and valley shape (Fig. 4).

Fig. 3
figure 3

Examples of the results for different valley shapes. Water pressure (a,d,g), shear strain (b.e.h) and the number of damage events (c,f,i) is showed for typical V-, UV- and U-valley shape simulations

Fig. 4
figure 4

Influence of valley shape (a,c,e) and glacier thickness (b,d,f) on landslide area, height of the valley toe on the slope bottom and landslide depth

In V-shaped valleys, slope instability shows the characters of deep-seated gravitational slope deformation (DSGSD), affecting the whole slope from the toe to the crest and beyond (Fig. 3). They are characterized by double-crested ridges, counterscarps and bulging in the lower sectors. In most simulations, horizontal and vertical displacement patterns outline compound mechanisms, with rotational kinematics in the upper part and more translational movements in the middle-lower sectors. Larger vertical displacements are observed in the upper sector of the slope and higher horizontal displacements in the lower portion. Topographic profiles show lowered upper-part and bulging at rockslide toe. The kinematics of slope instability is also significantly influenced by initial glacier thickness, since more rotational mechanisms are recognized with a thinner glacier while a more pronounced transition to shallower transitional instabilities can be noted with increasing ice thicknesses (Fig. 4f).

In UV-shaped valleys, rotational basal shear zones develop during the simulations (Fig. 3b). They daylight below the valley shoulders (i.e. slope breaks) and extend upslope behind the slope crest. Higher vertical displacements are located in the upper portion of the slope and greater horizontal displacements in the lower part of the slope. Simulations performed with glacier thickness equal to 80% of the valley depth result in shallower instability, daylighting at higher elevations; Fig. 4d and f).

In U-shaped valleys, slope instabilities are generally smaller with respect to those in V shaped valleys (Fig. 3c). In this case, all the movements can be classified as rotational. The depth of the failure surface depends from the slope material and the glacier depth.

We defined “unstable area”, the total area of all the elements with total displacement greater than 10 m (i.e. 10% of the maximum vertical displacement). We calculated the maximum and mean depth of the basal shear zone by extracting the number of ruptures (i.e. damage events) in each mesh element along different vertical sections along the slope profile and fitted them using stable statistical distribution. The mean value of the distribution represents the depth of the failure surface for each section. We then derived the mean value of the difference between the slope profile and the derived shear band, representing the mean value of the landslide depth. Finally, initial and final GSI distribution for each simulation were extracted for each scenario.

Results show that the size and depth of the simulated instabilities is negatively correlated with the glacier thickness while, for a given initial glacier thickness, the mean depth of the slope failure is generally smaller in UV- and U-shaped valleys (Fig. 4). Simulations with thinner glaciers are more similar to the results of simulations performed without deglaciation, which showed instabilities involving the entire slope and more rotational kinematics, while the rock mass within the unstable volume is characterized by little internal damage and deformation. GSI distributions show a lower frequency of damaged elements for almost all the simulations accounting for ticker glaciers, meaning more localized damage or smaller instabilities developing in these contexts. Simulations with thicker glaciers generally lead to a higher elevation landslide toe. Regardless of the ice thickness, smaller and deeper instabilities suspended on the valley bottom, develop in U- and UV-shaped valleys.

Slope failure timing

Model time was scaled to real time (i.e. years) using the deglaciation duration as a reference constrain. The time (t100) at which the 100-m vertical displacement threshold of the slope crest is reached was used to assess the temporal development of the instability. Lower values of t100 indicate a faster development of slope instability in response to deglaciation. Within an overall range of t100 between 750 and about 430 kyr from the beginning of deglaciation, we recognized five temporal cluster of instability development (Fig. 5) in particular: (a) Cluster 1, 350–450 kyr after the beginning of deglaciation: slope instabilities affecting V-shaped valleys carved in strong rocks (i.e. granitoids) and hosting a thin glacier (20–50%). These instabilities are generally deep and characterized by rotational kinematics, which implies a low influence of deglaciation processes on the instability development; (b) Cluster 2, around 200 kyr after the beginning of deglaciation: instabilities of V-shaped valleys in strong rocks, hosting a thick glacier (80%). These instabilities are shallower and smaller than in Cluster 1, and are characterized by more translational kinematics; (c) Cluster 3, 150–100 kyr after the beginning of deglaciation: V-shaped valleys carved in medium-strong rocks (i.e. gneiss) with a thin glacier lead to the development of almost rotational, deep instabilities in V-shaped valleys; (d) Cluster 4, 80–30 kyr after the beginning of deglaciation: instabilities developed in U and UV-shaped valleys in strong rocks (i.e. granitoids) and filled by thin glaciers (mostly 20–50%): and (e) Cluster 5, 50–0 kyr after the beginning of deglaciation: slope instabilities occurring in U- and UV-shaped valleys in granitoid rocks hosting thick glaciers, as well as instabilities affecting valleys with any shapes (V, UV, U) carved in gneiss and schist and hosting glaciers of different thickness.

Fig. 5
figure 5

Temporal clustering of simulated paraglacial slope instabilities. t100 extracted from the simulated scenarios. Horizontal coloured bars identify the recognized temporal instability clusters

We grouped the simulations based on rock type, initial glacier thickness and timing of deglaciation, to highlight possible relationships between these factors and the time-lag between the deglaciation and the instability development (Fig. 6). A general trend that can be observed in almost all the model results is the link between the slope material and the t100. Stronger rock masses results in longer slope response time. For any given slope material, the thicker the glacier the lower t100. Greater t100 are associated to V-shaped valleys, which is likely due to the different geometries of the instability (deeper and larger phenomena) with respect to U- and UV-shaped valleys. For most of the simulated scenarios, faster response can be recognized in correspondence with faster valley deglaciation (1000–3000 years).

Fig. 6
figure 6

Boxplot representing the influence of glacier thickness (a) and valley shape (b) on t100

Rockslide differentiation

Large rockslides are different in size, affected slope sectors, damage patterns and style of activity, and can evolve slowly over thousands of years or accelerate towards catastrophic failure (Eisbacher and Clague 1984; Crosta et al. 2013). The transition from slow deformation to fast movements through the hydro-mechanical differentiation of a rockslide, i.e. rockslide differentiation, occurs only in some cases and with very different timing with respect to deglaciation, allowing to distinguish “early paraglacial” from long-term paraglacial rock slope failures. The occurring of this transition depends from the acting mechanisms and damage evolution and it’s fundamental in a risk-oriented perspective. Differentiation occurs through the development of a continue basal cataclastic shear zone, with a dilatant character that evolves into a compacting one. This leads to permeability reduction and to the formation of perched aquifers within the rockslide mass. This increases the sensitivity to external hydrological forcing and the susceptibility to acceleration and catastrophic collapse (Agliardi and Crosta 2014; Crosta et al. 2014; Preisig et al. 2016; Alberti et al. 2019).

We define the rockslide differentiation as the formation of a continue localized failure surface. To understand the evolutionary trends and individuate the scenarios most prone to slow-to-fast transition, we concentrated on the indicators of the presence of this failure surface (Fig. 7), namely:

  1. a)

    Perched aquifer formation above the basal shear zone

  2. b)

    Displacement rate increase

  3. c)

    Damage evolution and creep behaviour

Fig. 7
figure 7

Rockslide differentiation outlined by the evolution of pore pressure (a, b, c), shear strain (d, e, f) and damage (g, h,i). Model states are extracted at specific times (1,2,3 dashed lines in the plot in inset j). Results represented in a, d, and g refer to time 1; images b, e and h were extracted at time 2; and final results c, f and i correspond to time 3. The cumulative damage rate and the vertical displacement of the slope crest and middle point through time are plotted in the lower graph

Perched aquifers

For each performed simulation, we quantified the occurrence and distribution of perched aquifers by calculating the percentage of unstable area (i.e. area above the basal shear zone) occupied by continuous clusters of wet finite elements (Fig. 8). As the failure surface propagates upwards, perched aquifer in the lower portion of the slope starts developing and extends as the shear zone reach maturity (i.e. lateral continuity and a compacting character). Thus, the location and distribution of perched aquifers within the landslide mass provide information about the degree of rockslide differentiation. Generally, simulations with larger initial glacier thickness result in more developed permeable clusters, due to the higher amount of distributed damage within the slope resulting from the removal of larger ice loads and to the evolution style of slopes initially covered by larger glaciers, characterized by shallower instability phenomena.

Fig. 8
figure 8

Percentage of landslide area occupied by perched aquifer for each simulation

Displacement rate increase

For each simulation, we extracted cumulative (horizontal and vertical) displacement-time curves at three reference points located at the top, middle and bottom of the slope profiles, respectively (Fig. 2). We identified knee-points (i.e. sharp increases of the displacement-time curve), analysing the variations in displacement-time curves inclination. These can be considered proxies of rockslide differentiation, assuming that sharp changes in displacement rates are related to strain localization and perched aquifer formation. The knee-point angle was also calculated comparing the mean slope of the displacement-time curve before and after the knee-point itself. It represents the sharpness of the behaviour change.

We recognize different knee-points clusters in the simulated instabilities developing between 0 and 30 kyr after the beginning of the deglaciation (Fig. 9), i.e. the timescale of the Alpine LGM (Ivy-Ochs 2015), namely: (a) knee-points between 12 and 17 kyr after the beginning of deglaciation, i.e. failures occurring today or in the recent past are associated with gneiss slopes in V-shaped valleys formerly occupied by a thick glacier, or in U-/UV-shaped valleys that hosted thin glaciers; (b) knee-points between 9 and 11 kyr after the beginning of deglaciation, i.e. Early Holocene rockslide differentiation, mainly observed in in UV-shaped gneiss slopes formerly buttressed by thin glaciers; and (c) knee-points between 4 and 6 kyr and between 3 and 0.5 kyr after the beginning of deglaciation, i.e. Lateglacial rockslide differentiation, associated with the instability of slopes with lower mechanical properties (schists and gneiss). In the latter case, the duration of deglaciation seems to be a key discriminant between rapid and slow instability development.

Fig. 9
figure 9

Knee-points extracted from the simulated scenarios, inset from 0 to 30 kyr after the beginning of the deglaciation. See 3.3.3 for “Jump” and Classical creep definition

Valleys affected by a thicker glacier and formed within more resistant materials were generally found to have a lower knee-point angle. This confirms a major external influence on the instability evolution when glacier thickness is higher. Strong slope materials lead to a smoother behaviour change, while for gneiss and most of all for schist, the displacement curves changes were found to be sharper.

Damage evolution and creep behaviour

Further insights about the instability evolution derive from the temporal evolution of damage. Different damage patterns could be identified depending on the simulated valley geometry. (a) In V-shaped valleys, damage initially affects the lower slope sector and then propagates upward, resulting in strain localization along a single shear zone and increased displacement rates. In the simulations showing a more transitional behaviour (glacier thickness = 80%), multiple shear zones start developing progressively at the bottom of the slope and migrate upwards with a retrogressive behaviour during the simulation. (b) In UV-shaped valleys, damage is initially distributed between the bottom and the shoulder of the valley, and then it localises and propagates upward resulting in a rotational kinematics. In some cases, especially those involving a thicker glacier and longer deglaciation time, an initial relatively small rockslide forms around the shoulder of the valley and subsequently a larger deep-seated phenomenon develops. (c) In U-shaped valleys, damage generally initiates in the middle part of the slope. When shallower instabilities develop, smaller failure surfaces form in the middle part of the slope and then propagate upward, with a retrogressive behaviour.

The temporal evolution of damage accumulation is similar in all the simulations, with a first stage of distributed damage associated with slow displacements, followed by an acceleration stage when damaged areas connect and damage localizes. The first stage corresponds to primary creep, while the second one involves the onset of transition to tertiary creep. To systematically discriminate between these two phases, the cumulative damage curve and the rate of the cumulative damage versus time were extracted from each simulation, outlining different typical creep styles. We recognized a classical creep behaviour with primary creep evolving into tertiary creep. Here the changing point is estimated by pinpointing the time at which the minimum rate of cumulative number of ruptures in time is reached (grey symbols in Figs. 9 and 10b). This most frequently characterizes simulations of U-or UV-shaped valleys, with fast deglaciation and a thick glacier on schists and gneiss. A more complex stepped behaviour (white symbols in Figs. 9 and 10a) is typical of V- and UV-shaped valleys with a thick glacier, in agreement with the sequential development of multiple slope failures at different scales. The timing of the curve steps correlates with sharp increases in displacement rates in the middle- and low-slope sectors. Finally, some simulations did not reach tertiary creep behaviour at the end of the simulations. This means that more time is required to fully differentiate the instability. In these cases, the slope instability is related to a large amount of diffuse damages rather than to a continue failure surface.

Fig. 10
figure 10

Conceptual illustration of the failure behaviour individuated observing the cumulative damage rate curves: (a) developing through multiple failure surfaces and (b) along one single failure surface. Crest and medium slope vertical displacements and cumulative damage rate for different simulations showing the “jump” creep behaviour (a) and the classical creep behaviour (b)

Combining all the gathered information, glacier thickness and slope material seem to have a strong influence on the transition to faster differentiated movements. Usually the transition from primary to tertiary creep anticipates a sharp increase of the registered amount of displacements. This is valid both when looking at instabilities affecting the whole slope, e.g. Fig. 10b, and when smaller retrogressive failures (Fig. 10a) or secondary phenomena are involved. Gneiss and schist rock slopes seem to be more prone to transition from slow-moving to catastrophic phenomena.

Discussion

Understanding the preparatory and triggering factors of paraglacial large rock slope instabilities is key to predict their future onset and evolution in alpine areas characterized by an increasing anthropogenic pressure. These factors have been previously recognized on both the local (Eberhardt et al. 2004; Bigot-cormier et al. 2005; Agliardi et al. 2009; Schwartz et al. 2017; Riva et al. 2018) and regional scale (Jarman 2006; Cossart et al. 2008; Prager et al. 2008; Agliardi et al. 2013; Crosta et al. 2013; Pedrazzini et al. 2016). Numerical modelling has been widely used to simulate the influence of glacier downwasting on rock slope stability (Agliardi et al. 2001; Eberhardt et al. 2004; Fischer et al. 2010; Grämiger et al. 2017). However, the site-specific character of these studies and the oversimplification of deglaciation modelling (Eberhardt et al. 2004; Baroni et al. 2014; Bacenetti et al. 2015; Grämiger et al. 2017) hamper a more general understanding of the interplay and relative importance of different factors related to the glaciation and deglaciation history on paraglacial slope stability. In this perspective, we performed a systematic, parametric evaluation of paraglacial rock slope instability, including mechanisms and timing, by coupling realistic deglaciation histories derived by a physically based ice dynamic model (Herman et al. 2011; Sternai et al. 2013) and an explicitly time-dependent model of progressive rock slope failure.

Our parametric study allows evaluating the influence of rock mass strength, initial glacier thickness, valley shape and deglaciation rate on the mechanisms and temporal evolution of paraglacial rock slope failures. In order to reduce the number of considered variables and achieve a consistent understanding of paraglacial controls on slope stability, we omitted the influence of site-specific stratigraphic, inherited structural or active tectonic controls (e.g. seismicity) that can significantly affect the deformation and failure mechanisms and timing on the local scale (Agliardi et al. 2019). For the same reason, ridge geometry was simplified, even if Kinakin and Stead (2005) highlighted using numerical modelling that different processes are responsible for the formation of slope deformations in different type of ridge. Compared to findings relative to a single-case study, this allowed to draw more general and robust conclusions about the factors influencing the instability development.

The modelled timing of landslide differentiation shows a good agreement with real case studies developing in similar scenarios. For instance, comparable timing with respect to cluster 1 in Fig. 9 (i.e. between 12 and 17 kyr after the beginning of the deglaciation) was derived for the Lavini di Marco rock avalanche (NE Italian Alps) and the 1991 Randa rockslide (Swiss Alps, Sartori et al. 2003). These failures occur in valleys covered by thick glaciers during the LGM. Our results corroborate conclusions by Eberhardt et al. (2004) regarding the latter case study, suggesting how the failure process could have been initiated by the deglaciation. Rock avalanches in Val Viola (Italian Alps), dated between 4000 and 8600 years after the complete valley deglaciation (Hormes et al. 2008) shows similarities with cluster 2 in Fig. 9 (i.e. between 9 and 11 kyr after the beginning of the deglaciation). Glacial over-deepening of the valley was not excluded as preparatory factor for instability development. The Sechilienne landslide, located in a micashist U-shaped valley in the wester Alps entirely filled by the Romanche Glacier at the LGM, can be included in cluster 3 (i.e. between 4 and 6 and between 3 and 0.5 kyr after the beginning of the deglaciation). In this case, the instability developed about 6–8 kyr BP, thus a minimum of 5 kyr after glacier retreat (Schwartz et al. 2017). Landslide kinematics, involving the entire slope, with the toe suspended on the valley bottom and a depth of the basal surface of about 150 m (Le Roux et al. 2011), agree with the results of our simulations. The initial destabilization of La Clapiere landslide (Argentera massif, SE France), developing in gneissic material, was dated with a delay of more than 3 kyr after the end of the deglaciation of the Tinee valley at approximately 18.5 kyr BP (Bigot-cormier et al. 2005). The Chironico landslide (southern Swiss Alps), mainly developing in gneiss material, was dated about 4 kyr after valley deglaciation. During the Last Glacial Maximum (LGM), the Valle Leventina was almost entirely filled with ice except for nunataks. Valley steepening, stress redistribution during glacial unloading and subsequent damage propagation were indicated among the predisposing factor for the Chironico failure (Claude et al. 2014), consistently with our results. Similar delays (about 6 kyr after deglaciation) were recognized also in the Canadian Cascade Mountains (Orwin and Clague 2004). The Cheam rock avalanche (Fraser Valley, Canada), involving about 175 million m3 of relatively weak materials, took place on the steep flank of a glacier-shaped valley, covered by high glacier thickness during the Fraser Glaciation (30,000–10,000 kyr BP) and underwent to long deglaciation. The Cheam landslide is the largest known catastrophic landslide in western Canada. The simulation of similar scenarios in our analyses was identified as prone to instability differentiation.

In general, rock mass strength is positively correlated with the timing of slope response to deglaciation, with weaker slopes failing earlier other conditions being equal. This is consistent with previous works relating rock strength with the extent and geometry of glacial erosion (Augustinus 1995), the amount of in situ stress relaxation (Engelder and Sbar 1977) and the regional distribution of slope instability development (Holm et al. 2004). Our results further highlight how discontinuity propagation due to debuttressing and stress-release is more effective in already damaged materials. Slope undercutting leads to more unstable scenarios when weaker or jointed rock masses are involved.

Melting thicker glaciers results in relatively shallow rock slope instabilities. They daylight at higher elevation and are characterize by a higher susceptibility to differentiate into rockslides with a shorter time lags between deglaciation and landslide enucleation. Generally, the thicker the glacier, the greater the influence of its downwastage on slope stability. Debutressing also influences rock mass permeability and groundwater distribution, which depends on rock mass damage in the different slope sectors. High glacier thickness was found to be correlated to landslide differentiation, being one of the main factors controlling slope evolution together with slope lithology. Shear stress and damage propagation often show a retrogressive behaviour, with the initial formation of smaller failure surfaces lower portion of the slope and subsequent upward migration. Furthermore, the removal of thick glaciers resulted in the development of secondary instabilities and, in general, high risk scenarios involving e.g. landslide differentiation and suspended rockslides on the valley bottom. Consistently with our results, the LGM ice thickness was found to be one of the parameters that control the occurrence of DSGSD in the European Alps (Crosta et al. 2013). In addition, Cossart et al. (2008) suggested that rock slope failures in the Upper Durance catchment (SE France) are concentrated in areas with inferred high glacial loading stress.

The slope lifecycle is strongly influenced by a less direct controlling factor as the long-term erosional imprint of glaciation. This is represented by the valley shape and slope profile. Slopes carved by glaciers on different extent (i.e. UV- and U-shaped valleys) are more prone to instability development other factors being equal, and their response time to deglaciation is significantly shorter than with V-shaped valleys. Rock slope instabilities phenomena developing in UV- and U-shaped valleys are usually smaller and deeper. Rock slope instabilities tend to cluster in correspondence of the shoulders of UV-shaped valleys, and for both U- and UV-shaped topographies, basal shear zones mostly daylight suspended over the valley floor, increasing the susceptibility to massive collapses potentially evolving into high-mobility rock avalanches.

Our results show that the development of secondary (nested) instabilities strongly depends on valley shape, with UV-shaped valley more prone to the development of smaller rockslides within a larger deep-seated gravitational slope deformation, which are key to be considered in a practical risk management perspective. Nested instabilities, highlighted by increased amount of damage, nucleate in correspondence of major slope breaks (valley shoulders). A similar evolution was recognized for the Ruinon rockslide, located in the middle-lower portion of Cresta di Saline slope, affected by the Saline DSGSD (Italian Alps, Riva 2017).

Faster valley deglaciation results in faster slope response. Indeed, perturbations after gradual load removal are expected to be lower with respect to ones occurring with faster unloading (Tao et al. 2013; Zhao et al. 2014). This finding appears to be of major interest in a context of climate change and accelerated glacial downwasting. We have to expect not only more instabilities due to recent glacier melting but these would probably develop with a shorter time lag after deglaciation. On the other hand, landslide differentiation seems to be poorly dependent on deglaciation velocity.

The presence of suspended aquifer inside the modelled landslide body confirms the enhanced sensibility of deglaciated slopes to hydrogeological inputs. This is especially evident with thicker glaciers, implying a larger decrease of confining pressure during the deglaciation, and an increase of rock mass damage. Failure events characterized by the transition from “slow to fast” are extremely destructive (Crosta et al. 2007). Despite a number of landslide inventories (Crosta et al. 2013; Agliardi et al. 2013) have been produced, there is a lack of knowledge regarding the factors and conditions leading to such transition (Hermanns et al. 2014). Analysis of the indicators of the presence of the failure surface, i.e. perched aquifers, changes in creep behaviour and slope displacement accelerations, provided the opportunity to investigate the scenarios most prone to landslide differentiation and catastrophic failures. These were linked to ticker glacier and strong glacial imprint. One of the most destructive and costly natural disaster occurred in Italy in the last decades, the Val Pola landslide, took place on the 28 July 1987, in Upper Valtellina, in a valley developed in gneiss, diorite, gabbro and ortho-quartzite. The valley is characterized by glacial trough morphology, strong relief and very steep slopes. It was covered by thick glaciers during the LGM. Anyway, in case of steep slopes and high glacier influence, also granitic slopes can undergo to disastrous failures. This is, for example, the case of the Miralago landslide (Valtellina, Italy), dated about 10–12 kyr BP.

Modern climate change and associated global temperature rise are causing increasing rates of glacial melting and thinning. Between 1850 and 1975, the European Alps have lost about half of their total glacial volume, with an additional 35–40% of the remaining amount lost since 1975 (Haeberli et al. 2007). Such a dramatic deviation from past conditions will have strong consequences on Alpine rockslides and DSGSD, which pose serious issues for human settlements. Our results suggest that faster deglaciation implies faster slope response. We expect that accelerating deglaciation rates will lead to widespread potential unstable slopes. In this context, the identification of the factors that will most likely lead current landslides to a “slow-to-fast” transition is fundamental to effectively address and define risk assessment strategies. Even if landslide differentiation seems to be poorly correlated with deglaciation rates, our parametric study defined the scenarios most prone to catastrophic instability evolution (i.e. high glacier thickness, strong glacial imprint) and gave indications about where and when slope instabilities are more probable to develop.

Conclusions

Since glacial-interglacial cycles have a strong impact on slope stability in the European Alps, paraglacial failures related to glacier downwasting is an expected outcome of current climate warming. Understanding the factors influencing probable instability evolution and failure timing is fundamental to manage slope-failure hazard in a context of accelerate glacier melting.

Using parametric numerical modelling to derive realistic deglaciation curves and to simulate their effects on different rock slope scenarios, we explore the relation between glacier thickness, deglaciation rate, valley shape, slope material and landslide development in paraglacial environments.

Our parametric study enables to link the melting of thicker glaciers to shallower rockslides daylighting at higher elevation, with a shorter slope response time following deglaciation. Shallower, suspended rockslides dominate in UV- and U-shaped valleys, where the stronger glacial imprint influences also the slope lifecycle. Weaker slopes and faster deglaciations were correlated to faster response.

In a risk-reduction perspective, the most relevant scenarios are those leading to large rock instability differentiation and involving secondary instability phenomena. These scenarios, where the collapse of the whole slope is unlikely, but the activation of the nested rockslide can cause significant damages, were recognized in valleys showing a strong glacial imprint, made of schists or gneiss covered by thick LGM glaciers. Time lags between the beginning of deglaciation and the instability evolution for these scenarios range between 30,000 and 500 years, with reduced slope response time in case of fast deglaciation.