Introduction

Determination of seismic slope stability is a crucial topic for Geotechnical Engineer. The traditional concept of slope stability is to evaluate the minimum factor of safety. The slope may be natural or man-made. The aim of the slope stability is to determine the safe and economic design of embankments, excavations, landfills etc. Several studies are published related to the slope stability and its analyses.

The traditional and well-established limit equilibrium method is widely used for slope stability analysis by engineers and researchers. The most common limit equilibrium techniques are Swedish circle method or method of slice in which the failure wedges are divided into the number of slices. Fellenius [1] established a line which is called Fellenius line to locate the center of the critical circle which gives a minimum factor of safety. Bishop [2] upgraded the method of slice in which inter-slice forces are taken into account. Consequently, the same method with noncircular failure surfaces was upgraded by Morgenstern and Price [3], Spencer method [4], Janbu method [5, 6] and Zhu et al. [7] based on different assumptions of inter slice forces.

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.

Stability analysis of slope under earthquake loadings

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 PP1.

Fig. 1
figure 1

Model of slope considered in the present study

Full size image

The generalized equation for the log-spiral curve is given by

$$r = r_{0} e^{{\theta_{1} \tan \varphi }}$$
(1)

The initial radius \(r_{0}\) is given by

$$r_{0} = \frac{{H\cos ec(\theta_{1} + \theta_{2} )}}{{e^{{\theta_{1} \tan \theta }} - \frac{{\sin \theta_{2} }}{{\sin (\theta_{1} + \theta {}_{2})}}}}$$
(2)

The width of the failed mass of the slope at the top surface,

$${\text{BE}} = b_{s} = \frac{{r_{0} \sin \theta_{1} }}{{\sin (\theta_{1} + \theta {}_{2})}} - H\cot (\theta_{1} + \theta_{2} ) - H\cot \beta$$
(3)

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 1st 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. ‘b1’ and ‘b2’ is the width of ith and jth slice respectively. The forces acting on each slice include weight (Wi, Wj) 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 (Ni, Nj) and tangential forces (Ti, Tj) acting on the lower boundary of the slice, Reaction (Ri, Rj) inclined at an angle φ with the normal of the slice.

Fig. 2
figure 2

Generalization of slice and forces considered on each slice

Full size image

The normal and tangential forces on ith slice are given by,

$$N_{i} = W_{i} \cos \alpha$$
(4)
$$T_{i} = W_{i} \sin \alpha_{i}$$
(5)

where \(\alpha_{i}\) is the angle between Ni and Wi

The angle between normal N and weight W for ith slice is given by,

$$\alpha_{i} = \sin^{ - 1} \frac{{\left( {2i - 1} \right)\frac{{b_{1} }}{2}}}{{r_{0} e^{{\theta_{1} \tan \varphi }} }}$$
(6)

The value of bo is given by,

$$b_{o} = AP = r\cos \left( {\theta_{1} + \theta_{2} } \right)$$
(7)

Height of ith slice in the zone ABD is given by

$$h_{i} = ib{}_{1}\left( {\tan \beta - \tan \xi } \right)$$
(8)

where

$$\xi = \tan^{ - 1} \left( {\frac{H}{{H\cot \beta + b_{s} }}} \right).$$
(9)

Weight of ith slice, is given by,

$$W_{i} = A_{i} \gamma$$
(10)

Mass of 1st slice in zone ADB,

$$m(z)_{1} = \frac{\gamma }{g}\frac{{b_{1} z_{1} }}{{h_{1} }}dz$$
(11)

Mass of ith slice in zone ADB

$$m(z)_{i = 2 \, to \, p} = \frac{{\gamma b_{1} }}{g}dz$$

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} = \sqrt {G/\rho }\) and \(v_{p} = \sqrt {2G(1 - \nu )/\rho (1 - 2\nu )}\) respectively where G,\(\rho ,\)\(\mu\) and \(\omega\) 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 khg and kvg, the accelerations at any depth z and time t below the top can be expressed as,

$$a_{h} \left( {z,t} \right) = k_{h} g\sin \omega \left( {t - \frac{h - z}{{v_{s} }}} \right)$$
(13)
$$a_{v} \left( {z,t} \right) = k_{v} g\sin \omega \left( {t - \frac{h - z}{{v_{p} }}} \right)$$
(14)

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 as,\(Q_{{h_{1} }} = \int\limits_{z = 0}^{{z_{1} = h_{1} }} {(q_{h} )_{1} } dz_{1} = \int\limits_{{z_{1} = 0}}^{{z_{1} = h_{1} }} {m(z)_{1} .a_{h} \left( {z_{1} ,t} \right)dz_{1} }\)[\(h_{1} =\)\(b_{1} \left( {\tan \beta - \tan \xi } \right)\) from Fig. 1]

$$= \frac{{\gamma b_{1} k_{h} \lambda }}{{4\pi^{2} h_{1} }}\left[ {\lambda \sin \left( {2\pi \frac{t}{T}} \right) - \lambda \sin \left\{ {2\pi \left( {\frac{t}{T} - \frac{{h_{1} }}{\lambda }} \right)} \right\} - 2\pi h_{1} \cos \left( {2\pi \frac{t}{T}} \right)} \right]$$
(15)

Again, the horizontal inertia force \(Q_{{h_{i} }}\) acting on the slice-2 to slice-p in the zone ABD can be expressed as,

$$(Q_{h} ){\mkern 1mu} _{{i_{{ = 2{\text{ to}}p}} }} = \int\limits_{{z_{i} = 0}}^{{z_{i} = h_{{a_{i} }} }} {(q_{h} )_{{i = {\kern 1pt} 2{\kern 1pt} to{\kern 1pt} p}} } dz_{i} = \int\limits_{{z_{i} = 0}}^{{z_{i} = h_{{a_{i} }} }} {m(z_{i} )_{{i = 2{\kern 1pt} to{\kern 1pt} p}} } a_{h} \left( {z_{i} ,t} \right)dz_{i} {\text{ }} = \frac{{\gamma b_{1} k_{h} \lambda }}{{2\pi }}\left[ {\cos \left\{ {2\pi \left( {\frac{t}{T} - \frac{{h_{{a_{i} }} }}{\lambda }} \right)} \right\} - \cos \left( {2\pi \frac{t}{T}} \right)} \right]$$
(16)

where,

$$h_{{a_{i} }} = \frac{{h_{i - 1} + h_{i} }}{2} = \frac{{(2i - 1)b_{1} }}{2}\left( {\tan \beta - \tan \xi } \right)$$
(17)

Similarly, the vertical inertia force \(Q_{{v_{1} }}\) acting on the slice-1 in the zone ABD can be expressed as,

$$Q_{{v_{1} }} = \frac{{\gamma b_{1} k_{v} \eta }}{{4\pi^{2} h_{1} }}\left[ {\eta \sin \left( {2\pi \frac{t}{T}} \right) - \eta \sin \left\{ {2\pi \left( {\frac{t}{T} - \frac{{h_{1} }}{\eta }} \right)} \right\} - 2\pi h_{1} \cos \left( {2\pi \frac{t}{T}} \right)} \right]$$
(18)

Also, the vertical inertia force \(Q_{{v_{i} }}\) acting on the slice-2 to slice-p in the zone ABD can be expressed as,

$$(Q_{v} )_{{i = 2\,{\text{to}, {p}}}} = \frac{{\gamma b_{1} k_{v} \eta }}{{2\pi }}\left[ {\cos \left\{ {2\pi \left( {\frac{t}{T} - \frac{{h_{{a_{i} }} }}{\eta }} \right)} \right\} - \cos \left( {2\pi \frac{t}{T}} \right)} \right]$$
(19)

On the other hand, height of slice in the zone BDE is given by,

$$h_{j} = \left\{ {b_{s} - \left( {j - p} \right)b_{2} } \right\}\tan \xi$$
(20)

For any jth slice, the angle between normal N and weight W is given by,

$$\alpha_{i} = \sin^{ - 1} \left[ {\frac{{pb_{1} + \left\{ {2\left( {j - p} \right) - 1} \right\}\frac{{b_{2} }}{2}}}{{r_{0} e^{{\theta_{1} \tan \varphi }} }}} \right]$$
(21)

Mass of jth slice in zone BDE,

$$m(z)_{{j = (p + 1)\,{\text{to}}\,{\text{(n}} - {\text{1)}}}} = \frac{{\gamma b_{2} }}{g}dz_{j}$$
(22)

Mass of nth slice in zone BDE,

$$m(z)_{n} = \frac{\gamma }{g}\frac{{b_{2} z_{n} }}{{h_{n} }}dz_{n}$$
(23)

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,

$$\begin{aligned} (Q_{h} )_{{j = \left( {p1 + } \right){\kern 1pt} to{\kern 1pt} \left( {n - 1} \right)}} = & \int\limits_{{z_{j} = 0}}^{{z_{j} = h_{{a_{j} }} }} {(q_{h} )_{{j = {\kern 1pt} \left( {p + 1} \right){\kern 1pt} to{\kern 1pt} \left( {n - 1} \right)}} } dz_{j} = \int\limits_{{z_{j} = 0}}^{{z_{j} = h_{{aj}} }} {m(z)_{{j = \left( {p + 1} \right){\kern 1pt} to{\kern 1pt} \left( {n - 1} \right)}} } a_{h} \left( {z_{j} ,t} \right)dz_{j} \\ = & \frac{{\gamma b_{2} k_{h} \lambda }}{{2\pi }}\left[ {\cos \left\{ {2\pi \left( {\frac{t}{T} - \frac{{h_{{a_{j} }} }}{\lambda }} \right)} \right\} - \cos \left( {2\pi \frac{t}{T}} \right)} \right] \\ \end{aligned}$$
(24)

where,

$$h_{{a_{j} }} = \frac{{h_{j - 1} + h_{j} }}{2} = \left[ {b_{s} - \frac{{\left\{ {2\left( {j - p} \right) - 1} \right\}b_{2} }}{2}} \right]\tan \xi$$
(25)

Besides, the horizontal inertia force \(Q_{{h_{n} }}\) acting on the slice-n in the zone BDE can be expressed as,

$$\begin{aligned} Q_{{h_{n} }} = & \int\limits_{{z_{n} = 0}}^{{z_{n} = h_{n} }} {(q_{h} )_{n} } dz_{j} = \int\limits_{{z_{n} = 0}}^{{z_{n} = h_{n} }} {m(z_{n} )a_{h} \left( {z_{j} ,t} \right)dz_{j} } \\ = & \frac{{\gamma b_{2} k_{h} \lambda }}{{4\pi^{2} h_{n} }}\left[ {\lambda \sin \left( {2\pi \frac{t}{T}} \right) - \lambda \sin \left\{ {2\pi \left( {\frac{t}{T} - \frac{{h_{n} }}{\lambda }} \right)} \right\} - 2\pi h_{n} \cos \left( {2\pi \frac{t}{T}} \right)} \right] \\ \end{aligned}$$
(26)

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,

$$Q_{{v_{j} }} = \frac{{\gamma b_{2} k_{v} \eta }}{2\pi }\left[ {\cos \left\{ {2\pi \left( {\frac{t}{T} - \frac{{h_{{a_{j} }} }}{\eta }} \right)} \right\} - \cos \left( {2\pi \frac{t}{T}} \right)} \right]$$
(27)

Also, the vertical inertia force \(Q_{v}\) acting on the slice -n can be expressed as,

$$Q_{{v_{n} }} = \frac{{\gamma b_{2} k_{v} \eta }}{{4\pi^{2} h_{n} }}\left[ {\eta \sin \left( {2\pi \frac{t}{T}} \right) - \eta \sin \left\{ {2\pi \left( {\frac{t}{T} - \frac{{h_{n} }}{\eta }} \right)} \right\} - 2\pi h_{n} \cos \left( {2\pi \frac{t}{T}} \right)} \right]$$
(28)

Expression of horizontal and vertical inertia forces due to surcharge

The horizontal inertia force due to the surcharge load (q) is given by,

$$Q_{hq} = \int\limits_{z = 0}^{{z = z_{q} }} {m_{q} \left( z \right)a_{hq} \left( {z,t} \right)dz} = \int\limits_{z = 0}^{{z = z_{q} }} {\frac{\gamma }{g}b_{s} .k_{h} g\sin \omega \left( {t - \frac{{H + z_{q} - z}}{{v_{s} }}} \right)dz}$$
(29)

where mq 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 = \(\frac{q}{\gamma }\) (Fig. 2).

Solving Eq. (29) the horizontal inertia force for any jth slice is given by,

$$\left( {Q_{hq} } \right)_{j} = \frac{{\gamma b_{2} k_{h} }}{2\pi }\left[ {\lambda \cos 2\pi \left( {\frac{t}{T} - \frac{{h_{j} + z_{q} }}{\lambda }} \right) - \cos 2\pi \left( {\frac{t}{T} - \frac{{z_{q} }}{\lambda }} \right)} \right]$$
(30)

Similarly, the vertical inertia forces are expressed as,

$$\left( {Q_{vq} } \right)_{j} = \frac{{\gamma b_{2} k_{v} }}{2\pi }\left[ {\eta \cos 2\pi \left( {\frac{t}{T} - \frac{{h_{j} + z_{q} }}{\eta }} \right) - \cos 2\pi \left( {\frac{t}{T} - \frac{{z_{q} }}{\eta }} \right)} \right]$$
(31)

Determination of factor of safety

The horizontal disturbing moment MH is expressed as,

$$M_{H} = k_{h} \left\{ {\sum\limits_{i = 1}^{p} {W_{i} \bar{y}_{i} } + \sum\limits_{i = p + 1}^{n} {\left( {W_{j} + Q_{{hq_{j} }} } \right)\bar{y}_{j} } } \right\}$$
(32)

The vertical disturbing moment MV is expressed as,

$$M_{V} = \left( {1 \pm k_{v} } \right)\left\{ {\sum\limits_{i = 1}^{p} {W_{i} \bar{x}_{i} } + \sum\limits_{i = p + 1}^{n} {\left( {W_{j} + Q_{{vq_{j} }} } \right)\bar{x}_{j} } } \right\}$$
(33)

\(\bar{x}_{1}\), \(\bar{x}_{2}\), \(\bar{x}_{i}\), \(\bar{x}_{j}\), \(\bar{x}_{n}\), \(\bar{y}_{1}\), \(\bar{y}_{2}\), \(\bar{y}_{i}\), \(\bar{y}_{j}\) and \(\bar{y}_{n}\) denotes the point of action w.r.t. center O for different slices.

The disturbing moment due to tangential component (T) is given by,

$$M_{T} = r\left[ {\sum\limits_{i = 1}^{p} {\left( {W_{i} \sin \alpha_{i} } \right) + \sum\limits_{j = p + 1}^{n} {\left( {W_{j} \sin \alpha_{j} } \right)} } } \right] = r_{0} e^{{\theta_{1} \tan \varphi }} \left[ {\sum\limits_{i = 1}^{p} {\left( {W_{i} \sin \alpha_{i} } \right) + \sum\limits_{j = p + 1}^{n} {\left( {W_{j} \sin \alpha_{j} } \right)} } } \right]$$
(34)

The resisting moment of soil slope is given by,

$$\begin{aligned} M_{R} = r\left[ {c_{{m_{f} }} \left\{ {\sum\limits_{i = 1}^{p} {\left( {b_{1} \sec \alpha_{i} } \right)} + \sum\limits_{j = p + 1}^{n} {\left( {b_{2} \sec \alpha_{j} } \right)} } \right\} + \left\{ {\sum\limits_{i = 1}^{p} {\left( {W_{i} \cos \alpha_{i} } \right)} + \sum\limits_{j = p + 1}^{n} {\left( {W_{j} \cos \alpha_{j} } \right)} } \right\}\tan \varphi_{{m_{f} }} } \right] \hfill \\ = r_{0} e^{{\theta_{1} \tan \varphi }} \left[ {c_{{m_{f} }} \left\{ {\sum\limits_{i = 1}^{p} {\left( {b_{1} \sec \alpha_{i} } \right)} + \sum\limits_{j = p + 1}^{n} {\left( {b_{2} \sec \alpha_{j} } \right)} } \right\} + \left\{ {\sum\limits_{i = 1}^{p} {\left( {W_{i} \cos \alpha_{i} } \right)} + \sum\limits_{j = p + 1}^{n} {\left( {W_{j} \cos \alpha_{j} } \right)} } \right\}\tan \varphi_{{m_{f} }} } \right] \hfill \\ \end{aligned}$$
(35)

where, mf denotes the mobilization factor.

There, factor of safety (FOS) is given by,

$${\text{FOS}} = \frac{{M_{R} }}{{M_{H} + M_{V} + M_{T} }}$$
(36)

Results and discussions

In this analysis, all the parameters of soil slope are constant except r0, θ1, t/T and mf. 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:

Fig. 3
figure 3

Flowchart for the evaluation of FOS

Full size image

  • β = 20°, 35° and 50°

  • φ = 20°, 30° and 40°

  • kh= 0.1, 0.2 and 0.3

  • kv= 0.0, kh/2 and kh

  • c/γH = 0.05, 0.1 and 0.15

  • q/γH = 0, 0.5 and 1

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.

Table 1 Value of factor of safety at the static condition at c/γH = 0.05 and q/γH = 0
Full size table
Table 2 Results of pseudo-dynamic factor of safety at c/γH = 0.05 and q/γH = 0
Full size table

Steps to use the results to solve field problem

The steps given below are to be followed.

  1. i

    Given data: Soil parameters (c, φ, γ); Surcharge loading, q; slope parameters (desired height, H and inclination of the slope, β); Seismic accelerations (kh and kv).

  2. j

    Evaluate c/γH and q/γH.

  3. k

    For the value of φ, β, kh, kv, 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].

Parametric study

Effect on factor of safety due to the variation of horizontal and vertical seismic acceleration coefficients

Figure 4 shows that variation of factor of safety with respect to horizontal seismic acceleration (kh) for kv = 0, β = 20°, φ = 40°, c/γH = 0.05 and q/γH = 0. It is observed that an increase in kh from 0.0 to 0.1, decrease the value of factor of safety by 21.5% over the static value and an increase of kh from 0.1 to 0.2 decreases the factor of safety by 9% on kh = 0.1. Similarly, from Fig. 5 it is seen that, increase in kv from 0 to kh/2 factor of safety increase by 3% for kh = 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.

Fig. 4
figure 4

The variations of factor of safety with kh for β = 20°, c/γH = 0.05 and q/γH = 0

Full size image
Fig. 5
figure 5

Variations of factor of safety with kv for kh = 0.1, β = 20°, φ = 40°, c/γH = 0.05 and q/γH = 0

Full size image

Effect on factor of safety due to the variation of angle of internal friction (φ)

Figure 6 shows the variation of factor of safety with respect to horizontal seismic acceleration (kh) 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 kh = 0.1, β = 20°, kv = kh/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.

Fig. 6
figure 6

The variations of factor of safety with at different angle of internal friction for kh = 0.1, kv = kh/2, c/γH = 0.05 and q/γH = 0

Full size image

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 φ = 40°, kv = kh/2, c/γH = 0.05 and q/γH = 0. It is observed that if slope angle increase from angle 20° to 35° factor of safety gets decrease by about 12.5% and increase from angle 20° to 50° decrease by about 30% at kh = 0.1, kv = kh/2, φ = 40°, c/γH = 0.05 and q/γH = 0. Due to the increase of the slope angle, the driving forces increases which causes a decrease in the factor of safety.

Fig. 7
figure 7

The variations of factor of safety with slope angle for kh = 0.1, kv = kh/2, φ = 40°, c/γH = 0.05 and q/γH = 0

Full size image

Influence of cohesion on factor of safety

Figure 8 demonstrates the variations of factor of safety with respect to horizontal seismic acceleration at a different c/γH ratio (c/γH = 0.05, 0.1 and 0.15) for φ =40°, kv = kh/2, β = 20° q/γH = 0. It is seen that factor of safety increases with an increase in c/γH ratio. For example, at φ = 40°, β = 20°, kh = 0.1, kv = 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.

Fig. 8
figure 8

The variations of factor of safety with cohesion for kv = 0, φ =40°, β = 20° q/γH = 0

Full size image

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°, kh = 0.1, kv = kh/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.

Fig. 9
figure 9

The variations of factor of safety with different surcharge for φ = 40°, β = 20°, kv = kh/2 and c/γH = 0.05

Full size image

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 pseudo-static and pseudo-dynamic analyses. Table 3 presents the comparison of the factor 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.

Table 3 Comparison of present results with the results obtained from pseudo-static analyses in available literature for c/γH = 0.025, q/γH = 0
Full size table

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°, kh = 0.1, kv = kh/2, the factor of safety obtained by present analysis are lower about 5% than values obtained considering circular rupture surfaces.

Table 4 Comparison of present results with the results obtained from pseudo-dynamic analyses in the available literature
Full size table

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.

Fig. 10
figure 10

Plot of geometry with loads and boundary conditions (all dimension are in m)

Full size image
Fig. 11
figure 11

Plot of mesh with significant nodes before analysis (β = 20°)

Full size image

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 (K0) 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.

Table 5 Input soil parameters for PLAXIS
Full size table

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.

    Fig. 12
    figure 12

    Input earthquake data for dynamic analysis in Plaxis as provided at the fixed end of slope mode

    Full size image

  • 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.

Fig. 13
figure 13

The plot of deformed mesh due to the loadings (exaggeration factor: 1)

Full size image

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/s2 and the corresponding kh 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/s2). 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.

Fig. 14
figure 14

Plot of acceleration vs depth a peak value at different depth, b at top and bottom surface

Full size image

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.

Fig. 15
figure 15

Plot of horizontal displacements at β = 20° a Shadings b Contour lines

Full size image
Fig. 16
figure 16

Plot of vertical displacements β = 20° a Shadings, b Contour lines

Full size image

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.

Fig. 17
figure 17

The plot of total stresses β = 20° a Shadings, b Contour lines. Horizontal Distance (m)

Full size image

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.

Table 6 Comparison of present analytical results with numerical analysis
Full size table
Fig. 18
figure 18

Potential failure surface at β = 20°, kh = 0.1, kv = 0, φ = 40°, c/γH = 0.05 and q/γH = 0

Full size image

According to Fig. 14, the peak ground acceleration found after the numerical analysis is 1.05 m/s2. In the analytical, the applied peak value is given by,

$$PGA = \left( {k_{h} \times g} \right) = \left( {0.1 \times 9.8} \right) = 0.98\quad {\text{m}}/{\text{s}}^{ 2}$$

Hence, the PGA value calculated in numerical analysis agrees well with the value obtained in the analytical solution.

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 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 non-homogeneous 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.