Slope stability analyses using limit equilibrium and strength reduction methods

This paper presents results of comparative slope stability analyses conducted by limit equilibrium and strength reduction methods. Several slopes, taken either from geotechnical practice or literature review, are considered. The influence of tension, distributed load, finite element size and model parameters is analysed in relation to the location, shape of the critical failure surface and the corresponding value of the factor of safety. Both methods provide similar results if they are correctly applied using appropriate software programs.


Introduction
In the engineering practice, slope stability analyses are most commonly performed using limit equilibrium methods (LEM).The most widely used are the methods of slices, where the soil mass above an assumed failure surface is divided into several, usually vertical, slices and, as a result, the value of the factor of safety is obtained: (1) where τ f is the actual available shear strength of the material, and τ m is the mobilised shear strength or the average shear stress on the hypothetical failure surface mobilized to maintain the body in equilibrium.Despite their inherent weaknesses [1][2][3], these methods have been developed and tested by means of actual case histories.In order to render the analysis statically determinate, some assumptions must be made.The assumption is made that the interslice shear force X is related to the interslice normal force E (total or effective), by mathematical expression: where l is the scaling constant representing the percentage of the function f(x) used for solving for the factor of safety equation, and f(x) is the functional relationship (interslice force function) that describes the manner in which the magnitude of X/E varies across the failure surface.The final result is usually not very sensitive to the choice of f(x) function [4][5][6], but examples where the solution is very sensitive to the assumed interslice force function have also been reported in the literature, e.g. by Krahn [2].Although limit equilibrium methods adopt the general philosophy of an upper bound solution [7], they do not meet all accuracy requirements.Izbicki [8] defined the results of limit equilibrium methods as a "reduced" upper bound, implying that the LEM factor of safety will be slightly lower than the upper bound factor of safety.Bishop's simplified method [9] is used to calculate the factor of safety of a circular failure surface, which has proven to provide the results similar to those obtained with some more rigorous formulations.As to the methods that satisfy all static equilibrium elements, the Maksimovic's method (1979), Spencer method (1967) and GLE method (Fredlund and Krahn, 1977) are used.The position and shape of the critical failure surface is determined by the use of the semi and fully automatic search techniques.The interactive optimization algorithm implemented in the software package BGSLOPE [10,11] is used to obtain the critical failure surface described by Bezier or a polygonal curve.This procedure is found to be the most versatile for comparison with the finite element method results, as it gives the user the possibility to successively manage the movements of control points of Bezier polygon or coordinate points of a polygonal curve.The analyses in LEM are performed by assuming the linear Mohr-Coulomb failure criterion.In this way, the input of no more than three parameters is needed, namely the total unit weight of material (γ), cohesion (c), and angle of shearing resistance (j).The finite element method (FEM) is a numerical technique for finding approximate solutions to boundary value problems for partial differential equations.It theoretically satisfies all requirements that must be met for a complete solution to a slope stability problem [12].The material behaviour in FEM is described by an elastic, perfectly plastic model complying with the Mohr-Coulomb failure criterion [13].This model takes into account shear strength and deformation parameters.Three additional parameters, along with the aforementioned ones, are the modulus of elasticity, E, Poisson's ratio, ν, and angle of dilatancy, ψ.Several authors [14][15][16] have shown that the deformation parameters E and ν, as well as the domain size, have an insignificant influence on the factor of safety value.The effect of dilatancy on the final result has also been investigated by several authors [16,17], and is reanalysed in this paper.In case of non-associated plasticity, the positive angle of dilatancy is calculated based on a simple, widely accepted, relationship (ψ = j-30°).In the case of j ≤ 30°, the value of ψ equals to 0°.Exceptionally, in the example No. 4, the angle of dilatancy is assumed to be the difference between the Coulomb angle and the basic angle used in the hyperbolic description of the non-linear failure envelope [18][19][20].There are two commonly used definitions of the factor of safety in FEM, namely the strength reserving definition, and the overloading definition [21][22][23].The location of the critical failure surface, as well as the factor of safety, is dependent on the chosen definition.In this paper, the value of the factor of safety is determined based on the strength reduction method (SRM), where the soil strength is artificially weakened until the slope fails.Numerically, this occurs when it is no longer possible to obtain a converged solution.This can be expressed by the equation: ( where SRF is the total multiplier (strength reduction factor) that is used to define the value of soil strength parameters at a given stage of the analysis, c, j are input values of shear strength parameters, and c m , j m are mobilized or reduced values used in the analysis.At the start of a calculation, SRF is set to 1.0, i.e. all material strengths are set to their input values.In the state of failure, the SRF defined by eqn.(3) corresponds to the value of the factor of safety given in eqn.(1).In FEM, no assumption needs to be made about the shape or location of the failure surface.Failure occurs through the zones within the soil mass in which the soil shear strength is unable to resist the applied shear stress [15].
All calculations are performed in static drained conditions assuming effective shear strength and deformation parameters, without considering seismic effects.The influence of groundwater level and distributed load is examined.
Slope stability analyses using limit equilibrium and strength reduction methods Various investigations [14,24,25] showed that it cannot be readily concluded which method gives the higher value of the factor of safety, as the analysis is dependent on the specific problem and software package used in the analysis.In order to investigate the influence the implemented numerical algorithm exerts on the final result, certain analyses have been performed using four commercial software packages.BGSLOPE and Slide are based on the limit equilibrium method, whereas Plaxis and Phase 2 are based on the finite element method.For comparison purposes, the shapes of critical failure surfaces are overlaid on the Plaxis finite element model.

Example 1: Side slope of the Beli potok tunnel on the Belgrade Bypass
According to the design [26], the side slopes of the cut and cover part of the Beli potok tunnel are 19 m high, with 1/2 inclination (v/h), and one berm at mid-height of the slope 5 m in width.The ground profile is composed of six horizontal layers, as shown in Figure 1.

Figure 1. Typical cross-section of side slope
Based on the monitoring of piezometer constructions, the groundwater level is situated at the depth of 18 m.Triaxial CU tests results (with pore pressure measurements), simple shear and oedometer results, and extensive in situ (CPT) measurements, were used to derive shear strength and deformation parameters.The parameters used for the analyses are summarized in Table 1.
For layers #3 and #4, laboratory tests showed variation of the cohesion and angle of shear resistance values, and so the analyses were performed for two different sets of parameters.
The BGSLOPE software was used to compute the factor of safety by means of limit equilibrium methods.In order to obtain the critical failure surface, seven control points of the smooth Bezier curve were successively moved in various directions, in 0.05 m increments, until the smallest factor of safety was obtained.The smallest factor of safety was calculated by assuming the half-sine interslice force function.In order to eliminate tension in the upper part of the slope, the tension crack consistent with the factor of safety was introduced.The influence of tension will be discussed in more detail in example No. 3. The Plaxis software was used to compute the factor of safety with SRM.The finite element model consists of 3011 15-noded triangular elements automatically generated based on a robust triangulation procedure.An average element size was 1.168 m.The critical failure surfaces are compared and shown in Figure 2. The values of the factors of safety are summarized in Table 2.

Example 2: embankment widening on M-19 road
The side embankment, approx.8 m in height, is to be constructed for the widening of the existing road M-19 on Belgrade-Ljig section (approximately at KM 4+100).The design proposed several alternative solutions as elaborated in [27].Only the solution with the lightweight fly-ash material is analysed in this paper.The slope geometry is shown in Figure 4. Soil parameters used in the analyses are derived from field investigations and laboratory testing, cf.Table 3. Shear strength parameters for fly-ash were determined by direct shear testing on compacted samples (according to Proctor test), consolidated for 42 hours [28].It is important to note that fly-ash exhibits pozzolanic activity, i.e. it has the self-binding ability due to chemical reaction with calcium hydroxide in the presence of water.This leads to an increase in shear strength during compaction and consolidation processes.The linear Mohr-Coulomb envelope was defined for the normal stress level of up to 150 kN/m 2 .Young's modulus was determined by correlation with the previously determined oedometer modulus (on compacted samples) [28] according to the theory of elasticity with the Poisson's ratio assumed to be equal to 0.3.
The BGSLOPE software was used to perform LEM analyses.Linear distribution of the interslice force function is assumed along the optimized Bezier curve.This distribution is defined by specifying one parameter (z 1 = 0.1), on the x-axis along failure surface, whereas one additional value of parameter  Slope stability analyses using limit equilibrium and strength reduction methods z is assumed to be equal to 1, at the middle of the failure surface.The finite element model in Plaxis is defined by 1716 15-noded triangular elements.An average element size is 0.942 m.Critical failure surfaces overlap, as shown in Figure 5.The corresponding factors of safety are summarized in Table 4.

Table 4. Calculation results
As can be seen in Figure 5, the location of the critical failure surface obtained by Plaxis (j = ψ) corresponds to the LEM solution, whereas the difference between the factors of safety obtained by Plaxis (j ≠ ψ) and LEM is insignificant.The critical circular failure surface overestimates the factor of safety by 6.7 %.The influence of load distribution on the final result is analysed in the following two examples.The occurrence of tension on top of the slope (with horizontal surface) is also discussed.

Example no. 3: NTNU slope
This relatively simple slope was analysed as a part of PhD thesis defended at Norwegian University of Science and Technology.In study [17], the software package Slope/W with Morgenstern-Price method was used to perform the LEM analyses.Tension was not considered.
In order to inspect the findings, an example is re-analysed.
The model consists of a two-layered slope, as shown in Figure 6.Parameters used to perform the analyses are summarized in  In the limit equilibrium calculations, the tension occurred in the upper part of the slope due to cohesion term in the soil strength description of layer #1, as shown in Figure 6.The tension is indicated as irregularity in the line of thrust, which appears outside the sliding body, rendering the solution physically inadmissible.As a result of tension, normal forces at bases of the slices and interslice forces are negative, meaning that the slices are in a state of buoyancy although in reality there will be no tendency for the mass to lift upwards.It was also established that the presence of negative forces could lead to numerical instability  Zoran Berisavljević, Dušan Berisavljević, Vladimir Čebašek, Dragoslav Rakić when computing the factor of safety [1,11].In order to avoid this, Spencer [29,30] suggested that a tension crack zone should be introduced at the top of the slope.The depth of the tension zone was assumed to be equal to the depth of the zero active effective stress, consistent with the factor of safety of the slope: where r u is the pore pressure ratio, and H c is the depth to the zero active effective stress.
In order to make FEM results comparable with LEM values, a zero tensile strength is assigned to all materials in FEM [25,31].After a number of trials, it was concluded that the depth of the tension point zone is in good agreement with the depth calculated according to eqn. ( 4), for the ordinary combinations of strength parameters.The deformation parameters E and ν do not have any influence on the tension zone depth.The angle of dilatancy increases the factor of safety, thus decreasing the depth of the tension zone.These findings are valid for a fully developed failure mechanism.The location of the critical failure surface given in Plaxis (j = ψ) corresponds to the BGSLOPE solution, as shown in Figure 7b.

Table 6. Calculation results
Failure mechanisms obtained from Phase 2 and Slide are identical to the Plaxis and BGSLOPE calculations, respectively, and are thus not included in Figure 7.The corresponding factors of safety are summarized in Table 6.
The calculated tension crack depth of 0.92 m is in good agreement with the depth of the tension zone obtained by Plaxis, as shown in Figure 7a.If compared to the Plaxis (ψ = 0°) solution, the critical failure circle overestimates the value of the factor of safety by 5 % and only by 2.7 % if compared to the Phase 2 (ψ = 0°) results.
In order to investigate the influence of distributed load on the stability of the slope, an arbitrary chosen value of 92 kPa is assumed to act on the top of the slope, as shown in Figure 8.The effect of distributed load eliminates the tension.Calculation results are summarized in Table 6.The failure surface in Plaxis shows sharp transition in the upper part of the slope due to formation of the active Rankine zone.The smallest factor of safety is obtained in BGSLOPE by optimization of the polygonal failure surface defined by nine points.The value of Fs = 1.387 is by 4 % higher than the value obtained by Plaxis for an assumed non-dilatant material behaviour.The shape of the polygonal failure surface corresponds well to the Plaxis solution, as shown in Figure 8.The same critical failure surface is obtained from Slide (GLE method) with the corresponding factor of safety, Fs = 1.386, which is practically the same value as the one computed in BGSLOPE.Slope stability analyses using limit equilibrium and strength reduction methods

Example 4: New Valley Project -stability of rock slopes
This example examines stability of the rock slope formed for construction of the foundation pit for the pumping station in the scope of the "New Valley Project" in Egypt.Plaxis calculations performed for software verification [20]

Conclusions
The comparative study was performed in order to investigate the applicability of FEM and LEM to calculate the smallest factor of safety and location of the critical failure surface.The influence of different parameters on the final result was investigated.
In the presented examples, the critical circular failure surface overestimates the factor of safety by up to 9 %.Cheng [32] gives a comprehensive review about the problems associated with the possibility of locating the critical failure surface in the LEM.The fully automatic optimization could sometimes lead to the solution corresponding to a local minimum.For this reason, the semi-automatic search technique, such as an interactive optimization algorithm, is suggested for finding the location of the critical failure surface.The smooth Bezier curve is in most cases more suitable for calculating the smallest factor of safety than the polygonal curve.However, sharp transitions, as in the case of distributed load acting at the top of the slope, are more readily described by polygonal curve.The shape of the interslice force function f(x) should be chosen so as to obtain the smallest factor of safety.Satisfactory results are obtained by commonly used distributions like Spencer, half-sine, or linear.
It is commonly accepted that the values E and ν and the domain size do not have any (or have a very small) influence on the factor of safety and the location of the critical failure surface in SRM.If realistic values of E and ν (determined by laboratory or in situ tests) are used, the results should be in close agreement with LEM results.Cheng et al. [14] show that the greater differences between the two methods could be expected in a special case of the slope with a thin soft band.Despite the generally low confinement environment, the angle of dilatancy exerts some influence on the SRM results.If the associated flow rule is assumed in SRM, the locations of critical failure surfaces obtained by the two methods are virtually the same.However, the best prediction of the factor of safety value compared to LEM is obtained for the non-associated flow rule, in case realistic values of positive dilatancy are incorporated into analysis [33,34].If, for example, the non-dilatant material behaviour is assumed in the first case given in example No. 4, the value of the factor of safety is 1.796.In case of the associated flow rule, this value is by 6 % higher.Another parameter influencing the SRM results is the finite element size.It is suggested that the analysis should be performed with the local mesh refinement, so that the time required for calculation is acceptable from the practical point of view.The advantage is given to the 15-noded over the 6-noded triangular finite elements.If zero tensile stresses are assumed in the tension cut-off option in Plaxis, a good correspondence is observed between the tension crack depth given in eqn.( 4), and the distribution of plastic tension points.Different numerical algorithms implemented in software packages such as Plaxis, Phase 2 , Slide, and BGSLOPE, and the varying choice of the number of allowed iterations, tolerance, number of slices, etc. have an influence on the final result in both FEM and LEM.From the practical point of view, both methods provide similar results If correctly applied with all available software options.

Figure 2 .
Figure 2. Results of analyses As shown in Figure2, the optimized critical failure surface obtained by LEM corresponds well with the finite element solution in case of the associated flow rule (j = ψ).The same safety factors were obtained for Plaxis (ψ = 0°) and an optimized Bezier curve.The critical circular failure surface overestimates the value of the factor of safety by 4.3 %.Safety factor values obtained with the second set of parameters are shown in Table2, while shapes of failure surfaces are given in Figure3.The critical circular failure surface is located at the upper part of the slope, overestimating the factor of safety by 3.9-4.3%.The critical failure surface in Plaxis is obtained for ψ = 0° and amounts to Fs = 1.128.It extends throughout the entire height of the slope.The critical LEM failure surface amounts to Fs = 1.137, and its position corresponds to the Plaxis j = ψ solution.If the Bezier curve is optimized to have the same shape as the Plaxis ψ = 0°° solution, a local minimum is obtained.If the smallest factors of safety are compared by the two methods, it can be seen that their difference is insignificant and amounts to 0.8 %.The smallest factor of safety of the circular failure surface that extends from the crest to the toe of the slope is Fs = 1.206.

Figure 4 .
Figure 4. Road embankment widening using fly-ash material

Figure 5 .
Figure 5. Results of analyses: a) LEM; b) LEM surfaces overlaid on a finite element model

Figure 8 .
Figure 8. Influence of distributed load

Figure 7 .
Figure 7. Results of the analyses: a) distribution of plastic points; b) factor of safety values

Figure 10 .Figure 9 .
Figure 10.Results of analyses a) tension zone depth; b) shape and position of failure surfaces

Figure 11 .Figure 12 .
Figure 11.Analysis for varying distribution of load

Table 2 . Calculation results
*Values of shear strength parameters c i φ correspond to an analysis shown in Figure3

Table 1 . Input parameters for analyses Method Safety factor
(figure 2) a Only the smallest factors of safety are presented Zoran Berisavljević, Dušan Berisavljević, Vladimir Čebašek, Dragoslav Rakić

Table 5 .
The BGSLOPE (Maksimovic's method) and Slide (GLE method) software are used to compute the LEM factor of safety.The half-sine interslice force function is assumed in calculations.The Plaxis and Phase 2 are used to obtain the SRM factor of safety.The calculation model used in Plaxis consists of 2249 15-noded triangular finite elements.An average element size is 0.516 m.Phase 2 model is represented by 7433 6-noded triangular finite elements.

Table 7 .
The analyses were performed assuming dry slope conditions.FEM results are presented for the model with 3459 15-noded triangular finite elements, with an average element size of 1.740 m.LEM calculations were performed with the software package Slide (GLE method).The interslice force function f(x) was assumed to have linear distribution defined along the horizontal axis by two parameters (z 1 = 0.4 and z 2 = 2) at the opposite sides of the failure surface.The location of the critical failure surface was found, as in the preceding example, by an automatic search technique based on the Monte Carlo method.Calculation results are shown in Figure10, and safety factors values are summarized in Table8.The second analysis considered the case with distributed loads acting on the top of the slope.