otsotropicfNew

1
Copyright © 2024 Seequent Limited, The Bentley Subsurface Company. All rights reserved.

No part of this work may be reproduced or transmitted in any form or by any means, electronic or
mechanical, including coping, recording, or by any information storage or retrieval system, without the
prior permission of Seequent (GeoStudio).

Trademarks: GeoStudio, SLOPE/W, SEEP/W, SIGMA/W, QUAKE/W, CTRAN/W, TEMP/W, AIR/W,

SLOPE3D, SEEP3D, TEMP3D, AIR3D, CTRAN3D, and other trademarks are the property or registered
trademarks of their respective owners.

Email: sales@geoslope.com

ii
Contents
Symbols.......................................................................................................................................................vii
1 Introduction..........................................................................................................................................1
2 Theory...................................................................................................................................................3
2.1 Background ...................................................................................................................................3
2.2 Factor of Safety Definition ............................................................................................................6
2.3 Column Forces ..............................................................................................................................7
2.3.1 Base Normal Force..............................................................................................................10
2.3.2 Intercolumn Normal Force..................................................................................................11
2.3.3 Intercolumn Shear Force ....................................................................................................12
2.4 Limit Equilibrium Formulation ....................................................................................................13
2.4.1 Moment Equilibrium...........................................................................................................14
2.4.2 Force equilibrium................................................................................................................15
2.5 Solution Algorithm......................................................................................................................15
2.5.1 2D Solution Workflow.........................................................................................................16
2.5.2 3D Solution Workflow.........................................................................................................17
2.6 Convergence and the Role of Slip Surface Shape .......................................................................19
2.6.1 Circular Slip Surface ............................................................................................................19
2.6.2 Planar Slip Surface ..............................................................................................................20
2.6.3 Composite Slip Surface .......................................................................................................21
2.6.4 Block Specified ....................................................................................................................22
2.6.5 Symmetrical Slip Surface in 3D ...........................................................................................22
2.6.6 Summary.............................................................................................................................23
2.7 Stress/Force Distributions ..........................................................................................................24
2.8 Commentary on the LE Method .................................................................................................28
3 Slip Surface Searching Techniques .....................................................................................................29
3.1 Introduction ................................................................................................................................29
3.2 Grid and Radius...........................................................................................................................29
3.3 Entry and Exit..............................................................................................................................32
3.3.1 Circular Slip Surfaces...........................................................................................................32
3.3.2 Ellipsoidal Slip Surfaces.......................................................................................................34
3.4 Fully Specified .............................................................................................................................36
3.4.1 2D Definition.......................................................................................................................36

iii
 3.4.2 3D Definition.......................................................................................................................37
3.5 Block Specified............................................................................................................................38
3.6 Cuckoo Search ............................................................................................................................40
3.7 Non-Circular Slip Surfaces...........................................................................................................42
3.7.1 Impenetrable Material Model ............................................................................................42
3.7.2 Weak Surfaces ....................................................................................................................43
3.8 Anisotropic Surfaces ...................................................................................................................44
3.9 Optimization ...............................................................................................................................46
3.9.1 2D Optimization..................................................................................................................47
3.9.2 3D Optimization..................................................................................................................50
3.10 Effect of Strength on Critical Slip Surface Location ....................................................................52
3.10.1 Frictional vs Undrained Strength ........................................................................................52
3.10.2 Minimum Depth Control.....................................................................................................55
3.10.3 Best Practice .......................................................................................................................55
3.11 Tension Cracks and Exit Projections ...........................................................................................55
3.12 Results Interpretation.................................................................................................................56
3.13 Commentary on Physical Admissibility .......................................................................................56
3.14 Error Message for Slip Surfaces ..................................................................................................57
4 Material Models .................................................................................................................................59
4.1 Introduction ................................................................................................................................59
4.2 Mohr-Coulomb ...........................................................................................................................59
4.3 Spatial Mohr-Coulomb................................................................................................................60
4.4 Undrained Strength ....................................................................................................................61
4.5 Impenetrable (Bedrock)..............................................................................................................61
4.6 Shear-Normal Function...............................................................................................................61
4.7 Anisotropic Function...................................................................................................................62
4.8 Anisotropic Strength...................................................................................................................63
4.9 SHANSEP .....................................................................................................................................64
4.10 Hoek-Brown ................................................................................................................................65
4.11 Compound Strength ...................................................................................................................69
4.11.1 Dip and Dip Direction..........................................................................................................70
4.12 High Strength ..............................................................................................................................73
4.13 Bilinear........................................................................................................................................73

iv
 4.14 Strength as a Function of Depth .................................................................................................74
4.15 Combined Models: Frictional-Undrained ...................................................................................75
4.16 Unsaturated Shear Strength .......................................................................................................76
4.17 Soil Unit Weight ..........................................................................................................................76
5 External Loads.....................................................................................................................................77
5.1 Ponded Fluid ...............................................................................................................................77
5.2 Surcharge Loads..........................................................................................................................78
5.3 Point Loads .................................................................................................................................81
6 Pore Water Pressure...........................................................................................................................81
6.1 Introduction ................................................................................................................................81
6.2 Piezometric Surfaces...................................................................................................................81
6.3 Ru Coefficients............................................................................................................................84
6.4 B-bar Coefficients .......................................................................................................................84
6.5 Spatial Function ..........................................................................................................................85
6.6 Finite Element Pore Water Pressures .........................................................................................86
7 Reinforcement and Structural Components.......................................................................................88
7.1 Reinforcement Load Calculation.................................................................................................88
7.1.1 Face Anchorage ..................................................................................................................89
7.2 Geosynthetics .............................................................................................................................90
7.3 Anchors and Nails .......................................................................................................................91
7.3.1 Grouted Friction Anchor .....................................................................................................92
7.4 Pile Reinforcement .....................................................................................................................93
7.5 User-defined Reinforcement ......................................................................................................93
7.6 Concentrated vs Distributed.......................................................................................................95
7.7 Manufacturer Library..................................................................................................................96
7.8 Factor of Safety Dependency......................................................................................................97
7.9 Results Visualization ...................................................................................................................97
8 Stress-based Stability Analysis............................................................................................................97
8.1.1 Calculating the Stability Factor ...........................................................................................97
8.1.2 Establishing Finite Element Stresses...................................................................................98
8.1.3 Commentary on Stress-Based Method...............................................................................99
9 Rapid Drawdown Analysis Methods ...................................................................................................99
9.1 Simple Effective Strength Method............................................................................................100

v
 9.2 Rigorous Effective Strength Method ........................................................................................100
9.3 Staged Drawdown (Duncan et al., 1990) ..................................................................................101
10 Seismic and Dynamic Stability Analysis ........................................................................................103
10.1 Pseudo-static Analysis ..............................................................................................................103
10.2 Undrained Response during Dynamic Loading .........................................................................105
10.2.1 Effective Stress Strength...................................................................................................105
10.2.2 Undrained Strengths (Duncan, Wright, and Wong, 1990) ................................................106
11 Newmark Permanent Deformation ..............................................................................................107
12 Probabilistic Analysis ....................................................................................................................112
12.1 Probability Density Functions ...................................................................................................114
12.2 c – 𝜙' Correlation......................................................................................................................121
12.3 Probability of failure and reliability index.................................................................................122
12.4 Spatial variability ......................................................................................................................123
12.5 Number of Monte Carlo Trials ..................................................................................................128
12.6 Multiple statistical parameters.................................................................................................129
13 Sensitivity Analysis........................................................................................................................129
14 Partial Factors ...............................................................................................................................130
14.1 Favorable or unfavorable actions .............................................................................................130
14.2 Pore water pressures and ponded water .................................................................................132
14.3 Shear strength ..........................................................................................................................133
14.4 Reinforcement Loads................................................................................................................134
15 References ....................................................................................................................................134

vi
Symbols
𝛼 Sliding direction, radians 𝑐' Effective cohesion, kPa
𝛼 Apparent dip, degrees '
in the horizontal direction, 𝑐ℎ
𝛽 Area of column base, m2
𝛽 Strike of the structure, degrees intercept of R undrained strength, 𝑐𝑅
'
𝛽 Reliability index in the vertical direction, 𝑐𝑣
𝛾 Unit weight 𝑐𝑢 Maximum undrained strength, kPa
of the soil, 𝛾 𝑑 Scale of fluctuation
of water, 𝛾𝑤 𝐷 External point loads, N
Γ Dimensionless variance function 𝐷 Rock mass disturbance
𝛿 True dip, degrees 𝐷 Bound diameter, m
𝛿 Interface shear angle, degrees 𝑑𝐾 Intercept of the isotropic consolidation
𝑐=1
𝜀 Normal standard deviation strength envelop
𝜃 Volumetric content, m3/m3 𝑒 Void ratio
𝐸 Intercolumn normal force, N
water content, 𝜃𝑤
𝑓 Intercolumn force function
saturated water content, 𝜃𝑠 𝐹 Factor of safety
residual water content, 𝜃𝑟 for force equilibrium, 𝐹𝑓
𝜃 Offset of lognormal mean
for moment equilibrium, 𝐹𝑚
𝜆 Scaling factor for equilibrium
in the 𝑥 direction, 𝐹𝑥
in the 𝑥 direction, 𝜆𝑥
in the 𝑧 direction, 𝐹𝑧
in the 𝑧 direction, 𝜆𝑧
𝐹 Force, N
𝜎 Stress, kPa
effective, 𝜎' horizontal, 𝐹ℎ
major principal, 𝜎1 vertical, 𝐹𝑣
𝐹𝑃𝑅 Factored pullout resistance, kPa
intermediate principal, 𝜎2
𝐹𝑇𝐶 Factored tensile capacity, N/m
minor principal, 𝜎3 𝑔 Gravitational acceleration, m/s2
horizontal, 𝜎ℎ 𝐺 Specific gravity of the grains
normal, 𝜎𝑛 𝐺𝑆𝐼 Geological strength index
𝐻 Intercolumn horizontal shear force, N
vertical, 𝜎𝑣 𝐻𝑠 Soil column height, m
𝜎 Strength, kPa ℎ𝑤 Height of water, m
tensile, 𝜎𝑡 𝑘 Seismic coefficient
uniaxial compressive, 𝜎𝑐𝑖 horizontal, 𝑘ℎ
𝜎 Standard deviation of statistical distribution
vertical, 𝑘𝑣
𝜏 Shear stress, kPa
𝑘 Ellipsoid aspect ratio
𝜏 Reinforcement-soil interface shear
𝑘 Correlation coefficient
strength, kPa
𝑙 Reinforcement length, m
𝜙' Effective angle of internal friction, degree
𝐿𝑎
' Effective anchorage length, m
in the horizontal direction, 𝜙ℎ 𝑚 Over-consolidation non-linearity exponent
intercept of R undrained strength, 𝜙𝑅 𝑚 Number of variables
' 𝑀 Lognormal mean
in the vertical direction, 𝜙𝑣 𝑚𝑏, 𝑠, 𝑎
𝜓𝐾 Hoek-Brown material constants
Slope of the isotropic consolidation
𝑐=1 𝑛 Porosity
strength envelop
𝑛 Hybrid ellipsoid exponent
𝑎, 𝑏, 𝑐 Principal semi-axis
𝑁 Total normal force (base of column), N
𝑎 Acceleration, m/s2
𝑁 Normalized random number
average, 𝑎𝑎𝑣𝑒
for 𝜙 , 𝑁1
'

yield, 𝑎𝑦
for 𝑐', 𝑁2
𝐵̅ Skempton’s B-bar coefficient
𝑐 Desired level of confidence adjusted for cohesion, 𝑁𝑎

vii
𝑁𝑚𝑐 Number of Monte Carlo trials
𝑂𝐶𝑅 Over-consolidation ratio
𝑃𝐹 Pullout force, N
𝑃𝑅 Pullout resistance, kPa
𝑅 Moment arm, m
𝑅𝑢 Pore water pressure ratio
𝑟 𝑥, 𝑟 𝑦 , 𝑟 𝑧 Ellipsoid radii

𝑅𝐹 Reduction factor
for creep, 𝑅𝐹𝐶𝑅
for chemical/environmental effects, 𝑅𝐹𝐶𝐻
for the extrapolation of data, 𝑅𝐹𝐷
for installation damage, 𝑅𝐹𝐼𝐷
overall, 𝑅𝐹𝑜
for weathering, 𝑅𝐹𝑊
𝑅𝑅𝐹 Resistance reduction factor
𝑠 Shear strength, kPa
𝑆 Shear force mobilized (base of column), N
𝑆 Degree of saturation
𝑆 Out-of-plane spacing, m
𝑆 Lognormal standard deviation
𝑆𝑒 Effective saturation
𝑆𝐼𝐴 Interface adhesion
𝑠𝑚 Mobilized shear forces, N
𝑆𝑁𝐶 Ratio of undrained strength
𝑠𝑟 Resisting shear forces, N
𝑠𝑢 Undrained shear strength, kPa
𝑆𝐴𝐹 Surface area factor
𝑆.𝐹. Stability factor
𝑡 Time, s
𝑇 Reinforcement force, N
𝑇 Strength, kPa
embedded end, 𝑇𝑒
reinforcement-facing connection, 𝑇𝑓
pullout resistance, 𝑇𝑠
reinforcement tensile, 𝑇𝑡
𝑇𝐶 Tensile capacity, N
𝑢 Pressure, kPa
of pore water, 𝑢𝑤
of pore air, 𝑢𝑎
𝑢 Mean of statistical distribution
𝑥', 𝑦', 𝑧' Local coordinates
𝑥0,𝑦0, 𝑧0 Coordinates of ellipsoid center
𝑉 Velocity, m/s
𝑊 Column weight, N
𝑋 Intercolumn shear force, N
𝑍 Distance between two sections, m

viii
1 Introduction
Analysing the stability of earth structures is the oldest type of numerical analysis in geotechnical
engineering. The idea of discretizing a potential sliding mass into columns was introduced early in the
20th Century. In 1916, Petterson (1955) presented the stability analysis of the Stigberg Quay in
Gothenberg, Sweden where the slip surface was taken to be circular, and the sliding mass was divided
into columns. During the next few decades, Fellenius (1936) introduced the Ordinary or Swedish method
of columns. In the mid-1950s Janbu (1954) and Bishop (1955) developed advances in the method. The
advent of computers in the 1960’s made it possible to handle the iterative procedures inherent to the
method, which led to mathematically rigorous formulations such as those developed by Morgenstern
and Price (1965) and by Spencer (1967). One of the reasons the limit equilibrium method was adopted
so readily, is that solutions could be obtained by hand-calculations. Simplifying assumptions had to be
adopted to obtain solutions, but the concept of numerically dividing a larger body into smaller pieces for
analysis purposes was rather novel at the time.

Modern limit equilibrium software has made it possible to handle increased complexity within an
analysis. It is now possible to deal with complex stratigraphy, highly irregular pore water pressure
conditions, various linear and nonlinear shear strength models, almost any kind of slip surface shape,
concentrated loads, structural reinforcement, and full three-dimensional (3D) analysis. Limit equilibrium
formulations based on the method of columns are also being applied increasingly to the stability analysis
of structures such as tie-back walls, slopes reinforced with soil nails or geosynthetics, and even the
sliding stability of structures subjected to high horizontal loading arising, for example, from ice flows.

While modern software is making it possible to analyse ever-increasingly complex problems, the same
tools are also making it possible to better understand the limit equilibrium method itself, particularly the
differences between the various analysis methods and the limitations in the formulation. Knowledge of
these limitations can be used to interpret the results with confidence.

It is important to note that the purpose of this ‘book’ is not to provide detailed instructions for
operating the software. The primary vehicle for that information is the learning section for GeoStudio on
the Seequent website (www.seequent.com), where the user can access example files, tutorial movies,
technical webinars and our online learning platform. In addition, help topics are available during
operation of the software in the Help menu (accessed by pressing F1). These resources provide valuable
information for those learning how to use GeoStudio.

This reference book covers the key elements to conduct a slope stability analysis using SLOPE/W or
SLOPE3D. The first section of this book reviews the background information and theory that underpin
these two GeoStudio products. A fundamental understanding of the limit equilibrium method of
columns is essential to conducting an effective slope stability analysis with GeoStudio. As such, an entire
chapter is devoted to a review of the fundamentals of limit equilibrium as a method of analysis. The
chapter reviews the consequences of a static equilibrium formulation, the differences between the
various methods, the importance of intercolumn forces, and the effect of slip surface shape on the

1
computed factor of safety. In addition, the chapter discusses the limitations of the method and ways of
overcoming these limitations.

Chapters 3 through 7 summarize the main components of stability analyses: slip surface search
methods, material definition, external loads, and pore water pressure definition. Additional analysis
options are discussed in Chapters 8 to 10, which cover seismic and dynamic analyses, probabilistic and
sensitivity analyses, and partial factor design. The ultimate objective of these chapters is to help the
reader understand the fundamentals of each component or analysis type to assist with setting up a
stability analysis and utilizing GeoStudio’s extended capabilities.

2
2 Theory

2.1 Background
There are many different solution techniques for the limit equilibrium method for analysing slope
stability. The methods are characterized by 1) the equations of statics that are satisfied; and, 2) the
assumed relationship between the intercolumn shear and normal forces. To better understand the
different methods, consider a 2D sliding mass discretized into 2D columns and a basic free-body diagram
for a single column, as illustrated in Figure 1. Normal and shear forces act on the base and sides of the
column, while the intercolumn shear and normal forces act on the vertical edges. The sliding mass is
considered to be at the state of limiting equilibrium (i.e. just at the point of failure) such that the sum of
all forces and/or moments acting on all columns are equal to zero.

Figure 1. Column discretization and column forces for a sliding mass.

The first method developed was known as the Ordinary, or Fellenius, method. The method ignored all
intercolumn forces and satisfied only moment equilibrium. Adopting these simplified assumptions made
it possible to compute a factor of safety using hand calculations. Later Bishop (1955) devised a scheme
that included intercolumn normal forces but ignored the intercolumn shear forces. Again, Bishop’s
Simplified method satisfies only moment equilibrium. Including the normal intercolumn forces yields a
non-linear factor of safety equation that requires an iterative solution procedure. The Janbu Simplified
method is like the Bishop Simplified method because it includes the normal intercolumn forces and
ignores the intercolumn shear forces. However, the Janbu method satisfies horizontal force equilibrium,
as opposed to moment equilibrium.

Computers made it possible to readily handle the iterative procedure inherent in the limit equilibrium
method, leading to mathematically rigorous formulations that included all intercolumn forces and
satisfied all equations of statics. Two such methods are the Morgenstern-Price and Spencer methods.
Table 1 lists the methods available in GeoStudio, along with the equations of statics that are satisfied for
each of the methods. Table 2 provides a summary of the intercolumn forces included in each method
and the assumed relationship between the intercolumn shear and normal forces.

3
Table 1. Equations of statics satisfied for each method.

Method 2D 3D Moment
Force Equilibrium
Equilibrium
Ordinary or Fellenius Yes No Yes No
Bishop Simplified Yes Yes Yes No
Janbu Simplified Yes Yes No Yes
Spencer Yes Yes Yes Yes
Morgenstern-Price Yes Yes Yes Yes
Corps of Engineers – 1 Yes No No Yes
Corps of Engineers – 2 Yes No No Yes
Lowe-Karafiath Yes No No Yes
Janbu Generalized Yes No Yes (by column) Yes
Sarma – vertical columns Yes No Yes Yes

Table 2. Intercolumn force characteristics and relationships.

Method Intercolumn Intercolumn X/E Relationship/

Normal (E) Shear (X) Inclination of Resultant
Ordinary or Fellenius No No No intercolumn forces
Bishop Simplified Yes No None/Horizontal
Janbu Simplified Yes No None/Horizontal
Spencer Yes Yes Prescribed/Constant
Morgenstern-Price Yes Yes User Function/Variable
Corps of Engineers – 1 Yes Yes Prescribed/Inclination of a
line from crest to toe
Corps of Engineers – 2 Yes Yes Prescribed /Inclination of
ground surface at top of
column
Lowe-Karafiath Yes Yes Prescribed /Average of
ground surface and column
base inclination
Janbu Generalized Yes Yes Prescribed /Applied line of
thrust and moment
equilibrium of column
Sarma – vertical columns Yes Yes Prescribed/
𝑋 = 𝑐 + 𝐸𝑡𝑎𝑛𝜙'

The relationship between the intercolumn shear and normal forces is one of the key differentiators for
the methods. The Spencer method assumes that 𝑋 𝐸 = 1.0, while Morgenstern-Price assumes a half-
sine relationship, which means the intercolumn shear is zero at the crest and toe of the slip surface,

4
while the shear and normal are equal at the centre of the slip surface. The Corps of Engineers method 1
uses the inclination (i.e., rise/run) of a line from crest to toe to determine the ratio, while the Corps of
Engineers method 2 uses the ground surface inclination (Figure 2 to Figure 4). The Lowe-Karafiath
method is similar, except that it is uses the average of ground surface and column base inclinations. The
Sarma – vertical column method uses the relationship presented in Table 2, making the intercolumn
shear force a function of the cohesion and friction angle.

1.144

Figure 2. Example slip surface.

Interslice Force Fn. vs. Distance

0.4

Applied Fn.
0.3
Interslice Force Fn.

0.2

0.1

Specified Fn.

0.0
0 10 20 30 40 50

Distance

Figure 3. Intercolumn force function for Corps. of Engineers 1.

5
 0.5

0.4 Applied Fn.

Interslice Force Fn.

0.3

0.2

0.1
Specified Fn.

0.0
0 10 20 30 40
Distance

Figure 4. Intercolumn force function for Corps. of Engineers 2.

2.2 Factor of Safety Definition

The factor of safety is defined as the ratio of shear strength of the soil/rock to the shear stress mobilized
in the system. The factor of safety is given by:

𝑠 Equation 1
𝐹=
𝑠𝑚

where 𝑠 is the shear strength and 𝑠𝑚 is the mobilized shear stress. In the context of the limit equilibrium
method, it is often more insightful to write the equation as follows:
𝑠 Equation 2
𝑠𝑚 =
𝐹

In this form, it is apparent that the factor of safety is the factor by which the shear strength of the
soil/rock must be reduced to bring the system into a state of limiting equilibrium. Section 2.4
demonstrates that this expression is the cornerstone of the formulation because it introduces 𝐹 into the
equations of equilibrium.

For an effective stress analysis, the shear strength of an unsaturated soil is defined as:
Equation 3
𝑠 = 𝑐' + (𝜎𝑛 ‒ 𝑢𝑎)𝑡𝑎𝑛𝜙' + 𝑆𝑒(𝑢𝑎 ‒ 𝑢𝑤)𝑡𝑎𝑛𝜙'

where 𝑐 is the effective cohesion, 𝜙 is the effective angle of internal friction, 𝜎𝑛 is the total normal
' '

stress, 𝑢𝑤 is the pore water pressure, 𝑢𝑎 is the pore-air pressure, and 𝑆𝑒 is the effective saturation given
as (van Genuchten, 1980):

6
 𝜃𝑤 ‒ 𝜃𝑟
Equation 4
𝑆𝑒 =
𝜃𝑠 ‒ 𝜃𝑟

where 𝜃𝑤, 𝜃𝑠, and 𝜃𝑟 are the volumetric water content, saturated volumetric water content, and
residual volumetric water content, respectively. The last term of Equation 3 represents the additional
strength derived from matric suction in the unsaturated zone. For a saturated soil, the strength equation
simplifies to:
Equation 5
𝑠 = 𝑐' + (𝜎𝑛 ‒ 𝑢𝑤)𝑡𝑎𝑛𝜙'

2.3 Column Forces

The sliding mass is discretized into columns for 2D and 3D analysis. GeoStudio adopts the coordinate
system shown in Figure 5; the y-direction is positive upwards, and North is aligned with the negative z-
axis (South being aligned with the positive z-axis).

Figure 5. Coordinate system in GeoStudio. North is aligned with the negative z-direction.

Figure 6 presents a plan view depiction of a 3D sliding mass discretized into columns. The sliding
direction is unknown in a 3D analysis and is opposite of the shear mobilized force direction (𝛼). For 2D
‒𝜋 𝜋
𝛼
analysis, the sliding direction is either left-to-right or right-to-left, in which case equals 2 or 2 ,
respectively. Figure 7 and Figure 8 top and side views of the forces acting on a 3D column. Table 3
defines the symbols.

7
Figure 6. Plan view of a 3D slip surface discretized into columns.

𝑦 Back 𝑥

𝐸𝑧�

𝐻𝑥�

𝐻𝑧� 𝑁�𝑔��
𝐸𝑥� 𝑁�𝑔�� 𝐸𝑥���
Left ∆𝑧� Right
𝑆�𝑓��
𝑆�𝑓�� 𝐻𝑧���

𝐻𝑥���

𝐸𝑧���
𝑧
∆𝑥�
Front

Figure 7. Plan view of the forces acting on a 3D column.

8
 ∆𝑥� ∆𝑧�

1 + 𝑘� 𝑊� 1 + 𝑘� 𝑊�

𝑦 𝑦

𝑋𝑥� 𝐸𝑥��� 𝑋𝑧� 𝐸𝑧���

∆𝐻𝑥� ∆𝐻𝑧�
𝐸𝑥� Right 𝐸𝑧� Front
𝑋𝑥��� 𝑋𝑧���
Left Back
𝑁�𝑔�� 𝑁�𝑔��
𝑆�𝑓�� 𝑆�𝑓��
𝑁�𝑔�� 𝑁�𝑔��

𝑆�𝑓�� 𝑆�𝑓��
𝑆��(�) 𝑆��(�)

𝑁� 𝑁�
𝑥 𝑧

Figure 8. Section views of the forces acting on a 3D column.

Table 3. Summary of column forces.

Parameter Definition
(1 + 𝑘 𝑣 )𝑊 Total weight of column, including vertical seismic load
𝑘 ℎ𝑊 Horizontal seismic load applied through the centroid of
the column
𝑁 Total normal force on the base of then column
𝑆 Shear force mobilized on the base of the column
𝐸 Intercolumn normal forces
𝑋 Intercolumn vertical shear forces
𝐻 Intercolumn horizontal shear forces
𝐷 External point loads

The normal and shear forces on the base of each column are resolved into x-, y-, and z- components:

𝑁𝑥𝑖 = 𝑔1 𝑁𝑖; 𝑁𝑦𝑖 = 𝑔2 𝑁𝑖; 𝑁𝑍𝑖 = 𝑔3 𝑁𝑖 Equation 6

𝑖 𝑖 𝑖

𝑆𝑥𝑖 = 𝑓1 𝑆𝑖; 𝑆𝑦𝑖 = 𝑓2 𝑆𝑖; 𝑆𝑍𝑖 = 𝑓3 𝑆𝑖 Equation 7

𝑖 𝑖 𝑖

where:

{𝑔 1 𝑖 𝑔2
𝑖
𝑔3
𝑖 }= unit vector in the direction of base normal force 𝑁𝑖

9
{𝑓 1 𝑖 𝑓2
𝑖
𝑓3
𝑖 } = unit vector in the direction of mobilized shear force 𝑆𝑖

The subscript 𝑖 denotes the column number, while the subscripts 1, 2, and 3 denote the components of
the unit vector relative to x-, y-, and z- axes, respectively. The unit vector for the shear force is opposite
of the sliding direction. An arbitrary inter-column force function 𝑓(𝑥,𝑧) is assumed to relate the
intercolumn shear and normal forces in the x- and z- directions:

𝑋𝑥𝑖 = 𝐸𝑥𝑖𝑓(𝑥,𝑧)𝜆𝑥; 𝑋𝑧𝑖 = 𝐸𝑧𝑖𝑓(𝑥,𝑧)𝜆𝑧 Equation 8

𝐻𝑥𝑖 = 𝐸𝑧𝑖𝑓(𝑥,𝑧)𝜆𝑥𝑧; 𝐻𝑧𝑖 = 𝐸𝑥𝑖𝑓(𝑥,𝑧)𝜆𝑧𝑥 Equation 9

where 𝐸 represents the intercolumn normal forces, 𝑋 is the intercolumn vertical shear force, and 𝐻 is
the intercolumn horizontal shear force.
𝐻𝑥 𝐻𝑧
In general, the intercolumn horizontal shear forces 𝑖 and 𝑖 are zero for the exterior columns. By
using the property of complementary shear, which is founded on moment equilibrium in the 𝑥-𝑧 plane,
𝐻𝑥 𝐻𝑧
and the specified 𝜆𝑥𝑧 or 𝜆𝑧𝑥, 𝑖 + 1 and 𝑖 + 1 can be determined sequentially as:
Δ𝑥𝑖
Equation 10
𝐻𝑥
𝑖+1
=
Δ𝑧𝑖
(𝐻𝑧 𝑖 + 1 + 𝐻𝑧 𝑖 ) ‒ 𝐻𝑥 𝑖

Δ𝑧𝑖
Equation 11
𝐻𝑧
𝑖+1
=
Δ𝑥𝑖
(𝐻𝑥 𝑖 + 1 + 𝐻𝑥 𝑖 ) ‒ 𝐻𝑧 𝑖

where, Δ𝑥𝑖 and Δ𝑧𝑖 are column spacing in the 𝑥-𝑧 plane. The column spacing is present because the
equations are derived from moment equilibrium about a y-axis located at the center of the column. If
the columns are spaced equally in both the x- and y- directions, then the first term is 1.0.

Only the left-most diagram in Figure 8 is applicable for 2D because the forces on the front and back of
𝐻𝑥 𝐻𝑧
the column are implicitly zero. Moreover, the intercolumn horizontal shear forces, 𝑖 and 𝑖, are zero
for a 2D analysis because there is no tendency for rotation about the y-axis. The weight of column 𝑖 is
represented as 𝑊𝑖. The shear force and normal force at the base of the column in plan view (i.e., z- 𝑥
𝑆𝑧𝑥 𝑁
plane) is given by 𝑖 and 𝑧𝑥𝑖, respectively. The base column forces are similarly shown in the two
𝑆𝑥𝑦 𝑁𝑥𝑦 𝑆𝑥𝑦 𝑁
section views as 𝑖, 𝑖, 𝑖 and 𝑥𝑦𝑖, where the subscript indicates the plane (e.g., 𝑥-y plane).
Intercolumn normal forces are always horizontal, are visible on both the left/right and front/back sides
and are aligned with the 𝑥-axis and z-axis, respectively.

2.3.1 Base Normal Force

The normal force at the base of a column is derived from the summation of forces in the vertical
direction as:

10
 𝑔2 𝑁𝑖 + 𝑓2 𝑆𝑖 ‒ (1 + 𝑘𝑣)𝑊𝑖 + ∆𝑋𝑥𝑖 + ∆𝑋𝑧𝑖 + 𝐷𝑦𝑖 = 0 Equation 12
𝑖 𝑖

where ∆𝑋𝑥𝑖 and ∆𝑋𝑧𝑖 represent the resultants of the intercolumn vertical shear forces. The shear force
mobilized at the base of a column (𝑆𝑖) can be written considering the area of the column base, 𝛽, and
ignoring the unsaturated strength as:

𝑠 𝑖𝛽 𝑖 𝛽𝑖 𝐶𝑖 + (𝑁𝑖 ‒ 𝑈𝑖)𝑡𝑎𝑛𝜙𝑖' Equation 13

𝑆𝑖 = 𝑠 𝑚 𝛽 𝑖 =
𝑖 𝐹
=
𝐹 ( '
( )
𝑐 𝑖 + 𝜎𝑛 ‒ 𝑢𝑤 𝑡𝑎𝑛𝜙 =
𝑖 𝑖
'
) 𝐹

where 𝐶𝑖 and 𝑈𝑖 are the cohesive force and pore water pressure force, respectively. After substituting
for the mobilized shear force 𝑆𝑖 into Equation 12, the base normal force can be expressed as:
𝑁 𝑖 = 𝐴𝑖 + 𝐵𝑖 𝑆𝑖 Equation 14

where

(1 + 𝑘𝑣)𝑊𝑖 ‒ ∆𝑋𝑥𝑖 ‒ ∆𝑋𝑧𝑖 ‒ 𝐷𝑦𝑖 Equation 15

𝐴𝑖 =
𝑔2
𝑖

‒ 𝑓2 Equation 16
𝑖
𝐵𝑖 =
𝑔2
𝑖

𝐶𝑖 + (𝐴𝑖 ‒ 𝑈𝑖)𝑡𝑎𝑛𝜙𝑖' Equation 17

𝑆𝑖 =
𝑡𝑎𝑛𝜙𝑖'
[ ( )]
𝐹 1 ‒ 𝐵𝑖
𝐹

The equation for the base normal force is nonlinear because it depends on the factor of safety, 𝐹. Since
𝐹 and the intercolumn shear forces, 𝑋𝑥 and 𝑋𝑧, are unknown, the base normal force cannot be solved
directly, and an iterative scheme is required (Section 2.5). Finally, it is noted that Equation 14 is also
𝑓2
valid for a 2D analysis given that 𝑖 is simply the unit vector component of the normal force in the
vertical direction.

2.3.2 Intercolumn Normal Force

The intercolumn normal forces, 𝐸𝑥𝑖 and 𝐸𝑧𝑖, are required for the calculation of intercolumn shear forces,
which are required for the calculation of the base normal force. The intercolumn normal forces are
solved using an integration procedure commencing from left-to-right for a 2D slip surface, and following
the patterns presented in Figure 9 for 3D. Summation of forces in the x-direction for each individual
column is given by:

11
 Equation 18
𝐸𝑥
𝑖+1 𝑖
(
= 𝐸 𝑥 + 𝐻𝑥
𝑖+1 𝑖
)
‒ 𝐻𝑥 + 𝑓1 𝑆𝑖 + 𝑔1 𝑁𝑖 + 𝐷𝑥𝑖
𝑖 𝑖

Similarly, the intercolumn normal force in the z-direction is given by:

Equation 19
𝐸𝑧
𝑖+1 𝑖
(
= 𝐸 𝑧 + 𝐻𝑧
𝑖+1 𝑖
)
‒ 𝐻𝑧 + 𝑓2 𝑆𝑖 + 𝑔2 𝑁𝑖 + 𝐷𝑧𝑖
𝑖 𝑖

Since the outer-most intercolumn forces are zero, the intercolumn normal forces of all columns can be
computed in a column-by-column and row-by-row sequence commencing from the exterior corner
columns (Figure 9). Note that the equations for computing the intercolumn normal forces are
dependent on the normal and mobilized shear forces at the column base and hence the computation of
intercolumn normal forces is part of the iterative solution process (Section 2.5). In a 2D analysis, only
Equation 18 is required, with the intercolumn horizontal shear forces once again zeroed.

Figure 9. Integration procedure for intercolumn normal force calculation in 3D. a) Summation of forces in the x-direction. b)
Summation of forces in the z-direction.

2.3.3 Intercolumn Shear Force

The intercolumn shear forces in the rigorous slope stability formulations are handled with Equation 8
and Equation 9 as proposed by Morgenstern and Price (1965). In these equations, 𝜆 is a scaling factor
required for equilibrium and 𝑓(𝑥, 𝑧) is the intercolumn force function for the x- and z- Cartesian
directions. The intercolumn force function depends on the slope stability method (see Table 2) and
represents the ratio between the intercolumn shear and normal forces following:
𝑋𝑥 𝑋𝑧
𝑓 (𝑥 ) = Equation 20
𝐸 𝑥 ; 𝑓 (𝑧 ) = 𝐸𝑧

The Spencer method assumes a constant intercolumn force function, while the Morgenstern and Price
method assumes a half-sine relationship. The scaling factor, 𝜆, is unique in the x- and z-directions for a
3D analysis, while there is only a single 𝜆 value for a 2D analysis commensurate with the x-direction.

For illustration, consider the typical half-sine force function vs column number shown in Figure 10 for a
2D analysis. The upper curve is the specified half-sine function. The lower curve is the function after the

12
factor of safety has been solved to convergence. The ratio between the two curves represents 𝜆, which
in this case is 0.43. At Column 10, 𝑓(𝑥) = 0.83, 𝐸 = 100 𝑘𝑁, so 𝑋 = (100)(0.43)(0.83) = 35.7 𝑘𝑁. The
‒1
resultant angle between the normal and shear force would therefore be 19.6tan [𝑋 𝐸]. This means
the intercolumn resultant force is inclined at 19.6 degrees from the horizontal at Column 10.

1.2

1.0
Interslice Force Functions

0.8

0.6

0.4

0.2

0.0
0 5 10 15 20 25 30

Slice #

Applied Fn. Specified Fn.

Figure 10. Intercolumn force function for a typical 2D analysis.

2.4 Limit Equilibrium Formulation

A general limit equilibrium formulation uses the following equations of statics in solving for the factor of
safety:

 The summation of forces in the vertical direction for each column is used to compute the normal
force, 𝑁 (Section 2.3.1).
 The summation of forces in the horizontal directions (𝑥 and z) for each column is used to
compute the intercolumn normal force, 𝐸. In 2D, summation only occurs in the x-direction.
 The summation of moments about a centre of rotation in both the x- and z- directions for all
columns is calculated. The equations can be rearranged and solved for the moment equilibrium
factors of safety, 𝐹𝑚𝑥.and 𝐹𝑚𝑧. In 2D, 𝐹𝑚𝑥 is the only component considered.
 The summation of horizontal forces in both the x- and z- directions for all columns produces the
force equilibrium factors of safety, 𝐹𝑓𝑥 and 𝐹𝑓𝑧. In 2D, only the 𝐹𝑓𝑥 equation is considered.

The 2D and 3D formulations are fundamentally the same. It is convenient to think of a 3D analysis as
two dependent 2D analyses. In a 2D analysis, the coordinate axis is x-y, there are no forces acting on the
back and front of the columns, and the shear and normal forces reside in the x-z plane (i.e., no z-force
components). In a 3D analysis, equilibrium is first satisfied in the x-y plane, then in the z-y plane as
shown in Figure 8. The intercolumn horizontal shear forces on the opposing faces of a column ‘bind’
these two directions together. For example, horizonal force equilibrium in the x-direction (left-right)

13
 𝐻𝑥 𝐻𝑥
requires the intercolumn shear acting on the back and front of the column, 𝑖 and 𝑖 + 1. Finally, it is
noted that forces in 3D are resolved into x, y, and z components, while only the x and y components are
active in 2D.

2.4.1 Moment Equilibrium

Consider the overall moment equilibrium in the x-direction including point loads, 𝐷, shown below. In
this case, the moments are taken about an axis parallel to the z-axis.

∑(𝑁𝑖𝑔2𝑖 + 𝑆𝑖𝑓2𝑖 ‒ (1 + 𝑘𝑣)𝑊𝑖)𝑑𝑥(𝑖) ‒ ∑(𝑁𝑖𝑔1𝑖 + 𝑆𝑖𝑓1𝑖)𝑑𝑦(𝑖) + ∑(𝐷𝑗𝑦𝑑𝑥(𝑗) ‒ 𝐷𝑗𝑥𝑑𝑦(𝑗)) = 0 Equation 21

where 𝑑𝑥(𝑖)/𝑑𝑥(𝑗) and 𝑑𝑦(𝑖)/𝑑𝑦(𝑗) are the moment arms for the y- and x-components of each force,
respectively. Note that the intercolumn normal and shear forces are not present in Equation 21 because
summation over all columns results in cancellation. After substituting for the mobilized shear force 𝑆𝑖
(Equation 13) and rearranging, the moment factor of safety in the x- equilibrium direction is:
Equation 22
∑{[𝐶𝑖 + (𝑁𝑖 ‒ 𝑈𝑖)𝑡𝑎𝑛𝜙𝑖'](𝑓2𝑖𝑑𝑥(𝑖) ‒ 𝑓1𝑖𝑑𝑦(𝑖))}
𝐹𝑚𝑥 =
∑((1 + 𝑘𝑣)𝑊𝑖 ‒ 𝑁𝑖𝑔2𝑖)𝑑𝑥(𝑖) + ∑𝑁𝑖𝑔1𝑖𝑑𝑦(𝑖) ‒ ∑(𝐷𝑗𝑦𝑑𝑥(𝑗) ‒ 𝐷𝑗𝑥𝑑𝑦(𝑗))

Equation 22 can be written to include reinforcement forces, 𝑇𝑘, as follows:

Equation 23

𝐹𝑚𝑥
∑{[𝐶𝑖 + (𝑁𝑖 ‒ 𝑈𝑖)𝑡𝑎𝑛𝜙𝑖'](𝑓2𝑖𝑑𝑥(𝑖) ‒ 𝑓1𝑖𝑑𝑦(𝑖))}
=
∑((1 + 𝑘𝑣)𝑊𝑖 ‒ 𝑁𝑖𝑔2𝑖)𝑑𝑥(𝑖) + ∑𝑁𝑖𝑔1𝑖𝑑𝑦(𝑖) ‒ ∑(𝐷𝑗𝑦𝑑𝑥(𝑗) ‒ 𝐷𝑗𝑥𝑑𝑦(𝑗)) ‒ ∑

In some cases, the factor of safety is applied to the reinforcement loads (Section 7.8), yielding the factor
of safety dependent equation:
Equation 24
∑{[𝐶𝑖 + (𝑁𝑖 ‒ 𝑈𝑖)𝑡𝑎𝑛𝜙𝑖'](𝑓2𝑖𝑑𝑥(𝑖) ‒ 𝑓1𝑖𝑑𝑦(𝑖))} ‒ ∑(𝑇𝑘𝑦𝑑𝑥(𝑘) ‒ 𝑇𝑘𝑥𝑑𝑦(𝑘))
𝐹𝑚𝑥 =
∑((1 + 𝑘𝑣)𝑊𝑖 ‒ 𝑁𝑖𝑔2𝑖)𝑑𝑥(𝑖) + ∑𝑁𝑖𝑔1𝑖𝑑𝑦(𝑖) ‒ ∑(𝐷𝑗𝑦𝑑𝑥(𝑗) ‒ 𝐷𝑗𝑥𝑑𝑦(𝑗))

For a 3D analysis, moment equilibrium must also be satisfied in the z-direction. In this case, the
moments are taken about an axis parallel to the x-axis. The moment factor of safety in the z-direction is:
Equation 25

14
 𝐹𝑚𝑧
∑{[𝐶𝑖 + (𝑁𝑖 ‒ 𝑈𝑖)𝑡𝑎𝑛𝜙𝑖'](𝑓3𝑖𝑑𝑦(𝑖) ‒ 𝑓2𝑖𝑑𝑧(𝑖))}
=
∑((1 + 𝑘𝑣)𝑊𝑖 ‒ 𝑁𝑖𝑔2𝑖)𝑑𝑧(𝑖) + ∑𝑁𝑖𝑔3𝑖𝑑𝑦(𝑖) + ∑(𝐷𝑗𝑧𝑑𝑦(𝑗) ‒ 𝐷𝑗𝑦𝑑𝑧(𝑗)) + ∑

Note that the moment arm 𝑑𝑦(𝑖) applies to the z- component of a force, while 𝑑𝑧(𝑖) applies to the y-
component of a force. The factor of safety dependent form of Equation 25 is given by:
Equation 26
∑{[𝐶𝑖 + (𝑁𝑖 ‒ 𝑈𝑖)𝑡𝑎𝑛𝜙𝑖'](𝑓3𝑖𝑑𝑦(𝑖) ‒ 𝑓2𝑖𝑑𝑧(𝑖))} + ∑(𝑇𝑘𝑧𝑑𝑦(𝑘) ‒ 𝑇𝑘𝑦𝑑𝑧(𝑘))
𝐹𝑚𝑧 =
∑((1 + 𝑘𝑣)𝑊𝑖 ‒ 𝑁𝑖𝑔2𝑖)𝑑𝑧(𝑖) + ∑𝑁𝑖𝑔3𝑖𝑑𝑦(𝑖) + ∑(𝐷𝑗𝑧𝑑𝑦(𝑗) ‒ 𝐷𝑗𝑦𝑑𝑧(𝑗))
Only Equation 23 or Equation 24 is required for a 2D analysis.

2.4.2 Force equilibrium

Consider the overall force equilibrium in the x-direction.
Equation 27
∑∆𝐻𝑥𝑖 + ∑𝐷𝑗𝑥 + ∑𝑁𝑖𝑔1𝑖 + ∑𝑆𝑖𝑓1𝑖 = 0

Here, ∆𝐻𝑥𝑖 represents the resultant of the intercolumn horizontal shear forces in the x-direction. After
substituting for the mobilized shear force, 𝑆𝑖, and rearranging the terms, the force factor of safety in the
x- equilibrium direction is:
Equation 28
∑[(𝑁𝑖 ‒ 𝑈𝑖)𝑡𝑎𝑛𝜙𝑖' + 𝐶𝑖]𝑓1𝑖
𝐹𝑓𝑥 =‒
∑𝐷𝑗𝑥 + ∑𝑁𝑖𝑔1𝑖 + ∑∆𝐻𝑥𝑖 + ∑𝑇𝑘𝑥
Reinforcement forces, 𝑇𝑘, can again be considered according to factor of safety dependency, yielding:
Equation 29
∑[(𝑁𝑖 ‒ 𝑈𝑖)𝑡𝑎𝑛𝜙𝑖' + 𝐶𝑖]𝑓1𝑖 + ∑𝑇𝑘𝑥
𝐹𝑓𝑥 =‒
∑𝐷𝑗𝑥 + ∑𝑁𝑖𝑔1𝑖 + ∑∆𝐻𝑥𝑖
Similarly, the equations of factor of safety for force equilibrium in the 𝑧-direction are:
Equation 30
∑[(𝑁𝑖 ‒ 𝑈𝑖)𝑡𝑎𝑛𝜙𝑖' + 𝐶𝑖]𝑓3𝑖
𝐹𝑓𝑧 =‒
∑𝐷𝑗𝑧 + ∑𝑁𝑖𝑔3𝑖 + ∑∆𝐻𝑧𝑖 + ∑𝑇𝑘𝑧
and factor of safety dependent:

15
 Equation 31
∑[(𝑁𝑖 ‒ 𝑈𝑖)𝑡𝑎𝑛𝜙𝑖' + 𝐶𝑖]𝑓3𝑖 + ∑𝑇𝑘𝑧
𝐹𝑓𝑧 =‒
∑𝐷𝑗𝑧 + ∑𝑁𝑖𝑔3𝑖 + ∑∆𝐻𝑧𝑖
𝐻𝑥
Only Equation 28 or Equation 29, with the horizontal intercolumn shear forces 𝑖 eliminated from the
equations, is required for a 2D stability analysis.

2.5 Solution Algorithm

The derivation of the force and moment factor of safety equations produces a set of non-linear
equations that must be solved using a Newton-Raphson solution procedure. The non-linearity is
manifest in the equations as follows:

 Normal force at the base of a column is a function of 𝐹;

 𝐹 is a function of the normal force at the base of a column;
 Intercolumn normal force is indirectly a function of 𝐹 via the normal force at the base of the
column and directly via the shear mobilized; and,
 Intercolumn shear force is a function of the intercolumn normal force, which is a function of 𝐹.

The solution schemes used to solve the 2D and 3D equations are discussed in the following sections.

2.5.1 2D Solution Workflow

Figure 11 presents the solution workflow for a 2D analysis. The solution scheme consists of the following
essential steps (𝑖 = 𝑖𝑡𝑒𝑟𝑎𝑡𝑖𝑜𝑛 𝑛𝑢𝑚𝑏𝑒𝑟):

1. Initialize the parameters 𝜆 = 0 and 𝐹𝑓 = 𝐹𝑚 = 1.0.

2. Search for 𝜆 that results in 𝐹𝑓 = 𝐹𝑚, iterating until 𝐹𝑓(𝑖) = 𝐹𝑓(𝑖 ‒ 1) and 𝐹𝑚(𝑖) = 𝐹𝑚(𝑖 ‒ 1) for each
value of 𝜆.
3. The calculation procedure stops when 𝐹𝑓 = 𝐹𝑚.

16
Figure 11. Solution workflow for a 2D analysis.

Figure 12 illustrates the plot used to judge convergence. At 𝜆 = 0 the intercolumn shear forces are zero
(Equation 8 and Equation 9) and 𝐹𝑚 and 𝐹𝑓 correspond to the Bishop and Janbu solutions, respectively.
Convergence of the inner loop (Figure 11) for any value of 𝜆 implies that the shear strength at the base
of each column was divided by a value 𝐹𝑓 or 𝐹𝑚 to bring each column, and therefore the entire sliding
mass, into force or moment equilibrium. Convergence of the outer loop implies that the specified inter-
column force function was scaled by the correct value of 𝜆 to obtain 𝐹𝑓 = 𝐹𝑚. This condition
corresponds to the intersection point of the two curves shown in Figure 12 and implies that the column
forces calculated using 𝐹𝑓 or 𝐹𝑚 would be equivalent. Force and moment equilibrium of each column,
and therefore the entire sliding mass, has been obtained.

1.15

Bishop
1.10
Factor of Safety

1.05

Morgenstern-Price
or Spencer
1.00

0.95

Janbu

0.90
0.0 0.1 0.2 0.3 0.4 0.5 0.6

Lambda

Moment Force

Figure 12. Factor of safety versus lambda for the critical slip surface.

2.5.2 3D Solution Workflow

Figure 13 presents the solution workflow for a 3D analysis. There are two loops for solving the x and z
equilibrium equations that are wrapped inside a 𝜆𝑥 ‒ 𝜆𝑧 dependency loop, which is wrapped in an 𝛼
loop. The 𝜆𝑥 ‒ 𝜆𝑧 lambda dependency loop is necessary because the solution for 𝜆𝑥 depends on the
value of 𝜆𝑧 and vice versa. The sliding direction 𝛼 loop is necessary to ensure that the force and moment
factors of safety in both equilibrium directions are equivalent. The solution scheme consists of the
following essential steps (𝑖 and 𝑗 are iteration counters):
𝑥 𝑥 𝑧 𝑧
1. Initialize the parameters 𝜆𝑧 = 0, 𝐹𝑓 = 𝐹𝑚 = 𝐹𝑓 = 𝐹𝑚 = 1.0, and the sliding direction 𝛼.

17
 𝑥 𝑥
2. For the current value of 𝜆𝑧 and 𝛼, search for 𝜆𝑥 that results in 𝐹𝑓 = 𝐹𝑚, iterating until
𝑥
𝐹𝑓(𝑖) = 𝐹𝑓(𝑖 𝑥‒ 1) and 𝐹𝑚(𝑖)
𝑥
= 𝐹𝑚(𝑖𝑥‒ 1) for each value of 𝜆𝑥.
𝑧 𝑧
3. For the current value of 𝜆𝑥 and 𝛼, search for 𝜆𝑧 that results in 𝐹𝑓 = 𝐹𝑚, iterating until
𝑧
𝐹𝑓(𝑖) = 𝐹𝑓(𝑖 𝑧‒ 1) and 𝐹𝑚(𝑖)
𝑧
= 𝐹𝑚(𝑖𝑧‒ 1) for each value of 𝜆𝑧.

4. Repeat steps 2 and 3 until 𝜆𝑥(𝑗) = 𝜆𝑥(𝑗 ‒ 1) and 𝜆𝑧(𝑗) = 𝜆𝑧(𝑗 ‒ 1), then proceed to the next step.
𝑥 𝑧 𝑥 𝑧
5. Calculate an estimate of the sliding direction 𝛼 that results in 𝐹𝑓 = 𝐹𝑓 and 𝐹𝑚 = 𝐹𝑚 and repeat
steps 2 onwards.
𝑥 𝑧 𝑥 𝑧
6. The calculation procedure stops when 𝐹𝑓 = 𝐹𝑓 and 𝐹𝑚 = 𝐹𝑚.

Figure 13. Solution workflow for a 3D analysis.

Figure 14 presents the convergence graphs from a 3D stability analysis, including a) 𝐹 vs 𝜆 in the x-
direction; b) 𝐹 vs 𝜆 in the z- direction; and, c) 𝐹 vs 𝛼. At convergence, the outer 𝜆 loops ensures that
𝐹𝑚𝑥 = 𝐹𝑓𝑥 and 𝐹𝑚𝑧 = 𝐹𝑓𝑧. The 𝛼 loop ensures that the 𝐹𝑥 = 𝐹𝑧, meaning that the entire sliding mass is
in equilibrium and therefore the correct sliding direction has been determined. In both loops,
convergence is acceptable if the 𝜆 cross-over(s) and 𝛼 cross-over are unambiguous.

18
Figure 14. Convergence plots for a 3D stability analysis.

2.6 Convergence and the Role of Slip Surface Shape

The effect of the intercolumn force function on the calculated factor of safety depends partly on the
shape of the slip surface. Intuitively, this makes sense given that the amount of rotation and/or
translation within a sliding mass would govern the amount of shear that develops between the columns.
This section discusses the effect of 𝜆 on the moment and force factor of safety values, and therefore the
calculated (converged) 𝐹. For simplicity, only 2D results are presented, but the concepts also apply for
3D stability analysis where the slip surface shapes are ellipsoidal or spherical.

2.6.1 Circular Slip Surface

Figure 15 presents a simple 2D circular slip surface together with the associated 𝐹 versus 𝜆 plot. In this
case the moment equilibrium is completely independent of the intercolumn shear forces, as indicated
by the horizontal moment equilibrium curve. The force equilibrium, however, is dependent on the
intercolumn shear forces.

19
 1.35

1.30
Factor of safety

1.25

1.20

1.15

1.10

1.05
0 0.1 0.2 0.3 0.4 0.5

Lambda

Moment Force

Figure 15. Factor of safety versus lambda for a circular slip surface.

The moment equilibrium is not influenced by the shear forces because the sliding mass as a free body
can rotate without any movement between the columns. However, substantial intercolumn slippage is
necessary for the sliding mass to move laterally. Consequently, the horizontal force equilibrium is
sensitive to intercolumn shear.

Any assumption regarding an intercolumn force function is irrelevant for this scenario because the
moment equilibrium is completely independent of intercolumn shear. The intercolumn shear can be
assumed to be zero, as in the Bishop’s Simplified method, and still obtain an acceptable factor of safety,
provided the method satisfies moment equilibrium. This is, of course, not true for a method based on
satisfying only horizontal force equilibrium such as the Janbu’s Simplified method. Ignoring the
intercolumn shear when only horizontal force equilibrium is satisfied for a curved slip surface results in a
factor of safety significantly different than when both force and moment equilibrium is satisfied.

2.6.2 Planar Slip Surface

Figure 16 illustrates a planar slip surface. The moment and force equilibrium curves have reverse
positions compared to a circular slip surface. The force equilibrium is completely independent of
intercolumn shear, while moment equilibrium is sensitive to the intercolumn shear. The soil wedge can
move or translate along the planar slip surface without any slippage between the columns. Considerable
slippage is required, however, for the wedge to undergo rotational movement (i.e. moment
equilibrium).

20
 1.20

1.15
Factor of safety

1.10

1.05

1.00

0.95
0.0 0.2 0.4 0.6 0.8 1.0

Lambda

Moment Force

Figure 16. Factor of safety versus lambda for a planar slip surface.

2.6.3 Composite Slip Surface

A 2D composite slip surface comprises an arc of a circle and a planar surface as illustrated in Figure 17.
The planar portion in this example follows a weak layer, which is a common situation in many
stratigraphic settings. In this case, both moment and force equilibrium are influenced by the
intercolumn shear forces. Force equilibrium factors of safety increase, while moment equilibrium factors
of safety decrease as the intercolumn shear forces increase (i.e., larger 𝜆 values).

1.25

1.20
Factor of safety

1.15

Water
1.10

1.05
0.0 0.1 0.2 0.3 0.4 0.5

Lambda

Moment Force

Figure 17. Factor of safety versus lambda for a composite slip surface.

This example illustrates that a Bishop’s Simplified analysis does not always produce a conservative
result. A rigorous formulation such as the Morgenstern-Price or Spencer method will give a lower factor
of safety than a Bishop Simplified factor of safety. This is not necessarily true for all composite slip
surfaces. For some composite slip surfaces, a mathematically more rigorous factor of safety may be
higher than the Bishop’s Simplified. It is not possible to know prior to solving when a more simplified
method will produce a lower factor of safety compared to a rigorous method. In summary, the effect of
intercolumn shear can be significant for composite slip surface as the sliding mass has the tendency to
undergo rotation and horizontal translation.

21
2.6.4 Block Specified
Figure 18 shows a 2D block-type slip surface. As with the previous composite slip surface, the moment
and force equilibrium are both influenced by the intercolumn shear. The force equilibrium is more
sensitive to the shear forces than the moment equilibrium, as indicated by the slope of the lines in
Figure 18. Significant slippage between columns is required for both horizontal translation and rotation,
giving rise to the importance of the shear forces.

1.60

1.55

1.50
Factor of safety

1.45

1.40

1.35

1.30

1.25

1.20
0.0 0.1 0.2 0.3 0.4 0.5
Lam bda

Moment Force

Figure 18. Factor of safety versus lambda for a block-specified slip surface.

2.6.5 Symmetrical Slip Surface in 3D

In 3D stability analyses, the slip surfaces will often be symmetric about the direction of movement for
geometry created by extruding a 2D cross-section. The factor of safety 𝐹 orthogonal to the direction of
movement is theoretically infinite because the driving forces and driving moments are zero.
Mathematically, this can produce noise in the x-direction and z-direction equations when the geometry
is aligned with the x-axis or z-axis. GeoStudio detects this scenario by checking if the sum of column base
normal forces and intercolumn normal forces, resolved into the direction orthogonal to movement, is
equal to zero. If symmetry is detected, with the sliding direction aligned with an axis, the solution
strategy switches to the ‘one-directional method’, in which force and moment equilibrium is calculated
only in the direction of movement, not in the x and z directions as discussed in Section 2.5.2.
Consequently, the 𝐹 versus 𝜆 plot in the direction orthogonal to sliding plot will display all 𝐹 values at a
𝜆 = 0 (e.g., of symmetry about the x-axis; Figure 19.

22
 1.25

1.2
Factor of Safety

1.15

1.1

1.05

1
-1 0 1

Z-Lam bda

Figure 19. Factor of safety versus lambda in the direction orthogonal to movement for the symmetric case.

2.6.6 Summary
Figure 20 presents the factor of safety versus 𝜆 plot for a 2D deep-seated slip surface behind a shoring
wall located in Calgary, Canada. The slip surface passes beneath the lower tip of the sheet piling. The 𝐹
versus 𝜆 curves for moment and force equilibrium are similar in this case; both are sensitive to the
intercolumn shear forces. Ignoring the intercolumn shear forces for this case results in a significant
underestimation of the factor of safety. Without including the intercolumn shear forces, the factor of
safety is less than 1.0 indicating an unstable situation. Including the shear forces increases the factor of
safety to 1.25. The difference is due to the tendency for rotational and lateral movement.

1.40

1.30
Factor of safety

1.20

1.10

1.00

0.90

0.80
0.0 0.1 0.2 0.3 0.4

Lambda

Moment Force

Figure 20. Factor of safety versus lambda for a deep-seated slip surface behind a shoring wall.

The previous examples demonstrate that the importance of the intercolumn force functions is strongly
related to the shape of the potential slip surface, which in turn is related to the tendency for the sliding
mass to undergo rotation and/or lateral movement. In 3D analysis, the intercolumn force function
becomes even more important as both the force and moment 𝐹 are sensitive to 𝜆 even for perfectly
spherical slip surfaces in a homogeneous slope. As such, a rigorous method such a Morgenstern-Price is

23
recommended because it satisfies both moment and force equilibrium subjected to an assumed
intercolumn force function.

2.7 Stress/Force Distributions

Limit equilibrium stability solutions only satisfy the requirement for equilibrium. In contrast, stress-strain
formulations also require stress-strain compatibility. The compatibility requirement can be expressed
mathematically by strain increments, which are linked to the equilibrium requirement by means of a
stress-strain constitutive law. A solution to the stress-strain equation therefore produces the stresses
within a domain. The same can not be said of the limit equilibrium solution. The computed stresses
(forces) are purely a by-product of the equilibrium requirement and the solution procedure and
therefore cannot be assumed to represent those generated by a static stress-strain formulation.

Consider the simple 45° slope shown in Figure 21 and Figure 22, where shallow and deep-seated slip
surfaces are respectively shown. The normal stress distribution along the slip surface from a limit
equilibrium Morgenstern-Price analysis with a constant intercolumn force function is compared with the
normal stress distribution from a linear-elastic finite element stress analysis. For the shallow slip surface,
the normal stresses deviate in the toe area. The normal stress distributions for the deeper slip surface
are closer, but still different for portions of the slip surface.

90

80

70

60
Normal Stress

50

40

30

20
L.E. F.E.
10

0
0 5 10 15 20 25 30
Slice number

Figure 21. Normal stress distributions along a slip surface exiting at the toe.

180

160

140
Normal stress

120

100

80

60

40
L.E. F.E.
20

0
0 5 10 15 20 25 30

Slice number

Figure 22. Normal stress distributions along a deep-seated slip surface.

Figure 23 presents the normal stress distribution for a reinforced excavation. The reinforcement loads
are applied at the point where the slip surface intersects the line of action. Again, there are significant

24
differences between the limit equilibrium normal stresses and the finite element stresses, particularly
for the columns which include the reinforcement loads.

200

150
Normal stress

100

50

0
L.E. F.E.

-50
0 5 10 15 20 25 30
Slice num ber

Figure 23. Normal stress distributions along a slip surface for a slope with reinforcement.

Figure 24 presents the case of a tie-back wall with two rows of anchors. The anchor forces are applied
where the slip surface intersects the anchor. The free body diagrams and force polygons for two
columns intersecting the upper and lower anchor are presented in Figure 25 and Figure 26, respectively.
Note that the intercolumn normal forces point away from the column on the right side. This indicates an
apparent tension between the columns, which is obviously not the case in the field. Plotting the
intercolumn forces in Figure 27 further highlights the difficulty. At the location of each anchor, the
intercolumn normal forces become negative and the intercolumn shear forces reverse direction. Of
significance, however, is the fact that the force polygons close, signifying that the columns are in
equilibrium. As such, the results fulfill the objective of the limit equilibrium formulation.

150 kN

150 kN

Figure 24. Anchor reinforced excavation.

25
 Slice 8 - GLE Method

26.512
36.004

34.093

36.257

150

38.29
139.62
81.175

Slice 13 - GLE Method

Figure 25. Free body diagram and force polygon for a column intersecting the upper anchor.

39.609
4.5301

4.2896
73.793

150

77.93

153.64 89.326

Figure 26. Free body diagram and force polygon for a column intersecting the lower anchor.

26
 40

20

0
Interslice force

-20

-40

-60

-80
0 5 10 15 20 25

Interslice number

Normal Force Shear Force

Figure 27. Intercolumn forces versus column number.

Figure 28 shows the intercolumn shear and normal forces when the anchor loads are applied on the face
of the wall. The normal force increases gradually except for the last two columns. The intercolumn shear
forces are negative and decrease gradually with column number. The negative intercolumn shear stress
reflects a negative 𝜆. Figure 29 presents the normal stress distributions along the slip surface for the two
cases. Applying the anchor load on the face of the wall yields a much smoother normal stress
distribution.

300

250

200

150
Interslice force

100

50

-50

-100

-150
0 5 10 15 20 25

Interslice number

Normal Force Shear Force

Figure 28. Intercolumn forces versus column number with anchor loads applied at the face of the wall.

27
 350
300
Normal stress - kPa

250
200
150
100
50
0
0 5 10 15 20
Slice number
On slip surface On wall

Figure 29. Column base normal stresses for the two cases.

Despite the vastly different stresses between the columns and along the slip surface, the factors of
safety are nearly identical for the two approaches of applying the anchor loads. Case 1 resulted in a
factor of safety of 1.075, while case 2 produced 1.076.

In summary, the stress distributions computed from a limit equilibrium analysis may be different from
finite element computed stresses. The finite element stresses are more realistic because the method
satisfies equilibrium while acknowledging the constitutive behaviour of the soil or rock. The implication
is that the limit equilibrium computed stresses are not representative of actual field conditions, but
rather are a result inherent to the method.

2.8 Commentary on the LE Method

Central to the limit equilibrium method is the assumption that the factor of safety is the same for all
columns. The limit equilibrium method requires an iterative technique to solve the nonlinear factor of
safety equations. In the Morgenstern-Price or Spencer methods, a second level of iterations is required
to find the column forces that result in the same moment and force factors of safety 𝐹𝑚and 𝐹𝑓.
Fundamentally, the iterations are required to meet two conditions:

 To find the forces acting on each column such that the column is in force and moment
equilibrium, and
 To find the forces acting on each column that will make the factor of safety the same for each
column.

The resulting forces/stresses of a limit equilibrium analysis do not represent the actual in situ stress
conditions because the method does not consider the constitutive behaviour of the material. Despite
this significant assumption, the method produces an acceptable factor of safety for one key reason:
global equilibrium is satisfied. Stated another way, the distribution of the local forces/stresses may be

28
incorrect, but the forces integrated over the whole are correct, local irregularities are smoothed out,
and the overall factor of safety for the entire sliding mass is correct.

Besides unrealistic stress distributions, limit equilibrium formulations have other limitations. One of the
main limitations is the difficulty with convergence under certain conditions. Most convergence problems
arise due large strength contrasts between adjacent materials, sharp corners in a slip surface, and
concentrated loads. Earth structures requiring reinforcement usually have a steep face and,
consequently, the critical slip surface position is inclined at a steep angle. The lateral forces, together
with a steep slip surface, can make it difficult to obtain a converged solution.

In addition, the critical slip surface (i.e., lowest 𝐹) is often directly adjacent to trial slip surfaces for which
a converged solution could not be computed. This casts doubt on the validity of the critical slip surface.
This problem is discussed further in the chapter on modeling reinforcement. Soil-structure interaction
also creates some difficulties in a limit equilibrium analysis. Some of these problems can be overcome
using the finite element stress-based stability analysis, or the shear strength reduction method (see the
Static Stress Strain Modeling reference book).

3 Slip Surface Searching Techniques

3.1 Introduction
Determining the position of the critical slip surface with the lowest factor of safety remains one of the
key issues in a stability analysis. During the solution process, a potential slip surface is generated, and
the associated factor of safety is calculated. This is repeated for many possible slip surfaces and the trial
slip surface with the lowest factor of safety is deemed to be the critical slip surface.

There are many approaches for defining the shape and positions of trial slip surfaces. This chapter
explains the procedures available in SLOPE/W and SLOPE3D and discusses the applicability of the
methods to various situations. Finding the critical slip surface requires considerable guidance from the
analyst despite the advanced capabilities of the software. The soil stratigraphy may influence the critical
mode of failure and the stratigraphy therefore must be considered in the selected shape of the trial slip
surfaces. In the case of a tie-back wall, it may be necessary to consider separately a toe failure and a
deep-seated failure. In an open pit mine, the issue may be stability of a slope bench versus overall high
wall stability. Generally, not all potential modes of failure can be investigated in a single analysis. In such
cases, the slip surface search needs to be specified and controlled by the analyst.

A general procedure for defining trial slips may result in some physically inadmissible trial slip surfaces;
that is, the trial slip surface has a shape that cannot exist in the physical system. Often it is not possible
to compute a factor of safety for such unrealistic situations, due to a lack of convergence. Sometimes,
however, safety factors can be computed for unrealistic slips, and then it is the responsibility of the
analyst to judge the validity of the computed factor of safety. The software cannot make this judgment.

3.2 Grid and Radius

The grid and radius technique is used to generate circular slip surfaces in a 2D analysis. Circular trial slip
surfaces were inherent in the earliest limit equilibrium formulations and the techniques of specifying

29
circular slip surfaces has become entrenched in these types of analyses. The trial slip surface is an arc of
circle. The arc is that portion of a circle that cuts through the slope. A circle can be defined by specifying
the 𝑥 -y coordinate of the centre and the radius. A wide variation of trial slip surfaces can be specified
with a defined grid of circle centers and a range of defined radii.

Figure 30 presents the configuration for a typical problem. The grid above the slope defines the center
of rotation for each trial slip surface. In this example there are thirty-six (6 x 6) grid points or circle
centers. The grid is defined by three points; the upper left (18), lower left (16) and lower right (12). The
radii for each centre of rotation are specified with the tangent lines. The lines are specified by the four
corners of a box. The four corners in Figure 30 are 19 (upper left), 21 (lower left), 22 (lower right) and 20
(upper right). The four points need to be drawn by starting at the upper left and proceeding in a
counter-clockwise direction around the box. The number of increments between the upper and lower
corners can be specified. In the above example there are five increments making the total number of
radius tangent lines equal to six.

18

16 12

19 20

21 22

Figure 30. The grid and radius method of specifying trial slip surfaces.

Trial slip surfaces are generated by first calculating the perpendicular distance between the radius line
and the grid centre (Figure 31). The perpendicular distance becomes the radius of the trial slip surface.
The specified radius lines are more correctly tangent lines; that is, lines that are tangent to the trial
circles. The trial slip surface is drawn and the portion of the arc below the ground surface is used to
define the sliding mass. In this case, the factor of safety will be calculated for 216 (36 x 6) trial slip
surfaces.

30
 18

16 12

Radius

19 20

21 22

Figure 31. Radius for a single trial slip surface using grid and radius.

The radius lines (points 19, 21, 22, 20) can be located at any convenient position and can form any
quadrilateral shape. The illustration in Figure 32 is entirely acceptable. It should also be noted that the
tangent lines are extended mathematically to the edge of the domain to ensure combability with the
grid. As such, the two configurations shown in Figure 33 are equivalent.

Figure 32. Specification of radius lines.

Figure 33. Automatic extension of the radius tangent lines.

31
It is often convenient to force all slip surface to pass through a single point to explore a specific mode of
failure (Figure 34). This is accomplished by collapsing the radius lines to a single point (click four times in
a single location). Similarly, grid of centers can be collapsed to a single point. This makes it possible to
analyze just one slip surface, which can be very useful for doing comparisons of various features or
options.

Radius point

Figure 34. Radius lines collapsed to a single point.

3.3 Entry and Exit

The Entry and Exit technique can be used in both 2D and 3D stability analyses. The main advantage of
this approach is the ability to visualize the extents and/or range of trial slip surfaces prior to solving. As
such, the analyst can easily isolate various modes of failure by re-positioning the entry and exit zones.
The following sections discuss the technique for circular slip surfaces in 2D and ellipsoidal slip surfaces in
3D.

3.3.1 Circular Slip Surfaces

Figure 35 presents an entry-exit zone for a 2D analysis. Each zone is drawn along the ground surface and
then sub-divided into specified increments. Each point in the entry zone is connected to a point along
the exit area to form a line, which is bisected by an orthogonal line. Radius points are created to form
the required third point of a circle. This radius point is used together with the entry-exit points to form
the equation of a circle. Controls are imposed on the locations of these radius points so that the circle
will not be a straight line (infinite radius), and the entry angle of the slip circle on the crest will not be
larger than 90 degrees (undercutting slip circle). The equation of a circle gives the center and radius of
the circle, the trial slip surface is then handled in the same manner as the conventional Grid and Radius
method.

32
 Entry points

Exit points
Radius points

Figure 35. Entry and exit zones for generating circular slip surfaces.

Figure 36 shows the valid slip surfaces when the entry increments, exit increments, and radius
increments are set equal to 5. A total of 216 (6 x 6 x 6) slip surfaces are generated. The entry-exit
method can also be combined with radius tangent lines to create control on the search zone (Figure 37),
which is useful for exploring different modes of failure. For example, the radius tangent lines can be
positioned lower in the soil profile to drive a deep-seated slip surface search that interacts with a weak
layer or impenetrable material surface.

Figure 36. Trial slip surfaces for an entry-exit search.

33
Figure 37. Entry-exit combined with radius tangent lines.

3.3.2 Ellipsoidal Slip Surfaces

An ellipsoid is a 3D surface defined by:
Equation 32

( ) ( ) ( )
𝑥' ‒ 𝑥0
𝑟
𝑥'
2
+
𝑦' ‒ 𝑦0
𝑟
𝑦'
2
+
𝑧' ‒ 𝑧0
𝑟
𝑧'
𝑛
=1

𝑟' 𝑟 ' 𝑟'

where (𝑥0, 𝑦0, 𝑧0) are the global coordinates of the centroid, 𝑥 , 𝑦 and 𝑧 are ellipsoid radii in 3
𝑟 ' = (𝑘 )𝑟 '
directions, and 𝑛 is the curvature exponent. The radius 𝑧 𝑥 , where 𝑘 is the aspect ratio. Setting
𝑘 = 2 and 𝑛 = 2 generates a standard ellipsoid (Figure 38). Setting 𝑘 = 1 with 𝑛 = 2 produces a sphere.
The exponent 𝑛 can be used to alter the tips of the ellipsoid (Figure 39).

Figure 38. Standard ellipsoid with 𝑛 = 2.

34
Figure 39. Ellipsoid with 𝑛 = 2 (left) and 𝑛 = 6 (right).

The entry and exit zones of a 3D analysis are defined as grids in plan view that are mapped down onto
the ground surface (Figure 40a). Each point in the entry grid is connected to a point in the exit grid and a
third radius point is generated along an orthogonal bisector line (Section 3.3.1). The coordinates of the 3
points are used to calculate the radii 𝑟 = 𝑟𝑥' = 𝑟𝑦' and the coordinates of the centroid (𝑥0,𝑦𝑜, 𝑧𝑜).

(a) (b)
Figure 40. Entry and Exit method in 3D stability analysis (a) 3D view (b) 3D ellipsoid in plan view.

35
3.4 Fully Specified

3.4.1 2D Definition
A trial slip surface can be specified with a series of data points in a 2D analysis. This allows for complete
flexibly in the position and shape of the slip surface (Figure 41). Note that the specified surface should
ideally start and end outside the geometry such that the intersection points with the ground surface can
be determined automatically by the software.

Computed intersection point

Computed intersection point

Specified points

Figure 41. Fully specified slip surface.

An Axis Point is then defined for the moment equilibrium calculations (Figure 42). In general, the Axis
Point should be located at the approximate center of rotation of all the specified slip surfaces. If the
point is not defined, the software automatically calculates one based on the lateral extents of the slip
surface. The factor of safety calculation is not sensitive to the position of the Axis point for methods that
satisfy both moment and force equilibrium (e.g., Spencer and Morgenstern-Price methods). However,
for simplified methods (e.g., Ordinary and Simplified Bishop), the factor of safety calculations can be
sensitive to the position of the Axis Point.
1.583

Computed intersection point Axis Point

Computed intersection point

Figure 42. Definition of an axis point for moment equilibrium.

The Fully Specified method is useful when the position of large portions of the slip surface is known from
slope inclinometer field measurements, geological stratigraphic units, and/or surface observations. The
option may also be useful for other cases such as the sliding stability of a gravity retaining wall (Figure

36
43). While the Fully Specified method is completely flexible with respect to trial slip surface shapes and
positions, it is limited in that every trial slip surface needs to be specified. The method is consequently
not suitable for doing many trials to find the critical slip surface.
6
23 22 21 20 19 18 17 16 15 14
1 2 13

Soil: 1
Soil: 2 Retaining Wall
Backfill

5 3 10 9 4
7 8
Soil: 3
11 Foundation Clay 12

Figure 43. Stability analysis of a gravity retaining wall.

3.4.2 3D Definition
A fully specified ellipsoid is parameterized by specifying the coordinates of the centroid (𝑥0, 𝑦0, 𝑧0)
𝑟' 𝑟 ' 𝑟'
ellipsoid radii in 3 directions ( 𝑥 , 𝑦 , 𝑧 ), and the curvature exponent 𝑛 (Equation 32). The ellipsoid can
also be rotated through three successive simple rotations by specifying the rotation angles (𝜃𝑥, 𝜃𝑦, 𝜃𝑧)
and the rotation sequence. The initial rotation can be taken about any of the three axes while successive
rotations can be taken about the other axes, including a repeated axis. Figure 44 shows the 𝑥'𝑧' planar
cross-section of an ellipsoid passing through various stages of a (𝑧' ‒ 𝑦' ‒ 𝑧') sequence of rotations. The
sequence starts by rotating the initial system of axes counterclockwise about the 𝑧’-axis by an angle 𝜙.
The axes are then rotated counterclockwise by an angle 𝜃 about the new 𝑦'-axis. Finally, the axes are
rotated counterclockwise about the new 𝑧’-axis by an angle 𝜓.

Figure 44. Rotation of the planar cross-section of an ellipsoid through a z-y-z sequence.

The ability to specify any sequence of rotations can result in illogical positing of the fully specified
ellipsoid; however, the generalization was necessary to facilitate transformation between different
coordinate systems. Figure 45 provides an example where the source (subscript 𝑠) East-North-Upwards
coordinate system has axis 𝑋𝑠, 𝑌𝑠, and 𝑍𝑠, respectively. The source axes map to the 𝑋, ‒ 𝑍, and 𝑌 axes,

37
respectively, in GeoStudio. A source axis rotation order of 𝑍𝑠-𝑌𝑠-𝑋𝑠 with angles 𝛾, 𝛽, and 𝛼, respectively,
would be specified in GeoStudio with an order of 𝑌-𝑍-𝑋 and angles 𝜃𝑦 = 𝛾, 𝜃𝑧 =‒ 𝛽, and 𝜃𝑥 = 𝛼,
respectively. Rotations are counterclockwise positive, as indicated by the rotation vectors, so 𝜃𝑧 =‒ 𝛽
rotates the navigation cube clockwise about the 𝑍 axis (Figure 45b).

(a) (b)
Figure 45. Transformation from a source coordinate system (i.e., a) to the GeoStudio coordinate system (i.e., b).

3.5 Block Specified

The Block Specified search technique is used to generate piece-wise linear slip surfaces in a 2D analysis.
Block shaped analyses can be performed by specifying two grids as shown in Figure 46. The grids are
referred to as the left block and the right block. Each grid is defined with an upper left point, a lower left
point, and a lower right point. Arrows are drawn at the upper left and right corners to graphically
portray the range of specified angles for the projection of the slip surface. The user-input projection
increments further subdivide the projection angles.

1 2
14

4
8 11

9 10 12 13

Figure 46. Block Specified slip surface definition.

The slip surface consists of three-line segments. The middle segment goes from each grid point on the
left to each grid point on the right. The other two segments are projections to the ground surface at a
range of specified angles. Figure 47 presents a single slip surface generated between two points for the
Block Specified search technique.

38
Figure 47. A single slip surface generated by Block Specified search.

Selection of the left and right projection angles (Figure 48) is often done based on earth pressure theory.
The situation in the toe area is similar to a passive earth pressure condition where the sliding mass is
being pushed outward and upward. In the crest area, the situation is analogous to active earth pressure
conditions. From lateral earth pressure considerations, the passive (toe) slip surface rises at an angle
' '
equal to (45 - 𝜙 /2), while the active slip line dips at an angle of (45+ 𝜙 /2). These considerations can be
used to guide the selection of the projection angles.

Figure 48. Projection angles for the Block Specified method.

Angles at the base of a 2D column are defined in a counterclockwise direction from the positive x-
coordinate axis. An angle of zero means a horizontal direction to the right, an angle of 90 degrees means
an upward vertical direction; an angle of 180 degrees means a horizontal direction in the negative x-
coordinate direction, and so forth. In the above example, the right toe (passive) projection angles vary
between 30 and 45 degrees, and the left crest (active) projection angles vary between 115 and 130
degrees.

The Block method also needs a defined Axis about which to take moments like the Fully Specified
method. If the Axis point is not defined, SLOPE/W will compute an Axis location based on the geometry
of the problem and the positions of the left and right blocks.

The Block method can often produce slip surfaces that fail to converge due to sharp corners (Figure 49).
This is particularly problematic when the grid blocks are close to each other. The Block method is most
suitable for cases where there is a significant distance between the two blocks.

39
Figure 49. Trial slip surface with a sharp corner.

Worth noting is that the two grid blocks can be collapsed to a line with points or to a single point. If the
two left specified points in the grid block are the same, the block will collapse to a line. If all three points
are the same, the grid block will collapse to a single point.

There are situations where it is preferable to have all the middle line segments of the trial slip surface
parallel. Take, for example, the case of a slope where the material is strongly bedded and the strength
along the bedding is less than the intact material (Figure 50). The grid blocks are placed so that the bases
are parallel to the bedding. By selecting the “No crossing” option, the middle segments of the trial slip
surfaces will not cross the bedding. This approach can be combined with the Anisotropic or Compound
Strength material models to ensure that the strength reduces when the slip surface is parallel to the
bedding.

Direction of bedding

Figure 50. Block search method in a slope with bedding.

3.6 Cuckoo Search

The Cuckoo search technique is used to generate circular and ellipsoidal slip surfaces in 2D and 3D
analyses, respectively. Contrary to Entry and Exit technique, where the parameters defining the slip
surface circle or ellipsoid are either prescribed or back calculated from the geometry, the Cuckoo search
randomly generates slip surfaces by altering the parameters that define the equation of the slip surface.
The Cuckoo technique is a nature-inspired metaheuristic optimization algorithm developed by Yang and
Deb (2009). It is based on the brood parasitic behavior of some Cuckoo species, which lay their eggs in
the nests of host birds of other species.

40
In the Cuckoo Search algorithm, each nest contains one slip surface (i.e., one Cuckoo’s egg) and the
Cuckoo routine randomly generates the first set of slip surfaces (i.e., the first generation of eggs). An
idealized model of Cuckoo breeding methods is then applied to develop the next generations.
Mathematically, it is a global optimization search technique for a mathematical function (i.e., an
‘objective function’ that attempts to find the global minimum, while avoiding local minima).

In a slope stability analysis, the objective function is ultimately the factor of safety plotted against the
parameters that define the slip surface. The final goal of the optimization algorithm is to minimize the
objective function and find the most critical slip surface. The objective function could be described as
follows:

 2D Circular: 𝐹 as a function of the center of rotation coordinates and radius of the circle.
 3D Ellipsoid: 𝐹 as a function of the coordinates of the center of symmetry (𝑥0, 𝑦0, 𝑧0), rotation
about the y-axis, and the radii defining the ellipsoid (𝑟𝑥, 𝑟𝑦, 𝑟𝑧).

The Cuckoo algorithm naturally includes a method that randomly generates new ‘nests’ of guesses that
drive the process toward a solution and/or replaces abandoned slip surface guesses. The procedure is
enhanced by the Lévy Flights random walk procedure.

In a 3D analysis, the maximum extents of the search zone can be defined when using the Cuckoo search
(Figure 51). Although this could encompass the entire domain, it is more computationally efficient and
best-practice to focus the search in certain areas based on engineering judgement. Any slip surface that
exits outside the zone is discarded from the analysis.

Figure 51. Defining the search zone extents for Cuckoo.

The total number of trial slip surfaces in Cuckoo search is the product of the following two Cuckoo
parameters:

 Number of nests: the number of nests where slip surfaces are produced.
 Number of iterations: The number of times a new generation of nests is reproduced.

41
The default parameters are generally acceptable for most problems but may need to be adjusted if users
deem that the generated slip surfaces do not adequately cover the search area. The distribution of slip
surfaces over the domain can be visualized via a slip surface color map. Since local search involves
random permutations of nest hosts, it is suggested to keep the default number of nests.

3.7 Non-Circular Slip Surfaces

The geology has a major influence on the location and shape of the critical slip surface. Multiple layers
with varying strength and stiffness, bedding layers, and joints/fractures often impose control on the
deformation pattern that develops. This section discusses various options for generating slip surfaces
that adhere to structural and/or geological controls.

3.7.1 Impenetrable Material Model

The impenetrable material model is used to force a composite slip surface shape in either a 2D or 3D
stability analysis. An ‘impenetrable’ material is assigned to all regions (2D) or solids (3D) beneath the
interface with the high strength material. After the trial slip surface is drawn, the intersection with the
impenetrable surface is determined. All columns intersecting the impenetrable surface are truncated
such that the slip surface follows the impenetrable contact (Figure 52). The advantage of this approach
is that the impenetrable material model can be used to capture any irregular contact (Figure 53).

Impenetrable (bedrock)

Figure 52. Composite slip surface controlled by an impenetrable layer.

Impenetrable (bedrock)

Figure 53. Irregularly shaped impenetrable layer.

The impenetrable material model can also be combined with geometry modifications and/or other
material models to capture weak layers within the stratigraphy. Consider the horizontal and sloping

42
weak layer configuration in Figure 54 and Figure 55, respectively. In both cases, a thin region was
created to model the weak layer. That stated, using a Weak Layer may be a better approach for these
scenarios (Section 3.7.2).

Figure 54. Horizontal weak layer.

Cover material

Impenetrable

Figure 55. Potential sliding at the interface of a geosynthetic liner.

3.7.2 Weak Surfaces

Weak surfaces can be used in 2D and 3D analyses to represent weak zones in soil or rock. A weak
surface could be used, for example, to represent a single a discontinuity, a shear zone at residual
strength, or an interface along a geomembrane. The use of a single weak surface is like using the
impenetrable material model (Section 3.7.1); however, it differs in that a unique material is associated
with the surface. Furthermore, multiple weak surfaces can be defined in a single analysis. Weak surfaces
are defined using a polyline in 2D and background mesh in 3D (Figure 57) and a material is associated
with each weak surface.

43
Figure 56. Example of a weak layer in a 3D stability analysis defined using a background mesh.

Trial slip surfaces that intersect a weak surface are truncated (Figure 57a). The top-most weak surface
governs (Figure 57b). The strength definition of columns intersected by the weak surface is defined by
the material associated with the weak surface. Weak surfaces should ideally extend across the domain
to mitigate the generation of poorly shaped slip surfaces (Figure 57c). Any slip surface in 2D or 3D with a
stairstep shape is deemed inadmissible.

Figure 57. Cross-section showing weak surfaces: a) single weak surface, b) two weak surfaces with different material
associations, and c) a discontinuous weak surface.

3.8 Anisotropic Surfaces

Figure 58 presents a hillslope with a folded geological unit containing bedding anisotropy. Anisotropic
surfaces can be used in 2D and 3D analyses to generate directionally dependent shear strength. Soil and
rock exhibit anisotropic strengths due to the presence of discontinuities, where the strength of the
discontinuities varies considerably from the surrounding (intact) material. Discontinuities are typically
the result of bedding, schistosity, joints, foliation, cleavage, fractures, or faults. The strength anisotropy

44
may occur ubiquitously in a geological unit (i.e., a discontinuity set) or as a single discontinuity. In
contrast to a Weak Layer, an Anisotropic Surface does not impose a direct control on the shape of a slip
surface. Anisotropic surfaces are defined using a polyline in 2D and a background mesh in 3D.

Figure 58. Anisotropic surface, search options, and detail for the angle of the anisotropy (𝛼) obtained by closest distance.

The input properties of an Anisotropic Surface, which are like those of the Compound Strength material
(Section 4.11), include: a base material representing the intact soil/rock matrix, a material representing
the discontinuity, and two angle ranges, A and B. The A and B parameters associated with each
anisotropic surface are angular tolerances around the discontinuity that control the transition of shear
strength from the anisotropic material to the intact material (Figure 59). The solver calculates the shear
strength based on the angle between the column base and the anisotropic feature (e.g., joint plane).
Within the range of A, the discontinuity shear strength is used. Outside the range of B, the shear
strength of the base material is used. For angles between A and B, the shear strength is computed based
on a linear interpolation between the shear strengths of the base and the discontinuity materials.

Figure 59. Definition of the A and B angle ranges.

45
The base angle of every column (𝜃) is known. The angle of the anisotropy (𝛼) corresponding to each
column base needs to be determined from a point on the anisotropic surface. There are three options to
locate the point on the anisotropic surface: 1) closest distance; 2) vertical; and 3) horizontal (Figure 58).
The horizontal option is only relevant in 2D. The closest distance option is shown in the inset of Figure
58. The strength at the base of the column 𝑆 is determined from the strength of the intact material 𝑆'
and/or anisotropy material 𝑆" as follows:

𝑖𝑓(|𝜃 ‒ 𝛼| ≤ 𝐴) {𝑆 = 𝑆"}
Equation 33
𝑒𝑙𝑠𝑒 𝑖𝑓(𝐴 ≤ |𝜃 ‒ 𝛼| ≤ 𝐵) {𝑆 = 𝑆'' + (𝑆' ‒ 𝑆")𝑅}
𝑒𝑙𝑠𝑒 {𝑆 = 𝑆'}
where

|𝜃 ‒ 𝛼 | ‒ 𝐴 Equation 34
𝑅=
𝐵‒𝐴
and angular difference |𝜃 ‒ 𝛼| is calculated as

𝑎∙𝑏 Equation 35
|𝜃 ‒ 𝛼| = acos
|𝑎||𝑏|
if the angular difference |𝜃 ‒ 𝛼| ≤ 𝜋 2 and otherwise as:

𝑎∙𝑏 Equation 36
|𝜃 ‒ 𝛼| = 𝜋 ‒ acos
|𝑎||𝑏|

3.9 Optimization
Physical observations of failed slopes demonstrate that the critical slip surface is often not perfectly
circular/elliptical (2D) or spherical/ellipsoidal (3D). The actual shape of the slip surface is governed by
several factors including spatial variations in strength and pore-water pressures, strength anisotropy,
and topographic variations. Optimization is used to search for a modified slip surface shape that
produces a lower factor of safety and a mode of failure that better honors the physical system (Figure 60
and Figure 61).

Figure 60. Circular (left) and optimized (right) critical slip surfaces.

46
Figure 61. Block specified (top) and optimized (bottom) slip surfaces.

3.9.1 2D Optimization
GeoStudio uses a modified Monte-Carlo random-walk procedure for 2D optimization (Greco, 1996;
Malkawi, Hassan and Sarma, 2001). The optimization procedure initiates by creating a piece wise linear
representation of the slip surface based on the specified number of starting points (Figure 62a). The
optimization algorithm proceeds as follows:

1. Calculate allowable x and y offsets for each point (Figure 62b).

2. Calculate a random set of coordinates for the current point within its allowable range.
3. Generate trial slip surfaces from the sets of coordinates (Figure 62c). Each slip surface is one
(optimization) iteration.
4. Calculate the factor of safety for each slip surface and identify the most critical.
5. Insert a point in the longest segment of the most critical slip surface if a) the current point is the
last point; and b) the number of inserted points has not exceeded the specified maximum
(Figure 62 d).
6. Advance to the next point and return to Step 1 until the maximum number of (optimization)
iterations has been reached.

47
Figure 62. 2D Optimization procedure: a) initial fit to critical slip surface; b) mutation boundaries for each point; c) trial slip
surfaces; and d) new critical slip surface with a point inserted.

Step 1 of the optimization procedure limits the mutation range of each point; however, this does not
ensure that all trial slip surfaces are kinematically admissible; consequently, there are user specified
angle limits on the concavity (angle) between columns. Trial slip surfaces that violate the concave angle
limits are deemed invalid (Figure 63; CCW+ angle):

𝑖𝑓(𝜃𝑛 ‒ 𝛼 > 𝜃𝑛 + 1) 𝑡ℎ𝑒𝑛 𝐼𝑁𝑉𝐴𝐿𝐼𝐷 Equation 37

where 𝛼 is the concavity angle tolerance on the resisting or driving side of the slip surface.

Figure 63. Concavity between columns exceeding the tolerance.

48
Optimized slip surfaces can still have sweeping concavities over larger distances that do not have
column-to-column concavity that violates the angle tolerances. These types of slip surface are often
mathematically converged but may be kinematically inadmissible. The Factor of Safety versus Lambda
convergence plots may indicate that convergence is tenuous or fortuitous; regardless, the optimized
solution must be validated against the physical reality.

The effectiveness of optimization is sometimes dependent on the starting position of the critical slip
surface because of the offset restrictions imposed in Step 1. Figure 64 shows the effect of optimization
on two Fully Specified slip surfaces passing through a domain comprising a thin weak layer. The
optimized slip surface only passes through the weak layer for the case where the fully specified slip
surface was within proximity at the onset.

Figure 64. Optimization of two Fully Specified slip surfaces.

Figure 65 shows a graph of Factor of Safety versus (optimization) iteration. This graph type can be used
to discern if the maximum number of optimization iterations needs to be increased or, conversely, if it
can be decreased to reduce solve time.

49
 Factor of Safety vs. Iteration Number

1.46

1.44

1.42
Factor of Safety

1.4

Slip 0

1.38

1.36

1.34

1.32
0 500 1,000 1,500 2,000

Iteration

Figure 65. Factor of safety versus optimization iteration.

3.9.2 3D Optimization
The optimization procedure initiates by fitting a Non-Uniform Rational Basis Spline (NURBS) to a
spherical or elliptical slip surface. The fitting procedure involves (Figure 66):

1. Calculating the extents of a bounding box over the slip surface in plan view.
2. Calculating the X and Z coordinates of the control points from the specified number of points
and bounding box extents.
3. Calculating the Y coordinates of the control points to best fit the NURBS to the original slip
surface.

50
Figure 66. Initial control grid before calculating the y-coordinate of the control points to fit the NURBS to the slip surface.

The slip surface shape optimization employs a Cuckoo Optimization Algorithm to modify the coordinates
of the NURBS control points (Figure 67). The optimization algorithm comprises the following stages:

1. The control points are uniformly translated in three directions and the grid is scaled in the XZ
plane. Each translation/scaling is one optimization iteration. The number of optimization
iterations completed during this stage is 𝑖1 = 𝑔 ∗ 𝑛, where 𝑔 is the number of generations and 𝑛
the number of nests. The control points that produce the lowest factor of safety are preserved
for the next stage of optimization.
2. The control points are variably translated in the Y direction. Each generation of y-coordinates is
one optimization iteration. The number of optimization iterations completed during this stage is
𝑖2 = 𝑔 ∗ 𝑛, where 𝑔 is the number of generations and 𝑛 the number of nests. The control points
that produce the lowest factor of safety are preserved for the next stage of optimization.
3. The control point coordinate histories are cleared, and Stage 2 is repeated a fixed number of
times (ℎ).
4. All stages of the optimization algorithm are repeated until the total number of iterations equals

the maximum iterations: ∑𝑖1 + ∑𝑖2ℎ = 𝑖𝑚𝑎𝑥 .

51
Figure 67. Section view of the control polygon and optimized slip surface (NURBS).

As with 2D optimization, the slip surface shapes generated by the optimization algorithm can sometimes
be mathematically converged but kinematically inadmissible. The control points defining the NURBS are
not located on the slip surface except by coincidence (Figure 67), making it harder to mitigate concavity.
Restrictions are placed on the relative vertical positioning of the control points, regardless, sweeping
concavities can still occur, making it necessary to validate the optimized slip surface against the physical
reality.

3.10 Effect of Strength on Critical Slip Surface Location

3.10.1 Frictional vs Undrained Strength

The shear strength of the soil often poses the greatest control on the location of the critical slip surface.
The two common cases that may present a challenge include: a) a frictional soil/rock (no cohesion) in
which the angle of repose approaches the friction angle; and, b) a slope comprising a single
homogenous material defined using undrained strength. In the first case, the critical slip surface will
always tend towards the infinite slope case where the factor of safety is given by:

𝑡𝑎𝑛𝜙' Equation 38
𝐹=
𝑡𝑎𝑛𝛼

where 𝛼 the inclination of the slope. The critical slip surface is parallel and immediately next to the slope
face as shown Figure 68 for a 2D analysis.

52
 2
1

1
Figure 68. Shallow critical slip surface for friction-only case.

The tendency to move towards the infinite slope case means the radius of the circle tends towards
infinity. The minimum factor of safety is therefore usually on the edge of the grid of rotation centers
(Figure 69). Making the grid larger does not resolve the problem. The Grid and Radius method appears
to break down under these conditions, and the concept that the minimum factor of safety should be
inside the grid is not achievable.

1.203
1.24
0
1.2
80

Figure 69. Center of rotation for the infinite slope case.

The opposite occurs when the soil strength is defined by a constant undrained strength; that is, the
friction angle is zero. In this case, the critical slip surface will tend to go as deep as possible as shown in
Figure 70. In this example the depth is controlled by the geometry. If the problem geometry was to be
extended, the critical slip surface would go even deeper.

53
Figure 70. Deep slip surface for homogeneous undrained case.

Figure 71 shows the factor of safety contour plot on the search grid. The minimum factor of safety is
always on the lower edge of the search grid. Once again, it is not possible to define a minimum center
inside the search grid for this homogenous undrained case.
1.200

1.300

1.088

Figure 71. Center of rotation for the undrained strength case.

Case 1 and 2 can be summarized as follows:

 The cohesion or ‘apparent cohesion’ is rarely zero near the ground surface due to desiccation or
the presence of a root zone. More importantly, the near surface soil/rock is usually within the
unsaturated zone, where the matric suction generates additional effective stress and therefore
additional strength.
 The undrained strength is rarely constant within a soil layer, even for very soft soils. If the
undrained strength is modeled more realistically with some increase with depth, the critical slip
surface position no longer tends to go as deep.

54
3.10.2 Minimum Depth Control
The minimum depth input can be used to filter out very shallow slip surface in both 2D and 3D analysis.
However, the critical slip surface will still tend towards the minimum for Case 1 unless the fundamental
issue is addressed.

3.10.3 Best Practice

The most realistic position of the critical slip surface is computed when effective strength parameters
are used along with a realistic pore water pressure definition. In fact, it could be argued that the number
one issue in most slope stability problems is pore water pressure, which naturally effects the effective
stress shear strength. Effective strength parameters can be readily defined with considerable accuracy
for most soils and rocks. Similarly, a reasonable definition of pore water pressure conditions be achieved
using a piezometric surface(s) for simple groundwater flow systems, or by using a finite element
seepage analysis for more complicated systems.

3.11 Tension Cracks and Exit Projections

The steepness of the of the slip surface near entry into the slope can have some adverse numerical
consequences. At the crest, the normal force at the base of the first column may point away from the
column, indicating the presence of tension instead of compression. Generally, this is considered
unrealistic, particularly for materials with little or no cohesion. Physically, it may suggest the presence of
a tension crack. In the exit area, steep slip surface angles can sometimes cause convergence difficulties.
These two issues can be handled using a tension crack zone and/or exit projection, respectively.

A tension crack can be defined for a 2D stability analysis using two options: a) tension crack angle; or, b)
tension crack line. In the case of a tension crack angle, any column in the active zone with a base angle
that exceeds the input value is discarded (Figure 72). The active earth pressure angle can be used as a
guide to determine this value. In the second option, a line representing the base of a tension crack is
drawn on the domain. All columns above this line are discarded. Finally, it is noted that the tension crack
can be assumed to be filled with water by specifying a value between 0 and 1 (i.e. 0 to 100% full). A
lateral load is then applied to the last column based on the depth of water and specified water unit
weight.
Tension crack

Figure 72. Example of tension crack specified via an angle.

The problem of a steep exit angle can be overcome by specifying a projection angle. If the column base
has an inclination greater than the specified angle, the trial slip surface is projected at the specified
angle (Figure 73). The passive earth pressure angle can be used as a guide to determine this value.

55
 Footing Load

Exit controlled by
specified angle

Figure 73. Exit angle controlled by specified angle.

3.12 Results Interpretation

The slope stability examples available on the website should be consulted for the best practice in
interpreting a stability analysis. The following features are noted:

 Convergence graphs: The 𝐹 versus 𝜆 plot(s) should always be reviewed to ensure that the cross-
over point of the moment and force equilibrium lines are unambiguous. For 3D, a convergence
graph if available for both the x and z directions.
 Graphing: To effectively interpret and understand the analysis results, plot graphs of
intercolumn data including normal stress, effective stress, pore water pressure, shear strength,
shear mobilized, strength parameters, and much more.
 Grid and Radius Contours (2D): Contour the 𝐹 in the grid to review the values around the critical
slip surface. Note that it is not always possible to have the critical slip surface centre in the
middle of the search grid. A more important consideration is that the correct mode of failure is
considered.
 Slip Surface Color Map with Filter: Use the slip surface color map to visualize all slip surfaces.
The filter in the slip surface docking window can be used to explore valid slip surface, invalid slip
surfaces, all slip surfaces within a range of factor of safety and more. This type of review builds
confidence in the analysis set-up and results.
 Safety Map: the slip surface safety map can be used to present all slip surfaces within a range of
factor of safety, helping to elucidate the mode of failure.
 Slip Surface Selection Tool: Manually select slip surfaces and review output data.
 View Column Information: Select columns to review the free body diagram and/or data.

3.13 Commentary on Physical Admissibility

The procedures used to generate trial slip surfaces can often produce slip surfaces that are not
physically admissible; that is, they cannot exist in reality, or movement cannot possibly take place along
the trial slip surface. Fortunately, in many cases it is not possible to obtain a solution for a physically
inadmissible trial slip surface due to lack of convergence. However, sometimes factors of safety are
computed and displayed for cases where it is highly improbable that the potential sliding mass can slide
or rotate.

56
This issue is complicated by the fact that there are no firm criteria to mathematically judge physical
inadmissibility. There are some general rules that prevent some unrealistic cases like, for example, the
inclination of the slip surface exit angle. If the slip surface exit is too steep, the results can become
unrealistic and so the exit angle is arbitrarily limited to a maximum value. It is, however, not possible to
develop similar rules for all potentially inadmissible cases. Consequently, it is necessary for the analyst
to make this judgment and ensure that the critical slip surface shape is in-keeping with the physical
system.

3.14 Error Message for Slip Surfaces

A typical analysis may involve many trial slip surfaces; however, some of the slip surfaces may not have a
valid solution. In such cases, a factor of safety with a value ranging from 983 - 999 is stored. For
example, a value of 999 is stored for slip surfaces that failed to converge. The following summarizes the
various error conditions:

983 – Slip Surface encounters missing stresses at the column base center in a stress-based stability
analysis. This happens when finite element stresses are used in the stability analysis and if the slip
surface is passes through a region without stresses defined.

984 – Slip Surface encounters an air gap or vertical/overhanging bedrock region (Figure 74). This
happens when adjacent regions are not connected properly or when a slip surface passes through a
vertical or overhanging bedrock, causing the slip surface to be discontinuous.

Bedrock

Figure 74. Slip surface intersecting an overhanging or vertical bedrock.

990 - Slip surface does not intersect the ground surface. All valid slip surfaces must enter and exit the
ground surface.

991 - Slip surface center defined below the slip surface exit points. This results in a column with base
angle larger than 90 degrees.

992 - Slip surface cannot be generated (Figure 75). This error covers a multitude of problems including a
slip surface that does not enter or exit on the ground surface and overhanding slip surfaces.

57
Figure 75. Slip surface cannot be generated.

993 - Slip surface is too shallow (Figure 4 48). This happens when the thickness of all columns are smaller
than the specified minimum slip surfaces thickness.

994 - No intersecting point is obtained in the factor of safety versus lambda plot (Figure 76). This
happens when the problem fails to converge in the lambda loop.

Factor of Safety vs. Lambda

1.4

1.3
Moment

1.2
Factor of Safety

1.1

1.0

0.9 Force

0.8
-1.5 -1.0 -0.5 0.0 0.5

Lambda

Figure 76. Example of a problem that failed to converge.

995 - Slip surface could not be analyzed. This happens when all other conditions appear to be fine, but
the solver encounters problems in the factor of safety computation. For example, the mobilized shear
resistance is in the opposite direction due to a large negative normal force. This error message could
also be related to extreme soil properties or large external forces relative to the sliding mass.

996 - Slip surface is not in the same direction as the specified direction of movement. This happens
when the direction of movement is not set properly in SLOPE/W or in a complex geometry where a slip
surface has an exit point higher than the entry point.

998 - Slip surface enters or exits beyond the slip surface limit. This happens when the slip surface
extends beyond the specified slip surface limits.

999 - Slip surface does not have a converged solution. This happens when the procedure fails to solve
the nonlinear factor of safety equation until the factor of safety converges within a specified tolerance.

58
4 Material Models

4.1 Introduction
The shear strength at the base of each column is central to the limit equilibrium formulation and
comprises three components: cohesion, frictional resistance, and unsaturated strength (Equation 3).
Referring to the force and moment 𝐹 equations in Section 2.4, either the cohesion and/or frictional
shear strength must be defined for every material. The unsaturated strength component is optional and
is excluded by default. It is therefore the responsibility of each material to return the necessary strength
components to calculate the 𝐹. Table summarizes the available material models for 2D and 3D analysis.
Table 4. Material models available in 2D and 3D.

Material Model 2D 3D
Mohr-Coulomb Yes Yes
Spatial Mohr-Coulomb Yes Yes
Undrained Strength Yes Yes
Impenetrable Yes Yes
Shear-Normal Function Yes Yes
Anisotropic Strength Yes No
Anisotropic Function Yes Yes
SHANSEP Yes Yes
Hoek-Brown Yes Yes
Compound Strength Yes Yes
High Strength Yes Yes
Bilinear Yes Yes
Strength as Function of Depth Yes Yes
Combined Models Yes Yes
Unsaturated Strength Yes Yes

4.2 Mohr-Coulomb
The Mohr-Coulomb strength equation is a linear relationship between the strength of a material and the
applied effective normal stress:
Equation 39
𝑠 = 𝑐' + (𝜎𝑛 ‒ 𝑢𝑤)𝑡𝑎𝑛𝜙'

where 𝑐 is the effective cohesion, 𝜑 is the effective angle of internal friction, 𝜎𝑛 is the total normal
' '

stress, and 𝑢𝑤 is the pore water pressure. This equation represents a straight line on a shear strength
versus normal stress plot (Figure 77). The intercept on the shear strength axis is the effective cohesion
and the slope of the line is the angle of internal friction.

59
Figure 77. Mohr-Coulomb shear strength envelope.

4.3 Spatial Mohr-Coulomb

The Spatial Mohr-Coulomb model provides the ability to vary the unit weight, effective friction angle,
and/or effective cohesion spatially in the domain in the x-y plane. There are two options: linear or
spatial. In the case of a linear variation, the material parameter is defined as a function of either x-
coordinate or y-coordinate (Figure 78). A spatial function requires three data points including x-
coordinate, y-coordinate, and the material parameter. Figure 79 presents a case in which the undrained
strength (section 4.4) is defined using a linear function in the top layer and spatial function in the
bottom unit.

Figure 78. Unit weight versus y-coordinate.

60
Figure 79. Undrained strength spatial variation in subsurface.

4.4 Undrained Strength

'
The Undrained Strength option provides a convenient approach for setting 𝜙 to zero in the Mohr-
Coulomb model. The shear strength of the material is only described by the cohesion (i.e., undrained
shear strength) and the pore water pressure has no effect on the shear strength calculated at the base
of the column.

4.5 Impenetrable (Bedrock)

The Impenetrable strength material model can be used to force the slip surface to traverse the contact
with a high strength layer (Figure 80). The columns inside the sliding mass are terminated at the
boundary with the impenetrable material. In contrast to Weak Layers (Section 3.7.2), the strength at the
base of the columns is defined by the material overlying the impenetrable material, that is, within the
sliding mass.

Figure 80. Composite slip surface generated by the Impenetrable material model.

4.6 Shear-Normal Function

The shear-normal material model can be used to define a general strength envelope (Figure 81). At solve
time, the effective normal stress at the base of each column is calculated. The effective stress is equal to
the total stress when pore water pressures are not defined. The effective normal stress is then used to

61
 determine the tangent slope of the spline function and the corresponding y-intercept, which represents
the friction angle and cohesion, respectively.


25
Shear Stress

20

15

10

5
C

0
0 20 40 60 80 100
Normal Stress

Figure 81. Shear strength versus normal stress function.

4.7 Anisotropic Function

The Anisotropic Function model is a general strength model for handling anisotropy that can be used in
' '
2D and 3D. The strength parameters 𝑐 and 𝜙 are multiplied by a Modifier Factor that is defined as a
function of inclination 𝛼, the inclination of the column base measured from a local axis 𝑋' aligned with
the sliding direction (Figure 82). Anisotropy dipping in the direction of sliding (Case 1) has a negative
inclination, while that dipping in a direction opposite to sliding (Case 2) has a positive inclination (Figure
82). The sign convention of the Modifier Function therefore holds regardless of the analysis being left-
to-right or right-to-left in a 2D analysis (Figure 83).

Figure 82. Definition of inclination 𝛼 relative to the local axis 𝑋' aligned with the sliding direction: left to right and right to left
sliding directions.

62
Figure 83. Modifier factor versus inclination for Case 1 and Case 2 in Figure 82.

As an example of Case 1, assume the friction angle decreases from 30° to 20° along anisotropy inclined
at 𝛼 = ‒ 45° ± 2.5°. The Modifier Function could be defined using a piece-wise linear function defined
by 4 points: (-50, 1.0), (-47.5, 0.67), (-42.5, 0.67), and (-40, 1.0), which assumes the transition from 30°
to 20° occurs over a 2.5° inclination range (Figure 84).

Figure 84. Example of a spline data point Modifier Factor function.

4.8 Anisotropic Strength

The Anisotropic Strength model calculates the cohesion and friction angle at the base of a column from
the horizontal and vertical input values as:
Equation 40
𝑐' = 𝑐ℎ' cos2 𝛼 + 𝑐𝑣' sin2 𝛼
Equation 41
𝜙' = 𝜙'ℎcos2 𝛼 + 𝜙'𝑣sin2 𝛼

where 𝛼 is the inclination of the column base (Figure 82). This material model can only be used for 2D
analysis.

63
4.9 SHANSEP
The SHANSEP model (Stress History and Normalized Soil Engineering Properties; Ladd and Foott, 1974;
Ladd, 1991) varies the undrained shear strength with vertical effective stress. There are three options
for parameterizing the model:
𝑠 ' '
 Specify a ratio 𝑆 = 𝑢 𝜎𝑣, where 𝑠𝑢 is undrained strength and 𝜎𝑣 is vertical effective stress;
 Specify a functional relationship between undrained shear strength and vertical effective stress
using a generalized spline function; and,
 Define the parameters of the closed-form equation:

Equation 42
𝑠𝑢 = 𝜎𝑣' 𝑆𝑁𝐶𝑂𝐶𝑅𝑚

where 𝑂𝐶𝑅 is over-consolidation ratio, 𝑚 an exponent controlling non-linearity, and 𝑆𝑁𝐶 the ratio of
undrained strength to vertical effective stress for the normally compressed state. The over-consolidation
𝑠 '
ratio can vary with y-coordinate, which causes the ratio 𝑢 𝜎𝑣 to vary nonlinearly with depth if 𝑚 < 1.

Consider the case of placing an embankment on soft ground as illustrated in Figure 85. The first 2m
below the ground is a desiccated layer and the water table is at the contact between the desiccated
layer and the underlying soft soil. The weight of the sand fill will increase the total overburden stress as
well as the pore water pressure in the soft clay layer. Since the effective overburden stress is required in
the SHANSEP model, special care is required to make sure that the pore water pressure in the soft clay
layer is modeled correctly.

Sand fill

Desiccated layer

Soft clay

Figure 85. Case with SHANSEP model applied to soft clay.

To get the correct pore water pressure and consequently the effective stress in the soft clay, it is
necessary to do the following:

 The pore water pressure is defined using a piezometric surface with B-Bar;
 The piezometric surface representing the water table is assigned only to the soft clay;
 A B-bar value of 1.0 should be assigned to the soft clay;
 The sand fill is flagged as the material that the “weight will be added” in the B-bar pore water
pressure calculation of the soft clay layer; and,

64
  The desiccated layer does not have a B-bar and does not have an assigned piezometric surface.

Using a B-bar of 1.0 ensures that the pore water pressure increase in the soft clay layer is commensurate
with the unit weight of the fill. In other words, the additional vertical stress from the sand fill is added to
the initial pore water pressure represented by the specified piezometric surface. Accordingly, the
vertical effective stress in the soft clay is not affected by the sand fill. This is the objective in the
SHANSEP model – the undrained shear strength is the same before and after the fill placement.

The pore water pressures used in the stability calculations should be thoroughly checked and verified by
plotting the pore pressure along the slip surface to make sure that appropriate shear strength is being
used. Surcharge loads as well as surface load due to ponded water on the ground surface are used in the
computation of the vertical stress and therefore are considered in the SHANSEP formulation. However,
individual point loads are not considered.

4.10 Hoek-Brown
The Hoek and Brown model is a nonlinear shear strength model for rock. Hoek, Carranza-Torres, and
Corkum (2002) present the following relationship for the principal stresses at failure:
𝜎3
𝜎1 = 𝜎3 + 𝜎𝑐𝑖 𝑚𝑏
( 𝜎𝑐𝑖
+𝑠
) 𝑎 Equation 43

where 𝜎1 is the major principal stress, 𝜎3 is the minor principal stress, 𝜎𝑐𝑖 is the uniaxial compressive
strength of the intact rock, and 𝑚𝑏, 𝑠, and 𝑎 are material constants. GeoStudio provides the option to
calculate the material constants based on the following expressions:

𝑚𝑏
( 𝐺𝑆𝐼 ‒ 100

= 𝑚 𝑒 28 ‒ 14𝐷
) Equation 44
𝑖

1 1 Equation 45
𝑎 = + (𝑒 ‒ 𝐺𝑆𝐼 15 ‒ 𝑒 ‒ 20 3)
2 6

( 𝐺𝑆𝐼 ‒ 100

𝑠 = 𝑒 9 ‒ 3𝐷
) Equation 46

where 𝐺𝑆𝐼 is the Geological Strength Index (0 – 100), 𝑚𝑖 is a property of the intact rock, and 𝐷 is the
rock mass disturbance factor (0-1). For intact rock, GSI is equal to 100 and s is calculated as 1.0. In an
earlier form of the failure criterion, the coefficient 𝑎 was assumed to be a constant value of 0.5. It is now
considered to be a variable dependent on GSI.

To establish the strength curve, 𝜎1 is computed for a range of 𝜎3 values given the material parameters.
The default range for 𝜎3 is from the tensile strength, which is a negative value, up to half of the uniaxial
compressive strength, where the tensile strength is computed as:
𝜎𝑐𝑖
Equation 47
𝜎𝑡 =‒ 𝑠
𝑚𝑏

65
Once the 𝜎1- 𝜎3 pairs at failure are known, a series of shear strength (𝜏) and normal stress (𝜎𝑛) data
points are computed as follows:

𝜎1 𝜎3
𝜎𝑟𝑎𝑡𝑖𝑜 =
𝜎3
= 1 + 𝑎𝑚𝑏 𝑚𝑏
( 𝜎𝑐𝑖
+𝑠
) 𝑎‒1 Equation 48

𝜎1 + 𝜎3 𝜎1 ‒ 𝜎3 𝜎𝑟𝑎𝑡𝑖𝑜 ‒ 1
𝜎𝑛 = ( 2 )(
‒
2 ) (
×
𝜎𝑟𝑎𝑡𝑖𝑜 + 1 ) Equation 49

𝜎𝑟𝑎𝑡𝑖𝑜
Equation 50
𝜏 = (𝜎 1 ‒ 𝜎 3 )
𝜎𝑟𝑎𝑡𝑖𝑜 + 1

Figure 86 presents the strength function generated for a 𝜎𝑐𝑖, 𝑚𝑖, GSI, and 𝐷, of 20,000 kPa, 10, 45, and
1.0, respectively. At solve time, the strength at the base of a column is calculated in the same manner as
a Shear-Normal material model.

1.4

1.2
Shear Stress (x 1000)

1.0

0.8

0.6

0.4

0.2

0.0
-0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5
Normal Stress (x 1000)

Figure 86. Strength function based on generalized Hoek-Brown parameters.

The following three tables provide some guidelines for selecting and specifying uniaxial compressive
strengths, 𝑚𝑖 constants and GSI values. The 𝐷 parameter is an indication of the amount of disturbance
the rock may have undergone during excavation or blasting. A 𝐷 value of zero means the rock has
undergone very little disturbance, while a 𝐷 value of 1.0 represents a condition where the rock has
undergone significant disturbance due to heavy production blasting that may exist in situations such as
an open pit mine (Hoek et al, 2002).

66
Table 5. Field estimates of uniaxial compressive strength. (after Hoek, 2000).

67
Table 6. Values of constant 𝑚𝑖 for intact rock by rock type; values in parentheses are estimates (after Hoek, 2000).

68
Table 7. Estimates of Geological Strength Index based on geological description (after Hoek, 2000).

4.11 Compound Strength

The Compound Strength model is used to simulate anisotropic (i.e., directionally dependent) strength in
soil or rock, where the strength of the discontinuities varies considerably from the surrounding (intact)
material. The discontinuity could represent bedding, schistosity, joint, foliation, cleavage, fractures, or
faults. The Compound Strength model is like the Anisotropic Surface (Section 3.8) except that the
discontinuity is defined in the model by means of dip and dip direction, which is ideal for ubiquitous
planar features such as those that exist in jointed rock masses.

The model requires a base material representing the intact soil/rock matrix, a material model for each
individual joint or bedding surface, and two angles A and B. The A and B parameters associated with

69
each anisotropic feature are angular tolerances around the discontinuity that control the transition of
shear strength from the anisotropic material to the intact material (Figure 59). The solver calculates the
shear strength based on the angle between the column base and the joint plane (Equation 33). Within
the range of A, the joint/bedding shear strength is used. Outside the range of B, the shear strength of
the base material is used. For angles between A and B, the shear strength is computed based on a linear
interpolation between the shear strengths of the base and the joint materials (Equation 33 through
Equation 35).

4.11.1 Dip and Dip Direction

Dip and dip direction is a measurement convention used to describe the orientation of a planar
geological feature. The dip is the steepest angle of descent of a tilted feature relative to a horizontal
plane. The dip direction is the azimuth (i.e., compass direction) of the dip line of the planar feature; that
is, the azimuth of the imagined line inclined downslope. Dip is always entered as a positive value
between 0 and 90 degrees. Dip direction, being an azimuth, is clockwise positive from North. In
GeoStudio, North is aligned with the negative z-axis. Figure 87 shows two joints dipping at 45 degrees.
Joint 1 and 2 have dip directions of 225 and 135, respectively. In a 2D analysis, the domain is in the x-y
plane, North is in the out-of-plane direction (i.e., into the page), and the joints are implicitly rotated
about the y-axis and viewed along the strike (Figure 87). As discussed subsequently, the dip direction
still plays an important role in a 2D analysis.

Figure 87. Dip and dip direction of two joints looking North along the negative z-axis.

Figure 88 shows pie charts for 3 joint sets with each chart depicting the dip direction, A and B angles,
and the materials for the base and anisotropic features. A plane with a dip direction between 0 and 180
degrees like Joint Set 2 in Figure 87, has the dip in the negative quadrant (despite the dip value being a
positive value). Conversely, a plane with a dip direction between 180 and 360 degrees like Joint Set 1 in
Figure 87, has the dip in the positive quadrant. The strength anisotropy is symmetrical around the center

70
of the pie chart. Hence, the angular range could appear to be discontinuous on the pie chart if the
anisotropy influence range crosses the vertical line as shown in the Figure 87. The pie chart should be
used to visually inspect that each joint set is correctly represented on the domain. Consider, for
example, analyzing the stability of the left and right sides of a geometrically symmetrical domain (Figure
89). The left-to-right analysis would produce a very different solution than the right-to-left solution if the
same compound strength material was used.

71
Figure 88. Strength anisotropy represented by a pie-chart.

Figure 89. Joint sets in a 2D analysis: (a) left-to-right (b) right-to-left.

Typically, a 2D analysis will be done with the assumption that the section is aligned with the dip
direction of the structures or bedding. As shown in Figure 90, the true dip (𝛿) should be used in this

72
case. If the section is at an oblique angle relative the strike of the structures (𝛽), then the apparent dip (
𝛼) could be used. The apparent dip can be calculated from the angle 𝛽 and 𝛿 as:

Equation 51
α = tan ‒ 1 (sin 𝛽tan 𝛿)

Conversely, the true dip can be calculated from the apparent dip as:

δ = tan ‒ 1 (tan 𝛼 sin 𝛽) Equation 52

Figure 90. Example of true dip and apparent dip.

4.12 High Strength

A High Strength material model provides a convenient approach for including materials with very high
strength such as a concrete retaining wall. This material model effectively acts as a slip surface filter: any
slip surface that passes through the material is discarded. This is based on the idea that the shear zone
cannot pass through the high strength material. The material unit weight is included in the calculation of
the column weight, thereby effecting the normal stress at the base of the column.

4.13 Bilinear
'
Figure 91 presents a bilinear Mohr-Coulomb material model, which is defined by a single input for 𝑐 ,
'
two inputs for 𝜙 , and a single input for the normal stress. For column base normal stresses greater than
'
the specified normal stress, the 𝜙 2line segment is used to compute an equivalent cohesion intercept.

73
Figure 91. Bi-Linear Mohr-Coulomb shear strength envelope.

4.14 Strength as a Function of Depth

Soft normally or slightly over-consolidated soil usually exhibits an increase in undrained shear strength
with depth. There are two options for describing this type of strength profile. One option is referenced
to the top of the material layer and the other is referenced to a specified datum. Figure 92 presents a
case in which the undrained strength is referenced to the top of a layer, with the rate of increase capped
by a maximum undrained strength input. The top of the soil layer does not need to be planar because
the depth is calculated for each column. Similarly, Figure 93 presents the option using a datum to
calculate depth.

Figure 92. Undrained shear strength as a function of depth.

74
Figure 93. Undrained shear strength as a function of depth below a datum.

4.15 Combined Models: Frictional-Undrained

The frictional-undrained shear strength model can be used to consider both the drained and undrained
Mohr-Coulomb strengths. The drained strength governs until the base normal stress produces a
strength that exceeds the undrained strength (𝑆𝑢; Figure 94). The effective cohesion (𝑐 ) and undrained
'

shear strength (𝑆𝑢) can vary with depth. A reference value for 𝑐' and 𝑆𝑢 can be specified at the top of the
stratigraphic unit or at an elevation datum together with a rate of change with depth. Alternatively, a
𝑐'
ratio of 𝑆𝑢 can be specified. If the value is non-zero, the effective cohesion (𝑐') is computed from the
specified ratio and the value of 𝑆𝑢.

Figure 94. Combined frictional-undrained strength model.

75
4.16 Unsaturated Shear Strength
In the same way that positive pore water pressure decreases the effective stress and thereby decreases
the shear strength, negative pore water pressure increases the effective stress and in turn increases the
strength. The shear strength of unsaturated soil is:
Equation 53
𝑠 = 𝑐' + (𝜎𝑛 ‒ 𝑢𝑎)𝑡𝑎𝑛𝜙' + 𝑆𝑒(𝑢𝑎 ‒ 𝑢𝑤)𝑡𝑎𝑛𝜙'

where 𝑐 is the effective cohesion, 𝜙 is the effective angle of internal friction, 𝜎𝑛 is the total normal
' '

stress, 𝑢𝑤 is the pore water pressure, 𝑢𝑎 is the pore-air pressure, and 𝑆𝑒 is the effective degree of
saturation given as (van Genuchten, 1980):
𝜃𝑤 ‒ 𝜃𝑟
Equation 54
𝑆𝑒 =
𝜃𝑠 ‒ 𝜃𝑟

where 𝜃𝑤, 𝜃𝑠, and 𝜃𝑟 are the volumetric water content, saturated volumetric water content, and
residual volumetric water content, respectively. The last term of Equation 3 represents the additional
strength derived from matric suction, (𝑢𝑎 ‒ 𝑢𝑤), in the unsaturated zone. The effective degree of
saturation in Equation 53 ensures that the role of matric suction in generating shear resistance is
weighted based on the water content of the soil.

4.17 Soil Unit Weight

As discussed in Section 2, the force and moment equilibrium equations require the weight of each
column. The weight of each column is calculated by numerical integration over the length of the column
to account for spatial variability in the unit weight. By default, the unit weight of the saturated and
unsaturated zones is assumed equivalent. However, the unit weight of the unsaturated soil can be
specified uniquely by either entering a constant value, or using a volumetric water content function to
calculate the unit weight according to:
𝐺 ‒ 𝑆𝑒 Equation 55
𝛾 = 𝛾𝑤
1+𝑒

where 𝐺 is the specific gravity, 𝑆 is the degree of saturation, and 𝑒 is the void ratio.

It should be noted that the application of a unique unit weight in the unsaturated zone has a marginal
effect on the results in most cases, particularly if the unit weight is correctly linked to the degree of
saturation. For illustration purposes, assume that 𝐺 = 2.7, 𝛾𝑤= 10 kN/m3, and e = 0.7. Below the water
table where the soil is saturated (S = 1), the unit weight is 20.0 kN/m3. Above the water table where the
soil is unsaturated, assume the degree of saturation is 80%, (i.e., S = 0.8) the unit weight would be 19.2
kN/m3, or about a 4% difference.

From a stability analysis point of view, the 4% difference in the unit weight is insignificant. First, the
degree of saturation in the capillary zone is 100%, so the total unit weight is unchanged. Secondly, factor
of safety calculations are not very sensitive to the unit weight. If the column weight increases, the
gravitational driving force increases, but the shear resistance also increases because of a higher normal

76
stress at the column base. Therefore, defining different unit weights above and below the water table
usually has little effect on the resulting factor of safety.

Figure 95 presents an example of a coarse rock fill in which the unit weight may very significantly in the
unsaturated zone. The air entry value of coarse rockfill is very low, so the capillary fringe would be small,
causing the transition toward the residual water content to occur almost immediately above the water
table. Assume that 𝑒 is 1.0, 𝐺 is 2.7 and 𝛾𝑤 is 10 kN/m3. Below the water line, where S is 1.0, the unit
weight is 18.5 kN/m3. Above the water line where S is zero, the unit weight is 13.5 kN/m3. This is a
significant difference and consequently worth the effort to specify two different unit weights.

Figure 95. Example of a rockfill with significant difference in saturated and unsaturated unit weights.

5 External Loads
External loads, 𝐷, enter into the equations for the column base normal force, intercolumn normal force,
and the force and moment factor of safety equations (Section 2). These external loads are often stresses
that must be converted to forces using the geometry of the columns (i.e. area of column faces).

5.1 Ponded Fluid

Figure 96 presents a 2D example of ponded fluid at the toe of a slope. The fluid provides a surcharge
load that acts normal to the ground surface, often providing a stabilizing load. An automatic fluid
surcharge load is included on the ground surface in 2D and 3D analyses if the pore water pressure is
defined using a piezometric surface or finite element analysis, provided that the pore water pressure is
positive. In the case of a piezometric surface, the pore water pressure is calculated based on the height
above the ground surface and the user-defined water unit weight as 𝛾𝑤ℎ𝑤. Note that piezometric
surfaces are associated with materials. As such, the automatic ponding surcharge will only apply if the
material exposed at the ground surface has been associated with a piezometric surface.

77
Figure 96. Exampled of ponded water in a 2D analysis.

A force representing the weight of the ponded fluid is added to the top of each column as a point load
(Figure 97). The force is equal to the water pressure at mid-column multiplied by the column area, which
is equal to the width in 2D. The computed surcharge forces on a column can be verified by viewing
Column Information. In the case of a vertical wall, the horizontal hydrostatic force is applied
automatically to the vertical side of the column.

Figure 97. Example of a force due to ponded water at the top of a column.

There are some cases in which the ponded fluid has a different unit weight than the groundwater. The
example in Figure 98 demonstrates a brine pond in which the groundwater piezometric surface applies
to all three materials. The additional unit weight of the brine is modeled using a surcharge load.
Accordingly, any trial slip surface passing through the light green material will use the piezometric
surface and water unit weight to calculate the pore water pressure at the base of the columns, while the
surcharge load at the ground surface will be calculated using the piezometric surface plus the surcharge
load.

78
Figure 98. Example of an additional surcharge load applied to account for brine.

5.2 Surcharge Loads

Figure 99 shows two surcharge loads applied to the ground surface representing an equipment or
building load and a soil berm. In 2D, a surcharge load is defined/drawn using a polyline and a user-
specified unit weight. The loading can be applied in a vertical direction or normal to the ground surface
line. The normal direction can be used when modeling water or other fluid types under hydrostatic
conditions. The surcharge load regions are cross-hatched and arrows are used to indicate the loading
direction. Surcharge load can be specified as positive or negative

Figure 99. Ground surface surcharge loads.

A force representing the surcharge is added to the top of each column as a point load. The force is
equal to the vertical distance at mid-column from the ground surface to the top of the surcharge load
times the column width, times the specified surcharge (or unit weight).

A useful method for applying a uniform pressure is to make the distance from the ground surface to the
surcharge load line unity (1.0). The specified surcharge will then equal the applied pressure. This can be
useful when simulating a footing pressure. The surcharge force is displayed at the top of the column at
the applied direction (Figure 100).

79
Figure 100. Example of a surcharge force at top of column.

In the example shown in Figure 101, the unit weight of water is specified as 10 kN/m3, but the fluid
material has a unit weight of 11 kN/m3. To model the additional weight of the fluid material, a
surcharge load of 1 kN/m3 is added to the pond with the load direction selected as normal.

Figure 101. Example of ponded water and surcharge load.

The detail column force of the column #15 is shown in Figure 102. The weight of the ponded fluid is
shown as a point load on the top of the column. The applied 17.916 kN vertical load is computed by the
unit weight of the fluid (11 kN/m3/m) multiplied by the depth of the pond (2m) and the width of the
column (0.81437m).

80
Figure 102. Example of a fluid surcharge force at top of column.

The detailed column force of the vertical column (column #14) is shown in Figure 103. The hydrostatic
force of the ponded water is shown with a horizontal force shown on the top right corner of the free
body diagram. The total is 22 kN which is equal to the total unit weight of the fluid (11 kN/m3/m)
multiplied by the square of the depth of the pond (4m2) divided by 2.

Figure 103. Example of a hydrostatic force on the side of column.

It is important to note that the “normal” direction of the surcharge load must be chosen in order to
model the hydrostatic nature of ponded waters. When the normal direction is chosen, the surcharge
load is applied to the top of each column in the direction normal to the ground surface. In the case of
ponding next to a vertical wall, the hydrostatic force is applied horizontally (normal) to the side of the
vertical column.

81
5.3 Point Loads
Point loads are applied to the ground surface or anywhere with the model domain. A point load is
included in the equilibrium equations if it sits within a particular column.

6 Pore Water Pressure

6.1 Introduction
GeoStudio has several methods for specifying the pore water pressure conditions in a slope stability
analysis including piezometric surfaces with or without B-bar or 𝑅𝑢 coefficients, spatial functions, and
finite element computed pore water pressures. It should be noted that the most realistic position of the
critical slip surface is obtained when effective strength parameters are used in the analysis. Effective
strength parameters, however, are only meaningful when they are used in conjunction with pore water
pressures given the strength relationship. It is therefore important to define the pore water pressures in
a manner that best reflects the actual conditions.

6.2 Piezometric Surfaces

A piezometric surface represents an imaginary surface that defines the level to which water would rise
in a confined aquifer based on piezometer measurements. In GeoStudio, a piezometric surface can also
be used to approximate a water table in an unconfined stratigraphic unit. If the piezometric surface lies
above the ground surface, then it represents an artesian condition for a confined stratigraphic unit; or, a
ponded condition for an unconfined aquifer or confining unit.

The pore water pressure at the base of each column is calculated by multiplying the vertical distance
from the column base mid-point to the piezometric surface by the unit weight of water (Figure 104). If
the column base mid-point is located above the piezometric surface, then the pore water pressure is
negative, but additional strength due to the matric suction is assumed to be zero unless specified (see
Section 4.16).

H (water)
1

u = H(water) x water unit weight

Figure 104. Pore water pressure from a piezometric surface.

In a 2D analysis, a piezometric surface is defined by either drawing a polyline or entering x-y coordinates
directly (Figure 105). In 3D, a background mesh is first tagged as a piezometric surface (Figure 106), then

82
added to the list of piezometric surfaces (Figure 107). The piezometric surface is then assigned to the
various materials in the domain in the same manner as 2D. Multiple piezometric surfaces can be added
and assigned to materials within a single analysis.

Figure 105. Definition of a piezometric surface in 2D.

Figure 106. Designation of a background mesh as a piezometric surface.

83
Figure 107. Adding a background mesh to the list of piezometric surfaces and applying the surface to materials.

There are cases in which the piezometric surface may not apply to a material in the domain. Consider
the construction of a sand embankment on soft clay soil as shown in Figure 108. Placement of the sand
embankment produces excess pore water pressure in the clay represented by the piezometric surface,
while the sand embankment remains full drained. As such, the piezometric surface is only assigned to
the clay unit.

Sand

Clay

Figure 108. Piezometric surface assigned only to the clay unit.

6.3 Ru Coefficients
The pore water pressure ratio 𝑅𝑢 is a coefficient that relates the pore water pressure to the vertical
overburden stress. The concept of 𝑅𝑢 was developed primarily for use with stability charts (Bishop and
Morgenstern, 1960). The coefficient is defined as:
𝑢𝑤
Equation 56
𝑅𝑢 =
𝛾𝐻𝑠

84
where 𝐻𝑠 is the height of the soil column. The denominator of the 𝑅𝑢 equation includes external loads
(Budhu, 2007). Accordingly, ponded water or surcharge loads are included in the calculation of pore
water pressure at the base of a column. By rearranging the variables, the pore water pressure becomes:
𝛾𝐻𝑠 Equation 57
𝑢𝑤𝑅 =
𝑅

Theoretically, the 𝑅𝑢 approach is only valid for hydrostatic flow systems in which the phreatic surface is
horizontal. In a multi-layered stratigraphic section, 𝑅𝑢 can be applied uniquely to each soil unit.
Moreover, Ru is not intended to model the generation of excess pore water pressures due to loading.

6.4 B-bar Coefficients

The Skempton B-bar coefficient (𝐵̅ ) is used to approximate excess pore water pressure due to loading
as:
∆𝑢𝑤
Equation 58
𝐵̅ =
∆𝜎𝑣

where ∆𝑢𝑤 is the change in pore water pressure, that is, the excess pore water pressure, and ∆𝜎𝑣 is the
change in vertical total stress. Consider the case of a sand embankment constructed over a clay
foundation as shown Figure 109. Prior to construction, the pore water pressure at the base of column 10
is 2.14m x 9.81 kN/m3 = 20.99 kPa. Assuming a 𝐵̅ = 0.2, the excess pore water from the sand load is
7.22m x 20.0 kN/m3 x 0.2 = 28.88 kPa. The total pore water pressure at the base of column 10 is
therefore 49.87 kPa.

Sand

10

H(s) = 7.22

H(w) = 2.14

Clay

Figure 109. Excess pore water pressure due to construction of an embankment.

The B-bar condition is defined along with a piezometric surface. There are three aspects to the definition
of a B-bar scenario: 1) assigning the piezometric surface to the materials in the domain (Section 6.2); 2)
ascribing a B-bar coefficient to the materials; and, 3) indicating if a material generates excess pore water
pressure via an ‘Add Weight’ toggle. In the example above, the piezometric surface would be assigned to
the In situ clay, but not the overlying sand, which is placed during construction. A B-bar coefficient

85
would be assigned to the clay, but not the sand. Finally, the sand embankment would be flagged with
‘Add Weight’ (Figure 110).

Figure 110. Configuration for a B-bar analysis.

6.5 Spatial Function

The pore water pressure conditions can be specified using the pore water pressure head defined at
discrete points in a 2D analysis (Figure 111). A smooth surface that passes through all the specified
points is constructed using kriging. Once the surface has been constructed, the pore water pressure can
be determined at any spatial coordinate in the domain. Spatial functions can be defined only the vertical
x-y plane and are therefore assumed constant in the z-direction for 3D analyses.

Figure 111. Grid of pore water pressure head.

The pressure head spatial function is useful when the pore water pressures are irregular or non
hydrostatic. A spatial function is also useful for directly including field piezometric measurements in an
analysis.

When using the pressure head with spatial function option, it is necessary to give careful thought to all
the locations the pore water pressure is known. For example, some arbitrary data points with pressure
head equal to zero may be required in the unsaturated zone if matric suction is to be ignored. If there is
ponded water, a pressure head equal to the depth of the water would need to be specified along the a
portion of the ground surface. On a seepage face where the pore water pressure is zero, the pressure
head would need to be specified as zero. In summary, the pore water pressure needs to be specified at
all locations where the pore water pressure condition is known with certainty.

86
6.6 Finite Element Pore Water Pressures
In GeoStudio, the pore water pressures in a stability analysis can be defined using a steady-state or
transient groundwater flow analysis, a transient stress-strain consolidation analysis, or a dynamic stress-
strain analysis with excess pore water pressure calculations. In keeping with piezometric surfaces, an
automatic surface load is applied along the ground surface if the pore water pressures are positive (see
Section 5.1).

To use the finite element computed pore water pressures, the coordinates of the column base centroid
are used to locate the element in the FE mesh. Next, the local coordinates within the element are
determined from the cartesian coordinates. The pore water pressures at the local coordinate is then
calculated via interpolation from nodal values. The advantage of this approach is that the pore water
pressures can have any irregular distribution and represent conditions at various times.

Figure 112 and Figure 113 present an example of using irregular finite element computed pore water
pressures in a 2D stability analysis. A low permeability layer in the slope causes a perched water table to
form due to the infiltration. This type of pore water pressure distribution can be modeled in SEEP/W or
SEEP3D and then used directly in the stability analysis (2D or 3D).

Infiltration
Flow path
Perched Water Table

Low permeability layer

Figure 112. Pore water pressures in a 2D seepage analysis.

Figure 113. Pore water pressures in a 2D stability analysis.

Figure 114 shows the distribution of pore water pressure along the slip surface. The pore water pressure
starts out negative at the crest, then becomes positive in the perched water zone, becomes negative
again in the unsaturated zone and then is positive through the base saturated zone.

87
 40

30
Pore-Water Pressure

20

10

-10

-20
0 10 20 30 40
Distance along slip surface

Figure 114. Pore water pressure distribution along the slip surface.

In the case of a transient finite element with multiple time steps, a specific time step may be selected
for the stability analysis. Alternatively, all time steps can be included, producing a FOS through time
(Figure 115). This option is particularly useful when the pore water pressures are changing significantly
due to change hydraulic conditions at the ground surface, such as a flood event for a levee, filling and/or
drawdown of earth dam, groundwater pumping around an open pit, and a multitude of other transient
scenarios.

Figure 115. Factor of safety versus time.

88
7 Reinforcement and Structural Components
There are several reinforcement types in GeoStudio for 2D analysis including anchors, nails,
geosynthetics, and piles. In a limit equilibrium analysis, the effect of the reinforcement is represented
using a point load. The point loads act on the free body, which is the potential sliding mass, and must
therefore be included in the moment and force equilibrium equations (see Section 2.4).

7.1 Reinforcement Load Calculation

Soil-structure interaction consists of an exchange of stress, whereby deformation within the ground
transfers load into the structure. Consider the grouted anchors shown in Figure 116. The reinforcement
system comprises face anchorage and a tendon. A portion of the tendon is bonded to the soil/rock,
leaving a length that is free to elongate. The active earth pressure acting on the wall leads to
deformation, causing elongation of the of the tendon. The length of tendon that is bonded transmits the
applied tensile load into the ground. The same concept applies for pre-stressing: the load is transferred
through the free length and into the length bonded to the soil/rock.

Figure 116. Example of grouted anchors.

Figure 117 presents a simplified conceptual model of the load transfer into the ground. The load realized
at the anchorage is transferred down the free length due to elongation, that is, strain, in the tendon. The
bonded length is therefore subject to a pull-out load, which mobilizes shear stress at the interface. This
shear stress could be converted to a force by multiplying by the area of the interface. More importantly,
the shear stress must balance the tensile load in the tendon. In a limit equilibrium analysis, a
reinforcement load, 𝑇, can be applied to the sliding mass to represent this interaction. For convenience,
the reinforcement load is applied at the intersection with the slip surface along the action.

89
Figure 117. Conceptual model of load transfer into ground.

In GeoStudio, pullout resistance (𝑃𝑅) is generally the parameter used to capture the stress transferred
into the ground. The reinforcement load is calculated as:
Equation 59
𝑇 = 𝑃𝑅(𝐴𝑟𝑒𝑎)
where the Area depends on the length of reinforcement behind the slip surface. The following sections
detail the customization for each reinforcement type.

7.1.1 Face Anchorage

Face anchorage is common for most reinforcement. There may be scenarios, however, where the
reinforcement is not anchored at the face of the slope, particularly for geosynthetic reinforcement.
Shear stress would therefore be mobilized at the reinforcement interface within the sliding mass. If the
face anchorage option is set to ‘no’, the reinforcement load is calculated for the portion of
reinforcement behind the sliding mass and the portion within the sliding mass. The lesser of these two
values is then compared to the tensile capacity (discussed below).

7.2 Geosynthetics
For geosynthetics, such as geotextiles and geogrids, two methods of obtaining the pullout resistance are
available. The selection is dependent on the stress transfer mechanism of the reinforcement. The
pullout resistance can be specified directly if passive resistance is the dominant stress transfer
mechanism. Passive resistance refers to the development of bearing type stresses on relatively stiff
members of the reinforcement that are situated normal to the direction of pullout (e.g. geogrid).
Alternatively, the pullout resistance can be calculated from the overburden stress if frictional resistance
dominates the stress transfer mechanism (e.g. geotextiles).

The pullout resistance is specified with the units of force per unit length of geosynthetic per unit width
in the out-of-plane direction. The calculated pullout resistance option requires the specification of:

 Interface adhesion (𝑆𝐼𝐴): apparent cohesion (adhesion) if effective drained soil strengths are
being considered. Alternatively, the parameter can be used to specify the undrained strength
at the interface between the geosynthetic and soil.
 Interface shear angle (𝛿): angle of interface shearing resistance if effective drained soil strengths
are being assumed.

90
  Surface area factor (𝑆𝐴𝐹): used to account for mobilized pullout resistance on the top and
bottom of the geosynthetic. The default parameter is 2 for resistance on both sides of the
geosynthetic; however, this input could be used to account for different soils on the top and
bottom of the geosynthetic.

The calculated pullout resistance PR is given by:

Equation 60
𝑃𝑅 = (𝑆𝐼𝐴 + 𝜎𝑣' 𝑡𝑎𝑛𝛿)𝑆𝐴𝐹

The effective overburden stress is calculated based on the height of soil above the intersection point
with the base of the column and includes the effects of surcharge loads. Point loads are not included in
the calculation. Regardless of the selected approach, the following inputs are required:

 Resistance reduction factor (RRF): can be used as a “scale effect correction factor” to account
for nonlinear stress reduction over the embedded length of highly extensible reinforcement.
 Tensile capacity (TC): tensile strength of the reinforcement.
 Reduction factor: accounts for the reduction in the ultimately tensile capacity of the
reinforcement due to physical processes such as installation damage, creep, and durability.
 Load orientation: a value ranging between 0 and 1 to control the orientation of the pullout force
on the free body. A value of zero 0 applies the load parallel to orientation of the geosynthetic; a
value of 1 applies the force parallel to the column base (Figure 118).

Slice 13 - Morgenstern-Price Method

88.926
5.8436

39.767 26.891

93.232
2.7344
20.794

33.924

Figure 118. Pullout force parallel to column base.

The factored pullout resistance 𝐹𝑃𝑅 per unit length of geosynthetic behind the slip surface is calculated
from the pullout resistance 𝑃𝑅 as:
Equation 61
𝐹𝑃𝑅 = 𝑃𝑅/𝑅𝑅𝐹

91
The pullout force (𝑃𝐹) is calculated as:
Equation 62
𝑃𝐹 = 𝐹𝑃𝑅(𝑙)

The maximum pullout force must not exceed the factored tensile capacity calculated as:
Equation 63
𝐹𝑇𝐶 = 𝑇𝐶/(𝑅𝐹)

The pullout force that is applied to the free body is the lesser of 𝑃𝐹 or 𝐹𝑇𝐶.

7.3 Anchors and Nails

Anchors and nails are essentially the same, except that the entire length of nail is bonded. The pullout
resistance (𝑃𝑅) is specified with the units of force per area of bonded length. The following additional
inputs are required for both the nail and anchor reinforcement types:

 Bond diameter (𝐷): the diameter of the bonded section in contact with the soil/rock.
 Resistance reduction factor (𝑅𝑅𝐹): can be used as a “scale effect correction factor” to account
for nonlinear stress reduction over the embedded length.
 Tensile capacity (𝑇𝐶): tensile strength of the reinforcement.
 Reduction factor (𝑅𝐹): accounts for the reduction of the ultimate tensile capacity due to
physical processes such as installation damage, creep, temperature, and durability.
 Shear force: to incorporate a force that represents the shear stress mobilized within the
reinforcement.
 Shear reduction factor: a factor to reduce the input shear force.
 Apply shear: a value ranging between 0 and 1 to control the orientation of the shear force on
the free body. A value of zero 0 applies the force parallel to the orientation of the
reinforcement; a value of 1 applies the force parallel to the column base.
 Factor of Safety Dependent – Yes or No: applies the computed factor of safety to the pullout
resistance and tensile capacity (see Section 7.8).
 Force Distribution: Concentrated or Distributed (see Section 7.6).
 Bond length: required for anchors.
 Bond diameter (𝐷): diameter of bonded length.
 Out-of-Plane Spacing (𝑆): the distance in the out-of-plane dimension between the anchors or
nails.
 Face Anchorage – Yes or No: indicates the presence of anchorage at the face of the structure.

The factored pullout resistance, 𝐹𝑃𝑅, per length of bonded section or nail behind the slip surface is
calculated from the specified pullout resistance 𝑃𝑅 as:
Equation 64
𝐹𝑃𝑅 = 𝑃𝑅(𝜋𝐷)/(𝑅𝑅𝐹𝑥𝑆)

92
The numerator converts the pull-resistance stress to a force, while the denominator normalizes the load
by the spacing for a 2D analysis. Once the length of reinforcement behind the slip surface is known (𝑙),
the pullout force (𝑃𝐹) is calculated as:
Equation 65
𝑃𝐹 = 𝐹𝑃𝑅(𝑙)

The maximum pullout force must not exceed the factored tensile capacity calculated as:
Equation 66
𝐹𝑇𝐶 = 𝑇𝐶/(𝑅𝐹𝑥𝑆)

Note that the tensile capacity must also be normalized by the spacing for a 2D analysis in order to be
commensurate with the 𝐹𝑃𝑅. The pullout force that is applied to the free body is the lesser of 𝑃𝐹 or
𝐹𝑇𝐶.

7.3.1 Grouted Friction Anchor

There are two options for calculating the pullout resistance of an anchor: 1) specify the pullout
resistance; or, 2) calculate the pullout resistance based on the effective overburden stress. Selection of
the latter makes the anchor behave like a grouted friction anchor. A grouted friction anchor accounts for
the stress-dependent frictional strength of the soil-grout interface. The pullout resistance is calculated
using Equation 60.

7.4 Pile Reinforcement

A pile reinforcement provides resistance to shearing in a limit equilibrium analysis. The inputs include
(Figure 119):

 Shear force: represents the shear stress mobilized within the reinforcement.
 Shear reduction factor: a factor to reduce the input shear force.
 Out-of-Plane Spacing (𝑆): the distance in the out-of-plane dimension between the piles.
 Direction of Application: parallel or perpendicular to the slip surface.

Figure 119. Slope reinforced with a pile wall.

93
7.5 User-defined Reinforcement
Reinforcement loads can be incorporated into stability calculations based on a user defined
reinforcement force versus distance function. The reinforcement force can either be a pullout or a shear
force. The distance in the function is the length from the starting point of the reinforcement to the
intersecting point with the slip surface. Other input parameters are:

 Out-of-Plane Spacing (𝑆)

 Reduction factor
 Force orientation: relative to the column base (0 to 1: 0 is Axial, 1 is parallel to column base) or
perpendicular to the reinforcement

The user defined reinforcement allows for consideration of various modes of failure such as pullout,
tensile capacity (of the reinforcement), plate capacity (or connector failure) and so on. The user defined
reinforcement can also be used model variations in the mobilized shear force over the length of a pile.
The application of the user defined reinforcement for simulating a pile requires that the force be applied
perpendicular to the reinforcement. The determination of a shear force function is relatively
straightforward. Only the pullout force along the reinforcement is discussed.

Consider the following reinforcement load/force components:

 𝑇𝑓 – reinforcement-facing connection strength

 𝑇𝑡 – reinforcement tensile strength

 𝑇𝑒 – embedded end strength

 𝑇𝑠 – pullout resistance, 𝑇𝑠 = 𝜏(𝐿𝑎), here 𝜏 is the reinforcement-soil interface shear strength

(F/L2) and 𝐿𝑎 is the effective anchorage length of reinforcement.

Relative to various potential failure modes (Figure 120), the maximum reinforcement forces at the
intersecting point can be determined as follows:

 Stripping (𝐹1) - stripping occurs if the pullout force exceeds the combined capacity from facing
connection and soil-reinforcement adhesive/frictional shear resistance in the sliding mass.
 Pullout (𝐹2) - reinforcement pullout from behind the slip surface when the pullout force exceeds
the combined capacity from embedded end bearing and soil-reinforcement adhesive/frictional
shear resistance in retaining zone.
 Overstressing (𝐹3) – the pullout force cannot exceed the reinforcement tensile capacity.

The reinforcement force versus distance function can be determined by 𝐹 = 𝑚𝑖𝑛⁡(𝐹1, 𝐹2, 𝐹3)

94
Figure 120. Reinforcement forces considering (a) stripping, (b) pullout, and (c) overstressing.

Figure 121 illustrates a reinforcement function for a 23 m long anchor with plate capacity of 40 kN, end
bearing capacity of 100 kN, reinforcement tensile capacity of 200 kN, and soil-reinforcement interface
shear force of 20 kN/m. Interpreting the graph from right-to-left, envisage a slip surface immediately
inside the terminal point of the anchor. The reinforcement load is governed solely by the embedded end
bearing capacity of 40 kN. The bonded length behind the slip surface increases as the slip surface gets
shallower. The pull-out force increases from 40 kN to 200 kN at a rate of 20 kN/m due to the
contribution from pullout resistance. Eventually, the length of reinforcement behind the slip surface
becomes large enough that the tensile capacity is reached. Once the depth of the slip surface is less than
8 m from the face, the reinforcement load is governed by stripping, that is, the contribution of pullout
resistance within the sliding mass and the face connection. It should be noted that if the face connection
capacity exceeds the tensile capacity, then stripping is not possible, and the function would be constant
at 200 kN between 0 m and 18 m.

95
Figure 121. An example reinforcement function.

7.6 Concentrated vs Distributed

Figure 122 presents a free body diagram of a column with a 200 kN reinforcement load and a 20 kN
shear force. The reinforcement load is ‘concentrated’ on a single column. It should be noted, however,
that the ‘Distributed’ option is default for all reinforcement. This option distributes the reinforcement
load equally to all columns below the elevation at which the reinforcement intersects the slip surface.
This helps to promote convergence by smoothing the normal force distribution along the slip surface,
thereby minimizing sharp contrasts in effective stress and shear strength between adjacent columns.

Slice 11 - Morgenstern-Price Method

68.251
3.6654

4.9452
92.09

20
200

58.572

192.52 90.264

Figure 122. Column free body diagram with a reinforcement load.

96
7.7 Manufacturer Library
Several manufacturer geosynthetic reinforcements are incorporated into GeoStudio. The reinforcement
parameters are automatically calculated based on installation and site-specific conditions that influence
the overall reduction factor on the tensile capacity. The pullout resistance (𝑃𝑅) for all manufacturer
reinforcements is calculated according to Equation 60. For manufacturer reinforcement, the interface
shear angle is calculated as:
Equation 67
𝛿 = cot ‒ 1 (𝑎' tan 𝜑')

' '
where 𝑎 is the pull-out coefficient and 𝜑 is the friction angle of the backfill. The pullout coefficient is
either prescribed by the manufacturer based on the backfill material or it is a user-input value.

The inputs for each manufacturer reinforcement are dependent on the product type selected. The
inputs are used to calculate the overall reduction factor (𝑅𝐹𝑂) for calculation of the factored tensile
capacity (FTC) as:
𝑇𝐶 𝑇𝐶 Equation 68
𝐹𝑇𝐶 = =
𝑅𝐹𝑂(𝐹𝑆) 𝑅𝐹𝐼𝐷 × 𝑅𝐹𝐶𝑅 × 𝑅𝐹𝐶𝐻 × 𝑅𝐹𝑊 × 𝑅𝐹𝐷

where 𝐹𝑇𝐶 is the factored tensile capacity, or the long-term tensile strength, 𝑇𝐶 is the tensile capacity,
or the short-term tensile strength, 𝑅𝐹𝑂 is the overall reduction factor, which is calculated based on the
other reduction factors. The factors include:

 𝑅𝐹𝐼𝐷 is the reduction factor for installation damage;

 𝑅𝐹𝐶𝑅 is the reduction factor for creep;

 𝑅𝐹𝐶𝐻 is the reduction factor for chemical/environmental effects;

 𝑅𝐹𝑊 is the reduction factor for weathering; and,

 𝑅𝐹𝐷 is the reduction factor for the extrapolation of data.

Each of the reduction factors are influenced by site-specific and environmental conditions. Required
inputs that determine the individual reduction factors may include soil type or particle size, fill material,
temperature, design life, and/or soil pH levels.

7.8 Factor of Safety Dependency

The effect of reinforcement can conceptually act immediately or develop with some strain. A pre-
stressed anchor, for example, acts immediately. The force is induced by the pre-stressing. The force in a
geosynthetic, on the other hand, may develop over time during construction and during stress re-
distribution of stress upon completion of the construction (Hoek and Bray, 1974). In other words, the
reinforcement forces are mobilized in response to straining in the same way that the shear stress is
mobilized as ground deforms/strains. Applying the factor of safety to the reinforcement load in the
equilibrium equations results in the factor of safety dependent equations presented in Sections 2.4.1
and 2.4.2

97
A point of significance is that the soil strength and the shear resistance arising from the reinforcement
are both divided by the same overall global factor of safety. The implication is that the soil shear
resistance and reinforcement shear resistance are developed and mobilized at the same rate. This is
likely not entirely correct, but that is the assumption with this approach.

The “F of S Dependent” Yes option is considered appropriate for ductile reinforcement such as some
geosynthetics. In this case, it is important to set the Reduction Factors to 1.0. Finally, it is noted that all
manufacturer reinforcements for geosynthetics assume ‘No’ factor of safety dependency.

7.9 Results Visualization

The detailed example files should be consulted for further insights into the functionality of
reinforcement and the graphical features available to assist with interpretation of the results.
Specifically, visual features are available to indicate (a) if the pullout force applied to the free body was
governed by the factored pullout resistance or factored tensile capacity; and (b) the length over which
pullout resistance was mobilized. In addition, the View | Object information can be used to explore both
the reinforcement inputs and values used in the factor of safety calculations.

8 Stress-based Stability Analysis

8.1.1 Calculating the Stability Factor
Section 2.7 discusses stress distributions obtained from a limit equilibrium analysis and shows that these
stress distributions are not necessarily representative of the actual field stresses. This discrepancy
occurs because the formulation only considers equilibrium, and it requires that the factor of safety is the
same for each column. The missing piece of physics is a stress-strain relationship. Including such a
relationship means displacement compatibility is satisfied, which in turn leads to much more realistic
stress distributions.

The stress-based stability analysis technique first establishes the stress distribution in the ground using a
finite element analysis and then use these stresses in the stability analysis. This method is available in
2D, in which the stresses can be defined using either a static or dynamic stress-strain analysis. The
method no longer seeks to find a factor by which the shear strength must be reduced to bring the
system into equilibrium, but rather, it calculates the stability factor – which is equivalent to factor of
safety - directly as:

𝑆.𝐹. =
∑ 𝑆𝑟 Equation 69

∑ 𝑆𝑚
where 𝑆𝑟 and 𝑆𝑚 are the resisting and mobilized shear forces, respectively. Summation of these forces
over all columns yields the total resisting shear force and total mobilized shear force. The resisting force
at the base of each column, ignoring the unsaturated strength component, is calculated as (see Equation
13):

98
 𝑆𝑟 = 𝛽𝑖𝑠 = 𝛽𝑖(𝑐' + (𝜎𝑛 ‒ 𝑢𝑤)𝑡𝑎𝑛𝜙') Equation 70

where 𝛽𝑖 is the column area. Similarly, the mobilized shear force at the base of a column is calculated as:
𝑆𝑚 = 𝜏 𝑚𝛽 𝑖 Equation 71

where 𝜏𝑚 is the mobilized shear stress. Of significance is the fact that both the normal stress (𝜎𝑛) and
the mobilized shear stress are determined from a stress-strain analysis. Therefore, the equations for
computing the stability factors are linear; that is, no iteration is required to establish the stability factors
as in the limit equilibrium method.

In a stress-strain analysis, the cartesian stress stresses (𝜎𝑥, 𝜎𝑦, 𝜎𝑥) are computed at Gauss regions within
an element. Using the interpolating functions for the element, the stresses can also be computed at any
local point within an element. In this manner, the stresses are calculated at the centroid of every column
base and the Mohr-Circle equations are used to calculate the normal stress and shear mobilized in 2D
as:
𝜎𝑥 + 𝜎𝑦 𝜎𝑥 ‒ 𝜎𝑦 Equation 72
𝜎𝑛 = + cos 2𝜃 + 𝜏𝑥𝑦sin 2𝜃
2 2

Equation 73
𝜎𝑥 ‒ 𝜎𝑦
𝜏𝑚 = 𝜏𝑥𝑦cos 2𝜃 ‒ sin 2𝜃
2

where 𝜃 is the inclination of the column base. The line of application of the normal stress is
perpendicular to the column base, whereas the line of application of the mobilized shear stress is
parallel to the column base.

8.1.2 Establishing Finite Element Stresses

The stresses from any GeoStudio finite element analysis can be used in a stress-based stability analysis.
The most simplistic approach is to conduct a Gravity Activation analysis, which assumes all materials
isotropic and elastic. That stated, a more rigorous approach is required for most practical problems and
involves the following (see “Stress-Strain Modeling with GeoStudio” reference book):

1. Conduct an In Situ analysis to establish the initial stresses and pore water pressure in the ground
using the Gravity Activation, K0, or Field Stress procedure.
2. Conduct a Stress Correction to return stresses within a domain to legal stress space; that is, to a
location in stress space that is on or below the yield surface. These illegal stresses can occur
when an Isotropic Elastic material model is used to establish the initial stresses.
3. Conduct a load-deformation or consolidation analysis to model staged construction, soil-
structure interaction, and/or, consolidation.

99
A stress-based stability analysis can you use the stresses at any point in the modeling sequence,
including at multiple points in a staged construction analysis or at any time step within a consolidation
analysis.

8.1.3 Commentary on Stress-Based Method

The use of finite element computed stresses inside a limit equilibrium framework to assess stability has
many advantages including:

• There is no need to make assumptions about intercolumn forces.

• The stability factor is deterministic once the stresses have been computed, and
consequently, there are no iterative convergence problems.
• The issue of displacement compatibility is satisfied.
• The computed ground stresses better capture the physical system.
• Stress concentrations are indirectly considered in the stability analysis.
• Soil-structure interaction effects are readily handled in the stability analysis

In most cases, the limit equilibrium and stress-based method will produce similar factor of safety values.
As such, the stress-based approach is often convenient if a stress-strain analysis is already being used to
consider deformations and/or pore water pressures, the problem involves soil-structure interaction, or
the limit equilibrium analysis is suffering from convergence problems. A key benefit of the method is
that it is anchored to the familiarity of limit equilibrium methods routinely used in current practice.

9 Rapid Drawdown Analysis Methods

Stability analysis during rapid drawdown is an important consideration in the design of embankment
dams. During rapid drawdown, the stabilizing effect of the water on the upstream face is lost, but the
pore water pressures within the embankment may remain high. As a result, the stability of the upstream
face of the dam can be much reduced. The dissipation of pore water pressure in the embankment is
largely influenced by the permeability and the storage characteristic of the embankment materials.
Highly permeable materials drain quickly during rapid drawdown, but low permeability materials take a
long time to drain.

Stability during rapid drawdown can be modeled using two approaches, namely the “effective strength”
approach (simple and rigorous) and the “staged undrained strength” approach. In the effective strength
approach, the effective stress parameters are used to compute the effective shear strength of the
embankment during rapid drawdown. The advantage of this approach is that effective shear strength
parameters are used which are fundamentally better. The disadvantage of this approach is the extra
work required to estimate the pore water pressures within materials during the drawdown process.

In the staged undrained strength method, the shear strength is evaluated based on the total stress
parameters and consequently, the requirements associated with the estimation of the pore water
pressure conditions for effective stress analysis can be avoided. The disadvantage of this approach is
that total stress parameters are used, and these total stress parameters are not commonly available.
Furthermore, the hydraulic properties of the materials are ignored and the stability of the embankment
dam at different times during the drawdown process cannot be evaluated.

100
9.1 Simple Effective Strength Method
One simple and conservative approach is to assume that the rapid drawdown process occurs
instantaneously. The stabilizing force due to the upstream water is assumed to immediately vanish while
the piezometric surface follows the upstream slope and remains unchanged within the embankment
dam (Figure 123). It is the piezometric surface that is used to compute the pore water pressures along
the slip surface.

Piezometric Line
Drain
Clay Core

Embankment
Bedrock

Figure 123. Piezometric surface for a simple rapid drawdown analysis.

In general, the simple effective strength approach represents the worst-case scenario. This rapid
drawdown seldom occurs instantaneously and pore water pressures in the upstream material may
dissipate quite readily during the drawdown process.

9.2 Rigorous Effective Strength Method

A more sophisticated approach is to model the rapid drawdown and evaluate the pore water pressure
conditions by doing a finite element transient seepage analysis. The advantage of this approach is that
the hydraulic properties of the materials are considered, and a time component can be included in the
analysis. With this approach, rapid drawdown is not just instantaneous but can be modeled as a process.
The factor of safety of the embankment dam at different times during the entire drawdown process can
be evaluated as shown in Figure 124.

2.5

2
Factor of Safety

1.5

0.5

0
0 2 4 6 8 10 12 14 16 18

T ime Steps
.
Figure 124. Factor of safety versus time during rapid drawdown.

101
9.3 Staged Drawdown (Duncan et al., 1990)
Duncan, Wright and Wong (1990) proposed a three-stage procedure for analyzing the rapid drawdown
case. The procedure assures that the undrained strength, calculated prior to drawdown, does not
exceed the drained strength calculated post drawdown. The following inputs are required to complete
the analysis:
' '
 Effective stress strength parameters 𝑐 and 𝜙 (via Mohr-Coulomb or spatial Mohr-Coulomb soil
strength models).
 The intercept 𝑐𝑅and slope 𝜙𝑅 of the R (total stress) undrained strength envelope taken as the
tangent to Mohr circles plotted using the minor principal stress at consolidation and principal
stress difference at failure (diameter of circle).
 A piezometric surface defining the pore-water pressure condition prior to drawdown.
 A piezometric surface defining the pore-water pressure condition after drawdown.
 Drained materials should have 𝑐𝑅 = 0 and 𝜙𝑅 = 0. A verification check ensures that 𝑐𝑅 is greater
than 𝑐 and 𝜙𝑅 is less than 𝜙 for undrained materials.
' '

Two piezometric surface must be defined to complete a staged rapid drawdown analysis. The higher
surface characterizes the pore-water pressure conditions prior to drawdown. The lower surface
characterizes the pore-water pressure conditions after drawdown.

The publications should be consulted for the details of the procedure. The following is a synopsis of the
calculation stages:

Stage 1: The first stage involves a conventional limit equilibrium stability analysis using the long-term
effective stress strengths and pore-water pressure conditions prior to rapid drawdown. The computed
effective normal and shear stresses along the slip surface are retained. These effective stresses are
assumed to represent those that exist on a potential failure plane at the time of consolidation and are
referred to as the consolidation stresses. The consolidation stresses are used to estimate the undrained
shear strengths for materials that do not drain freely.

Stage 2: The second stage involves a stability analysis after drawdown when the water level is low. The
effective strength parameters are used for the freely drained materials. Undrained strengths are
calculated for the poorly drained materials from the computed effective normal and shear stresses from
Stage 1. Parameters that define the isotropic consolidation strength envelope are calculated from the
effective (𝑐 and 𝜙 ) and total (𝑐𝑅and 𝜙𝑅) strength parameters by:
' '

cos 𝜙𝑅𝑐𝑜𝑠𝜙' Equation 74

𝑑𝐾 = 𝑐𝑅
𝑐=1 1 ‒ sin 𝜙𝑅

sin 𝜙𝑅𝑐𝑜𝑠𝜙' Equation 75

𝜓𝐾 = 𝑎𝑟𝑐𝑡𝑎𝑛
𝑐=1 1 ‒ sin 𝜙𝑅

102
The isotropic consolidation strength envelope relates the strength on the failure plane to the effective
consolidation pressure; that is, the isotropic pressure on a triaxial specimen prior to shearing.
Subsequently, a procedure is used to interpolate the undrained shear strength between the effective
stress and isotropic consolidation shear strength envelopes. The procedure attempts to account for the
fact that anisotropically consolidated undrained strengths might be lower than isotropically
consolidated undrained strengths. The undrained strength at the base of each column is compared with
the drained strength, and the smaller strength is used in the third stage.

Stage 3: The third stage involves a stability analysis using the lesser of the effective stress strengths
corresponding to the low water level or the total stress strengths from Stage 2. The factor of safety
computed from the third stage analysis is used to represent the stability after rapid drawdown. The
factor of safety computed from Stage 3 will be the same as the factor of safety from Stage 2 if the
undrained strengths were less than the effective stress strengths.

A total of four rapid drawdown examples published by Duncan, Wright and Wong (1990) and by the
Corps of Engineers (2003) have been repeated GeoStudio (see table below). Figure 125 present the
GeoStudio result of the Walter Bouldin Dam.
Table 8. Comparison to publications.

Example GeoStudio Publication Reference

Walter Bouldin Dam 1.03 1.04 Duncan, Wright and Wong (1990)
Pumped Storage Project 1.55 1.56 Duncan, Wright and Wong (1990)
Dam
Pilarcitos Dam 1.05 1.05 Duncan, Wright and Wong (1990)
USACE benchmark 1.46 1.44 Corps of Engineers (2003) EM
1110-2-1902

103
Figure 125. Walter Bouldin Dam verification after rapid drawdown.

10 Seismic and Dynamic Stability Analysis

Dynamic forces due to seismic activity are usually oscillatory, multi-directional, and act only for
moments in time. Despite this complex response, static forces are sometimes used to represent the
effect of the dynamic loading in a limit equilibrium analysis. The dynamic effects can be considered in
several ways. The simplest is a pseudo-static type of analysis, which can be conducted in both 2D and 3D
analyses. A more complex approach is to use QUAKE/W finite element dynamic stresses and pore water
pressures together with SLOPE/W for a 2D analysis.

10.1 Pseudo-static Analysis

A pseudo-static analysis represents the effects of earthquake shaking by accelerations that create
inertial forces. These forces act in the horizontal and vertical directions at the centroid of each column.
The horizontal forces (𝐹ℎ) are defined as:
𝑎ℎ Equation 76
𝐹ℎ = 𝑊 = 𝑘 ℎ𝑊
𝑔

and the vertical forces (𝐹𝑣) are calculated as:

𝑎𝑣 Equation 77
𝐹𝑣 = 𝑊 = 𝑘𝑣𝑊
𝑔

where 𝑎ℎ and 𝑎𝑣 are the horizontal and vertical pseudo-static accelerations, 𝑔 is the gravitational

104
acceleration constant, 𝑊 is the weight of the column, and 𝑘ℎ and 𝑘𝑣 are the horizontal and vertical
seismic coefficients, respectively. In GeoStudio, a seismic load is defined by specifying the seismic
coefficient. A 𝑘ℎ of 0.2, for example, would mean that the horizontal acceleration is 0.2 𝑔. The surcharge
load from ponded water acting on a slope is not included in the calculation of the inertial forces. Water
has no shear strength and therefore inertial forces acting on the water do not contribute to destabilizing
the slope.

In a 3D stability analysis, the horizontal seismic force is applied on each column in the direction of
sliding. In both 2D and 3D, the vertical seismic load is added to the column weight. Consider the column
free body diagram presented in Figure 126 and Figure 127. In the first case, a 𝑘ℎ of 0.2 yields a
horizontal force 22.7 kN on a column with a weight of 113.48 kN. In the second case, a 𝑘𝑣 of 0.1 is
included, producing a vertical seismic load of 11.35 kN, which yields a total column weight of 124.83 kN
(113.48 kN + 11.35kN).

Slice 10 - Morgenstern-Price Method

113.48
148.43

309.69 336.22
22.697

173.45
41.097

78.447

Figure 126. Column FBD with a horizontal seismic load.

Slice 10 - Morgenstern-Price Method

124.83
159.33

336.41 364.63
22.697

185.89
44.589

87.766

Figure 127. Column FBD with a horizontal and vertical seismic load.

105
The application of vertical seismic coefficients often has little impact on the factor of safety because the
vertical inertial forces alter the column weight, which in turn alters the column base normal and shear
resistance. The added mobilized shear arising from the added weight tends to be offset by the increase
in shear strength. This is only true for frictional strength components and not for cohesive or undrained
strength.

Horizontal inertial seismic forces can have a dramatic effect on the stability of a slope. Even relatively
small seismic coefficients can lower the factor of safety greatly, and if the coefficients are too large, it
becomes impossible to obtain a converged solution. It is consequently good practice to apply the seismic
forces incrementally to gain an understanding of the sensitivity of the factor of safety to this parameter.
It is often useful to create a graph such as in Figure 128. As the seismic coefficient increases, there
should be smooth gradual decrease in the factor of safety.

1.70

1.60

1.50

1.40
Factor of Safety

1.30

1.20

1.10

1.00

0.90

0.80
0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35
Horizontal K coefficient

Figure 128. Factor of safety vs the horizontal seismic coefficient.

10.2 Undrained Response during Dynamic Loading

In some cases, dynamic loading is so rapid that the soil behaves in an undrained manner. This may be
the case for saturated clayey soils, for example, but not for unsaturated gravels or rock fill. As discussed
in the theory section, a limit equilibrium formulation determines the forces acting on a column such that
the column is in force and moment equilibrium and the factor of safety is the same for every column.
The application of a dynamic force to a column naturally alters the entire equilibrium and therefore the
normal force and shear strength at the column base. In GeoStudio, there are two Staged Pseudo-Static
options that calculate the strength along the slip surface prior to including seismic forces: the Effective
Stress Strength approach and the Undrained Strength method by Duncan et al. (1990). Both these
techniques are only applicable to 2D analysis.

106
10.2.1 Effective Stress Strength
In this case of an Effective Stress Strength Pseudo-static analysis, the analysis is first solved without any
seismic forces to establish the effective stress and shear strength along the slip surface (stage 1). The
shear strength at the base of each column is then converted into an equivalent undrained (cohesive)
strength. The analysis is then repeated with the seismic loading. The undrained shear strength is not a
function of the normal stress and, consequently, the seismic loads do not alter the shear strength. This is
illustrated in Figure 129. The graph on the left shows the frictional strength before the seismic forces are
applied. The graph on the right shows that the strength was converted into an equivalent cohesive
strength. Note that the strengths are the same in both graphs.

80 80

Cohesive Cohesive
60 60
Strength

Strength

40 40

20 20
Frictional Frictional

0 0
0 10 20 30 40 0 10 20 30 40
Distance Distance

Figure 129. Conversion of frictional strength to cohesive strength in a seismic analysis.

10.2.2 Undrained Strengths (Duncan, Wright, and Wong, 1990)

As above, the analysis is first solved without any seismic forces to establish the effective stress and shear
strength along the slip surface (stage 1). The undrained strength is then calculated by the equations
proposed by Duncan et al. (1990), which make use of the R-envelope strength properties in the material
model definition (See Section 9.3). The effective stress strengths are computed during the second stage
of the analysis when the seismic forces are applied. Finally, the smaller of the undrained or effective
stress strength at the base of a column is used to calculate the factor of safety.

Seismic K = 0.20

Figure 130. Illustrative example with seismic forces.

107
 70

60

50
Shear Strength

40

30

20

Without seismic forces With seismic forces

10

0
0 5 10 15 20 25 30 35

Distance along slip surface

Figure 131. Effect of seismic forces on altering the shear strength.

11 Newmark Permanent Deformation

The Newmark method approximates the deformation of a 2D sliding mass due to earthquake loading.
Figure 132 presents an example of a 10-meter slope. The slope is subjected to 10 seconds of shaking
resulting from the earthquake record shown in Figure 133. Figure 134 shows the deformed mesh at two
steps in the earthquake history.
25

Example Earthquake Analysis

20
Elevation (m)

15

10

0
0 10 20 30 40 50 60 70

Distance (m)

Figure 132. Slope subjected to earthquake shaking.

108
 0.30
0.25
0.20
0.15
Acceleration ( g )

0.10
0.05
0.00
-0.05
-0.10
-0.15
-0.20
-0.25
0 2 4 6 8 10
Time (sec)

Figure 133. Example earthquake record.

Figure 134. Deformation at two steps in the earthquake cycle.

Section 8.1.1 presents the stress-based stability procedure, which is the basis of the Newmark method.
The Newmark procedure inherently assumes that there is no change in the material strength during
shaking; therefore, the dynamic and static shear resistances are equal and Equation 69 becomes:

𝑆.𝐹. =
∑𝑆𝑟(𝑠𝑡𝑎𝑡𝑖𝑐) Equation 78

∑ 𝑆𝑚
where 𝑆𝑟(𝑠𝑡𝑎𝑡𝑖𝑐) is the resisting shear force calculated from the initial static stresses {𝜎𝑠𝑡𝑎𝑡𝑖𝑐}. The
mobilized shear force 𝑆𝑚 at any time step is calculated from the accumulated stresses {𝜎} at the base of
the column. Figure 135 presents the factor of safety versus time for the slip surface that was deemed
most critical. The graph vividly illustrates the oscillation in the factor of safety due to the shaking. Note
that the factor of safety momentarily falls below 1.0 at several times during the shaking. The Newmark
method assumes that irreversible deformation occurs during those periods.

109
 1.4

1.3

1.2
Factor of Safety

1.1

1.0

0.9

0.8
0 2 4 6 8 10
Time

Figure 135. Factor of safety versus time during shaking.

The average acceleration of a sliding mass can be calculated from Newton’s second law of motion:

𝑎𝑎𝑣𝑒 =
∑𝑆𝑚 ‒ ∑𝑆𝑚(𝑠𝑡𝑎𝑡𝑖𝑐) Equation 79

∑𝑊
where the numerator comprises the difference between the current and initial (static) mobilized shear
force and the denominator is the summation of the column weights; that is, the total mass. Figure 9
presents the average acceleration versus time for the slip surface that was deemed most critical.
0.15

0.10
Average Acceleration

0.05

0.00

-0.05

-0.10

-0.15
0 2 4 6 8 10
Time

Figure 136. Average accelerations versus time for a slip surface.

Figure 137 presents the stress-based factor of safety (Figure 135) plotted against the average
acceleration at each time step of the dynamic analysis. The average acceleration corresponding to a
factor of safety of 1.0 is called the yield acceleration, 𝑎𝑦. It is the overall average acceleration that will
cause the slide mass to move. In this example, 𝑎𝑦 is 0.048.

110
 1.4

1.3

1.2
Factor of Safety

1.1

1.0

0.9

0.8
-0.15 -0.10 -0.05 0.00 0.05 0.10 0.15
Average Acceleration

Figure 137. Factor of safety versus average acceleration for a slip surface.

The velocity of the sliding mass is calculated over each time interval, starting when 𝑎𝑎𝑣𝑒 > 𝑎𝑦 by:

𝑡 + 𝑑𝑡
Equation 80
𝑉= ∫ (𝑎𝑎𝑣𝑒 ‒ 𝑎𝑦)𝑑𝑡
𝑡

where the end of the time interval (𝑡 + 𝑑𝑡) corresponds to the moment when the velocity 𝑉 decreases
back to zero, and the sliding mass stops moving. This means the integration continues over a portion of
the time interval where 𝑎𝑎𝑣𝑒 < 𝑎𝑦. Integrating the area under the velocity versus time curve gives the
cumulative displacement (Figure 139).
0.20

0.15
Velocity

0.10

0.05

0.00
0 2 4 6 8 10
Time

Figure 138. Velocity versus time of a slip surface.

111
 0.10

0.08
Deformation

0.06

0.04

0.02

0.00
0 2 4 6 8 10
Time

Figure 139. Cumulative displacement versus time.

The results of a 2D dynamic stress-strain analyses are independent of the orientation of the domain
relative to the sign on the earthquake accelerations. In contrast, the Newmark solution will change if the
domain is mirrored or if the earthquake acceleration history is inverted (i.e., multiply each value by -1).
Consider, for example, that the positive peak acceleration would displace a left-to right sliding mass
downslope. Inverting the record, and therefore changing the sign on the peak acceleration, would
displace the sliding mass upslope at the same instant in time. To partially compensate for this nuanced
deficiency, the software inverts the simulated average acceleration history (Figure 136), repeats the
double integration with the same 𝑎𝑦, and retains the results corresponding to the largest displacement.
There is no guarantee that this procedure will obtain the largest possible displacement, so a better
practice would involve cloning the analyses and solving with an inverted earthquake acceleration record.

GeoStudio calculates the displacement for each trial slip surface, making it possible to find the trial slip
surface with the largest potential movement, which may not correspond to the slip surface with the
lowest factor of safety. The trial slip surfaces can be sorted by displacement rather than by factor of
safety, to identify the one with the largest movement.

The Newmark procedure for estimating the permanent deformation is based on the concept of a sliding
block analysis. A Newmark type of analysis is applicable only in certain specific circumstances. The
procedure is ideal for cases where there is little or no degradation in the soil shear strength during the
dynamic shaking. This may be the case for the rock-fill shells of an earth embankment, steep natural
rock slopes or unsaturated, heavily over-consolidated clayey slopes. This procedure is not applicable to
cases where there is a significant potential for a large loss in shear strength due to either the generation
of excess pore-water pressures or the collapse of the soil grain structure as may be the case for loose
silty and sandy soils. Kramer (1996, p. 437) points out that these types of analyses can be unreliable for
situations where there is a large build up of excess pore-water pressures or where there is more than
about a 15% degradation in the shear strength during the earthquake shaking.

112
12 Probabilistic Analysis
A probabilistic analysis requires a probability density function (PDF) for one or more parameters (Figure
140). The abscissa of the PDF is the random offset from a user input parameter (e.g., cohesion).
Numerical integration of the PDF produces the cumulative distribution function (CDF; Figure 140). The
Monte Carlo sampling method randomly generates a number between [0, 1] and obtains the
corresponding offset from the CDF (Figure 140). Figure 140 shows a sample with 5 data points (i.e., 5
trials). Random samples are similarly generated for each parameter having a PDF and the offsets are
added to the deterministic input values to obtain the combination of parameters used for each trial run
(i.e., each analysis).
1.6 1

0.9
1.4

0.8
1.2
0.7
Probability Density

1
0.6

Probability
0.8 0.5

0.4
0.6

0.3
0.4
0.2

0.2
0.1

0 0
-1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1
Offset Offset

Figure 140. Probability density function (left) and cumulative distribution function (right) with 5 random Monte-Carlo
samples.

The Monte-Carlo sampling method tends to cluster the samples around the mean; a problem that is
accentuated by a steep CDF. This deficiency can be mitigated by using a larger sample size (i.e., number
of trials). Alternatively, the Latin Hypercube sampling method can be used to sample more evenly across
the CDF. The Latin Hypercube method is a stratified Monte Carlo method that involves dividing the
cumulative density function (CDF) into equal partitions and then calculating a data point in each
partition. The data points in each partition can be chosen randomly; however, GeoStudio adopts a
systemic mid-point sampling. For example, 5 data points sample discussed previously involves dividing
the CDF into 5 equal partitions and calculating one point in each partition (Figure 141). The sampling of a
single variable is referred to as one-dimensional Latin hypercube sampling. Multi-dimensional Latin
hypercube sampling involves generating one-dimensional samples for each variable and then combining
the lists, randomly, into an array. For example, Figure 142 shows an example of randomized partitions
from which 10 samples are obtained for 4 parameters (i.e., a, b, c, and d). The number of samples is
equal to the number of trials specified by the user.

113
 1.6 1

0.9
1.4

0.8
1.2
0.7
Probability Density

1
0.6

Probability
0.8 0.5

0.4
0.6

0.3
0.4
0.2

0.2
0.1

0 0
-1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1
Offset Offset

Figure 141. Probability density function (left) and cumulative distribution function (right) with 5 mid-point Latin Hypercube
samples.

Paramater (Partition)
Sample a b c d
1 9 1 1 2
2 3 9 10 6
3 6 5 7 1
4 7 2 2 8
5 10 4 9 5
6 1 3 3 4
7 5 6 8 10
8 2 7 4 3
9 4 10 5 7
10 8 8 6 9

Figure 142. Randomized partitions for 10 samples of parameters a, b, c, and d.

GeoStudio labels the Probability axis of the CDF ‘Random Number’ to clearly communicate its role in a
probabilistic analysis. In addition, the axes are swapped to communicate that a number X between [0, 1]
is used to obtain an offset 𝑌 = 𝑓(𝑋) (Figure 143).

114
Figure 143. Example of an offset sample function.

12.1 Probability Density Functions

The types of probability density functions (PDF) that can be specified include Normal, Lognormal,
Uniform, Triangular, and Generalized Spline.

The normal probability density function is given by:

‒ [(𝑥 ‒ 𝜇)2/2𝜎2]
𝑒 Equation 81
𝑓(𝑥) =
𝜎 2𝜋

where 𝑥 is the value for which the probability density is required, 𝜇 is the mean, 𝜎 is the standard
deviation. The probability density function 𝑓(𝑥) is plotted against the offset 𝛿𝑥, which is calculated as:
Equation 82
𝛿𝑥 = 𝑥 ‒ 𝜇 + 𝛿

where 𝛿 is the offset mean. The user inputs are 𝜎 and the offset mean 𝛿. The mean (𝜇) is assumed equal
to the deterministic input value. Figure 144 shows a typical normal distribution for a parameter 𝑥 with
𝛿 = 0, 𝜎 = 2 , and the deterministic input value, which is assumed equal to the means, of 𝜇 = 10.

115
Figure 144. Normal probability density function.

By default, the PDF extends to ± 5 standard deviations from the mean offset. Given a standard of
deviation of 2, the offset will range from -10 to +10. If, for example, the parameter 𝑥 was effective
friction angle with a mean of 30, then the randomly sampled friction angles would vary between 20° to
40°. Alternatively, the minimum and maximum values of the function can be specified. Figure 145 shows
the probability density function when clipped with a minimum/maximum offset of ± 5.

116
Figure 145. Clipped normal probability density function.

A variable 𝑥 is lognormally distributed if the natural log of 𝑥 (i.e., ln 𝑥) is normally distributed. The
general lognormal probability density function is given by:

‒ ((ln (𝑥) ‒ 𝑀)2/2𝑆2)

𝑒 Equation 83
𝑓(𝑥) =
𝑥𝑆 2𝜋

where 𝑆 is the lognormal standard deviation and 𝑀 is the lognormal mean, which are calculated from
the input standard deviation (𝜎) and mean (𝜇) by:

𝑆 = 𝑙𝑛 ((𝜎𝜇) + 1)
2 Equation 84

1 Equation 85
𝑀 = 𝑙𝑛(𝜇) ‒ 𝑆2
2

The probability density function 𝑓(𝑥) is plotted against the offset 𝛿𝑥, which is calculated as:
Equation 86
𝛿𝑥 = 𝑥 + 𝛿

where 𝛿 is the offset mean. The user inputs are 𝜇, 𝜎, and the offset mean 𝛿. Contrary to Equation 82, the
offset mean should be set equal to, or close to, the mean. Figure 146 shows a lognormal distribution for
𝜇 = 10, 𝜎 = 2, and the offset mean 𝛿 =‒ 10 (compare with Figure 144). The corresponding sampling
function is shown in Figure 147.

117
Figure 146. Lognormal probability density function.

Figure 147. Lognormal sampling function.

118
A uniform function simply describes a range of values as illustrated in Figure 148. The characteristic of a
uniform distribution is that all values have an equal probability of occurrence. This is further illustrated
by the fact that the sampling function is a straight line as in Figure 149.

Figure 148. A uniform probability density function.

Figure 149. A uniform sampling function.

119
A triangular probability density function can be specified with three points: the minimum and maximum
offsets and the offset of the apex value. Figure 150 shows a typical triangular probability density
function with the minimum and maximum offset at -5 and 25, respectively, and the apex located at an
offset of 0. Figure 151 shows the associated sampling function.

Figure 150. A triangular probability density function.

120
Figure 151. A triangular sampling function.

The probability distribution offset function can also be defined using a generalized spline function. As
with the normal and lognormal distributions, the abscissa of the PDF must be defined using an offset.
Figure 152 shows a bimodal PDF define by 5 data points that varies the parameter from its deterministic
value by as much as ± 10. The corresponding sampling function is shown Figure 153.

Figure 152. A probability density function defined using spline data points.

121
Figure 153. A generalized spline sampling function.

'
12.2 c – 𝜙 Correlation
A correlation coefficient expresses the relative strength of the association between two parameters.
Laboratory tests on a wide variety of soils (Lumb, 1970; Grivas, 1981 and Wolff, 1985) show that the
'
shear strength parameters 𝑐 and 𝜙 are often negatively correlated with correlation coefficient ranges
from -0.72 to 0.35. Correlation between strength parameters may affect the probability distribution of a
'
slope. GeoStudio allows the specification of 𝑐 and 𝜙 correlation coefficients for all soil models using 𝑐
'
and 𝜙 parameters.

Correlation coefficients will always fall between -1 and 1. When the correlation coefficient is positive, 𝑐
'
and 𝜙 are positively correlated, implying that larger values of 𝑐 are more likely to occur with larger
' '
values of 𝜙 . Similarly, when the correlation coefficient is negative, 𝑐 and 𝜙 are negatively correlated
'
which reflects the tendency of a larger value of 𝑐 to occur with a smaller value of 𝜙 . A zero-correlation
'
coefficient implies that 𝑐 and 𝜙 are independent parameters.
'
The correlation coefficient can be specified for all material models using 𝑐 and 𝜙 . The input is visible
when cohesion is selected as the probabilistic parameter. The correlation is considered via the random
number generated to sample the cohesion function using:

𝑁𝑎 = 𝑁1𝑘 + (1 ‒ |𝑘|)𝑁2 Equation 87

where 𝑁𝑎 is the adjusted normalized random number for the cohesion, 𝑘 is the correlation coefficient
between 𝑐 and 𝜙 , 𝑁1 is the normalized random number for 𝜙 , and 𝑁2 is the normalized random
' '

number for the 𝑐.

122
12.3 Probability of failure and reliability index
A factor of safety is really an index indicating the relative stability of a slope. It does not imply the actual
risk level of the slope due to the variability of input parameters. With a probabilistic analysis, two useful
indices are available to quantify the stability or the risk level of a slope. These two indices are known as
the probability of failure and the reliability index.

As illustrated in Figure 154, the probability of failure is the probability of obtaining a factor of safety less
than 1.0. The probability of failure is determined by counting the number of safety factors below 1.0 and
then taking this number as a percentage of the total number of converged Monte Carlo trials. For
example, if there are 1000 Monte Carlo trials with 980 converged safety factors and 98 of them are
below 1.0, then the probabilistic of failure is 10%.
Probability Distribution Function
100

80 P (F of S < x)
Probability (%)

60

40

20 P (Failure)

0
0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7

Factor of Safety

Figure 154. Factor of safety probability of failure.

The probability of failure can be interpreted in two ways (:

 If a slope were to be constructed many times, what percentage of such slopes would fail; or
 The level of confidence that can be placed in a design.

The first interpretation may be relevant in projects where the same slope is constructed many times,
while the second interpretation is more relevant in projects where a given design is only constructed
once and it either fails or it does not. Nevertheless, the probability of failure is a good index showing the
level of risk of instability.

There is no direct relationship between the deterministic factor of safety and probability of failure. In
other words, a slope with a higher factor of safety may not be more stable than a slope with a lower
factor of safety (. For example, a slope with factor of safety of 1.5 and a standard deviation of 0.5 may

123
have a much higher probability of failure than a slope with factor of safety of 1.2 and a standard
deviation of 0.1.

Another way of looking at the risk of instability is with what is known as a reliability index. The reliability
index, 𝛽, is defined in terms of the mean and the standard deviation of the trial factors of safety as
shown in the following equation:

(𝜇 ‒ 1.0) Equation 88
𝛽=
𝜎

The reliability index describes the stability by the number of standard deviations separating the mean
factor of safety from its defined failure value of 1.0. It can also be considered as a way of normalizing the
factor of safety with respect to its uncertainty. It is applicable to normal distributions.

When the shape of the probability distribution is known, the reliability index can be related directly to
the probability of failure. Figure 155 illustrates the relationship of the reliability index to the probability
of failure for a normally distributed factor of safety.

Figure 155. Relationship between reliability index and the probability of failure.

12.4 Spatial variability

Spatial variability can be an important issue if a potential slip surface has a significant length within one
material. Consider the case in Figure 156. A significant portion of the slip surface is in the clay and
consequentially it is reasonable to assume that there is some statistical variability in the clay along the
slip surface. The length of the potential slip surface in the clay is around 25 m. The horizontal distance is
24 m.

124
 Sand

Clay
24 m

Figure 156. Case with statistical variability in the clay.

The clay has a mean undrained strength (𝑐) of 15 kPa with a standard deviation of ±3. A clipped normal
distribution having an offset of -10 kPa and +10 kPa is assumed as shown in Figure 157.

There are three options for considering the spatial variability along the slip surface. The options deal
with the number of times the soil properties are sampled along the slip surface. The options are:

 Sample the properties only once for each soil for the entire slip surface.
 Sample the properties for each column along the slip surface.
 Sample the properties at a specified distance along the slip surface.

If the clay strength is sampled only once for each Monte Carlo trial, the strength is constant along the
slip surface, as illustrated in Figure 158. The graph shows possible strengths for four different Monte
Carlo trials.

Figure 157. Probability density distribution for the clay.

125
 25.0

23.0

21.0

19.0

17.0 Trial 1
Cohesion

Trial 2
15.0
Trial 3
13.0 Trial 4

11.0

9.0

7.0

5.0
0.0 5.0 10.0 15.0 20.0 25.0

Distance along slip surface

Figure 158. Strength variations in the clay when sampled once for each trial run.

Statistically, the strength in the clay could be as low as 5 kPa or as high as 25 kPa. Intuitively, it seems
unlikely that the strength along the entire slip surface could be at these outer limits, particularly if the
slip surface has a significant length. Sampling more than once is likely more appropriate.

Another option is to sample the soil strength for each column. The variation along the slip surface then
could be as illustrated in Figure 159 for two trial runs. There are 10 columns in the clay and,
consequently, there are 10 different strengths along the slip surface.

21.0

19.0

17.0

15.0
Cohesion

Trial 1
13.0
Trial 2

11.0

9.0

7.0

5.0
0.00 5.00 10.00 15.00 20.00 25.00

Distance along slip surface

Figure 159. Strength variation along the slip surface when soil is sampled for every column.

126
As indicated on Figure 159, to sample the soil strength for each column in every Monte Carlo trial may
result in significant fluctuation of soil strength along the slip surface even within the same material.
Although this is statistically possible, it seems quite unlikely in the field.

The third option is to sample the soil at a specified distance. For example, if the sampling distance is 10
m, then the strength variation along the slip surface for two different Monte Carlo runs could be as
shown in Figure 160.

23.0

21.0

19.0

17.0
Cohesion

15.0
Trial 1
Trial 2
13.0

11.0

9.0

7.0

5.0
0.00 5.00 10.00 15.00 20.00 25.00

Distance along slip surface

Figure 160. Strength variation along the slip surface with a sampling distance of 10 m.

At this sampling distance, there are three different strengths along the slip surface. There are two
complete 10-m sampling distances and a third partial sampling distance. A partial sampling distance is
considered to be correlated with the preceding sampling distance. The coefficient of correlation
between two soil sections can be computed with the equation as proposed by Vanmarcke (1983):

𝑍02Γ(𝑍0) ‒ 𝑍12Γ(𝑍1) + 𝑍22Γ(𝑍2) ‒ 𝑍32Γ(𝑍3) Equation 89

𝜌(∆𝑍,∆𝑍') =
0.5
2Δ𝑍Δ𝑍'[Γ(Δ𝑍)Γ(Δ𝑍')]

where:
∆𝑍,∆𝑍’ is the length between two sections;

Z0 the distance between the two sections;

Z1 is ∆𝑍 + 𝑍0;

Z2 is ∆𝑍 + 𝑍0 + ∆𝑍';

Z3 is ∆𝑍' + 𝑍0; and

127
Γ is a dimensionless variance function.

Vanmarcke (1983) showed that the dimensionless variance function can be approximated by:
Γ(𝑍) = 1.0 when 𝑍 ≤ 𝑑; and
Γ(𝑍) = 𝑑/𝑍 when 𝑍 > 𝑑.

The symbol 𝑑 is called the scale of fluctuation. It is about twice the autocorrelation distance. When
spatial variation is considered, the value of 𝑑 is related to the desired sampling distance along the slip
surface.

As is evident from the above equations, the correlation is strong when the partial distance is short
relative to the specified sampling distance, and the correlation is weak when the partial distance
approaches the specified sampling distance.

The probability of failure is related to the sampling distance, especially if there are only a few statistical
parameters. The probability of failure is the highest when the soil is sampled only once for each trial run
and is the lowest when the soil is sampled for each column. For the simple case under discussion here,
the probabilities of failure for various sampling frequencies are listed in Table 9.
Table 9. Probability of failure for various sampling frequencies.

Sampling Frequency Probability of Failure (%)

Only once 19
Every 15 m 10
Every 10 m 6
Every column <1

There are more trial runs with a factor of safety less than 1.0 when the strength is sampled only once
and therefore a high probability of failure. Sampling the soil for every column tends towards an average
mean strength and therefore a low probability of failure.

Selecting an appropriate sampling distance is not trivial. Generally, a thorough understanding of the
statistical characteristics of the soil stratum is required, and spatial variability is a factor that can only be
assessed considering the details at each project and the amount of data available.

The two options discussed above – sampling the soil only once or sampling the soil strength for each
column – are perhaps unrealistic for some situations; however, these two options are nonetheless
useful. They provide a convenient means of bracketing the possible range of probabilities of failure.
Knowing the limits can then help with selecting and evaluating a proper sampling distance.

If sufficient data is available, a variogram can be created as discussed by Soulié et al. (1990). The form of
a variogram is illustrated in Figure 161. The differences in strength, for example, increase with distance
between sampling points up to a certain distance after which the difference remains relatively constant.
The point at which the variogram becomes horizontal is called the range, and it is the distance beyond
which there is no correlation between strength measurements. The range from a variogram would be an
appropriate sampling distance for a probabilistic analysis.

128
 Range
between sampling points
Strength difference

Distance between sampling points

Figure 161. The form of a typical variogram.

12.5 Number of Monte Carlo Trials

The number of Monte Carlo trials is dependent on the number of variable input parameters and the
expected probability of failure. In general, the number of required trials increases as the number of
variable input increases, or the expected probability of failure becomes smaller. It is not unusual to do
thousands of trials to achieve an acceptable level of confidence in a Monte Carlo probabilistic slope
stability analysis.

It has been suggested that the number of required Monte Carlo trials is dependent on the desired level
of confidence in the solution, as well as the number of variables being considered. Statistically, the
following equation can be used (Harr, 1987):

𝑐2
𝑁𝑚𝑐 =
[ 4 (1 ‒ 𝜀 ) 2 ] 𝑚 Equation 90

where 𝑁𝑚𝑐 is the number of Monte Carol trials, 𝑐 is the desired level of confident (0 to 100%) expressed
as a decimal, 𝜀 the normal standard deviation corresponding to the level of confidence, and 𝑚 the
number of variables. The number of Monte Carlo trials increases geometrically with the level of
confidence and the number of variables. For example, if the desired level of confidence is 80%, the
normal standard deviation will be 1.28; the number of Monte Carlo trials will be 10 for 1 variable, 100
for 2 variables and 1,000 for 3 variables. For a 90% level of confidence, the normal standard deviation
will be 1.64; the number of Monte Carlo trials will be 67 for 1 variable, 4,489 for 2 variables and 300,763
for 3 variables. In fact, for a 100% level of confidence, an infinite number of trials will be required.

For practical purposes, the number of Monte Carlo trials to be conducted is usually in the order of
thousands. This may not be sufficient for a high level of confidence with multiple variables; fortunately,
in most cases, the solution is not very sensitive to the number of trials after around 3000 trials have
been run. Mean and standard deviation of the simulated factor of safety can be plotted again the run
number to discern if the number of trials was adequate.

129
12.6 Multiple statistical parameters
A wide range of parameters can be assigned a statistical dispersion and included in a Monte Carlo
simulation. This includes soil strength parameters, concentrated loads, pore water pressures, pseudo-
static seismic coefficients and so forth. Each of the parameters is sampled for each Monte Carlo run.
Using too many statistical parameters together in the same run tends to reduce the distribution of the
computed safety factors. That is, the combination of multiple statistical parameters tends to smooth out
the extreme of the variability and the probability of failure diminishes. Perhaps this is an appropriate
response. The main point here is that it is important to be cognizant of this and to evaluate the results
carefully when many parameters are included in a probability analysis.

13 Sensitivity Analysis
A sensitivity analysis is somewhat analogous to a probabilistic analysis. Instead of selecting the variable
parameters randomly, the parameters are selected in an ordered fashion using a Uniform Probability
Distribution function. Consider the simple example in Figure 162. The goal of this analysis may be to
determine if the stability is more sensitive to the frictional strength of the sand or the undrained
strength of the clay.

Sand

Clay

Figure 162. Illustrative case for a sensitivity analysis.

Consider an example in which the friction angle for the sand is 30° and the variation in the friction angle
will range from 25° to 35° degrees. The undrained strength of the clay is 20 kPa and the range will be
from 15 to 25 kPa. First, the cohesive strength of the clay is held constant at 20 kPa and the factor of
safety is computed for each friction angle value. Next the friction angle of the sand is kept constant and
factors of safety is computed for each cohesion value.

For presentation purposes, all the strengths are normalized to a value ranging between 0.0 and 1.0. Zero
means the lowest value and 1.0 means the highest value. The point where the two sensitivity curves
'
cross is the deterministic factor of safety. In this case, zero means a 𝜙 of 25° and a 𝑐 of 15 kPa. A value
'
of 1.0 means a 𝜙 of 35 degrees and a 𝑐 of 25 kPa. The results are presented as a sensitivity plot (Figure
163). For this case, the stability is much more sensitive to changes in the clay cohesion than to changes
in the sand friction angle. This is intuitively correct since a major portion of the slip surface is in the clay
layer.

130
 1.5

1.4
Phi

1.3
Factor of Safety

1.2

1.1

1.0 C

0.9
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Sensitivity Range

Figure 163. Illustration of a sensitivity plot.

14 Partial Factors
Partial factors are valid for 2D and 3D analyses. These factors can be applied to characteristic material
strength properties, actions, or resistances in accordance with a limit state design approach such as
Eurocode 7, Norwegian Standard NS 3480, and British Standard 8006. The factored values are in turn
used in the stability analysis to check if an ultimate limit state has been reached. The ultimate limit state
is defined as a state associated with complete collapse or failure. The calculated factor of safety is
interpreted as an over-design factor when partial factors are applied.

Material strength properties include effective stress or undrained strength properties. Actions refer to
gravity forces, seismic forces, and point and surcharge loads. Gravity forces are factored via the unit
weight of the soil and seismic forces via the seismic coefficients. Point and surcharge loads must be
declared by the user to be permanent or variable. The software determines if an action is favorable or
unfavorable as discussed subsequently. Resistances include soil strengths (as opposed to material
strength parameters) and reinforcement resistances.

14.1 Favorable or unfavorable actions

An action is deemed to be favorable or unfavorable by considering the angle of the action relative to the
angle of the column base and relative to the direction of movement of the sliding mass. Figure 164
illustrates an action on an individual 2D column. Actions that resist and promote movement of the
column are considered favorable and unfavorable, respectively. The action shown would be deemed
favorable for a right-facing slope and unfavorable for a left-facing slope. Actions oriented at 𝛿±90° are
deemed favorable if acting to increase the base normal force and unfavorable if otherwise.

131
Figure 164. Favorable/unfavorable action applied to a column.

Gravity forces are vertical and are therefore always deemed unfavorable for columns in the active zone.
Figure 165 compares the results of an analysis completed using a) characteristic values of unit weight;
and b) derived values of unit weight calculated by multiplying the unit weight of unfavorable columns by
a partial factor of 1.25. Columns 11 to 25 were deemed unfavorable, with the transition occurring
between Column 10 and 11 where the column base angle changed sign. The increased column weights
in the active zone resulted in an over-design factor of 1.29 as opposed to the overall factor of safety of
1.44. Incidentally, the increased gravity forces in the active zone cause an increase in the base normal
force, which in turn increases the frictional shear resistance if effective stress strengths are being used.
It is therefore possible for the over-design factor to exceed the overall factor of safety for some cases.

Figure 165. Result of applying a partial factor to the unit weight of unfavorable columns.

132
14.2 Pore water pressures and ponded water
Limit state design approaches sometimes allow for calculation of derived pore water pressures and
water loads by application of partial factors. GeoStudio, however, does not apply partial factors to pore
water pressures or to water loads because of the possibility for illogical force distributions. Design
values of pore water pressures and water loads should be derived by manually applying a safety factor
to the piezometric line(s) or to the boundary conditions that control the porewater pressures in a finite
element analysis. Consideration should be given to the permanent and variable components of the
hydraulic regime when applying the safety factor.

Pore water pressures are influenced by partial factors if the Ru or B-bar methods are used because of
the implicit relationship to the unit weight of the soil (See Sections 8.3 and 8.4). Figure 166 compares
the results of an analysis using a piezometric line with B-bar and completed using a) characteristic values
of unit weight; and b) derived values of unit weight calculated by multiplying the unit weight of
unfavorable columns by a partial factor of 1.25. The upper soil layer adds the weight, and the lower soil
layer has a B-bar of 1.0. Columns 4 through 14 are in the active zone and were therefore deemed
unfavorable. Factoring the column weights in the active zone by 1.25 resulted in the excess pore water
pressures being 1.25 greater than those produced using characteristic values. The pore water pressures
calculated by the Ru and B-bar methods can be thought of as being factored in the same manner as the
unit weights.

Figure 166. Excess pore water pressures resulting applying a partial factor to the unit weight of unfavorable columns.

133
14.3 Shear strength
Table 10 summarizes the way material model strengths are factored. The categories refer to: (1)
effective cohesion, 𝑐'; (2) effective coefficient of friction, 𝑡𝑎𝑛𝜙 ; (3) undrained strength, 𝑠𝑢; and, (4)
'

'
overall strength. For example, 𝑐 and 𝑡𝑎𝑛𝜙 are factored for Mohr-Coulomb, while the overall strength is
factored for Hoek-Brown. The effective stress strength models shown in the table also comprise an
optional suction strength component. The suction strength is calculated as a function of a coefficient of
friction value; consequently, suction strength is implicitly part of category (2). An Earth Resistance
partial factor can also be specified separately of the 4 material parameter factors, resulting in division of
the entire numerator of the force and moment factor of safety equations by a single value. An
equivalent result could be obtained by multiplying the factors in categories 1 through 4 by the required
resistance factor while leaving the Earth Resistance factor specified as 1.0.
Table 10. Application of partial factors to material model parameters and/or strength.

Material Model Material Parameter Category

Mohr-Coulomb 1,2
Undrained 3
S=f(depth) 3
S=f(datum) 3
High strength N.A.
Bedrock N.A.
Bilinear 4
Anisotropic strength 1,2
Anisotropic function 1,2
Shear-Normal function 4
Combined, S=f(depth) 1,2,3
Combined, S=f(datum) 1,2,3
S=f(overburden) 3
Spatial M-C 1,2
Generalized Hoek-Brown 4

Staged rapid drawdown and staged pseudo-static analyses involve the calculation of both drained
effective stress strengths and undrained strengths. Partial factors are applied to the effective stress
strength properties and undrained shear strengths in all stages of the analysis and in a manner
consistent with material models in categories 1 and 2.

134
14.4 Reinforcement Loads
Partial factors are applied to the reinforcement resistance parameters: pullout resistance, shear force,
and tensile capacity. The partial factor is applied to the pullout resistance even if it is calculated from
overburden stress and interface strength properties.

15 References
Abramson, L.W., Lee, T.S., Sharma, S., and Boyce, G.M., 2002. Slope Stability and Stabilization Methods.
John Wiley & Sons Inc. pp.712.

Budhu, M. 2007. Soil Mechanics and Foundations, 2nd Edition. John Wiley & Sons.

Bishop, A.W. and Morgenstern, N., 1960. Stability coefficients for earth slopes. Géotechnique, Vol. 10,
No. 4, pp. 164 169.

Christian, J.T., Ladd, C.C. and Baecher, G.B., 1994. Reliability Applied to Slope Stability Analysis. Journal
of Geotechnical Engineering, Vol. 120, No. 12. pp. 2180-2207.

Corps of Engineers (2003) “Appendix G - Procedures and Examples for Rapid Drawdown”. Engineering
Manual, EM 1110-2-1902. Department of the U.S Army Corps of Engineers.

Duncan, J.M.,Wright S.G. and Wong, K.S. (1990). “Slope Stability during Rapid Drawdown”. Proceedings
of H. Bolton Seed Memorial Symposium. Vol. 2.

El-Ramly, H., Morgenstern, N.R., Cruden, D.M. 2002. Probabilistic Slope Stability Analysis for Practice.
Canadian Geotechnical Journal. Vol. 39. pp. 665-683.

Fellenius, W., 1936. Calculation of the Stability of Earth Dams. Proceedings of the Second Congress on
Large Dams, Vol. 4, pp. 445-463.

Fredlund, D.G., 1974. Slope Stability Analysis. User's Manual CD-4, Department of Civil Engineering,
University of Saskatchewan, Saskatoon, Canada.

Fredlund, D.G., and Krahn, J., 1977. Comparison of slope stability methods of analysis. Canadian
Geotechnical Journal, Vol. 14, No. 3, pp. 429 439.

Fredlund, D.G., Krahn, J., and Pufahl, D. 1981. The relationship between limit equilibrium slope stability
methods. Proc. 10th Int. Conf. Soil Mech. Found. Eng. (Stockholm, Sweden), Vol. 3, pp. 409-416.

Fredlund, D.G., Zhang, Z.M. and Lam, L., 1992. Effect of the Axis of Moment Equlibrium in Slope Stability
Analysis. Canadian Geotechnical Journal, Vol. 29, No. 3.

Fredlund, D.G., Xing, A. Fredlund M.D., and Barbour, S.L., 1996. The Relationship of the Unsaturated Soil
Shear Strength to the Soil-water Characteristic Curve. Canadian Geotechnical Journal, Vol. 33, pp. 440-
448.

135
Greco, V.R. 1996. Efficient Monte Carlo Technique for Locating Critical Slip Surface. Journal of
Geotechincal Engineering. Vol 122, No. 7. ASCE, pp.517-525.

Grivas, D.A., 1981. How Reliable are the Present Slope Failure Prediction Methods? Proceedings of the
Tenth International Conference of Soil Mechanics and Foundation Engineering, Stockholm, Sweden, Vol.
3, pp.427-430.

Harr, M.E., 1987. Reliability-Based Design in Civil Engineering. McGraw-Hill Book Company. pp. 290.

Higdon A., Ohlsen, E.H., Stiles, W.B., Weese, J.A. and Riley, W.F., 1978. Mechanics of Materials. John
Wiley & Sons. pp.752.

Hoek, E. 2000. Course Notes entitled Practical Rock Engineering (2000 ed., Chapter 11).

Hoek, E., Bray, J.W., 1974. Rock Slope Engineering – Appendix 3. Published for the Institute of Mining
and Metallurgy., pp. 352 – 354.

Hoek, E., Carranza-Torres, C. and Corkum, B., 2002. Hoek-Brown Failure Criterion – 2002 Edition. Proc.
North American Rock Mechanics Society meeting in Toronto in July 2002.

Janbu, N., Bjerrum, J. and Kjaernsli, B., 1956. Stabilitetsberegning for Fyllinger Skjaeringer og Naturlige
Skraninger. Norwegian Geotechnical Publications, No. 16, Oslo.

Janbu, N. 1954. Applications of Composite Slip Surfaces for Stability Analysis. In Proceedings of the
European Conference on the Stability of Earth Slopes, Stockholm, Vol. 3, p. 39-43.

Krahn, J., 2003. The 2001 R.M. Hardy Lecture: The Limits of Limit Equilibrium Analyses. Canadian
Geotechnical Journal, Vol. 40. pp.643-660.

Krahn, J., Price, V.E., and Morgenstern, N.R., 1971. Slope Stability Computer Program for Morgenstern-
Price Method of Analysis. User's Manual No. 14, University of Alberta, Edmonton, Canada.

Kramer, S.L. 1996. Geotechnical Earthquake Engineering. Prentice Hall, pp. 437.

Ladd, C.C., 1991. Stability evaluation during staged construction: 22nd Terzaghi Lecture. Journal of
Geotechnical Engineering, ASCE, 117(4), pp. 537-615

Ladd, C.C. and Foott, R. 1974. New Design Procedure for Stability of Soft Clays. Journal of the
Geotechnical Engineering Division, ASCE, Vol. 100, (GT7), pp. 763-786.

Lam, L., and Fredlund D.G., 1993. A general Limit Equilibrium Model for Three-Dimensional Slope
Stability Analysis. Canadian Geotechnical Journal. Vol. 30. pp. 905-919.

Lambe, T.W. and Whitman, R.V., 1969. Soil Mechanics. John Wiley and Sons, pp. 359-365.

Lapin, L.L., 1983. Probability and Statistics for Modern Engineering. PWS Publishers. pp. 624.

136
Li, K.S. and Lumb, P., 1987. Probabilistic Design of Slopes. Canadian Geotechnical Journal, Vol. 24, No. 4,
pp. 520-535.

Lumb, P., 1966. The Variability of Natural Soils. Canadian Geotechnical Journal, Vol. 3, No. 2, pp. 74-97.

Lumb, P., 1970. Safety Factors and the Probability Distribution of Soil Strength . Canadian Geotechnical
Journal, Vol. 7, No. 3, pp. 225-242.

Malkawi, A.I.H, Hassan, W.F. and Sarma, S.K. 2001. Global Serach Method for Locating General Slip
Surface Using Monte Carlo Technique. Journal of Geotechnical and Geoenvironmental Engineering. Vol
127, No. 8, pp. 688-698.

Morgenstern, N.R., and Price, V.E., 1965. The Analysis of the Stability of General Slip Surfaces.
Géotechnique, Vol. 15 (1), pp. 79-93.

Mostyn, G.R. and Li, K.S., 1993. Probabilistic Slope Stability Analysis - State-of-Play, Proceedings of the
Conference on Probabilistic Methods in Geotechnical Engineering, Canberra, Australia. pp. 281-290.

Petterson, K.E. 1955. The Early History of Circular Sliding Surfaces, Géotechnique, Vol. 5, p. 275-296.

Sarma, S.K., 1973. Stability Analysis of Embankments and Slopes. Géotechnique, Vol. 23 (3), pp. 423-433.

Sarma, S.K., 1979. Stability Analysis of Embankments and Slopes. Journal of the Geotechnical
Engineering Division. ASCE. Vol. 105, No. GT12. pp. 1511-1524.

Smith, G.N., 1974. Elements of Soil Mechanics for Civil and Mining Engineers. Crosby Lockwood Staples,
London. pp. 126-144.

Spencer, E. 1967. A Method of Analysis of Embankments assuming Parallel Inter-Slice Forces.

Géotechnique, Vol 17 (1), pp. 11-26.

Soulié, M., Montes, P. and Silvestri, V. 1990. Modelling Spatial Variability of Soil Parameters, Canadian
Geotechnical Journal, Vol. 27, No. 5, pp. 617-630.

van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulic conductivity of
unsaturated soils. Soil Science Society of America Journal 44(5): 892-898.

Vanapalli, S.K., Fredlund D.G., Pufahl, D.E. and Clifton, A.W., 1996. Model for the Prediction of Shear
Strength with respect to Soil Suction. Canadian Geotechnical Journal, Vol. 33, pp. 379-392.

Vanmarcke, E.H., 1983. Random Fields: Analysis and Synthesis. MIT Press, Cambridge, Mass.

Whitman, R.V. and Bailey, W.A., 1967. Use of Computer for Slope Stability Analysis. Journal of the Soil
Mechanics and Foundation Division of ASCE, Vol. 93, No. SM4.

Wolff, T.F, 1985. Analysis and Design of Embankment Dams: A Probabilistic Approach. Ph.D. Thesis,
Purdue University, West Lafayette, IN.

137
Yang X-S and Deb S. Engineering optimisation by cuckoo search. PhD thesis. International Journal of
Mathematical Modelling and Numerical Optimisation. 2010; 1:330–343.

138