Numerical Modeling and Simulation of the Stability of Earth Slopes 
 
A Thesis Submitted to the Department of Nuclear Engineering 
 
DEPARTMENT OF NUCLEAR ENGINEERING 
SCHOOL OF NUCLEAR AND ALLIED SCIENCES 
COLLEGE OF BASIC AND APPLIED SCIENCES 
UNIVERSITY OF GHANA 
 
BY 
 
BRENDAN DAGEMANYIMA ATARIGIYA, 10507155 
BSc (KNUST, Kumasi), 2012 
 
 
In Partial Fulfilment of the Requirements for the Degree of 
 
 
MASTER OF PHILOSOPY 
IN 
COMPUTATIONAL NUCLEAR SCIENCE AND ENGINEERING 
 
July 2016 
1 
 
DECLARATION 
I hereby declare that with the exception of references to other people’s work which 
have duly been acknowledged, this Thesis is the result of my own research work and 
no part of it has been presented for another degree in this University and elsewhere. 
 
 
……………………………    ………………………….. 
Brendan Dagemanyima Atarigiya      Date 
(Candidate) 
 
I hereby declare that the preparation of this project was supervised in accordance with 
the guidelines of the supervision of Thesis work laid down by the University of Ghana 
 
 
…………………………    ………………………… 
Dr. Nii Kwashie Allotey   Nana (Prof.) A. Ayensu Gyeabour I 
(Principal Supervisor)      (Co-Supervisor) 
 
…………………………..   ………………………………….. 
Date       Date 
  
ii 
 
DEDICATION 
I dedicate this work to God Almighty, my family and friends. 
  
iii 
 
ACKNOWLEDGMENT 
Firstly I want to thank and praise God for the good health and strength given me during 
all this period of schooling.  
My sincere and utmost gratitude goes to my Principal Supervisor, Dr. Nii Kwashie 
Allotey of the Ghana Atomic Energy Commission (GAEC) and the University of Ghana 
for his expertise in the field of research and exemplary guidance towards the progress 
of this research. I am also grateful to Prof. (Nana) A. Ayensu Gyeabour my Co-
Supervisor for his creative suggestions and motivation through this research project. 
My thanks also goes to Ms. Rita Awura Abena Appiah of GAEC, thank you for pushing 
me to make this possible. You helped me improve my programming skills a lot.  
To my big family, the Atarigiya, Aguyire, and Allotey families, thank you for the 
support. I could not have done it without you.  
To my amazing course mate, Linda Sarpong, thank you for your words of 
encouragement during out period in school; you are a strong woman!  Furthermore, to 
the nuclear engineering department (Samiru, Efia, Matilda, Henry…), thanks for the 
wonderful time we had together. 
Last and not the least, to the woman who stood firmly behind me from day one, from 
when this journey began, Mercy Selina Somhayin Namateng, I LOVE YOU. 
GOD RICHLY BLESS YOU ALL. 
 
  
iv 
 
ABSTRACT 
Ghana, as most other countries, has a considerable variation in its topography.  In an 
attempt to build cheaper, but yet the safe structures (i.e., roads, apartments, etc.), we 
are most often times faced with building on hill-sides and in valleys.  This then calls 
for the need to correctly assess the stability of any adjacent slopes.   
In recent times, due to the extensive need for stability analysis in engineering practice, 
slope stability analysis programs have been developed.  It is noted that these 
commercial slope stability programs are used extensively in the industry but are very 
expensive and require purchasing yearly licenses.  As a result of this, slope stability 
analysis is not routinely conducted in local geotechnical engineering practice.  The need 
for cheaper more accessible options is thus considered needful.   
This research initiative uses MATLAB, a commercially available, user-friendly and 
easy to access computing platform to develop a slope stability analysis program.  The 
method used is the General Limit Equilibrium Method (GLE) with the adoption of the 
Morgenstern-Price (M-P) factor of safety (FoS) approach to develop a cheap, efficient, 
and yet effective model for slope stability analysis and design.  The results of the 
program are validated by comparing with the results of SLOPE/W, a commercial slope 
stability program.   
The results show four model outputs from the developed program and SLOPE/W for a 
homogeneous material. Two different failure mechanisms are shown (i.e., toe and base 
failures).  It is noted that the percentage error in the M-P FoS is less than 5%.  
It is anticipated that with the availability of this computer code, Ghanaian Engineers 
can more readily assess the safety of slopes in routine design works. 
  
v 
 
Contents 
DECLARATION ........................................................................................................... ii 
DEDICATION .............................................................................................................. iii 
ACKNOWLEDGMENT............................................................................................... iv 
ABSTRACT ................................................................................................................... v 
LIST OF FIGURES ...................................................................................................... ix 
LIST OF TABLES ........................................................................................................ xi 
CHAPTER ONE: INTRODUCTION ............................................................................ 1 
1.1 Background ..................................................................................................... 1 
1.2 Problem Statement .......................................................................................... 3 
1.3 Relevance and Justification of Study .............................................................. 4 
1.4 Research Goal ................................................................................................. 5 
1.5 Research Objectives ........................................................................................ 5 
1.6 Scope ............................................................................................................... 5 
1.7 Format of the Thesis ........................................................................................ 6 
CHAPTER TWO: LITERATURE REVIEW ................................................................ 7 
2.1 Factors Causing Instability .............................................................................. 7 
2.2 Types of Slip Surfaces..................................................................................... 8 
2.3 Definition of FoS ........................................................................................... 10 
2.4 Slope Stability Analysis Methods ................................................................. 12 
2.4.1 Limit Analysis Method .......................................................................... 12 
2.4.2 Variational Calculus Method ................................................................. 13 
2.4.3 Strength Reduction Method ................................................................... 14 
2.4.4 General Discussion on Limit Equilibrium Method ................................ 15 
2.5 General Limit Equilibrium Method of Slices (GLE method) ....................... 19 
vi 
 
2.6 The Ordinary or Fellenius Method (OMS) ................................................... 20 
2.7 Simplified Bishop’s Method ......................................................................... 21 
2.8 Janbu’s Simplified Method ........................................................................... 22 
2.9 Morgenstern-Price (M-P) Method ................................................................. 23 
CHAPTER THREE: RESEARCH METHODOLOGY AND SOFTWARE USED ... 27 
3.1 Selection of Factor of safety method............................................................. 27 
3.2 Morgenstern-Price Method ........................................................................... 28 
3.3 Assumptions .................................................................................................. 29 
3.4 Numerical Method Development .................................................................. 29 
3.5 Derivation of Equations ................................................................................ 30 
3.6 Structured Program ....................................................................................... 32 
3.7 Numerical Algorithm .................................................................................... 33 
3.8 Software and Programs Used ........................................................................ 34 
3.8.1  SLOPE/W ............................................................................................. 34 
3.8.2 MATLAB ............................................................................................... 36 
CHAPTER FOUR: RESULTS AND DISCUSSION .................................................. 37 
4.1 Introduction ................................................................................................... 37 
4.2 Programme Test Examples............................................................................ 37 
4.2.1 Toe Failure ............................................................................................. 37 
4.2.2 Base Failure ........................................................................................... 42 
4.3 Comparison with Cases from Literature ....................................................... 46 
4.3.1 Case 1 ..................................................................................................... 46 
4.3.2 Case 2 ..................................................................................................... 50 
CHAPTER FIVE: CONCLUSIONS ........................................................................... 53 
5.1 Conclusions ................................................................................................... 53 
5.2 Recommendations ......................................................................................... 54 
REFERENCES ............................................................................................................ 55 
vii 
 
APPENDICES ............................................................................................................. 60 
Appendix A: Structured Programme ........................................................................ 60 
Appendix B: MATLAB Code for Solving FoS ....................................................... 63 
 
  
viii 
 
LIST OF FIGURES 
Figure 1.1: Over 40 m high slope on the Ayi-Mensah-Aburi road   3 
Figure 2.1: Types of circular slip failure surface     9 
Figure 2.2: Typical non-circular slip surfaces     10 
Figure 2.3: Various definitions for FoS      11 
Figure 2.4: Swedish Slip Circle Method      16 
Figure 2.5: Slice Discretization and Slice Forces in a Sliding Mass   17 
Figure 2.6: Forces considered in the Ordinary Method of Slices   21 
Figure 2.6: Forces considered in the Ordinary Method of Slices   22 
Figure 2.8: Forces considered in the M-P method     24 
Figure 2.9: Inter-slice force function types      25 
Figure 2.10: Variation of FoS with respect to Fm and Ff vs. λ for the M-P method 26 
Figure 3.1: Sketch of a Slope Section       30 
Figure 3.2: Forces acting on a Single Slice from a Mass Slope   31 
Figure 3.3: Algorithm flowchart for solving for the FoS    33 
Figure 3.4: SLOPE/W KeyIn Analyses Page      35 
Figure 3.5: SLOPE/W KeyIn Material Page      35 
Figure 3.6: SLOPE/W KeyIn Entry and Exit Range Page    36 
Figure 4.1: SLOPE/W Output of Toe Failure: Case 1    38 
Figure 4.2: MATLAB Output of Toe Failure: Case 1    39 
Figure 4.3: SLOPE/W Output for Toe Failure: Case 2    40 
ix 
 
Figure 4.4: MATLAB Output for Toe Failure: Case 2    41 
Figure 4.5: SLOPE/W Output for Base Failure: Case     42 
Figure 4.6: MATLAB Output for Base Failure: Case 1    43 
Figure 4.7: SLOPE/W Output for Base Failure: Case 2    44 
Figure 4.8: MATLAB Output for Base Failure: Case 2    45 
Figure 4.9: Homogeneous Slope without Foundation    46 
Figure 4.10: Analysis using SLOPE/W - FoS = 1.385    47 
Figure.4.11: FoS based on M-P approach for Toe failure - FoS = 1.451  48 
Figure 4.12: FoS based on M-P approach for Base Failure: FoS = 1.375  49 
Figure 4.13:  Slope Model Geometry from Slide 3.     50 
Figure 4.14: FoS based on M-P Approach for ACAD Problem   51 
 
  
x 
 
LIST OF TABLES 
Table 2.1: The Main Limit Equilibrium Methods      18 
Table 3.1: Brief Comparison of Limit Equilibrium Methods    27 
Table 4.1: Slope Dimensions and Material Properties for Toe Failure: Case 1 38 
Table 4.2: Slope Dimensions and Material Properties for Toe Failure: Case 2 40 
Table 4.3: Slope Dimensions and Material Properties for Base Failure: Case 1 42 
Table 4.4: Slope Dimensions and Material Properties for Base Failure: Case 2 44 
Table 4.5: Slope dimensions and material properties     46 
Table 4.6: Slope Dimensions and Material Properties for ACAD Problem  50 
Table 5.1: Summary of FoS Outputs for all Case Studies.    53 
 
  
xi 
 
LIST OF SYMBOLS AND ABBREVIATIONS  
A  cross-sectional area of slice 
dx  width of slice  
l  length of the bottom of the slice 
c  cohesion of soil 
G  total unit weight of soil 
τf shear strength 
τ shear stress 
γw unit weight of water 
 angle of internal friction of soil 
α  inclination from horizontal of the bottom of the slice (degrees) 
cal cos(α) 
sal sin(α) 
tph tan() 
H Height of slope 
N  total normal force on the bottom of the slice 
S  shear force on the bottom of the slice 
W  weight of the slice 
havg = average height of slice  
u  pore water pressure 
xii 
 
∆X shear component of the inter-slice force 
∆E  inter-slice force on the downslope of the slice 
ns number of slice faces 
FoS factor of safety 
F assumed factor of safety 
Ff force factor of safety 
Fm  moment factor of safety 
Xc x coordinate of centre of slip circle 
Yc y coordinate of centre of slip circle 
R Radius of slip circle 
ytop Ground surface  
ybot Slope surface 
FSom Factor of Safety for Ordinary Method of Slices 
FFm Moment Factor of Safety for Morgentern-Price Method of Slices 
FFf Force Factor of Safety for Morgentern-Price Method of Slices 
F(x) Interslice force function 
λ Scale factor of the assumed f(x) 
 
xiii 
 
CHAPTER ONE: INTRODUCTION 
1.1 Background  
Slope stability is the potential for ground slopes to resist movement [1]. Slope 
instability has been the subject of continued concern because of the tremendous loss of 
life, property and infrastructure caused annually in many places in the world [2].  In the 
field of construction, slope instability can occur due to rainfall, increasing the water 
table, and the change in stress conditions. Similarly, tracks of land that have been stable 
for years may suddenly fail due to changes in the geometry, external forces and loss of 
shear strength [3]. 
Slope failures, also called slides or landslides, whether sudden or gradual, are due to 
the increased stress of slope materials or foundations compared with their mobilized 
strength [3]. 
The majority of the slope stability analyses performed in practice still use traditional 
limit equilibrium approaches involving the method of slices, and has remained virtually 
unchanged for decades. 
Analysis of the slope stability is carried out to assess the safety of artificial or natural 
slopes (e.g., dams, road cuts, mining open pit excavations and landfills). For human 
made slopes, analysis of slope stability is used to evaluate various design options, which 
then provides a basis for a form of engineering design with associated costing 
comparisons.  The efficient engineering of natural and artificial slopes has therefore 
become a common challenge faced by both researchers and practitioners. 
Slope stability assessment mainly involves the use of the factor of safety (FoS) method 
to determine how close a given slope is to the onset of instability, or to what extent the 
state of the slope is from failure. 
1 
 
When this ratio is well above 1, the resistance to shear failure is generally higher than 
the driving shear stress, and the slope is considered stable.  When this ratio is near to 1, 
the shear strength is almost equal to the shear stress, and the slope is close to failure.  If 
the FoS is less than 1, the slope is considered to have failed, or considered to be trigger-
point ready [4]. 
Ghana is not noted to be a frequent serious victim of mass movement (slope failures).  
It is, however, noteworthy that Ghana has not been without slope failures.  In Ref. [5], 
reference is made of a slope failure which occurred in 1968 near Jamasi in the Sekyere 
District involving about 1,500 cubic meters of rock, soil and vegetation. The failure 
blocked the main Kumasi-Mampong truck road for a total of ten days. 
Reference is also made in Refs. [6 & 7] on potential slope failures on the stretch of the 
Accra-Aburi road, when rocks began to fall unhindered onto the road in 2014.  
Furthermore, Ref. [8] notes that in 2013 after a heavy downpour of rain, loose parts of 
the Akoaso Mountain fell and blocked the Kasoa-Weija portions of the road, resulting 
in a significant traffic jam for hours.  
These recent records of slope instability in the country have served as a wakeup call for 
researchers and practicing engineers to take a critical look at this issue.   
2 
 
 
Figure 1.1: Over 40 m high slope on the Ayi-Mensah-Aburi road 
 
1.2 Problem Statement 
In an attempt to build cheap and yet the safe structures (i.e., roads, living apartments, 
etc.) for man-kind, we are most times faced with building in valleys and on mountains.  
Either way, we are faced with the problem of slope instability.   
In the past decades, computer software for slope stability analysis and design have been 
developed and marketed extensively.  These commercial software, which have been 
developed over many years, are able to perform rigorous stability calculations, and give 
fast and accurate answers to complex slope stability problems.  These software have 
become widely accepted in industry, and are now part of most large design engineering 
offices.  These software are, however, expensive and normally require the annual 
renewal of licenses.  Notwithstanding, their wide acceptance in industry, most 
Ghanaian geotechnical engineers still resort to the use of the rule of thumb slope 
stability analysis methods, and old charts for their daily slope stability analysis. This is 
3 
 
due to the relatively high cost of these commercial software, and the limited financial 
capacity of local engineering firms. 
This has created a gap between local and international engineers, and has resulted in 
major geotechnical engineering projects involving complex slope instability problems 
being awarded to international firms, rather than local engineering firms. 
 
1.3 Relevance and Justification of Study  
Landslides, rock falls, and mass movement of any kind, are undoubtedly, one of the 
oldest natural disasters that have resulted in huge damages, loss of lives, and a great 
deal of pain to mankind.  Like other mountainous countries, Ghana has large variations 
in its topography.  The impending threat of landslides in the case of [5 or 6], or rock 
falls in the case of [7] is now accepted as life threating, and the need for these slopes to 
be properly engineered is critical.   
It has been already noted that the available commercial slope stability programs that 
are extensive slope stability analysis are very expensive and require purchasing yearly 
licenses.  This has necessitated the need to develop a simple, yet efficient slope stability 
program, that can be easily accessed by local engineers to appraise local slope stability 
problems using the most rigorous and accurate methods available. 
This research initiative uses MATLAB, (a commercially available, and  easily 
accessible computing platform with great user-friendly interface) using General Limit 
Equilibrium method (GLE) and adopting the Morgenstern-Price FoS to develop a cheap 
and efficient, yet effective model for slope stability analysis and design.  
 
4 
 
1.4 Research Goal  
The goal of this thesis is to develop and implement a slope stability numerical model 
that will aid local geotechnical engineers to readily appraise the stability of slopes in 
local practice. 
 
1.5 Research Objectives 
The objectives of the research are to: 
 Develop a physical model to represent the problem 
 Develop mathematical equations to solve the problem 
 Develop a numerical algorithm and write a code to solve for the FoS of earth 
slopes. 
 Verify and validate the code using Geoslope International’s SLOPE/W 
commercial slope stability programme. 
 
1.6 Scope 
For the purpose of this study, this thesis is limited to the development of a slope stability 
programme for homogenous soil and rock media.  In this regard, the goal of the study 
is to develop the generic algorithm for slope stability analysis. Furthermore, similar to 
the existing commercial programmes, the study is limited to two-dimensional slope 
stability problems.  
 
 
5 
 
1.7 Format of the Thesis  
Following the introduction to slope stability problems in Chapter 1, a detailed literature 
review of the methods of slope stability analysis is provided in Chapter 2.  Chapter 3 
then presents the proposed solution method for calculating the factor of safety, in which 
the GLE method is explained, and the solution algorithm developed.   
Chapter 4 presents sample results from the developed MATLAB code.  It presents 
comparisons between the results of the developed code and results from the commercial 
slope stability programme, SLOPE/W.  Chapter 5 finally presents the conclusions of 
the study, and also provides recommendations for further studies.   
  
6 
 
CHAPTER TWO: LITERATURE REVIEW 
Analysis of slope stability problem is an important subject area in geotechnical 
engineering. The phenomenon of landslides and related slope instability is a problem 
in many parts of the world.  Slope failure mechanisms and the geological history of a 
slope can be very complicated problems and required complex forms of analyses. 
 
2.1 Factors Causing Instability 
The failure of the slopes occurs when the downward movements of soil or rock material 
because of gravity and other factors, creates shear stresses that exceed the inherent shear 
strength of the material. Therefore, factors that tend to increase the shear stress or 
decrease the shear strength of a material increases the risk of the failure of a slope. 
Various processes can lead to a reduction of the shear strength of a soil/rock mass.  
These include: increased pore pressure, cracking, swelling, decomposition of 
argillaceous rock fills, creep under sustained loads, leaching, softening, weather and 
cyclic loading, among others.   
On the other hand, shear stress within a rock/soil mass may increase due to additional 
loads on top of the slope, and increase in water pressure due to cracks at the top of the 
slope, an increase in the weight of soil due to increasing water content, the excavation 
of the base of the slope, and seismic effects.  Furthermore, additional factors that 
contribute to the failure of a slope include the rock/soil mass properties, slope geometry, 
state of stress, temperature, erosion, etc.   
The presence of water is the most critical factor that affects the stability of slopes.  This 
is because it increases both the driving shear stress, and also decreases the soil/rock 
mass’ shear strength.  The speed of sliding movement in a slope failure can vary from 
7 
 
a few millimeters per hour to very fast slides where great changes have occurred within 
seconds.  Slow slides occur in soils with a plastic stress-strain characteristics, where 
there is no loss of strength with increased strain.  Fast slides occurs in situations where 
there is a sudden loss of strength, as in the liquefaction of sensitive clay or fine sand. 
Increase in shear stresses across the soil mass result in movement only when the shear 
strength mobilized on given possible failure surface in the ground is less than the 
driving shear stresses along that surface.   
 
2.2 Types of Slip Surfaces 
To calculate the FoS of a slope, it is always assumed that the slope is failing in some 
shape, normally in a circular or non-circular shape. For computational simplicity, the 
slide surface is often seen as circular or composed of several straight lines [9]. Different 
sliding surfaces are normally assumed with the computation of a corresponding FoS.  
The sliding surface with the minimum FoS is then selected as the FoS of the slope in 
question. 
A circular sliding surface, like that shown in Figure 2.1, is often used because it is 
suitable to sum the moments about a centre. The use of a circle also simplifies the 
calculations. Wedge-like surfaces have their failure mechanisms defined by three or 
more straight line segments defining an active area, central block, and the passive area 
as shown in Figure 2.2. This type of sliding surface can be used for analysis of slopes 
where the critical potential sliding surface comprises a relatively long linear sector 
through low material bounded by a stronger material. 
As noted above, the critical slip surface is the surface with the lowest factor of safety.  
The critical slip surface for a given problem analysed by a given method, is found by a 
8 
 
systematic procedure to generate sliding test surfaces until the one with the minimum 
safety factor is found [9]. 
 
 
Figure 2.1: Types of circular slip failure surface [3] 
  
9 
 
 
Figure 2.2: Typical non-circular slip surfaces [9] 
 
2.3 Definition of FoS 
FoS is usually defined as the ratio of the ultimate shear strength to the shear stress 
mobilized at imminent failure. There are several ways to formulate the FoS, The most 
common formulation assumes the safety factor to be constant along the sliding surface, 
and it is defined in relation to limit equilibrium, force and moment equilibrium [4]. 
10 
 
These definitions are given in Figure 2.3 below. As will be developed further in the 
limit equilibrium method, the first definition is based on the shear strength which can 
be obtained in two ways: a total stress approach (su‐analysis) and an effective stress 
approach (c’- φ’ −analysis).  The type of strength consideration depends on the soil 
type, the loading conditions and the time elapsed after excavation.  The total stress 
strength method is used for short–term conditions in cohesive soils, whereas, the 
effective stress method is used in long- term conditions in all soil types, or in short-term 
conditions in cohesive soils where the pore pressure is known [3].  
 
Figure 2.3: Various definitions for FoS [3] 
11 
 
2.4 Slope Stability Analysis Methods 
Slope stability problems deal with the condition of ultimate failure of a soil or rock 
mass. Analyses of slope stability, bearing capacity and earth pressure problems, all fall 
into this area. The stability of a slope can be analysed by a number of methods, among 
others of which are the:  
 Limit analysis method,  
 Variational calculus method,  
 Strength reduction method and 
 Limit equilibrium method,  
 
2.4.1 Limit Analysis Method 
The limit analysis method theory is based on a rigid-perfectly plastic model material.  
Drucker and Prager [10] first formulated and introduced the upper and lower bound 
plasticity theorems for soil/rock masses. 
The general analysis process includes construction of a statically admissible stress field 
for the lower-bound analysis, or a kinematically admissible velocity field for the upper-
bound analysis. 
For both upper- and lower-bound analysis, one of the following two conditions has to 
be satisfied: 
 Geometrical compatibility between internal and external displacements or 
strains.  This is usually concerned with kinetic conditions – velocities must be 
compatible to ensure no gain or loss of material at any point. 
 Stress equilibrium, i.e., the internal stress fields must balance the externally 
applied stresses (forces). 
12 
 
The basis of limit analysis rests upon two theorems, which can be proved 
mathematically. In simple terms, these theorems are: 
 Lower Bound: any stress system in which the applied forces are just sufficient 
to cause yielding. 
 Upper Bound: Any velocity field that can operate is associated with an upper 
bound solution. 
The lower- bound approach has been used in 2D slope stability analysis by Zhang [11] 
and Kim et al. [12].  The upper-bound approach was first used in 2D slope stability 
analysis by Drucker and Prager [10] to determine the critical height of a slope.  
Subsequently, Refs. [13, 14 & 15] also applied and extended the upper-bound 
approaches in 2D slope analysis. Michalowski [16] proposed an upper-bound approach 
based on a translational failure mechanism.  The vertical slice techniques, which are 
often used in traditional limit equilibrium approaches, were employed to satisfy the 
force equilibrium condition for all individual slices. Two extreme kinematical solutions 
neglecting the inter-slice strength, or fully utilizing the inter-slice strength of the soil 
were then obtained. The traditional limit equilibrium solutions for slices with a proper 
implicit assumption of failure mechanism can fall into the range of these two extremes.  
Donald and Chen [17] also presented an upper-bound method on the basis of a multi-
wedge failure mechanism, in which the sliding body was divided into a small number 
of discrete blocks. 
 
2.4.2 Variational Calculus Method 
The variational calculus approach does not require assumptions on the inter-slice 
forces.  It was first used for 2D stability analysis by Baker and Garber [18].  This 
13 
 
approach was subsequently employed by Jong [19] for vertical-cut analysis in cohesive 
directionless soils.  Cheng, et al. [20] also used it in their research where they developed 
a numerical algorithm based on the extremum principle by Pan [21].  The formulation 
which relies on the use of a modern try and error optimization method and can be 
viewed as an equivalent form of the variational method in a discretized form, but is 
applicable for a complicated real life problem.  
 
2.4.3 Strength Reduction Method 
In recent decades, there have been great developments in the area of the strength 
reduction method (SRM) for slope stability analysis.  The general procedure of the 
SRM analysis is the reduction of the strength parameters by the FoS, while the body 
forces (due to the weight of soil and other external loads) are applied until the system 
cannot maintain a stable state. This procedure can determine the FoS within a single 
framework, for both two and three-dimensional slopes. 
The main advantages of the SRM are as follows; 
 The critical failure surface is automatically determined from the 
application of gravity loads and/or the reduction of shear strength, 
 It requires no assumption about the distribution of the inter-slice shear 
forces, 
 It is applicable to many complex conditions, and,  
 It can give information such as stresses, movements (deformations) and 
pore pressures. 
One of the main disadvantages of the SRM is the long time required to develop the 
computer model and to perform the analysis to arrive at a solution. With the 
14 
 
development of computer hardware and software, 2D -SRM can now be done in a 
reasonable amount of time suitable for routine analysis and design. This technique is 
also adopted in several well-known commercially available geotechnical finite element 
or finite difference programs. In strength reduction analysis, the convergence criterion 
is the most critical factor in the assessment of the FoS. 
Investigation results show that; the FoS obtained and the corresponding slip surface 
determined by the SRM, demonstrate good agreement with the results of the Limit 
Equilibrium Method (LEM).  
 
2.4.4 General Discussion on Limit Equilibrium Method 
2.4.4.1  Swedish Slip Circle Method of Analysis 
The Swedish Slip Circle method assumes that the friction angle of the soil or rock is 
zero.  In other words, when angle of friction is considered to be zero, the effective stress 
term tends to zero, which is therefore equivalent to the shear cohesion parameter of the 
given soil.  The Swedish slip circle method assumes a circular failure interface, and 
analyses stress and strength parameters using circular geometry and statics as shown in 
Figure 2.4. The moment caused by the internal driving forces of a slope is compared 
with the moment caused by resisting forces in the sliding mass. If forces resisting 
movement are greater than the forces tending to cause movement, then the slope is 
assumed to be stable. 
15 
 
 
Figure 2.4: Swedish Slip Circle Method [3] 
 
2.4.4.2  Method of Slices 
Despite all the above methods, limit equilibrium methods are by far the most used form 
of analysis for slope stability studies. They are the oldest best-known numerical 
technique in geotechnical engineering.  These methods involve cutting the slope into 
fine slices so that their base can be comparable with a straight line.  The governing 
equilibrium equations equilibrium of the forces and/or moments, Figure 2.5 are then 
developed. According to the assumptions made on the efforts between the slices and 
the equilibrium equations considered, many alternatives have been proposed in Table 
2.1. They give, in most cases, quite similar results. The differences between the values 
of the FoS obtained with the various methods are generally below 6% [22]. 
 
 
16 
 
 
 
 
 
 
 
 
 
 
 
 
Figure 2.5: Slice Discretization and Slice Forces in a Sliding Mass [23.24] 
The idea of dividing a potential sliding mass into slices dates back to the early 1900’s.  
The first documented use of the method of slices is the analysis of the 1916 failure at 
the Stigberg Quay in Gothenburg, Sweden [25]. 
This Limit Equilibrium method is well known to be a statically indeterminate problem 
and assumptions about the inter-slice shear forces are needed to make the problem 
statically determinate.  On the basis of the assumptions on the internal forces and the 
force and/or moment equilibrium, there are more than ten methods developed for 
analysis of slope stability problems [26]. Famous methods include those by Fellenius 
[27], Bishop [28], Janbu [29, 30], Spencer [31], Morgenstern-Price [32] and the General 
Limit Equilibrium (GLE) [25]. 
Table 2.1 shows the differences between the various methods of stability analysis, on 
the basis of forces and moments equilibrium.   
17 
 
Table 2.1: The Main Limit Equilibrium Methods [5] 
Inter Inter 
Moment Force 
Moment Force Slice Slice Inter Slice 
No  Methods Factor of Factor of 
Equilibrium Equilibrium Normal Shear Force Function 
Safety Safety 
Forces Forces 
Fellenius, Swedish Circle or 
1  Yes No No No Yes No No 
Ordinary Method (1936) 
2  Bishop Simplified (1955)  Yes No Yes No Yes No No 
3  Janbu Simplified (1954)  No Yes Yes No No Yes No 
4  Spencer Method (1967)  Yes Yes Yes Yes Yes Yes Constant 
Constant Half- 
Sine 
Morgenstern-Price Method 
5  Yes Yes Yes Yes Yes Yes Clipped-Sine 
(1965)  
Trapezoidal 
Specified 
 
18 
 
All limit equilibrium methods utilise the Mohr‐Coulomb criterion to determine the 
shear strength (τf) along the sliding surface.  The mobilised shear stress at which a soil 
fails in shear is defined as the shear strength of the soil. According to Janbu [29], a state 
of limit equilibrium exists when the mobilised shear stress (τ) is expressed as a fraction 
of the shear strength. Nash [33] states that, at the moment of failure, the shear strength 
is fully mobilised along the failure surface when the critical state conditions are 
reached. The shear strength is usually expressed by the Mohr‐ Coulomb linear 
relationship as:  
𝑓 = 𝑐 + tan(⁡ )        (2.1) 
and  
⁡ = ⁡ 𝜏𝑓/𝐹𝑜𝑆         (2.2) 
where  
c is the soil cohesion and  is the soil frictional angle. 
The strength of the soil available depends partially on the type of soil and the normal 
stresses acting on it. The mobilized shear strength on the other hand depends on the 
external forces acting on the soil mass. This defines the FoS as a ratio of the τf to τ in a 
limit equilibrium analysis [30], as defined in Equation 2.2. 
 
2.5 General Limit Equilibrium Method of Slices (GLE method)  
The GLE formulation was developed by Fredlund at the University of Saskatechewan 
[34, 35]. This method encompasses the key elements of all the other methods of slices 
(the Ordinary, Bishop’s, Janbu and M-P methods). The GLE formulation is based on 
two factors of safety equations. One equation gives the factor of safety with respect to 
19 
 
moment equilibrium, Fm (Equation 2.3) while the other equation gives the factor of 
safety with respect to horizontal force equilibrium, Ff (Equation 2.4). 
𝑐∗𝑙∗𝑅∗𝑊∗cos[𝛼𝑖]−μl)∗𝑅∗tan[𝜙]𝐹𝑚 = ∑       (2.3) 𝑊∗𝑥
𝑐∗𝑙∗cos[𝛼𝑖]+(𝑁𝑖−μ𝑙)∗tan[𝜙]∗cos[𝛼𝑖]𝐹𝑓 = ∑      (2.4) 𝑁𝑖∗sin[𝛼𝑖]
Where, 
l is the length of the bottom of the slice, R is the radius of the slip circle, W is the weight 
of slice, Ni is the total normal force on the bottom of the slice, c  is the soil cohesion, 
µ⁡is the pore water pressure, α is the inclination of the slip surface at the middle of the 
slice. 
2.6 The Ordinary or Fellenius Method (OMS) 
The ordinary method is considered the simplest of the methods of slices since it is the 
only procedure that results in a linear factor of safety equation. It is generally stated that 
the inter-slice forces can be neglected because they are parallel to the base of each slice 
[26]. This notwithstanding, the Newton's principle of 'action equals reaction' is not 
satisfied between slices. The change in direction of the resultant inter-slice forces from 
one slice to the next results in factor of safety errors that may be as much as 60% [37]. 
The normal force on the base of each slice is derived either from summation of forces 
perpendicular to the base or from the summation of forces in the vertical and horizontal 
directions. The forces considered to act in this method are represented in Figure 2.6 
below. The FoS is based on moment equilibrium and computed as: 
20 
 
 
Figure 2.6: Forces considered in the Ordinary Method of Slices 
𝑐∗𝑙+𝑁∗tan⁡()
𝐹𝑚 =         (2.5) 𝑊∗sin⁡(𝛼)
N⁡ = ⁡W ∗ cos(α)⁡– ⁡µ*l        (2.6) 
where 
W is the weight of each slice 
 µ⁡is the pore water pressure, 
l is the base length of the slice, 
α is the inclination of the slip surface at the middle of the slice. 
 
2.7 Simplified Bishop’s Method 
The simplified Bishop method neglects the inter-slice shear forces and thus assumes 
that a normal or horizontal force adequately defines the inter-slice forces as shown in 
Figure 2.7 below. The normal force on the base of each slice is derived by summing 
forces in a vertical direction.  The factor of safety is derived from the summation of 
moments about a common point as in OMS since the inter-slice forces cancel out. 
Therefore, the factor of safety equation is the same as for the ordinary method [28]. 
21 
 
However, the definition of the normal force is different. The FoS is based on moment 
equilibrium and computed as: 
 
Figure 2.7: Forces considered in the Simplified Bishop’s Method 
𝑐∗𝑙+𝑁∗tan⁡()
𝐹𝑚 =         (2.7) 𝑊∗x
𝑐∗𝑙∗sin[𝛼𝑖]−𝑢∗𝑙∗sin[𝛼𝑖]∗tan[𝜙]𝑊−[ ]
𝑁 = 𝐹𝑚𝑖 sin[𝛼𝑖]∗tan[𝜙]       (2.8) Cos[𝛼𝑖]+ 𝐹𝑚
where ‘x’ is the horizontal distance from the mid-base of the slice to the centre of 
rotation. 
As can be seen in the equation above, the expression is nonlinear due to the appearance 
of Fm in both sides of the equation and will require an iterative procedure to reach a 
solution. 
 
2.8 Janbu’s Simplified Method 
In Janbu's simplified method, the normal force in each slice is derived from the 
summation of vertical forces, with the inter-slice shear forces ignored. The horizontal 
force equilibrium equation is used to derive the factor of safety.  The sum of the inter-
slice forces must cancel and FoS equation becomes;  
22 
 
𝑐∗𝑙∗sin[𝛼𝑖]−𝑢∗𝑙∗sin[𝛼𝑖]∗tan[𝜙]𝑊−[ ]
𝑁 𝐹𝑓𝑖 = sin[𝛼𝑖]∗tan[𝜙]       (2.9) Cos[𝛼𝑖]+ 𝐹𝑓
𝑐∗𝑙∗cos(𝛼)+(𝑁𝑖−⁡μ∗𝑙)∗cos⁡(𝛼)∗tan⁡()
𝐹𝑓 =⁡∑       (2.10) 𝑁𝑖∗sin⁡(α)
 
2.9 Morgenstern-Price (M-P) Method 
This has become the most widely used method developed for analysing generalized 
failure surfaces.  The method was initially described by Morgenstern and Price [36].  
The Morgenstern‐Price method also satisfies both force and moment equilibriums and 
the overall problem is made determinate by assuming a functional relationship between 
the inter-slice shear force and the inter-slice normal force. According to Morgentern 
and price [32], the inter-slice force inclination can vary with an arbitrary function f(x) 
as: 
𝑋𝑖 = 𝐸𝑖 ∗ ⁡λ ∗ f(x)⁡        (2.11) 
where,  
f(x) = inter-slice force function that varies continuously along the slip surface and  
λ = scale factor of the assumed function. 
The method suggests assuming any type of force function, f(x), for example half‐sine, 
trapezoidal or user defined as shown in Figure 2.9. The relationships for the base normal 
force (N) and inter-slice forces (E, X) are the same as given in Janbu’s generalised 
method. The forces considered to act in this method are represented in figure 2.8 below. 
For a given force function, the inter-slice forces and the Newton-Raphson numerical 
technique can be used to solve the moment and force equations for the FoS and  until, 
Ff is equals to Fm in equations (2.12) and (2.13) [33] below.  
23 
 
 
 
 
 
 
 
Figure 2.8: Forces considered in the M-P method [24] 
𝑐∗𝑙∗cos[𝛼𝑖]+(𝑁𝑖−μ𝑙)∗tan[𝜙]∗cos[𝛼𝑖]𝐹𝑓 = ∑       (2.12) 𝑁𝑖∗sin[𝛼𝑖]
𝑐∗𝑙∗𝑅∗𝑊∗cos[𝛼𝑖]−μl)∗𝑅∗tan[𝜙]𝐹𝑚 = ∑        (2.13) 𝑊∗𝑥
where  
𝑐∗𝑙∗sin[𝛼 ]−𝑢∗𝑙∗sin[𝛼 ]∗tan[𝜙]
𝑊−∆𝑋𝑖−[ 𝑖 𝑖 ]
𝑁 = 𝐹𝑖 sin[𝛼𝑖]∗tan[𝜙]        (2.14) Cos[𝛼𝑖]+ 𝐹
F is Fm or Ff depending on which equilibrium equation is being solved. 
and  
∆𝑋𝑖 =⁡∆𝐸𝑖 ∗ λ ∗ f(x)         (2.15) 
24 
 
 
 
Figure 2.9: Inter-slice force function types 
An alternative derivation for the Morgenstern-Price method was proposed by Fredlund 
and Krahn in [34]. It presents a complete description of the variation of the factor of 
safety with respect to λ. On the first iteration, the vertical shear forces are set to zero. 
On subsequent iterations, the horizontal inter-slice forces are first computed and then 
the vertical shear forces are computed using an assumed λ value and side force function. 
Fm and Ff are solved for a range of λ values and a specified side force function. These 
25 
 
FoS are plotted on a graph as shown in Figure 2.10. The FoS vs. λ are fit by a second 
order polynomial regression and the point of intersection satisfies both force and 
moment equilibrium. 
 
 
Figure 2.10: Variation of FoS with respect to Fm and Ff vs. λ for the M-P method [34] 
 
  
26 
 
CHAPTER THREE: RESEARCH METHODOLOGY AND SOFTWARE 
USED 
 
3.1 Selection of Factor of safety method 
There are different limit equilibrium methods having varying superiorities over each 
other as discussed in Chapter 2.  Methods are based on different assumptions on 
equilibrium conditions to be satisfied.  Duncan and Wright [24] summarized some of 
limit equilibrium methods with respect to their limitations, assumptions and 
equilibrium conditions to be satisfied as shown in Table 3.1. 
Table 3.1: Brief Comparison of Limit Equilibrium Methods [24]. 
Procedure Use 
Swedish circle method Applicable to slopes where = 0.  
Applicable to non-homogeneous slopes soils where slip surface can be 
Ordinary method of slices approximated by a circle. Very convenient for hand calculations. In 
accurate for effective stress analysis with high pore water pressure. 
Applicable to non-homogeneous slopes soils where slip surface can be 
approximated by a circle. More accurate than ordinary method of slices, 
Simplified Bishop method 
especially for analysis with high pore water pressures. Calculations 
feasible by hand or spreadsheet. 
An accurate procedure applicable to virtually all slope geometries and soil 
Spencer's method profiles. The simplest complete equilibrium procedure for computing 
FoS. 
An accurate procedure applicable to virtually all slope geometries and soil 
Morgenstern and Price's 
profiles. Rigorous, well-established complete equilibrium procedure.  
method 
Requires solution of nonlinear equations with an iterative procedure. 
The analysis and design of failing slopes requires an in-depth understanding of the 
failure mechanism in order to choose the right slope stability analysis method.  
In this research, as stated above, the slice approach for GLE procedure is used and the 
FoS according to the Morgenstern-Price’s procedure, which satisfies all the 
27 
 
requirements for static equilibrium is adopted.  Regardless of whether equilibrium is 
considered for a single free body or a series of individual vertical slices, there are more 
unknowns (forces, locations of forces, factor of safety, etc.) than the number of 
equilibrium equations; the problem of computing a FoS is thus statically indeterminate. 
Therefore, assumptions must be made to achieve a balance of equations and unknowns. 
This method allows for analysis of any failure shape (circular, non-circular or 
compound). The solution for the (FoS) is derived from the summation of forces 
tangential and normal to the base of a slice and the summation of moments about the 
centre of the base of base slice. The steps required to provide the input data for 
performing the slope stability analysis include [37]. 
 A survey of the elevation of the ground surface on a section perpendicular to 
the slope. 
 Estimation of ground stratigraphy from borehole logs and soil/rock properties 
from engineering soil/rock tests.   
 The determination of ground water level from piezometer readings to estimate 
ground water pore-water pressures. 
 
3.2 Morgenstern-Price Method 
Morgenstern and Price [32] developed the method for FoS similar to the Spencer 
method [31]. The method considers both normal inter slice force and shear forces. 
Therefore it satisfies both moment and force equilibrium.  
 
 
28 
 
3.3 Assumptions 
To begin with the generalized formulation, the following assumptions are made for a 
body of mass slope: 
1. The failure surface is assumed to be circular. 
2. The soil/rock is a homogeneous. 
3. The soil/rock is an isotropic material. 
4. The failure mass is a rigid body. 
5. The base normal force acts at the middle of each slice. 
6. The Mohr-Coulomb failure criterion is used. 
 
3.4 Numerical Method Development  
In this work, the geotechnical generalized method of slices approach, which is a form 
of the Finite Strip Method would be used for the computation of the Factor of Safety 
(FoS), which is a measure of the degree of safety of the slope. The Finite Strip Method 
is a variant form of the known differential equation numerical discretization 
approaches.  It discretizes the equilibrium differential equation (based on slices) and 
imposes equilibrium boundary conditions at the ends to solve for the required 
unknowns. 
The method of slices method is a numerical approach used to solve the equilibrium 
differential equation at the limiting condition, i.e., at the point of slope failure.  The 
method consists of dividing the slope into a number of fine slices so that their bases can 
be comparable with a straight line.  The equilibrium differential equations are then 
developed for each slice, and then the global force and moment equilibrium problem 
solved numerically to obtain the FOS [24]. 
29 
 
This study will make available a locally developed slope stability program that would 
simulate the behaviour of a simple slope undergoing circular failure.  The results from 
this work, which will include both factor of safety for moment and force equilibrium, 
will be visualized in the command window of MATLAB. 
 
3.5 Derivation of Equations  
Figure 3.1. is a sketch of a slope section and Figure 3.2 is a sketch of a slice of mass 
from Figure 3.1 with the forces acting on it at the point of failure. By taking moment 
about the midpoint of the base of the slice in Figure 3.1; 
′ ′ dy ′ ′ ′ ′ dy dx𝐸 [(𝑦 − 𝑦𝑡) − (− )] − (𝐸 + dE ) [𝑦 + dy − 𝑦𝑡 − dy𝑡 + ( )] − 𝑋 ( ) − (𝑋 +2 2 2
dx
dX) ( ) = 0⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡⁡  (3.1) 
2
After simplification and proceeding to the limit as dx→0, it can be shown readily that; 
𝑑(𝐸′.𝑦′𝑡) dE
′
𝑋 = − 𝑦( )        (3.2) 
dx dx
 
Figure 3.1: Sketch of a Slope Section [23]. 
30 
 
 
 
 
 
 
 
Figure 3.2: Forces acting on a Single Slice from a Mass Slope [24]. 
For equilibrium in the N direction, we have, from Figure 3.2; 
dN′ = dWcos[𝛼] − dXcos[𝛼] − dE′sin[𝛼]⁡     (3.3) 
For equilibrium in the S direction, we have from Figure 3.2; 
dS = dE′cos[𝛼] − dXsin[𝛼] + dWsin[𝛼]⁡⁡⁡⁡⁡⁡    (3.4) 
The Mohr-Coulomb failure criterion in terms of effective stresses may be expressed as; 
1
dS = {𝑐′dxsec[𝛼] + (dN′)tan[𝜙′]}⁡     (3.5) 
𝐹
Equation (3.5) defines the factor of safety in terms of shear strength. 
Eliminating dS from Equations (3.4) and (3.5); 
1
[𝑐′dx⁡sec[𝛼] + (dN′)sin[𝜙′]] = dE′cos[𝛼] − dXsin[𝛼] + dWsin[𝛼] (3.6) 
𝐹
By eliminating dN′ from Equations (3.3) and (3.6) and dividing by dxcos[𝛼], we have; 
𝑐′ 2 tan[𝜙
′] dW dX dE′ dE′ dX dW
sec [𝛼] + [ − − tan[𝛼]] = − + tan[𝛼]  (3.7) 
𝐹 𝐹 dx dx dx dx dx dx
dy
In the specified co-ordinate system, tan[𝛼]=-  and equation (3.7) becomes; 
dx
31 
 
𝑐′ dy 2 tan[𝜙′] dW dX dE′ dy dE′ dX dW dy
[1 + ( ) ] + [ − + ( )] = − − ( )  (3.8) 
𝐹 dx 𝐹 dx dx dx dx dx dx dx dx
Therefore, the governing differential equation becomes using equations (3.2); 
dE′ tan[𝜙′] dy dX tan[𝜙′] dy 𝑐′ dy ′]
[1 − ( )] + [ + ] = [1 + ( )2
dW tan[𝜙 dy
] + [ + ](3.9) 
dx 𝑓 dx dx 𝐹 dx 𝐹 dx dx 𝐹 dx
E’ and E’+dE’ are the horizontal inter-slice forces, 
X and X+dX are the vertical inter-slice forces, 
𝜙′ is the effective friction angle, 
W is the weight of ith slice, 
dN and dS are resultants of the normal and tangential forces on the  slice base of length 
li  
α is the inclination of the slice with respect to the horizontal. 
F is the factor of safety. 
The governing differential equation is developed for each slice, and summed up to 
develop the equilibrium equation for entire mass.  Based on the assumption on the inter-
slice forces as provided in Chapter 2, the various expressions for Fm, Ff, N, X & E for 
the different slices are computed. 
 
3.6 Structured Program 
Below (Appendix A) is an algorithm developed to solve for the force and moment factor 
of safety’s (Ff and Fm) using the M-P method.  
 
 
32 
 
3.7 Numerical Algorithm 
Figure 3.3 shows the numerical algorithm used to solve the problem. Equation numbers 
referred to here are that from appendix A. 
 
Generate ytop 
 Input H, and ybot eq 
Xc, Yc, R  
1&2 
 
Compute lb, W and alpha-eq 
 10, 7& 6  
 Input c and phi 
Compute FSom-eq 15 &16 
 
Set lambda and 
 initial FoS for the 
first iteration 
 Compute N, FSm, FSf and 
∆E-eq 17, 15, 16 &11. 
No  
Is FSm-FoS<tol  
Is FSf-FoS <tol else 
NFSm=FSm, NFSf=FSf 
 
Yes 
 
Output 
FSm(i), Increase Output 
FSf(i) FSm, FSf   
 
 
Figure 3.3: Algorithm flowchart for solving for the FoS 
 
 
 
33 
 
3.8 Software and Programs Used 
3.8.1  SLOPE/W  
SLOPE/W, developed by GeoStudio International Canada, is part of a professional 
geotechnical software suite of programmes.  SLOPE/W is the programme used for slope 
stability analysis.  This software is based on the theories and principles of the LE 
methods discussed in the previous sections.  To check the accuracy of the program 
written in MATLAB and its output, SLOPE/W is employed.  For this study, a full 
licensed version of Geostudio 2007 has been used.  
For each model, using the drawing tools, the geometry of the slope is entered. Then by 
using the slip surface dialogue box, the slip surface is specified by using a range. 
Then by using the “Materials” dialogue box, soil parameters will be entered and the 
selected soil will be assigned to the drawing in the software. 
After entering all of the input data into the software by hitting the “Start” button under 
the “Solve Manager”, the program starts to analyze the slope and find the minimum 
factor of safety and its related failure surface. Figures 3.4 to 3.6 show some of the stages 
of developing the model using SLOPE/W  
34 
 
 
Figure 3.4: SLOPE/W KeyIn Analyses Page 
 
Figure 3.5: SLOPE/W KeyIn Material Page 
35 
 
 
Figure 3.6: SLOPE/W KeyIn Entry and Exit Range Page 
 
3.8.2 MATLAB 
MATLAB (Matrix Laboratory) is a multi-paradigm numerical computing environment 
and fourth-generation programming language. A proprietary programming language 
developed by MathWorks, MATLAB allows matrix manipulations, plotting of 
functions and data, implementation of algorithms, creation of user interfaces, and 
interfacing with programs written in other languages, including C, C++, Java, Fortran 
and Python [38].  The programme is also easily accessible and is available in most 
design engineering offices.  These features of MATLAB make it suitable to be used for 
this research work.   
  
36 
 
CHAPTER FOUR: RESULTS AND DISCUSSION 
 
4.1 Introduction 
This chapter presents outputs from the MATLAB code and SLOPE/W programmes and 
compares the results of the two. First, the analysis is run for assumed models (the same 
geometry, soil properties) of homogeneous soil mass and modelled for two failure 
scenarios; i.e., for toe and base failure conditions. 
The second set of comparison is done for existing case studies. Case one is the study 
from Griffiths and Lane [39] and Rocscience [40] and case two is from the study by 
ACADS [41].  
The MATLAB script is presented in appendix B. 
 
4.2 Programme Test Examples 
The program was tested with two failure mechanisms, toe and base failures. In slope 
stability design and analysis, and for most design projects, the worst case scenario is 
always adopted. The worst case when it comes to FoS in slope stability analysis is the 
least FoS. This is adopted in choosing the FoS in this MATLAB program in situations 
when more than one FoS is outputted. 
 
4.2.1 Toe Failure 
Two models were presented for the toe failure mechanism. The slope angles for both 
models are 45° or 1H: 1V (H: horizontal distance; V: vertical distance). The rest of the 
soil parameters and geometry are presented in Tables 4.1 and 4.2. 
37 
 
Case 1 
Table 4.1: Slope Dimensions and Material Properties for Toe Failure: Case 1 
 
c´ ’ γ H 
 (kPa) (ͦ) (kN/m3) (m) 
Case 1 15 38 18 5 
 
 
 
Figure 4.1: SLOPE/W Output of Toe Failure: Case 1 
  
38 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Figure 4.2: MATLAB Output of Toe Failure: Case 1 
  
39 
 
Case 2 
Table 4.2: Slope Dimensions and Material Properties for Toe Failure: Case 2 
 
c´ ’ γ H 
 (kPa) ()ͦ  (kN/m3) (m) 
Case 2 5 28 16 5 
 
 
Figure 4.3: SLOPE/W Output for Toe Failure: Case 2 
  
40 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
  
 
Figure 4.4: MATLAB Output for Toe Failure: Case 2 
  
41 
 
It is noted from the outputs that the MATLAB code give reasonable outputs of FoS 
compared to SLOPE/W. In both cases, the percentage error is found to be 4%. The 
small error margin is an indication that the developed MATLAB code is working to an 
acceptable degree. 
 
4.2.2 Base Failure 
Two models were presented for the base failure mechanism. The slope angles for both 
models are 1H: 1V. The rest of the soil parameters and geometry are presented in Tables 
4.3 and 4.4. 
Case 1 
Table 4.3: Slope Dimensions and Material Properties for Base Failure: Case 1 
 
c´ ’ γ H 
 (kPa) ()ͦ  (kN/m3) (m) 
Case 1 40 20 20 5 
Figure 4.5: SLOPE/W Output for Base Failure: Case 1  
42 
 
 
  
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
  
 
Figure 4.6: MATLAB Output for Base Failure: Case 1 
  
43 
 
Case 2 
Table 4.4: Slope Dimensions and Material Properties for Base Failure: Case 2 
 
c´ ’ γ H 
 (kPa) ()ͦ  (kN/m3) (m) 
Case 2 0 15 18 5 
 
Figure 4.7: SLOPE/W Output for Base Failure: Case 2 
  
44 
 
 
 
 
 
 
 
 
 
  
 
 
 
 
 
 
 
 
 
 
 
Figure 4.8: MATLAB Output for Base Failure: Case 2 
  
45 
 
4.3 Comparison with Cases from Literature 
 
4.3.1 Case 1 
The geometry of a homogeneous slope without foundation is shown in Figure 4.9. This 
case follows the analyses performed by [39] and [40], which were used as benchmark 
cases to study the applicability of 3-D FE analyses to slope stability. The slope angle is 
26.25° or 2H: 1V for this case. According to Griffiths and Lane, the adopted parameters 
are based on cʹ/γH = 0.05. The height of the slope is assumed to be 40 m, thus the 
corresponding parameters are summarized in Table 4.5. 
 
Figure 4.9: Homogeneous Slope without Foundation 
Table 4.5: Slope dimensions and material properties 
c´ ’ γ H 
(kPa) (ͦ) (kN/m3) (m) 
40 20 20 40 
The slope stability analyses using the computer program, SLOPE/W, is shown in Figure 
4.10. The FoS based on the GLE and the adopted M-P approach is conducted for both 
toe and slope failures. The FoS for toe and slope failures are found to be 1.385 and 
1.373 respectively. 
46 
 
The same analysis was conducted using the MATLAB code, the FoS for toe and slope 
failures are found to be 1.431and 1.375 respectively and shown in Figure 4.10-4.11. 
 
Figure 4.10: Analysis using SLOPE/W - FoS = 1.385 
 
 
 
 
 
 
 
 
 
 
 
 
47 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
  
 
Figure.4.11: FoS based on M-P approach for Toe failure - FoS = 1.431  
48 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
  
Figure 4.12: FoS based on M-P approach for Base Failure: FoS = 1.375 
  
49 
 
4.3.2 Case 2 
In 1988 a set of 5 basic slope stability problems, together with 5 variants, was 
distributed both in the Australian Geomechanics profession and overseas as part of a 
survey sponsored by ACADS [41].  
This problem is taken from the verification manual of Slide 3.0 (Verification #1). It was 
distributed to the Australian Geomechanics profession and overseas in 1988 with some 
other models to verify the efficiency and accuracy of Slide 3.0. The slope model 
geometry is presented in Figure 4.13.  The slope material properties are shown in Table 
4.6. 
 
Figure 4.13:  Slope Model Geometry from Slide 3.0 [41] 
Table 4.6: Slope Dimensions and Material Properties for ACAD Problem 
c´  γ H 
(kPa) ()ͦ  (kN/m3) (m) 
3 19.6 20.2 10 
 
  
50 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Figure 4.14: FoS based on M-P Approach for ACAD Problem. 
  
51 
 
The FoS from this case study when modelled with Slide 3.0 and the GLE (M-P) 
approach adopted, was found to be 0.986. Upon modelling the same problem with the 
MATLAB code, as shown in Figure 4.14 above, the FoS is found to be 1.035. Again, 
this can be said to be a reasonably accurate estimate of the FoS by the MATLAB code. 
The percentage error is 4.7%. 
 
  
52 
 
CHAPTER FIVE: CONCLUSIONS 
 
5.1 Conclusions 
The primary focus of this research was to:  
 Develop a MATLAB code for solving the FoS of homogeneous earth slopes  
 Verify and validate the code with SLOPE/W models and also with cases from 
literature  
Seven cases were studied in this research as summarized below in Table 5.1 
Table 5.1: Summary of FoS Outputs for all Case Studies. 
Percentage 
Case 1 Case 2 
  Error (%) 
  SLOPE/W MATLAB SLOPE/W MATLAB Case 1 Case 2 
Toe Failure 5.649 5.406 3.659 3.514 4 4 
Base 
Failure 3.152 3.168 0.611 0.614 0.5 0.5 
Comparison with Literature 
 Griffiths Slide 
and Lane MATLAB (Rocscience) MATLAB Case 1 Case 2 
Toe Failure 1.385 1.431   3.2  
Base 
Failure 1.373 1.375 0.986 1.035 0.1 4.7 
From the Table 5.1 above, it is seen that the MATLAB program gives a good estimate 
of the FoS when compared with SLOPE/W and some other problem evaluated with 
different programs for a homogeneous soil material.  Two different failure mechanisms 
are shown (i.e., toe and base failures).  It is noted that the maximum percentage error 
in the M-P FoS is 5%.  
It is anticipated that with the availability of this computer code, Ghanaian Engineers 
can more readily assess the safety of slopes in routine design works. 
53 
 
5.2 Recommendations 
The current work could be extended in the future to include the following.   
 The provision for a more "user friendly" interface and development of the 
programme into a stand-alone interactive program. 
 The direct extension of the model to cater for heterogeneous soil and rock 
material 
 The direct extension of the model to include earthquake effect, surcharge 
loading and increase in pore water pressures.   
It is hoped that at the end of these improvements, this work that has started during this 
research study, would be found in most geotechnical design offices in Ghana and other 
less-developed countries, who cannot afford the existing expensive commercial 
software available.  
  
54 
 
REFERENCES 
[1] Wikepedia, https://en.wikipedia.org/wiki/Slope_stability, accessed June 2016. 
[2] Shioi Y. and Sutoh S., Collapse of high embankment in the 1994 far-off Sanriku 
Earthquake. Slope Stability Engineering, Eds. Yagi, Yamagami & Jiang, Balkema, 
Rotterdam, pp.559-564, 1999. 
 [3] Abramson, L. W.; Lee, T. S., Sharma, S., and Boyce, G. M., Slope Stability and 
Stabilisation Method, 2nd edition, John Wiley & Sons, New York, 2002. 
[4] Rouaiguia, A., Dahim, M. A., Numerical modeling of slope stability analysis, 
IJESIT, Vol. 2, No 3, May 2013. 
[5] Ayetey, J. K., The extent and effects of mass wasting in Ghana, Bulletin of the 
International Association of Engineering Geology, No. 43, Paris, 1991. 
[6] CitifmOnline, http://citifmonline.com/2014/06/09/death-trap-rocks-hang-
dangerously-along-ayi-mensah-peduase-road, accessed May, 2016. 
[7] GraphicOnline, http://graphic.com.gh/news/general-news/22493-danger-
looms-on-aburi-accra-road-as-rocks-fall-from-cliff.html, accessed, May 2016 
[8] GraphicOnline, http://graphic.com.gh/news/general-news/8384-landslide-can-
occur-around-akoasa-mountain.html, accessed, May 2016 
[9] DeNatale, J. S., Rapid identification of critical slip surfaces: structure, Journal 
of Geotechnical Engineering, ASCE, Vol. 117, No. 10, pp. 1568-1589, 1991. 
 [10] Drucker, D. C. and Prager, W., Soil mechanics and plastic analysis or limit 
design. Quarterly Journal of Applied Mathematics, Vol. 10, pp. 157-165, 1952 
[11]  Zhang, X., Slope stability analysis based on the rigid finite element method. 
Geotechnique, Vol. 49, No. 5, pp. 585-593, 1999. 
55 
 
[12] Kim, J., Salgado, R., and Lee, J., Stability analysis of complex soil slopes using 
limit analysis. Journal of Geotechnical and Geoenvironmental Engineering, ASCE, 
Vol. 128, No. 7, pp. 546-557, 2002. 
[13] Chen, W. F., Limit Analysis and Soil Plasticity, Elsevier Scientific Publishing 
Co., New York, 1975. 
[14] Karal, K., Energy method for soil stability analyses. Journal of the Geotechnical 
Engineering Division, Vol. 103, No. 5, pp. 431-445, 1977. 
[15] Izbichki R. J., Limit plasticity approach to slope stability problems. Journal of 
the Geotechnical Engineering Division, Vol. 107, No. 2, pp. 228-233, 1981.  
[16] Michalowski, R. L., Slope stability analysis: a kinematical approach. 
Geotechnique, Vol. 45, No. 2, pp. 283-293, 1995. 
[17] Donald, I. and Chen, Z. Y., Slope stability analysis by the upper bound 
approach: fundamentals and methods. Canadian Geotechnical Journal, 34, 853-862, 
1997.  
[18] Baker, G., and Garber, M., Theoretical analysis of the stability of slopes, 
Geotechnique, Vol. 28, No. 4, pp. 395-411, 1978. 
[19] Jong, G. D., Application of the calculus of variations to the vertical cut off in 
cohesive frictionless soil, Geotechnique, Vol. 30, No. 1, pp. 1-16, 1980. 
[20] Cheng, Y. M., Zhao, Z. H., and Wang, J. A., Realization of Pan Jiazheng’s 
extremum principle with optimization methods. Chinese Journal of Rock Mechanics 
and Engineering, 27(4), pp. 782-788 (in Chinese), 2008. 
[21] Pan, J. L., Goh, A. T. C., Wong, K. S., and Selby, A. R. 2002. Three-
dimensional analysis of single pile response to lateral soil movements. International 
56 
 
Journal for Numerical and Analytical Methods in Geomechanics, 26(8), pp. 747-758, 
2002. 
[22] Duncan, J.M., State of art: limit equilibrium and finite-element analysis of 
slopes, Journal of Geotechnical Engineering, Vol. 122, No. 7, pp. 557-596, 1996. 
[23] Duncan, J.M. and Wright, S.G., Soil Strength and Slope Stability. John Wiley 
& Sons Inc., 2005. 
[24] Zhao, Y., Tong, Z., Lü, Q., Slope Stability analysis using slice-wise factor of 
safety, Mathematical Problems in Engineering, Volume 2014, Paper No. 712146, 6 
pages, 2014. 
[25] Petterson, K.E., The early history of circular sliding surfaces, Geotechnique, 
Vol. 5:pp. 275-296, 1955. 
[26] Krahn J., Stability Modeling with SLOPE/W: An Engineering Methodology, 
July 2012 Edition, Geo-Slope International, 2012. 
[27] Fellenius, W., Calculation of stability of earth dams. Proceedings of the 2nd 
Congress Large Dams pp. 445-463, Washington, D. C., 1936. 
[28]  Bishop, A. W., The use of the slip circle in the stability analysis of slopes. 
Geotechnique, Vol. 5, No. 1, pp. 7-17, 1955. 
[29] Janbu, N., Slope Stability Computations, Embankment-Dam Engineering, 
Casagrande Volume, Ed. R. C. Hirschfield and S. J. Poulos, John Wiley and Sons, New 
York, pp. 47-86, 1973. 
[30] Janbu, N., Application of composite slip surface for stability analysis, European 
Conference on Stability of Earth Slopes, Vol. 3, pp.39-43, Stockholm, 1954. 
57 
 
[31] Spencer, E., A method of analysis of embankments assuming parallel inter-slice 
forces-Geotechnique, Vol. 17(1), pp. 11-26, 1967. 
[32] Morgenstern, N. R. and Price, V. E., The analysis of stability of general slip 
surfaces. Geotechnique, Vol. 15, No. 1, pp. 77-93, 1965. 
[33] Nash, D., A comparative review of limit equilibrium methods of stability 
analysis, Slope Stability, Ed. M. G. Anderson and K. S. Richards, pp. 11 -75, 1987. 
[34] Fredlund, D.G, and Krahn, J., Comparison of slope stability methods of 
analysis, Canadian Geotechnical Journal, Vol. 14, No. 3, pp. 429-439, 1977.  
[35] Fredlund, D.G, Krahn, J., and Pufahl, D.E., The relationship between limit 
equilibrium slope stability methods, Proceedings of the International Conference on 
Soil Mechanics and Foundation Engineering, Vol. 3, pp.409-416, 1981. 
[36] Whitman, R. M., and Bailey, W. A., Use of computers for slope stability 
analysis, Journal of the Soil Mechanics and Foundation Division, ASCE, Vol. 93, No. 
SM4, pp. 475-498, 1967. 
[37] Baba, K., Bahi, L., Ouadif, L., and Akhssas, A., Slope stability evaluations by 
limit equilibrium and finite element methods applied to a railway in the Moroccan Rif, 
Open Journal of Civil Engineering, Vol. 2, pp. 27-32, 2012. 
[38] Matlab, https://en.wikipedia.org/wiki/MATLAB, accessed April 3, 2016 
[39] Griffiths, D.V., and Lane P.A., Slope stability analysis by finite elements, 
Geotechnique, Vol. 49, No. 3, pp. 387-403, 1999. 
[40] Rocscience Inc., Application of the Finite Element Method to Slope Stability. 
Rocscience, 2004. 
58 
 
[41] ACADS, 2D elasto-plastic finite element program for slope and excavation 
stability analyses - Slope Stability Verification Manual, Rocscience, 2011.  
59 
 
APPENDICES 
Appendix A: Structured Programme 
Below is an algorithm developed to solve for the force and moment factor of safety’s 
(Ff and Fm) using the M-P method.  
1. By using the equation of a straight line, define the slope surface. 
𝑦top = 𝑚 ∗ 𝑥 + 𝑐        (1) 
2. By using the equation of a circle, define the slip surface. 
𝑦 = 𝑦 − √(abs(𝑅2 − (𝑥 − 𝑥 )2bot 𝑐 𝑐       (2) 
3. By using the equation of a straight line, define the piezometric line. 
𝑦w = 𝑚 ∗ 𝑥 + 𝑐        (3) 
4. Divide the slope into n number of slices by vertical lines. 
5. For each slice, the width, dx, bottom inclination, α, and average height, havg., are 
determined. 
ℎ +ℎ
ℎavg =
𝑖 𝑖+1         (4) 
2
𝑋
dx =          (5) 
𝑛𝑠−1
(ℎ −ℎ )
𝛼 = tan−1( 𝑖+1 𝑖 )        (6) 
dx
6. The area of the slice, A, is computed by multiplying the width of the slice (dx) 
by the average height, havg.  
7. The weight, W, of the slice is computed by multiplying the area of the slice by 
the total unit weight of soil:  
8. W = γA         (7) 
9. If piezometric height is above slip surface, continue with step 9, else, skip to 
step 17. 
60 
 
10. The piezometric height is determined at the downslope boundary, centre and 
upslope boundary of each slice. The piezometric height at the downslope and 
upslope boundaries of the slice, hb and ht, respectively, are used to compute the 
forces from water pressures on the sides of the slice. Here, a triangular 
hydrostatic distribution of pressures is assumed on the sides of the slice. The 
piezometric height at the centre of the slice, hp, represents the pressure head for 
pore water pressures at the base of the slice. 
11. Hydrostatic forces from water pressures on the sides of the slice are computed 
from, 
ℽ𝑤∗ℎ
2
µ = 𝑏𝑏          (8) 2
ℽ𝑤∗ℎ
2
µ 𝑡𝑡 =          (9) 2
12. The pore water pressure is computed by multiplying the piezometric head at the 
centre of the base of the slice by the unit weight of water: µ = Ɣw hp  
13. The length of the bottom of the slice, l, is determined; the length can be 
computed from the width, dx, and base inclination, 
𝑑𝑥
𝑙 =         (10) 
cosα
14. A trial FoS (F) is assumed and a λ ranging from 0 to 1 set. 
15. Beginning with the first slice at the toe; 
tan[𝜙]cos[𝛼 ] cl
𝑊[sin[𝛼 ]− 𝑖𝑖 ]−𝐹 𝐹
∆𝐸𝑖 =       (11) tan[𝜙]sin[𝛼 ]
cos[𝛼𝑖]+
𝑖
𝐹
16. The inter-slice shear force, X, for each slice is calculated from equation…. 
∆𝑋𝑖 = ∆𝐸𝑖λf(𝑥)       (12) 
61 
 
Where  λ is the percentage (in decimal form) of the function used and f(x) is 
inter-slice force function representing the relative direction of the resultant 
inter-slice force given by 
𝑋 −𝑋
𝑓(𝑥) = sin[( 𝑖+1 𝑖) ∗ 𝜋]      (13) 
𝑋𝑛−𝑋𝑖
17. The normal force, N, at the base of each slice is calculated from, 
𝑐∗𝑙∗sin[𝛼𝑖]−𝜇∗𝑙∗sin[𝛼𝑖]∗tan[𝜙]𝑊−(∆X)−[ ]
𝑁 = 𝐹𝑖 sin[𝛼 ]∗tan[𝜙]      (14) 
cos[𝛼𝑖]+
𝑖
𝐹
18. The new factor of safety’s, Ff and Fm , are calculated from 
𝑐∗𝑙∗cos[𝛼𝑖]+(𝑁𝑖−μ𝑙)∗tan[𝜙]∗cos[𝛼 ]F 𝑖f⁡= ∑      (15) 
𝑁𝑖∗sin[𝛼𝑖]
𝑐∗𝑙∗𝑅∗𝑊∗cos[𝛼𝑖]−μl)∗𝑅∗tan[𝜙]            𝐹𝑚 = ∑      (16) 𝑊∗𝑥
19. Iterate for a number of slip surfaces and determine the slip surface with the least 
FoS. 
62 
 
Appendix B: MATLAB Code for Solving FoS 
%%%this code calculates the factor of safety of a homogeneous soil mass 
%%%using the method of slices and the M-P interslice force function 
Clear all; close all; clearvars; clc;  
%% Slope and soil parameters 
H   = 5; %inpu('Enter the height of the slope');                % height of 
slope 
Xc =1; %inpu('Enter the x coordinate of center of slip circle: Xc should be 
>0');  
% x coordinate of center of circle 
Yc = 7; %input ('Enter the y coordinate of center of slip circle: Yc should 
be greater than H');               % y coordinate of center of circle 
R   = 8; %:5: Yc+10      % radius of center of slip circle 
fr = 45;              % slope angle in degrees 
ns = 31;              % number of slice faces 
G   = 20;               % unit weigh of soil 
c   = 40;              % cohesion of soil 
phi = 20;               % frictional angle of soil (degrees) 
tph = tand(phi);        % tangent of the frictional angle of soil 
tfr = tand(fr);         % slope angle in gradient 
%% Slip surface generation 
xmin = Xc - sqrt(abs(R^2 - Yc^2))             % Exit point of slip circle 
xmax = Xc + sqrt(abs(R^2 - (Yc-H)^2))         % Entry point of slip circle 
x    = linspace (xmin,xmax,ns)';                 % Positions where the forces 
will be analyzed. 
ybot = Yc - sqrt(abs(R^2-(x-Xc).^2));           % Equation for generating slip 
circle 
%% Slope surface generation 
ytop = (x>=0).*(x<(H/tfr)).*(x*tfr) + (x>=(H/tfr)).*H; % Slope surface. 
equation of a line is used 
figure(1) 
hold on 
plot (x,ybot,'-r',x,ytop,'-r') 
xlabel ('Distance (m)') 
ylabel ('Elevation (m)') 
plot(reshape([x,x,x].',1,[]),reshape([ytop,ybot,ytop*NaN].',1,[])) 
hold off 
%% Other parameters 
hs    = ytop-ybot;                  % Height at each node 
havg = (hs(1:end-1)+hs(2:end))/2;  % 5) Average height of each slice 
xavg  = abs(((x(1:end-1)+x(2:end))/2)-Xc);    %    midpoint distance of each 
slice. 
dx    = (xmax-xmin)/(ns-1);                 %    Width of each slice 
yb    = (diff(ybot)); 
alpha = rad2deg(atan2(yb,dx));        %    Slice angle. With atan2 the sign 
is clear. 
sal = sind(alpha); 
cal = cosd(alpha); 
A     = dx.*havg;                   % 6) Area of each slice. 
W     = G*A;                        % 7) Weight of each slice. 
s= R*sal;                    % offset f for non-circular slip surfaces 
l     = dx. /cal;            % 12) Length of the bottom of the slice, assuming 
straight border. 
 u     = 0;  
 %% factor of safety computation 
upO = sum((c.*l) + (W.*cal*tph) - (u.*l*tph)); 
downO = sum(W.*sal); 
FSom = upO. /downO 
OFSm = FSom; 
NFSm = 1.2*FSom; %initial guessed F 
OFSf = FSom; 
NFSf = 1.2*FSom; 
Tol =0.001; 
t =-1:0.2:0; 
y=t'; 
63 
 
i=1; 
k=1; 
n=1; %: ns-1; 
for lbd=t; 
        E (1) =0; 
        d=1: ns-1; 
        f=sind (((xavg (d)-xmin)/ (xmax-xmin))*180); 
        X=E.*lbd.*f; 
    while abs(OFSm-NFSm)>tol 
        OFSm = NFSm; 
        Nbu = ((W-X) - ((c.*l.*sal) - (u.*l.*tph.*sal)). /OFSm); 
        Nbd = cal + (tph.*sal). /OFSm; 
        Nob = Nbu. /Nbd;        
        upb = sum((c.*l*R) + (Nob.*tph*R) - (u.*l*tph*R)); 
        downb = sum(W.*xavg); 
        FSm = upb. /downb; 
        FSmm = ((FSm)) %./(ns-1)); 
        NFSm = FSmm; 
        FSM (i)=FSmm; 
     end  
  
    while abs(OFSf-NFSf)>tol           
        OFSf = NFSf;         
        Nju = ((W-X) - ((c.*l.*sal) - (u.*l.*tph.*sal)). /OFSf); 
        Njd = cal + (tph.*sal). /OFSf; 
        Noj = Nju. /Njd;     
        upj = sum( (c.*l.*cal) + (tph*cal.*Noj) - (u.*l.*tph.*cal) ); 
        downJ = sum(Noj.*sal); 
        FSf = (upj. /downJ); 
        FSf = (FSf) %./(ns-1)); 
        NFSf = FSf; 
        FSF (i)= FSf; 
    end    
    E = (((c.*l-u.*l.*tph).*cal). /OFSm) + (sal.*Nob)-((tph.*cal).*Nob./OFSm); 
    E = (((c.*l-u.*l.*tph).*cal). /OFSm) + (sal.*Noj)-((tph.*cal).*Noj./OFSm); 
  
    OFSm = FSom; 
    OFSf = FSom; 
    i=i+1; 
     
end 
fSM=FSM' 
fSF=FSF' 
figure(2) 
plot(y,fSM,'-r',y,fSF,'-b') 
%%    finding intersection point 
[xi, yi] = polyxpoly(y, fSM, y, fSF); 
xlabel('lambda') 
ylabel('FoS') 
legend('FSm', 'FSf') 
mapshow(xi,yi,'DisplayType','point','Marker','o', 'MarkerEdgeColor','k'); 
axis normal 
  
strValues = strtrim(cellstr(num2str([xi yi],'(%0.3f,%0.3f)'))); 
text (xi,yi,strValues,'VerticalAlignment','bottom'); 
[xi yi] % Display intersection points 
 
 
64