Rock Glacier Dynamics by a Thermo-Elastic-Viscoplastic Constitutive Relationship

Abstract

As a result of mountain permafrost creep, rock glaciers are common features in high-altitude periglacial areas. From a practical point of view, beyond their localization and inventorying, both the monitoring and prediction of their evolution due to climate changes are crucial. One of the effects of climate change is the thickening of the basal shear zone (the portion of the rock glacier where most deformations are localized), eventually leading to the development of unexpected and unprecedented (in terms of location, magnitude, frequency, and timing) instability phenomena. These phenomena bear consequences for the understanding of landscape evolution, natural hazards, and the safe and sustainable operation of high-mountain infrastructures. Most of the studies about active rock glaciers are focused on the analysis of monitoring data, while just a few studies are focused on modeling their behavior to understand their possible further evolution. The active rock glacier response is characterized by a viscous (rate-dependent) behavior, influenced by seasonal temperature oscillations, and characterized by a seasonal transition from slow to fast. In this work, a new thermo-mechanical model based on the delayed plasticity theory and calibrated on experimental results is proposed. The model is employed to evaluate the influence of geometry and forcing (air temperature) on a real rock glacier (Murtèl-Corvatsch rock glacier) creep behavior.

1. Introduction

The main processes related to climate change in high-mountain environments involve the role of permafrost degradation in promoting slope instability (amongst others, [1,2]). Rock glaciers are globally recognized as an expression of creeping permafrost in the periglacial zones of high-latitude and high-altitude environments (e.g., [3,4]). According to [5], rock glaciers are quaternary forms of perennially frozen debris supersaturated with interstitial ice and ice lenses or even with bodies of massive ice moving downslope.

The long-term monitoring of rock glacier movement in the European Alps indicates that, in recent decades, the increase in ground temperatures [6,7,8] induced rock glacier velocity increases, on both seasonal [9,10,11] and decadal [12,13,14,15,16,17] timescales.

For instance, as shown in [18], in the European Alps, annual terrestrial geodetic surveys suggest that the average rock glacier surface velocity has increments of about +52% each year during the monitoring period (2010–2014), with most of the rock glaciers reaching surface velocity maxima in 2014 [6].

Temperatures inside a rock glacier are governed by heat conduction/convection processes (e.g., [19]), whose boundary conditions are the mean annual air temperature (MAAT) and the permafrost base temperature. Temperature increments may induce local ice melting (e.g., [20]), associated with both variations in water pressure and flow rate in the rock glacier body and a reduction in ice intergranular bonds strength, potentially inducing local instability phenomena [21]. Hence, to correctly analyze the mechanical response of rock glaciers, hydro/thermo/mechanical coupled problems have to be solved. To further increase the complexity of the problem, a severe spatial heterogeneity in ice content and debris grain size, ranging from silty sand to large boulders [22], is usually observed in rock glacier bodies. Smaller volumetric ice content and debris grain size are usually observed in the shear horizons [23,24,25], the parts of the rock glaciers, usually characterized by a thickness of a few meters (e.g., [25]), where most shear strains (up to 90%) and strain rates are concentrated [26]. In addition, local ice melting within different parts of the rock glacier body may induce differential rates of downslope creep and consequently geomorphology variations. However, only a few studies have measured internal movements in rock glaciers [26,27,28], which are essential for investigating rock glacier instability.

The first step for the solution of hydro/thermo/mechanical coupled problems is the definition of constitutive relationships capable of reproducing the material behavior. To this aim, experimental laboratory tests on samples directly collected from rock glaciers (e.g., [21,28,29,30]) were performed. The test results highlighted that the strength of frozen materials essentially depends on their composition, in terms of mixtures of ice and rock elements, and on temperature. The experimental results are, however, quite rare due to difficulties in accessing these areas, and in collecting samples, especially at depth. It is worth mentioning that, with regard to the dimension of the rock elements of a rock glacier body, the response of small-scale laboratory samples is not necessarily representative of large-scale behavior. In the future, to overcome laboratory test limitations, “numerical experimental tests” (for instance, by using the discrete element method, as is shown in [31,32]) can be employed, in the framework of multiscale approaches, for better understanding the global system response.

The constitutive models currently adopted to reproduce the mechanical behavior of this material are an adaptation of Glen’s [33] or Nye’s [34] laws, originally introduced for describing glacier movements. These laws are empirical power law functions, depending on the state of stress, temperature, and ice content. One limitation of these empirical laws is that they cannot satisfactorily reproduce the response of the shear horizons and the rest of the rock glacier body with a single set of constitutive parameters [20]. In other words, they cannot exhaustively describe the heterogeneity of rock glaciers characterized by debris and debris ice-riched layers (e.g., [25,28]). Moreover, as will be shown in Section 2.1, these laws do not always satisfactorily reproduce the previously cited laboratory test results. Several authors (e.g., [12,28,35,36]) demonstrated a strict link between the rising MAAT and accelerating creep velocity of alpine rock glaciers. This acceleration is usually attributed to a variation in the thermal regime of the ground, induced by climate change [14,22,37,38,39].

The final goals of the research carried out are: (i) a better understanding of the currently observed dominant controls on geomorphological changes and (ii) the introduction of reliable and simple tools to model the response of rock glaciers, capable of reproducing both short- and long-term adaptation of these systems to changing environmental factors.

As a first step toward this direction, in this paper, the authors intend to introduce a simplified approach for modeling the thermo-mechanical coupled response of rock glaciers (Section 2). This is based on: (i) the introduction of a new constitutive law, conceived according to the delayed plasticity theory proposed in [40], for describing the mechanical behavior of ice-rock mass mixtures (Section 2.1), (ii) the solution of the heat diffusion equation (Section 2.2) and (iii) the schematization of the rock glacier into a one-dimensional problem (Section 2.3).

To exemplify the practical application of this simplified approach, the case of Corvatsch rock glacier, located in the Swiss Alps is presented and the in-situ measurements compared with the numerical simulation results in Section 3. Finally, in Section 4 the results of a sensitivity analysis allowed us to put in evidence the main geometrical factors influencing rock glacier response.

2. A Simplified Approach to Model Rock Glacier Response

The numerical modeling of the mechanical response of rock glaciers is usually based on the employment of a continuum approach. According to this approach, to solve the problem, a system of partial differential equations describing (i) the balance of momentum, (ii) the compatibility condition and (iii) the material constitutive law (Section 2.1) has to be solved. In addition, since temperature plays an important role, the heat diffusion equation (Section 2.2) also has to be solved and the coupling between mechanical and thermal processes has to be taken into account.

In general, numerical (e.g., finite element or finite difference) models can account for the correct geometry of the system and the heterogeneity/spatial variability of the rock glacier internal structure (e.g., [21,39,41,42]). However, in most of the cases, a reliable description of geometry and heterogeneity/variability is missing. Hereafter, a simplified one-dimensional geometry (Section 2.3) is considered, and the material is assumed to be homogeneous (the material mechanical properties/density/thermal diffusivity are assumed constant along depth).

2.1. Constitutive Relationship

The constitutive relationship proposed by the authors is conceived according to the delayed plasticity theory proposed in [40].

According to this theory, strain rate tensor is decomposed into an instantaneous/reversible and a delayed/irreversible part (${\dot{\epsilon }}_{ij}^{irr}$). With respect to standard plasticity theory, the consistency condition is abolished (implying that the effective stress state can be either inside or outside the yield locus f) and the flow rule is expressed as it follows:

$${\dot{\epsilon }}_{ij}^{irr}=\mathsf{\Phi }\frac{\partial g}{\partial {\sigma }_{ij}}$$

(1)

where g is the plastic potential, ${\sigma }_{ij}$ is the stress tensor and Φ is the viscous nucleus. Φ can be defined as the distance (measured through a suitably defined mapping rule) between the state of stress and the yield locus position. In most of the cases, Φ is a non-negative and increasing function of f. Under this assumption, f may be interpreted as a scalar measure of the probability of occurrence with time of irreversible strain accumulation.

The elastic law is assumed to be linear and isotropic. The yield function (f), plastic potential (g) and the viscous nucleus (Φ) are defined by starting from the experimental results reported in [29]. According to these results:

• The material mechanical response is practically unaffected by the confining pressure;
• The material shear strength (both peak and residual) depends on temperature;
• During creep tests, the logarithm of the deviatoric strain rate linearly increases with the applied stress deviator.

For the sake of simplicity, the yield function is defined according to the Tresca criterion, where the cohesion (c′) is assumed to depend on temperature (T). Since rock glaciers accumulate large irreversible strains with time, only the residual conditions are accounted for. The yield function is defined in a non-dimensional form:

$$f=\frac{{\tau }_{eq}}{c^{\prime }}-1$$

(2)

being

$${\tau }_{eq}=\frac{{\mathsf{\sigma }}_I-{\sigma }_{III}}{2}$$

(3)

where ${\mathsf{\sigma }}_I$ and ${\sigma }_{III}$ are the maximum and the minimum principal stress, respectively. In practice, Equation (3) represents the radius of the Mohr circle associated with the state of stress. The evolution law of c’ with T, defined and calibrated based on the experimental results of [29] (points of Figure 1) and performed on ice–silty sand mixtures, is expressed as:

$$c^{\prime }=c_{1 } \left(1-e^{\frac{T}{c_2} }\right)$$

(4)

where c₁ and c₂ are, in general, expected to be dependent on the ice content. The line reported in Figure 1 is obtained by imposing c₁ = 500 kPa and c₂ = 1.2 °C, and it is valid for ice content values ranging between 75 and 90%. For T = 0°, when the ice intergranular bonds are completely thawed, cohesion is expected to be nil, since, owing to the large strains developing in the shear zone, the unbonded granular material is expected to be at a critical state (implying the cohesion to be nil). The unbonded material shear strength is expected to be associated with a residual friction angle (an additional parameter that is difficult estimate) and effective stresses (depending on pore water pressure, this can also be difficult estimate). The introduction of this additional contribution surely enriches the model, but it makes the calibration significantly more complex. To provide safe side estimations, the contribution of friction to the material shear strength is disregarded.

The plastic potential is defined as:

$$g={\tau }_{eq}-c^{\ast }$$

(5)

where c* is a dummy variable.

The employment of a non-dimensional yield function and of a dimensional plastic potential is particularly convenient, since, under this assumption, the viscous nucleus parameter values are expected not to depend on the applied stresses.

To obtain a linear relationship between deviatoric stress and the logarithm of deviatoric strain rate, an exponential viscous nucleus is employed:

$$\mathsf{\Phi }\left(f\right)=\eta e^{\alpha f}$$

(6)

where η and α are constitutive parameters.

To calibrate η and α, the experimental results reported in [29] (triangles of Figure 2) obtained for T values ranging between −1.92 and −1.98 °C were considered. The dashed red line is obtained by integrating the constitutive relationship under triaxial creep stress paths and by imposing T = −1.95 °C and η = 3.5 × 10⁻⁷ 1/s and α = 5.2 (

Table 1. Summary of imposed and calculated parameters.

Parameters	Value
C1	500 kPa
C2	1.2 °C
η	3.51–7 1/s
α	5.2

).

To validate the model, the experimental results reported in [29] (circles of Figure 2) for T ranging between −0.8 and −0.83 °C (not employed in the model calibration) are compared with the constitutive model prediction (dashed-dotted line of Figure 2). Even in this case, the agreement is very satisfactory.

For the sake of completeness, in Figure 2, the line obtained by employing the Glen’s law with the parameter values suggested by [12] for T = −0.81 °C is also reported (solid line). This relationship cannot satisfactorily reproduce the experimental results.

2.2. Heat Transmission

By following what has been suggested by other authors [12,20], the heat transmission is assumed to be only governed by conduction. The differential equation describing the variation of temperature in time (t) and space is the well-known heat equation:

$$\kappa \mathsf{\Delta }T =\frac{\partial T}{\partial t}$$

(7)

where Δ is the Laplace operator and κ is thermal conductivity. As suggested by [12,13,14,15,16,17,18,19,20], the κ value is assumed to be coincident with the one of ice. For the sake of simplicity, the heat generation is associated with the accumulation of irreversible strains. As shown in the following (Section 3), these simplifying assumptions introduced for heat conduction do not compromise the capability of the model of reproducing the variation in time of temperature.

It is worth mentioning that the problem hereafter studied is only “partially coupled”: temperature influences the material response (Equation (2)), but the mechanical problem does not influence the thermal one. This implies that the thermal problem can be solved independently, and its results can be employed as input data for the mechanical one.

2.3. One-Dimensional Problem Schematization

Analogously to what was proposed in [12,20], the rock glacier is assimilated to an infinite slope h-thick, θ-inclined (Figure 3) and resting on a rigid bedrock. The rock glacier body is subdivided into three subdomains (Figure 3): (i) an upper part where ice is present only in cold seasons (subdomain 1 of Figure 3), (ii) the shear horizon (subdomain 3 of Figure 3), where most of the strains develop and the particle size distribution is analogous to the one employed by [29] for the experimental tests, and (iii) the portion between the two (subdomain 2 of Figure 3), where strains are negligible owing to the interaction between large rock boulders.

For the sake of simplicity, for the material in subdomains 1 and 2, a rigid behavior is assumed, whereas the constitutive law introduced in Section 2.1 is employed to reproduce the mechanical behavior of material belonging to subdomain 3.

The whole domain is under simple shear conditions. Therefore, normal (σ_n) and tangential stresses (τ) can be assessed by simply imposing the balance of momentum (di Prisco and Pisanò, 2011). Since rock glacier accumulates large strains in the shear horizon, the authors assume the evolution of stress σ_t (Figure 3), associated with the static redundancy of the problem [43,44], to be ended and, according to the employed (nil) dilatancy angle value (Equation (5)), to be equal to σ_n [44]. It is worth mentioning that, owing to the employed yield function and plastic potential, the response is not influenced by the value of the intermediate (out-of-plane) stress.

As far as the thermal problem is concerned, only subdomains 2 and 3 are considered. To solve the thermal problem, (i) the initial temperature along the whole stratum and (ii) the variation in time of temperature at the top of subdomain 2 and at the base of the shear horizon have to be known.

In the following section (Section 3), this simplified geometrical scheme will be employed to reproduce the response of Murtèl-Corvatsch rock glacier. This case was chosen since (i) after the first investigations of this area started in of the 1970s [5], it is now one of the best investigated permafrost areas in the Alps and part of the PERMOS network [6,23,45,46] and (ii) the stratigraphies of two boreholes located within a distance of 30 m shows significant small-scale heterogeneities in the rock glacier [28,41].

3. Case Study: Murtél-Corvatsch Rock Glacier

The Murtèl-Corvatsch rock glacier (amongst others: [22,23,28,46]) is situated in the upper Engadine (Figure 4), facing Piz Corvatsch (3300 m a.s.l.), in the southeastern part of Switzerland (UTM, 563 131, 5 142 001; zone 32T). It extends from sea level from 2620 to 2850 m a.s.l. It is 400 m long and 200 m wide (total area of about 0.08 km²), characterized by a mean slope of about 10° and marked furrows and ridges on its surface. Within the area, the Murtél-Corvatsch rock glacier is one of the dominant periglacial features, but further rock glaciers and talus slopes are present to the west (Figure 4). In 1987, a 60 m deep borehole was drilled and equipped with thermistors that have been manually logged in 1 month intervals until 1992, and since then data have been stored automatically by a logger collecting temperature data in 6 h intervals [46]. The site is characterized by a coarse blocky (the diameters of the boulders are in the range of decimeters up to several meters) surface layer of approximately 3–3.5 m in thickness above a massive ice core down to 28–30 m and a frozen blocky layer underneath, reaching from 28 to 50 m, probably adjacent to the bedrock [21] (Figure 5). The ice core has an average temperature of −2 °C at 10 m depth and −1.4 °C at 25 m depth (Figure 6). The active layer has an average thickness of about 2 m. The involved lithologies mainly consist of granodiorite and granite (Schneider et al. 2012), and studies on back-weathering at this site have shown back-weathering rates of 2 mma⁻¹ (Müller et al. 2014). For Corvatsch rock glacier, likely one of the most investigated quaternary forms (e.g., [23,46]), there exists a long-term and continuous data acquisition process (i.e., ground and depth temperatures, superficial displacements). Attempts to determine the age of this rock glacier [47,48] obtained an age of 5000 to 6000 as a minimum value [49]). The rock glacier is characterized by a rather slow creep velocity (10–17 cm/year) and is considered a thick and ice-rich landform with an average volumetric ice content of 60% [20,21,47].

Annual precipitation at the site is about 900 mm (982 mm St Moritz 1951–1980; Segl-Maria 978 mm 1979–2021; 900.1 mm Piz Corvatsch 1979–2021) [6]. Typical maximum snow cover thickness is between 1 and 3 m. The mean ground annual surface air temperature is 1.8 °C for the observation period of 2011 to 2019, while the mean annual air temperature since 1979 at Piz Corvatsch weather station (46.41° N 9.82° E, 3305 m a.s.l.) is −5.4 °C.

4. Model Application and Results

The Corvatsch rock glacier case study is discussed by employing the simplified approach presented in Section 2.

The slope inclination is assumed to be equal to 10° and the shear horizon is assumed to be characterized by a thickness equal to 2 m, ranging from a depth from 28 to 30 m (Figure 3).

For the solution of the thermal problem, the in-situ measurements [6] of variation in time of temperature from 2016 to 2020 are employed. In particular, the distribution of temperature on 1 January 2016 is employed as initial conditions (Figure 7a), whereas variation in time of temperature at the base of the shear horizon and variation at the top of subdomain 2 of Figure 3 are employed as boundary conditions (Figure 7b). The numerical results in terms of variation in time of the temperature along depth will be compared with the in-situ measurements in Figure 8.

The numerical results in terms of temperature time history are in satisfactory agreement with the experimental measures (Figure 8), confirming the suitability of the simplifying assumptions (Section 2.2). This also suggests that during the considered time-period (50 months) significant variations in the thermal properties of the system (e.g., variations in the ice content or development of melted water flow paths in the rock glacier body) did not take place. Moreover, the results also highlight that those seasonal variations in temperature are observed only in the superficial layers: for the considered period, for depths larger than 14 m, the variation is lower than 1 °C.

By following what was suggested by [12], for calculating stresses acting in the rock glacier body, the ice unit weight is employed. The numerical results are reported as variation in the time of displacement profile along depth in the shear horizon (Figure 9a) and variation in the time of surface velocity (Figure 9b, obtained by GPS measurements). The displacement distributions (Figure 9a) along depth in the shear zone is practically linear, implying that strains are almost constant. The ground surface velocity is almost constant in the considered time period. Its average value (28 cm/year) is in reasonably good agreement with the in-situ measurements (10–17 cm/year).

A multiscale approach is essential to benefit from physical data in laboratory tests and physical models, as validation and parameter calibration opportunities.

4.1. Sensitivity Analysis

In the previous section, it was shown that the model proposed can reproduce the mechanical response of rock glaciers. In this section, a sensitivity analysis was carried out to explore the influence on the surface velocity of different input geometries and parameters. The difference in shear zone thickness, slope inclination and the depth of shear zone values were considered. For all the cases considered, the thermal boundary and initial conditions are the ones of Figure 7.

The results are summarized in Figure 10. Surface velocity linearly increases with the shear zone thickness, testifying that for all the cases considered almost constant strain distributions in the shear zone were obtained.

As expected, surface velocity also significantly increases with the inclination.

Considering an increment on shear zone inclination up to three times (5°, 10° and 15°), the results show an increment up to 300% (Figure 10). Absolute mean velocities are strongly affected by variations in geometry, due to changing stress conditions (Equation (3)). However, over the timescales considered in this study (seasonal to multi-annual), the geometry of a rock glacier is not expected to change substantially.

Larger velocity values are also observed for larger shear zone depths. This implies that, for the time-period considered and for the imposed thermal boundary conditions, the influence of the state of stress (stresses increase linearly with depth) is significantly larger with respect to the one of temperature (larger for more superficial layers).

However, this conclusion is not general, since a considerable increase in temperature due to global warming would lead to faster flow, as also reflected in the observational data sets presented in [12,20].

4.2. Model Assumptions and Limitations

For the sake of clarity, in this paragraph the limitations of the proposed model, clarifying the assumptions made, are discussed.

First, we considered the rock glacier dimensions, i.e., a thickness of 30 m circa for the body and 2 m for the shear zone, which are significantly smaller with respect to the planar dimension (400 m). As is also done for landslides in these cases, the assumption of schematizing the geometry as one-dimensional is very common. In principle, the approach proposed here by the authors can also be extended to two- or three-dimensional problems. Nevertheless, in these cases an accurate definition of both ground surface and shear zone geometry is crucial. In all the cases in which an accurate description of geometry is missing, implying that two- and three-dimensional solutions are not necessarily reliable, one-dimensional solutions have to be preferred, at least for these preliminary evaluations.

Furthermore, despite the heterogeneity shown by quaternary forms such as rock glaciers, a simplification on the homogeneity of the shear band is necessary in all the cases in which the information regarding the geometry and the composition/mechanical behavior of the material in the shear zone is very limited. Even if the rock glacier considered here is one of the most studied, an accurate description of the shear horizon, allowing one to properly characterize the heterogeneity, is missing. For this reason, the authors assumed the shear horizon to be homogeneous.

Concerning the thermo-mechanical model, the main assumptions are as follows:

• The thermal problem is characterized by a complex time-varying boundary condition (on the top). In very simple, but unrealistic, cases (e.g., linearly increasing temperature or harmonically varying temperature), a closed form solution can be derived. Unfortunately, for the imposed boundary conditions (Figure 7b), an analytical solution for the thermal problem cannot be derived. Moreover, due to both the thermo-mechanical coupling and to the non-linearity of the constitutive relationship, which are also for the mechanical problem, a closed-form solution cannot be derived.
• The authors assumed that the presence of ice bonds provides a cohesion to the material. When the ice is completely melted, cohesion is expected to be nil. The melted material (a rock fragments–water mixture) is expected to behave as a standard soil at a critical state (i.e., without cohesion). However, since the lack of measurements of pore water pressure, a reliable evaluation of the effective stresses and, therefore, of material (frictional) shear resistance is not possible. For the sake of simplicity, the authors disregarded this contribution.
• From the experimental data of [29] (Figure 2), we derive that during creep tests, the logarithm of the deviatoric strain rate seems to linearly increase with the applied stress deviator. To reproduce this experimental trend, the authors introduced an exponential viscous nucleus.
• When irreversible processes take place (in this case, the accumulation of irreversible strains) energy is dissipated (the dissipation can be measured as the product of stresses and irreversible strain rates). In principle, this dissipated energy is expected to induce a temperature increase (heat generation). For the sake of simplicity, in the proposed model, the authors disregarded this aspect. This assumption, however, does not compromise the model capability of reproducing the experimental results (Figure 8).

5. Conclusions

In this paper, a thermo-elastic–viscoplastic 1D model is proposed. The numerical model couples heat conduction to a temperature-dependent elastic–viscoplastic constitutive relationship, phenomenologically calibrated on a series of experimental laboratory test results. The proposed constitutive law is suitable for reproducing the response of the material belonging to the shear horizon. In principle, the same model can also be employed to reproduce the response of the other parts of the rock glacier body, but, up to now, laboratory test results are not available. In the future, to overcome the limitation of laboratory apparatus, “numerical experimental tests” (for instance, by using the discrete element method) may be employed for understanding the local (at the intergranular ice bonds scale) thermo-mechanical processes and introducing more complex constitutive laws, conceived by following suitable upscaling (micro-macro) approaches.

By comparing the model prediction with in-situ measurements of the tested rock glacier (Murtèl-Corvatsch rock glacier complex), the main findings of this study can be summarized as follows:

• Numerical results in terms of temperature time history are in satisfactory agreement with the experimental measures, confirming the suitability of the proposed model;
• In the considered time period (about 4 years), significant variations in the thermal properties of the system do not take place;
• The proposed model reproduces the correct order of magnitude of mean surface velocities.

Additional sensitivity analysis shows a major correlation between the geometry of the rock glacier body (i.e., characteristics of shear zone) and velocity.

In spite of these limitations, this work represents the first attempt to model in detail, through an elastic–viscoplastic constitutive relationship, the rock glacier response, by considering the intrinsic spatial variability of a rock glacier body and to use laboratory experimental test results as constraints. The approach is very promising but has to be adjusted in some details. The ice bonds provide cohesion to the material, i.e., when the material is melted the cohesion is nil and the material behaves as a soil at a critical state. Therefore, to properly reproduce the behavior, a reliable evaluation of effective stresses is necessary, implying that in situ local measurements of pore water pressure are needed. In the future, the model may be improved by introducing a more sophisticated constitutive relationship considering both intergranular friction and local thermo-mechanical processes taking place on the ice bond scale.

To validate respective models, however, more knowledge on the thermo-mechanical processes in the material has to be gained by performing suitable laboratory tests.