WATER RESOURCES RESEARCH, VOL. 48, W09524, doi:10.1029/2012WR011830, 2012

Analysis of rainfall-induced slope instability using a field of local
factor of safety
Ning Lu,1 Başak Şener-Kaya,1 Alexandra Wayllace,1 and Jonathan W. Godt2
Received 4 January 2012; revised 24 July 2012; accepted 25 July 2012; published 18 September 2012.

[1] Slope-stability analyses are mostly conducted by identifying or assuming a potential
failure surface and assessing the factor of safety (FS) of that surface. This approach ofassigning a single FS to a potentially unstable slope provides little insight on where thefailure initiates or the ultimate geometry and location of a landslide rupture surface. We
describe a method to quantify a scalar field of FS based on the concept of the Coulomb
stress and the shift in the state of stress toward failure that results from rainfall infiltration.
The FS at each point within a hillslope is called the local factor of safety (LFS) and isdefined as the ratio of the Coulomb stress at the current state of stress to the Coulomb stressof the potential failure state under the Mohr-Coulomb criterion. Comparative assessmentwith limit-equilibrium and hybrid finite element limit-equilibrium methods show thatthe proposed LFS is consistent with these approaches and yields additional insight into
the geometry and location of the potential failure surface and how instability may initiate
and evolve with changes in pore water conditions. Quantitative assessments applying
the new LFS field method to slopes under infiltration conditions demonstrate that the LFS
has the potential to overcome several major limitations in the classical FS methodologies
such as the shape of the failure surface and the inherent underestimation of slope
instability. Comparison with infinite-slope methods, including a recent extension tovariably saturated conditions, shows further enhancement in assessing shallow landslide
occurrence using the LFS methodology. Although we use only a linear elastic solution
for the state of stress with no post-failure analysis that require more sophisticatedelastoplastic or other theories, the LFS provides a new means to quantify the potential
instability zones in hillslopes under variably saturated conditions using stress-field based
methods.
Citation: Lu, N., B. Şener-Kaya, A. Wayllace, and J. W. Godt (2012), Analysis of rainfall-induced slope instability using a fieldof local factor of safety, Water Resour. Res., 48, W09524, doi:10.1029/2012WR011830.

1. Introduction[2] Slope-stability analyses are commonly performed using
limit-equilibrium methods because of their proven effectiveness and reliability [e.g., Fellenius, 1936; Bishop, 1955;
Morgenstern and Price, 1965; Janbu, 1973; Lade, 1992;
Michalowski, 2002; Duncan and Wright, 2005]. For suchmethods, the geometry and sometimes the location of thepotential failure surface are approximated or predetermined.
All two-dimensional limit-equilibrium methods consist of
discretizing the mass of a potential failure slope into smaller
vertical slices. All limit-equilibrium methods assume thateach individual slice is treated as a unique sliding block1
Department of Civil and Environmental Engineering, Colorado Schoolof Mines, Golden, Colorado, USA.
2
Geologic Hazards Science Center, U.S. Geological Survey, Denver,
Colorado, USA.

Corresponding author: N. Lu, Department of Civil and Environmental
Engineering, Colorado School of Mines, Golden, CO 80401, USA.
(ninglu@mines.edu)©2012. American Geophysical Union. All Rights Reserved.
0043-1397/12/2012WR011830(as shown in Figure 1a) and shear strength is mobilized withthe same ratio of shear stress to shear strength for all slices.
[3] By applying the principles of force and/or moment
equilibrium, all limit-equilibrium methods seek an indicatorof slope stability called the “factor of safety” or FS, which istypically defined as the ratio of the available shear strength tothe shear stress required for equilibrium along the prescribedfailure surface. Depending on how interslice forces are handled and the equilibrium principles that are involved, many
variations of limit-equilibrium methods have been developed.
For later comparison with the new method proposed in thiswork, the commonalities and differences among these classical methods are illustrated in Figure 1 and summarized
below.
[4] The Ordinary Method (Figure 1b) [Fellenius, 1936] isthe simplest and the oldest method of slices. It ignores all theinter-slice forces because of the assumption that the forces areparallel to the base of each slice. The method satisfies only
moment equilibrium and assumes a circular slip surface. Dueto its simplicity, it is possible to compute the factor of safety
by hand calculations. The ordinary method is generally lessaccurate compared to the other methods of slices. In somestudies [e.g., Duncan and Wright, 2005; Fredlund and Krahn,
1977; Abramson et al., 1996], it has been shown that theW09524

1 of 14

W09524

LU ET AL.: RAINFALL-INDUCED LANDSLIDE ANALYSIS

W09524

Figure 1. (a) Slice discretization of a potential sliding mass and forces acting on a typical slice used inlimit-equilibrium methods. Illustration of the commonalities and differences among the classical methodsof slices for (b) Ordinary (Fellenius) Method, (c) Bishop’s Simplified Method (vertical force equilibrium),(d) Janbu’s Simplified Method (horizontal force equilibrium), and (e) Morgenstern-Price Method.

accuracy is decreased for effective stress analyses and underhigh pore water pressure conditions.
[5] Bishops’s simplified method (Figure 1c) [Bishop,
1955] satisfies vertical force equilibrium for each slice andoverall moment equilibrium at the center of the circular trial
surface (Figure 1a). The method also ignores inter-slice shear
forces (X). Bishop’s simplified method is considered to be
more accurate than the Ordinary method, especially for
effective stress analyses with high pore water pressures [e.g.,
Duncan and Wright, 2005; Abramson et al., 1996].
[6] The Janbu Simplified Method (Figure 1d) satisfies
only overall horizontal force equilibrium and also assumes
that there are no interslice shear forces (X). The assumptionthat all the interslice forces are horizontal often leads tosmaller factor of safety when compared to more rigorous
methods that satisfy complete force and moment equilibrium.
To account for the effect of the shear forces, a correction
factor related to cohesion, angle of friction and the shape ofthe failure surface were presented by Janbu [e.g., Duncan andWright, 2005; Abramson et al., 1996].
[7] The Morgenstern and Price Method (Figure 1e, 1965)
considers both shear (X) and normal (E) interslice forces andsatisfies both moment and force equilibrium. The method
assumes an arbitrary mathematical function to describe thedirection of the interslice forces [e.g., Fredlund and Krahn,1977].
[8] Some recent developments of methods to calculate thefactor of safety include more accurate computation of theinterslice stress distributions [e.g., Duncan, 1996; Yu et al.,
1998; Swan and Seo, 1999] and numerical algorithms forshear strength reduction analysis from finite element methods
for multidimensional analysis under a continuum mechanics
framework [e.g., Matsui and San, 1992; Smith and Griffiths,
2004; Krahn, 2003]. The classical limit-equilibrium methods
described above and a hybrid FEM-limit equilibrium method[Krahn, 2003; GEO-SLOPE International, 2007] represent

2 of 14

W09524

LU ET AL.: RAINFALL-INDUCED LANDSLIDE ANALYSIS

W09524

Figure 2. Illustration of the concept of a scalar field of local factor of safety. The current state of stressat any point in hillslope (shown as the solid-line circle) has a shear stress t and shear strength t*. The
quantity t*/t can be used as an indicator of how far the current state of stress is from failure and isdefined as the local factor of safety (LFS). If suction stress ss is reduced at this point, the circle shifts
leftward and the LFS is reduced, indicating that the state of stress is closer to the failure condition or
LFS = 1.0.

the state-of-the-art in slope stability analysis and will be used
in what follows for comparison with the new LFS method.
[9] Because limit-equilibrium and stress-field based methods typically seek a single stability indicator for the entire
slope, they provide little insight into the location of initial
failure and how the region of potential instability evolves withchanges in pore water conditions that accompany rainfall
infiltration. We describe a method to calculate a scalar field offactor of safety to overcome some of these limitations andtake full advantage of modern numerical solutions for variably saturated flow and the state of stress in hillslopes.

2. Definition of the Scalar Field (or Local)of Factor of Safety[10] The basic idea is to seek a scalar quantity at each point
in a hillslope that is indicative of stability at that point. Our
method extends the approach described by Iverson and Reid
[1992] to cohesive, variably saturated porous material anduniquely provides a physically meaningful quantity to assess
the stability of partially saturated hillslopes. As illustrated in
two-dimensional stress space in Figure 2, if the current stateof stress at a point of interest is the solid circle, the stressquantity t indicates of the state of failure. This stress t can be
uniquely defined by shifting the Mohr circle leftward until it
touches the Mohr-Coulomb (M-C) failure envelope. Kinget al. [1994] called the shear stress on the Mohr circle along
the direction perpendicular to the M-C failure criterion“Coulomb stress.” Scalar fields of Coulomb stress calculated
in this manner near faults have been used to understand thestate of stress in the Earth’s crust and the occurrence ofearthquake aftershocks [e.g., Toda et al., 1998]. The pathsof changes in effective stress in slopes that result fromchanges in pore water conditions follow a similar leftward
shift. The size of the Mohr circle remains nearly constantunder changing pore water conditions because the size of
the Mohr circle is determined by the difference in principal
total stresses. The principal total stresses are mainly affectedby the geometry of the slope and self-weight of slope
materials. The former does not change for a given slope
over timescales relevant to rainfall infiltration and changes
in the latter are relatively small compared with changes inpore water pressure and effective stress during rainfall
infiltration and seepage that lead to a leftward shift of theMohr circle.
[11] The pattern of the shift in the stress path implies thatthe direction of failure under variably saturated soil conditions would be similar to the direction of Coulomb stress. If
the reference point for soil strength at the current state ofstress (the solid Mohr circle) can be estimated by the interceptof the extension of the Coulomb stress with the M-C envelope, shown as point B in Figure 2, a local factor of safety(LFS) at each point within the entire slope can be defined as
LFS ¼

t*
t

ð1Þ

where t* is the potential Coulomb stress and t is the current
state of Coulomb stress. By the similarity of triangles ACD
and ABE, the LFS can be expressed conveniently in terms of
the ratio of the adjusted mean stress of the current state ofstress to the adjusted mean stress of the potential failure state
under the Mohr-Coulomb criterion. The adjustment herecomes from the fact that hillslope materials are typicallycohesive and tensile stresses could occur in the area near theslope crest. To include such a situation, a stress value of AD,
in lieu of the stress invariant OD (Figure 2) is used. Similarly,
a stress value of AE in lieu of the stress invariant OE is used.
Graphically, the LFS defined by equation (1) is equal to

3 of 14LU ET AL.: RAINFALL-INDUCED LANDSLIDE ANALYSIS

W09524

W09524

where ua is the prevailing air pressure. Suction stress ss is acharacteristic scalar function of a soil and varies with soil
suction or saturation. It is equal to the pore water pressure
when soil is fully saturated but generally smaller than thepore water pressure in partially saturated soil. Suction stressfollows different patterns in clay, silt, and sand but can be
unified in a single closed form for all soils [Lu and Likos,
2004, Lu and Griffiths, 2004, Lu et al., 2010b] using thesame set of parameters to describe soil-water characteristic
curves in van Genuchten’s [1980] model, i.e.,s s ¼ ðua  uw Þ

ss ¼ 

ua  uw ≤ 0 ð5aÞðua  uw Þ
ð1 þ ½aðua  uw Þn Þ

ua  uw ≥ 0 ð5bÞðn1Þ=nwhere a and n are unsaturated soil parameters identical tothose in van Genuchten’s [1980] soil-water characteristic
curve model that relates the equivalent moisture content qe(equals to the volumetric water content q normalized by theresidual water content qr) to soil suction (ua  uw):
qe ¼

Figure 3. Illustration of the time-matching scheme inthe one-way coupled hydro-mechanical framework in
HILLSLOPE (FS)2.
AE over AD. If the Mohr-Coulomb failure criterion is used,
it can be shown that
LFS ¼

2 cos f′ðc′ þ p′ tan f′Þq′ð2Þ

where c′ is the drained cohesion of the slope material, f′ is thedrained friction angle of the slope material, p′ and q′ are themean and deviator effective stresses in two-dimensional
space, respectively, and are defined byp′ ¼

s 1 ′ þ s3 ′ s1 þ s3
 ss
¼
2
2

q′ ¼ s1 ′  s3 ′ ¼ s1  s3

ð3aÞs1 ′ ¼ s1  ua  s s

ð4aÞ

s3 ′ ¼ s3  ua  s s

ð4bÞ



1
1 þ ½aðua  uw Þn11=nð6Þ[12] Thus, in slopes the LFS so defined by equations (1)–(5) is a scalar quantity that varies in both space and time. As
illustrated in Figure 2, if the state of stress in a slope represented by the Mohr circle is below the material’s failureenvelope, the LFS is greater than unity. Whereas for thosezones in a slope with LFSs equal to unity, the states of stresshave reached their limit states and local failure occurs.
Therefore, LFS has the potential to delineate zones of stability or failure in slopes.

3. A Coupled Hydromechanical Framework[13] To compute the fields of stresses, soil suction andsaturation, and the LFS in a two-dimensional slope, a transient finite element framework was used. The code simulates
one-way coupled variably saturated flow and stress problems using finite elements (FEM) for hillslope flow andstress and factor of safety (HILLSLOPE FS2). The code wasdeveloped and validated by the authors [Lu et al., 2010a]
and was built from two previously uncoupled FEM codes(Figure 3). Both of the component codes were tested andvalidated by the original developers. The solution for transient unsaturated flow (Richards’ equation) is based on the
U.S. Department of Agriculture model SWMS_2D [Simunek
et al., 1994]. The Richards equation is [e.g., Freeze andCherry, 1979; Lu and Likos, 2004]

ð3bÞ

with the generalized effective stresses for variably saturated
porous media defined by suction stress ss [Lu and Likos,
2004, 2006] asq  qr¼
1  qrr  KðhÞrH þ W ¼∂ q ðhÞ∂t

ð7Þ

where K(h) is hydraulic conductivity function (HCF), h ispressure or suction head, H is total head (H = h + z), W isflux from a source or to a sink, and q (h) is volumetricmoisture content function commonly called the soil-water
retention curve (SWRC). In general, K(h) and q (h) are
nonlinear functions of pressure head h. Van Genuchten’s

4 of 14LU ET AL.: RAINFALL-INDUCED LANDSLIDE ANALYSIS

W09524[1980] model shown in equation (6) is used for SWRC.
Mualem’s [1976] is used for HCF.
n
k ¼ ks1

1  ðajhjÞn1 ½1 þ ðajhjÞn n1

1

f1 þ ½ajhjn g2 2n
1

o2
ð8Þ

where ks is saturated hydraulic conductivity.
[14] Solutions of the transient fields of suction headh (or matric suction (ua  uw)/g), moist unit weight g, andmoisture content q are used as the coupling variables for thesubsequent computation of the fields of total stress, suction
stress, and LFS. The solutions for stress and displacement are
based on FEM2D [Reddy, 1985] and use the field of moist
unit weight as the gravity term. The total stresses in the slope
are mainly affected by the geometry of the slope and thehillslope material’s self-weight. The governing equations
describing the total stresses in linear-elastic materials arefrom balance of linear momentum:r  ðsÞ þ g b ¼ 0

ð9Þ

where s is a stress tensor with three independent stressvariables in two-dimensional space, g is the bulk unit weightof hillslope materials, and b is the unit vector of body forces
with two components.
[15] The fields of total stress and suction head are then
combined to compute a field of suction stress (Figure 3), andfinally the total stress field and suction stress field are used
to compute the field of LFS.
[16] Because plane-strain linear elasticity theory is implemented in FEM2D [Reddy, 1985], only two elastic properties,
namely the elastic modulus E and Poisson’s ratio v, areneeded for solutions of both stress and displacement. Two
unsaturated hydromechanical properties, namely a and n, arecommonly defined in equations (5), (6), or (8). Additional
modifications to the code include improved computational
capacity and efficiency of the equation solver, addition of an
automatic mesh generator, and graphical input and output.
[17] Comparisons of simulation results of suction headand stress fields with those obtained from a commercial
slope-stability software package GeoStudio [GEO-SLOPE
International, 2007] were made and close agreement
between the two codes were obtained. For all the hybrid
FEM-limit equilibrium analyses and the classical limitequilibrium methods described in the following sections,
SLOPE/W [GEO-SLOPE International, 2007] was also used.
[18] We point out that the hydromechanical framework
based on linear elasticity solves statically admissible fields ofstress and displacement without invoking failure or the
redistribution of postfailure stress or displacement. A statically admissible stress distribution is defined only as one thatsatisfies the equilibrium differential equations, in the interiorof the body, and the boundary conditions wherever any
boundary tractions are prescribed [Malvern, 1969]. Thus LFS
less than unity only provide an indication of hypothetical
failure or the stress state in comparison with the hypothetical
failure parameters, namely drained cohesion c′ and angle ofinternal friction f′. This simple computational strategy isused to illustrate the location of potential failure zones without involving more sophisticated statically indeterminate
elastoplasticity theories. For this reason, in figures showing

W09524

contours of the LFS, no distinction is made for computed
LFS less than unity.
[19] In what follows, we first compare the LFS withseveral limit-equilibrium methods for stability analysis ofslopes of various geometries under either dry or saturated andfully drained conditions. The potential applicability of the
LFS to slope-stability analysis under variably saturated conditions is then explored for an example slope under transient
rainfall conditions and followed with comparisons to somerecent developments of infinite slope-stability theory undertransient infiltration conditions.

4. Comparisons With the Classical Factor
of Safety Methodologies
4.1. FS and LFS in Vertical Cuts
[20] We first compute the stress distribution in a verticalcut of either dry or saturated, but fully drained soil, undergravity for four different stability numbers m (0.05, 0.10,
0.15, and 0.20). The stability number is defined as m =c′/g H with g being the unit weight of the slope material,and H the height of the slope. The internal friction angle f′ isassumed to be 30 in all cases here. A small m value indicatesthe slope is less stable compared to those with large m values.
Based on the fields of stresses and strength parameters of theslope, the field of LSF is computed by equation (2) andplotted in Figure 4 (left). It can be seen that there are three
weak or potentially unstable zones in a vertical cut slope from
Figure 4a; near the vertical face, beneath the toe, and far fieldfrom the slope at the same elevation as the toe on horizontal
ground. The low LFS zone behind the vertical face is characterized by the principal maximum stress in the verticaldirection, whereas the low LFS zones beneath the toe and far
field from the slope are characterized by the principal maximum stress in the horizontal direction. For all the m values
examined here except for m = 0.05 (Figure 4a), the LFS in thefar field area away from the slope is never below 1.0, whereas
the LFS in the other two zones could be below 1.0 and thearea expands as the m value decreases (Figure 4, left). The
most pronounced weak zone is behind the vertical face withan up-leftward trend approximately 60 (or p4 þ f′2 = 60 ) fromhorizontal as the m value decreases (Figure 4a).
[21] To better illustrate this trend, contours of LFS equal to
1.0 are plotted in the right column in Figure 4 (thick solidlines) together with other potential failure surfaces identified
by the four limit-equilibrium stability methods from the
SLOPE/W software. To broaden the comparison and assessment, potential failure surfaces predicted by the Culmannmethod [Culmann, 1875; Lohnes and Handy, 1968] (analytical solution for finite slopes with a plane failure surface),
and a global FS calculated using an FEM stress solution
combined with a method of slices and iterative-search procedure to identify the least stable failure surface (Global
method from SLOPE/W) are also plotted in Figure 4 (right).
[22] For a relatively strong vertical cut (m = 0.20 andshown in Figure 4d), all four limit-equilibrium methods(Ordinary, Bishop, Janbu, and Morgenstern-Price), as well as
the Culmann and Global methods, predict that the vertical cut
is stable or the FS is greater than 1.0, whereas the LSF
method delineates a zone of potential failure behind the vertical face and near the toe. The potential failure surfaces for

5 of 14

W09524

LU ET AL.: RAINFALL-INDUCED LANDSLIDE ANALYSIS

W09524

Figure 4. (left) Distribution of the LFS for different stability numbers m = c′/gH calculated using
HILLSLOPE (FS)2 for vertical cuts. (right) Comparison of LFS with FS from various limit-equilibrium
methods. The black zone in the left column indicates the region where factors of safety are less than
1.0 from the linear elastic analyses. (a) m = 0.05, (b) m = 0.10, (c) m = 0.15, and (d) m = 0.20 (g = 20
kN/m3, f′ = 30 ).
Janbu, Ordinary, and Morgenstern-Price are identical (thinsolid curves), and very close to that given by the global
FS in the lower portion of the cut, but different than thatby the Bishop method (long dashed curve) in the middle
portion of the cut. The potential failure surface predictedby the Culmann method is quite different than the rest of
the limit-equilibrium methods. It predicts failure of the cut(short dashed line with FS = 1.0). Note that the lower partof the LFS = 1.0 contour follows closely the potential
failure surfaces predicted by the other methods. Overall,
Bishop’s method yields the lowest FS value (1.054) and isclose to the FS given by the Ordinary method (1.06). Both
Janbu and Morgenstern’s methods yield similar FS values(1.27 versus 1.26), much higher than the FS by Bishop andOrdinary, yet much smaller than the global FS method (1.59).
[23] The fundamental difference between the LFS = 1.0
and the failure surfaces calculated by the different limitequilibrium methods is that all points along a LFS = 1.0

contour are at the failure state, whereas the state of stress at
any point along a FS = 1.0 failure surface from the limitequilibrium methods may not reach the material strength atthat point.
[24] For vertical cuts with a stability number m = 0.15 all
the limit-equilibrium methods predict the same potential
failure surface with nearly identical FS values (1.05 – 1.06).
Compared to the case where m = 0.20, the Global FS method
yields much lower FS value (1.32) and the potential failuresurface is farther away from that given by the four limitequilibrium methods (Figure 4d). Again, the lower part of
the LFS = 1.0 contour follows closely the potential failuresurfaces predicted by the other methods. Furthermore,
the LFS method identifies a failure zone beneath the toe
where the LFS is less than unity. As the stability number mdecreases (shown respectively for m values of 0.10 and
0.05 in Figures 4b and 4a), all four limit-equilibrium methods
predict instability of the vertical cut with the identical FS

6 of 14

W09524

LU ET AL.: RAINFALL-INDUCED LANDSLIDE ANALYSIS

W09524

Figure 5. (left) Distribution of the LFS for different stability numbers m = c′/gH calculated using
HILLSLOPE (FS)2 for 60 slopes. (right) Comparison of LFS with FS from various limit-equilibrium
methods. The black zone in the left column indicates the region where factors of safety are less than 1.0 fromthe linear elastic analyses. (a) m = 0.05, (b) m = 0.10, (c) m = 0.15, and (d) m = 0.20 (g = 20 kN/m3, f′ = 30 ).
value (0.85 for m = 0.10 and 0.65 for m = 0.05) and similar
failure surface. For limit-equilibrium analysis, zones where
the FS is less than unity can also be identified and are shownas shaded areas in Figures 4a and 4b. Once again the lower
part of the LFS = 1.0 contour follows very closely with thepotential failure surfaces predicted by the limit-equilibrium
methods. Compared to the potential failure surfaces withhigher m values (Figures 4c), both the Global FS and Culmannmethods predict nearly identical failure surfaces but awayfrom those by either the limit-equilibrium methods or the
LFS method.
[25] In light of the above analysis for the stability ofvertical cuts it appears that the failure surface initiates fromthe toe region behind the vertical face and progresses up-andleft toward the top of the cuts. This trend can only be predicted by the new LFS methodology as the limit-equilibrium
methodology explicitly excludes detailed information on the
stress distribution and prescribes the failure surface along thegeneral direction of the Rankine’s active failure planes, i.e.,

f′
p
4þ 2

= 60 where the vertical is the maximum principal
stress direction and horizontal the minimum principal stressdirection.
4.2. FS and LFS in 60 Slopes[26] The distributions of LFS are shown in Figure 5 (left)
for an example 60 slope. All limit-equilibrium methods
produce similar (potential) failure surfaces except for thecase of stability number m = 0.05 (Figure 5a), where
differences near the top of the slope between the limitequilibrium and the Global FS methods appear. Only one
zone with LFS less than unity is identified for all m values.
This zone initiates near but above the toe and propagates
upward throughout the entire face of the slope as the stability
number m decreases. For the stability number m = 0.20
(Figure 5d), a zone of LFS less than unity coincides with thepotential failure plane predicted by the Culmann method but
is located above the potential failure surface predicted by all
the limit-equilibrium methods. For the stability number m =7 of 14

W09524

LU ET AL.: RAINFALL-INDUCED LANDSLIDE ANALYSIS

W09524

Figure 6. (left) Distribution of the LFS for different stability numbers m = c′/gH calculated using
HILLSLOPE (FS)2 for 30 slopes. (right) Comparison of LFS with FS by various limit-equilibrium
methods. The black zone in the left column indicates the region where factors of safety are less than 1.0 fromthe linear elastic analyses. (a) m = 0.05, (b) m = 0.10, (c) m = 0.15, and (d) m = 0.20 (g = 20 kN/m3, f′ = 30 ).
0.20 (Figure 5d), m = 0.15 (Figure 5c), and m = 0.10
(Figure 5b), all the limit-equilibrium methods predict that theslope is stable with similar FS among the four methods.
The FS (1.97) is slightly higher for the Morgenstern-Pricemethod for the stability number m = 0.20 case and is closer tothat given by the Global FS method (FS = 2.02). For the stability number m = 0.05 (Figure 5a), the FS predicted by all thelimit-equilibrium methods are similar and predict that theslope is unstable (FS < 1.0). However, the global FS method
predicts stability (FS = 1.03). The contour of unity for the LFS
follows the potential failure surfaces given by the limitequilibrium methods near the toe for the stability numbers
m = 0.15 and m = 0.10, but is somewhat different for thestability number m = 0.05. The contour of unity for the LFS is
substantially different compared to the failure surfaces predicted by the limit-equilibrium methods near the crest of theslope. The LFS unity contour exits the domain on the slope
face whereas the equilibrium methods predict that the failuresurface will intercept the ground surface behind the slope

crest. The LFS method predicts failure under conditions
where the other methods predict stability providing additional insight into the shape of the failure surface and sensitivity to the distribution of stress. The LFS method predicts
small unstable zones on the slope face for the m = 0.15 andm = 0.10 cases, whereas the limit-equilibrium methods indicate the slope is stable. This implies that the LFS provides anincrease in precision and sensitivity in predicting failure
geometry over that of the limit-equilibrium methods.
4.3. FS and LFS in 30 Slopes[27] The trend of larger radii of curvature of the failuresurface with decreasing slope gradient predicted by the LFS
methodology is examined further in the example showing
results for slopes inclined at 30 (Figure 6). The LFS
methodology predicts a low FS zone of shallow depth nearlyparallel to the slope surface for the slopes with stability
numbers of m = 0.10 and m = 0.15 (Figures 6b and 6c). As
the stability number decreases (from Figures 6d to 6a),8 of 14

W09524

LU ET AL.: RAINFALL-INDUCED LANDSLIDE ANALYSIS

this zone progressively elongates and eventually a quasitranslational zone with LFS < 1.0 that encompasses nearly
the entire slope surface develops. This pattern of failure issimilarly predicted by Savage’s [1994] exact solutions forstresses in finite slopes. In contrast, none of the limitequilibrium methods predict this shape of failure surface
or the potential instability for a slope with a stability
number m = 0.05. The smallest FS calculated by the limitequilibrium methods is 1.64 (Janbu) and all potential
failure surfaces are curved, in contrast to the nearlytranslational surface predicted by the LFS method. Shallow
landslides, typically translational slope failures a few
meters thick of unlithified soil mantle or regolith, may dominate mass-movement processes in hillslope environments[Cruden and Varnes, 1996; Trustrum et al., 1999; Sidle andOchiai, 2006]. We show that the limit equilibrium-methods
do not adequately describe the shape of such failures,
whereas the LFS methodology can predict the initiation,
evolution, and geometry of such slope failures.
[28] Slope-stability analysis of translational landslides hascommonly relied on the infinite-slope model [e.g., Wu andSidle, 1995; Dietrich et al., 1995; Baum et al., 2010],
which requires several strict assumptions, including slopeparallel failure surfaces, saturated materials, a sharp contrast
between the properties of overlying materials and substrate,and simultaneous failure at given depths. The assumptionsof saturated materials, sharp contrast between layered materials, and fixed failure depth have been removed using variably saturated theory for effective stress analysis [Lu andGodt, 2008]. Comparison of the LFS methodology with theinfinite-slope stability models is presented in the next section
for slopes under infiltration conditions.

5. Infiltration-Induced Slope Stability Analysis
Using LFS
5.1. Hydrologic Response[29] Here we use a 30 slope composed of homogeneous
silty soil as an example to illustrate the effects of rainfall
infiltration on slope stability. The domain and mechanicaland hydrologic properties of the hillslope materials areshown in Figure 7a. The surface infiltration boundary is
imposed using a 5-month rainfall record for September 2005
to January 2006 collected in the Seattle, Washington, area
[Godt et al., 2009] and shown in Figure 7b. This rainfall
record is representative of winter season precipitation forcoastal regions of the northwestern United States. Soil water
retention and suction stress characteristic curves of the silty
hillslope are shown in Figure 7c. The water table is set at 5 m
below the toe and forms the lower boundary. An initialhydrostatic condition above the water table is imposed,
which results in linear variation of soil suction with distanceabove the water table (not shown) and a nonlinear variationin volumetric soil moisture above the water table (shown in
Figure 8a). Both left and right sides of the domain are
assumed to be no-flow boundaries. The water table location
at the bottom of the domain is fixed.
[30] The variation in soil moisture content resulting fromthe 5-month rainfall is shown in Figure 8 and Figure 9a.
Because the soil near the slope surface is initially quite dry(q < 20%) under hydrostatic conditions and the prescribed
flux is nearly the same as the saturated hydraulic conductivity, a nearly saturated zone preceded by a sharp wetting front

W09524can be observed during the 5-month period (Figures 8b–8d).
The progress of the wetting front is also illustrated in
Figure 9a as profiles of moisture content at the middle of theslope (A-location in Figure 7a) at different times from the
beginning of the rainfall. The wetting front becomes morediffuse as time elapses. Soil suction in the slope changes in asimilar pattern to that of the moisture content, as illustrated by
the suction profiles at different times in Figure 9b.
[31] The 5-month rainfall episode results in 480 mm water
infiltrating into the slope and causes moisture content
redistribution throughout the entire slope. The infiltrated
water reaches the water table below the toe area in about
112 days, as shown in Figure 8f.
5.2. Mechanical Response[32] The dynamics of the mechanical response of theslope to the 5-month rainfall period can be quantitatively
investigated by examining the suction stress field. Under theframework of the unified effective-stress principle defined inequations (4)–(5), the field of suction stress in the slope canbe determined once either the moisture content or soil
suction field is known. Figure 10 shows the suction stressfield in the slope at different times. Before the rainfall
commences, as shown in Figure 9c, suction stress varies
nearly linearly with depth below the ground surface andhas much smaller absolute values than that of soil suction(shown in Figure 9b) at the same location. For example,
near the slope surface at the location A, soil suction is98.0 kPa (Figure 9b) whereas the corresponding suction
stress is 34.2 kPa (see Figure 9c). After the rainfall begins,suction stress quickly increases behind the wetting front asshown in Figure 9c and Figures 10b–10f. At the end of
December (90 days), a zone with suction stress >10 kPa
is about 0.5 m, as shown in Figure 9c and Figure 10e. After
that, suction stress continues to increase in near the slope
surface, (Figure 9c and Figures 10f). At the end of the
5-month simulation period (112 days), a zone with suction
stress >10 kPa is about 1.1 m (Figure 9c).
[33] The dynamics of slope stability now can be quantified
using the LFS in light of varying field of effective stress in theslope. The field of LFS is calculated simply by using
equations (2)–(4) with the field of total stress under thegravity computed from the FEM simulations. The fields of
LFS for the slope at different simulation times are shown in
Figure 11. Before the rainfall episode, the slope is stable(Figure 11a) and LFS near the slope surface are between 1.50
and 2.00. At 10 days after the rainfall begins, a thin zone
nearly parallel to the slope surface with LFS between 1.25
and 1.50 develops (Figure 9d and Figure 11b). At 90 days, athin zone nearly parallel to the slope surface with LFS lessthan 1.0 develops. This zone can be identified near the upper
middle portion of the slope (Figure 11e). At 112 days, thiszone expands to about a maximum thickness of 20 cm near
the middle of the slope (Figure 11f and Figure 9d). This zoneof LFS < 1.0 spatially and temporally coincides with the large
increase in suction stress shown in Figures 10e–10f. Within
this zone, as a result of increases in suction stress ss ordecreases in effective stress (equation (4)), the states of stressshift leftward in in the Mohr diagram (as shown in Figure 2
moving from the solid Mohr circle to the dashed Mohr
circle) leading to failure of the soil. It is noteworthy that while
the geometry of the moisture content and suction stress contours are more or less parallel to the slope surface (Figures 8

9 of 14

W09524

LU ET AL.: RAINFALL-INDUCED LANDSLIDE ANALYSIS

W09524

Figure 7. (a) The geometry, material properties and FEM mesh, (b) rainfall scenario from Seattle area
recorded during 2005–2006 [Godt et al., 2009], and (c) soil water retention curve (SWRC) and suction
stress characteristic curve (SSCC) of the silty soil used in the simulations. Values for residual volumetricwater content qr, saturated volumetric water content qs, van Genuchten [1980] parameters, a, n, saturated
hydraulic conductivity, Ks, soil cohesion c′, moist unit weigh g, angle of internal friction f′, Poisson’s
ratio n, and Young’s modulus E, are shown in the inset table in Figure 7a.
and 10), the contours of LFS near the slope surface have alens shape (Figures 11e-11f) that is consistent with fieldobservations of shallow landslides [e.g., Cruden and Varnes
1996]. The LFS, as illustrated here, provides quantitative
information regarding the initial failure location and apotentially rigorous way to compute the timing, geometry,and progression of potential instability driven by changes in
the pore water conditions.
[34] The LFS also provides a quantitative way to examine
the suitability of the widely used infinite-slope model for
shallow landslides. In recent work by Lu and Godt [2008],
the infinite-slope stability model was extended to accommodate variably saturated conditions by using the unified

10 of 14

LU ET AL.: RAINFALL-INDUCED LANDSLIDE ANALYSIS

W09524

W09524

Figure 8. Contours of simulated moisture content at different times for a silty slope under rainfall
conditions.
effective stress equation (4) in lieu of Terzaghi’s effective
stress equation for saturated soil. The FS is [Lu and Godt,
2008]
Fðz; tÞ ¼

tan f′
2c′
ss ðz; tÞþð tan b þ cot bÞ tan f′ ð10Þtan b g z sin 2 bgz

where b is the slope angle, g is the moist unit weight of soil
within the distance z from the slope surface, and t is time.
Equation (10) is used to calculate profiles of factor of safety(FS) at location A shown in Figure 7a, using the suction
stress profiles at different times shown in Figure 9c and thesame material strength properties. The calculated FS profiles
are shown in Figure 9e and Figure 9f (same with Figure 9e
but with different scales in FS values). The FS for thehydrostatic conditions (t = 0) calculated from the infiniteslope stability model increases monotonically from about
2.0 at 3 m below the slope surface to more than 20 near theslope surface (Figure 9e). The monotonic variation pattern isalso obtained for cases in which the soil is either dry or
saturated in the infinite-slope models. Under unsaturated
seepage conditions, the FS profiles are no longer monotonic
and the location of minimum FS can occur within unsaturated zone, as illustrated in Figure 9f. Lu and Godt [2008]
and Godt et al. [2009] show that FS could be less than 1.0
in sandy soil on steep slopes under heavy rainfall conditions.
For the silty slope after the 5-month rainfall episode, FS canbe reduced to less than 3.0 but never reaching 1.0, implyingthat this slope is stable if infinite-slope theory were used toassess stability. Comparing Figure 9d and Figure 9f, both

LFS and the infinite-slope models produce about the samevalues of FS near 2.0 at relative large distances from theslope surface, where z > 2.5 m. However, using the effective
stress calculation by FEM (in Figure 9d and Figure 11), the
LFS is shown to be less than 1.0 near the slope surface fora certain time period. The overestimation of stability near
the slope surface by the infinite slope model compared to the
LFS method results from the appearance of distance z in
the denominator in equation (10) and is overcome by the
new LFS method.

6. Summary and Conclusions[35] Slope-stability analysis is typically performed by
identifying or assuming the potential failure surface andassessing the factor of safety (FS) of that surface. The
approach of calculating a single FS for a slope does notprovide vital information on how the failure surface initiatesand evolves with changes in pore water conditions. A method
to calculate a scalar field of FS is proposed based on theconcept of Coulomb stress and the shift of stress paths toward
the failure state in slopes under variably saturated infiltrationand seepage. The FS at each point within a slope is called
local factor of safety (LFS) and is defined as the ratio of
the adjusted mean stress of the current state of stress to theadjusted mean stress of the potential failure state under the
Mohr-Coulomb failure criterion. This LFS includes bothcompressive and extensional stress failure regimes.
[36] Comparative assessments with the hybrid FEM-limit
equilibrium slope stability analyses and stress field-based

11 of 14

W09524

LU ET AL.: RAINFALL-INDUCED LANDSLIDE ANALYSIS

W09524

Figure 9. Profiles of (a) moisture content, (b) suction, (c) suction stress, (d) local factor of safety, (e) and(f) factor of safety calculated using an infinite-slope stability analysis in the middle of the silty slope at
different times.
FEM-LFS method show that the proposed LFS is consistent
with the classic methodologies and yields further insights
into the geometry of the potential failure surface and how
the failure surface initiates and evolves. From sensitivity

analysis of the stability of slopes with various inclined
angles and materials, it is shown that the limit-equilibrium
methods provide an estimate of stability that exceeds thatfrom the LFS method by a factor up to 2.39 as shown in

12 of 14

LU ET AL.: RAINFALL-INDUCED LANDSLIDE ANALYSIS

W09524

Figure 10.

Figure 11.

Contours of simulated suction stress at different times for a silty slope under rainfall conditions.

Contours of simulated LFS at different times for a silty slope under rainfall conditions.
13 of 14

W09524

W09524

LU ET AL.: RAINFALL-INDUCED LANDSLIDE ANALYSIS

Figure 6a. That is, local failure could occur in a slope even if
the global FS for the slope is calculated to be 2.39. For steep
slopes, while the limit-equilibrium methods produce similargeometries of failure surface near the toe region, they do not
capture the geometry of failure surfaces in the middle part ofthe slope and at higher regions near the crest. For gentle
slopes, the limit-equilibrium methods fail to capture both thegeometry of the failure surfaces and their occurrence in thenear-surface environments.
[37] Quantitative assessments applying the new LFS field
method to slopes under infiltration conditions demonstratethat it has a potential to overcome several major limitationsin the commonly used global FS methodologies, including
accurate prediction of the rupture surface geometry andtheir occurrence near the slope surface. Comparison with theinfinite-slope methods including some recent expansion forvariably saturated slope conditions shows further enhancementin predicting shallow landslides by the LFS methodology.
Although this work only employs linear elasticity with nopostfailure elastoplastic analysis that require more sophisticated elastoplasticity theories, the LFS can provide newperspectives to the stability of slopes for stress field-based
methodologies.
[38] Acknowledgments. This research is partially supported by a grantfrom the National Science Foundation (NSF CMMI-0855783) to NL and JWG
and a grant from the U.S. Geological Survey (USGS G09AC00085) to NL.

References
Abramson, L. W., T. S. Lee, S. Sharma, and G. M. Boyce (1996), Slope
Stability and Stabilization Methods, 629 pp., John Wiley, Hoboken, N. J.
Baum, R. L., J. W. Godt, and W. Z. Savage (2010), Estimating the timing
and location of shallow rainfall-induced landslides using a model fortransient, unsaturated infiltration, J. Geophys. Res., 115, F03013,
doi:10.1029/2009JF001321.
Bishop, A. W. (1955), The use of slip circle in the stability analysis ofslopes, Geotechnique, 5(1), 7–17, doi:10.1680/geot.1955.5.1.7.
Cruden, D. M., and D. J. Varnes (1996), Landslide types and processes, inLandslides: Investigation and Mitigation, Transportation Research
Board Special Report, vol. 247, edited by A. K. Turner and R. L.
Schuster, pp. 36–75, Natl. Acad. Press, Washington, D.C.
Culmann, C. (1875), Die Graphische Statik (Graphic Statics), 644 pp.,
Meyer and Zeller, Zurich.
Dietrich, W. E., R. Reiss, M. L. Hsu, and D. R. Montgomery (1995),
A process-based model for colluvial soil depth and shallow landsliding
using digital elevation data, Hydrol. Processes, 9(3–4), 383–400, doi:10.
1002/hyp.3360090311.
Duncan, J. M. (1996), State of the art: Limit equilibrium and finite-element
analysis of slopes, J. Geotech. Geoenviron. Eng., 122(7), 577–596.
Duncan, J. M., and S. G. Wright (2005), Soil Strength and Slope Stability,
309 pp., John Wiley, Hoboken, N. J.
Fellenius, W. (1936), Calculation of the stability of earth dams, inTransactions of the 2nd Congress on Large Dams, vol. 4, pp. 445–463,
Int. Comm. on Large Dams, Washington, D.C.
Fredlund, D. G., and J. Krahn (1977), Comparison of slope stability
methods of analysis, Can. Geotech. J., 14(3), 429–439, doi:10.1139/
t77-045.
Freeze, R. A., and J. A. Cherry (1979), Groundwater, 604 pp., Prentice Hall,Englewood Cliffs, New Jersey.
GEO-SLOPE International (2007), GeoStudio, Calgary, Alberta, Canada.
Godt, J. W., R. L. Baum, and N. Lu (2009), Landsliding in partially saturated materials, Geophys. Res. Lett., 36, L02403, doi:10.1029/
2008GL035996.
Iverson, R. M., and M. E. Reid (1992), Gravity-driven groundwater flow
and slope failure potential 1. Elastic effective-stress model, Water
Resour. Res., 28(3), 925–938, doi:10.1029/91WR02694.
Janbu, N. (1973), Slope stability computations. Embankment-dam engineering, in Casagrande, vol. R, edited by C. Hirschfeld and S. J. Poulos,
pp. 49–86, John Wiley, New York.

W09524

King, G. C. P., R. S. Stein, and J. Lin (1994), Static stress changes and thetriggering of earthquakes, Seismol. Soc. Am. Bull., 84, 935–953.
Krahn, J. (2003), The 2001 R.M. Hardy Lecture: the limits of limit equilibrium analyses, Can. Geotech. J., 40, 643–660, doi:10.1139/t03-024.
Lade, P. V. (1992), Static instability and liquefaction of loose fine sandy
slopes, J. Geotech. Eng., 118(1), 51–71, doi:10.1061/(ASCE)07339410(1992)118:1(51).
Lohnes, R. A., and R. L. Handy (1968), Slope angle in friable loess,J. Geol., 76(3), 247–258, doi:10.1086/627327.
Lu, N., and J. W. Godt (2008), Infinite slope stability under unsaturated
seepage conditions, Water Resour. Res., 44, W11404, doi:10.1029/
2008WR006976.
Lu, N., and D. V. Griffiths (2004), Profiles of steady-state suction stress inunsaturated soils, J. Geotech. Geoenviron. Eng., 130(10), 1063–1076,
doi:10.1061/(ASCE)1090-0241(2004)130:10(1063).
Lu, N., and W. J. Likos (2004), Unsaturated Soil Mechanics, 584 pp., John
Wiley, New York.
Lu, N., and W. J. Likos (2006), Suction stress characteristic curve for unsaturated soils, J. Geotech. Geoenviron. Eng., 132(2), 131–142, doi:10.1061/(ASCE)1090-0241(2006)132:2(131).
Lu, N., A. Wayllace, and J. W. Godt (2010a), A hydro-mechanical model
for predicting infiltration-induced landslides, paper presented at 2010
Modeling Conference, U. S. Geol. Soc., Bloomfield, Colo.
Lu, N., J. Godt, and D. T. Wu (2010b), A closed-form equation for effective
stress in unsaturated soil, Water Resour. Res., 46, W05515, doi:10.1029/
2009WR008646.
Malvern, L. E. (1969), Introduction to the Mechanics of a Continuous
Medium, Prentice-Hall, Upper Saddle River, N. J.
Matsui, T., and K. C. San (1992), Finite element slope stability analysis byshear strength reduction technique, Soil Found., 32(1), 59–70, doi:10.
3208/sandf1972.32.59.
Michalowski, R. L. (2002), Stability charts for uniform slopes, J. Geotech.
Geoenviron. Eng., 128(4), 351–355, doi:10.1061/(ASCE)1090-0241(2002)128:4(351).
Morgenstern, N. R., and V. E. Price (1965), The analysis of the stability ofgeneral slip surface, Geotechnique, 15(4), 289–290.
Mualem, Y. (1976), Hysteretical models for prediction of the hydraulicconductivity of unsaturated porous media, Water Resour. Res., 12(6),1248–1254, doi:10.1029/WR012i006p01248.
Reddy, J. N. (1985), Introduction to the Finite Element Method, McGrawHill, New York.
Savage, W. Z. (1994), Gravity-induced stress in finite slopes, Int. J. Mech.
Min. Sci. Geomech. Abstr., 31(5), 471–483, doi:10.1016/0148-9062(94)90150-3.
Sidle, R. C., and H. Ochiai (Eds.) (2006), Landslides: Processes, Prediction,and Land Use, in Water Resour. Monogr., vol. 18, 312 pp., AGU,
Washington D.C.
Simunek, J., T. Vogel, and M. T. van Genuchten (1994), The SWMS2D code for simulating water flow and solute transport in two-dimensionalvariably saturated media, Version 1.21, Res. Rep. 132, USDA-ARS
U.S. Salinity Lab., Riverside, Calif.
Smith, I. M., and D. V. Griffiths (2004), Programming the Finite Element,
4th ed., John Wiley, West Sussex, U. K.
Swan, C. C., and Y. K. Seo (1999), Limit state analysis of earthen slopes
using dual continuum/FEM approaches, Int. J. Numer. Anal. Methods
Geomech., 23(12), 1359–1371, doi:10.1002/(SICI)1096-9853(199910)
23:12<1359::AID-NAG39>3.0.CO;2-Y.
Toda, S., R. S. Stein, P. A. Reasenberg, J. H. Dieterich, and A. Yoshida(1998), Stress transferred by the 1995 Mw = 6.9 Kobe, Japan, shock:Effect on aftershocks and future earthquake probabilities, J. Geophys.
Res., 103(B10), 24,543–24,565, doi:10.1029/98JB00765.
Trustrum, N. A., B. Gomez, M. J. Page, L. M. Reid, and D. M. Hicks(1999), Sediment production, storage and output: The relative role oflarge magnitude events in steepland catchments, Z. Geomorphol.
Suppl., 115, 71–86.
van Genuchten, M. T. (1980), A closed form equation for predicting thehydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J., 44,
892–898, doi:10.2136/sssaj1980.03615995004400050002x.
Wu, W., and R. C. Sidle (1995), A distributed slope stability model forsteep forested basins, Water Resour. Res., 31(8), 2097–2110, doi:10.
1029/95WR01136.
Yu, H. S., R. Salgado, S. Sloan, and J. Kim (1998), Limit analysis versus
limit equilibrium for slope stability assessment, J. Geotech. Geoenviron.
Eng., 124(1), 1–11, doi:10.1061/(ASCE)1090-0241(1998)124:1(1).

14 of 14