This paper has been published by the Canadian Geotechnical Journal as an e-first article on 18 June 2012, and is currently available here:
I did some work on progressive failure in large landslides in sensitive clay, drawing on principles of fracture mechanics, which have been used, but not commonly, in geotechnique to describe progressive failure processes. That work is documented in my doctoral thesis, which one can find here: Quinn doctoral thesis, 2009, and partly elaborated in a journal paper here: Quinn et al, 2011.
The journal paper essentially re-frames one of the chapters of the thesis. There is a subsequent chapter in the thesis that takes the model a step further. That chapter was submitted to the Canadian Geotechnical Journal late 2010 and was accepted with comments early 2011, but I didn’t find time to clean it up for re-submission until just recently. This post will present the draft manuscript, as submitted. Apologies for the formatting issues… there may be some symbols missing although I will try to clean them up…
Development of Progressive Failure in Sensitive Clay Slopes
P.E. Quinn (1), M.S. Diederichs (2), R.K. Rowe (3), D.J. Hutchinson (2)
(1) BGC Engineering Inc., Victoria, BC
(2) Department of Geological Sciences and Geological Engineering
(3) Department of Civil Engineering
GeoEngineering Centre at Queen’s-RMC
Queen’s University, Kingston, Ontario, Canada
The development of progressive failure in sensitive clay slopes is analysed using principles from fracture mechanics, where failure stress depends on a material property (toughness) and length of an existing weakness, unlike in traditional soil mechanics where failure stress is generally assumed to be only a function of material strength. An existing weakness, or partially developed failure surface, can propagate due to sudden loads, as might be induced by seasonal erosion cycles, small local landslides at the river bank, or earthquake shaking. The paper examines the growth of the failure surface over time due to such loads. The analysis shows that a large landslide in sensitive clay can occur after a number of seasonal erosion cycles with no obvious trigger. It also shows that a larger load due either to a small landslide or earthquake can cause a developing failure surface to propagate freely toward general collapse. Large landslides in sensitive clay are often preceded by a smaller landslide at the river bank, but can be triggered by large earthquakes, or may occur for no obvious reason. This agreement between predictions and real behaviour suggests that the model for progressive failure of large landslides in sensitive clay deserves further consideration.
Sensitive clay, progressive failure, fracture mechanics, retrogressive landslides, brittleness
Large, natural landslides in sensitive clay tend to originate at steep river banks and subsequently involve significant areas of very gently sloping terrain beyond the river bank. A typical example is shown in Figure 1, which shows the 1971 South Nation River slide of eastern Ontario. In Canadian sensitive clay slopes, failure is most often triggered by erosion (Lebuis et al. 1983), with many of the largest landslides being triggered by earthquake shaking (e.g. Aylsworth and Lawrence, 2003, Desjardins, 1980). These landslides are often described as retrogressive, developing as a series of individual rotational slumps (e.g. Bjerrum, 1955). An alternative hypothesis has been proposed for some of these landslides by several authors, including Bernander (2000), Quinn et al. (2007), and Locat et al. (2008), whereby failure is suggested to develop progressively upward from the river bank, with the complete continuous failure surface forming before the initiation of significant, noticeable movement and disruption in the ground mass above this incipient failure plane. This mechanism is explained in detail in Quinn (2009) and Quinn et al. (2011), and is illustrated in Figure 2. Phenomenological evidence for the model has been detailed in Quinn (2009) and Quinn et al. (2007).
Figure 1. South Nation River slide of 1971. (Reproduced with the permission of Natural Resources Canada 2008, courtesy of the National Air Photo Library)
Figure 2. Upward propagation of progressive failure.
The problem of progressive failure in natural slopes is a difficult challenge, defying simple analytic treatment. Limit equilibrium methods for slope stability analysis have limited applicability when geological materials exhibit strain-weakening behaviour: use of the peak strength in analyses will surely overestimate the factor of safety, and use of the residual strength will give a lower bound, but no clear idea of the time scale over which failure may occur. Similarly, finite element and finite difference models, which have revolutionized the analysis of geotechnical problems involving ductile materials, have limited applicability: release of energy into the model with strain-weakening cannot be handled deterministically, therefore model behavior after the initiation of collapse has no meaning, and development of progressive failure cannot be modeled. Linear elastic fracture mechanics has been used to develop analytic models of progressive failure in other brittle materials, in particular snow involved in dry slab avalanches. This work follows from those endeavours, and attempts to develop an analytic framework for understanding progressive failure in sensitive clay, with an emphasis on exploring triggering mechanisms.
Sensitive clays are quasi-brittle under undrained loading conditions, with a pronounced loss in strength after an initial peak for most stress paths (Lefebvre and LaRochelle, 1974). Quasi-brittle materials can fail by propagation of localised fracture (Anderson, 1995). This paper uses the principles of linear elastic fracture mechanics (LEFM) to illustrate how progressive failure can propagate in sensitive clay slopes over time, due either to an accumulation of annual erosion cycles or the incidence of a single larger perturbation, such as a moderate local landslide or an earthquake. This serves to provide a unifying analytic basis for two distinctly different landslide triggers, both of which are understood to cause large landslides in sensitive clay in eastern Canada: most landslides are associated with erosion (Lebuis et al. 1983); however, an important subset, comprising perhaps most of the largest landslides, occur less frequently and are associated with major earthquakes (Aylsworth and Lawrence, 2003). This paper suggests a mechanism consistent with both of these triggers.
1.3 Scope and limitations
LEFM concepts are not commonly used to examine failure in clay slopes. In LEFM, the analysis of the propagation of a shear fracture (i.e. shear band) depends on the length of the process zone (i.e. plastic zone, or zone within which shear stress transitions from peak strength to residual strength) ahead of the propagating shear band. This length is a material property, and is related to fracture toughness, which is a measure of the energy required to propagate fracture (Anderson, 1995). Fracture toughness is a material property that is not typically measured for soils, but the length of the process zone can be estimated from work by Tavenas et al. (1983), as shown by Quinn et al. (2011). Note, however, that the toughness, or brittleness, of sensitive clay appears to vary by at least two orders of magnitude for the Canadian sensitive clays tested by Tavenas et al. (1983). The present work would therefore benefit from specific laboratory testing of fracture toughness of sensitive clay.
This paper examines the propagation of fracture (i.e. shear bands) approximately parallel to the ground surface adjacent to a natural cut slope (e.g. river bank), and determines the conditions for unlimited propagation (i.e. general collapse due to propagation of fracture). The work first presents the conditions for collapse, as interpreted from an LEFM analysis of progressive failure of overconsolidated clay slopes by Palmer and Rice (1973). This analysis establishes the stress conditions ahead of a developing shear band necessary to propagate failure. Next, the stress conditions ahead of a shear band due to annual erosion cycles at the river bank and due to potential large earthquakes are determined. These estimates of stress are then used to determine how quickly failure can develop under several typical slope scenarios with varying rates of river erosion and different levels of earthquake shaking. The findings of this paper are believed to be generally applicable to the progressive failure of sensitive clay slopes, and may be broadly applicable to other types of slopes with different quasi-brittle materials other than sensitive clay.
This paper draws heavily on concepts developed by Rice (1968) and Palmer and Rice (1973), both of which make several simplifying assumptions. In particular, both of those references assume the propagating shear band is long with respect to problem geometry, and also in relation to the nominal displacement of the brittle material. This paper includes examination of shorter shear bands early in the development/initiation of progressive failure. That phase of failure, as described herein, may be considered to be resting on a less stable theoretical foundation than the subsequent propagation of failure at the end of a longer shear band.
2 FRAMEWORK FOR EVALUATING THE POTENTIAL FOR SHEAR BAND PROPAGATION
In fracture mechanics, the potential for failure is a function of three factors: driving stresses, fracture toughness (i.e. a measure of the energy required to propagate fracture), and the length of an existing weakness, or degree of accumulated damage (Anderson, 1995). Failure propagates over time. This approach differs from classic soil mechanics approaches which tend to assume limiting equilibrium conditions are reached at the point of collapse, with failure occurring simultaneously along the full failure surface. Classic soil mechanics problems are generally solved by comparing working stresses with yield stress, whereas in fracture mechanics, both toughness and defect size must be considered in relation to working stresses.
Fracture toughness of the material determines how applicable LEFM concepts are for analysing failure (Anderson, 1995). Very brittle materials (i.e. materials with low fracture toughness) will fail by brittle fracture, and LEFM applies. Very tough materials are subject to instantaneous collapse rather than progressive failure. Conventional limit equilibrium analysis applies to these materials, and LEFM concepts are not applicable. Materials with intermediate toughness will also fail progressively over time, but non-linear fracture mechanics methods may be necessary. The materials with intermediate or low toughness can fail with average driving stresses lower than the peak strength of the material. This phenomenon is typically observed in slopes that fail progressively. Bishop (1971) suggested that the difference between the actual safety factor back calculated based on the mean failure stress along the failure surface and that calculated from the measured peak strength depended on the brittleness of the clay, defined in terms of the amount of strength loss after the peak. Sensitive clays can be very brittle, with significant loss of strength over relatively small strains or displacements. It therefore follows that progressive failure is possible in these materials, and that fracture mechanics approaches should be valid for their analysis.
Figure 3. Conceptual model for progressive failure of clay slopes, modified from Palmer and Rice (1973).
Palmer and Rice (1973) used a method developed by Rice (1968) to analyse the propagation of progressive failure in overconsolidated London clay slopes, using the simplified model of a step excavated into an infinite slope to represent a road or railway cut in a natural clay slope. Their model is illustrated in Figure 3. It considers the existence of a shear band of length, L, propagating upward from the toe of the cut, with height, h. If the slope angle, alpha, is made very small, this same model approximates the erosion of a river channel in the gentle Champlain clay plains of eastern Canada.
One can obtain relationships for critical conditions leading to continued propagation of a shear band, using the approach of Palmer and Rice (1973), for a slope under static loads with the geometry given in Figure 3, as follows:
 Lcr = [sqrt(2E’h(del_bar)(tao_p – tao_r))]/(tao_D – tao_r)
Lcr is the critical length. Once the shear band reaches this length, it will propagate with no additional external load;
E’ is the Young’s modulus of the clay;
del_bar is the nominal displacement of the clay. This is a material property defined by:
del_bar(tao_p – tao_r) = integral(tao – tao_r)d(del)
This value is some fraction of the shear displacement along a shear surface required to fully remould the clay, del_r, and it may be considered an indirect measure of fracture toughness;
tao_p and tao_r are the peak and residual shear strengths;
po is the initial at rest lateral earth pressure in the slope, averaged over the height h, before the introduction of the cut; and
tao_D is the shear stress along the line of the potential shear band due to gravity or other driving forces, as shown in Figure 3.
All of these variables can be determined with some certainty from standard laboratory tests and slope geometry except for the nominal displacement, del_bar. That value requires an understanding of the complete stress-displacement curve to the point where the fully remoulded strength has been reached. High quality, large displacement shear tests on sensitive clay are scarce in the literature. Stark and Contreras (1996) performed ring shear testing on Drammen clay, but their total displacement of nearly 60 mm was likely insufficient to reach the fully remoulded strength, since the observed residual strength was only roughly half the peak strength. Sensitive Drammen clay would be expected to have a much lower remoulded strength, so additional strength loss may have been observed if the test had been continued to much larger displacements. A careful examination of novel test data presented by Tavenas et al. (1983) suggests that values ranging between 50 mm and 5 m might be an appropriate estimate of nominal displacement for a range of Canadian sensitive clays (Quinn et al. 2011).
The nominal displacement, del_bar, varies with brittleness (or with the magnitude of displacement/strain required to obtain the residual strength) as illustrated in Figure 4. This value is higher for materials with more ductile behaviour, and decreases with increasing brittleness.
Figure 4. Nominal displacement, (del_bar), and brittleness.
The utility of Equation  can be explored by considering a typical river bank carved in the gentle sensitive clay plains of eastern Canada. Consider a bank height of 30 m with the clay plains above the crest sloping at 0.5 degrees toward the river bank. If the nominal displacement of the clay, del_bar, is 0.5 m (e.g. Saint Alban clay, from Tavenas et al. 1983, as interpreted by Quinn et al. 2011), bulk density is 18 kN/m3, and peak and remoulded strengths are 80 and 1 kPa (note that sensitivity of Saint Alban clay tested by Tavenas et al. 1983 was reported as ranging between 24 and > 600, with ~ 80 being perhaps a typical value), respectively, then the critical length, Lcr, can be calculated as 220 m. Note that this assumes E’ = 10,000 kPa and a friction angle of 30 degrees, leading to Ko = 0.5 and associated average lateral earth pressure, po = 135 kPa, assuming normally consolidated conditions and gently sloping ground. This means that if a shear band grows, through some mechanism, to a length of 220 m, it will begin to propagate under no additional external load, resulting in general collapse of the slope, according to LEFM principles and the analysis of Palmer and Rice (1973).
The presence of additional driving stresses due to some external load would serve to reduce the critical length. For example, an additional shear stress of 50 kPa reduces the critical length to 15 m. Therefore, if an earthquake, or some other perturbation, generates additional shear stresses along the line of the potential shear band on the order of 50 kPa, an existing shear band of any length greater than 15 m will begin to propagate to failure. Similarly, a load at the river bank due to erosion may generate additional shear stresses at the end of the 15 m long shear band on the order of 50 kPa. In this case, the shear band can also begin to propagate, but only to the extent within which the additional shear stresses exceed 50 kPa, and only while such additional stresses exist. The critical lengths for various additional shear stresses for the example slope are summarized in Table 1.
Table 1. Reduced critical lengths with additional driving stresses.
The following sections of this paper examine the magnitude of shear stresses that would develop along the line of a potential shear band due to erosion cycles and large earthquakes. This will serve to establish the conditions under which an existing sensitive clay slope might suffer a large landslide by progressive failure propagating upward from the river bank.
3 ADDITIONAL STRESSES AT THE TIP OF A PROPAGATING SHEAR BAND DUE TO EROSION
Consider a semi-infinite solid loaded at the free end, with an existing discontinuity (i.e. crack, or shear band) as illustrated in Figure 5. It has been shown (e.g. Rice 1968) that this loading condition will result in a stress concentration at the end of the crack, and finite element analyses confirms the result. The inset in Figure 5 shows a typical distribution of shear stresses near and ahead of the end of the crack under these loading conditions. The crack can either be empty or filled with a softer material than the surrounding plate, and in either case a similar result will be obtained. It follows that a slope with an existing weak zone, or propagating shear band, might have similar stress concentrations at the end of the weak zone. Finite element analyses were used to estimate the magnitude of such stresses under typical loading conditions. Note that the use of finite element analysis follows from a suggestion by Palmer and Rice (1973), as they noted that a number of simplifying assumptions were necessary in developing their theoretical framework. The use of a small, finite thickness weak zone in the present analysis, instead of an arbitrarily thin shear band, differs from the original models of Rice (1968) and Palmer and Rice (1973). The validity of its use is investigated in the preliminary finite element analysis, which shows similar behavior for both very thin empty cracks and thicker, but still relatively thin, weaker zones filled with low stiffness material. Use of the finite thickness weak zone is felt to better represent the physical situation under study in this paper, as the liquefaction of sensitive clay is not necessarily expected to lead to propagation of arbitrarily thin shear bands.
Figure 6. Typical slope model for FEM analysis, 45 degree slope, 30 m high, with
150 m long existing discontinuity. Main image shows global boundary conditions, bottom image shows mesh details near the slope and weak zone, and inset shows details of 10 m excavation or small landslide.
Figure 6 shows a typical model slope, with a sample loading condition, in this case an excavation representing a small erosion episode or slope failure 10 m by 10 m into the river bank. The model boundaries were located 150 m to the left of the slope, 1000 m to the right, and 250 m below the toe of the slope. Lateral boundaries were restrained by rollers, and the lower boundary was fixed. Each model had approximately 12,000 to 15,000 elements. It should be noted that a number of simplifying assumptions have been made in this analysis, and therefore the results are meant to be illustrative of general behavior, rather than predictive of precise stress-deformation conditions.
The tip of the weak zone is shown in the inset of Figure 6. The slope is 30 m high, standing at 45 degrees, fairly typical values for river banks in Canadian sensitive clay deposits. The gentle plains above the river bank slope at 0.5 degrees toward the river. An existing weak zone 0.3 m thick extends 150 m inward from the toe of the slope, extending at 0.5 degrees, parallel with the ground surface. An excavation is made at the river bank after gravity loading, and in this case extends 10 m inward, resulting in an excavation volume of 50 m3. The upper 10 m of the model below the ground surface is comprised of a stiff clay crust, with a Young’s modulus, E, of 20,000 kPa, reflecting an undrained shear strength of approximately 160 kPa and Poisson’s ratio of 0.49. The unit below the stiff clay crust is a softer clay with E = 10,000 kPa. The weak zone is 100 times softer, with E = 100 kPa, corresponding to fully remoulded conditions along the propagating shear band for a remoulded strength of 1 kPa.
The model was analysed for weak zones with lengths, L, of 50, 80, 100, 150 and 200 m, and for excavations with length, Lex, of 3, 5, 10, 20 and 21.5 m. These excavation lengths correspond with excavation volumes, Vex, of 4.5, 12.5, 50, 200 and 231 m3, respectively, since the assumed excavations are triangular with height approximately equal to Lex.
The dimensions, L and Lex, are illustrated in Figure 6. Variation of these parameters allowed the interpretation of relationships between stress at the end of the weak zone and excavation size or length of weak zone. All analyses were conducted assuming elastic material behaviour, which is appropriate for small strain undrained loading of clay, and is considered suitable for the purpose of estimating the stress concentration associated with sudden loads at the river bank, such as annual erosion cycles associated with the spring freshet, and secondary erosion cycles associated with the wet fall season (Lebuis et al. 1983).
Finite element analyses show a bulb of elevated shear stress at the end of the weak zone for all excavations at the river bank, with peak stress occurring about 0.1 m ahead of the tip, and stresses dropping with additional distance. Peak stress ahead of the tip increases with the volume of excavation (i.e. 0.5[Lex]^2, where Lex is illustrated in Figure 6), and decreases with length of the weak zone (i.e. L, as illustrated in Figure 6) in non-linear fashion. The rate of decrease of shear stress away from the tip past the peak value at 0.1 m is independent of both excavation size and length of the weak zone. Figure 7 shows typical distributions of shear stress ahead of the tip of the weak zone. All stress distribution curves have the same general shape, with stresses varying according to magnitude of driving stress (i.e. size of the excavation) and length of the weak zone, and tapering off away from the tip according to the same general function of distance.
Figure 7. Shear stress beyond the end of the weak zone (100, 150 and 200 m long weak zones, 231 m3/m excavation, Young’s modulus 10,000 kPa in soft clay, 100 kPa in weak zone).
Figure 8. Peak shear stress ahead of the tip versus excavation volume, Vex (150 m long weak zone), with interpreted best fit power law trend.
Figure 8 shows the relationship between peak shear stress ahead of the tip and volume of excavation at the river bank, Vex. Figure 9 shows the relationship between peak shear stress ahead of the tip and length of the weak zone, L, with peak shear stress varying according to a power law. Figure 10 shows the relationship between shear stress and distance ahead of the weak zone for different length weak zones, L.
Figure 9. Peak shear stress ahead of the weak zone versus length of the weak zone, L (21.5 m excavation length), with interpreted best fit power law trend.
Figure 10. Shear stress distribution beyond the end of the weak zone for different length (L) weak zones (Lex = 21.5 m excavation).
Synthesis of the results shown in Figure 8 to Figure 10 allows the interpretation of the following relationship for shear stress (in kPa) ahead of the tip of the weak zone:
 tao_xy(x) ~ 2700(Vex^0.56)(L^-1.5)(x^-.45)
where Vex = volume (m3/m) of excavation at the river bank (i.e. erosion) associated with one erosion cycle;
L = length (m) of the weak zone; and
x = distance (m) ahead of the tip of the weak zone.
This relationship is valid for x > 0.1 m. Shear stresses drop closer to the tip of the shear band, and so the relationship is not valid at shorter distances. Equation  can be reorganized to obtain a relationship for the distance within which a specified shear stress will be observed ahead of the weak zone:
 x_tao = [(2700(Vex^0.56))/(tao_xy(L^1.5))]^2.2
Equation  can be used to estimate the shear stresses at any point directly ahead of the shear band, and Equation  can be used to estimate how far a given value of excess shear stress is observed ahead of the shear band, for any combination of shear band length and excavation size. These will be used subsequently to explore the potential for shear band propagation over time due to an accumulation of small load cycles or the sudden occurrence of one large perturbation.
4 ADDITIONAL STRESSES IN A SLOPE DUE TO EARTHQUAKE LOADING
The Saint Lawrence Lowlands of eastern Canada, which contain extensive deposits of sensitive Champlain clay, are confined largely within the Saint Lawrence and Ottawa-Bonnechere grabens, which are rift valleys bounded by zones of normal faulting (Tremblay et al. 2003). These areas are seismically active, so that the areas of sensitive clay in eastern Canada are subject to elevated seismic hazard in comparison with much of the rest of Canada. Figure 11 shows design spectral acceleration values for a period of 1 second for eastern Canada from Halchuk and Adams (2008). The locations of a number of landslides believed to have been triggered by earthquakes are also shown in that Figure.
Figure 11. Spectral acceleration for a period of 1.0 seconds in southeastern Canada at a probability of 2%/50 years for firm ground conditions (NBCC soil class C), redrawn from Halchuk and Adams (2008), with locations of landslides known to have been triggered by earthquakes.
Many of the largest landslides in sensitive clay in eastern Canada are thought to have been triggered by large earthquakes (e.g. Aylsworth and Lawrence, 2003 and Desjardins 1980). This includes a number of extremely large ancient landslides east of Ottawa dated to a probable single event some 4500 years before present (YBP) according to Aylsworth and Lawrence (2003), several very large features near Shawinigan dated by Desjardins (1980), and a number of large landslides thought to be associated with the 1663 Charlevoix earthquake, including a very large landslide within which the 1971 Saint-Jean-Vianney slide occurred (Leggett and LaSalle, 1978). The 1663 earthquake is believed to have been magnitude 7 (M7), whereas a 7000 YBP event that triggered a small number of additional landslides east of Ottawa has been tentatively estimated as M6.5 or greater (Aylsworth and Lawrence, 2003). Aylsworth and Lawrence (2003) have estimated that M5.9-6.0 represents an approximate lower threshold for the occurrence of large earthquake-triggered landslides based on observations of landslides associated with a number of moderate (M5.3 to M6.3) historic eastern Canadian earthquakes. For example, Lefebvre et al. (1992) identified a number of small landslides in southern Quebec associated with the 1988 Saguenay earthquake (M5.9). That event did not trigger any known large retrogressive earthquakes.
The 2005 version of the National Building Code of Canada (NBCC, National Research Council 2005) identifies uniform seismic hazard spectra for Canadian sites, specifying an annual probability of 2% in 50 years (1/2500 per annum) for structural design. Atkinson and Beresnev (1998) developed simulated ground motion time histories that could be used for dynamic modelling purposes based on the design criterion of 10% in 50 years (1/500 per annum) in the previous (e.g. 1995) version of the NBCC. These time histories have appropriate amplitude and frequency content for the old uniform hazard spectra. In the Ottawa area, for example, an M5.5 earthquake 30 km away (M5.5R30) scaled at 0.8 (i.e. amplitude of the acceleration time history is reduced to 0.8 times that of the full record) addressed the short period component of the design uniform hazard spectrum, and M7R150 at 0.7 addressed the long period component. Additional time histories updated by Atkinson and Beresnev in 1999 for the 2005 Code were obtained from Atkinson (personal communication 2008). For Ottawa, for example, M6R30 x 0.8 and M7R70 x 0.8 combine to address the long and short period components of the 2005 design hazard.
The earthquake events thought to have triggered the large landslides mentioned in the preceding paragraphs were likely larger than these design events, or rather were perhaps similar to these events without reduction by the scaling factor. The earthquake records (i.e. M6R30 and M7R70), both with and without scaling factor reductions, are therefore believed to represent a reasonable first estimate of the lower threshold of earthquake shaking likely to lead to widespread occurrence of large landslides. These events are therefore used for much of the analyses reported in this paper.
The simulated earthquake records were used to estimate the amplification of ground shaking upward through the soil column, with a particular emphasis on obtaining the magnitude of cyclic shear stress at the elevation of a potential shear band propagating inward from the toe of a river bank. Several different bank heights, and corresponding depths of shear band, were considered, and several different soil profiles were examined. Four different cases were considered for most scenarios: M6R30 x 0.80, M6R30 x 1.0, M7R70 x 0.80, and M7R70 x 1.0.
Additional cases of strong motion were considered to represent probable shaking at known landslide locations due to the 1663 Charlevoix earthquake (M7). Landslides at Saint-Jean-Vianney and Shawinigan occurred about 105 and 225 km, respectively, from the earthquake. M7R105 and M7R225 records were not available from Atkinson (pers. comm. 2008), so the M7R100 record was used, with scaling factors of 0.85 and 0.32, respectively, to obtain the approximately correct peak ground acceleration, as determined by comparing peak ground acceleration with distance for the various distance scenarios for the M7 event. Note that this approximation should be reasonable for Saint-Jean-Vianney, but the frequency content of the record might not be appropriate for the Shawinigan location. Nevertheless, it serves as a useful case for first examination.
Numerous stratigraphic profiles were examined to reflect the variability in soil conditions throughout the Saint Lawrence Lowlands. Actual overburden thicknesses at the Aylsworth and Lawrence (2003) landslides is approximately 24 +/- 14 m near the headscarps, as interpreted from thickness data provided by Natural Resources Canada (2004). Overburden thickness at the Saint-Jean-Vianney and Shawinigan landslides has been similarly estimated as 44 +/- 5 m and 40 +/- 6 m, respectively. Total thickness of the modelled soil column was therefore varied from 20 to 60 m to cover this range of soil thicknesses. The soil column was either soft clay directly over bedrock, or an upper soft clay unit with stiffer soil (e.g. stiff clay or till) between the upper soft clay and underlying bedrock.
Shear wave velocities for soft clay were estimated based on the following relationship from Hunter et al. (2007):
 v_s = 110.3 + 5.179z^0.871 +/- 55.4
where z is depth below the ground surface, in metres, and shear wave velocity, v_s, is in m/s. Bedrock was assumed to have a shear wave velocity of 2500 m/s. Stiffer soils between the soft clay and bedrock, when included, were assumed to have a shear wave velocity of 500 m/s at the top of the stiffer unit, increasing by 5-10 m/s per m of thickness.
The reduction in shear modulus with strain (i.e. G/Gmax) was based on the relationship obtained for Leda clay by Anderson and Richart (1976). The relationship between damping and strain was based on the average relationship for clay reported by Seed and Idriss (1970).
All analyses were conducted using SHAKE91 (Idriss and Sun, 1991), accessed through the SHAKE2000 interface (Ordonez, 2008). SHAKE91 computes the response of a semi-infinite horizontally layered soil deposit overlying a uniform half-space to vertically propagating shear waves, using a one-dimensional wave propagation method (Schnabel et al. 1972). The program derives a continuous solution of the wave equation adapted for fast Fourier transform techniques. The analysis is performed in the frequency domain, and the analysis is linear. An iterative procedure is used to account for non-linear behaviour, as would be encountered for real soils under strong ground motions. Non-linearity is captured in the model by using strain-dependent values for shear modulus and damping. Initial values of shear modulus and damping are assumed in the first iteration, and the analysis is repeated until strain-compatible modulus and damping values are reached. A strain-compatible solution is generally reached within five to eight iterations. Additional detail is available in Idriss and Sun (1991) and Schnabel et al. (1972).
Earthquake motions propagating upward from bedrock through a column of very soft soil may amplify, resulting in higher peak ground accelerations at the ground surface, in comparison with the motion at the underlying bedrock. Amplification of peak ground acceleration for different thicknesses of soft clay over bedrock are shown in Figure 12 for the long period (i.e. M7R70) and short period (i.e. M6R30) events at full strength and scaled to 0.8. Note that the variation with depth is not a uniform function, but rather has peaks at specific depths, and the trends are different for the two different events. Response to upward propagating shear waves depends on the frequency and amplitude content of the input motion as well as the resonant characteristics of the soil column. For a given input motion with certain dominant wavelengths, the response will peak for different soil column lengths as a function of both input wavelength and soil stiffness.
Figure 12. Amplification of peak ground acceleration versus thickness of soft clay for various earthquake magnitudes (e.g. M6 and M7) at various distances (e.g. R = 30 or 70 km) and different scaling factors (e.g. 0.8 or 1.0).
Ground amplification may result in large shear stresses in the soil column, varying with depth, strength and frequency content of ground shaking, and characteristics of the soil column (e.g. thickness, stiffness). A typical profile of earthquake-induced shear stress is presented in Figure 13, which shows peak shear stress with depth for a M6R30 event scaled to 0.80 applied to bedrock below 40 m of soft sensitive Champlain clay. Shear stresses in the region of interest, namely 20 to 40 m, which corresponds to typical river bank heights, range from about 60 to 90 kPa. Similar values are observed for depths of 10, 20 and 30 m for the four primary design events (i.e. M6R30 and M7R70 x 0.8 and 1.0), but are not shown. Figure 14 shows shear stress at 20 m depth for different thicknesses of soil column for the four design events. The shear stress relationships for depths of 10 and 30 m have similar magnitudes. Note that shear stress generally increases with depth (i.e. from 10 to 30 m), but the relationship with thickness of the soil column is more complex, varying both with thickness and the nature of ground shaking (i.e. the long period versus short period design events). This is a result of the changing resonant characteristics of the soil column with stiffness and thickness, as discussed earlier in relation to ground amplification.
Ground amplification and the resulting shear stress profile also vary with the nature of the soil profile, in addition to its thickness. Figure 15 shows shear stress at 20 m depth for 30 m of soft clay that is either directly over bedrock, or with up to 30 m of stiffer soil (stiff clay or till) between the soft clay and bedrock. It can be seen that shear stresses are higher – by close to 20 % – for the case of 30 m soft clay over 30 m of stiff soil and the M6R30 event. The difference is less significant for the M7R70 event, which has a greater proportion of low frequency content. Figure 16 shows shear stress with depth for two different profiles of 60 m thickness: 30 m soft clay over 30 m stiff soil, and 60 m of soft clay directly on rock. Shear stresses are greater for the mixed profile in comparison with 60 m of soft clay, and the difference is substantially greater for the M6R30 event, which has a greater proportion of high frequency content than the M7R70 event.
Figure 15. Shear stress at 20 m depth with stiff soil of different thickness between soft clay and bedrock.
Figure 16. Shear stress with depth for 60 m soft clay over rock and 30 m soft clay plus 30 m stiff soil over rock.
Consider again the natural slope in sensitive clay described in section 2, comprised of a river bank carved in otherwise gentle clay plains. The river bank is 30 m high and stands at about 45 degrees. The clay plains above the river bank slope toward the river at 0.5 degrees. At the elevation of the toe of the slope, peak and residual shear strengths are 80 kPa and 1 kPa respectively. The nominal displacement of the clay, del_bar, is 0.5 m, which is similar to the value interpreted by Quinn (2009) for Saint-Alban clay from data from Tavenas et al. (1983). Equation  yields a critical length, Lcr, of 220 m. This implies that an existing shear band of that length would propagate further under no additional external load. Using Equation  again, one can calculate the critical length for the same slope under additional stresses near the shear band tip, due either to local loads at the river bank (e.g. erosion cycles) or an earthquake, for example. A range of values are presented in Table 1, and these are shown graphically in Figure 17. A power law best fit curve suggests that for this scenario:
 tao_D = 1150L^-1.23 – 1.46(L/Lcr)
where tao_D is the additional driving stress (in kPa) required to initiate failure propagation for a given length, L (in m), of the shear band. The second term is included to force the driving stress to zero in the best fit power law when length, L, reaches the critical length for the slope. That is, the calculated additional driving stress required to propagate failure should be zero when the shear band is 220 m long.
It is now possible to examine shear band propagation under different loading conditions. Consider the same slope described above, but assume the existence of a 10 m long pre-existing weakness, as a first example. The occurrence of a sudden erosion event will cause an undrained load to be transmitted to the end of the shear band. Assume the erosion event has removed a volume of soil from the slope extending 1 m into the river bank, which equates to 0.5 m3 of soil removed per metre length along the river. Given the prior existence of the 10 m shear band, Equation  shows that an additional shear stress of 67 kPa is sufficient to propagate failure. Equation  indicates that this shear stress is exceeded for 1.9 m ahead of the shear band tip, given the current shear band length and the magnitude of the erosion cycle, during the period of earthquake shaking. This erosion event will therefore cause the shear band to grow by 1.9 m as a result of this sudden load. Note that shear band propagation will only occur when the additional stresses are present, thus continuation of propagation to general collapse is not necessarily guaranteed. In this case, a relatively small erosion event (0.5 m3) has caused an existing weakness of moderate length (10 m) to grow by nearly 20 %. Similar analyses for different loading conditions and different lengths of weak zone show that a weak zone of any length will grow incrementally by some distance provided the load is sufficiently large.
Now consider the same slope, but starting with a 1 m long shear band, to examine growth of an existing weakness of negligible length, as might exist at the onset of river erosion. This slope will be subject to a large number of load cycles over time. Assume these load cycles occur with a random distribution of size, Lex (i.e. horizontal length of erosion into the river bank for a triangular excavation, as shown in Figure 6), corresponding to a Poisson distribution with maximum value 3 m and mean value 0.83 m. This corresponds with an average volume of 0.34 m3/m per erosion cycle, and maximum volume of 4.5 m3/m.
It is possible to calculate the incremental growth of the shear band due to each load cycle, and plot the progress of its growth, following the same logic described for the 10 m weakness previously. Figure 18 shows the progressive growth of the shear band with an accumulation of random load cycles over time. In this specific case, the shear band propagates to failure after about 260 load cycles. Figure 18 also shows the additional shear stress due to external forces required to propagate failure, and this value decreases as the shear band grows. If a single erosion cycle occurs each year, with a random size according to the assumed distribution, then this slope would be expected to fail suddenly about 260 years after the initiation of this pattern of erosion.
It can be seen that the growth of the shear band occurs in three distinct stages. Initially, when the shear band is short, propagation occurs very quickly. Subsequently, it slows, stabilizing to a roughly constant rate. Finally, as the shear band approaches its critical length for conditions under no external loads (i.e. 220 m), the growth rate accelerates until failure occurs suddenly after a final small load increment. At this point, the volume of accumulated erosion (112 m3) is the equivalent of approximately 3.7 m of lateral migration of the river bank (i.e. assuming the total eroded volume forms a parallelogram, and migration distance is 112 m3/30 m), if one assumes the river bank maintains the same approximate profile over time as it advances inward with the accumulation of individual erosion cycles, which may affect different parts of the bank each year. Note that since the river bank will have moved over time due to the accumulated erosion, the actual length of the weak zone would be only about 220 m after about 260 load cycles. Actual failure might therefore be delayed somewhat beyond 260 load cycles. However, this method is intended to estimate approximate length of the weak zone and approximate number of cycles to failure, and the potential error does not affect the validity of the general concept being illustrated here.
Consider now three different rates of erosion for the same slope and same initial weakness 1 m long: mild, intermediate, and aggressive. The mild erosion case corresponds to that described above, and the intermediate and aggressive erosion cases correspond to average annual erosion loads equivalent to excavating volumes of 0.6 and 1.0 m3, respectively. Figure 19 shows the distribution of random load cycles for these mild, intermediate and aggressive erosion cases applied to the model. Figure 20 shows the interpreted growth of the shear band with accumulated erosion cycles for these three erosion scenarios. The intermediate and aggressive erosion cases result in failure after 140 and 90 load cycles, respectively, corresponding to total erosion of about 105 m3/m (being slightly higher for the intermediate case), or the equivalent of about 3.5 m of bank migration. Failure therefore appears to depend on the total volume of erosion at failure, which is nearly independent of rate of erosion, decreasing slightly with a moderate increase in erosion rate (i.e. 112, 107 and 104 m3 for the mild, moderate and aggressive erosion cases). Failure occurs more quickly for the faster erosion rates, varying with L_in^-2.2, where L_in is the mean annual inward advance of the river bank due to an accumulation of erosion cycles, or varying nearly linearly with the inverse of the mean annual volume of erosion.
Consider now the occurrence of a small landslide, equivalent to a single excavation removing 50 m3 of soil per m along the river (i.e. landslide extending 10 m into the river bank), occurring after 80 load cycles under otherwise mild erosion conditions. The growth of the shear band under this scenario is shown in Figure 21, in comparison with the aggressive erosion conditions described previously. This event causes a sudden growth of the shear band, and reduces the duration to ultimate failure from about 260 to 110 load cycles. Note that if this landslide involved 112.5 m3/m instead of 50 m3/m (i.e. landslide extending 15 m into the riverbank instead of 10 m), the induced shear stress would be sufficient to trigger ultimate failure, given the length of the shear band at that time.
Careful examination of Figure 18 shows that shear stresses extending throughout the slope, well beyond a small stress concentration bulb at the end of the shear band, due to earthquake shaking, could trigger failure at any time, provided the induced shear stresses are high enough. Further, the threshold shear stress for earthquake triggering of large landslides reduces over time as the shear band grows due to annual erosion cycles. Figure 21 shows the growth of a shear band under intermediate erosion conditions, but in this case, the slope is subjected to strong shaking after 125 erosion cycles due to a M6 earthquake located 30 km away. The induced cyclic shear stress is sufficient to propagate failure, and in fact an M6R30 earthquake, a very significant event generating 94 kPa shear stress at 30 m depth, would be sufficient to trigger failure after only four erosional load cycles for the example slope. By contrast, an M7R100 x 0.32 event, corresponding to the magnitude of shaking during the 1663 Charlevoix earthquake at Shawinigan, 225 km from the epicentre, would result in 27 kPa for the example slope, assuming 40 m of soft clay over bedrock. This stress would be sufficient to trigger failure of this example slope after 12 erosion cycles had occurred.
This paper has used a linear elastic fracture mechanics approach to describe progressive failure in sensitive clay slopes. The potential for propagation of failure depends on a number of factors as outlined in Equation , including those relating to slope geometry (i.e. slope height, active earth pressure and driving stress) and soil properties (i.e. stiffness, peak and residual strength, and nominal displacement). This model involves consideration of brittleness, or toughness, of the soil, as expressed by the nominal displacement, , introducing a different dimension to the analysis of failure beyond those typically considered in soil mechanics problems (i.e. available strength and working stresses only).
Elastic finite element modelling confirms the prediction of analytical models that a load applied at a river bank will be felt as a concentrated load at the end of a weak zone, or shear band, extending inward from the toe of the slope. The shear stress ahead of the shear band is an important consideration in the evaluation of potential for continuing propagation of the shear band. This stress decreases with distance ahead of the tip as an inverse power law, and is a linear function of the peak stress observed just ahead of the shear band tip. This peak stress varies according to a power law with the volume of a single erosion event into the river bank, and with the length of the shear band as an inverse power law. These relationships can be combined to develop a unifying relationship between shear stress and distance ahead of the shear band given the magnitude of erosion and length of the shear band, as presented in Equation . This relationship can be rearranged (Equation ) to determine the distance ahead of a propagating shear band within which some specific stress will be felt, given an erosion event of known magnitude.
The model has been applied to a typical slope to determine the rate of shear band propagation under different erosion scenarios. The additional stresses at the end of the shear band due to an annual erosion cycle can force the propagation of failure by some incremental distance. The shear band growth rate, expressed in incremental magnitude of growth per load cycle, is seen to be a function of the erosion rate, varying nearly linearly with the inverse of annual volume of erosion. The rate of shear band growth is also very sensitive to the brittleness of the clay, being much higher for lower values of del_bar. A single small landslide larger than the typical annual erosion event may also result in either a sudden increase in shear band length, or could serve as the final trigger for catastrophic failure.
Strong ground shaking due to earthquakes can cause substantial additional shear stress in a slope. The magnitude of this stress depends on frequency and amplitude content of the earthquake motion, thickness of soil over bedrock (i.e. generally increasing with increasing thickness), and potential existence of a stiffer soil unit between the upper soft clay and underlying bedrock (i.e. generally increasing with the presence of a stiffer interface). Slopes with an existing weakness, or shear band, of finite length may be subject to sudden failure due to the occurrence of moderate to large earthquakes. The potential for failure depends on the stage of development of the existing shear band, and the brittleness of the soil. The potential for a given earthquake to trigger large landslides therefore depends on its magnitude as well as the current state of erosion along rivers within the range of influence of the earthquake.
These general findings lead to some interesting observations. The pattern of shear band growth illustrated in Figure 18 and Figure 20 shows three distinct stages of development, and resembles the three stages of creep. Sensitive clay slopes are known to creep to failure by local rotational landslides, as documented by Lefebvre and LaRochelle (1974) and Lefebvre (1981). The pattern of shear band growth over time also resembles the growth of fatigue cracks under cyclic loads, as documented by, for example, Weibull (1963). The process of accumulating load cycles at a river bank adjacent to an existing shear band or weakness is a reasonable analog for a metal test specimen with an existing crack subjected to repetitive loading, so this similarity in predicted behaviour is very interesting. The final stage of shear band growth deserves additional consideration. The model predicts a gradual acceleration in shear band growth, with catastrophic failure occurring some time after shear band growth begins to accelerate. This is a common phenomenon with slow moving large landslides, suggesting the model may be broadly applicable for other types of large progressive landslides.
A model for progressive failure of sensitive clay slopes has been presented, based on principles of linear elastic fracture mechanics. Application of the model, using typical slope conditions for river banks carved in the sensitive Champlain clay of eastern Canada, suggests a number of key findings:
• The incidence of large landslides should be more frequent along rivers with greater erosion rates;
• Large landslides should be more common for more brittle clays;
• Large landslides can occur:
o suddenly with no obvious trigger, as a result of an accumulation of small annual load cycles due primarily to erosion;
o soon after a small local landslide at the river bank; or
o following an earthquake or other large transient load (i.e. due to other natural or anthropogenic causes, such as pile driving or blasting).
These phenomena are all commonly observed in relation to large landslides in sensitive clay, suggesting the model is worthy of further study. The nature of the soil stratigraphy, considering both thickness and the presence of a stiffer intermediate unit between bedrock and soft clay, has a strong influence on the potential for large earthquake-triggered landslides. Further, the pattern of growth of the failure surface over time, preceding ultimate failure, bears a strong resemblance to other natural phenomena, including creep failure of natural slopes, and the sudden acceleration to catastrophic collapse observed in slow moving large landslides. This suggests the model may have broader applications beyond sensitive clay slopes.
The first author gratefully acknowledges funding through NSERC PGS and OGSST scholarships. This work was supported financially by the GEOIDE network and the Railway Ground Hazards Research Program, funded by NSERC, CN Rail, CP Rail and Transport Canada. The work benefited substantially from the leadership and interest of Mario Ruel at CN Rail. Initial inspiration for this work came from discussions with Stig Bernander, private consultant, Göteborg, Sweden who has written extensively about progressive failure in sensitive clay slopes, and with Louis Delmas, a graduate student studying snow avalanches at the Norwegian Technical University in Trondheim. Consideration of earthquake triggering in this paper was done at the suggestion of Jim Hunter, Geological Survey of Canada. John Adams of the Geological Survey of Canada provided advice regarding seismic design parameters, and Gail Atkinson provided copies of her simulated earthquake records that conform to the NBCC design criteria.
Anderson, T.L. 1995. Fracture mechanics: fundamentals and applications. CRC Press LLC, Boca Raton, FL.
Anderson, D.G. and Richart, F.E. 1976. Effects of shearing on shear modulus of clays. ASCE Journal of the Geotechnical Engineering Division, 102(GT9): 975-987.
Aylsworth, J.M., and Lawrence, D.E. 2003. Earthquake-induced landsliding east of Ottawa: A contribution to the Ottawa Valley Landslide Project, In: Proceedings of Geohazards 2003, 57-64.
Atkinson, G.M., and Beresnev, I.A. 1998. Compatible ground-motion time histories for new national seismic hazard maps. Canadian Journal of Civil Engineering, 25: 305-318.
Bazant, Z.P. and Planas, J. 1998. Fracture and size effect in concrete and other quasibrittle materials. CRC Press LLC, Boca Raton, FL.
Bernander, S. 2000. Progressive failure in long natural slopes: formation, potential extension and configuration of finished slides in strain-softening soils. Licentiate Thesis, Lulea University of Technology.
Bishop, A.W. 1971. The influence of progressive failure on the choice of the method of stability analysis. Geotechnique, 21: 168-172.
Bjerrum, L. 1955. Stability of natural slopes in quick clay. Geotechnique, 5: 101-119.
Desjardins, R. 1980. Tremblements de terre et glissements de terrain: Corrélation entre des datations au 14C et des données historiques à Shawinigan, Québec. Géographie physique Quaternaire, 34(3): 359-362.
Halchuk, S., and Adams, J. 2008. Fourth generation seismic hazard maps of Canada: maps and grid values to be used with the 2005 National Building Code of Canada. Geol. Survey of Can. Open File 5813, 32 pgs.
Hunter, J.A., Burns, R.A., Good, R.L., Aylsworth, J.M., Pullan, S.E., Perret, D., and Douma, M. 2007. Borehole shear wave velocity measurements of Champlain Sea sediments in the Ottawa-Montréal region. Geological Survey of Canada, Open File 5345.
Idriss, I.M. and Sun, J.I. 1991. Users’ Manual for SHAKE91: A modified version of SHAKE for conducting equivalent linear seismic response analyses of horizontally layered soil deposits. Centre for Geotechnical Modelling, University of California, Berkley.
Lebuis, J., J.-M. Robert, and P. Rissmann. 1983. Regional mapping of landslide hazard in Quebec. In: Symposium on Slopes on Soft Clays, Swedish Geotechnical Institute Report No. 17, Linkoping, 205-262.
Lefebvre, G., and La Rochelle, P. 1974. The analysis of two slope failures in cemented Champlain clays. Canadian Geotechnical Journal, 11: 89-108.
Leggett, R.F., and P. LaSalle. 1978. Soil studies at Shipshaw, Quebec: 1941 and 1969. Canadian Geotechnical Journal, 15: 556-564.
Locat, A., Leroueil, S., Bernander, S., Demers, S., Locat, J., and Ouehb, L. 2008. Study of a lateral spread failure in an eastern Canada clay deposit in relation with progressive failure: the Saint-Barnabé-Nord slide. GeoHazards 2008, Quebec City, 89-96.
National Research Council of Canada. 2005. National Building Code of Canada.
Natural Resources Canada. 2004. Urban Geology of the National Capital Area – Drift Thickness, Geological Survey of Canada, Terrain Sciences Division. Available http://gsc.nrcan.gc.ca/urbgeo/ [Accessed September 2006].
Ordonez, G. 2008. Users’ Manual for SHAKE2000: A computer program for the 1-D analysis of geotechnical earthquake engineering problems. June 2008 revision.
Palmer, A.C., and Rice, J.R. 1973. The growth of slip surfaces in the progressive failure of over-consolidated clay. Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 332(1591): 527-548.
Quinn, P.E. 2009. Large landslides in sensitive clay in eastern Canada and the associated hazard and risk to linear infrastructure. PhD Thesis, Department of Geological Sciences and Geological Engineering, Queen’s University.
Quinn, P.E., Diederichs, M.D., Hutchinson, D.J., and Rowe, R.K. 2007. An exploration of the mechanics of retrogressive landslides in sensitive clay. Proc. of the Canadian Geotechnical Conf, Ottawa, 721-727.
Quinn, P.E., Diederichs, M.D., Rowe, R.K. and Hutchinson, D.J. 2011. A new model for large landslides in sensitive clay using a fracture mechanics approach. Canadian Geotechnical Journal, 48 (8): 1151-1162.
Rice, J.R. 1968. A path independent integral and the approximate analysis of strain concentration by notches and cracks. Journal of Applied Mechanics, 35: 379-386.
Schnabel, P.B., Lysmer, J., and Seed, H.B. 1972. SHAKE: a computer program for earthquake response analysis of horizontally layered sites. UCB/EERC-72/12, Earthquake Engineering Research Center, University of California, Berkeley, 92 pages.
Seed, H.B. and Idriss, I.M. 1970. Soil moduli and damping factors for dynamic response analysis. Report No. UCB/EERC-70/10, University of California, Berkley.
Stark, T., and Contreras, I.A. 1996. Constant volume ring shear apparatus. Geotechnical Testing Journal, 19(1): 3-11.
Tavenas, F., Flon, P., Leroueil, S., and Lebuis, J. 1983. Remolding energy and risk of slide retrogression in sensitive clays. In: Symposium on Slopes on Soft Clays, Swedish Geotechnical Institute Report No. 17, Linkoping, 423-454.
Tremblay, A., B. Long, and Massé, M. 2003. Supracrustal faults of the St. Lawrence rift system, Québec: kinematics and geometry as revealed by field mapping and marine seismic reflection data. Tectonophysics, 369: 231-252.
Weibull, W. 1963. A theory of fatigue crack propagation in sheet specimens. Acta Metallurgica, 11: 745-752.