Seismic analysis of slope considering log-spiral failure surface with numerical validation

This paper presents the solution of c-φ soil slope stability problem using the pseudo-dynamic method assuming the log-spiral failure mechanism. The vertical slice analysis methodology is introduced to analyze the stability of the slope under vertical and horizontal seismic loads. The adopted pseudo-dynamic method in the present analysis considers the effect of phase difference and primary and shear wave velocity which travel through the soil medium during an earthquake. Considering the limit equilibrium method, the moment equilibrium condition for each slice was established and a factor of safety of the slope is evaluated. Results are presented in a tabular and graphical form to understand the effect of different parameters of slope soil. Comparisons of the present analysis are highlighted with the existing solutions. A numerical solution is also done by PLAXIS 2D to compare the present study under earthquake loading conditions.

The slope stability analysis problem and its solution under earthquake loading conditions have been widely published in the literature for many years. The seismic behavior of soil slope by pseudo-static approaches has been reported by Terzaghi [8]. In the pseudo-static analyses, the loads coming from the earthquakes is represented by the static inertia forces which act at the center of the failure soil mass. Several investigations [9][10][11][12][13][14][15] has been done to analyze the stability of slope by means of limit equilibrium technique at earthquake loading conditions. However, a pseudo-static method considers the dynamic behavior of an earthquake in a very approximate way and does not consider the effects of time and phase difference of seismic waves.
The displacement-based analysis has been investigated by Newmark [16] on dams and embankments due to earthquake shaking. Sarma [17] proposed a displacement analysis based on a limit equilibrium principle using Newmark model. Ling et al. [18] developed a procedure based on pseudo-static limit equilibrium which incorporates the permanent displacement limit. Ling et al. [19] analyzed the slope of planer surfaces considering the influence of horizontal and vertical seismic accelerations for the calculation of the permanent displacement. Lee et al. [20] proposed a analysis to estimate the permanent displacement of slope due to earthquake using Newmark's sliding block method and response-history analysis. However, Newmark sliding-block method adopted assumptions similar to the pseudo-static case as the effect of time and phase difference are not taken into account for the stability analysis of slopes.
To overcome the deficit of pseudo-static analysis, Steedman and Zeng [21] introduced a pseudo-dynamic method on retaining wall considering the effects of time and finite shear wave velocity. Further, Choudhury and Nimbalkar [22] extended the pseudo-dynamic approach taking into account both shear and primary wave velocities. The pseudo-dynamic approach has been applied by several researchers [23][24][25][26][27][28][29] for different structures to overcome the drawback of the pseudo-static approach. Basha and Babu [30] reported a reliability-based optimization technique of reinforced soil structure supporting φ backfill based on the pseudo-dynamic approach. Chanda et al. [31] introduced this pseudo-dynamic approach to solve the slope stability problem considering circular failure mechanism for homogeneous slope in which Fellenius line is used to locate the center of the most critical circle. On the other hand, Xu et al. [32] developed a pseudo-dynamic method to analyze the slope by wedge-failure mechanism during earthquake loading conditions. A pseudo-dynamic approach [33] is given using the discretization-based kinematic analysis to investigate the seismic stability of the slope. A gravity retaining wall is analyzed by the pseudo-dynamic method of analysis against sliding and overturning modes of failure [34]. A pseudo-dynamic method is suggested for the slope considering a circular failure mechanism [35]. However, evaluation of the stability of slopes considering log-spiral rupture surface under earthquake loading condition has received little attention. In this study, an attempt is made to develop a methodology to solve homogeneous c-φ soil slope considering log-spiral failure surface using the pseudo-dynamic approach. The effects of both the horizontal and vertical seismic accelerations acting on the potential failure mass are taken into account. The weight zone is divided into a number of vertical slices and the limit equillibrium condition is applied for both individual and for all the slices as a whole. Consequently, a numerical solution has been suggested by PLAXIS 2D [36] to compare the results obtained from analytical solutions. Finally, a factor of safety of slope under seismic loading condition is determined and the values as obtained from the analytical solution are compared with the values obtained from the numerical solution.

Assumptions
• The soil mass is considered rigid. • The soil is a homogeneous particulate material.
• The failure surface of the slope is assumed as a logarithmic spiral. • The interactions between the slices are neglected as the resultant force is parallel to the base of each slice. • The seismic forces are considered as sinusoidal harmonic motion.

Method of analysis
Let us consider an arbitrary soil slope of height H, inclined at an angle β with the horizontal made up of c-φ soil as shown in Fig. 1. Assuming rotational failure passing through the toe which generates logarithmic spiral failure arc. The failed soil mass is divided into a number of vertical slices. This logarithmic spiral surface is governed by the height of slope H and the location of the center of the logarithmic spiral arc 'O' . OE and OA are the initial and final radius respectively of the logarithmic spiral curve. The location of the log-spiral arc can be defined by the subtended angle θ 1 . ' AE' is the straight length of the log spiral wedge which makes an angle ξ with the horizontal line PP 1 .
The generalized equation for the log-spiral curve is given by The initial radius r 0 is given by The width of the failed mass of the slope at the top surface, The failure wedge is divided into 'n' number of slices. The wedge is divided into two zones in which the zone ADB is divided into 'p' number of slices and zone BDE is divided into 'n-m' number of slices. The 1 st and nth slices are triangular in shape whereas, the shapes of the other slices are trapezoidal.
The forces acting on ith slice and jth slice during equilibrium is shown in Fig. 2. 'b 1 ' and 'b 2 ' is the width of ith and jth slice respectively. The forces acting on each slice include weight (W i , W j ) acting vertically at the midpoint of the slice in the direction normal to the arc, cohesive force (c) acting along failure surface in the direction opposite to the movement of the failure mass, the normal (N i , N j ) and tangential forces (T i , T j ) acting on the lower boundary of the slice, Reaction (R i , R j ) inclined at an angle φ with the normal of the slice.
The normal and tangential forces on ith slice are given by, where α i is the angle between N i and W i The angle between normal N and weight W for ith slice is given by,

Calculation of pseudo-dynamic inertia forces
During earthquake condition, the shear and primary waves transmitting through the soil mass at a velocity are v s = √ G/ρ and v p = √ 2G(1 − ν)/ρ(1 − 2ν) respectively where G, ρ,µ and ω are shear modulus, density, Poisson's ratio and angular frequency respectively.
According to Choudhury and Nimbalkar [22], if the base of the slope is subjected to harmonic horizontal and vertical seismic accelerations of amplitude k h g and k v g, the accelerations at any depth z and time t below the top can be expressed as,

Derivation of horizontal and vertical inertia forces
The horizontal inertia force Q h 1 acting on the slice-1 in the zone ABD can be expressed Again, the horizontal inertia force Q h i acting on the slice-2 to slice-p in the zone ABD can be expressed as, where, Similarly, the vertical inertia force Q v 1 acting on the slice-1 in the zone ABD can be expressed as, Also, the vertical inertia force Q v i acting on the slice-2 to slice-p in the zone ABD can be expressed as, On the other hand, height of slice in the zone BDE is given by, For any jth slice, the angle between normal N and weight W is given by,

Mass of jth slice in zone BDE,
Mass of nth slice in zone BDE, Now, the horizontal inertia force Q hj acting on the slice-(p + 1) to slice-(n-1) in the zone BDE can be expressed as, Besides, the horizontal inertia force Q h n acting on the slice-n in the zone BDE can be expressed as, In a similar way, the vertical inertia force Q v j acting on the slice-(p + 1) to slice-(n-1) in the zone BDE is given by, Also, the vertical inertia force Q v acting on the slice -n can be expressed as,

Expression of horizontal and vertical inertia forces due to surcharge
The horizontal inertia force due to the surcharge load (q) is given by, where m q is the mass of the elemental strip in the surcharge load, a hq denotes horizontal inertia force due to surcharge, z q = equivalent surcharge height = q γ (Fig. 2). Solving Eq. (29) the horizontal inertia force for any jth slice is given by, Similarly, the vertical inertia forces are expressed as, (24)

Determination of factor of safety
The horizontal disturbing moment M H is expressed as, The vertical disturbing moment M V is expressed as, The disturbing moment due to tangential component (T) is given by, The resisting moment of soil slope is given by, where, m f denotes the mobilization factor. There, factor of safety (FOS) is given by,

Results and discussions
In this analysis, all the parameters of soil slope are constant except r 0 , θ 1 , t/T and m f . The optimization is done by MATLAB [37] and to evaluate the factor of safety with respect to θ 1 , θ 2 and t/T. Figure 3 shows the flowchart for obtaining the factor of safety. From the global optimization curve, the minimum factor of safety is taken as the optimized value of factor of safety. To study the effect of pseudo-dynamic seismic inertia forces on the slope, the parametric study is conducted for the variation of different parameters. The results obtained from the present analysis are presented in the form of a table and graph. The variations of parameters of the typical model slope are as follows:  Tables 1 and 2 show the optimized factor of safety at both static and seismic loading conditions. The following subsections present the variation of factor of safety of slope due to the variation of different parameters.    11:12 Steps to use the results to solve field problem

Start
The steps given below are to be followed.
i Given data: Soil parameters (c, φ, γ); Surcharge loading, q; slope parameters (desired height, H and inclination of the slope, β); Seismic accelerations (k h and k v ). j Evaluate c/γH and q/γH. k For the value of φ, β, k h , k v , c/γH and q/γH, find out the value of FOS from Tables Tables 1 and 2.
For intermediate portion, linear interpolation is suggested. If the FOS is more than one, then the slope is safe. On the other hand if the FOS is less than one, then either the slope parameters H, β are to be changed untill unless FOS reaches the value ~ one or reinforcements are to be provided following Ghosh and Paul [38]. Figure 4 shows that variation of factor of safety with respect to horizontal seismic acceleration (k h ) for k v = 0, β = 20°, φ = 40°, c/γH = 0.05 and q/γH = 0. It is observed that an increase in k h from 0.0 to 0.1, decrease the value of factor of safety by 21.5% Table 2 Results of pseudo-dynamic factor of safety at c/γH = 0.05 and q/γH = 0    Fig. 5 it is seen that, increase in k v from 0 to k h /2 factor of safety increase by 3% for k h = 0.1, β = 20°, φ = 40°, c/γH = 0.05 and q/γH = 0. Due to an increase in seismic acceleration, the disturbance of the soil increases which results in a reduction in the value of factor of safety. Figure 6 shows the variation of factor of safety with respect to horizontal seismic acceleration (k h ) at different angle of internal friction (20°, 30° and 40°) of soil. From the plot, it is seen that increasing angle of internal friction from 20° to 30° factor of safety gets increased by 16% at k h = 0.1, β = 20°, k v = k h /2, c/γH = 0.05 and q/γH = 0. It is described that factor of safety increases with the increase of angle of internal friction and increases the internal resistance of the soil particle which allows an increase in factor of safety.

Effect on factor of safety due to the variation of slope angle
Factor of safety decreases with the increase in the value of slope angle. Figure 7 shows the variation of factor of safety with different slope angle (β = 20°, 35° and 50°) for

Influence of cohesion on factor of safety
It is seen that factor of safety increases with an increase in c/γH ratio. For example, at φ = 40°, β = 20°, k h = 0.1, k v = 0 and q/γH = 0 factor of safety increases approximately 14% and 22% due to increase of c/γH for 0.05 to 0.1 and 0.05 to 0.15 respectively. Cohesion increases the intermolecular attraction between the soil particles and hence the factor of safety of slope increases.

Influence of surcharge on the factor of safety
The variations of factor of safety with various surcharge loading are plotted in Fig. 9. From the plot, it is observed that factor of safety decreases with the increase of surcharge loading. As an example, for φ = 40°, β = 20°, k h = 0.1, k v = k h /2 and c/γH = 0.05, the factor of safety increases about 4% when surcharge coefficient changes from 0 to 1. With the increase in the surcharge, the driving forces acting on the slope increase and ultimately factor of safety decreases.

Comparison of results
The results of factor of safety obtained by the present pseudo-dynamic analysis are compared with the factor of safety as obtained by existing solutions of pseudostatic and pseudo-dynamic analyses. Table 3  of safety obtained by the present method with pseudo-static analyses results evaluated by Newmark [16] and Choudhury et al. [11]. It is noticed that results of factor of safety obtained from present study lower than those obtained by Newmark [16] and Choudhury et al. [11] as present study assumed log-spiral failure surface whereas, Newmark [16] and Choudhury et al. [11] assumed planer and circular failure surfaces respectively.
On the other hand, Table 4 shows the values of factor of safety obtained by the present study and the values reported in Chanda et al. [31]. It is observed that factor of safety obtained from the present analysis is marginally lower than Chanda et al. [31]. The marginal difference is that Chanda et al. [31] analyses pseudo-dynamic analysis considering circular rupture surfaces whereas the present study considers log-spiral rupture surfaces. For example, at φ = 40°, β = 35°, k h = 0.1, k v = k h /2, the factor of safety obtained by present analysis are lower about 5% than values obtained considering circular rupture surfaces.

Numerical modelling
Slope stability analysis by finite element techniques has been widely used in the last few decades. The advantages of the finite element method over other methods are that no failure mechanism has to be assumed a priori. Pasternack and Gao [39] reported a numerical procedure of slope for the determination of factor of safety. Mastoni and San [40] developed a shear strength technique for the finite element slope stability problem. Ugai and Leshichinsky [41] reported a comparison between the three-dimensional limit equilibrium approach and finite element method of the slope. Duncun [42] developed a hyperbolic non-linear elastic stress-strain relationship to express the dependence of soil behavior on the shear strength parameters. Several studied [43][44][45] of slope stability analysis based on finite element method are reported Hammouri et al. [46] reported a slope stability analysis using limit equilibrium method and results obtained from the limit equilibrium approach compared with the results obtained from finite element analyses.

Development of a finite element model
In the present study, a two-dimensional finite element model is developed by PLAXIS [36] to analyzed the slope under earthquake loading conditions. The geometry of the finite element model slope is illustrated in Fig. 10. The standard boundary of the slope is considered at the base at a depth 15 m below the top level of slope and roller at the two vertical sides; one vertical side passes through the center line of slope and another at a distance of 40 m away from the toe of the slope. Three slopes with different slope angles are modeled in the present analysis. A load of 50kN/m is applied as a surcharge load on the top surface of the slope. Figure 11 shows the arrangement of mesh with significant nodes before the analysis. The slope of the present study is modeled with axisymmetric plane strain and 6-noded triangular elements to determine the minimum factor of safety. The initial horizontal to vertical effective stress ratio (K 0 ) is assumed as 0.5 based on Jaky's formula.

Material model and soil properties
The linear elastic perfectly-plastic Mohr-Coulomb soil model is adopted in the present numerical analysis to model the behavior of the soil. A total five number of parameters are required to model the Mohr-Columb criteria. The parameters for Mohr-Coulomb model used in the present analysis are listed in Table 5.

Calculation procedure
The calculation of the present analysis involves three phases: • In the first phase, the only static load is applied.
• The dynamic analysis has been performed in the second phase in which an earthquake load (Fig. 12) is introduced in addition to the static load of the first phase. • In the third phase, the phi-c reduction technique is applied to evaluate the global factor of safety of the model slope considering plastic equilibrium criteria are to be satisfied.

Results and discussions of numerical analysis
The results obtained from the numerical approach and sensitivity analysis are presented in the following sections.

Study of deformed mesh
The generated deformed mesh is shown in Fig. 13. From the arrangement of the mesh, it is observed that the mesh is deformed in a vertical direction under the surcharge loading. As the distance from the surcharge loading increases, the vertical deformation due to the surcharge loading gradually going to be decreased. Also, the maximum horizontal deformation occurs at the sloping side of the slope as shown in the plot.   11:12 Response of slope due to acceleration Figure 14a shows the peak acceleration value at different depth of slope at a distance of 6.5 m from AB face (Fig. 10) whereas Fig. 14b presents the acceleration vs time curve of slope at top and bottom surface of the slope after the application of acceleration at the fixed end of the slope. The peak value of applied acceleration is 1.073 m/s 2 and the corresponding k h value is 0.1. From the plot, it is seen that the peak value of acceleration increase towards the direction of the surface and at the top surface of the slope it becomes maximum (1.027 m/s 2 ). It is obvious as the soil at the top surface is  free to move. It is also seen from the plot at the fixed boundary of the slope there is no acceleration.

Horizontal and Vertical deformation of the slope
The horizontal deformation and corresponding contour lines are shown in Fig. 15a, b respectively due to the application of the earthquake loading. From the plot, it is seen that the soil at the sloping face displaced at a larger amount indicating slope type of rotational failure of the soil mass. The range of horizontal displacement for slope failure is greater than 22.5 mm. When the displacement contour passes through toe it generates 15 mm displacement and if displacement contour passes through base the related displacement is 12.5 mm or less. On the other hand, Fig. 16a, b shows the variation of vertical displacement and contour lines respectively. From the plot, it is noticed that the maximum vertical deformation occurs just below the top surface of the slope at which point the load is applied. It is also seen that the vertical deformation gradually decreases with the depth and nearly after two-third of the depth, the vertical deformation completely disappears.

Variation of total stress
The variation of total stresses and corresponding contour lines due to the application of earthquake is presented in Fig. 17a, b respectively. From the diagram, it is observed that due to the applied force maximum compressive stresses occur at the bottom middle part of the slope and it gradually decreases towards the top and sloping side of the slope.

Verification of the present analysis
A two-dimensional finite element analysis has been performed to compare the present analytical solutions. Table 6 presents the comparison of results obtained by the analytical solution with results obtained from finite element analysis. From Table, it is seen that a good match is observed between the results obtained from analytical solutions and finite element analysis. The comparison of the potential failure surfaces between the analytical approach and numerical analysis are shown in Fig. 18.
According to Fig. 14, the peak ground acceleration found after the numerical analysis is 1.05 m/s 2 . In the analytical, the applied peak value is given by,

Conclusions
An attempt is made to evaluate the factor of safety of soil slope using a pseudo-dynamic limit equilibrium method and considering a log-spiral failure surface. An attempt is also made to compare the analytical solution with the results obtained from the PLAXIS 2D solution. A detailed parametric study is done for the variation of different parameters in the case of both analytical and numerical solutions. It is seen that factor of safety decreases with the increase of horizontal and vertical seismic accelerations, slope angle and surcharge loading. On the other hand, with the increase of soil friction angle and cohesion factor of safety increases. As the increase of depth, the acceleration of the slope decrease and the acceleration at the top surface becomes maximum. It is also observed that maximum horizontal deformation occurs at the sloping side of the slope, whereas slope deformed maximum vertical deformation occurs just below the load is applied. Slope experiences maximum compressive stress at the bottom middle part of the slope. The results obtained from the present analysis are compared with the previous pseudo-static and pseudo-dynamic analyses of slope and found to be a lower value than previous analyses. After optimization, the center PGA = k h × g = (0.1 × 9.8) = 0.98 m/s 2 Fig. 16 Plot of vertical displacements β = 20° a Shadings, b Contour lines of rotation is determined and finally, the initial and final radius are known. Based on those values, the potential failure surface is plotted to show the shape of the slip surface under critical condition and a comparison with the numerical and experimental failure surfaces proves that the logarithmic failure surface is feasible in comparison to linear or circular failure surface. The results as obtained from the present study may be used for stability analysis of slope under seismic loading condition. As a future scope of the work, it is mentioned here that the present suggested analytical methodology can easily be extended for nonhomogeneous soil slope even if there are an arbitrary number of non-homogeneous blocks of soil within the failure wedge.
The present methodology uses an analytical method where earthquake forces are considered in an approximate way without writing down the displacement equation under a real dynamic environment which ultimately moves to numerical analysis.