The ski resort town of St Moritz, Switzerland, is partially constructed on a large creeping landslide, which has been causing damage to buildings and infrastructure. At the town centre, the landslide is constrained by a rock outcrop, creating a compression zone in the sliding mass. After decades of gradual slowing down,s in the beginning of 1990s the landslide started to accelerate, in spite of the fact that the average yearly precipitation and the pore water pressure on the sliding surface remained fairly constant. The paper shows that a constrained creeping landslide experiences progressive failure caused by the propagation of a zone of intense shearing along the slip surface resulting in significant earth pressure increase and visco-plastic yielding of soil in the compression zone. This basic physical mechanism, relying on extensive laboratory and field tests and long-term displacement monitoring, explains the paradox of the St Moritz landslide acceleration. Although the model predicts that the landslide could eventually slow down, its displacements may become excessive for some buildings, requiring an early warning system and further stabilization of the historic Leaning Tower. In general, by predicting the onset of yielding, the model can provide an important timeframe for stabilization of constrained landslides.
Understanding of failure mechanisms of creeping landslides is of critical importance for assessment and mitigation of their hazard (Terzaghi 1950; Chandler & Pook 1971; Glastonbury & Fell 2008; Schulz et al. 2009). In general, when a creeping landslide is constrained by an obstacle, either natural (a rock outcrop) or artificial (a retaining wall), it slows down, creating an impression of being stabilized (Bernander & Olofsson 1981; Wiberg et al. 1990; Puzrin & Sterba 2006). Sometimes, however, the stabilization phase is succeeded by acceleration (e.g. Puzrin & Schmid 2007), which in the case of the St Moritz-Brattas slide (figure 1) resulted in the last 19 years in more than threefold increase in its displacement rates (figure 2). The onset of landslide acceleration can be observed around September 1991, when the time–displacement curve changes its curvature.
The St Moritz slide (figure 1a) has a length of 1500 m, width of 600–800 m and an average inclination of about 20° (Schluechter 1988). The landslide consists of two parts (figure 1b), with some geological evidence of a rock outcrop at the boundary between them. The upper 800–900 m are a rockfall of boulders reaching 2–3 m in diameter; the lower part of 600–700 m is the 17–23 m thick sliding layer built of a silty soil matrix with boulder inclusions. Inclinometer and Increx measurements taken in 2006–2010 in the lower part show that the shear strain is largely localized within the slip surface (figure 3a), while the compression along the slide leads to the heave in the perpendicular direction, which is fairly uniformly distributed within the sliding layer thickness (figure 3b). At the bottom of the lower part the landslide is constrained by a rock outcrop, creating a compression zone in the sliding mass.
In the lower 600–700 m of the landslide, the average displacement rates (measured in 2006–2010 close to the latitudinal centre of the landslide) increase uphill from the rock outcrop (figure 4) and can reach up to 0.5 m yr−1. Further landslide acceleration may, therefore, have potentially devastating consequences, in particular for sensitive structures, such as the twelfth century Leaning Tower of St Moritz (figure 1a), located close to the compression zone of the landslide. The Tower was stabilized in 1986 (before the landslide acceleration) and its inclination was partially corrected by lifting it with hydraulic jacks and placing on a new foundation on three Teflon-bearing pads (Sterba et al. 2000). Additional correction took place in 2005 (Alonso et al. 2010).
The ‘usual suspect’ for the creeping landslide acceleration is an increase in precipitation (e.g. Iverson et al. 2000). Indeed, the landslide displacements exhibit sensitivity to precipitation rate changes (e.g. in 2001–2006, when deviation of the precipitation rates from the average value is observed in figure 2). Around the onset of acceleration, however, no precipitation increase could be observed (figure 2), which was also confirmed by fairly constant ground-water table (figure 5), measured using a piezometer pipe installed in a borehole about 37 m above the rock outcrop at a depth of 14 m. No rising in the ground-water table has been observed around the onset of acceleration in 1991. On the contrary, the ground-water table around 1991 seems to drop below the average level (with an average depth of 9.92 m), which appears to be the same both in the period of 1984–1996 (around the onset of acceleration) and during the entire measurement period of 1978–2010.
A more likely reason for the landslide acceleration is the passive failure of soil at the bottom of the landslide, in the zone characterized by high compression strain rates (figure 6), measured by Tschudi & Angst (1999) within the constructed area in the lower 200 m of the landslide in 1988–1998, close to the latitudinal centre of the landslide. Location of the boundary of the compression zone can be estimated at l≈21 m, below which the highest compression strain rates were observed.
It remains, however, to be explained how the earth pressure could (i) reach the high level of passive pressure when the landslide was slowing down, and (ii) keep increasing even after the failure. In fact, a small but steady pressure increase of about 0.8–1.0 kPa yr−1 has been measured in the compression zone in 2008–2010 using a novel inclinodeformometer (Schwager et al. 2009). While the rate dependency of the residual shear strength on the sliding surface (2% per log cycle of the shear strain rate measured in ring shear tests, figure 7) could, in principle, explain some pressure increase during stabilization (Puzrin & Sterba 2006; Alonso et al. 2010), it is not sufficient for the passive pressure to develop. Moreover, when landslide accelerates, this rate dependency should lead to the pressure decrease, contrary to the observations.
Inability to explain the St Moritz landslide evolution using conventional approaches became the main motivation behind the field and laboratory testing programme carried out in 2006–2010 at the ETH Zurich. This experimental study was combined with basic analytical modelling of the data from many decades of the St Moritz-Brattas landslide monitoring in an effort to provide a mechanism explaining both the stabilization and acceleration phases of the constrained creeping landslide evolution.
2. The mechanism
Such a mechanism is indeed available and is based on the phenomenon of the progressive propagation of a zone of intense shearing, which causes earth pressure increase in the compression zone of the landslide, followed by visco-plastic yielding of soil in this zone (figure 8). Strain localization and propagation of shear bands (few millimetre thick zones of intense shearing, where shear strength decreases to residual value) take place in soils exhibiting strain softening (Palmer & Rice 1973; Vardoulakis et al. 1981; Saada et al. 1994; Puzrin & Germanovich 2005; Saurer & Puzrin 2011). In a strain-softening material (figure 9), shear strength τ* on the slip surface drops from the peak τp to residual value τr when the shear strain γ in the shear band reaches the critical value of γr.
Progressive shear band propagation has been shown to cause delayed failure of slopes cut in overconsolidated clays (Skempton 1964; Skempton & LaRochelle 1965; Bjerrum 1967). Catastrophic shear band propagation has helped to explain the mechanisms of gigantic tsunamigenic landslides (Nisbet & Piper 1998; Puzrin et al. 2004, 2010). Fewer attempts have been made, however, to investigate this mechanism in relation to creeping landslides, in particular to those constrained by natural or artificial obstacles (Bernander & Olofsson 1981; Wiberg et al. 1990).
Furthermore, the zone of intense shearing, propagating along the slip surface in the lower part of the St Moritz landslide, is not necessarily a thin classical shear band. The slip surface is confined there to a weaker, about 20 cm thick silty clay layer with gravel and organic inclusions located at the depth 17–23 m according to inclinometer measurements and core drilling. The thickness of this layer may not be sufficient for a single continuous thin shear band to propagate in it (Palmer & Rice 1973), as indirectly confirmed by the lack of slickenside surfaces in the core samples. Therefore, the shear deformation is likely to involve the entire clayey layer, either being smeared over its thickness d≈20 cm or localized in a number of parallel discontinuous thin shear bands, as can be also concluded from a relatively broad inclinometer trace of the band (figure 3a). In any case, this shear zone retains the important property of the shear band: shear strength τ* in it reduces to residual value τr when the shear strain there reaches the critical value of γr (figure 9).
A simple conceptual model of the St Moritz landslide is shown in figure 8. The sliding layer of thickness h, with the depth of phreatic surface hw and thickness of the sliding surface d, is subjected to the uniform effective pressure p′t acting at its uphill boundary, and has zero displacement δ on the downhill boundary. Where the landslide displacement has exceeded the value of δr=γrd, the shear resistance τ* along the sliding surface has dropped down to the residual shear strength τr indicating formation of the shear zone. Beyond the tip of the shear zone there is a process zone, where the shear resistance grows to the peak value τp and then decreases to the gravitational shear stress τg in the compression zone. Rate dependency of the residual shear strength has been neglected.
In contrast, the soil in the sliding layer has elasto-visco-plastic properties (figure 10). Pre-yielding compression in soil is visco-elastic, characterized by the first Kelvin element in the series, with Young's modulus E and viscosity η. When the pressure reaches the yield stress p′y, the second Kelvin element with Young's modulus ξE and viscosity ξη is mobilized, causing a reduction in the equivalent stiffness and viscosity during yielding. In general, different yielding factors ξ could be applied to elastic and viscous parameters, though adopting them equal simplifies the model considerably.
Creeping and yielding of the soil in the compression zone cause the tip of the shear zone, where shear strain has reached γr and the shear resistance fell to its residual value τr (figure 9), to move down along the slip surface. This leads to the progressive propagation of the shear zone and the drop of the shear strength along the slip surface, which has to be compensated by an increase of the effective earth pressure in the compression zone. After this pressure reaches the yield stress p′y (figure 10), the elastic and viscous resistances decrease significantly, causing increase in the landslide displacement rates. They do not, however, vanish entirely (as would happen at the classical perfectly plastic failure), allowing for continuing earth pressure increase in the compression zone.
4. Analytical model
The above assumptions lead to the following simplified analytical model allowing for the proposed mechanism to be quantified and validated against the monitoring and test data. Combination of equilibrium and geometric equations with constitutive relationships within the compression zone produces a differential equation, which can be resolved in a closed form.
(a) Equilibrium equation
The equilibrium equation for the sliding layer in figure 8 is in general given by 4.1and 4.2where p is the total uniform normal earth pressure in the direction of sliding; τ* the shear resistance along the slip surface; τg the gravitational shear stress; γt the total unit weight of soil in the sliding layer; h the thickness of the sliding layer and α the inclination of the slope and the slip surface.
In the compression zone 0≤x≤l(t), however, we assume τ*=τg, so that the earth pressure does not vary along the zone p(x,t)=p(t) and its evolution in time is defined by the equilibrium of the sliding layer above the compression zone: 4.3and 4.4where pt is the uniform total pressure acting at the uphill boundary of the sliding layer; τr the residual shear strength along the slip surface (figure 9); L the length of the sliding layer; l the length of the compression zone; γw the total unit weight of water; hw the depth of the phreatic surface; φ′r the residual angle of internal friction on the slip surface.
The above assumption neglects the ‘process zone’, where the shear resistance along the slip surface drops from its peak value τp to the residual value τr (figure 8). In reality, depending on the value of γr, this zone could appear to be of a non-negligible length within the scale of the problem. Considering, however, the general indeterminacy with respect to the shear resistance along the slip surface in the compression zone, and the fact that it is going to be partially smaller than τg and partially larger, it seems reasonable to average the resistance as τ*=τg. The validity of this assumption will be assessed later as a part of the overall ability of the model to simulate landslide displacements and pressures.
Water flow in an infinite slope produces identical pore water pressure profiles along the slope, so that p=p′+u and pt=p′t+u, and the equilibrium equation (4.3) can be rewritten in effective stresses: 4.5It follows, that the earth pressure in the compression zone increases in time proportionally to the decrease in the compression zone length.
(b) Constitutive relationship
The elasto-visco-plastic behaviour (figure 10) of the soil in the sliding layer is different for pre-yielding and yielding phases. Pre-yielding compression in soil is visco-elastic, characterized by the first Kelvin element in the series is given by 4.6where ε is the total strain in the model; E is Young's modulus of the soil; η the viscosity of the soil and p′y the yield stress.
When the pressure reaches the yield stress p′y, the second Kelvin element is mobilized. Since the two Kelvin elements work in series, the stresses in both elements are equal to the total stress p′: 4.7while the strains are different and have to be summed up to define the total strain and its rate: 4.8where ξ is the yielding factor. Multiplying the first equation (4.7) by this yielding factor ξ, summing it up with the second equation (4.7) and substituting equations (4.8) into them, gives the following constitutive relationship for the yielding phase: 4.9In general, different yielding factors ξ could be applied to elastic and viscous parameters. This, however, would lead to a more complex, second-order differential constitutive relationship (4.9). Adopting the yielding factors equal simplifies the model considerably and leads to a closed form solution, while still allowing for a broad class of material behaviour to be described within a thermomechanically consistent framework (Houlsby & Puzrin 2002; Puzrin & Rabaiotti 2010).
(c) Geometric expressions in the compression zone
The strain in the sliding layer in the direction of sliding is given by 4.10where δ(x,t) is the displacement of the sliding layer in the direction of the slope. From the equilibrium equation (4.5), it follows that the effective earth pressure does not vary along the compression zone p′(x,t)=p′(t). Therefore, from the constitutive equations (4.6) and (4.9), it can be concluded that the strain is also constant. In this case, equation (4.10) can be rewritten as 4.11Because the boundary of the compression zone is defined at the tip of the shear band, the relative displacement there is always the same and equal 4.12at which the shear strain in the slip surface reaches the critical value of γr and the shear resistance drops to its residual value τr (figure 8). Substitution of the equation (4.12) into equation (4.11) gives 4.13It follows that the strain in the compression zone is inversely proportional to the length of the compression zone.
(d) Differential equation
The differential equation for the strain ε in the compression zone during yielding is obtained by combining the equilibrium equation for the sliding layer (4.5) with the constitutive relationship (4.9) and the geometric expression for the length of the compression zone (4.13): 4.14where 4.15The equation for the pre-yielding phase is also given by equation (4.14) with parameters obtained from equations (4.15) by setting : 4.16
For the yielding phase, a closed form solution can be obtained in parametric representation via a free parameter y: 4.17where 4.18with ly, εy and being the length of the compression zone, the strain and the strain rate in it at the onset of yielding ty, respectively. For certain combinations of parameters, equation (4.17) describes initial acceleration of the landslide, followed by further stabilization.
When the compression zone becomes short (l≪L), equation (4.14) has an asymptotic solution. For the pre-yielding phase, this solution describes the landslide stabilization: 4.19where is the strain rate in the compression zone at the onset of yielding ty. If the yield stress p′y is not reached, the strain in the compression zone will asymptotically approach a finite value, resulting, from equation (4.13) in a finite value of the length of the compression zone. That is, propagation of the shear zone along the slip surface will asymptotically come to a halt and the landslide will become stable.
5. Determination of the model parameters
The landslide length, L≈650 m, was defined from geodetic measurements performed in 2006–2010 (figure 4). Indeed, in the lower 600–700 m of the landslide, displacement rates increase uphill, while in the upper part no significant uphill increase of displacements can be observed. Location of the boundary between the two parts of the landslide can be estimated at about L≈650 m, where displacement rates below are higher than those above. This indicates that the two parts of the slide are moving apart relative to each other, as can be confirmed by the observed tension cracks and some geological evidence of a rock outcrop in this area. As a consequence, the earth pressure p′t could be assumed equal to the average active failure pressure p′a=94 kPa, calculated after Chu (1991) using the peak angle of internal friction, φ′p=33.7°, estimated in the sliding layer from the Marchetti dilatometer tests (figure 11a) with standard deviation 2.0°. The average passive pressure is also calculated using the same peak angle of internal friction as p′p=558 kPa (Chu 1991). For these calculations, the total unit weight was taken as γ=20 kN m−3; the average depth of the sliding surface, h=20 m, was defined from the inclinometer measurements; and the average slope inclination was adopted as α=20°; unit weight of ground water is γw=9.81 kN m−3, and the depth of phreatic surface, m, was determined from the ground-water table depth measurements (figure 5) with s.d. 0.4 m.
The length of the compression zone at the onset of yielding in 1991, ly=21 m, was estimated from geodetical measurements carried out by Tschudi & Angst (1999) in 1988–1998 (figure 6). Within this zone, high compression and heaving rates of the ground and structures can be observed, significantly larger than those measured outside of the compression zone in figures 6 and 3b, respectively.
The gravitational shear stress, τg=137 kPa, was calculated using equation (4.2) with γ=20 kN m−3, h=20 m and α=20°. The residual angle of internal friction on the slip surface, φ′r=23.7°, with rate dependency of 0.4°–0.5° per log cycle of shear strain and s.d. 0.4°, was determined (figure 7) using the ring shear apparatus at loads 100–200 kPa, bounding the pressure range at the sliding surface. The residual shear strength τr=122 kPa was then calculated using equation (4.4) with γw=9.81 kN m−3 and hw=9.4 m.
Young's modulus, E=10−50 MPa, was determined using in situ Cambridge and Marcchetti dilatometer tests (figure 11b; Puzrin et al. 2008). The relatively broad range of obtained moduli values is due to the fact that Marchetti and Cambridge tests provide lower and upper limits for the stiffness, respectively.
The yearly increase in earth pressure in the compression zone of 0.8–1.0 kPa was back-calculated from the changes in the dimensions of inclinometer pipes measured in 2008–2010 using the novel inclinodeformometer (Schwager et al. 2009).
All the above physical and geometric parameters are summarized in table 1.
6. Model calibration and validation
In general, for large-scale time-dependent geotechnical problems, reliable determination of some material parameters is not possible from laboratory and field tests owing to limitations of the available experimental techniques (Van Asch et al. 2007). This problem can be solved with the help of inverse analysis, which uses, as an input, displacement monitoring data together with more reliable model parameters, in order to back-calculate the remaining parameters and use them in the model predictions (Tacher et al. 2005; Puzrin & Sterba 2006; François et al. 2007).
Applying this methodology to the St Moritz landslide, after estimating the date of the onset of yielding (around ty=4 September 1991) and the corresponding strain rate from figure 2, we observe that equation (4.19) provides a good fit (figure 12) to the stabilization branch of the time-displacement curve using just one parameter: a=0.081yr−1. Subsequently, equation (4.17) gives a reasonable fit to the acceleration branch of the curve using two additional parameters: yy=−0.34 and λ=0.37. (The deviation of the predicted curve in figure 12 from the displacement data measured in 2001–2006 correlates temporally with the fluctuations in the rate of precipitation.) This fitting allows for back-calculation of the model parameters E, p′y and ξ by substituting a,yy and λ into equations (4.6), (4.16) and (4.18) together with parameters γr, ly, L, h, d, τg, τr and p′t estimated above.
The resulting back-calculated parameters happen to satisfy some important constraints and provide an insight into the mechanism of the St Moritz landslide. For a very broad range of the critical shear strains γr=100–1000%, corresponding to realistic strains at the onset of yielding of about εy=1–10%, the back-calculated yield stress falls within a remarkably narrow range –559 kPa, very close to the estimated (passive) failure pressure p′p=558 kPa, justifying the significant drop in stiffness and viscous resistance described by the similarly narrow range of the back-calculated yielding parameter ξ=0.0147−0.0193. This remaining resistance is, however, essential for explaining relatively moderate landslide accelerations and positive pressure increments in the compression zone. Indeed, the predicted pressure increase rates for the years 2008–2010: 6.1compare well with the observed values of 0.8−1.0 kPa yr−1. Finally, the range of the back-calculated Young's moduli, while being rather broad: E=5.4–41.2 MPa, also overlaps with the 10–50 MPa (figure 11b) obtained from the field tests.
All the above model parameters and their validation are summarized in table 2. In the above analysis, the fact has been neglected that the boundary of the compression zone, which was at the onset of yielding in 1991 was at ly=21 m, has been moving downhill ever since. For the yielding phase, displacements for x=21 m point in figure 12 were calculated taking the strains in the (shrinking) compression zone and assuming they also occur in the zone between the current tip of the compression zone and the fixed position at x=21 m. This is a slight approximation (since p′ varies in the region above the compression zone), but is within the level of accuracy required, because all the soil below x=21 m has reached the yield stress p′y at the onset of yielding, and no significant pressure changes beyond this value are expected during the further landslide evolution.
7. Discussion and concluding remarks
The mechanism of the progressive propagation of the zone of intense shearing and subsequent yielding of soil in the compression zone seems to be capable of explaining the ‘paradox’ of the St Moritz landslide acceleration in spite of the absence of the ‘usual suspects’: increase in the yearly precipitation and pore water pressure, as well as rate dependency on the sliding surface. As a result, this mechanism allows for improved assessment of its hazard and suggests measures for its mitigation. The model predicts that the acceleration phase could continue until about 2025 and may cause an additional 30 per cent increase in the current displacement rates. After that the slow stabilization of the landslide could be expected, provided no increase in the yearly precipitation and pore water pressure (figures 2 and 5) and no drop in the yielding factor ξ (figure 10) take place.
Even for this favourable scenario, however, the final displacements can become excessive for some buildings, requiring an early warning system, which is currently being installed in the Leaning Tower and within the inclinometer pipes directly in the landslide body. Further stabilization of the Leaning Tower is likely to be necessary within the next 10 years. The displacement rates, however, can only be controlled by lowering the phreatic surface with the help of a drainage system. As suggested by the proposed mechanism, however, such a system would be more effective if it had been constructed in 1980s, i.e. before the earth pressure in the compression zone reached the yield stress. Thus, by predicting the time of the onset of yielding, the model can provide an important timeframe for stabilization of constrained landslides.
The contributions of I. Sterba, S. Messerklinger, S. Annen and M. Schwager, all of the ETH Zurich, to this study are highly appreciated, as well as a valuable discussion with Prof. Eduardo Alonso, of the UPC Barcelona. The work has been partially supported by the ASTRA/VSS (grant no. VSS 2005/502) ‘Landslide-Road-Interaction’.
- Received January 24, 2011.
- Accepted March 1, 2011.
- This journal is © 2011 The Royal Society