Stability Analysis of Three-Dimensional Slopes Considering the Earthquake Force Direction

Simplified 3D Spencer’smethod is proposed to evaluate the stability of slopes under earthquake force and a random search algorithm is used to determine the critical slip surface. A computer code was written in Matlab to determine the critical slip surface and calculate its safety factor and direction of sliding (DOS). The applicability of the presented code is verified by reanalyzing three slope stability problems. Then, a comparison is made to investigate the effects of the direction of earthquake force on the factor of safety and the DOS of the critical slip surface. Accounting for the effects of the direction of earthquake force, the difference in safety factors, and the DOS can be up to 40% and 40, respectively. Finally, a parametric study is also conducted to determine the effects of the soil strength parameters (c, φ). It is found that the effective cohesion c and friction angle φ of the soil have nearly no influence on the unique DOS of the critical surface. Application of the presented approach into 3D slopes under quasistatic load is straightforward.


Introduction
Estimating the slope stability remains one of the most important problems for geotechnical engineers.Slope stability analyses have usually been performed using two-dimensional (2D) simulation [1][2][3][4].However, in most cases, the geometry of the slopes has three-dimensional (3D) characteristics, such as corners, conical heaps, and dams in narrow valley.Traditional 2D methods do not consider 3D effects resulting conservatism in the analysis of slope stability under complex condition [5].
Recently, various 3D slope stability methods have been proposed to obtain more reasonable solutions.These 3D analysis methods can be divided into three categories: limit equilibrium method, limit analysis method, and finite element method.The limit equilibrium method is most widely used in practice to analyze the stability of slopes.A number of 3D limit equilibrium methods, based on simple extension of convention 2D methods, have been proposed to estimate the 3D effects on the safety of slope [6][7][8].All these methods assumed the sliding body to move along the plane of symmetry.Such an assumption, however, can lead to miscalculation of the factor of safety for slopes under asymmetrical conditionals, such as the complex geometry, variable soil stratigraphy, external loading, seismic forces, or other factors [9,10].
Huang and Tsai [9] introduced the DOS that defines the direction of movement of the sliding mass in the horizontal plane and proposed a new 3D Bishop's method that satisfies two-direction moment equilibrium.Based on this new concept of 'two-direction equilibrium' , a number of researchers proposed various 3D asymmetrical slope stabilities analysis methods [10,[13][14][15].Among these methods, simplified Spencer's method proposed by Wan et al. [15] satisfies all the force equilibrium and two-direction moment equilibrium.This method has relatively rigorous theoretical basis.The method presented by Wan et al. [15] is further exploited to investigate the stability of 3D slopes under quasistatic load.
Slope failure induced by earthquake loads can lead to serious damage and loss of lives.For example, the Wenchuan earthquake (Ms=8.0) on May 12th, 2008, directly caused more than 15000 landslides which resulted in more than 20000 deaths.Since the mid-1970s, increasing attention has been directed toward the stability analysis of slopes under earthquake forces [17][18][19].At present, there are four methods for seismic stability analysis of slopes: pseudo-static method, Newmark sliding block analysis method, numerical analysis method, and testing method.The conventional pseudo-static method is still widely applied for evaluating the effects of earthquake forces on the stability of a slope [20,21].Under three-dimensional conditions, earthquake forces can act in all three directions-along two horizontal directions (x and y) and the vertical direction (z) as shown in Figure 1.Most of the pseudostatic analyses have involved the consideration of the horizontal earthquake force Q h, and vertical earthquake force Q v .However, previous studies assumed that the horizontal earthquake force Q h is parallel to the x-axis and neglected the effects of the direction of earthquake force  e on the slope safety.To improve the accuracy of 3D slope stability analysis, the effects of the direction of the earthquake force on 3D slope stability analysis should be taken into account.
In this paper, the method presented by Wan et al. [15] is further exploited to evaluate the stability of 3D slopes under quasistatic load.An efficient optimization method proposed by Chen [22] is attempted to determine the critical slip surface of slopes.The effect of the direction of the earthquake force on the factor of safety and the DOS of the critical slip surface is also investigated.The difference in factor of safety and the DOS is compared to investigate the effects of the direction of earthquake force on the safety of 3D slopes.A parametric study is also conducted to show the effects of the soil strength parameters (effective cohesion   , effective friction angle   ) on the DOS of the critical slip surface.

Formulation of 3D Spencer's Method
As with other 3D limit equilibrium methods, the potential failure mass of a slope is discretized into a number of columns with vertical interface.Figure 2 illustrates the internal and external forces acting on the various faces of the column ith.To establish the force and moment equilibrium equations, the following assumptions are included: (1) The conventional definition for factor of safety F reduces the available shear strength parameters   and   by the following equations to bring the slope to a limiting state.
where   and   are the effective cohesion and friction angle of the soil, respectively;    and    are soil strength parameters necessary to maintain the structure in limit equilibrium, respectively.
(2) The horizontal shear forces, H xi and H zi , are assumed to 0. Such an assumption was also adopted by Hungr [6], Huang and Tsai [9] and Cheng and Yip [10].
The following relationships between the intercolumn vertical shear forces (X xi and X yi ) and the normal force (E xi and E yi ) of ith column are assumed: where  x and  y are shear force factors in the x-and ydirections, respectively.
(3) Soil columns are assumed to move in the same direction, i.e. the DOS is unique on the x-y plane for all columns.The unique DOS, a, is equal to the angle between positive x-axis and the projection of the shear force T i on the base of each soil column in xy plane (measured counterclockwise from the positive x-axis).
(4) The pseudostatic method is adopted by the writers in the present formulation and by many other researchers [20,21].The horizontal earthquake force, Q hi , is assumed to act towards the center of the column.The components of the horizontal earthquake forces, Q hxi and Q hyi , can be expressed as where K h is the pseudostatic acceleration coefficient and  e is the angle between positive x-axis and the earthquake force Q hi (measured counterclockwise from the positive x-axis).
The unit vector, (m xi , m yi , and m zi ), for the shear force T i of the ith column can be determined by where (n xi , n yi , n zi ) is unit vector for the normal force N i of the ith column.
Establish the vertical and horizontal force equilibrium for the ith column in x, y, and z-directions, as By substituting (3), ( 4), ( 9), (10), and Mohr-Coulomb's failure criterion 12 into (8), the effective normal force, N i, , can be expressed as where A i is the area of the column base.
The shear force T i can be determined by application of Mohr-Coulomb's failure criterion.
Considering the overall force equilibrium in the x-and y-directions Establishing the overall moment equilibrium equation about the x-and y-axis, respectively, where x i , y i , and z i are coordinate values of the center of a column base; z ei is the z-coordinate of the center of a column.The unknowns, namely, F,  x ,  z , and a, can be obtained by solving the system of equations consisting of ( 13), ( 14), (15), and (16).
The presented Spencer's method can be simplified to Bishop's method by considering only moment equilibrium equations ( 15) and ( 16) about axes parallel to x-and y-axis direction, respectively.Similarly, presented Spencer's method can be simplified to Janbu's simplified method by considering only force equilibrium equations ( 13) and ( 14) in x-and yaxis direction, respectively.One more assumption that the intercolumn vertical shear force (X xi and X yi ) is neglected is made in Bishop's and Janbu's simplified method.The normal force N i can be expressed as

Optimization Procedure
The studies dealing with slope stability analysis within the framework of limit equilibrium involves two sequential steps: calculating the factor of safety for a given slip surface and determining the critical slip surface and its corresponding minimum factor of safety.Factor of safety is calculated by the proposed limit equilibrium method.Optimization process is used to determining the shape and location of the critical slip surface.The slip surface is commonly assumed to be a particular shape, such as spherical, ellipsoidal, and nonuniform rational B-splines.This paper introduces a general ellipsoidal shape as the shape of failure in 3D analysis.This ellipsoid can rotate on the x-y plane and the rotated ellipsoid is presented as equation.
[cos  * ( −   ) − sin  * ( −   )] where  is the rotation angle of the ellipsoid in the x-y plane; X c , Y c , and Z c , are coordinates of the center of the ellipsoid in the x-, y-, and z-directions, respectively; R x , R y , and R z are semiradiuses of the ellipsoid in the x-, y-, and z -directions, respectively.An optimization technique proposed by Chen [22] that was found to be effective in previous work by Gao et al. [23] is used to find the minimum factor of safety.In this technique, a random search algorithm is used to find an initial estimate slip surface as the starting point in search of the global minimum, and then a minimization procedure is carried out until the bandwidths of the search variables become less than predefined values.The procedure is proposed to consist of the following steps.
Step 2. Generate a new slip surface Z 1 and its factor of safety F 1 .The coordinates of Z 1 are determined by where Z 1 is a new slip surface; l  is the bandwidths of the search variables; u(0,1) is a random number ranging between -1 and 1.
Step 4. Repeat Steps 2-3 and Z 0 and F 0 are renewed until the number of trials reached a specified value N, then l  = l  /2.

Results and Discussion
Example 1. Example 1 is a homogeneous, laterally symmetrical slope.The geometry and material properties of Example 1 are shown in Figure 3.This example was analyzed by Zhang [16] using Bishop's, Janbu's simplified and Spencer's method.
In these methods, a special ellipsoid slip surface (R y =R x , =0) was employed to determine the critical slip surface and its corresponding minimum factor of safety.The same example was performed to verify the performance of random search method in determining the minimum factor of safety.
Based on Bishop's, Janbu's simplified, and Spencer's method, random search method was applied to determine the critical slip surface and the corresponding minimum factor of safety.The results obtained by the presented method, compared to those provided by Zhang [16] using Bishop's, Janbu's simplified, and Spencer's method are listed in Tables 1, 2, and 3, respectively.The calculated DOS are equal to 0 ∘ , which is the same as expect.The critical slip surface obtained by the presented method and those provided by Zhang [16] are illustrated in Figure 4.It can be seen from Tables 1, 2, and 3 that the minimum factors of safety obtained by presented method is slightly smaller than those provided by Zhang [16].
Example 2. Example 2 is performed to validate the performance of the random research in determining the minimum factor of safety and its corresponding critical slip surface.This example was first analyzed by Alkasawneh et al. [11] to study the effects of different search techniques on the minimum factor of safety and was reanalyzed by Kalatehjari et al. [12] to verify the application of particle swarm optimization (PSO).Figure 5 illustrates the geometry, strength parameters, and unit weight of Example 2. For validation, a comparison of the minimum safety factors is summarized in Table 4.The critical slip surface obtained by presented method is illustrated in Figure 6.The results obtained by the presented study are close to those of Alkasawneh et al. [11] and Kalatehjari et al. [12].The difference between the result of this study and that of Kalatehjari et al. [12] is 1.6%.Moreover, the difference between the results of current study and the results proposed by Alkasawneh et al. [11] for 2D and 3D analysis are 0.18% and        0.06%, respectively.The calculated unique DOS is equal to 0 ∘ , which is the same as Kalatehjari et al. [12].
Example 3. To further investigate the ability of the present method to determining the 3D shape and unique DOS of the critical slip surface, a hypothetical 3D slope problem is designed and analyzed.The geometry of Example 3 and its geomechanical properties are shown in Figure 7.This example involves four different 3D models, as shown in Figure 8.The unique DOS of the four cases can be predicted because the boundaries of slopes were constant and the face was rotated 10 ∘ clockwise in x-y plane in each step of the example.The minimum factors of safety should be equal to each other for all the models.Table 5 shows the minimum factors of safety and the unique DOS for all conditions.Figure 9 shows the critical slip surfaces for Example 3. The maximum difference between the unique DOS obtained by the present method and those logically estimated is less than 0.2 ∘ , which demonstrates the ability of the proposed method to determining the unique DOS of the critical slip surface.
The minimum safety factors of the four models are not exactly equal to each other, but the maximum difference is less than 0.01%.This difference may be attributed to the discretization of the sliding mass and the different grid widths.This problem can be minimized by sufficiently small grid width.These results demonstrate the ability of the presented method to determining the 3D shape, minimum factor of safety, and unique DOS of the critical slip surface.
Example 4. Example 4 is selected to investigate the effects of the earthquake force and its direction on the minimum factor of safety and the unique DOS of the critical slip surface.The geometry of this example is symmetry about   F (2D) F (3D) Alkasawneh et al. [11] 1.70 (PLAXIS) 1.80 (PLAXIS) Kalatehjari et al. [12] 1.77 (Bishop) current(Ry=50) 1.697 (Spencer) 1.799 (Spencer)  the neutral plane, but the earthquake force acting on this slope is asymmetry as shown in Figure 10.The pseudo-static acceleration coefficient K h is set equal to 0.1, 0.2, and 0.3, and the direction of the earthquake force  e ranges from 0 ∘ to 90 ∘ .
Figure 11 presents the difference in the factor of safety F/F 0 (where F is the factor of safety with variable  e , and F 0 is the factor of safety with  e =0) between considering and neglecting the direction of earthquake force.It can be seen from ( 5) that the components of seismic force in the direction of x-axis Q hxi decrease the increase of  e .A smaller Q hxi generates a failure surface with a larger factor of safety.Thus, the factor of safety increases with the increase of  e , as is shown in Figure 11.It can be seen that the difference in factor of safety (F/F 0 ) slightly increases as the direction of earthquake force increases when  e ≤30 ∘ and then tends to increase obviously.The direction of earthquake force ( e ) has significant effects on the difference in factor of safety.Neglecting the effects of the direction of earthquake force on the factor of safety can lead to overestimation in 3D slope stability analysis.Typically, the difference can reach a maximum value of 40.5% when the direction of earthquake force  e =90 ∘ and the pseudo-static acceleration coefficient Figure 12 shows the effects of the direction of earthquake force ( e ) on the unique DOS of the critical slip surface.It can be seen from ( 6) that the components of seismic force in the direction of y-axis Q hyi increase with the increase of  e , which results in an increase of the DOS as is shown in Figure 12.When the earthquake force acts parallel to the x-axis ( e =0, Q hyi =0), the DOS is equal to 0 ∘ which is the same as predicated.It can be seen that as  e and K h increased, the value of the DOS increased linearly.Neglecting the effects of the direction of earthquake force ( e ) can lead to miscalculation of the DOS of the critical slip surface.Typically, when  e =90 ∘ and K h =0.3, the miscalculation of the DOS can be up to over 40 ∘ .
Figure 13 illustrates the difference in the unique DOS of critical slip surface between different geomechanical properties (  ,   ).The results of Figure 9 point toward the same trends that the unique DOS increases as the direction of earthquake force ( e ) increases.It also found out that all the values DOS are consistent with each other, the maximum difference between which are less than 0.5 ∘ .It means that the geomechanical properties have little effects on the unique DOS of critical slip surface.
The calculated results indicated that the 3D analysis yields a conservative estimate for the factor of safety when the direction of the earthquake force is assumed to parallel to the x-axis.It is conservative because the effects of the direction of earthquake force are not included in the analysis.The 3D analysis neglecting the direction of earthquake force is appropriate for slope design because it yields a conservative estimate for the factor of safety.A 3D analysis considering the effects of direction of earthquake force is recommended for the backanalysis of 3D slope failures so that the backcalculated shear strength reflects the effects of direction of earthquake force.The backcalculated shear strength then can be used in remedial measures for failed slopes or slope design at sites with similar conditions.

Conclusion
Simplified 3D Spencer's method with a random search technique was proposed to assess the stability of slopes under quasistatic load.A computer code was developed in Matlab to carry out the calculations.In addition, a meaningful comparison is made to investigate the effects of the direction of the earthquake force on the factor of safety and the DOS.
Based on the results presented, the following conclusions can be made: (1) The differences in factor of safety F/F 0 increase with an increase in the pseudostatic acceleration coefficient K h and the direction of earthquake force  e values.Neglecting the effects of  e on the factor of safety can lead to overestimation in 3D slope stability analysis.
(2) The direction of earthquake force  e has obvious influence on the unique DOS of the critical slip surface.The DOS of the critical slip surface increases linearly with an increase in  e and K h values.(3) The effective cohesion   and friction angle   of the soil have nearly no influence on the unique DOS of the critical surface.

Figure 3 :
Figure 3: Geometry and properties of Example 1.

Figure 4 :
Figure 4: The calculated critical slip surface compared with those provided by Zhang [16].

Figure 5 :
Figure 5: Geometry and properties of Example 2.

Figure 7 :
Figure 7: Geometry and properties of Example 3.

Table 1 :
Comparisons of F and critical slip surface for Example 1 (Bishop's method).
Note.The same special ellipsoid slip surface (R y =R x , =0) used by Zhang is adopted in this paper.

Table 2 :
Comparisons of F and critical slip surface for Example 1 (Janbu's simplified method).Note.The same special ellipsoid slip surface (R y =R x , =0) used by Zhang is adopted in this paper.

Table 3 :
Comparisons of F and critical slip surface for Example 1 (Spencer's method).
Note.The same special ellipsoid slip surface (R y =R x , =0) used by Zhang is adopted in this paper.

Table 4 :
Comparison of F for Example 2.

Table 5 :
Comparison of F and DOS for Example 3.