The Role of Temperature in the Stress–Strain Evolution of Alpine Rock-Slopes: Thermo-Mechanical Modelling of the Cimaganda Rockslide

In this work, the thermo-mechanical stress–strain history of an Alpine slope is analyzed, with particular focus on the historical Cimaganda large landslide (Sondrio Province, Italy), which mobilized an estimated volume of 7.5 mm3 of rock material. Accurate geomorphological and geomechanical characterization involving field surveys and laboratory testing was carried out, leading to the development of a conceptual model of the slope. A thermo-mechanical semi-coupled approach was developed, considering both glacial debuttressing and thermo-mechanical effects due to gradual exposure of the slope to atmospheric conditions and paleo-temperature redistribution resulting from the Last Glacial Maximum deglaciation. A 2D distinct-element numerical approach was adopted, supported by a 2D finite-element analysis to simulate heat diffusion over the Valley cross-section. Modelling results allow to simulate the general evolution of the Cimaganda rock-slope and to highlight the significance of thermal processes in preparing rock-slope instabilities. While the mechanical effect of ice thickness reduction alone brings about moderate rock mass damage, the introduction of temperature couplings results in a substantial increase of damage, representing a significant factor controlling the stress–strain evolution of the slope. Simulated displacement and the development of a deep region of shear strain localization at a depth roughly corresponding to that of the detected Cimaganda sliding surface, allow to highlight the significance of temperature influence in preparing the rock-slope to instability. The effects of surface temperature variations on the stress-strain evolution of an Alpine rock-slope are examined through numerical simulation. This study contributes to the understanding of Thermo-Mechanical processes driving damage in rock masses. Results showed how Thermo-Mechanical stresses induced by temperature variations at both long- and short-term timescales, can represent a significant preparatory factor to rock slope instability. The factors introduced in the analysis allowed to simulate the general evolution of the Cimaganda rock-slope, which is corroborated by field survey data. The effects of surface temperature variations on the stress-strain evolution of an Alpine rock-slope are examined through numerical simulation. This study contributes to the understanding of Thermo-Mechanical processes driving damage in rock masses. Results showed how Thermo-Mechanical stresses induced by temperature variations at both long- and short-term timescales, can represent a significant preparatory factor to rock slope instability. The factors introduced in the analysis allowed to simulate the general evolution of the Cimaganda rock-slope, which is corroborated by field survey data.


Introduction
Rock slope instabilities result from progressive rock mass damage, involving slip propagation along existing discontinuities, failure of intact rock bridges and the creation of new fractures. These processes, collectively referred to as progressive failure mechanisms (e.g., Terzaghi 1962; Erismann 1 3 and Abele 2001; Eberhardt et al. 2004), are brought about by the evolution of stress-strain state and strength variation within the rock mass. Factors leading to a long-term change in the stress state and resisting forces are known as 'preparatory factors' (Gunzburger et al. 2005) and mainly include climate-driven phenomena such as water pressure changes, long-term glacial evolution, surface temperature variations, chemical and mechanical weathering processes. Although the single increment of damage may be small in magnitude, their cumulative effect can lead to irreversible deformations. Moreover, factors with a multiple cyclic action over time can induce rock microstructure damage, causing progressive weakening of the rock mass mechanical properties and forcing slopes into a critical disequilibrium (Li 2016).
When limit conditions for slope stability are reached, the rock-slope collapse can occur in conjunction with an ultimate event, known as 'triggering factor'. Commonly reported cases invoke water pressure changes induced by rainfall events as the main cause of rock mass damage (Preisig et al. 2016;Loew et al. 2017;Grämiger 2020), especially as ultimate factor for rock-slope failures (Davies 2014;Liu 2021). Another significant trigger for slope instabilities can be represented by earthquake shaking, as discussed in Strom (2015).
Considering different natural slopes, preparatory and triggering factors do not act identically: geological and geomorphological background factors such as lithology, fracture geometry and topography (angle and height of slopes), define a different stress-strain behavior of slopes. Factors controlling the susceptibility to gravitational-instability events are generally referred to as 'predisposing factors'.
Focusing on the Alpine environment, glacial unloading is often suggested to be the predominant mechanism preparing slope instabilities, leading to the development of tensile rock mass failures and large-scale stress release (Holm et al. 2004;Cossart et al. 2008). Numerical analyses showed that failure processes could initiate through the development of tensile rock mass damage following glacial debuttressing (Eberhardt et al. 2004). However, a large lag time typically separates the end of deglaciation and the occurrence of large-scale slope instabilities, implying the unlikelihood for glacial unloading to act as a direct failure triggering factor (Prager et al. 2008;Ivy-Ochs et al. 2009;Zerathe et al. 2014). The complex stress-strain evolution affecting Alpine slopes after deglaciation and the occurrence of deep deformations may be explained using time-dependent constitutive models (Apuani et al. 2007;Agliardi et al. 2013), implying the occurrence of deformation under a constant mechanical stress state along the slope.
The Alpine environment, however, is also affected by significant variations of climate conditions on both short-and long-term periods, such as changes in subsurface temperature and groundwater circulation, which may in turn induce considerable stress field variations. The role of thermomechanical (TM) and thermo-hydro-mechanical (THM) couplings in driving rock-slope instabilities, despite having been widely acknowledged, is frequently overlooked due to the numerical computational costs, and only recently has been investigated in more detail (Wu et al. 2019).
The importance of including long-term temperature changes supporting deglaciation was highlighted by Baroni et al. (2014), who performed a numerical analysis with a continuum approach applied to the Adamello valley. TM effects can induce slope deformation and fracture propagation as a result of volumetric expansion and contraction caused by superficial temperature changes (Pasten 2013;Bakun-Mazor et al. 2013). Gischig et al. (2011) using a simplified elastic rock-slope model, showed that TM stresses can induce rock mass deformation and joint failures even at greater depths below the thermo-active-layer. Based on this observations Grämiger et al. (2018) performed a numerical analysis showing that TM stresses could represent a significant factor driving rock mass damage in paraglacial environments.
Moreover, several studies based on field data measurement and physical models at the rock-block scale, showed the capability of periodic temperature oscillations to induce plastic damage on discontinuous elements (Gunzburger et al. 2005, Bakun-Mazor et al. 2020. Recently Marmoni et al. (2020), studied the spatial distribution of daily near-surface temperature variations and the heat flux propagation through a jointed rock-mass in a rocky quarry wall, pointing out the importance of cyclic thermal loads in the propagation of rock mass failures along a natural slope.
Other processes that are relevant to a more detailed THM analysis and remain rather unexplored in rock-slope modelling literature, are represented by freeze-thaw cycles and permafrost degradation. If partially saturated conditions apply, freeze-thaw cycles can affect the strength of rock masses through both ice segregation-melting processes (i.e., repeated expansion and contraction within a joint inducing fracture propagation and failure of intact rock bridges) and thermal chemio-mechanical weathering. Freeze-thaw cycles have been recognized among the climate factors driving the occurrence of rockfall events in the Alpine region (Nigrelli et al., 2018;Draebing et al. 2019;). The triggering of many large rockslides in the periglacial environment is attributed to the melting of permafrost and interstitial ice (i.e., reducing rock-mass cohesion) due to an increase in the average atmospheric temperature; examples in the Alps are the Thurwieser (Sosio et al. 2004), the Morteratsch (Fisher et al. 2010) and the Cengalo rockslides (Walter et al. 2020).
In this work, the effects of surface temperature variations will be examined through a numerical simulation, to evaluate whether temperature may serve as a preparatory factor for large rock-slope instabilities. Considering (i) glacial debuttressing resulting from the Last Glacial Maximum (LGM) deglaciation, (ii) paleo-temperature redistribution and (iii) seasonal temperature fluctuations, a TM semi-coupled analysis is applied to investigate the stress-strain evolution of a slope, along which a massive rockslide event occurred (the Cimaganda rockslide-Italian Alps).
This study aims at giving a contribution in the understanding of TM processes driving damage in rock masses, highlighting the influence of climatic factors in the evolution of Alpine slopes. These aspects are fundamental to help in determining future scenarios of slope stress-strain evolution under climate changes.
The paper is structured as follows. In Sects. 2 and 3, the study area is described focusing on the main geological and morphological features connected to the Cimaganda instability event. In Sect. 4, the geomechanical conceptual model of the slope is built and the modelling procedure is discussed: a 2D distinct-element (DEM) approach was adopted, supported by a 2D finite-element (FEM) analysis to simulate transient heat diffusion over the valley crosssection. A thermo-mechanical semi-coupled analysis was performed, considering both glacial debuttressing and thermo-mechanical effects due to the gradual exposure of the slope to atmospheric conditions and paleo-temperature redistribution resulting from the Last Glacial Maximum deglaciation. Results are presented in Sect. 5 focusing on the simulated stress-strain fields in the rock-mass domain. Discussions and conclusions of the performed analyses are provided in Sects. 6 and 7, respectively.

Study Site: The San Giacomo Valley
The study area, known as the Valchiavenna region, is located in the Central Italian Alps between the village of Chiavenna and the Splügen Pass, connecting Italy to Switzerland (Fig. 1a, b). The modelled rock-slope is sited on the left flank of the San Giacomo Valley, in correspondence of the Cimaganda village. The valley, furrowed by the Liro Torrent, follows a N-S-striking tectonic lineament and represents the natural divide between the Lepontine and Rhaetian Alps.
The geological framework of the area is related to the Alpine Pennidic Nappe arrangement, characterized by the emplacement of two sub-horizontal gneissic bodies: 'Suretta' and 'Tambò' units. They consist of polycyclic Fig. 1 a Location of the study area; b distribution of the principal DSGSD bodies along San Giacomo and Bregaglia Valley and poly-metamorphic basement mainly composed by paragneiss and orthogneiss lithologies (Montrasio and Sciesa 1988). The tectonic contact between the two units gently dips to the E-NE and is characterized by the presence of a meta-sedimentary cover unit (Spluga-syncline) composed by quartzites, marbles and schists. The contact widely outcrops on the east slope of the valley at an elevation of about 1900-2200 m a.s.l., showing significant deformation and thickness variations.
In the Valchiavenna region, the main structural alignments have the following directions: WNW-ESE, NW-SE and NE-SW (Apuani et al. 2009). They reflect the orientation of main regional tectonic systems as the Insubric Line, Gruf Line, Forcola Line and Engadine Line.
The morphology of the valley is characterized by high sub-vertical rock cliffs where gravitational processes, favored by the orientation of discontinuity sets, define the actual slope dynamics . Post-glacial rockfall talus cone and landslides deposits are particularly frequent along the entire valley floor, with blocks even greater than 10 m 3 . The last relevant rockslide event, in May 2018, involved a rock volume of about 5.000 m 3 over the Gallivaggio sanctuary, a few kilometers south of the Cimaganda site (Carlà et al. 2019, Dei Cas et al. 2021. Deep-seated gravitational slope deformations (DSGSD), favored by the combination of high-relief landscapes and lithostructural framework, frequently affect both valley flanks with different states of activity (Fig. 1b). Morphological evidence of DSGSDs (trenches, counter-slopes and splitting of ridges) are particularly frequent along the tectonic contact between the Tambò and Suretta units (Aldighieri and Mazzoleni 2011).

Glacial History
During the Pleistocene age, in the last 2.6 million years, European Alpine valleys were affected by several glacial and interglacial periods (Bini, 1996).
Ice extension during the Last Glacial Maximum (LGM) in the Valchiavenna region has been reconstructed by Tantardini et al. (2013) in agreement with Bini et al. (2009). During the LGM the whole San Giacomo Valley was occupied by ice and only the highest peaks and ridges were ice-free. The ice thickness between the valley bottom and the trimline was 450 m in the upper valley, at the Montespluga Lake, and increased downstream up to 1550 m at San Giacomo Filippo Village due to the topographical and morphological constrains (Tantardini et al. 2013). In correspondence of the Cimaganda slope the ice level reached an elevation of 2150 m a.s.l, resulting in a glacial thickness of 1230 m. The San Giacomo glacier joined the Engadin-Bregaglia glacier in the Chiavenna area, making a single body which flowed southwards until it joined the Adda glacier.
Despite the maximum expansion of ice during the LGM has been completely mapped, no exposure ages of surfaces and glacial deposits are available for the study area. The LGM in the Alps is dated at ~ 28 to 18 kyr BP (Ivy-Ochs et al. 2015). Strong retreat of the LGM systems occurred about 19-18 kyr BP, when glaciers had progressively receded back within their mountain front (Ivy-Ochs et al. 2015;Wirsig et al. 2016). Successive and limited readvances occurred during the Late Pleistocene and Holocene period (Darnault et al. 2012;Ivy-Ochs et al. 2008).
As described by Tantardini et al. (2013), after LGM expansion, the San Giacomo glacier passed from a single and wide body to a smaller tongue forced to flow into the valley floor. As melting proceeded, the glacier thickness decreased and progressively receded back to the upper valley. In agreement with studies conducted in the Alps and Engadin region (Olhendorf 1998;Florineth 1998;Maisch 2003;Bini et al. 2009, Ivy-Ochs et al. 2015, for the study area at hand it is reasonable to assume glacial maximum conditions until 19-18 kyr BP, followed by a general melting phase. Once deglaciation in the principal San Giacomo Valley was completed and ice-free conditions were reached, valley flanks have been no longer subject to glacier advances.

Climate Regime
During late Pleistocene and Holocene age, deglaciation in Alpine regions was supported by a change in climate conditions with a general long-term air warming trend. Studying isotope records and pollen data, different authors provide paleo-temperature evolution in European areas (Ivy-Ochs et al. 2006;Vinther et al. 2009). According to Ivy Ochs et al.
(2008) the paleo-temperature evolution in the Alpine region can be outlined as follows: during LGM, mean temperatures were about 12 °C colder than today and began to rise following the deglaciation process; at the end of Pleistocene, mean air temperatures had risen to values of 3.5 °C colder than today and stabilized around the current values during Holocene age, varying only by ± 1 °C as shown by Davis (2003) based on pollen data. Current climate conditions exhibit typical features of an Alpine environment. Elevations in the San Giacomo Valley range from 333 m a.s.l. of the Chiavenna Village up to 3279 m a.s.l. of the Pizzo Tambò. Air temperature data relative to the last decades recorded by the three ARPA (Agenzia Regionale per la Protezione dell'Ambiente) stations located along the valley were analyzed to evaluate their annual evolution. Daily and monthly average temperature series are shown in Fig. 2. Annual fluctuations with negative average values during winter and positive values during summer were observed. Fitting the time series with a periodic sine function, a mean annual amplitude of 10.2 °C was identified. While mean annual temperatures depends on altitude, the amplitude of seasonal oscillation is not strongly affected by elevation (Fig. 2).

The Cimaganda Rock-Slope
The east flank of the San Giacomo Valley in correspondence of the Cimaganda Village, is dominated by significant morphological evidence of an ancient rockslide event, the so-called Cimaganda rockslide (Fig. 3).
The massive rockslide event mobilized an estimated volume of 7.5 Mm 3 of rock material and represents the failure evolution of a DSGSD bounded upstream by the structural terrace of 'Bondeno' and on the lateral sides by the 'Stua' and 'Avero' Valleys. Typical features of large-scale gravitational stress release, with trenches, counter-slopes and sub-vertical tensile fracturing have been identified along the slope (Morcioni et al. 2020, Morcioni et al. 2021. Rock masses are mainly composed by the gneissic bodies of the Tambò nappe unit: a spatial variation in textural and compositional features can be observed, with intercalated orthogneiss, paragneiss and amphibolite bodies of the Corbet zone. The tectonic contact between Tambò and Suretta units is located above the failure area and is highlighted by a bundle of meta-sedimentary rocks. Rock masses show clear evidence of glacial cycles, with striae and scoured rocks that indicate an ice flow direction corresponding to the main valley axis. Morphologies and glacial deposits (moraines and erratics) related to the last retreat phases of the San Giacomo Glacier, are quite widespread along the 'Bondeno' terrace, where flat morphologies allow their preservation (Fig. 3a).
The landslide detachment zone is wedge-shaped and is characterized by high subvertical rock cliffs (Fig. 3b). It is bounded by brittle geological elements consisting of (i) a system of vertical fractures with a N-S direction parallel to the slope orientation, (ii) a persistent discontinuity set merging towards NW defining its left flank and (iii) a WSW-ENE subvertical fault defining its right flank.
The landslide deposit is cone-shaped extending from the head to the valley bottom. Due to the high energy of the event, it climbed the opposite slope, damming the Liro torrent and partly modifying its path, leading to the deposition of thin and lacustrine deposits upstream of the accumulation zone. The deposit is composed by angular blocks of gneissic

3
rock heterogeneously distributed along the slope. The largest blocks are found in its lower part with mean volumes of several tens of m 3 and huge boulders with volume larger than 100 m 3 . Block sizes decrease towards the top of the deposit, where vegetation partially colonized the rock material.
The age of the event is difficult to estimate as no information or citation data have been found in the local historical archives. Lichenometric dating was performed by Mazzoccola (1993) on the rockslide deposit, indicating that the main massive sliding episode occurred before 900 A.D. Other more recent, secondary events probably involved the main detachment zone, leading to the current slope morphology. These evidences, also supported by numerical analyses (Mazzoccola 1993), suggested a two-stage failure evolution of the slope, favored by the presence of water pressure within the subvertical fractures representing the upstream release elements (Cancelli et al. 1994).
Along the Cimaganda slope, the presence of highly persistent and open fracture systems parallel to the valley axis (Fig. 1e, f), leads to occasional slope failure phenomena which, in most cases, involve limited volumes of rock material. However, they represent a serious threat for the Cimaganda Village and for the SS36 national road, which is the only connection between Chiavenna and the Splügen Pass. The right flank of the Cimaganda rockslide is particularly active, due to the presence of high rock cliffs. The latest significant failure event occurred in September 2012, mobilizing 20.000 m 3 of rock material, blocking the SS36 national road and isolating the upper valley for a few days (Morcioni et al. 2020).

Geomechanical Characterization
Following the ISRM suggested methods (ISRM 1981), geomechanical field surveys were performed in the Cimaganda area (Fig. 3a). Focusing on the fracturing conditions of the rock mass, at least three systems of discontinuity (K1, K2 and K3) were identified along the slope. Joint geometrical parameters as orientation (dip direction and dip), spacing (distance normal to joint tracks) and trace (length of joint segment) were collected and analyzed with a statistical approach, as follows.
• K1, representing pervasive foliation of gneissic rocks, gently dips towards NNE with a low inclination angle (dip direction/dip: 34° ± 10°/27° ± 3°; spacing: 1.5 ± 0.7 m; trace: 9 ± 3.5 m). It can be related to the regional tectonic contact between Tambò Bieniawski, 1989) and the evaluated GSI (Geological Strength Index; Hoek and Marinos, 2000) ranging from 60 to 70, reveal a good rock mass mechanical quality. At the crown of the Cimaganda rockslide, in the upper slope, lower GSI values (50 ± 5) were detected, due to a more intense fracturing and a higher weathering degree of the rock masses. As shown by kinematic analysis, gravitational instability along the slope is favored by the orientation of the joint systems relative to the slope direction. High persistent K2 shear surfaces may act as a basal sliding plane of significant volumes of rock material, which can be isolated upstream by the K3 tensile system. Moreover, in the upper part of the slope, the presence of disarticulated rock masses with highly persistent and opened joints, leads to rockfalls events, mainly as toppling or sliding evolution of blocks and rock-towers.
To investigate the mechanical behavior of both intact rock and discontinuities, uniaxial compressive tests on rock samples and direct-shear tests of joint surfaces were conducted according to the ASTM standards (D7012-14 and D5607-08). Gneissic samples collected along the rock scarp of the 2012 rockslide event, exhibit a mean uniaxial compressive strength (UCS) of 105 MPa with a Young modulus (E) of 35 GPa and a Poisson ratio (υ) of 0.3. Dominant discontinuity shear strength values obtained from direct shear tests showed a strong connection with the surface morphological features, with values of peak friction angle varying from 28° to 48°. The lowest values are referred to the smooth and weathered surfaces of the K2 set, whereas the largest ones to the rough joints of the K3 system. Multiple shear loading cycles were also applied to each specimen, to obtain a residual friction angle of 25° (K2) and 33° (K3). According to the ISRM Suggested Methods (ISRM 2015), normal and shear joint stiffness range, respectively, from 1.5 to 4.08 GPa/m and from 1.01 GPa/m to 1.45 GPa/m (Table 1), following the relation with morphological features of joint surfaces as described above for the strength parameters. 4 Numerical Analysis

Modelling Approach
Numerical analysis of the Cimaganda rock-slope was performed using the 2D DEM code UDEC (version 7.0-Itasca Consulting Group 2018), suited to analyze the behavior of discontinuous media. In the adopted DEM code, the jointed rock mass is modeled as an assemblage of distinct deformable blocks, that are in mechanical interaction through contacts representing discontinuities (Cundall and Hart 1992). Deformable blocks are defined by a continuum mesh of triangular finite-difference zones, where each zone behave according to an assigned linear or nonlinear stress-strain law; blocks can thus deform, translate, rotate and interact with each other. Rock joints are treated as interfaces between intact blocks, along which contact forces and stresses are examined. Deformability and strength properties of interfaces are represented by spring-slider elements located at contact points between a block corner and an adjacent block edge (Lorig and Cundall, 1989). If the prescribed tension or shear strength limit is exceeded, failure of the interface is assumed. Stresses can accumulate at slip fronts and additional failure or slip propagation can occur until the stress state drops below the strength limit. The DEM method is thus particularly suited to evaluate the stress-strain evolution of a discontinuous medium like the rock mass at hand.
In a distinct numerical approach, the process of crack propagation within a jointed rock mass can be simulated via the so-called Voronoi tessellation, which creates randomly sized polygonal blocks (Bishop 2007). Their sides represent randomly oriented elements which allow failure along new pathways. Assigning them mechanical properties of intact rock, the process of fracture propagation and the failure of rock bridges can be simulated: 'fracturing' occurs when the joint strength between Voronoi blocks is exceeded. The Voronoi tessellation procedure in UDEC is controlled by two factors: the "Average edge length" and the "Iteration number". The Voronoi algorithm begins by distributing points randomly within the region to divide. Then, an iteration procedure moves the points depending on the iteration number: the higher the number of iterations, the more uniform the spacing between points will be. Next, triangles are created between all points and the Voronoi polygons are defined by constructing perpendicular bisectors of all triangles that share a common side (Itasca 2018).
In this work, the effects of temperature changes in the stress-strain evolution of the Cimaganda rock-slope are analyzed with a thermo-mechanical semi-coupled approach. Starting from the Last Glacial Maximum (LGM) condition, the combined effects of (i) glacial debuttressing (ii) paleotemperature redistribution and (iii) seasonal temperature fluctuations are explored. Transient long-term temperature fields were calculated using the FEM code COMSOL Multiphysics (2018), that allows greater flexibility compared to UDEC in defining boundary conditions and material thermal properties. Moreover, its multi-physics structure will be useful for further planned developments including additional couplings, such as considering heat advection by groundwater or air circulation.
To analyze the effects of temperature changes in the stress-strain field along the slope, a linear thermo-elastic behavior of the rock mass was assumed. Variations in the stress state are thus induced by thermal strains resulting from isotropic volumetric expansion or contraction of rock blocks, governed by a constant thermal expansion coefficient. For each block, considering the constraint exerted by the presence of the surrounding blocks at all locations but the slope surface, thermally induced volumetric strains will cause stress state changes, which can even propagate at depths greater than the so-called thermally active layer (maximum depth at which a surface temperature change is propagated). Changes in the stress state within the rock mass domain will affect interfaces between intact blocks (discontinuities and Voronoi sides) for which an elasto-plastic behavior is assumed. As the stress state evolves, the shear stress at interfaces may reach failure conditions, thus inducing plastic straining.
In the above outlined framework, the effects of both glacial debuttressing and TM loading were analyzed, focusing on the propagation of fractures within the rock-mass domain which is bound to represent a preparatory factor to instability along the Cimaganda slope. The adopted modelling strategy, involving the use of two software was developed following the modelling strategy outlined in Fig. 4 and discussed in the next Sects. 4.2 and 4.3.

Geomechanical Conceptual Model
The cross-section represented in the numerical model (Fig. 5a), taken along the Cimaganda rockslide's direction of movement, extends from the west to the east slope of the San Giacomo Valley (Fig. 1b, 3a) and is based on the 2015 Digital Terrain Model of the Lombardy Region, modified by reproducing the pre-failure topography conditions, assuming a continuity with the lateral slopes. The geomechanical conceptual model includes two rock mass basic components: discontinuities and intact rock (Fig. 5a). Discontinuities are represented by a multiple joint network based on the three different systems identified with geomechanical surveys (Sect. 3.1). The UDEC builtin joint generator was used, to reproduce with a statistical approach similar joint patterns to those observed in the field. Relevant parameters for the joint generation procedure are represented by the dip angle (angle between joint track and x-axis), trace (length of joint segment) and spacing (distance normal to joint tracks). They were determined based on geomechanical field measurements performed on rock masses outcropping along the slope (Sect. 3.1). A statistical analysis of collected parameters was performed and representative values were obtained: for each of them, a mean value and a standard deviation representing random fluctuation around the mean were provided (Sect. 3.1). Considering scale effects, the geometrical values directly measured at outcropping rock masses were scaled up by a factor of 100 in the model. Also, to account for a confinement effect, joint spacing was increased at depths greater than 1000 m by a factor of 2 and at depths greater than 2500 m by a factor of 4 (Fig. 5a). Discontinuities were assigned an elasto-plastic Mohr-Coulomb constitutive law including slip-weakening of friction and cohesion (Table 1). The joint model aims at simulating displacement-weakening by loss of frictional, cohesive and tensile strength: once shear or tensile failure is reached, strength properties instantaneously drop from peak to shear residual values and null tensile strength.
Joint strength parameters were derived from both directshear laboratory tests (K2 and K3 discontinuity sets; Sect. 3.1) and from Barton-Bandis model (K1 set), developed using JRC and JCS data evaluated during geomechanical surveys (Sect. 3.1) and corrected in relation with the scale effect (Barton and Bandis 1982). So-called Voronoi tessellation was introduced to simulate intact rock, while allowing for the possibility of newly formed discontinuities (Wu 2019;Rwechungula 2021). Polygons were built with an average edge length of 100 m and iteration number of 5 (software default value). Considering the extension of the slope and the geometrical properties of discontinuities, these values allowed a proper geometrical representation with a good computational cost.
Rock strength parameters were calculated by applying the GSI approach (Hoek et al. 2002) considering results from geomechanical surveys and uniaxial compression tests (Sect. 3.1). This resulted in estimated values of cohesion (5 MPa), friction angle (53°) and tensile strength (1 MPa) that were assigned to Voronoi contacts assuming a Mohr-Coulomb constitutive law (Table 1).
Rock blocks, resulting from the intersection of discontinuities and Voronoi polygons, were set as deformable and discretized into a triangular finite-difference mesh. Blocks were assigned thermo-elastic properties considering a homogeneous and isotropic material (Table 1), where mechanical properties derived from site-specific laboratory tests (Sect. 3.1) while thermal expansion properties were taken from the literature (Cermak et al. 1982;Robertson et al. 1988;Grämiger et al. 2018) (Table 1).
Considering that TM analyses are computationally intensive, models require an as simplified as possible geometry. Hence, the modelled area was restricted to the east slope of the San Giacomo Valley, where the Cimaganda rockslide developed, and Voronoi tessellation was introduced only in the uppermost 500 m of this area (Fig. 5a). Roller boundaries were applied to the lateral sides of the model, whereas the bottom boundary was fixed. Domain borders were placed far enough from the area of interest to avoid boundary effects (Fig. 5a).

Initialization
The initial elastic stress field due to the gravity load and topography constraints under ice-free conditions was calculated by a mechanical initialization simulation until force-equilibrium was reached (Fig. 4). In this phase, the strength of discontinuities was set sufficiently high (i.e., cohesion and friction angle were increased by one order of magnitude) to prevent failure. Then, Mohr-Coulomb elasto-plastic properties were activated thus allowing joints to yield, and equilibrium stresses were recalculated (Fig. 4). After this two-step initialization procedure, the extent of mechanical damage and stress-strain distribution in the slope are deemed to represent the ice-free, pre-LGM mechanical conditions in the San Giacomo Valley. Next, the glacial load corresponding to the LGM ice level was introduced as an isotropic load boundary condition, at the slope surface (Fig. 5a). Assuming an ice density of 900 kg m −3 , ice loads were defined by the following simplified relation: where Z is the ice thickness, ICE the ice density and g the gravity acceleration. At LGM, the ice level in correspondence of the Cimaganda village reached an elevation of 2150 m a.s.l, resulting in a glacial thickness of 1230 m (Tantardini et al. 2013).
LGM ice boundary stresses were thus applied to the slope portion covered by ice, and another elasto-plastic calculation step was performed until mechanical equilibrium. This was set to represent the starting point of the analysis, hence resulting accumulated displacements were reset to zero.

Glacier Retreat
Once the initialization procedure was completed, the effect of glacial unloading was explored by reducing the ice thickness acting on the Cimaganda rock-slope following Eq. 1 (Figs. 4, 5a). To simulate the progressive and continuous deglaciation process, the numerical analysis was performed starting from LGM conditions and representing the glacial retreat with thirteen stages of ice surface lowering, until the ice-free condition was reached. Each stage corresponded to the removal of a 100-m thick ice layer, and, for each stage, mechanical equilibrium was calculated (Figs. 4, 5a).
Stepwise debuttressing is a valid approach to minimize unbalanced forces and reach numerical mechanical equilibrium without excessive computational cost. Moreover, for longterm mechanical studies on glacial time scales, modelling ice as a hydrostatic stress boundary condition leads to more realistic results compared to explicitly representing the ice body in the integration domain as an elastic solid material (Leith et al. 2014). In fact, the latter assumption would lead to significant overestimation of the damage associated with glacier retreat (Grämiger et al. 2017).

Long-Term Temperature Effects
The effects of long-term thermal changes related to the paleo-temperature evolution during LGM deglaciation, were numerically simulated and superposed to the ice-unloading effects.
Temperature distributions along the valley cross-section were simulated in COMSOL Multiphysics considering the rock mass as a continuous medium. The process of heat transport is governed by thermal conduction defined by thermal diffusivity (D), which depends on thermal conductivity (λ), specific heat capacity (C P ) and rock density (ρ) following the relation D = λ(ρC P ) −1 . Values of λ = 2.9 W m −1 K −1 and C P = 780 J Kg −1 K −1 were chosen for the problem at hand, corresponding to generalized values for gneissic lithologies found in the literature (Cermak et al. 1982;Robertson et al. 1988;Gramiger et al. 2018).
Adiabatic boundary conditions define the lateral sides of the thermal model, while a constant vertical geothermal heat flux of 65 mW m −2 (Swisstopo, 1982) is applied along the lower boundary (Fig. 5b). A prescribed temperature history was assigned to the slope surface, representing the main driver of temperature evolution within the rock mass governed by the ice level (i.e., ice-air interface) position. Specifically, below the glacier cover, rock surface temperatures were held constant at 0 °C, while above the ice level, they were set to depend on environmental conditions, as a function of altitude (y) and time (t). Rock surface temperatures were thus defined by: where 15.3 °C is the mean rock surface temperature at the reference elevation of y = 0 m; 0.005 °C m −1 represents the rate of temperature variation with altitude (y); ΔT paleo (t) is the factor that accounts for the paleo-temperature change relative to the present (Fig. 5d) (2) T � (t, y) = 0, for y ≤ glacier elevation T � (t,y) = 15.3 − (0.005 * y) + ΔT paleo (t) for y > glacier elevation Initial temperature conditions were obtained upon calculating steady-state thermal equilibrium under the LGM ice extension (Fig. 4). Transient temperature evolution was evaluated for each mechanical step of deglaciation defined during the glacial retreat analysis (i.e., every 100 m of ice lowering; Figs. 4, 5a). Considering a melting rate of 0.6 m year −1 , deglaciation is completed at 16 kyr BP (Fig. 5d). Once icefree conditions were reached, temperature distributions were evaluated every 1 kyr, since paleo-temperature stabilized around the current mean values (ΔT paleo (t) = 0) at 6 kyr BP (Figs. 4, 5d).
The transient temperature fields calculated by FEM were then imported into the DEM code, by assigning a temperature value to each grid-point as input data. UDEC was then run to a thermo-mechanical equilibrium, considering thermo-elastically induced strains, bringing about changes in the stress field (Fig. 4). This represents a thermo-mechanical semi-coupled approach, where temperature variations may result in stress changes, but mechanical changes do not result in temperature variations. This restriction is not believed to be of great significance for quasi-static practical rock engineering problems (Itasca 2011).
Adopting a standard thermo-elastic isotropic constitutive law, thermally induced stresses in the rock mass domain are thus given by: where Δσ ij is the change in stress, δ ij is Kronecker delta (δ ij = 1 for i = j and δ ij = 0 for i ≠ j), K is the bulk modulus, α is the linear thermal expansion coefficient, and ΔT is temperature change.
The above relation is applied by the DEM code assuming a constant temperature in each triangular zone of the mesh, which is obtained by an interpolation from the surrounding grid-points. Thermo-mechanical coupling is thus provided by the influence of temperature change on the volumetric strain of a zone, through the thermal expansion coefficient α. Thermally induced stresses are superposed to the previously calculated mechanical stress state.

Short-Term Temperature Effects
The annual seasonal temperature fluctuations were finally introduced in the simulation and their effect superposed to that of the long-term temperature changes and ice unloading. Due to the significant number of involved cycles, short-term temperature distributions were computed directly within UDEC with a transient thermal analysis, to simplify the coupling procedures (Fig. 4).
Blocks between discontinuities were assigned the same diffusive thermal properties as those defined for the FEM long-term analysis (Sect. 4.3.3). Below the glacier cover, temperatures were held constant at 0 °C, while, above the ice level, rock surface temperatures were set to depend on environmental conditions defined by altitude, time (paleotemperature changes relative to present) and seasonal fluctuations.
Including short-term temperature effects, Eq.
(2) that defines the rock surface temperatures, now becomes (Grämiger et al. 2018): The term Asin(2 ft) represents the sinusoidal temperature annual signal. The amplitude A was set to 10 °C, as discussed in Sect. 2.2 when analyzing the time series of temperature data recorded in the study area.
Starting from the long-term temperature distributions, a sinusoidal temperature history with 20 °C = 2A peak-to-peak amplitude and one-year period was imposed at the ice-free slope surface. Boundary conditions in the rest of the domain assumed zero heat flux at both bottom and lateral sides of the model, so that the short-term thermal effect can be consistently superposed to the long-term temperature distributions (Fig. 5c).
Modelling annual cycles throughout the entire late Pleistocene and Holocene period was deemed practically unfeasible, due to the extremely high computational cost with standard calculation facilities. Hence, a sensitivity analysis was performed to assess the impact of annual temperature oscillations on mechanical damage in the rock-mass domain, quantified through the magnitude of shear and normal displacement rates and the occurrence of plastic straining recorded along discontinuities with elasto-plastic properties. An initial significant damage signal, reflecting the effects of suddenly applying a sinusoidal temperature fluctuation to the slope surface, was observed. The magnitude of the damage signal progressively decreases to a stable level and becomes negligible after about thirty peak-to-peak cycles, when the system reached a quasi-static equilibrium. These results are in accordance with Gischig et al. (2011) who observed a consistent decrease in the damage signal after five years using a simplified rock-slope model with elastic properties.
Results from the sensitivity analysis allowed to significantly reduce the number of annual cycles to be modeled. For each step of modelling (i.e., for every change in ice level and long-term temperature distribution), a number of annual cycles Nc = 50 was assumed to be sufficient to accommodate the effect of the new sinusoidal thermal boundary conditions (this corresponds to parameter Nc reported in the modelling flowchart of Fig. 4). The conductive heat equation was solved by the UDEC code via a Finite Differences implicit scheme with a timestep of 1 day. A thermo-mechanical equilibrium

Results
The output of the above outlined numerical analysis is evaluated in terms of stress-strain field redistribution in the rock mass domain and the development of plastic shear strain and tensile failures along discontinuities. The impact of temperature in the stress-strain field and in the development and propagation of fractures within the rock mass, is assessed by comparing results from the purely mechanical glacial unloading analysis with those from the long-term and short-term thermo-mechanical analyses. To evaluate localized stress changes, some monitoring points were defined at different depths along the slope at which the evolution of simulated normal and shear stress and displacement rates were recorded. Considering that block interface elements (discontinuities and Voronoi polygons) behave as an elastoplastic medium, the stress path at the selected points was plotted and compared to the corresponding strength envelope, analyzing the role of thermal loads in the occurrence and propagation of plastic failures.

Glacial Retreat Mechanical Analysis
As a result of topographical constrains, during LGM the glacial overburden was maximum at the bottom of the valley and progressively reduced moving upwards along the slope. The ice thickness acting on the valley floor was 1230 m, resulting in a maximum normal stress of about 10.5 MPa. Debuttressing leads to a general stress reduction within the rock-mass domain and a progressive stress redistribution along the slope. Volumetric stresses (σ kk ) decrease with a maximum variation of 20 MPa at the bottom of the valley, following the glacial loading pattern (Fig. 6a-a'). Differential stresses (σ max -σ min ) acting on the valley floor at a depth of 100 m increased from 5 to 15 MPa, as a result of reduction in the minimum principal stress component (σ min ). In contrast, the valley flanks experienced only a slight Fig. 6 Volumetric stresses (a-a') and differential stresses (b-b') simulated along the slope at two different steps of the analysis: LGM and Ice-free. In plots (a-a') and (b-b') extensive stresses are positive, while compressive ones are negative. The effects of ice unloading are shown by the spatiotemporal damage distribution (c) and the normal displacements (d) of discontinuities, the shear strain increment (e) and the total displacements (f) of blocks. The stress path recorded at point 'A' is presented in (g) differential stress variation, consistently with the smaller ice loading action (Fig. 6b-b').
Stress redistribution along the slope was thus not homogeneous, leading to a variable damage pattern. Simulated new tensile and shear failures mainly occurred as a propagation of discontinuities that had already undergone failure at LGM. Further damage mainly occurred at the valley bottom and at the midportion of the slope where the magnitude of stress changes was largest (Fig. 6c). It can be observed that debuttressing also induced opening of the subvertical discontinuities, along which normal displacements reach values over 5 × 10 -3 m as shown in Fig. 5d.
Expectedly, total displacements show a general uplift, as a result of elastic rebound due to unloading, reaching the maximum values at the valley floor (Fig. 6e, f). Overall, considering a purely mechanical analysis, the effect of ice removal does not seem to create significant rock mass damage, as it mainly results in a post-glacial elastic volume increase. Figure 6g shows the stress path simulated at point A, located along a mid-slope subvertical joint at a depth of 100 m. During initialization at ice-free equilibration, the stress state lies (in absolute value) below the peak failure strength envelope (stress-path point 1). Additional LGM ice loading caused an increase in the normal stresses, thus the stress state moved further away from the envelope (point 2). Then, progressive ice unloading induces a normal stress decrease and a shear stress increase, and the unloading path is roughly coincident with the loading one, implying reversible (elastic) loading conditions. However, when the ice level reached an elevation of 1350 m a.s.l and the slope surface above point A became ice-free, the stress state moved closer to the failure envelope (point 4). As deglaciation proceeded, normal stress in the considered joint further decreased, so that by the end of deglaciation, the joint ended up more critically stressed compared to the initial ice-free conditions, although failure was not attained.

Long-Term Thermo-Mechanical Analysis
The transient thermal analysis performed with the FEM approach allowed to evaluate the temperature redistribution along the slope during deglaciation, considering both ice level lowering and climate conditions changing. The initial temperature field for the transient thermal analysis was defined under LGM conditions with a stationary solution (Fig. 7a). Thermal gradients are shown to be highly influenced by the extent of the glacier cover, holding the ground surface temperature at 0° C. During deglaciation, the slope is progressively exposed to atmospheric conditions, which depend on the altitude and on the paleo-temperature factor (Sect. 4.3.3). Considering that the 0 °C isotherm for the mean ground temperature during LGM is located at 500 m a.s.l. (Eq. 2), during deglaciation ice-free rock surface portions are gradually exposed to negative temperatures, resulting in a progressive cooling of the slope. Ice-free valley conditions are reached at 16 kyr BP, when the 0 °C isotherm has risen to 1400 m a.s.l. (Eq. 2). Figure 7b shows the differences between LGM and ice-free temperature distributions. It can be observed that the slopes portion located between the LGM ice level and the ice-free 0 °C isotherm experienced a net cooling process (blue scale colors), while, in the rest of the slope, temperature globally increased (red scale colors). Simulated temperature variations are restricted to the slope subsurface within a maximum depth of thermal diffusion of about 300 m (Fig. 7b). The deglaciation-induced temperature difference distribution suggests that the slope surface experienced significant thermal shocks due to sharp temperature gradients, with special reference to (i) the top and (ii) the lower-mid portion of the scarp. The latter area, shown in Sect. 5.1 to be affected by most of the debuttressing-induced mechanical damage, is bound to especially suffer from thermally induced degradation.
Once deglaciation is completed, temperature distribution is no longer affected by the ice level lowering and depends only by the ΔT paleo factor (Eq. 2). This results in a progressive heating of the slope until 6 kyr BP, when temperature changes relative to present become negligible (ΔT paleo = 0; Fig. 4d). Considering that the ΔT paleo factor is independent of altitude, the heating process in this phase can be considered homogeneous along the entire slope, causing the maximum depth of temperature diffusion to increase to values of about 1000 m (Fig. 7c). At the end of the analysis, the 0 °C ground isotherm has risen to 3000 m a.s.l. (Eq. 2) and the slope has experienced an overall heating process (red scale colors in Fig. 7c).
Temperature fields were then imported in the DEM model as input data for the long-term TM analysis, where the effects of paleo-temperature changes, associated with the deglaciation process, were explored. Compared to a purely mechanical model, the TM one accounts for additional stresses due to the temperature variations along the slope (Eq. 3).
In Fig. 8, the variation in volumetric stresses (σ kk ) and differential stresses (σ max -σ min ) between the LGM, ice-free and Holocene thermal stabilization conditions are shown. During deglaciation (from 18 up to 16 kyr BP), the effects of glacial unloading and slope surface cooling, which causes a rock-mass volume contraction, leads to a general stress reduction within the rock-mass domain. Focusing on the mid-slope area, volumetric stresses decrease (Fig. 8a), while differential stresses increase (Fig. 8b) as a result of a reduction in the minimum principal stress component (σ min ).
The subsequent temperature and stress redistribution between the ice-free and the Holocene thermal stabilization conditions (from 16 up to 6 Kyr BP), is characterized by an extensive, overall heating of the slope. The rise in temperature induced a general volumetric expansion of the rock blocks, causing increments in the volumetric stresses along the entire slope, with maximum values of Δσ kk = 10 MPa (Fig. 8a'). Differential stresses progressively increased, starting from the valley floor and extending upslope, identifying a region of high differential stresses at the toe of the slope (Fig. 8b'). Monitoring the out-of-plane stresses, a noticeable variation in their magnitude could be observed, resulting in some cases in assuming the maximum principal stress value, especially at shallow depths where the slope morphology highly affects the stress distribution.
During deglaciation, the progressive cooling and contraction of rock blocks promoted both tensile and shear failures along the subsurface, as well as the opening of sub-vertical discontinuities (Fig. 8c, d). The subsequent heating phase promoted damage at greater depths where the confining stress level is higher, mainly as shear failures of subvertical discontinuities, with significant additional slip propagation Fig. 8 Simulation results show the increment of volumetric stresses (a-a') and differential stresses (b-b') along the slope at two different steps of the TM analysis: 'ice-free' and 'thermal stabilization'. In plots (a-a') and (b-b') extensive stresses are positive, while compressive ones are negative. Increments are referred to the LGM stress distribution presented in Fig. 6. The effects of ice unloading and long-term temperature change at 'ice-free' and 'thermal stabilization' conditions, are shown by plotting the spatiotemporal damage distribution (c-c') and the normal displacements (d-d') of discontinuities, the shear strain increment (e-e') and the total displacements (f-f') of blocks. The stress path recorded at point 'A' is presented in (g) as shown in Fig. 8c'. The volumetric expansion of blocks also induced a closure of discontinuities, especially at depth, where the stress increment is larger (Fig. 8d').
Focusing on the mid-slope area, where the Cimaganda rockslide occurred, a shear strain localization region is observed to begin to develop at the end of the deglaciation (Fig. 8e) and to progressively expand as a result of the heating process (Fig. 8e'). This region corresponds to the area of most significant differential and volumetric stress increase. Intensive shear straining within the elastic rock blocks gradually caused plastic yield along discontinuities and boundaries of the Voronoi polygons: a significant increment in the extent of joint failure and slip propagation along sub-vertical discontinuities was observed (Fig. 8c'-e'). As a result, simulated total displacements concentrate at this portion of the slope (Fig. 8f').
Examining the stress-path recorded at point A (Fig. 8g), introducing the long-term temperature effects a different stress state is attained compared with the purely mechanical analysis (Fig. 6g). Until the joint remains below the ice level and slope temperature is kept constant by the glacier cover (points 2-4), the stress path is coincident with the mechanical one. However, when point A becomes ice-free, temperature changes penetrate to depths of 300 m (Fig. 7b) and the rock mass is subjected to a cooling process, resulting in ΔT = -1 °C. When deglaciation is completed (point 5), stress conditions move closer to the failure envelope. The subsequent phase (points 5-6), characterized by a global heating trend (ΔT ~ 7 °C), leads to an increase in shear stress: the stress path thus reaches the failure envelope causing joint failure, and the instant attainment of residual strength conditions, according to the adopted joint constitutive law (point 6).

Short-Term Thermo-Mechanical Analysis
Introducing short-term thermal effects, as the slope surface is exposed to atmospheric conditions it is suddenly affected by temperature cycling on a seasonal basis, that are superposed to the long-term temperature changes. Considering that the diffusion depth of seasonal temperature changes is limited to the first tens of meters, the yearly temperature oscillations affect only the shallowest blocks along the slope surface, causing their cyclic isotropic expansion and contraction. However, thermally induced volumetric strains lead to stress variations which can propagate at depths well below the thermally active layer, as is shown analyzing the stress path recorded at point A and plotted in Fig. 9a.
Until the monitored joint remains covered by ice, no changes in the stress path were observed compared to the previous analyses. When ice level reached an elevation of 1450 m a.s.l., the slope surface above point A started being subjected to seasonal temperature fluctuations, and TM induced cyclic stresses are recorded in the stress path (point 3). However, at this stage, the discontinuity and the surrounding blocks are not already completely ice free, hence stress variations exhibit a little magnitude, of about 50 kPa (Fig. 9b).
TM cyclic stresses become highly significant when ice level reaches an elevation of 1350 m a.s.l., and the monitored joint is completely ice-free. In fact, heating phases cause an expansion of blocks, inducing in turn an increment in normal stresses with an amplitude of ~ 0.1 MPa at point A; on the other hand, cooling phases cause a contraction of blocks inducing a decrease of normal stress and a shear stress increase (hence an overall decrease of shear strength, as the failure envelope is approached) with an amplitude of ~ 0.1 MPa. After 21 peak-to-peak annual thermal cycles, failure is attained (Fig. 9b). It is important to note that temperature at point A is not affected by seasonal temperature fluctuations, as shown in Fig. 9c. In fact, point A is located at a depth greater than the thermo-active layer, so temperature values are only affected by the long-term changes. Variations in the stress state are thus induced only by contractive/ expansive thermal volumetric strains of rock blocks located above point A, closer to the slope surface.
It can be observed that the development of TM stresses within the elastic blocks brings about rupture propagation along the elasto-plastic interfaces. Slip dislocation along discontinuities follow the thermally induced cyclic pattern, involving both directions of shear: although a small amount of permanent dislocation remains after each full cycle, most of the slip occurring during the summer period is recovered in the subsequent winter one (Fig. 9d). These features apply as long as the discontinuity remains within the elastic domain; when plastic failure conditions are reached, irreversible slip occurs, and will accumulate during the entire process (Fig. 9d). The transition from peak to residual resistance results also in a stress redistribution in the rock mass domain, as highlighted by the increment of differential stresses (Δ(σ max -σ min )) shown in Fig. 9f. Thus, the stress redistribution induced by failure may promote critical loading of neighboring discontinuities, supporting the process of fracture propagation.
Compared with the purely mechanical and long-term TM analyses, when additional short-term temperature fluctuations are considered, the extent of unstable rock mass increases, extending further toward the slope top. New damage occurs mainly along the slope surface, represented by (i) tensile failures of discontinuous elements that represent the sides of blocks directly affected by the thermally induced volumetric changes, and (ii) plastic shear failures along subvertical joints, driving the process of TM stress propagation ( Fig. 10a-a').
Focusing on the mid-slope area where the Cimaganda rockslide occurred, the shear strain localization region that 1 3 begins to develop at the end of the deglaciation and progressively propagates as a result of paleo-temperature heating, is now more extended (Fig. 10b). Irreversible displacements take place, as a result of the combination of elastic glacial rebound and long-and short-term thermal effects, especially along the shallower blocks that are directly affected by shortterm temperature fluctuations (Fig. 10b).

Discussion
Results presented above are in accordance with other TM studies that identified the importance of including longterm temperature changes supporting Alpine deglaciation (Baroni et al. 2014;Grämiger et al. 2018) and the role of periodical thermal loads in the propagation of rock mass failures along a natural slope (Gunzburger et al. 2003, Pasten 2013Bakun-Mazor et al. 2013. In this work, the combined effect of long-and shortterm thermal loading was explored overlapping the warming trend resulting from the Last Glacial Maximum (LGM) conditions to the Holocene age and the seasonal temperature oscillation. By comparing ice-unloading analysis with the TM ones, the role of temperature in the stress-strain evolution of rock-slopes was highlighted, demonstrating that induced damage and displacements are strongly increased in the presence of TM stresses.
Considering a purely mechanical model, the effect of ice removal does not create significant rock mass damage, and it mainly results in a post-glacial elastic rebound.
Although subvertical joints open and become more critically stressed, the occurrence of plastic failure is limited to the propagation of yield along pre-LGM failed elements. At the end of deglaciation, no critical conditions for slope stability were observed along the slope. The subsequent introduction of both long-and short-term TM stresses leads to a substantial increase in the occurrence of plastic yield and slip propagation along discontinuities. During deglaciation, as soon as the slope is exposed to atmospheric conditions, seasonal temperature oscillations induce significant new tensile failures along the slope subsurface. Even if the effects of temperature fluctuations are progressively accommodated by the system, normal and shear displacements continue to accumulate along critically stressed and failed discontinuities. This results in accumulated inelastic deformation along joints, able to induce surface block instability. Considering the modelling approach, these TM effects can be interpreted as a large-scale fatigue process, involving incremental slip along critically stressed discontinuities following seasonal cyclic loading. Unlike the short-term thermal loading, the long-term one (i.e., paleo-temperature changes from LGM conditions up to present) is monotonic, it occurs over a much longer time scale and includes larger values of temperature variations. Hence, the diffusion of long-term signals reaches larger depths enhancing stress-strain field variations, as emphasized by the occurrence of deep plastic yield during the heating phase that follows deglaciation. While the short-term temperature variations mainly induce instability at the slope surface, long-term thermal effects can represent an important factor in the development of large-slope deformations.
The application of the analysis on a slope along which a large instability event occurred, allowed to calibrate the modelling approach and to validate results. As corroborated by geomechanical field surveys (see Sect. 3), the Cimaganda rock-slope evolution is mainly controlled by the presence of opened joints and highly persistent sub-vertical discontinuities, which promoted instability events along the entire slope (Fig. 11). This field evidence supports the numerical simulation results, confirming the occurrence of new plastic shear and tensile failures, leading to superficial rock-mass instabilities, especially if seasonal temperature fluctuations are considered. In accordance with field observations, damage is particularly evident at the top portion of the slope in correspondence of the Bondeno structural terrace, where highly disarticulated rock-masses outcrop (Fig. 11b).
During the deglaciation process, sub-vertical discontinuities become critically stressed, and the subsequent heating phase promoted plastic yield as well as the propagation of slip. At the end of the analysis the predisposition to instability along these elements is quite evident. Shallow rock sliding such as the documented 2012 Cimaganda event, could be easily initiated if destabilizing driving forces are introduced in the model (Fig. 11c, d).
Focusing on the central part of the slope where the historical rockslide event occurred, in the numerical model it can be observed that a shear straining region begins to localize at the end of deglaciation and progressively propagates as a result of the heating process. The presence of such strain localization area does not bring about critical Fig. 10 The combined effects of ice unloading, long-and short-term thermal loading at 'ice-free' and 'thermal stabilization' conditions, are shown by plotting the spatiotemporal damage distribution (a-a') and the normal displacements of discontinuities (b-b'), the shear strain increment (c-c') and the total displacements of blocks (d-d') deep-seated instability conditions, but it clearly defines a preferential region for the development of a failure zone, at a depth roughly corresponding to that of the detected Cimaganda sliding surface. Moreover, this area corresponds to the portion of the slope where the magnitude of simulated displacements (of the order of 10 -1 m) takes the maximum values. The resulting displacements and the development of a deep shear strain localization zone are in agreement with the evidence of the historical Cimaganda collapse. Even if a critical condition for slope stability is not reached in the simulation, a preferential zone where damages concentrate was observed. This zone is bound to develop further, once other destabilizing factors are considered and introduced in the model.
To bring about critical instability conditions, it would be necessary to induce damage and propagation of fracturing along the sub-planar elements parallel to the slope surface, which are not critically stressed by thermal loading due to their unfavorable orientation that prevents accumulation of tangential stresses.
As pointed out by Grämiger (2018), if a strength reduction analysis is applied to the model, sufficient rock-mass damage to destabilize the slope and initiate landsliding will be generated. This process, however, would force the rockslope to instability and is thus not really suited to back-analyze a failure event. For that purpose, triggering factors such as groundwater circulation and surface hydrology dynamics should be introduced in future development of the analysis. Implementing a coupled thermo-hydro-mechanical (THM) analysis, in addition to exploring the role of water pressure, will also allow to examine the advective component of heat transport.
Considering numerical studies conducted in the Alpine environment, strains and total displacement resulting from the presented TM analysis are comparable with the effects induced by other widely acknowledge preparatory factors, pointing out the significance of considering atmospheric temperature in the evolution of rock-slopes. As an example, numerical studies conducted in the upper San Giacomo Valley on the Vamlera DSGSD with similar geological frameworks and geometries, have shown a similar trend and magnitude of displacements considering the effects of deglaciation using a time-dependent constitutive law (Apuani et al. 2007). In oncoming developments of the analysis, the introduction of a constitutive model able to account for (i) a time-dependent visco-elasto-plastic behavior and (ii) the evolution and degradation of mechanical properties with thermal loads, will allow to enhance the understanding of time-dependent deformation processes and the progressive achievement of ultimate strength conditions that can affect natural slopes.
Finally, as discussed in the Sect. 5.2, a noticeable variation in the magnitude of the out-of-plane stresses could be observed with the proposed plane strain analysis, especially at shallow depths where the slope morphology highly affects the stress distribution. This suggests that a more realistic TM analysis will be achieved if a fully 3D model is developed.

Conclusions
In this work, a thermo-mechanical semi-coupled numerical approach was presented and applied to investigate the stress-strain evolution of an Alpine rock-slope where a massive rockslide event occurred. By coupling a thermal FEM with a mechanical DEM analysis, TM effects associated with ice redistribution, paleo-temperature evolution and seasonal temperature fluctuations were analyzed. Results showed a clear relation between TM stresses and the occurrence and evolution of damage within the rock mass.
The factors introduced in the analysis allowed to simulate the general evolution of the Cimaganda rock-slope, which is corroborated by field survey data. The predisposition to instability is highlighted by the occurrence of tensile failure of subvertical elements, along which irreversible slip can accumulate and propagate with depth. The resulting displacements and the development of a deep shear-strain localization zone are in agreement with the evidence of the historical collapse event, clearly defining a preferential region for the development of a failure zone. To bring about critical conditions for a deep-seated and large-scale instability event, additional destabilizing factors should be introduced in the analysis. To this end, groundwater circulation and hydrological effects are bound to play the most important role.
Considering the modelling assumptions, the numerical results suggest that: • The only action of unloading due to glacier melting (socalled debuttressing) appears insufficient to generate deep critical damage along discontinuities. However, it could represent an important factor in the development of critically stressed and open joints which may boost the effects of other coupled factors.
• Paleotemperature changes between LGM and Holocene thermal stabilization induce important variations in the stress state as a result of thermo-elastic expansion of rock blocks, causing failures and deep deformations. • Short-term temperature fluctuations induce significant new damage along the slope subsurface, causing irreversible instability of shallower blocks. TM stresses resulting from the periodical expansion and contraction of blocks can propagate at depths that are not directly affected by the thermal diffusion process, causing cyclic shear and normal displacements until the limiting strength of discontinuities is reached. When failure occurs, irreversible slip can accumulate, preparing the rock mass to global instability.
Hence, TM stresses induced by temperature variations along the slope surface at both long-and short-term timescales, can represent a significant preparatory factor to rockslope instability in the Alpine environment. Temperature factors should not be neglected when the stress-strain evolution of rock slopes is analyzed, as well as during the assessment of rock fall susceptibility. This is bound to become an increasingly important aspect, in light of forthcoming climate change-related temperature exacerbation.