New pseudo-dynamic analysis of two-layered cohesive-friction soil slope and its numerical validation

Abstract

Natural slopes consist of non-homogeneous soil profiles with distinct characteristics from slopes made of homogeneous soil. In this study, the limit equilibrium modified pseudo-dynamic method is used to analyze the stability of two-layered c-φ soil slopes in which the failure surface is assumed to be a logarithmic spiral. The zero-stress boundary condition at the ground surface under the seismic loading condition is satisfied. New formulations derived from an analytical method are proposed for the predicting the seismic response in two-layered soil. A detailed parametric study was performed in which various parameters (seismic accelerations, damping, cohesion, and angle of internal friction) were varied. The results of the present method were compared with those in the available literature. The present analytical analysis was also verified against the finite element analysis results.

Keywords

layered soil / limit equilibrium method / seismic analysis / damping / PLAXIS

Introduction

In geotechnical earthquake engineering, checking the stability of a slope is very important and the capacity of the slope to sustain the additional load due to earthquake shaking should be ensured during the design of the soil slope. Several researchers have developed various methods to analyze the stability of slopes under earthquake loading conditions. The slope is subjected to an additional force due to the seismic load during an earthquake. Hence, the evaluation of slope stability under seismic loading conditions is an important topic of research. The limit equilibrium method, in which the failure surface is assumed to be circular, has been widely adopted to assess slope stability [¹–³]. This approach has been applied by several researchers to analyze the slope stability in nonlinear failure surfaces [⁴,⁵].

An earthquake has significant effects on the stability of a slope. To introduce the effect of the earthquake on the slope, the pseudo-static approach has been adopted to analyze the slope stability. In this approach, the earthquake loadings are considered as equivalent inertia forces, and the dynamic behavior of the slope is determined on the basis of static equilibrium considerations. The earliest method for analyzing the slope stability under earthquake loading conditions was developed by Terzaghi [⁶]. Terzaghi pointed out that landslides affect the slope stability. Several researchers have adopted the method by assuming circular failure surfaces to analyze the stability of slopes [⁷–¹⁰]. Using the same pseudo-static approach, the slope stability analysis was extended to consider a logarithmic spiral rupture surface [¹¹]. However, the seismic behavior of earthquakes is very roughly approximated in the pseudo-static approach, which may not reproduce the actual seismic effects on the slope. To overcome this drawback, a pseudo-dynamic approach was introduced [¹²] to account for the effects of time and phase differences due to earthquake loading. This approach has been applied by several researchers [¹³–¹⁶] on some common geotechnical structures. However, this pseudo-dynamic approach has some limitations. (a) The method regards the wave as propagating upwards through the soil medium. This does not satisfy the zero-stress boundary condition at the free surface. (b) The method does not consider damping in the soil. (c) This method requires an amplification factor to be introduced. To overcome these limitations, a modified pseudo-dynamic approach was suggested [¹⁷]. This approach was later extended for predicting the active and passive earth pressure on the retaining wall [¹⁸]. A similar method was used to predict the stability of the slope considering logarithmic spiral failure surfaces [¹⁹].

It should be noted that the above-mentioned works investigated geotechnical structures with homogeneous (i.e., single-layered) soil profiles. However, slopes are sometimes made up of two or three different layers depending on the requirements. In fact, both the soil properties and distribution of seismic loads vary non-uniformly with the depth. A non-homogeneous slope stability analysis has been performed [²⁰] using the limit equilibrium framework. A number of researchers have used the Morgenstern and Price method [⁴] and the Spencer method [⁵] to model slopes made up of non-homogeneous soil. An expression was given [²¹] for the safety factor of a c-f soil slope based on limit analysis considering the log-spiral failure mechanism and soil non-homogeneity. A stability analysis of a non-homogeneous slope was performed [²²] for purely cohesive soils exhibiting a linear variation of the cohesion with depth. The stability of the slope was studied [²³] in homogeneous and layered soil considering the effects of rapid drawdown, crack location, and undrained clay soils. A three-dimensional slope stability analysis was presented [²⁴] for anisotropic and non-homogeneous slopes using the upper bound limit analysis approach. A stability chart was proposed [²⁵] for two-layered cohesive slopes using finite element upper and lower bound limit analysis methods. Finite element upper and lower bound limit analysis methods were used [²⁶] to investigate the three-dimensional slope stability of two-layered undrained clay slopes. A series of laboratory tests were conducted [²⁷] on two-layered model slopes subjected to surcharge loads. A stability analysis of a non-homogeneous slope [²⁸] made up of different soil combinations was performed using limit equilibrium and finite element methods.

There are very few studies in the literature [²⁹,³⁰] that analyzed the stability of two-layered slopes using static and pseudo-static methods. As the pseudo-static method provides only very approximate solutions for earthquake-related problems, the pseudo-dynamic method was introduced and was further improved by introducing the effects of significant factor damping and imposing the zero-stress boundary condition at the ground surface. The solution for a two-layered soil slope using this new pseudo-dynamic method has yet to be determined. Therefore, the pseudo-dynamic limit equilibrium method is used to solve two-layered soil slopes in the present study, and the zero-stress boundary condition and soil damping are also considered. The failure surface is considered to be a logarithmic spiral failure surface. Several relevant parameters are optimized through MATLAB, and the variation of these parameters due to seismic waves at different soil layers is also considered in this analysis. To verify the proposed model, the results obtained from the method are compared with those obtained by numerical analysis and earlier studies.

Method of analysis

Assumptions

1) Each soil layer is homogeneous, isotropic, semi-infinite, and dry with constant values of c and f.

2) Both the horizontal and vertical seismic inertial forces are sinusoidal and act at the center of gravity of each strip.

3) The slope fails due to the forces acting on the slope under seismic loading conditions. It is assumed that the logarithmic spiral failure surface generated passes through the toe of the slope.

Problem definition

The failure surface in field conditions is neither linear nor circular, and may be nonlinear. To better model this situation, a more realistic logarithmic spiral failure surface is considered for the analysis [²⁹–³¹]. Let us consider an arbitrary c-f soil slope with two different horizontal layers, as shown in Fig. 1. The top layer is defined by the thickness H₁, cohesion c₁, angle of internal friction f₁, and unit weight g₁. Similarly, the bottom layer is defined by the thickness H₂, cohesion c₂, angle of internal friction f₂, and unit weight g₂.The logarithmic sloping surface is inclined at an angle b to the horizontal.

Collapse Mechanism

The failure surface is assumed to be a combination of logarithmic spiral arcs (shown in Fig. 1) passing through the toe of the slope with the focus O. From the plot, it can be seen that the log-spiral curve portions ODPE and OAME are governed by the respective heights of the soil slopes H₁ and H₂. q₁ and q₂ are the angles, respectively, subtended by the logarithmic spiral curves ODPE and OAME at the center, and q_b is the angle at the center subtended by a line passing through the center and point D of the curve and a horizontal line OO₁ passing through the center.

Problem analysis

The generalized expressions for the logarithmic spiral curve based on the friction angle in the layer are given byr_1=r_0{\mathrm{e}}^{{\theta }_1\tan {\varphi }_1},r_2=r_1{\mathrm{e}}^{{\theta }_2\tan {\varphi }_2},where r_{\mathrm{0}} is the initial radius OD, r_{\mathrm{1}}is the final radius OE at the top layer, andr_{\mathrm{2}}is the final radius OA at the bottom layer.

The initial radius r_{\mathrm{0}}is given by

(3)r_{\mathrm{0}}={\displaystyle \frac{H\cos \mathrm{ec}({\theta }_{\mathrm{a}}+{\theta }_{\mathrm{b}})}{e^{({\theta }_{\mathrm{1}}\tan {\phi }_{\mathrm{1}}+{\theta }_{\mathrm{2}}\tan {\phi }_{\mathrm{2}})}-{\displaystyle \frac{\sin {\theta }_{\mathrm{b}}}{\sin ({\theta }_{\mathrm{a}}+\theta {}_{\mathrm{b}})}}}},

where{\theta }_a={\theta }_{\mathrm{1}}+{\theta }_{\mathrm{2}}.

The width of the soil surface,

(4)CD=b_{\mathrm{s}}={\displaystyle \frac{r_{\mathrm{0}}\sin {\theta }_{\mathrm{a}}}{\sin ({\theta }_{\mathrm{a}}+\theta {}_{\mathrm{b}})}}+H\cot ({\theta }_{\mathrm{a}}+{\theta }_{\mathrm{b}})-H\cot \beta .

Mass of the wedges

For the top layer, the total mass of the potential wedge is calculated by dividing the potential wedge into n strips (Fig. 2).

Therefore, the area of the ith strip in the top layer is given by

(5)A_i={\displaystyle \frac{H_{\mathrm{1}}}{n}}\left[\left\{\begin{array}{l}{\displaystyle \frac{r_{\mathrm{0}}{e}^{(i-\mathrm{1})\frac{{\theta }_{\mathrm{1}}}{n}\tan {\phi }_{\mathrm{1}}}}{\sin \left[{\theta }_{\mathrm{a}}+{\theta }_{\mathrm{b}}+(i-\mathrm{1}){\displaystyle \frac{{\theta }_{\mathrm{1}}}{n}}\right]}}\sin \left[{\theta }_{\mathrm{a}}-(i-\mathrm{1}){\displaystyle \frac{{\theta }_{\mathrm{1}}}{n}}\right]+\\ \left[H-(i-\mathrm{1}){\displaystyle \frac{H_{\mathrm{1}}}{n}}\right]\cot ({\theta }_{\mathrm{a}}+{\theta }_{\mathrm{b}})-\\ \left[H-(i-\mathrm{1}){\displaystyle \frac{H_{\mathrm{1}}}{n}}\right]\cot (\beta )\end{array}\right\}+\left\{\begin{array}{l}{\displaystyle \frac{r_{\mathrm{0}}{e}^{(\mathrm{2}i-\mathrm{1})\frac{{\theta }_{\mathrm{1}}}{n}\tan {\phi }_{\mathrm{1}}}}{\sin \left[{\theta }_{\mathrm{a}}+{\theta }_{\mathrm{b}}+(\mathrm{2}i-\mathrm{1}){\displaystyle \frac{{\theta }_{\mathrm{1}}}{n}}\right]}}\sin \left[{\theta }_{\mathrm{a}}-(\mathrm{2}i-\mathrm{1}){\displaystyle \frac{{\theta }_{\mathrm{1}}}{n}}\right]+\\ \left[H-(\mathrm{2}i-\mathrm{1}){\displaystyle \frac{H_{\mathrm{1}}}{n}}\right]\cot ({\theta }_{\mathrm{a}}+{\theta }_{\mathrm{b}})-\\ \left[H-(\mathrm{2}i-\mathrm{1}){\displaystyle \frac{H_{\mathrm{1}}}{n}}\right]\cot (\beta )\end{array}\right\}\right].

The mass of the potential wedge in the top layer is therefore given by

6)m_I(z_{\mathrm{1}},t)={\displaystyle \frac{{\gamma }_{\mathrm{1}}}{g}}{\displaystyle \sum _{i=\mathrm{1}}^n{A}_i}.

Similarly, the area of the jth strip in the bottom layer is given by

(7)A_j={\displaystyle \frac{H_2}{n}}\left[\left\{\begin{array}{l}{\displaystyle \frac{r_{\mathrm{0}}{\mathrm{e}}^{[{\theta }_{\mathrm{1}}\tan {\phi }_{\mathrm{1}}+(j-\mathrm{1})\frac{{\theta }_2}{n}\tan {\phi }_2]}}{\sin \left[{\theta }_{\mathrm{a}}+{\theta }_{\mathrm{b}}+(j-\mathrm{1}){\displaystyle \frac{{\theta }_2}{n}}\right]}}\sin \left[{\theta }_2-(j-\mathrm{1}){\displaystyle \frac{{\theta }_2}{n}}\right]+\\ \left[H-(j-\mathrm{1}){\displaystyle \frac{H_2}{n}}\right]\cot ({\theta }_{\mathrm{a}}+{\theta }_{\mathrm{b}})-\\ \left[H-(j-\mathrm{1}){\displaystyle \frac{H_2}{n}}\right]\cot (\beta )\end{array}\right\}+\left\{\begin{array}{l}{\displaystyle \frac{r_{\mathrm{0}}{\mathrm{e}}^{[{\theta }_{\mathrm{1}}\tan {\phi }_{\mathrm{1}}+(2j-\mathrm{1})\frac{{\theta }_2}{n}\tan {\phi }_2]}}{\sin \left[{\theta }_{\mathrm{a}}+{\theta }_{\mathrm{b}}+(\mathrm{2}j-\mathrm{1}){\displaystyle \frac{{\theta }_{\mathrm{1}}}{n}}\right]}}\sin \left[{\theta }_2-(\mathrm{2}j-\mathrm{1}){\displaystyle \frac{{\theta }_2}{n}}\right]\\ \left[H-(\mathrm{2}j-\mathrm{1}){\displaystyle \frac{H_{\mathrm{1}}}{n}}\right]\cot ({\theta }_{\mathrm{a}}+{\theta }_{\mathrm{b}})-\\ \left[H-(\mathrm{2}j-\mathrm{1}){\displaystyle \frac{H_{\mathrm{1}}}{n}}\right]\cot (\beta )\end{array}\right\}\right],

where j denotes the strip number on the bottom layer.

Accordingly, the mass of the potential wedge in the bottom layer is given by

(8)m_{\mathrm{II}}(z_{\mathrm{2}},t)={\displaystyle \frac{{\gamma }_{\mathrm{2}}}{g}}{\displaystyle \sum _{j=\mathrm{1}}^n{A}_j}.

Formulation of visco-elastic soil media

Consider a harmonic stress wave traveling through two-layer soil deposits and approaching an interface between the two layers with different soil properties, as shown in Fig. 2. Practical ground response problems usually involve soil deposit layers with different damping characteristics. Because the wave is traveling toward the interface, it is referred to as the incident wave (A_I). When the incident wave reaches the interface, part of its energy is transmitted through the interface and continues to travel through the top layer to form the transmitted wave (A_t). The remaining wave energy is reflected at the interface and travels back through the bottom layer as the reflected wave (A_r). The dynamic response of the soil is described by the behavior of viscoelastic materials, which exhibit both elastic and viscous behavior. In the soil, part of the absorbed energy is converted to the equivalent damping ratio (D), which is a dimensional measure that describes how the oscillation in a system decays after a disturbance. It is assumed that each layer of soil behaves as a Kelvin-Voigt model. This model is represented by a purely viscous damper and a purely elastic spring connected in parallel.

The equation of motion of the Kelvin-Voigt visco-elastic medium in vectorial form [³²] is given by

(9)\rho {\displaystyle \frac{{\partial }^{\mathrm{2}}U}{\partial t^{\mathrm{2}}}}=\left\{(\lambda +G)+({\eta }_{\mathrm{1}}+{\eta }_{\mathrm{s}}){\displaystyle \frac{\partial }{\partial t}}\right\}\mathrm{grad}(\theta )+\left(G+{\eta }_{\mathrm{s}}{\displaystyle \frac{\partial }{\partial t}}\right){\nabla }^{\mathrm{2}}U,

where \rhois the density of the soil material, \lambdais the Lame elastic modulus, G is the shear modulus,{\eta }_{\mathrm{1}}and{\eta }_sare the viscosities of the slope soil, U is the displacement with components u_x,u_y, and u_z, and\theta =div(U).

Considering wave propagation along the z-axis in the Kelvin-Voigt model, the solution of Eq. (9) is given by

(10)\rho {\displaystyle \frac{{\partial }^{\mathrm{2}}{u}_h}{\partial t^{\mathrm{2}}}}=G{\displaystyle \frac{{\partial }^{\mathrm{2}}{u}_h}{\partial z^{\mathrm{2}}}}+{\eta }_s{\displaystyle \frac{{\partial }^{\mathrm{3}}{u}_h}{\partial z^{\mathrm{2}}\partial t}},

(11)\rho {\displaystyle \frac{{\partial }^{\mathrm{2}}{u}_v}{\partial t^{\mathrm{2}}}}=(\lambda +\mathrm{2}G){\displaystyle \frac{{\partial }^{\mathrm{2}}{u}_v}{\partial z^{\mathrm{2}}}}+({\eta }_{\mathrm{1}}+\mathrm{2}{\eta }_s){\displaystyle \frac{{\partial }^{\mathrm{3}}{u}_v}{\partial z^{\mathrm{2}}\partial t}}.

In the present analysis, the base of the slope is considered to be rigid. Therefore, any downward-traveling waves in the soil are completely reflected back toward the ground surface by the rigid base.

Top layer

The equation of motion for harmonic waves may be written as

(12)u_I(z_{\mathrm{1}},t)=A_{\mathrm{t}}{e}^{\stackrel{*}{i}\left(\omega t+\stackrel{*}{k_{\mathrm{s}}}{z}_{\mathrm{1}}\right)}+A_{\mathrm{r}}{e}^{\stackrel{*}{i}\left(\omega t-\stackrel{*}{k_{\mathrm{s}}}{z}_{\mathrm{1}}\right)},

where i* is an imaginary number.

To satisfy the zero-stress boundary condition at the ground surface, A_t = A_r.

The equation of motion for harmonic waves in the top layer thus yields

(13)u_I(z_{\mathrm{1}},t)=A_{\mathrm{t}}\left\{e^{\stackrel{*}{i}\left(\omega t+\stackrel{*}{k_{\mathrm{s}}}{z}_{\mathrm{1}}\right)}+e^{\stackrel{*}{i}\left(\omega t-\stackrel{*}{k_{\mathrm{s}}}{z}_{\mathrm{1}}\right)}\right\}.

The displacement amplitude of the transmitted wave is given by

(14)A_{\mathrm{t}}={\displaystyle \frac{\mathrm{2}}{(\mathrm{1}+{\alpha }_{\mathrm{z}})}}{A}_I,

where the impedance ratio{\alpha }_z={\displaystyle \frac{{\rho }_{\mathrm{2}}{v}_{\mathrm{2}}}{{\rho }_{\mathrm{1}}{v}_{\mathrm{1}}}}. r₁ and r₂ are the densities of the top and bottom soil layers, respectively, and v_s1 and v_s2 the corresponding shear wave velocities.

The complex wave number can be written as

(15)\stackrel{*}{k_{\mathrm{s}}}=k_{\mathrm{s1}}+i{k}_{\mathrm{s2}},

where k_{\mathrm{s1}}={\displaystyle \frac{{\omega }_{\mathrm{s1}}}{v_{\mathrm{s1}}}}{\left\{{\displaystyle \frac{{(\mathrm{1}+\mathrm{4}{D}_{\mathrm{s1}}^{\mathrm{2}})}^{\raisebox{1ex}{$\mathrm{1}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{2}$}\right.}+\mathrm{1}}{\mathrm{2}{(\mathrm{1}+\mathrm{4}{D}_{\mathrm{s1}}^{\mathrm{2}})}^{\raisebox{1ex}{$\mathrm{1}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{2}$}\right.}}}\right\}}^{\raisebox{1ex}{$\mathrm{1}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{2}$}\right.}and k_{\mathrm{s2}}={\displaystyle \frac{{\omega }_{\mathrm{s1}}}{v_{\mathrm{s1}}}}{\left\{{\displaystyle \frac{{(\mathrm{1}+\mathrm{4}{D}_{\mathrm{s1}}^{\mathrm{2}})}^{\raisebox{1ex}{$\mathrm{1}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{2}$}\right.}-\mathrm{1}}{\mathrm{2}{(\mathrm{1}+\mathrm{4}{D}_{\mathrm{s1}}^{\mathrm{2}})}^{\raisebox{1ex}{$\mathrm{1}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{2}$}\right.}}}\right\}}^{\raisebox{1ex}{$\mathrm{1}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{2}$}\right.}.

Considering the damping as a function of D_{\mathrm{s}}={\displaystyle \frac{{\eta }_{\mathrm{s}}{\omega }_{\mathrm{s}}}{\mathrm{2}G}}, the real and imaginary parts of the horizontal displacement at any depth z₁ can be obtained as

(16)u_{{\mathrm{h}}_{\mathrm{I}}}(z_{\mathrm{1}},t)=A_{\mathrm{t}}\{\cosh (k_{\mathrm{s2}}{z}_{\mathrm{1}})+\sinh (k_{\mathrm{s2}}{z}_{\mathrm{1}})\}\left\{\cos \left(\mathrm{2}\pi {\displaystyle \frac{t}{T}}-k_{\mathrm{s1}}{z}_{\mathrm{1}}\right)-\cos \left(\mathrm{2}\pi {\displaystyle \frac{t}{T}}+k_{\mathrm{s1}}{z}_{\mathrm{1}}\right)\right\},(\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{l}\mathrm{p}\mathrm{a}\mathrm{r}\mathrm{t})

(17)u_{{\mathrm{h}}_{\mathrm{II}}}(z_{\mathrm{1}},t)=A_{\mathrm{t}}\{\cosh (k_{\mathrm{s2}}{z}_{\mathrm{1}})+\sinh (k_{\mathrm{s2}}z)\}\left\{\sin \left(\mathrm{2}\pi {\displaystyle \frac{t}{T}}-k_{\mathrm{s1}}{z}_{\mathrm{1}}\right)-\sin \left(\mathrm{2}\pi {\displaystyle \frac{t}{T}}+k_{\mathrm{s1}}{z}_{\mathrm{1}}\right)\right\}.(\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{r}\mathrm{y}\mathrm{p}\mathrm{a}\mathrm{r}\mathrm{t})

Similarly, considering the damping as a function ofD_{\mathrm{p}}={\displaystyle \frac{{\eta }_{\mathrm{p}}{\omega }_{\mathrm{p}}}{\mathrm{2}(\lambda +\mathrm{2}G)}} where{\eta }_{\mathrm{p}}=({\eta }_{\mathrm{1}}+\mathrm{2}{\eta }_s), the real and imaginary parts of the vertical displacement can be expressed as

(18)u_{{\mathrm{v}}_{\mathrm{I}}}(z_{\mathrm{1}},t)=A_t\{\cosh (k_{\mathrm{p2}}{z}_{\mathrm{1}})+\sinh (k_{\mathrm{p2}}{z}_{\mathrm{1}})\}\left\{\cos \left(\mathrm{2}\pi {\displaystyle \frac{t}{T}}-k_{\mathrm{p1}}{z}_{\mathrm{1}}\right)-\cos \left(\mathrm{2}\pi {\displaystyle \frac{t}{T}}+k_{\mathrm{p1}}{z}_{\mathrm{1}}\right)\right\},(\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{l}\mathrm{p}\mathrm{a}\mathrm{r}\mathrm{t})

(19)u_{{\mathrm{v}}_{\mathrm{II}}}(z_{\mathrm{1}},t)=A_t\{\cosh (k_{\mathrm{p2}}{z}_{\mathrm{1}})+\sinh (k_{\mathrm{p2}}z)\}\left\{\sin \left(\mathrm{2}\pi {\displaystyle \frac{t}{T}}-k_{\mathrm{p1}}{z}_{\mathrm{1}}\right)-\sin \left(\mathrm{2}\pi {\displaystyle \frac{t}{T}}+k_{\mathrm{p1}}{z}_{\mathrm{1}}\right)\right\}.(\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{r}\mathrm{y}\mathrm{p}\mathrm{a}\mathrm{r}\mathrm{t})

In this analysis, only the real part is considered. Differentiating Eq. (18) twice, the horizontal acceleration in the top layer at any depth z_i can be obtained as

(20)a_{{\mathrm{h}}_{\mathrm{I}}}(z_i,t)=k_{\mathrm{h}}g{A}_{\mathrm{t}}\left[\begin{array}{l}\mathrm{2}{k}_{\mathrm{s1}}{k}_{\mathrm{s2}}[\{\cosh (k_{\mathrm{s2}}{z}_i)-\sinh (k_{\mathrm{s2}}{z}_i)\}\sin (\omega t+k_{\mathrm{s1}}{z}_i)]+\\ \left[\begin{array}{l}\{\cosh (k_{\mathrm{s2}}{z}_i)+\sinh (k_{\mathrm{s2}}{z}_i)\}\\ \sin (\omega t-k_{\mathrm{1}}{z}_i)\end{array}\right]-(k_{\mathrm{s1}}^{\mathrm{2}}-k_{\mathrm{s2}}^{\mathrm{2}})\left[\begin{array}{l}\{\cosh (k_{\mathrm{s2}}{z}_i)-\sinh (k_{\mathrm{s2}}{z}_i)\}\\ \cos (\omega t+k_{\mathrm{s1}}{z}_i)\end{array}\right]\\ +[\{\cosh (k_{\mathrm{s2}}{z}_i)+\sinh (k_{\mathrm{s2}}{z}_i)\}\cos (\omega t-k_{\mathrm{s1}}{z}_i)]\end{array}\right].

Accordingly, the vertical acceleration in the top layer is expressed as

(21)a_{{\mathrm{v}}_{\mathrm{I}}}(z_i,t)=k_{\mathrm{v}}g{A}_{\mathrm{t}}\left[\begin{array}{l}\mathrm{2}{k}_{\mathrm{p1}}{k}_{\mathrm{p2}}[\{\cosh (k_{\mathrm{p2}}{z}_i)-\sinh (k_{\mathrm{p2}}{z}_i)\}\sin (\omega t+k_{\mathrm{p1}}{z}_i)]+\\ \left[\begin{array}{l}\{\cosh (k_{\mathrm{p2}}{z}_i)+\sinh (k_{\mathrm{p2}}{z}_i)\}\\ \sin (\omega t-k_{\mathrm{1}}{z}_i)\end{array}\right]-(k_{s\mathrm{1}}^{\mathrm{2}}-k_{\mathrm{p2}}^{\mathrm{2}})\left[\begin{array}{l}\{\cosh (k_{\mathrm{p2}}{z}_i)-\sinh (k_{\mathrm{p2}}{z}_i)\}\\ \cos (\omega t+k_{\mathrm{p1}}{z}_i)\end{array}\right]\\ +[\{\cosh (k_{\mathrm{p2}}{z}_i)+\sinh (k_{\mathrm{p2}}{z}_i)\}\cos (\omega t-k_{\mathrm{p1}}{z}_i)]\end{array}\right],

wherez_i=\left[\left({\displaystyle \frac{\mathrm{2}i-\mathrm{1}}{\mathrm{2}}}\right){\displaystyle \frac{H_{\mathrm{1}}}{n}}\right].

Bottom layer

Note that the displacements at the layer boundaries must be compatible, that is, the displacement at the top of the bottom layer must be equal to the displacement at the bottom of the top layer. Applying the compatibility requirement at the interface, the real and imaginary parts of the horizontal displacement in the bottom layer can be expressed as

(22)u_{{\mathrm{h}}_{\mathrm{II}}}={\displaystyle \frac{A_{\mathrm{t}}}{\mathrm{2}}}[(\mathrm{1}+{\alpha }_{\mathrm{z}})I_{\mathrm{1}}+(\mathrm{1}-{\alpha }_{\mathrm{z}})I_{\mathrm{2}}],(\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{l}\mathrm{p}\mathrm{a}\mathrm{r}\mathrm{t})

(23)u_{{\mathrm{h}}_{\mathrm{II}}}={\displaystyle \frac{A_{\mathrm{t}}}{\mathrm{2}}}[(\mathrm{1}+{\alpha }_{\mathrm{z}})I_{\mathrm{3}}+(\mathrm{1}-{\alpha }_{\mathrm{z}})I_{\mathrm{4}}].(\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{r}\mathrm{y}\mathrm{p}\mathrm{a}\mathrm{r}\mathrm{t})

The detailed expressions of I₁, I₂, I₃, and I₄ are given in the appendix.

The horizontal acceleration in the bottom layer can be obtained as

(24)a_{{\mathrm{h}}_{\mathrm{II}}}(z_j,t)=k_{\mathrm{h}}g{A}_{\mathrm{t}}(\mathrm{1}+{\alpha }_{\mathrm{z}})\left[\begin{array}{l}[\{(k_{\mathrm{s1}}^{\mathrm{2}}-k_{\mathrm{s2}}^{\mathrm{2}})(y_{\mathrm{s2}}-y_{\mathrm{s4}})-\mathrm{2}{k}_{\mathrm{s1}}{k}_{\mathrm{s2}}(y_{\mathrm{s1}}-y_{\mathrm{s3}})\}\{\sinh (k_{\mathrm{s2}}{z}_j)-\cosh (k_{\mathrm{s2}}{z}_j)\}]\sin (\omega t+k_{\mathrm{s1}}{z}_j)\\ -[\{(k_{\mathrm{s1}}^{\mathrm{2}}-k_{\mathrm{s2}}^{\mathrm{2}})(y_{\mathrm{s2}}+y_{\mathrm{s4}})-\mathrm{2}{k}_{\mathrm{s1}}{k}_{\mathrm{s2}}(y_{\mathrm{s1}}+y_{\mathrm{s3}})\}\{\sinh (k_{\mathrm{s2}}{z}_j)+\cosh (k_{\mathrm{s2}}{z}_j)\}]\sin (\omega t-k_{\mathrm{s1}}{z}_j)\\ +[\{(k_{\mathrm{s1}}^{\mathrm{2}}-k_{\mathrm{s2}}^{\mathrm{2}})(y_{\mathrm{s1}}-y_{\mathrm{s3}})+\mathrm{2}{k}_{\mathrm{s1}}{k}_{\mathrm{s2}}(y_{\mathrm{s2}}-y_{\mathrm{s4}})\}\{\sinh (k_{\mathrm{s2}}{z}_j)-\cosh (k_{\mathrm{s2}}{z}_j)\}]\cos (\omega t+k_{\mathrm{s1}}{z}_j)\\ -[\{(k_{\mathrm{s1}}^{\mathrm{2}}-k_{\mathrm{s2}}^{\mathrm{2}})(y_{\mathrm{s1}}+y_{\mathrm{s3}})+\mathrm{2}{k}_{\mathrm{s1}}{k}_{\mathrm{s2}}(y_{\mathrm{s2}}+y_{\mathrm{s4}})\}\{\sinh (k_{\mathrm{s2}}{z}_j)+\cosh (k_{\mathrm{s2}}{z}_j)\}]\cos (\omega t-k_{\mathrm{s1}}{z}_j)\end{array}\right],

wherey_{\mathrm{s1}}=\cos (k_{\mathrm{s1}}{H}_{\mathrm{1}})\cosh (k_{\mathrm{s1}}{H}_{\mathrm{1}}),

y_{\mathrm{s2}}=\sin (k_{\mathrm{s1}}{H}_{\mathrm{1}})\sinh (k_{\mathrm{s1}}{H}_{\mathrm{1}}),

y_{\mathrm{s}3}=\cos (k_{\mathrm{s1}}{H}_{\mathrm{1}})\sinh (k_{\mathrm{s1}}{H}_{\mathrm{1}}),

y_{\mathrm{s}4}=\sin (k_{\mathrm{s1}}{H}_{\mathrm{1}})\sinh (k_{\mathrm{s1}}{H}_{\mathrm{1}}),

Similarly, the vertical acceleration in layer II at any depth z_j can be expressed as

a_{{\mathrm{v}}_{\mathrm{II}}}(z_j,t)=k_{\mathrm{v}}g{A}_{\mathrm{t}}(\mathrm{1}+{\alpha }_z)\left[\begin{array}{l}[\{(k_{\mathrm{p1}}^{\mathrm{2}}-k_{\mathrm{p2}}^{\mathrm{2}})(y_{\mathrm{p2}}-y_{\mathrm{p4}})-\mathrm{2}{k}_{\mathrm{p1}}{k}_{\mathrm{p2}}(y_{\mathrm{p1}}-y_{\mathrm{p3}})\}\{\sinh (k_{\mathrm{p2}}{z}_j)-\cosh (k_{\mathrm{p2}}{z}_j)\}]\sin (\omega t+k_{\mathrm{p1}}{z}_j)\\ -[\{(k_{\mathrm{p1}}^{\mathrm{2}}-k_{\mathrm{p2}}^{\mathrm{2}})(y_{\mathrm{p2}}+y_{\mathrm{p4}})-\mathrm{2}{k}_{\mathrm{p1}}{k}_{\mathrm{p2}}(y_{\mathrm{p1}}+y_{\mathrm{p3}})\}\{\sinh (k_{\mathrm{p2}}{z}_j)+\cosh (k_{\mathrm{p2}}{z}_j)\}]\sin (\omega t-k_{\mathrm{p1}}{z}_j)\\ +[\{(k_{\mathrm{p1}}^{\mathrm{2}}-k_{\mathrm{p2}}^{\mathrm{2}})(y_{\mathrm{p1}}-y_{\mathrm{p3}})+\mathrm{2}{k}_{\mathrm{p1}}{k}_{\mathrm{p2}}(y_{\mathrm{p2}}-y_{\mathrm{p4}})\}\{\sinh (k_{\mathrm{p2}}{z}_j)-\cosh (k_{\mathrm{p2}}{z}_j)\}]\cos (\omega t+k_{\mathrm{p1}}{z}_j)\\ -[\{(k_{\mathrm{p1}}^{\mathrm{2}}-k_{\mathrm{p2}}^{\mathrm{2}})(y_{\mathrm{p1}}+y_{\mathrm{p3}})+\mathrm{2}{k}_{\mathrm{p1}}{k}_{\mathrm{p2}}(y_{\mathrm{p2}}+y_{\mathrm{p4}})\}\{\sinh (k_{\mathrm{p2}}{z}_j)+\cosh (k_{\mathrm{p2}}{z}_j)\}]\cos (\omega t-k_{\mathrm{p1}}{z}_j)\end{array}\right]

+A_{\mathrm{t}}(\mathrm{1}-{\alpha }_{\mathrm{z}})\left[\begin{array}{l}[\{(k_{\mathrm{p1}}^{\mathrm{2}}-k_{\mathrm{p2}}^{\mathrm{2}})(y_{\mathrm{p2}}+y_{\mathrm{p4}})-\mathrm{2}{k}_{\mathrm{p1}}{k}_{\mathrm{p2}}(y_{\mathrm{p1}}+y_{\mathrm{p3}})\}\{\sinh (k_{\mathrm{p2}}{z}_j)-\cosh (k_{\mathrm{p2}}{z}_j)\}]\sin (\omega t+k_{\mathrm{p1}}{z}_j)\\ -[\{(k_{\mathrm{p1}}^{\mathrm{2}}-k_{\mathrm{p2}}^{\mathrm{2}})(y_{\mathrm{p2}}-y_{\mathrm{p4}})-\mathrm{2}{k}_{\mathrm{p1}}{k}_{\mathrm{p2}}(y_{\mathrm{p1}}-y_{\mathrm{p3}})\}\{\sinh (k_{\mathrm{p2}}{z}_j)+\cosh (k_{\mathrm{p2}}{z}_j)\}]\sin (\omega t-k_{\mathrm{p1}}{z}_j)\\ +[\{(k_{\mathrm{p1}}^{\mathrm{2}}-k_{\mathrm{p2}}^{\mathrm{2}})(y_{\mathrm{p1}}+y_{\mathrm{p3}})+\mathrm{2}{k}_{\mathrm{p1}}{k}_{\mathrm{p2}}(y_{\mathrm{p2}}+y_{\mathrm{p4}})\}\{\sinh (k_{\mathrm{p2}}{z}_j)-\cosh (k_{\mathrm{p2}}{z}_j)\}]\cos (\omega t+k_{\mathrm{p1}}{z}_j)\\ -[\{(k_{\mathrm{p1}}^{\mathrm{2}}-k_{\mathrm{p2}}^{\mathrm{2}})(y_{\mathrm{p1}}-y_{\mathrm{p3}})+\mathrm{2}{k}_{\mathrm{p1}}{k}_{\mathrm{p2}}(y_{\mathrm{p2}}-y_{\mathrm{p4}})\}\{\sinh (k_{\mathrm{p2}}{z}_j)+\cosh (k_{\mathrm{p2}}{z}_j)\}]\cos (\omega t-k_{\mathrm{p1}}{z}_j)\end{array}\right]

where z_j=H_{\mathrm{1}}+\left[\left({\displaystyle \frac{\mathrm{2}j-\mathrm{1}}{\mathrm{2}}}\right){\displaystyle \frac{H_{\mathrm{2}}}{n}}\right],

y_{\mathrm{p1}}=\cos (k_{\mathrm{p1}}{H}_{\mathrm{1}})\cosh (k_{\mathrm{p1}}{H}_{\mathrm{1}}),

y_{\mathrm{p2}}=\sin (k_{\mathrm{p1}}{H}_{\mathrm{1}})\sinh (k_{\mathrm{p1}}{H}_{\mathrm{1}}),

y_{\mathrm{p3}}=\cos (k_{\mathrm{p1}}{H}_{\mathrm{1}})\sinh (k_{\mathrm{p1}}{H}_{\mathrm{1}}),

y_{\mathrm{p4}}=\sin (k_{\mathrm{p1}}{H}_{\mathrm{1}})\cosh (k_{\mathrm{p1}}{H}_{\mathrm{1}}).

Calculation of horizontal and vertical inertia forces

The horizontal inertia force Q_{{\mathrm{h}}_{\mathrm{I}}}acting on the top layer can be expressed as

(26)Q_{{\mathrm{h}}_{\mathrm{I}}}={\gamma }_{\mathrm{1}}{k}_{\mathrm{h}}{\displaystyle \sum _{i=\mathrm{1}}^n{m}_{\mathrm{I}}(z_i)a_{h_I}(z_i,t)}.

On the other hand, the horizontal inertia force Q_{h_{\mathrm{II}}}acting on the bottom layer is given by

(27)Q_{h_{\mathrm{II}}}={\gamma }_{\mathrm{2}}{k}_{\mathrm{h}}{\displaystyle \sum _{j=\mathrm{1}}^n{m}_{\mathrm{II}}(z_j)a_{h_{\mathrm{II}}}(z_j,t)}.

Therefore, the total horizontal inertia force Q_{\mathrm{h}}acting on the failure wedge is given by

(28)Q_{\mathrm{h}}=k_{\mathrm{h}}\left\{{\gamma }_{\mathrm{1}}{\displaystyle \sum _{i=\mathrm{1}}^n{m}_{\mathrm{I}}(z_i)a_{{\mathrm{h}}_{\mathrm{I}}}(z_i,t)+{\gamma }_{\mathrm{2}}{\displaystyle \sum _{j=\mathrm{1}}^n{m}_{\mathrm{II}}(z_j)a_{{\mathrm{h}}_{\mathrm{II}}}(z_j,t)}}\right\}.

Similarly, the total vertical inertia force Q_{\mathrm{v}}acting on the failure wedge is given by

(29)Q_{\mathrm{v}}=k_{\mathrm{v}}\left\{{\gamma }_{\mathrm{1}}{\displaystyle \sum _{i=\mathrm{1}}^n{m}_{\mathrm{I}}(z_i)a_{{\mathrm{v}}_{\mathrm{I}}}(z_i,t)+{\gamma }_{\mathrm{2}}{\displaystyle \sum _{j=\mathrm{1}}^n{m}_{\mathrm{II}}(z_j)a_{{\mathrm{v}}_{\mathrm{II}}}(z_j,t)}}\right\}.

The disturbing moments at the top and bottom layers due to horizontal forces are given by

(30){(M_{\mathrm{h}})}_{\mathrm{1}}=Q_{{\mathrm{h}}_{\mathrm{I}}}\times {\overline{y}}_{\mathrm{1}},

(31){(M_{\mathrm{h}})}_{\mathrm{2}}=Q_{{\mathrm{h}}_{\mathrm{II}}}\times {\overline{y}}_{\mathrm{2}}.

The disturbing moments at the top and bottom layers due to vertical forces are given by

(32){(M_{\mathrm{v}})}_{\mathrm{1}}=Q_{{\mathrm{v}}_I}\times {\overline{x}}_{\mathrm{1}},

(33){(M_{\mathrm{v}})}_{\mathrm{2}}=Q_{{\mathrm{v}}_{\mathrm{II}}}\times {\overline{x}}_{\mathrm{2}},

where {\bar{x}}_{\mathrm{1}},{\bar{y}}_{\mathrm{1}}and {\bar{x}}_{\mathrm{2}},{\bar{y}}_{\mathrm{2}}are respectively the horizontal and vertical distances of the centers of gravity from the centers of the logarithmic arcs for the two layers, and can be expressed as

(34){\overline{x}}_{\mathrm{1}}={\displaystyle \frac{M_{\mathrm{ODPE}}+M_{{\mathrm{BEE}}_1}-M_{{\mathrm{ODD}}_1}-M_{{\mathrm{CD}}_1{\mathrm{E}}_1}}{w_{\mathrm{ODPE}}+w_{{\mathrm{BEE}}_1}-w_{{\mathrm{ODD}}_1}-w_{{\mathrm{CD}}_1{\mathrm{E}}_1}}},

(35){\overline{y}}_{\mathrm{1}}={\displaystyle \frac{{M^{\prime }}_{\mathrm{ODPE}}+{M^{\prime }}_{{\mathrm{BEE}}_1}-{M^{\prime }}_{{\mathrm{ODD}}_1}-M_{{\mathrm{CD}}_1{\mathrm{E}}_1}}{w_{\mathrm{ODPE}}+w_{{\mathrm{BEE}}_1}-w_{{\mathrm{ODD}}_1}-w_{{\mathrm{CD}}_1{\mathrm{E}}_1}}},

(36){\overline{x}}_{\mathrm{2}}={\displaystyle \frac{M_{\mathrm{OAME}}-M_{{\mathrm{BEE}}_1}-M_{\mathrm{ACJ}}-M_{{\mathrm{OD}}_1\mathrm{J}}}{w_{\mathrm{OAME}}-w_{{\mathrm{BEE}}_1}-w_{\mathrm{ACJ}}-w_{{\mathrm{OD}}_1\mathrm{J}}}},

(37){\overline{y}}_{\mathrm{2}}={\displaystyle \frac{{M^{\prime }}_{\mathrm{OAME}}-{M^{\prime }}_{{\mathrm{BEE}}_1}-{M^{\prime }}_{\mathrm{ACJ}}-{M^{\prime }}_{{\mathrm{OD}}_1\mathrm{J}}}{w_{\mathrm{OAME}}-w_{{\mathrm{BEE}}_1}-w_{\mathrm{ACJ}}-w_{{\mathrm{OD}}_1\mathrm{J}}}},

w_{\mathrm{ODPE}},w_{{\mathrm{BEE}}_1},w_{{\mathrm{ODD}}_1},w_{{\mathrm{CD}}_1{\mathrm{E}}_1},w_{\mathrm{OAME}},w_{\mathrm{ACJ}}, and w_{{\mathrm{OD}}_1\mathrm{J}} are the weights of portions ODPE, BEE₁, ODD₁, CD₁E₁, OAME, ACJ, and OD₁J, respectively. M_{\mathrm{ODPE}},M_{{\mathrm{BEE}}_1},M_{{\mathrm{ODD}}_1},M_{{\mathrm{CD}}_1{\mathrm{E}}_1},M_{\mathrm{OAME}},M_{\mathrm{ACJ}},M_{{\mathrm{OD}}_1\mathrm{J}},{M^{\prime }}_{\mathrm{ODPE}},{M^{\prime }}_{{\mathrm{BEE}}_1},{M^{\prime }}_{{\mathrm{ODD}}_1},{M^{\prime }}_{{\mathrm{CD}}_1{\mathrm{E}}_1},{M^{\prime }}_{\mathrm{OAME}},{M^{\prime }}_{\mathrm{ACJ}}, and {M^{\prime }}_{{\mathrm{OD}}_1\mathrm{J}} are the moments of portions ODPE, BEE₁, ODD₁, OAME, ACJ, and OD₁J, respectively. The detailed expressions of these terms are given in the appendix.

The moment of force due to the cohesive force acting on the top layer is given by

(38){(M_{{\mathrm{c}}_{\mathrm{m}}})}_I={\displaystyle \frac{c_{{\mathrm{m}}_1}}{\mathrm{2}\tan {\phi }_{{\mathrm{m}}_1}}}(r_{\mathrm{1}}^{\mathrm{2}}-r_{\mathrm{0}}^{\mathrm{2}}).

The moment of force due to the cohesive force acting on the bottom layer is given by

(39){(M_{{\mathrm{c}}_{\mathrm{m}}})}_{\mathrm{II}}={\displaystyle \frac{c_{{\mathrm{m}}_{2_{}}}}{\mathrm{2}\tan {\phi }_{{\mathrm{m}}_2}}}(r_{\mathrm{2}}^{\mathrm{2}}-r_{\mathrm{1}}^{\mathrm{2}}).

The moment of force due to the reaction force acting on the top layer is given by

(40){(M_{{\mathrm{R}}_{\mathrm{m}}})}_{\mathrm{I}}={\displaystyle \frac{w_{\mathrm{BCDEB}}{r}_{\mathrm{0}}\sin {\phi }_{\mathrm{m}}}{\mathrm{1}+{\tan }^{\mathrm{2}}{\phi }_{\mathrm{m}}}}\left[\begin{array}{l}\{(\cos ({\theta }_{\mathrm{b}}+{\theta }_{\mathrm{1}})-\sin ({\theta }_{\mathrm{b}}+{\theta }_{\mathrm{1}})\tan {\phi }_{\mathrm{m}})\}-e^{{\theta }_{\mathrm{2}}\tan {\phi }_{\mathrm{m}}}\\ \{\cos ({\theta }_{\mathrm{b}}+{\theta }_{\mathrm{1}}+{\theta }_{\mathrm{2}})-\sin ({\theta }_{\mathrm{b}}+{\theta }_{\mathrm{1}}+{\theta }_{\mathrm{2}})\tan {\phi }_{\mathrm{m}}\}\end{array}\right].

The moment of force due to the reaction force acting on the bottom layer is given by

(41){(M_{{\mathrm{R}}_{\mathrm{m}}})}_{\mathrm{II}}={\displaystyle \frac{w_{\mathrm{ABEMA}}{r}_{\mathrm{0}}\sin {\phi }_{\mathrm{m}}}{\mathrm{1}+{\tan }^{\mathrm{2}}{\phi }_{\mathrm{m}}}}[(\cos {\theta }_{\mathrm{b}}-\sin {\theta }_{\mathrm{b}}\tan {\phi }_{\mathrm{m}})-e^{{\theta }_{\mathrm{1}}\tan {\phi }_{\mathrm{m}}}\{\cos ({\theta }_{\mathrm{b}}+{\theta }_{\mathrm{1}})-\sin ({\theta }_{\mathrm{b}}+{\theta }_{\mathrm{1}})\tan {\phi }_{\mathrm{m}}\}],

where m denotes the mobilization factor.

The factor of safety (FOS) of the slope is thus

(42)FOS={\displaystyle \frac{{(M_{{\mathrm{c}}_{\mathrm{m}}})}_{\mathrm{I}}+{(M_{{\mathrm{R}}_{\mathrm{m}}})}_{\mathrm{I}}+{(M_{{\mathrm{c}}_{\mathrm{m}}})}_{\mathrm{II}}+{(M_{{\mathrm{R}}_{\mathrm{m}}})}_{\mathrm{II}}}{{(M_{\mathrm{h}}+M_{\mathrm{v}})}_{\mathrm{I}}+{(M_{\mathrm{h}}+M_{\mathrm{v}})}_{\mathrm{II}}}}.

Results and discussions

From Eq. (42), it can be seen that for a particular layered slope under a particular seismic loading condition, all the variables are constant except for q₁, q₂, q_b, m, and t/T. To locate the critical logarithmic spiral failure surfaces, an optimization technique was implemented in MATLAB [³³] to obtain the corresponding values of the FOS with respect to q₁, q₂, q_b, m, and t/T. The details of the optimization algorithm are shown in Fig. 3. The results of the analytical solution are given in tabular form in Table 1.

For the two-layered soil slopes, the following typical parameters are used: g₁ = 16 kN/m³; g₂ = 16 kN/m³; f₁ = 30°; f₂ = 18°; c₁ = 10 kN/m³; c₂ = 25 kN/m³; k_h = 0.1, 0.2, and 0.3; k_v = 0, k_h/2, and k_h; b = 10° –50°; D₁/D₂ = 0.25, 0.5 and 1.0; A_i = 100 kN/m² and v_p/v_s = 1.87 for m = 0.3.

For simplicity, it is assumed that D= D_s= D_p. The influence of the damping ratio on the slope stability is demonstrated in Fig. 4. Here, D₁ is the damping ratio of the bottom layer and D₂ is the damping ratio of the bottom layer. Keeping the damping ratio of the bottom layer constant, the damping ratio of the top layer is increased such that the ratio D₁/D₂ increases. So, it is obvious that when the damping ratio increases, the FOS will also increase. For example, at f₂/f₁ = 0.6, b = 30°, k_h = 0.1, k_v = k_h/2, c₁/c₂ = 0.4, g₁/g₂ = 1.0, and H₁/H₂ = 0.6, the FOS increases by 32.9% when D increases from 10% to 20%.

The effects of the horizontal and vertical seismic accelerations on the slope stability are plotted in Fig. 5 and 6, respectively. The results indicate that the FOS decreases with the increase of both the horizontal and vertical seismic accelerations. It is noted that the effect of horizontal seismic accelerationsis more prominant than that of vertical seismic accelerations. It is well established that the seismic accelerations increase with the earthquake magnitudes, which in turn increase the shear strains or the amplification of the soil layers, and hence decrease the FOS. For example, at f₂/f₁ = 0.6, b = 30°, D₁/D₂ = 0.5, k_v= k_h/2, c₁/c₂ = 0.4, g₁/g₂ = 1.0, and H₁/H₂ = 0.6, the FOS decreases by 16.23% at k_h = 0.1 when k_h increases from 0.1 to 0.3. Similarly, for the same slope and k_h = 0.1, if k_v increases from k_h/2 to k_h, the FOS decreases by approximately 5.89% at k_v = k_h/2.

As discussed above, it is also observed that with the increase of the slope angle, the FOS decreases because of the increase in the driving force in the slope. The slope may become unstable if the FOS<1.

Figure 7 shows the variation in the FOS for different ratios of c₁/c₂at b = 30°, f₂/f₁ = 0.6, k_h = 0.1, k_v = k_h/2, D₁/D₂ = 0.5, g₁/g₂ = 1.0, and H₁/H₂ = 0.6. It is observed that the FOS increases with the ratio of c₁/c₂ because the cohesion increases the attraction between the soil particles, which in turn increases the shear strength of the soil. For example, the FOS increases by approximately 11% and 19.5% when the ratio of c₁/c₂ increases from 0.4 to 0.6 and from 0.4 to 1.0, respectively, for c₁/c₂ = 0.4. Here, the ratio of c₁/c₂ is increased by keeping the value of c₂ constant.

Figure 8 presents the influence of f₂/f₁ on the slope stability. Here, the ratio of f₂/f₁ is increased by keeping f₁ constant. It can be seen that the FOS increases with an increase in the ratio of f₂/f₁. This clearly suggests that the slope stability is significantly affected by the soil friction angle and that a better compacting effect in the slope leads to a higher soil friction angle. For example, the FOS increases by approximately 34% when f₂/f₁ increases from 0.6 to 0.8 atD₁/D₂ = 0.5, b = 30°, k_h = 0.1, k_v = k_h/2, c₁/c₂ = 0.4, g₁/g₂ = 1.0, and H₁/H₂ = 0.6.

Figure 9 shows the variation in the FOS with respect to the slope angle for different ratios of g₁/g₂ at D₁/D₂ = 0.5, b = 30°, k_h = 0.1, k_v = k_h/2, c₁/c₂ = 0.4, f₂/f₁ = 0.6, and H₁/H₂ = 0.6. From the plot, it is observed that the FOS increases with the g₁/g₂ ratio. For example, the FOS increases by approximately 8.5% when the g₁/g₂ ratio varies from 0.6 to1.0. An increase in the g₁/g₂ ratio implies the tightening of the soil particles, which results in an increase in the shear strength of the soil and thus, the FOS. Here, the ratio of g₁/g₂ is increased by keeping g₂ constant.

Figure 10 depicts the variation in the FOS with the H₁/H₂ ratio at D₁/D₂ = 0.5, b = 30°, k_h = 0.1, k_v = k_h/2, c₁/c₂ = 0.4, f₂/f₁ = 0.6, and g₁/g₂ = 1.0. Here, H₁ is increased while H₂ is kept constant. From the plot, it is noticed that the FOS increases with the ratio of H₁/H₂. For example, the FOS decreases by approximately 12.5% when the H₁/H₂ ratio changes from 0.6 to1.0.The increase in the H₁/H₂ ratio decreases the self-weight of the slope, which may decrease the driving forces acting on the slope and hence increase the FOS.

Numerical analysis

To perform the analysis, a two-dimensional finite element model of the theoretical slope model during an earthquake was developed and analyzed using PLAXIS 2D [³⁴]. Slope stability analysis using the finite element method has been reported by several researchers [³⁵–³⁹]. For the standard boundaries of the slope, the base was set at a depth 13m below the top level of the slope and rollers were placed on the two vertical sides; one vertical side passes through the center line of the slope and the other at a distance of 25 m from the toe of the slope. The rollers are such that they only allow vertical movement of the slope material. A full fixed boundary condition was considered at the base of the slope. The slope within these boundaries was divided into a number of 6-node triangular elements. The geometry of the finite element slope model and the mesh arrangement are depicted in Fig. 11. A sensitivity analysis was performed to understand how the size of the mesh affected the results obtained by the numerical analysis. Figure 12 shows that with the decrease in the element size, the FOS of the slope also decreased up to an element size of 0.5 m. When the element size was further reduced to less than 0.5 m, no considerable FOS reduction was observed. Therefore, a mesh size of 0.5 m was used in the present numerical analysis and the soils were discretized with 378 triangular elements.

The numerical study considered the propagated shear and primary waves transmitted through the soil medium. The maximum size of the finite element mesh also satisfies the condition that the size of the element should not be greater than 20% of the wavelength of a shear wave [⁴⁰]. According to Bakr and Ahmad [⁴¹], if l_min and v_s denote the wavelength and shear wave velocity, respectively, and f_max is the maximum frequency of the seismic input motion, the height of the finite element mesh should be less than or equal to he_max:

(43)h{e}_{\max }={\displaystyle \frac{{\lambda }_{\min }}{\mathrm{5}}}={\displaystyle \frac{v_{\mathrm{s}}}{\mathrm{5}{f}_{\max }}}.

To ascertain the generation of the free failure surface without any obstruction, the boundary of the model was extended beyond than the size of the actual model in both the horizontal and vertical directions, that is, along the x and y-axes. The slope was initially considered to be at rest. Hence, the coefficient of initial effective stress (K₀), which is equal to the coefficient of earth pressure at rest condition, was calculated using Jaky’s formula K₀ = 1-sinf, where f is the angle of internal friction of the soil and the initial stresses of the top and bottom layers were taken to be 0.837 and 0.636, respectively.The analysis was performed for dry soil.

Soil models

For finite element analysis, the behavior of soils may be modeled using various models, such as the Mohr–Coulomb model, hardening-soil model, soft-soil-creep model, and soft soil model. In the present analysis, the soft-soil-creep model (SSC) was adopted because of the analytic consideration of the viscous effect. The SSC model is a slight extension of the Mohr-Coulomb model. When soils are modeled with the SSC model, the initial conditions for the Mohr-Coulomb model are also required. The input parameters used in the finite element analysis are listed in Table 2.

Calculation phases

The numerical analysis was performed in three phases.

1) In the first phase, a static load was applied, which would act along with the dynamic load to satisfy plastic equilibrium.

2) The second phase included dynamic analysis in which accelerogram earthquake data was applied at the base in the horizontal direction (Fig. 13) to predict the possible failure of the slope during the earthquake.

3) For the design of the slope, it is important to consider not only the stability but also the final stability during construction. Therefore, the global safety factor of the slope was computed after the first and second phases. The global safety factor of the model slope was calculated in the last phase using the shear strength reduction technique or phi-c reduction analysis to assess the safety factor of the model slope after plastic equilibrium was attained. In this phase, the cohesion and friction angle of the soil mass were reduced until a collapse mechanism and the corresponding FOS were determined.

Results comparison

Figure 14 shows that the peak ground acceleration found is 0.93 m/s². In the analytical method, the peak applied value is given by

PGA=(k_{\mathrm{h}}\times g)=(\mathrm{0.1}\times \mathrm{9.81})=\mathrm{0.981}\mathrm{m}/{\mathrm{s}}^{\mathrm{2}}.

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

A comparison of the failure surfaces obtained from the present analytical method and the numerical analysis are shown in Figs. 15(a) and 15(b). From the plot, it can be seen that the patterns of the failure surfaces obtained from both the analytical and numerical methods have the same trends, and the shape of the failure surface is nonlinear. Figure 16 shows a comparison of both the analytical and numerical failure surfaces obtained from the present analysis at H₁/H = 0 with the failure surface obtained from a 1D slope on a shake Table [⁴²]. From the comparison, it can be seen that the failure surface obtained from the analytical study lies between the experimental and numerical failure surfaces. A comparison of the results obtained from the proposed analytical method and numerical analysis is shown in Table 3. The FOS results obtained from the analytical solution are in good agreement with the FOS obtained from the numerical solution. A limitation of the numerical solution is that it does not show any value for the FOS if the the FOS is less than one but only indicates that the FOS is less than one if the soil fails.

To verify the present analytical approach, the values obtained from the present analysis were compared with those reported by Chatterjee and Krishna [²⁸]. Table 4 compares the values obtained from the present analysis with those obtained from the available literature. It is noticed that at the static condition, the present FOS values are in good agreement with those reported by Chatterjee and Krishna [²⁸]. However, under seismic loading conditions, the FOS values obtained from the present study are much lower than those obtained by Chatterjee and Krishna [²⁸]. This is because the slopes are more vulnerable to instability when subjected to earthquake loading, and hence, the FOS decreases.

The stability number values obtained from the present pseudo-dynamic analysis were also compared with the values obtained by Kumar and Samui [²⁹] and Sarkar and Chakraborty [³¹]. Table 5 presents a comparison of the values of stability number at k_v = 0 and c₁/c₂ = 1. From the comparison, it is noted that the stability number obtained from the present pseudo-dynamic method are much greater than the stability numbers obtained by Kumar and Samui [²⁹] and Sarkar and Chakraborty [³¹]. Kumar and Samui [²⁹] have used static and pseudo-static analysis, whereas Sarkar and Chakraboty [³¹] used the pseudo-static approach for the slope stability analysis.

However, the present analysis is performed using the pseudo-dynamic method which also considers the damping of the slope soil. It can be seen that stability number values obtained from the analytical solution are very close to the values obtained from the numerical analysis. Thecomparison therefore demonstrates the accuracy of the analytical solution.

Conclusions

In this study, two-layered cohesive-frictional soil slopes were analyzed using a modified pseudo-dynamic method in which the damping of the soil is also considered. The logarithmic spiral failure surface was considered for the analysis. The same problem was solved numerically using the PLAXIS 2D model, and the results were compared. From the analysis, it was seen that increasing the damping ratio, cohesion, angle of internal friction, unit weight, and height of the top soil layer increase the FOS of the slope, whereas the FOS of the slope decreases with an increase in the seismic acceleration and the slope inclination. The results obtained from the analytical solution are in good agreement with both the numerical and experimental results. Hence, the FOS obtained from the analytical solution can be easily used for practical problems.

The assumptions made in the problem solution result in a number of limitations. The erratic earthquake forces were assumed to be sinusoidal and the effects of waves other than the shear and primary waves ignored. The assumption of a logarithmic spiral surface passing through the toe of the slope may also result in deviations of the solution from actual field situations. Future studies will be needed to address these limitations.

In addition, in future works, the phase-field model (PFM) will be useful for numerically quantifying the gradually evolving failure surface and damage region in the two-layered slope, and for considering the effect of seepage on the slope stability [⁴³–⁴⁷]. A more detailed study of the artificial neural network method [⁴⁸], deep collocation method [⁴⁹], peridynamics modeling [⁵⁰], and meshfree method [⁵¹] are possible avenues for research on the two-layered slope problem.

References

[1] Fellenius W. Calculation of the stability of earth dams. In: Proceedings of the Second Congress of Large Dams. Washington, D.C., 1936, 4: 445–463

[2] Taylor D W. Stability of earth slopes. Journal of the Boston Society of Civil Engineers, 1937, 24: 197246

[3] Bishop A W. The use of slip circle in the stability analysis of earth slopes. Geotechnique, 1955, 5(1): 7–17

[4] Morgenstern N R, Price V E. The analyses of the stability of general slip surfaces. Geotechnique, 1965, 15(1): 79–93

[5] Spencer E. A method of analysis of the stability of embankments assuming parallel interslice forces. Geotechnique, 1967, 17(1): 11–26

[6] Terzaghi K. Mechanisms of Land Slides. Engineering Geology (Berkeley) Volume. New York: Geological Society of America, 1950

[7] Sarma S K. Stability analysis of embankments and slopes. Geotechnique, 1973, 23(3): 423–433

[8] Sarma S K. Stability analysis of embankments and slopes. Journal of Geotechnical Engineering, 1979, 105(12): 1511–1524

[9] Koppula S D. Pseudo-static analysis of clay slopes subjected to earthquakes. Geotechnique, 1984, 34(1): 71–79

[10] Leshchinsky D, San K C. Pseudostatic seismic stability of slopes: Design charts. Journal of Geotechnical Engineering, 1994, 120(9): 1514–1532

[11] Hazari S, Ghosh S, Sharma R P. Pseudo-static analysis of slope considering log spiral failure mechanism. In: Latha Gali M, Raghuveer Rao P, eds. Indian Geotechnical Conference. Guwahati: IIT Guwahati, 2017

[12] Steedman R S, Zeng X. The influence of phase on the calculation of pseudo-static earth pressure on a retaining wall. Geotechnique, 1990, 40(1): 103–112

[13] Choudhury D, Nimbalkar S. Seismic passive resistance by pseudo-dynamic method. Geotechnique, 2005, 55(9): 699–702

[14] Ghosh P. Seismic passive pressure behind a non-vertical retaining wall usingpseudo-dynamic method. Geotechnical and Geological Engineering, 2007, 25(6): 693–703

[15] Ghosh S, Sharma R P. Pseudo-dynamic active response of non-vertical retaining wall supporting c-F backfill. Geotechnical and Geological Engineering, 2010, 28(5): 633–641

[16] Saha A, Ghosh S. Pseudo-dynamic analysis for bearing capacity of foundation resting on c–f soil. Journal of Geotechnical Engineering, 2015, 9(4): 379–387

[17] Bellezza I. A new pseudo-dynamic approach for seismic active soil thrust. Geotechnical and Geological Engineering, 2014, 32(2): 561–576

[18] Bellezza I. Seismic active earth pressure on walls using a new pseudo-dynamic approach. Geotechnical and Geological Engineering, 2015, 33(4): 795–812

[19] [].Chanda N.Seismic stability analysis of slope assuming log-spiral rupture surface using modified pseudo-dynamic method. International Journal of Geotechnical Engineering, 2018in press)

[20] Fredlund D G, Krahn J. Comparison of slope stability methods of analysis. Canadian Geotechnical Journal, 1977, 14(3): 429–439

[21] Chen W F, Snitbhan N, Fang H Y. Stability of slopes in anisotropic, nonhomogeneous soils. Canadian Geotechnical Journal, 1975, 12(1): 146–152

[22] Koppula S D. Pseudo-static analysis of clay slopes subjected toearthquakes. Geotechnique, 1984, 34(71): 1–79

[23] Hammouri A, Malkawi A I H, Yamin M M A. Stability analysis of slopes using the finite element method and limit equilibrium approach. Bulletin of Engineering Geology and the Environment, 2008, 67(4): 471–478

[24] Chang-yu H C, Jin-jian X. Stability analysis of slopes using the finite element method and limiting equilibrium approach. Journal of Central South University of Technology, 2013, 21: 1142–1147

[25] Qian Z G, Li A J, Merifield R S, Lyamin A V. Slope stability charts for two-layered purely cohesive soils based on finite-element limit analysis methods. International Journal of Geomechanics, 2015, 15(3): 06014022

[26] Lin Y, Leng W, Yang G, Li L, Yang J. Seismic response of embankment slopes with different reinforcing measures in shaking table tests. Natural Hazards, 2015, 76(2): 791–810

[27] Ni P, Wang S, Zhang S, Mei L. Response of heterogeneous slopes to increased surcharge load. Computers and Geotechnics, 2016, 78: 99–109

[28] Chatterjee D, Krishna A M. Stability analysis of two-layered nonhomogeneous slopes. International Journal of Geotechnical Engineering, 2018,in press)

[29] Kumar J, Samui P. Stability determination for layered soil slopes using the upper bound limit analysis. Geotechnical and Geological Engineering, 2006, 24(6): 1803–1819

[30] Qin C, Chen Chian S. Kinematic stability of a two-stage slope in layered soils. International Journal of Geomechanics, 2017, 17(9): 06017006

[31] Sarkar S, Chakraborty M. Pseudostatic slope stability analysis in two-layered soil by using variational method. In: Proceedings of the 7th International Conference on Earthquake Geotechnical Engineering. Rome, Italy: Taylor & Francis, 2019, 4857–4864

[32] Kramer S L. Geotechnical Earthquake Engineering. New Jersey: Prentice Hall, 1996

[33] MATLAB: The MathWorks, Inc. Version R2013a.Natick, MA: The MathWorks, Inc., 2014

[34] Plaxis: Plaxis 2D.Reference Manual. Delft, Netherlands, 2012

[35] Pasternack S C, Gao S. Numerical methods in the stability analysis of slopes. Computers & Structures, 1988, 30(3): 573–579

[36] Matsui T, San K C. Finite element slope stability analysis by shear strength reduction technique. Soil and Foundation, 1992, 32(1): 59–70

[37] Duncan J M. State of the art: Limit equilibrium and finite element analysis of slopes. Journal of Geotechnical Engineering, 1996, 122(7): 577–596

[38] Griffiths D V, Lane P A. Slope stability analysis by finite elements. Geotechnique, 1999, 49(3): 387–403

[39] Chugh A K. On the boundary conditions in slope stability analysis. International Journal for Numerical and Analytical Methods in Geomechanics, 2003, 27(11): 905–926

[40] Kuhlemeyer R L, Lysmer J. Finite element method accuracy for wave propagation problems. Journal of the Soil Mechanics and Foundations Division, 1973, 99: 421–427

[41] Bakr J, Ahmad S M. A finite element performance-based approach to correlate movement of arigid retaining wall with seismic earth pressure. Soil Dynamics and Earthquake Engineering, 2018, 114: 460–479

[42] Hazari S, Ghosh S, Sharma R P. Experimental and numerical study of soil slopes at varying water content under dynamic loading condition. International Journal of Civil Engineering, 2019, 18: 215–229:

[43] Zhou S, Zhuang X, Rabczuk T. Phase-field modeling of fluid-driven dynamic cracking in porous media. Computer Methods in Applied Mechanics and Engineering, 2019, 350: 169–198

[44] Zhou S, Zhuang X, Rabczuk T. Phase field modeling of brittle compressive-shear fractures inrock-like materials: A new driving force and a hybrid formulation. Computer Methods in Applied Mechanics and Engineering, 2019, 355: 729–752

[45] Zhou S, Rabczuk T, Zhuang X. Phase field modeling of quasi-static and dynamic crack propagation: COMSOL implementation and case studies. Advances in Engineering Software, 2018, 122: 31–49

[46] Zhou S, Zhuang X, Rabczuk T. A phase-field modeling approach of fracture propagation in poroelastic media. Engineering Geology, 2018, 240: 189–203

[47] Zhou S, Zhuang X, Zhu H, Rabczuk T. Phase field modelling of crack propagation, branching and coalescence in rocks. Theoretical and Applied Fracture Mechanics, 2018, 96: 174–192

[48] Anitescu C, Atroshchenko E, Alajlan N, Rabczuk T. Artificial Neural Network methods for the solution of second order boundary value problems. Computers, Materials & Continua, 2019, 59(1): 345–359

[49] Guo H, Zhuang X, Rabczuk T. A deep collocation method for the bending analysis of Kirchhoff plate. Computers, Materials & Continua, 2019, 59(2): 433–456

[50] Ren H, Zhuang X, Rabczuk T. Dual-horizon peridynamics: A stable solution to varying horizons. International Journal for Numerical Methods in Engineering, 2004, 61(13): 762–782

[51] Rabczuk T, Belytschko T. Cracking particles: A simplified meshfree method for arbitrary evolving cracks. International Journal for Numerical Methods in Engineering, 2004, 61(13): 2316–2343

RIGHTS & PERMISSIONS

Higher Education Press