Stability Analysis of Earth Slope Using Combined Numerical Analysis Method Based on DEM and LEM

Based on a 2D discrete element method (PFC2D code) and limit equilibrium method, a new method is proposed to obtain critical slip surface and to determine the factor of safety of the slope model. After the calibration process, the excavation process on the rectangular assembly of particles can make the slope model of the DEM. As the trial acceleration of gravity increases, the critical slip surface tends to change from shallow to deep, which indicates more serious failures in the slope. With higher trial accelerations of gravity, the factor of safety calculated by this method continually increases as well. Therefore, the critical trial acceleration of gravity that places the slope in a critical state of failure is essential to ensure the results are correct. The gravity increase method is applied to render the critical state of failure for slope model. As soon as the continuous failure surface forms, the influence parameters of particles on the critical slip surface also can be determined. Connecting the centre of these particles forms the critical slip surface, which is more consistent with the actual situation in nature due to the absence of assumptions of the critical slip surface in this method. Then, the factor of safety is calculated by the Spencer method that has the ability of dealing with noncircular, complex slip surface. This method proves to be a promising tool to analyse slope stability with the DEM. The character of progressive failure process of slope can be well modelled by DEM by investigating the development of failure in the slope body. By studying temporal and spatial risk of slope instability by the total failure process, it is possible for managers and engineers to detect temporal and spatial risks so robust suggestions and mitigation measures can be made.


INTRODUCTION
In the field of slope engineering, limit equilibrium method (LEM) is, currently, the most popular method due to its simplicity and applicability.However, it is limited to overly simplistic problems because of adopting many assumptions, such as circular critical slip surface, and providing little insight into slope failure mechanisms [17,4].As a matter of fact, slope failure involves complex mechanisms (e.g., nonlinear constitutive models, progressive failure and strain-softening behaviour) [7,8,19,13,22,23,24,27]. LEM is inadequate when addressing these challenges.
Numerical modelling techniques based on continuum mechanics were initially proposed to conduct research related to the stress fields and displacement fields in the slope [6,25,28].Then, the advance of strength reduction method and the gravity increase method makes it possible to find the critical slip surface and to compute the factor of safety (FS) either by reducing shear strength parameters of the slope or by increasing the gravity loading respectively [18].The numerical results-normally the finite element method (FEM) with the strength reduction method is adopted, are close to the limit equilibrium method, which has been proved by the literature [9,12,29].It requires fewer prior assumptions in searching for the critical slip surface, but the model is restricted to oversimplified stressstrain constitutive relationship.
A tool that appears to be promising for modelling the behaviour of soil and rock mass, particularly large-scale damage and non-linear behaviours [14], is the discrete element method (DEM) originally developed by [5].Unlike FEM, Particles Flow Code (PFC) can simulate the physical micro-mechanisms directly and complex constitutive relation of numerical model is bypassed.The advantages of this code, a well-known discrete element method, are that the deformation and failure of objects stem from the interaction between the micro-particles; therefore, even block breakage and separation can be considered [11].Significant developments in slope stability analysis have been led by the use of this code for numerical analysis.Wang et al. [21] first created a particles assembly to represent heavily jointed rock slopes and analyzed the impact of joint persistence on slope stability.Tang et al. [20] investigated the kinematics and mechanical behaviour of the Tsaoling landslide with the DEM.Scholtès and Donzé [16] demonstrated the application of this technique in studying the step-path failure of joint slope stability involving a complex interaction between pre-existing discontinuities and intact rock bridges.Camones et al. [2] studied stability of a rock slope using the gravity increase method with PFC, in which crack propagation and coalescence were observed.
Because of the repeating calibration process of the DEM for its properties, unlike the FEM, it may not be convenient and direct to apply the strength reduction method (SRM) to DEM.Further, the factor of safety obtained by gravity increase method (GIM) is not precise in some situations [30].The objective of this paper is to develop a new method combining concepts of the limit equilibrium method and the DEM to determine the critical slip surface and to calculate factor of safety of the slope.The slope model of the DEM can automatically form the critical slip surface from the calculated results at the stage of the incipient slope failure.The factor of safety can be further obtained by LEM--the Spencer method in this study.The DEM of the slope also shows its capability of studying progressive slope failure processes, which can improve our knowledge on total slope failure processes.Finally, the proposed method is verified by two slopes.It is shown that numerical results agree well with the laboratory test, numerical results of classical method -the LEM method and strength reduction method.

DEM FOR THE SLOPE STABILITY ANALYSIS 2.1 Bonded-Particle Model
In this study, a discrete granular simulation technique (PFC 2D) is adopted, and an assembly of bonded particles is used to represent the solid (e.g., rock and soil).There are two basic bond models: contact bond model and parallel bond model.The contact bond model can be treated as the simplified parallel bond model because the contact bond between two bonded particles behaves as parallel bond of radius zero.While parallel bond resists a bending moment or opposes rolling, contact bond can only resist force at the contact point.The parallel bond model was applied to all the PFC models in this paper.
For the parallel bonded-particle model, the total force and moment acting at each bonded contact consist of the particle-based force and the force and moment carried by the cement.To the basic particle-based contact, it is a linear elastic relation between the particle-based force and particle-particle overlap, given by n s s s where F n is the normal contact force, K n is the contact normal stiffness, U n is the contact overlap, increment of shear force is ΔF s , contact shear stiffness is K s and shear displacement is U s .The cement-based portion of the force and moment is computed in an incremental fashion.When the parallel bond is formed, the force and moment are initialized to zero.In such case, the elastic force and moment is written by: When the normal tensile stress or the shear stress of the cement surpass its respective parallel bond strength, the bond is broken and removed, by which cracks are produced in the model.The micro-cracks initiation, propagation and coalescence are automatically determined without further process (e.g., re-meshing).More details on explanations of parallel bonded-particle model can be founded in the literature [14].
Unlike other numerical methods, the above input micro-properties of PFC are usually unknown from experiments in laboratory and field investigations.To obtain these parameters, the calibration process-the socalled trial-and-error process-is required for the model.If a particular set of model micro-parameters can reproduce the desired physical material, these micro-parameters are treated as effective input parameters and can be utilized to perform further research.This calibration process consists mainly of comparisons between the numerical simulation of tests and corresponding physical tests (e.g., uniaxial compressive test and direct shear test).

Formation of the Slope Model
In this study, an example is used to demonstrate the procedure of this method.The mechanical properties of the slope model are listed in Tab. 1.The rectangular assembly of particles should be first generated prior to the formation of the slope [21].The particles were arbitrarily generated in a square domain confined with walls, and minimum radii and radii ratios (Rmax/Rmin) of the particles mainly control the particle size distribution.Micro-properties (i.e., contact modulus, particle normal/shear stiffness and friction coefficient) are required for the deformability of the particles.After the specified assembly of particles was generated, parallel bonds were installed throughout the assembly between all of the particles to reproduce the behaviour of soil and rock mass.For the parallel bond model, the micro-properties that define its deformability and strength are the parallel bond modulus, the ratio of parallel bond normal stiffness and shear stiffness and the parallel bond strengths.The uniaxial compression strength (UCS) of samples is 75 kPa.
To obtain micro-properties of the model described above, a set of biaxial compression tests was applied to calibration process as shown in Fig. 1.In simulating a series of biaxial compression test, top and bottom walls as load face were at the velocity, v p = 0.01 m/s.To maintain specify confining pressure (σ 3 ), lateral walls' velocity was decided by servo control function, which is described in detail by Itasca [10] and we will not repeat it here.After the calibration process, micro-properties used in the model were obtained and shown in Tab. 2 and the associated stress-strain behaviour was presented in Fig. 1(b).
The maximum vertical stress at the bottom should be 0.18 MPa because the unit weight of the slope is 19.6 kN/m 3 and its depth is 9 m.However, the inherent porosity of the packed particles decreases the overall density of the entity.To compensate for this decrease, another calibration process is required to reach the objective stress at the bottom, which is performed by simply increasing the density of particles [21].
After that process, a slope model could be created by deleting specified particles above the slope face, which can be achieved by the customized FISH function.This procedure, which models the actual excavation process in the field, was divided into five steps in this case (Fig. 2) to minimize the dynamic impact of deleting the particles [21].During the simulation, the boundary is fixed using walls in the DEM model.Therefore, the displacements of particles are not allowed in the normal directions of the left, right and bottom boundaries.

COMPUTATION OF FACTOR OF SAFETY 3.1 Gravity Increase Method
There are two effective ways for numerical modelling to render the critical state of failure for the slope when conducting slope stability analysis: the strength reduction method (SRM) and the gravity increase method (GIM) [18].The gravity increase method was adopted in this study.The calibration process of the DEM to obtain appropriate micro-properties is essential to the reliable model results but it is time-consuming and monotonous as well.The strength reduction method is performed by reducing the rock or soil shear strength (i.e., internal friction angle and cohesion) in stages until the slope fails; therefore, a series of shear strength properties are required to the slope model.Consequently, it seems impossible to implement strength reduction method to PFC.
Based on the gravity increase method, we only need to take the trial acceleration of gravity into consideration while holding the material properties constant.Highly increased gravitational acceleration was adopted to determine the collapse of slopes, and the critical gravitational acceleration of the model needs to be carefully calibrated until the continuous failure surface of the slope can be formed.The bisection method was applied in this process to find critical gravitational acceleration that iteratively narrows the interval until sufficient accuracy of the solution is met.
The factor of safety can be defined by the ratio of the gravitational acceleration as follows, according to gravity increase method: where g 0 is the acceleration of gravity in the initial state (9.81 m/s 2 ) and g trial is the acceleration of gravity at failure (m/s 2 ).
However, in this method, the factor of safety is not determined by the formula used in the gravity increase method but by the concepts of the limit equilibrium method.The reason is that the factor of safety in the limit equilibrium method, which is defined as the ratio of the sum of the resisting forces divided by the sum of the driving forces, is more reasonable and acceptable.Further, in some extreme cases, the g trial has to be much larger to render the slope model of the DEM unstable [25], such that the factor of safety by the gravity increase method may be unreliable.
The fixed maximum horizontal displacement for the slide mass that can indicate the formation of continuous failure surfaces is treated as the criteria for slope failure.By contrast, if the slope does not have the obvious deformation, we adopted average unbalanced contact force acting on all particles to diagnose the state of the model.Theoretically, the average unbalanced contact force is almost zero at equilibrium.Taking the magnitude of the average contact forces of assembly of particles as reference, the specific ratio of average unbalanced force divided by the average force over all particles is treated as the criteria for the determination of slope stability.If the ratio is less than a given value, the slope is considered stable and we would stop the calculation.In this paper, the given value of the ratio is 0.01.As a matter of fact, if the slope is unstable and calculation steps are enough, the initiation of failure and the large-scale movement of the slide mass in the transportation stage during the slope failure process are clearly observed.

Location of the Critical Slip Surface of Slope
Finding the critical slip surface within the slope body continues to be a controversial and unsolved issue.For the DEM, the continuous failure surface of the slope in the initiation failure (the critical slip surface in this study) can automatically form without any assumption.
For this slope model, once the critical g trial of 15.9 m/s 2 was reached, the continuous failure surface was automatically formed, as shown in Fig. 3(a).Tensile cracks in the top of the slope were just connected at this stage.After the initiation failure stage, the transportation and comminution stage of slope failure can be seen in Fig. 3(b) and the slip surface tends to be smooth due to the slipping of the slide mass.The critical slip surface of the slope model is visualized by deleting the particles of the slide mass (Fig. 4).The slide mass of the slope is defined by the horizontal displacement in this study; if the horizontal displacement of a particle is greater than a certain value (e.g., 0.01 m in this case), it can be classified into the slide body.Technical Gazette 25, 5(2018), 1265-1273 The exact location of the critical slip surface can be obtained automatically from the coordinates of the particles right above the failure surface (black particles in Fig. 4).Data of these particles were output from the assembly of particles, and then, the entire curve of the failure surface would be obtained by connecting the centre of these particles, as shown in Fig. 5.The more accurate location of the critical slip surface is dependent on the smaller radius of the particles in the slope model.However, this requires the generation of more particles in the model, which will highly increase the computational burden.The simplified Bishop method [1], a classical method for slope engineering, was also used in this study for comparison, as shown in Fig. 5. From the comparison of the critical slip surfaces, their positions are very close.Furthermore, the surface found in this study is not a circular arc, and the tension fracture at the top of the slope can be indicated.

Determination of the Factor of Safety
The Spencer method was used to calculate the factor of safety for the above critical slip surface.It is one type of the limit equilibrium method, commonly adopted by many engineers to compute factor of safety due to its widely acceptable results and its ability to deal with noncircular, complex slip surfaces.
Acting forces on a typical slice are shown in Fig. 6.This method assumes that the ratio λ between the vertical inter-slice forces and horizontal inter-slice forces is constant.θ indicates the angular of the resultant forces, i.e., tan θ = λ; W i -the weight of slice i; α i -the inclination angle of slip surface at the bottom of slice i; b i -the horizontal width of slice i; length of slice, l = b i secα i .
Based on force and moment equilibrium, Eq. ( 5) can be obtained.By using the Newton-Raphson method, the factor of safety, FS, and λ are determined as follows: ( , ) 0 ( , ) 0 the first equation with E n is referred to force balance while the moment balance is related to the second equation of M n .FS and λ are solved using an iterative algorithm.For further details, readers are referred to the reference [3].The factor of safety calculated by this method is 1.706, and other methods (i.e., GIM, the simplified Bishop method, the simplified Janbu method and SRM by FLAC) were also utilized in this section for comparison, as listed in Tab. 3. It is necessary to note that all limit equilibrium analyses are conducted using the code, Slide [15].In terms of the calculated results, the factors of safety found using different methods agree with each other very well, which proves that the results calculated by this method are effective and acceptable.

EFFECT OF TRIAL ACCELERATION OF GRAVITY
To render the slope model critical state of slope failure in this study, the trial acceleration of gravity g trial was loaded according to the gravity increase method.Therefore, g trial , as an important factor for slope failure, was studied in this section.
The g trial varied from 15.9 to 17.0 (i.e., g trial = 15.9, 16.0, 16.5 and 17.0, respectively) in this study.Obviously, when g trial was less than the critical g trial (15.9), the deformation of the slope could be observed, but the slope was stable.From computation results of the same slope model with different trial accelerations of gravity, the Fig. of slope failure and the critical slip surface are shown in Fig. 7, and the relationship between g trial and the factor of safety calculated by this method is listed in Tab. 4. When g trial was close to the critical g trial , 15.9, the progressive failure process was actually the same, and it is difficult to distinguish them from the critical slip surface without exact measure.As the g trial increases, the critical slip surface tends to change from shallow to deep in the slope, which means more serious failure happens in the slope by higher g trial , as shown in Fig. 8.Moreover, a secondary failure surface was formed in the inner slope body when the g trial reaches 17.0, even though the continuous failure surface was first formed (Fig. 7(d)).

PROGRESSIVE FAILURE ANALYSIS
The progressive failure of the above slope with the critical g trial , 15.9, was also studied (Fig. 9).The  The simulated process agrees well with the actual failure process usually developed in the slopes.At the first stage of slope failure, the failure cannot be observed solely from the shape of slope (Fig. 9(a)).However, from the Fig.
of the bonds and the cracks, the failure in the toe region is formed initially, and cracks also scatter in the slope body, even in the bottom due to the heterogeneity of the DEM slope model.As the time step increases, failure in the slope progresses up towards the top, and cracks gradually propagate to there as well.When the time step reaches 110 thousand steps, the crest of the slope fractures due to coalesced tensile cracks and, right at this time, the continuous failure surface is formed, thus demonstrating that the critical slip surface can be determined using this method.After the transportation stage of the slope failure process, the slope remains stable again when the time step reaches 820 thousand steps.
With the DEM, the progressive failure process of the slope can be simulated, and most importantly, the critical failure surface can be generated automatically without assumptions for the location, shape or optimization methods for determining the critical failure surface.Furthermore, this method can produce the total slope failure, which can be used to study temporal and spatial risk associated with massive slope failure [17].
Consequently, robust recommendations and mitigation measures can be studied and made.

VERFICATION AND EXAMPLICATION
To demonstrate the capacity of the method, two slope examples were used for analysis.Example 1.The example originated from a model experiment [26].After the calibration process with biaxial compression tests, the micro-parameters were determined (Tab.5).The minimal radius of the particles in this slope model was 0.8 mm, in which 18,506 particles were generated.Fig. 10 shows the DEM of 24 cm high and 85° slope with the unit weight of 25 kN/m 3 , Young's modulus of 200 MPa, Poisson's ratio of 0.35, an internal friction angle of 32° and a cohesion of 40 KPa.
The calculated result of the DEM is shown in Fig. 11(a).The continuous failure surface is obvious, which is similar to the calculated result from displacement vector by the FEM, as shown in Fig. 11(b).Further, the progressive failure is the same, from the toe region to the top where tensile cracks can be observed.Finally, the factor of safety calculated by this method was compared with those found using other methods and the experiment [26] as shown in Tab. 6.In the experiment, as soon as g trial reached 45g 0 , the critical failure surface was formed, and the slope failure instantaneously happened.Therefore, the factor of safety is 45, which is very similar to the FS of 43.9 determined by this study, and the difference is small when compared with other methods.
Example 2.   The weak layer was made among the assembly of the particles by specified FISH functions before creating the slope model.During the excavation process for building slope model, more calculation steps were needed to decrease the unbalanced force among the particles.The aim is to protect the intactness of the bonds from the unloading effect of the excavation process because the strength of this slope is relatively weak.Fig. 13(a) shows the slope failure of the DEM after sufficient computational steps.Shear failure of the slope expanded along the weak layer and was associated with tension fractures in the upper slope, which is exactly described by the distribution of the failure zone in FLAC3D with strength reduction method (Fig. 14(c), where block state-of-none refers to the elastic state, and block state-of-shear-n, tension-n, as well as shear-n tension-n, represents the shear failure now, the tension failure now and mixed shear and tension failure now, respectively).Although the numerical modelling of the continuum methods with appropriate constitutive model, such as FLAC3D in this example, can also study the failure mechanism of the slope at the incipient stage, it cannot model the actual movement of the slide mass and the postfailure behaviour of slope, which is the unique advantage of discrete element method.For the same slope, the slope stability was analysed using the simplified Bishop method, Spencer method, as shown in Fig. 14.The critical slip surface was obtained from the DEM when a continuous failure surface was automatically formed, as shown in Fig. 13 (b).The position of the critical slip surface of the other classical methods is similar to this method.The factor of safety of the slope computed by the proposed method is in good agreement with the other methods, as shown in Fig. 13 and Fig. 14.This method does not require assumptions on interstice forces and experiences on locations of critical slip surface.The critical slip surface and slide mass are automatically formed as the parallel bonds progressive breakage, which means that this method is more direct approach to investigate the slope stability and it has the potential in further development in rock slope of complex geologic structures involving bedding, foliations, joints etc.Unlike finite element methods, it can also study the failure mechanism of the slope at the incipient stage to the postfailure behaviour of slopes.

CONCLUSION
The combination of DEM and LEM is a promising method to calculate the factor of safety for slope stability analysis because it has the advantage of automatically locating the critical slip surface.This method expands the application of the discrete element method into slope stability analysis to determine the factor of safety for formed critical slip surfaces by adopting the limit equilibrium method.Its calculated results were verified by other well accepted methods, which proves this method is effective.
The slope model of the DEM (particle flow code) can be made using the excavation process on a rectangular assembly of particles.The gravity increase method is utilized to render the slope model critical state of slope failure.As soon as the continuous failure surface is formed, the critical slip surface of the slope also can be automatically determined and output, with no assumptions for the critical slip surface in this process.Then, the factor of safety of this critical slip surface for the slope is calculated by the Spencer method, which has the capability of dealing with noncircular, complex slip surfaces.
As the trial acceleration of gravity increases, the critical slip surface tends to change from shallow to deep, which indicates more serious failures in the slope.With higher trial accelerations of gravity, the factor of safety calculated by this method continually increases as well.Therefore, the critical g trial that places the slope in a critical state of failure is essential to ensure the results are correct.
The character of progressive failure process of slope can be well modelled by DEM by investigating the development of failure in the slope body.By studying temporal and spatial risk of slope instability by the total failure process, it is possible for managers and engineers to detect temporal and spatial risks (e.g., mapping hazard zone) so robust suggestions and mitigation measures can be made.

M
are the normal-and shear-directed forces and moments of the cement, n k and s k are the are normal and shear bond stiffness per unit area, A, J, I are the area, polar moment of inertia and moment of inertia and θ n and θ s are the components of rotation angle.

Figure 1
Figure 1 Numerical biaxial compression test: (a) Test sample with principal stresses exerted, (b) Stress-strain relationships.The uniaxial compression strength (UCS) of samples is 75 kPa.

Figure 2
Figure 2 Geometry of the slope and excavation by 5 steps

Figure 3
Horizontal displacement diagram of slope failure.(Unit: m)

Figure 4 Figure 5
Figure 4 Critical slip surface of the slope

Figure 6
Figure 6 Force diagram of slice for Spencer method.

Figure 7 Figure 8
Figure 7 Slope failure modes under different trial accelerations of gravity.Horizontal displacement under the g trial 15.9 (a), 16.0 (b), 16.5 (c), 17.0 (d) Fig. of the parallel bonds can indicate the stability of the slope as the assembly of scattered particles is bonded together to have strength (Fig. 9(b)).The failure process in the slope, therefore, can be seen directly from the Fig. of the parallel bonds.Each bond-break event in the slope was recorded as a segment linked to the centre of the parent particles to represent the crack in the slope body, as shown in Fig. 9(c).

Figure 9
Figure 9 Progressive process of above slope model: (a) Slope model of the DEM (b) Visualization of the parallel bonds (c) Distribution of cracks in the model

Figure 10
Figure 10 Schematic view of the DEM slope model in Example 1 Fig. 12 shows a slope model that includes two types of soil, in which the black particles indicate a weak intercalated layer with the thickness being 0.5 m in the slope.The macro-parameters are shown in Tab. 7, and Tab. 8 displays the main micro-parameters.The calibration process for the weak layer was performed by uniaxial compression and direct shear tests.

Figure 11
Figure 11 Calculated results of displacement in Example 1: (a) The horizontal displacement of slope failure (b) Displacement vector of FEM (Li 2009)

Figure 12
Figure 12 DEM in Example 2 (a) Horizontal displacement of PFC in Example 2 (b) Factor of safety of this study Figure 13 Calculation result of this method (a) Simplified Bishop method (b) Spencer method (c) Distribution of failure zone by SRM Figure 14 Results of the classical methods in Example 2

Table 1
Properties used for the simulated slope

Table 2
Input micro-properties of the model by the results of a series of biaxial mechanical tests

Table 3
Calculated results of the factor of safety

Table 4
Relationship between the trial acceleration of gravity and the factor of safety

Table 5
Main micro-parameter of the DEM in Example 1

Table 6
Comparison of factors of safety from numerical and experimental results

Table 7
Parameters for soils in Example 2

Table 8
Main micro-parameter of the DEM in Example 2