A temperature-dependent mechanical model to assess the stability of degrading permafrost rock slopes

. In the last two decades, permafrost degradation has been observed to be a major driver of enhanced rock slope instability and associated hazards in high mountains. While the thermal regime of permafrost degradation in high mountains is addressed in several modelling approaches, no mechanical models have been developed which thoroughly explain rock slope destabilisation controls in degrading permafrost. Meanwhile, recent laboratory studies show that degrading permafrost affects both, rock and ice-mechanical strength parameters as well as the strength of rock-ice interfaces. This study presents the 10 first general approach for a temperature-dependent numerical stability model that simulates the mechanical response of a warming and thawing permafrost rock slope. The proposed procedure is exemplified with a rockslide at the permafrost-affected Zugspitze summit crest. Laboratory tests on frozen and unfrozen rock joint and intact rock properties provide material parameters for discontinuum models developed with the Universal Distinct Element Code (UDEC). Geophysical and geotechnical field surveys reveal information on permafrost distribution and the fracture network. This model can demonstrate 15 how warming decreases rock slope stability to a critical level and why thawing initiates failure. A generalised sensitivity analysis of the model with a simplified geometry and warming trajectory below 0 °C shows that progressive warming close to the melting point initiates instability above a critical slope angle of 50–62°, depending on the orientation of the fracture network. The increase in displacements intensifies for warming steps closer to zero degree. The simplified and generalised model can be applied to permafrost rock slopes (i) which warm above -4 °C, (ii), with ice-filled joints, (iii) with fractured 20 limestone or probably most of the rock types relevant for permafrost rock slope failure, (iv) with a wide range of slope angles (30–70°) and orientations of the fracture network (consisting of three joint sets). Here we present a benchmark model capable of assessing the future destabilisation of degrading permafrost rock slopes


Introduction
Rock slope failures in high mountain areas endanger human lives, settlements and infrastructure.The stability of alpine rock slopes can be considerably affected by the climate-induced degradation of bedrock permafrost (Fischer et al., 2006;Gruber et al., 2004;Gruber and Haeberli, 2007;Ravanel and Deline, 2015).Best empirical trends are provided by the well-monitored rock slope failures in the Mont Blanc Massif that have been linked to the degradation of bedrock permafrost and ice-filled joints (Ravanel et al., 2010;Ravanel and Deline, 2008;Ravanel and Deline, 2011).
Mountain permafrost warmed globally by 0.19 °C between 20070.19 °C between and 20160.19 °C between (Biskaborn et al., 2019)).Simulations of long-term permafrost evolution in Swiss, German or Norwegian mountain ranges show an overall warming of permafrost and thaw layer deepening up to the end of the century (Gallemann et al., 2017;Hipp et al., 2012;Marmy et al., 2016).Rock slope failures influenced by permafrost degradation are expected to become more frequent as a result of the warming climate (Gobiet et al., 2014;Huggel et al., 2012).
The warming of permafrost in rock slopes reduces the shear resistance along rock joints by (i) reducing the fracture toughness of cohesive rock bridges, (ii) by lowering the friction along rock-rock contacts, (iii) by altering the creep of ice infillings and (iv) by reducing the fracture toughness of ice fillings and of rock-ice contacts (Krautblatter et al., 2013).Ice fillings contribute to higher rock joint shear strengths in terms of interlocking and adhesion along the rock-ice interface and increase the stability of a rock slope.The strength of ice in fractures decreases when the temperature approaches the melting point (Davies et al., 2001;Gruber and Haeberli, 2007;Mamot et al., 2018).Several laboratory test series demonstrate the influence of warming or thawing on rock-ice-mechanical properties: The ductile temperature-and stress-dependent creep of ice and ice-rich soils has been investigated by Arenson and Springman (2005) and Sanderson (1988).Whilst the brittle failure and creep of ice-filled rock joints have been observed by Davies et al. (2000) and Günzel (2008), Mamot et al. (2018) proposed a new failure criterion that demonstrates the loss of ice friction and cohesion in rock fractures dependent on normal stress and temperature.The mechanics of frozen and unfrozen intact rock have been studied by Dwivedi et al. (2000), Inada and Yokota (1984), Kodama et al. (2013) and Mellor (1973).The friction along frozen and unfrozen rock joints was measured by Krautblatter et al. (2013).
Discontinuum discrete-element codes simulate the deformation behaviour of fractured rock masses and consider the influence of jointing (Stead et al., 2006).Among these codes, UDEC (Universal Distinct Element Code) is one of the most-used modelling tools for distinct-element modelling of jointed and blocky material in 2D and has been applied to numerous published stability analyses in unfrozen rock slopes (Bhasin and Kaynia, 2004;Fischer et al., 2010;Gischig et al., 2011aGischig et al., , 2011bGischig et al., , 2011c;;Glamheden and Lindblom, 2002;Kveldsvik et al., 2008;Welkner et al., 2010).
In this paper, we present the first numerical model that calculates the influence of warming and thawing of frozen rock and rock joints on the mechanical strength and deformation behaviour of a rock slope.The proposed model is applied to a rockslide at the permafrost-affected Zugspitze summit which is one of the highest peaks of the Northern Calcareous Alps, located in the Eastern Alps.In this context, we address the following research questions: 1) How can mechanical and thermal data from the field and from the laboratory be combined to develop a numerical stability model for a warming and thawing permafrost rock slope?
2) Is the numerical stability model of the Zugspitze summit crest capable of simulating the warming-dependent changes in stability observed at the laboratory scale?
3) How is the generic stability of a warming permafrost rock slope dependent on the slope angle and the orientation of the fracture network?

A general approach applied to the Zugspitze summit crest
A slowly moving rockslide at the Zugspitze summit crest (2900 m a.s.l.), Germany, with extensive monitoring data on kinematics and degrading permafrost bedrock was chosen to develop a close to nature rock slope instability model for degrading permafrost rock slopes (Fig. 1; see also Mamot et al., 2018).Homogeneous lithology, well known permafrost degradation history for the Zugspitze summit area and the accessibility to an ice-filled rock sliding surface (Fig. 1c) made it an ideal benchmark for the mechanical model development.The model profile and the characteristics of the fracture network for the numerical model were taken from the study site.Rock specimens collected closeby to the rockslide were used for a unique set of mechanical laboratory tests on a broad spectrum of frozen and unfrozen intact rock and rock joint properties.
Intact rock properties were converted to rock mass characteristics using internationally standardised mathematical equations.
Electrical resistivity tomography (ERT) was applied to identify current frozen and unfrozen slope sections.These were implemented in the numerical model to spatially assign frozen or unfrozen rock mass and joint parameters to the corresponding model sections and to derive a rough spatial warming pattern.Considering the numerous laboratory-based observations of the strength reduction of rock-ice-mechanical properties at warming/thawing (Table 1), we expect to model a significant effect on the slope stability.For this, we set up a discontinuum model in 2D with UDEC 7. In a first step, we modelled the numerical impact of warming and thawing on the stability of the Zugspitze summit crest with ice-filled joints.In a second step, we applied the same model to a frozen rock slope with simplified geometry to study the numerical impact of rising subzero temperature on the rock slope stability for varying slope angles and orientations of the joint set configuration.Fig. 2 demonstrates a general procedure to develop a numerical discontinuum model to calculate the influence of warming and thawing on the mechanical strength and deformation behaviour of a (partially) frozen rock slope.Here, we present one exemplary way how to apply the set-up of the model and the modelling procedure to a real-world example of alpine rock slope instability potentially affected by permafrost degradation.For instance, the detection of bedrock permafrost can also be performed by other techniques like seismic refraction tomography; similarly, the conversion of the intact rock properties to those of the rock mass can also be done by the Q-value (and not by the GSI scheme, as presented in this article).
Figure 2: General procedure for a temperature-dependent mechanical discontinuum model for the stability of a warming or thawing permafrost rock slope.The procedure is embedded as methodological approach in this manuscript and, hence, is linked to specific sections.
2.1 Characterisation of the mechanical and thermal setting

Model profile
The numerical model was set up in 2D and, thus, required a cross section of the summit crest.This profile covers a distance of 100 m, it runs from the north-face, across the crestline, to the south-face, and it crosses one of the main shear zones of the rockslide (Fig. 1a, Fig. 1b).135

Rock joint geometry and kinematic analysis
Four discontinuity sets (K1, K2, K3, K4) were identified in scan lines and field mapping (Fig. 3a).The profile for the mechanical model strikes at 146° which is in the direction of the assumed movement and follows the dip of the southern slope face (45/160).K2 (33/063) was excluded from the numerical analysis as the dip direction deviated by 66.5° from the model profile.The remaining joint sets deviating by 20-33° were implemented in the model as the standard deviations of their dip 140 directions ranged between 5-20°, falling within the tolerable range of 30° deviation (Table 2).Joint set K1 represents the bedding planes and daylights in the south-face at an angle of 24° (Fig. 3b).Note: Standard deviations are only given for measured parameters.K1-K4 = Joint sets.SZ = Shear zone.a The JRC was measured along joints not included in the scanlines.
The possible failure mechanism and shear zones which delimit the currently observed unstable rock mass were assessed manually, with a 3D point cloud of the Zugspitze summit crestline derived by UAV-based (Unmanned Aerial Vehicle) photogrammetry (Fig. 1e), and with a kinematic analysis of a potential plane and wedge failure (Fig. S1, supplementary material).The latter was conducted with DIPS 7.0 (Rocscience) for the southern slope-face and along the main shear zones identified due to field mapping and the preceding analysis of the point cloud.The instability is driven by a complex combination of a plane and a wedge failure, though planar sliding is the dominant kinematic failure mode which affects the major left part of the unstable rock mass (Fig. 1a, Fig. 1e, Fig. S1).Planar sliding mainly occurs along shear zone SZ1, while wedge failure supports the displacement along SZ1/SZ3 including a tension crack SZ2 in the upper part of the slope instability.
Two further shear zones (SZ4 and SZ5) form the assumed downslope boundary of the unstable rock mass in the lower slope sections.Planar sliding can occur along SZ4 and SZ5, while wedge failure potentially enhances the failure process along SZ3 and the planes SZ4 and SZ5.
The geophysical transect and the numerical model profile run along the stepped planar failure part of the unstable rock mass (shear zones SZ1, SZ4, SZ5; left part in Fig. 1e).Shear zone SZ1 (equivalent to K4) and shear zone SZ4 (corresponds to K1) were included in the numerical model (Fig. 3b) to allow simplification and a representation in 2D.We assume the slope instability to be significantly dominated by shear zone SZ1.Observed maximum displacements along the upper part of SZ1 lie in the range of decimetres and point to a high level of joint persistence.Further, three carst dolines with a depth of several metres were mapped along SZ1 and at the intersection of shear zones SZ1 with SZ2 (Fig. 1a, Fig. S2).SZ1 has a trace length of approximately 70 m, runs in a maximum depth of 10-15 m and dips to the southeast.It is highly fractured and displays joint opening of several decimetres filled with fine material which ranges from clay-to gravel-size (Fig. 1d).
The 3D point cloud analysis was used to assess the unstable rock mass at the south-face of the summit crest to a volume of approximately 2.9 * 10 4 m³ (Fig. 1e).Measured fracture displacements between 2013 and 2019 show that the rock mass creeps slowly at an average of 2.1 mm yr -1 (Fig. S3).Between October and May, deformation rates are 84 % lower than in the remaining summer months.Details on the volume estimation and the measurement of the crack displacements are given in the supporting material.

Spatial permafrost distribution
In the 1960s and 1990s, permafrost at the Zugspitze summit was found in deep and permanently ice-filled fractures with a maximum aperture of 0.1 m (Körner and Ulrich, 1965;Ulrich and King, 1993).Current borehole temperatures at the peak of the Zugspitze average -1.1 °C within the permafrost core area (20-25 m away from the rock walls).At the margin to the active layer (ca. 5 m away from the north-face), the temperatures approach minima of -6 °C during winter (Böckli et al., 2011;Gallemann et al., 2017;Krautblatter et al., 2010;Noetzli et al., 2010).A model for the permafrost evolution at the Zugspitze peak by Gallemann et al. (2017) shows that between 1915 and 2015 the permafrost decreased by 2 m at the north-face and by 7 m at the south-face, with an increase in the permafrost temperature of approximately 1 °C.
To accurately assess the current distribution of permafrost along the model profile at the study site, we applied laboratorycalibrated electrical resistivity tomography (ERT).This technique has been used to characterise and monitor the spatial variability and evolution of mountain bedrock permafrost in steep rock walls (Keuschnig et al., 2017;Krautblatter et al., 2010;Murton et al., 2016).The outlined procedure of the measuring setup in the field was similarly applied to the Zugspitze summit crest.Two surveys were conducted at the end of August 2014 and 2015, which corresponds to the time the active layer at the Zugspitze peak usually reaches the maximum extent (Gallemann et al., 2017).Data processing and inversions were performed with Res2Dinv (see supplement material for measuring setup, data acquisition and analysis of the field surveys).
To assign the measured ERT of the summit crest to frozen or unfrozen slope sections, we used the resistivity of saturated frozen and unfrozen Wetterstein limestone in the laboratory, defined (i) by Krautblatter et al. (2010) in the same lithology, and (ii) in the context of this study based on two rock blocks from the field site.For the laboratory tests, we followed the suggested procedure by Krautblatter et al. (2010), though limited to a single cooling and subsequent freezing trajectory of the samples from 10 down to -6 °C.The rock samples (both with a size of 0.02 m³) were kept in a self-built, temperature-controlled cooling box for the duration of the tests.Rock temperature was measured simultaneous to resistivity with two Pt100 sensors (Greisinger GMH3750, with a 0.03 °C precision) inserted in the rock samples (Fig. 4a).Measurements were performed along 3-4 different Wenner arrays installed parallel to each other.
The laboratory ER clearly showed values lower than 19 kΩm for unfrozen rock and values higher than 28-29 kΩm for frozen rock (Fig. 4a).The equilibrium freezing point at -0.2 to -0.4 °C is indicated by a 21 to 37-fold increase of the frozen temperature-resistivity gradient.This pattern of electrical resistivity versus rock temperature is similar to the one demonstrated by Krautblatter et al. (2010): here, the ER measures 30 ± 3 kΩm at the equilibrium freezing point (-0.5 °C).Overall, we defined ER-values lower than 19 kΩm and higher than 33 kΩm to be unfrozen or frozen, and the range of values lain in between as possibly frozen rock.The ERT results of 2014 and 2015 are interpreted as follows (Fig. 4b and Fig. 4c): -The shallow rock layer is almost continuously unfrozen on both aspects of the crest.This assumption is confirmed by (i) borehole rock temperatures at the Zugspitze peak which show a 2-4 m thick active layer on the north-face (Gallemann et al., 2017), and (ii) thermistor-based rock temperatures measured along the ERT transect in a depth of 10-80 cm below the ground surface.The corresponding thermistors (specially modified iButtons in rods) show that the rock temperatures do not remain below zero degrees throughout the year and are positive at the moment of the ERT surveys (Fig. 1b, Fig. S4).
-Due to this evidence, the current boundary to permafrost on the north-face and at the upper part of the south-face was set to a continuous depth of approximately 5 m (black line in Fig. 4b and Fig. 4c).
-Most near-surface interstitial patches of resistivity > 33 kΩm relate to massive ice visible in the meter deep carst doline close to the crestline and surficial drying of rock usually limited to the upper decimeter (Hauck and Kneisel, 2008;Sass, 2005).
-On the lower south-slope, the unfrozen bedrock becomes thicker with decreasing altitude and can be found down to a depth of 10 m and more indicated by resistivities < 19 kΩm.
-Other individual patches of low electrical resistivity (< 19 kΩm) refer to major shear zones or to the highly porous talus slope below the northern rock wall filled with air or with unfrozen, saturated and fine material (Hauck and Kneisel, 2008;Supper et al., 2014).Zones of high conductivity can be promoted by precipitation or water from snow melt which easily infiltrate along unfrozen fine fillings or fractures (Keuschnig et al., 2017).
-The resulting spatial extent of permafrost relates to the core of the crest and rock layers below 5 m at the north-face.

Mechanical properties
The Zugspitze summit crest consists of Triassic Wetterstein limestone which has a thickness of 1000 m and constitutes the majority of the Zugspitze massif and the Wetterstein mountains (Gallemann et al., 2017;Miller, 1962).The limestone is massy, fine-grained and dolomised and shows little heterogeneity in terms of lithological properties (Krautblatter et al., 2010).Its porosity ranges between 1.9 % (Draebing and Krautblatter, 2012) and 4.4 % (Krautblatter et al., 2010).
The distinct element method in UDEC enables to simulate the mechanical behaviour of a discontinuous medium: the assemblage of discrete rock blocks separated from each other by discontinuities.We extensively investigated frozen and unfrozen rock and joint properties in the laboratory to assess the effect of thawing on the mechanical strength and deformability, and to determine the required material parameters for the model.All tests were performed with fully water saturated Wetterstein limestone samples collected at the study site.Frozen rock samples had a mean temperature of -5 °C, while unfrozen rock specimens were tested at room temperature (~22 °C).Detailed information on the preparation of the rock samples, the test setups and conditions can be looked up in the supporting material (see Sect.S6).
The strength and deformability of the model blocks and zones within the blocks were represented by the Mohr-Coulomb plasticity constitutive model with a tension cut-off.At the slope scale, material properties of the rock mass instead of intact rock are needed as discontinuities reduce rock mass strength and induce non-linearities and anisotropy in the stress-strain behaviour (Hoek et al., 2002;Wyllie, 2018).Hence, blocks and zones were allocated rock mass properties.The intact rock data were transferred to rock mass characteristics and switched to field-scale properties following the GSI (Geological Strength Index) rock mass classification system by Hoek and Brown (1997).The GSI system is directly linked to Mohr-Coulomb or Hoek-Brown strength parameters and rock mass deformation moduli (Cai et al., 2004).The deformation modulus of the rock mass Em can be described based on a field survey of the rock mass including description and measurement of block volumes and joint characteristics (Hoek et al., 2002): � �σ  /100 * 10 ((−10)/40) Eq. ( 1) where GSI is the geological strength index for the rock mass, D is the disturbance factor and   is the uniaxial compressive strength of the intact rock.The latter was tested in the laboratory and is a function of temperature (Table 3).Following Cai et al (2004), a GSI value of 65 was derived which relates to an interlocked, mostly blocky to very blocky rock mass formed by four joint sets with a mean spacing of 56 cm (Table 2), and rough slightly weathered surfaces.The disturbance factor D was set to 0 as the near surroundings of the profile are not affected by significant blast damage or excavating actions.
Values for the elastic rock mass shear and bulk moduli Gm and Km were derived according to the equations presented by Tipler and Mosca (2004) for intact rock: Eq. ( 2) and Eq. ( 3) where the Poisson's ratio υ, measured in this study (Table 3), describes the ratio of the transverse strain to the longitudinal strain under conditions of uniaxial stress (Jaeger et al., 2007).For Gm and Km, the dynamic Young´s modulus Edyn was substituted by Em (see Eq. ( 1)).The mass density ρ (g/cm³) and dilatational wave velocity VD (m/s) were measured in additional tests (Table 3) and utilised to calculate the dynamic Young's modulus   for intact rock (Rentsch and Krompholz, 1961): Eq. ( 4) Edyn and G were later used to determine the joint stiffness parameters.The rock mass internal friction angle φm and cohesion cm were estimated with the GSI and the Hoek-Brown constant mi according to Cai et al. (2004).Consider, while rock mass cohesion could be determined for both thermal conditions (as it is based on the temperature dependent uniaxial compressive strength), the unfrozen friction angle value was used for all thermal conditions (as the respective relation is not temperature dependent).The uniaxial tensile strength σtm was calculated with the following relation Eq. ( 5) where mb and s are constants for the rock mass (Hoek and Brown, 1997).
All mechanical properties of the intact rock and the rock mass, either tested in the laboratory or calculated, decrease temperature-dependently with 0.3 to 0.8 % °C-1 (absolutely 8-22 %) upon thawing (see Table 3 and Table 4); this pattern corresponds with other mechanical studies (see Table 1).While the increase in deformability is more pronounced for intact rock (0.8 % °C-1 ; 20.4-20.5 %) than for rock mass (0.3 % °C-1 ; 8.3-8.9 %), the tensile strength decreases by a similar amount when comparing intact rock with rock mass.This difference can be explained by the small variation in frozen and unfrozen Em-values caused by the similar frozen and unfrozen intact rock uniaxial compressive strength (see Eq. ( 1)).Note: a The frozen rock mass friction angle was given the same value as for the unfrozen friction angle, resulting in a decrease in thawing of zero.
In UDEC, the deformability and strength of rock discontinuities were simulated by the Mohr-Coulomb area contact constitutive model.Rock joint parameters for ice-free and ice-filled joints were considered in the presented stability analyses.The deformability of the joints is described by the joint normal stiffness kn: Eq. ( 6) and the joint shear stiffness ks: Eq. ( 7) where s is the joint spacing (m), Em and E are the Young's moduli (GPa) for rock mass and intact rock respectively (Barton, 1972) and Gm and G are the shear moduli (GPa) for rock mass and intact rock respectively (Glamheden and Lindblom, 2002;Itasca Consulting Group, 2019).Finally, temperature () dependent cohesion c  [] = 104.5 − 143.8 *  Eq. ( 8) and friction φ tan [] = 0.19 − 0.1 *  Eq. ( 9) of frozen ice-filled rock joints are determined based on the failure criterion presented by Mamot et al. (2018) which is valid for temperatures ranging from -4 to -0.5 °C and normal stresses between 100 and 800 kPa.Therefore, the cohesion and friction angle were calculated for temperatures -4, -3, -2, -1 and -0.5 °C, which are currently measured in the frozen sections of the Zugspitze summit (Gallemann et al., 2017).The stress range of 100-800 kPa corresponds to the rock overburden represented by the model of the Zugspitze summit crest (Fig. 3b).
Values for the cohesion of frozen and unfrozen, ice-free rock joints were roughly estimated from (Krautblatter et al., 2013).
As initial displacements lead to the destruction of asperities and smoothing of joint surfaces, the frictional strength of rockrock contacts was represented by the residual friction angle φr, estimated following Barton and Choubey (1977): where φb is the basic friction angle, r is the Schmidt hammer rebound value of weathered surfaces, and R is the Schmidt hammer rebound value of unweathered, sawn surfaces.Rebound hardness values of weathered frozen and unfrozen Wetterstein limestone surfaces were collected with the Schmidt hammer (N-type) following the proposed method by Aydin et al. (2005) and Ulusay (2015).In contrast, the basic friction angle was determined by tilt tests of frozen and unfrozen joint surfaces as recommended by Barton (2013).
The investigated and calculated joint deformation and mechanical strength properties reduce absolutely by 2 to 22 % upon thawing (0.1-0.8 % °C-1 ; Table 5).The small variation between frozen and unfrozen Em-values results in a small difference between frozen and unfrozen kn and ks (0.1 % °C-1 ; 2.6-2.8%).In contrast, the stiffness values vary more strongly between the various joint sets which are also dependent on the joint spacing.In general, the estimated unfrozen kn and ks lie well within the range of values proposed by the UDEC database (Itasca Consulting Group, 2019) and Kulatilake et al. (1992), or measured by Barton (1972) and Bandis et al. (1983).Similarly, the measured unfrozen φb corresponds well to the values listed e.g. in Barton and Choubey (1977).10), calculated after Barton and Choubey (1977) x

Model setup
The topography of the rock slope used in the model was derived from a digital terrain model of the Zugspitze summit area.
Three fully persistent joint sets (K1, K3, K4) and the dominant shear zone were included (Fig. 3b).For simplification and shorter computation time, the joint spacing was scaled up in the model by factor five.The presence of intact rock bridges was accounted for according to the approach by Jennings (1970).After this, the joint cohesion and friction angle were proportionately increased by the respective rock mass cohesion and friction angle dependent on the estimated percentages of the surface areas of rock bridges and joints within the rock slope.
The blocks in the model are supposed to be deformable and were subdivided into a mesh of finite-difference elements with a size of 1.5 m.The stress field was initialised according to the varying density of overburden which depends on the topography of the slopes.We assumed a horizontal to vertical stress ratio of 0.5.Roller boundary conditions were implemented at the sides of the model (i.e.vertical movements were allowed), whereas vertical and horizontal displacements were suppressed along the base.
The combined boundary of bedrock permafrost within the summit ridge, derived from ERT measurements of 2014 and 2015 (Fig. 4b), was used to define the current frozen and unfrozen sections of the numerical model.Then, we introduced six subsurface layers to simulate a stepwise warming or thawing from the slope surface to the core of the crest by adjusting the temperature-dependent material parameters for a specific temperature level and subsequent numeric cycling.The layers are oriented parallel to the derived permafrost boundary and account for a stronger warming signal directed from the south-slope (Fig. 5).The defined current permafrost boundary and the estimated spatial pattern of layers for warming are both in accordance with modelled current and future thermal fields of arbitrary mountain ridge geometries (Noetzli et al., 2007) and, in particular, of the Zugspitze (Böckli et al., 2011;Noetzli, 2008).(current state of permafrost distribution in the crest).

Procedure for the numerical analysis
The general modelling procedure consists of four steps: the calculation of an initial equilibrium and three successive steps of potential destabilisation which include the progressive destruction of intact rock bridges, the warming and the thawing of a permafrost rock slope: 1) Initial state: The numerical analysis was started by developing an initial equilibrium with a low joint persistence of 30 % (Model stage 1 in Table 6).The summit crest is assumed to be fully frozen representing an undefined moment in the past.
2) The progressive loss of rock bridges: The destabilisation of a permafrost-affected rock slope is initiated by the progressive destruction of cohesive intact rock bridges which prepare new shear planes along which displacement can take place (Krautblatter et al., 2013).The current overall joint persistence was assessed to 90-100 % as we observed a high trace length of the discontinuities and a displacement of several decimetres along the main shear zone indicating that most of the rock bridges are lost.The past destruction of rock bridges was simulated by a reduction of the apparent joint cohesion and friction angle in three steps which represents a progressive increase of the joint persistence from 30 to 90 % (Stages 1-4 in Table 6, Table 7; for more details see supplement, Sect.S7).
Joints may be partly filled with ice during the loss of rock bridges.However, joints were assigned properties of icefree joints as, according to Krautblatter et al. (2013), deformation and shear strength are mainly controlled by rockmechanical processes at this early stage of destabilisation.As long as a certain part of intact rock bridges and asperities along joint surfaces are present, stresses are supposed to concentrate at these locations while the softer ice fillings may be squeezed away (in particular, in greater depths).
Following Krautblatter et al. (2013), the advanced stage of accelerated displacements in a permafrost-affected rock slope begins as soon as most of the rock bridges are broken.Then, ice mechanics (fracturing of rock-ice contacts and ice, as well as the creep of ice) will increasingly replace rock mechanics (friction of rock-rock contacts and fracture of rock bridges) in controlling displacements along discontinuities and their potential acceleration.To introduce this stage of destabilisation, the joint persistence was further increased to 100 % and rock joints were simulated to be icefilled (Stage 5 in Table 6).For simplification and due to data availability, the joint shear resistance was solely given by the fracture of ice or rock-ice contacts, while the creep of ice was neglected.
3) Warming permafrost: The frozen crest was gradually warmed from the slope surface to the core from -4 to -0.5 °C (Stages 6-9 in Table 6).This was applied to the six rock layers defined in Fig. 5. Firstly, the outermost (orange) section was warmed to -3 °C, while inner sections remained at -4 °C (Stage 6).After that, the outermost section was warmed to -2 °C and the adjacent inner section (yellow layer in Fig. 5) was warmed to -3 °C, while inner sections remained unchanged (Stage 7).This procedure was continued until the outermost layer was at -0.5 °C (Stage 9).Each warming step was characterised by adjusting the temperature-dependent material parameters to the warmer temperature and by subsequent numeric cycling.The simulated temperature range of the rock mass is consistent with currently monitored borehole temperatures in permafrost rock walls across the European Alps (Gallemann et al., 2017;Noetzli et al., 2019).
4) Thawing from the slope surface to the core was implemented in the same way as warming (Stages 10-15 in Table 6), until all subsurface layers of the summit ridge were unfrozen.The spatial pattern of thawing is illustrated in Fig. 5.
To assess the level of a potential slope stability loss due to warming or thawing, we calculated the factor of safety (FS) by using the common strength reduction technique (Wyllie, 2018).Thereby, the cohesion and frictional strength of the rock mass and the joints are reduced simultaneously and gradually by increasing trial factors in a series of numerical simulations until failure occurs.A bracketing solution approach was applied to progressively reduce the bracket between stable and unstable solutions until it falls below a specified threshold (Itasca Consulting Group, 2019).The resulting FS is a single indicator of minimum stability which globally refers to the entire slope.3 Sensitivity analysis for a simplified, warming permafrost rock slope The relation between the dip of the joint sets, especially the bedding, and the slope-face is a crucial factor for rock slope stability (Cruden, 2003;Wyllie, 2018).Thus, we performed a sensitivity analysis on the numerical impact of varying (i) slope angles and (ii) orientations of the fracture network on rock slope stability.For this, we transferred the Zugspitze model to a rock slope with a simplified topography and warming steps were applied to the entire model domain without spatial differentiation.The topography of the Zugspitze south-face was modified to a straight line.In accordance with the dimensions of the crest geometry, the height of the slope and the width of the upper-face were standardised to 84 m and 32 m, respectively (Fig. 7).
Modelling was applied to 12 different slope angles between 30 and 69°.For inclinations of 70° and higher, no initial equilibriums could be calculated as cycling exceeded the maximum computation time.The joint set orientations were kept constant.In a first step, an initial equilibrium was calculated for a frozen rock slope with ice-filled joints and a joint persistence of 90 %.In a second step, the joint persistence was increased to 100 %, and then the rock slope was warmed in four steps from -4 to -0.5 °C (by reducing the rock joint strength equivalent to the Zugspitze model).However, each warming step to the next degree centigrade was applied to the entire rock slope and not just to a single rock layer.Reducing the complexity of the model facilitates the transferability to other frozen rock slopes.
We also remodelled the stability for 12 different dip angles of the fracture network with reference to the slope.The positioning of the three model joint sets to each other was held constant, while rotating them counter-clockwise in steps of 15° for each of the 12 slope angles used for the first sensitivity test.The analysed orientations of the fracture network represented anaclinal and cataclinal slopes after Cruden (2003;Fig. 3).Hereafter, cataclinal slopes have a slope-face and bedding with the same dip direction whereas, for anaclinal slopes, the dips of both planes have the opposite direction.The orientation and dip of the fracture network was represented by the bedding and then applied to the scheme after Cruden (2003).The simulated temperature of the rock mass was constantly at -4 °C since a stepwise warming was not applied.

Stability of the warming/thawing Zugspitze crest
Figure 6 shows the cross section of the Zugspitze summit crest with a factor of safety for different model states.While an increase in joint persistence or warming hardly lead to a rising block displacement magnitude (indicated by colours, note that a-c and d-f use different scales), thawing causes a drastic increase in displacements.As soon as thawing starts, the south-face becomes unstable and maximum displacements increase to 0.14 m (Stage 10, Fig. 6d).Rock deformations rise to 0.42 m when the dimension of thawed bedrock reaches the current state (Stage 11, Fig. 6e), and measure 0.78 m for the whole rock slope thawed (Stage 15, Fig. 6f).The highest deformations (> 0.4 m) concentrate on the upper 5-7 m of rock in the most inclined lower part of the slope (with a mean angle of 49°).While warming brings the factor of safety close to unity, it lies below 1 for all Stages 10-15 which refer to thawing.

Sensitivity to the slope angle
The model results for exemplary slope inclinations between 60-64° are shown in Figure 7.The block displacements increase and the related factors of safety decrease with higher slope angles and temperatures below the melting point.The factor of safety falls below unity at temperatures above -1 °C when a rock slope with a slope angle of 62° or more is warmed up.A 60° steep rock slope becomes unstable at temperatures above -0.5 °C.The interaction between slope gradient, maximum displacements and temperature can be separated into three domains (Fig. 8a): 1) Domain 1 (30 to 50°): The maximum displacements increase slightly by 0.4 mm per degree slope angle.The correlations for all temperatures are significant with p-values ≤ 0.04.
This becomes obvious due to the abrupt increase in displacements and the p-value > 0.05.For rock slopes at -1 °C and below, the onset of instability is visible above an inclination of 55°.Note that all displacements are one magnitude higher than in Domain 1.
3) Domain 3 (64 to 69°): The displacements increase suddenly with a further increase in the p-value and a corresponding decrease in the R²-value for a slope angle of 64° and for all rock temperatures.This points to a second stage of accelerating slope destabilisation.The displacements are one magnitude higher than in Domain 2.
The displacements are sensitive to the sub-zero temperatures at all slope gradients.However, the intensity of rock slope destabilisation is increased by higher rock temperatures close to the melting point (Fig. 8b): A warming from -4 to -2 °C leads to an increase in absolute displacements, being low in Domain 1 (0.2 mm °C-1 ), but high in Domain 3 (7 mm °C-1 ).A respective warming from -2 to -0.5 °C shows also an increasing trend with higher slope gradients, but is more pronounced with a massive increase from Domain 1 (0.6 mm °C-1 ) to Domain 3 (45 mm °C-1 ).Nevertheless, the highest temperature-dependent relative increase in displacements is observed in Domain 2 (Fig. 8b): Values increase from 9.2 % °C-1 (Domain 1) to 38.5 % °C-1 (Domain 2) and fall to 25 % °C-1 for slopes steeper than 62° (Domain 3), which is valid for a warming from -2 to -0.5 °C.
Again, the same pattern is less pronounced for a warming from -4 to -2 °C.
Comparatively to the displacements, the factor of safety also depends on the slope angle and the temperature (Fig. 8c): For a rising steepness from 30 to 69°, the calculated FS decreases inversely and is always lower for higher temperatures.Warming from -4 to -0.5 °C reduces the FS of a permafrost rock slope by a mean factor of 3.3.The pattern of domains in Fig. 8a can be easily applied to the relation between FS, temperature and slope angle (Fig. 8c): The mean relative reduction of the FS per degree of the slope angle is low in Domain 1 (-1.8 % ± 0.4) and becomes higher over Domain 2 (-4.5 % ± 3.3), until it reaches a maximum in Domain 3 (-9.8% ± 8.3).Correspondingly, rock slopes are stable for all modelled temperatures in Domain 1.
When shifting to Domain 2, the FS falls below unity at -0.5 °C and at -1 °C which coincides with an abrupt increase in displacements.Rock slopes at -4 °C keep stable over the whole range of slope gradients.
As we did for the maximum block displacements, we tested if warming closer to the melting point led to a distinct diminishing pattern of the FS (Fig. 8d): The absolute decrease is very similar for both warming steps from -4 to -2 °C and from -2 to -0.5 °C.However, the relative decrease from -4 to -2 °C is always lower than the one at higher temperatures.

Sensitivity to the fracture network orientation
Modelling resulted in 8 out of 12 orientations of the fracture network that show a significant displacement acceleration and initiating slope instability above a critical slope inclination (Fig. 9a): For anaclinal slopes, the critical angles range between 50-62° and, for cataclinal slopes, they range between 55-62°.Cataclinal rock slopes with a bedding dip of 84 and 69° do not indicate a significant knickpoint in the displacement curves.As such, the corresponding critical slope angles were set to the maximum of the investigated values (69°).Fracture networks rotated counter-clockwise by 30° and 75° from the joint set of the Zugspitze summit crest showed very noisy model data and were not used for the analysis.

Discussion
This study successfully simulated the mechanical response of a warming and thawing permafrost rock slope using a 525 temperature-dependent numerical model.The following sections show the general limitations of the presented models and critically discuss the site specific model, the model of the generalised permafrost rock slope and its transferability to other high mountain conditions.

Limitations
Mechanical models are often sensitive to material parameters.However, the availability of these parameters and the technical opportunities to establish them are limited.Therefore, the presented mechanical model of the Zugspitze summit crest bases on a few simplifications.First, it is crucial to consider the temperature dependency of the material parameters as freezing can increase rock mass strength.However, the influence of freezing effects is mostly ignored in the GSI rock mass classification system (see also Gambino and Harrison, 2017).In our model, the friction angle after Cai et al. (2004;Fig. 5) applies the temperature-independent Hoek-Brown constant mi and the GSI value, defined for the unfrozen rock mass.To our knowledge, the friction angle of a frozen rock mass and its change due to thawing is hardly investigated and not yet described by a general law.Model test runs with very low to unrealistically high values of the mechanical rock mass properties demonstrated that the displacements in the rock slope mostly remained within the same order of magnitude (Fig. S9).Therefore, the rock mass internal friction angle was derived due to the GSI scheme, and its unfrozen value was assigned to both the frozen and the unfrozen rock mass.Second, the knowledge of the subsurface structure (i.e.joint spacing and joint persistence) is limited, and has to be estimated due to inspection of the rock surface.Third, mechanical parameters were not determined for the specific joint sets.Only the joint normal and shear stiffness varied between the joint sets due to their different joint spacing.Fourth, a standardised recommendation for the disturbance factor D for an unstable rock slope (without blast damage) does not exist; such a recommendation should be validated numerically and include a classification of various degrees of instability.Therefore, D was set to 0, but the potential influence of a higher D on the slope stability was analysed due to additional model runs with the geometry of the Zugspitze crest (Fig. S6); a modified D approaching 1 leads to higher displacements but the factor of safety remains unchanged.
Furthermore, the implementation of the numerical model implies a few limitations which are addressed with the following approaches: -When implementing non-persistent joints, UDEC discards those which do not split a block.An irregular joint and block pattern may be formed which does not correspond to reality.Therefore, discontinuities were constructed fully persistent, but joint strength parameters were increased according to Jennings (1970) to omit the potential result of favoured higher displacements and velocities along persistent joints.
-The size of the Zugspitze model domain was constrained to a local section of the rock slope ending in a steep part at the south-face.The model domain fully incorporates the unstable rock mass under investigation and thereby perfectly accounts for the simulated spatial conditions of the rock-ice-mechanical processes.To investigate the slope angle at the domain boundary, we rerun the model with a bigger domain size ending in a flatter section of the south-slope (Fig. S6): the results demonstrated that a change in the size of the model domain does not significantly change the overall stability of the slope.
-The initial temperature of the frozen rock slope does not affect the displacement magnitudes and the general pattern of dependency between maximum displacements and slope inclination.The potential effect of a lower strength on the calculation of the initial equilibrium and on the resulting maximum displacement values was proven by reiterating the generalised model with an initial temperature at -2 °C and comparing it with the original -4 °C (Fig. S7).
-Warming to the next degree centigrade was, by default, calculated with 3000 cycles which was a compromise between a representative number of cycles for potentially reaching a numerical equilibrium (for stable rock slopes), and a reasonable, not too time-consuming calculation.Computations for low/intermediate slope angles and temperatures reached an equilibrium prior to 3000 cycles, while calculations for higher slope angles and temperatures mostly failed to reach an equilibrium (for unstable rock slopes) or scarcely required > 3000 cycles.Hence, the numerical calculation was repeated with warming steps with 6000 cycles to test a potential cut-off effect of cycles required for the rock slope to react to a change in stress-strain due to a modification of material parameters.However, the results coincide well with the model runs with 3000 cycles since displacements significantly accelerate above slope angles of 50° and 62°, and they are higher for temperatures closer to melting (Fig. S8).

Stability assessment of the warming/thawing Zugspitze summit crest based on the factor of safety
The numerical model for the Zugspitze summit ridge simulates the strength reduction at warming and thawing based on laboratory tests of this and of previous studies (Table 1).Warming of ice-filled joints takes the slope close to a critical level of stability (FS = 1.3), but it does not initiate instability.As soon as thawing sets in and the ice is lost, the south-face becomes massively destabilised as displacements increase by two orders of magnitude and the factor of safety falls below unity.
Destabilisation is triggered by a slight reduction of the joint cohesion from 0.18 MPa (at -0.5 °C) to 0 MPa (Table 5), even though the joint friction angle is increased from 13.5 to 29.3°.Davies et al. (2001) postulate that the FS of a permafrost affected rock slope is higher for unfrozen joints than for ice-filled joints between approximately -1.5 °C and zero degrees.In contrast to that, our numerical analysis of the thawing Zugspitze summit crest shows that many ice-filled fractures reach a minimum stability only upon thaw.The reduction in joint stability is significantly influenced by a related loss of the joint cohesion.The differing observations on the stability upon thaw may be explained as follows: (i) Our model assumes a lower joint roughness and uses the residual friction angle which is lower than the peak friction angle taken for the centrifuge model.We suppose that the simulated rock slope has already been moving for a long time and joint surfaces are flattened due to the progressive destruction of asperities along joints.In contrast, the sliding surface in the centrifuge model is not affected by any previous displacements.(ii) The calculated FS by Davies et al. (2001) corresponds to an independent moment in time, while the FS in the Zugspitze model is path-dependent and affected by the history of the previous numerical calculations.
The permafrost evolution model by Gallemann et al. (2017) predicts permafrost at the Zugspitze peak to be lost within the second half of the 21 st century which may result in an elevated occurrence of rock slope instability.This scenario is related to a projected increase in the mean air temperature by 3.2 °C until the end of the century.Stages 11 to 15 of the mechanical model represent the same time period (current state until second half of the century) and simulate progressive thawing of the summit crest that leads to the continuation of increasing displacements and the preparation of final slope failure (Fig. 6).
The factors of safety for all stages with one or more thawed rock layers, including the current extent of permafrost (Stage 11, Fig. 6e), lie below unity which indicates slope instability.This is consistent with the observations of the current displacements and characteristics of the main shear zone at the field site.However, the high values calculated by the model are not comparable with the measured displacement rates at the field site (2.1 mm y -1 ) for the following reasons: (i) The numerically calculated displacements are maximum values and refer to local sections in the model.(ii) Cycling in UDEC was not linked to specific time-steps.As such, we could not assign specific periods of time to the distinct, non-linear steps of warming or thawing.This did not permit us to estimate the approximate time period required for the calculated displacements to occur.(iii) Absolute displacements are path-dependent and affected by the history of calculations, including input and alteration of parameters as well as the stress-strain behaviour in the system.
Based on the 3D kinematic analysis (Sect.2.1.2),we postulate that the most likely mechanism controlling the slope instability is a combination of a plane and a wedge failure.However, for the 2D numerical model and at the scale of the failure plane, we consider various possible complex mechanisms: (i) Sliding of several rock blocks along a polygonal shear zone SZ1/SZ4 that is constituted by joint sets K1, K3 and K4.This complex type of failure is provided by antithetic fractures which enable shear displacements between the moving blocks (Eberhardt et al., 2004;Poisel and Preh, 2004).(ii) Shear or tensile failure along joint sets K3 or K4 may favour shear displacements along K1 by supplying additional stress onto the blocks and inducing "step-path" failure (Huang et al., 2015;Mejía Camones et al., 2013).(iii) Failure at the scale of single blocks can be induced by forward block toppling with K1 and K3 or backward block toppling with K1 and K4 (Itasca Consulting Group, 2019).Both options will only work for higher columns subdivided by surfaces of K1 which stick together due to locally higher frictional or cohesive strength.(iv) The two joint sets K1 and K4 favour the displacements to concentrate on the south-face.The highest deformations concentrate on the lower slope section with a mean gradient of 49°, but they are not coupled to the prominent shear zone.This can be the result of (i) the dominating influence of the slope angle on the stability of the system and (ii) the lack of variation in the joint properties of the shear zone and the remaining joint sets (Table 5).

The stability of a simplified permafrost rock slope with rising temperature
The destabilisation of a permafrost rock slope with ice-filled joints and a fracture network and rock/joint properties similar to the Zugspitze south-face can be initiated at a slope gradient > 50 or 55° (Domain 2 in Fig. 8: "First onset of instability") and significantly favoured above a gradient > 62° (Domain 3 in Fig. 8: "Accelerating slope destabilisation").The transition from Domain 1 (Fig. 8: "Displacements increase slightly with the slope angle") to Domain 2 seems to be the most critical in terms of slope destabilisation, as it is characterised by the highest relative increase in displacements per degree of the slope angle.
Interestingly, this could be confirmed by the study of the critical slope angle for varying orientations of the fracture network (Fig. 9): 50 % of the studied orientations have a critical slope angle of 50 or 55°, which corresponds to the transition from Domain 1 to Domain 2. 30 % of the orientations lead to instability with angles of between 57.5-62°, and 20 % remain stable within the studied range of slope gradients.The relevance of slope angles at the described transition from Domain 1 to Domain 2 (Fig. 8), identified by the numerical model, is in accordance with documented rock slope failures in the Swiss Alps, adjacent areas in France and Italy (Fischer et al., 2012) and in the Mont Blanc Massif (Ravanel et al., 2010), which are mostly attributed to the degradation of bedrock permafrost: The mean slope gradient of the detachment zones in the Mont Blanc Massif was 54°, while the highest percentage of detachment zones in the Swiss Alps and adjacent areas had slope angles of between 40-60°, leading to a postulated lower slope gradient threshold of 40-45° for rock slope failures.
We further demonstrate that the factor of safety falls below 1 when a rock slope with a gradient ≥ 50° is warmed from -4 °C to a temperature between -3 and -0.5 °C.This sensitive temperature range corresponds well to (i) the temperatures currently monitored in boreholes in permafrost rock walls in the European Alps (Gallemann et al., 2017;Noetzli et al., 2019) and (ii) the temperatures which are characteristic for warm permafrost areas or the lower permafrost boundary (-5 to 0 °C).The latter posed the release zones for most of the rock slope failure events documented by Fischer et al. (2012) and Ravanel et al. (2010).
The model results also demonstrate that the increase in displacements becomes more pronounced when approaching the melting point which indicates a non-linear relation (Fig. 8b): the relative and the absolute increase are always higher for a warming from -2 to -0.5 °C than for a warming from -4 to -2 °C.The relation between rising subzero temperature and reduced stability was also observed by Davies et al. (2000Davies et al. ( , 2001)).However, on the basis of our model results, we propose to complement the conclusions by Davies et al. (2001) in the following way: (i) The temperature of critical stability is highly dependent on the inclination of the main fractures versus the slope angle.(ii) It is of crucial importance to consider the stress conditions in a rock slope when extrapolating the results from the laboratory scale to the field scale.(iii) The use of rock instead of concrete samples for the laboratory determination of joint parameters results in a closer reproduction of the real conditions along rock joints in the field.

Transferability and implications for other field sites
The general procedure for a temperature-dependent mechanical stability model (Fig. 1) can be applied to any warming or thawing permafrost rock slope.However, the results of the sensitivity analysis with input data from the Zugspitze summit ridge (Sect.4.2) are valid for rock slopes with increasing subzero temperatures, consisting of fractured limestone, with strength and deformability similar to the Wetterstein limestone tested in this study.A transfer to permafrost rock slopes with a different lithology is possible, but it requires more laboratory calibration tests and modelling.A provisional transfer may be possible with the following justification: Mamot et al. (2018Mamot et al. ( , 2020) ) studied the shear strength of ice-filled permafrost rock joints and developed a resilient temperatureand stress-dependent failure criterion which is valid for the majority of the rocks that are relevant for observed permafrost rock slope failures in the Alps.Since this failure criterion was used to calculate the cohesion and the friction of the ice-filled joints in the presented model, the numerical description of them may be applied to other lithologies, too.However, the joint stiffness and the mechanical parameters of the rock mass still vary among different rock types and may lead to different model results.
To get a first impression of this potential effect, we performed a couple of model test runs with varying values for the mechanical properties of the rock mass ranging from very low to unrealistically high.As the displacements remained within the same order of magnitude (Fig. S9), we infer that the model results will not be significantly influenced by other rock types with distinct properties of strength and deformability.Nevertheless, we strongly recommend to verify the effect of different rock types by thorough modelling in the future.The transferability of the numerical results is also limited to rock slopes with a similar joint set configuration as in this article since the rule-of-thumb for critical slope angles and fracture network orientations (Fig. 9) is valid for a specific constellation of three joint sets and not only for the bedding, as presented by Cruden (2003).Considering the constraints above, the simplified model, developed by the sensitivity analysis (Sect.4.2), can be transferred to warming permafrost rock slopes with ice-filled discontinuities and slope gradients smaller than 70°.In addition, these rock slopes consist of limestone (and probably mostly all rock types relevant for permafrost rock slope failures in the Alps), and they contain three joint sets separated from each other by an angle of 45°.

Conclusion
This study presents the first numerical model which is capable of performing comprehensive mechanical stability analyses of degrading permafrost rock slopes.In this context, we provide (i) a universal procedure for the input of thermal and mechanical data, the model setup and the modelling strategy, (ii) a numerical benchmark application to a specific test site at the Zugspitze peak, and (iii) the first numerically derived critical stability thresholds related to the slope angle, rock mass temperature and orientation of the fracture network.The related main findings are summarised as follows: (i) The temperature-dependent mechanical stability model can be used for any permafrost-affected rock slope across the globe which is subjected to climatic warming.
(ii) Laboratory tests and field reconnaissance of the benchmark site Zugspitze exemplify thermal, geometrical and mechanical input data for the numerical model.
(iii) Frozen and unfrozen bedrock material properties were assigned to specific sections in the model and changed due to warming/thawing.The modelling procedure was divided into three stages: rock bridge destruction, warming with ice-filled joints, and thawing.Process-specific and temperature-dependent input parameters were modified when switching from one stage to the next.(vi) The dependence between instability and the slope angle can be classified into a stable first domain (≤ 50°), a second domain with a first onset of instability (55-62°), and a third domain characterised by an accelerated slope destabilisation (≥ 64°).The greatest relative increase in displacements is observed in the second domain.
Destabilisation is more pronounced for warming closer to the melting point (from -2 to -0.5 °C) than for warming from -4 to -2 °C.This difference becomes greatest in the second domain.
(ix) The calculated critical slope angles and rock mass temperatures correspond well to the characteristics of documented rock slope failures in permafrost areas in the European Alps, which often showed large amounts of residual ice in their scars.(xi) The critical thresholds can be used to detect rock slopes which are susceptible to fail in the future and potentially endanger human life and mountain infrastructure.For this, it is a prerequisite to have data on the fracture network, lithology, geometry and the thermal field of the investigated rock slopes.In contrast, a detailed stability assessment of a single rock slope requires a number of further site-specific input data.
We are at the beginning of developing numerical models for the mechanical response of a rock mass to warming/thawing.
Therefore, there is a huge potential for improvements.For instance, more mechanical laboratory tests are necessary to improve process understanding and to extend the number of frozen and unfrozen material parameters for modelling.Future numerical models should be based on further sensitivity analyses, or include water pressure or the creep of ice at lower displacement rates within the advanced stage of rock slope destabilisation.Time-dependent cycling is not yet implemented in most of the numerical models, neither it is in the presented ones.Hence, they are not suited to forecast a precise time of rock slope failure.
Further, the spatial warming pattern for the Zugspitze model is simplified and based on a static temperature field derived from ERT, rock thermistors and published thermal models.For a more accurate characterisation of the warming pattern, future mechanical models should include a subsurface temperature field based on a combination of geophysics and a heat flow model for different time steps.

Note: a
The P-wave velocity of rock correlates highly to Mode I fracture toughness which correlates closely to Mode II fracture toughness (Chang et al., 2002).b P-wave velocities are parallel to cleavage/bedding.

Figure 1 :
Figure 1: (a) Left: Location of the study site (red box) in the Zugspitze summit area, ca.60 m below the peak (2963 m a.s.l.).The red arrow points to the sliding direction of the unstable rock mass.The hillshade is calculated from a digital elevation model with a cell 110

Figure 4 :
Figure 4: (a) Left: Laboratory-tested electrical resistivity of frozen and unfrozen Wetterstein limestone collected at the study site.210

Figure 5 :
Figure 5: Spatial warming/thawing pattern for the numerical model of the Zugspitze summit crest.For each warming/thawing step, the six rock layers (coloured areas) are assigned an individual temperature with a characteristic set of mechanical parameters.Warming or thawing proceeds from the outside to the inside of the crest.The dashed black line marks the current boundary between frozen and unfrozen rock, derived in Section 2.4.The presented thermal field is estimated and refers to model stage 11 in Table 6 (current state of permafrost distribution in the crest).

Figure 6 :
Figure 6: Calculated spatial distribution and magnitude of displacements for the Zugspitze summit ridge.(a)-(c) Factors of safety (FS) are given for Stages 1, 4 and 9 which refer to the loss of rock bridges and warming.(d)-(f) Stages 10, 11 and 15 correspond to thawing with FS below 1.The prominent shear zone is marked by a white-red dashed line.The permafrost boundary in (d) and (e) is highlighted by a white dashed line.

Figure 7 :
Figure 7: Numerically calculated block displacements (on the log scale) for a simplified rock slope geometry with exemplary inclinations of 60°, 62° and 64°, and five temperatures between -4 and -0.5 °C.

Figure 8
Figure 8 provides critical slope angles and rock temperatures beyond which instability is introduced.(a) Calculated maximum displacements against slope angle for temperatures between -4 and -0.5 °C.The curves are assigned to three distinct domains dependent on the slope gradient: Above 50°, and a second time above 62°, the curves show a sudden onset of instability.(b) Absolute and relative increase in displacements versus slope angle for a warming from -4 to -2 °C and from -2 to -0.5 °C.(c) Calculated factor of safety (FS) against slope angle for different sub-zero temperatures.(d) Absolute and relative decrease in FS versus slope angle for a warming from -4 to -2 °C and from -2 to -0.5 °C.

Figure 9 :
Figure 9: (a) A rule-of-thumb for critical slope gradients and fracture network orientations of anaclinal and cataclinal permafrost rock slopes.For cataclinal slopes, both the slope-face and the bedding have the same dip direction, but for anaclinal slopes, the dips of both are oriented oppositely (Cruden, 2003).The green area refers to stable rock slopes while the red area refers to unstable ones.For simplification, the fracture network is represented by the bedding dip.The coloured areas indicate the studied range of slope angles.(b) Five exemplary joint set constellations for cataclinal and anaclinal slopes with an inclination of 60°.The numbers of the rock slopes correspond to the numbers in (a).The joint spacing is ten-times bigger than in reality and is only applied for illustration.Example 5 refers to the Zugspitze summit crest.

(
iv) The Zugspitze model demonstrates a stability decrease towards a critical level as a result of (a) rock bridge destruction and (b) gradual warming of frozen rock and ice-filled joints from -4 to -0.5 °C.Surficial rock slope failure starts coincident to thawing of the outermost rock layer.Upon full thaw of the summit crest, expected within the next five decades, the model predicts an increase in displacements which potentially lead to final slope failure.(v)We developed a framework to generalise and upscale the Zugspitze model: Critical stability thresholds were calculated in the context of a sensitivity analysis of the Zugspitze model with simplified geometry and warming pattern.
(x) The critical thresholds can be applied to warming permafrost rock slopes with (a) ice-filled joints, (b) limestone equivalent to Wetterstein limestone, or probably most of the rock types relevant for permafrost rock slope failure in the Alps, (c) slope angles smaller than 70° and (d) various orientations of the fracture network consisting of three joint sets.