Evaluation of the stability of terraced slopes in clayey gravel soil using a novel numerical technique

Abstract

Conventional geotechnical software limits the use of the strength reduction method (SRM) based on the Mohr–Coulomb failure criterion to analyze the slope safety factor (SF). The use of this constitutive model is impractical for predicting the behavior of all soil types. In the present study, an innovative numerical technique based on SRM was developed to determine SF using the finite element method and considering the extended Cam–clay constitutive model for clayey gravel soil as opposed to the Mohr–Coulomb model. In this regard, a novel user subroutine code was employed in ABAQUS to reduce the stabilizing forces to determine the failure surfaces and resist and drive shear stresses on a slope. After validating the proposed technique, it was employed to investigate the performance of terraced slopes in the context of a case study. The impacts of geometric parameters and different water table elevations on the SF were examined. The results indicated that an increase in the upper and lower slope heights led to a decrease in SF, and a slight increase in the horizontal offset led to an increase in the SF. Moreover, when the water table elevation was lower than the toe of the terraced slope, the SF increased because of the increase in the uplift force as a resistant component.

Graphical abstract

Keywords

terraced slope / safety factor / finite element method / ABAQUS / extended Cam–clay

1 Introduction

Terraced slopes (also multi-level and multi-tier slopes) are groups of slopes constructed at different elevations [1], as shown in Fig.1. During slope construction, the cutting and filling operation costs are balanced against slope performance considerations. The widespread use of terraced slopes in roadway design has promoted research with respect to their stability.

Cheng et al. [2] compared the failure surface locations and safety factors (SFs) in different single-soil slope configurations obtained via traditional limit equilibrium methods (LEM) with the strength reduction method (SRM) considering the Mohr–Coulomb failure in FLAC software. Later, Bojorque et al. [3] examined Cheng et al.’s [2] results using PLAXIS geotechnical FEA software. For all slope configurations, the PLAXIS SFs were lower than those obtained by Cheng et al. [2], presumably because of the nonlinear solution algorithm. Additionally, the evaluation of a homogeneous and nonhomogeneous single slope by Xu et al. [4], using the reduction in cohesion and friction angle in the Mohr–Coulomb constitutive model in ABAQUS, demonstrated the ability and accuracy of this software in determining the sliding surface. Using the linear Mohr–Coulomb failure criterion, Zhao et al. [5] conducted numerical simulations of the stability of homogeneous slopes with benches based on upper bound limit analysis and SRM. They determined stability charts for estimating the stability status, failure modes, and formula for the slope SF in homogeneous slopes with benches. Zhao et al. [6] conducted numerical analyses using a finite element model (FEM) to determine the SF of homogeneous soil slopes by considering linear and nonlinear Mohr–Coulomb failure criteria. Their linear and nonlinear analysis results showed that the selection of the formulation technique in a nonlinear failure scenario significantly affected the calculated SF. Using the sliding block model, Zhang et al. [7,8] presented the permanent displacement of near-fault seismic slopes. As one of the main findings, they stated that the computed displacements were highly sensitive to the azimuth of the slope.

With a different focus, the effect of cracks and 3D shape of the slope were observed for single slopes using the Mohr–Coulomb failure criterion [9–12]. Li et al. [13] examined the simultaneous effect of existing cracks in a soil mass, along with seismic loads, on the SF of a three-dimensional slope. They demonstrated that the presence of cracks and seismic loads diminished slope stability. Additionally, they reported that the consideration of the 3D effects results in a higher SF. Yang and Zhang [14] later conducted analytical research on the stability of a three-dimensional slope reinforced with piles. They considered cracks in their model and used the kinematic approach of the limit analysis to calculate the SF. They reported that slopes with large inclination angles were sensitive to natural cracks, whereas slopes with mild inclination angles were less influenced. Zhang et al. [15] investigated the effect of the tension zone on the stability of a single slope under static and dynamic loads, considering the Mohr–Coulomb method using two approaches: a) introducing cracks in the slope; b) reducing the tensile strength of the soil. They showed that the consideration of tension zone effects can sharply reduce the stability factor of slopes, particularly for steep slopes and large horizontal seismic loads.

In the previous studies, researchers investigated the impacts of mechanical properties of slopes, along with other parameters, such as earthquakes, rainfall, or variations in groundwater elevation, on slope SF. They proposed methods to determine the SFs under various conditions based on the Mohr–Coulomb failure criterion. However, based on the type of soil, i.e., gravel, sandy, clayey, or a combination of the three, the utilization of different constitutive models has a substantial effect on the soil response with respect to mechanical, hydraulic, or thermal loads. For instance, in accordance with computational geotechnics theories [16-18], the Mohr−Coulomb constitutive model is suitable for granular materials such as sandy soils. However, Drucker–Prager and Cam–clay can be used for soils when the impact of volume variation on the bulk property and resistance to shear should be considered, as in the case of soils including clay [19]. Therefore, the selection of a proper constitutive model for soil is a significant factor for analyzing slope stability. Nevertheless, the common geotechnical software and analytical methods used by researchers rely on conventional Mohr–Coulomb failure criterion according to the SRM to calculate the SF in slopes. In the present study, we first introduce a novel numerical technique using the finite element method considering the extended Cam–clay constitutive model for soil, as opposed to the routine Mohr−Coulomb criterion, to analyze the stability and safety of slopes. After validating the proposed technique, it was employed to evaluate the stability and safety of terraced slopes in the context of a case study. Finally, in this study, we investigated the impact of slope geometry and water-table elevation on the SF of terraced slopes using the proposed method.

2 Methodology and objectives

In the present study, by using FEM and the geostatic stress reduction technique in soil elements, the SF of a two-dimensional terraced slope, consisting of two slopes, was evaluated with different dimensional parameters and water-table elevations. Based on the plain-strain condition, the ABAQUS software was used to model and analyze the displacements and plastic strains. First, by considering the type of soil at the slope location (clayey gravel), three constitutive models, namely extended Cam–clay, modified Drucker–Prager/Cap, and Mohr–Coulomb, were evaluated based on laboratory data. Then, a suitable constitutive model was selected to simulate the shear failure in soil owing to geostatic stress reduction. UFIELD (user subroutine to redefine field variables at nodal points) was implemented in ABAQUS to simulate stress reduction, and finally, calculate the slope’s SF by considering nonlinear deformations.

After verifying the proposed numerical technique with respect to LEM, the method was used to evaluate the stability of a terraced slope. The effects of geometric factors (dimensional variables) and groundwater elevation on the SF were investigated by defining different scenarios (33 scenarios). The dimensional variables were slope heights and angles as well as the horizontal offset between the two slopes. Additionally, by defining the boundary of the free water surface, the water table level was considered to vary between the bottom of the model and crown of the upper slope. After running all scenarios, the obtained SFs are presented in the diagrams. Subsequently, data points are interpreted and discussed, and conclusions and recommendations are presented based on the results.

2.1 Case study

In this study, to investigate the stability of terraced slopes, which are constructed by cutting and filling and adjacent to a valley with seasonal changes in water table elevations, a slope located east of Las Vegas, near the Colorado River, was selected as the case study. Based on the unified soil classification system and the sampling of six points (B-1 to B-6), GC (clayey gravel) was determined as the dominant soil group in this location (Tab.1). Thus, the GC soil type was considered for all models in this study.

2.2 Constitutive model

Three elasto-plastic constitutive models, including extended Cam–clay, modified Drucker–Prager/Cap, and Mohr–Coulomb, were evaluated in ABAQUS based on laboratory data to evaluate and select a suitable constitutive model of the soil in the study area.

Compared with uniaxial conditions, in which only a yield stress exists, in multiaxial conditions, the yield stress is defined as a function. The yield function, F, is defined as a scalar function of stress and state parameters {k} [20].

(1)F(\{{\sigma }^{\prime }\}, \{\mathit{k}\})=0.

This function separates the elastic from the elasto–plastic behavior. In general, the surface is a function of the stress state ({\sigma }^{\prime }), and its size changes as a function of the state parameters {k}, which can be related to the hardening−softening parameters. For perfect plasticity, {k} is constant and denotes the magnitude of the stresses at the yield. For hardening and softening plasticity, {k} varies with plastic strain to denote the magnitude of the stress state with respect to yield changes. The value of the yield function F is used to identify the type of material behavior. It should be noted that yield functions can be investigated in various spaces, which can be transformed using map functions.

2.2.1 Extended Cam–clay model

In the extended Cam–clay (Eq. (2)), an elasticity model, in which the bulk elastic stiffness increases as the material undergoes compression, was used to calculate the elastic strains [21–24]. The plastic strains were calculated using the theory of plasticity, i.e., a yield surface with the associated flow (the plastic potential is the same as the yield surface) along with a hardening rule that allows the yield surface to grow or shrink in the stress space. In Eq. (2), p' denotes the mean effective stress (Eq. (2a)); \chi is a constant, which is used to modify the shape of the yield surface on the wet-side of the critical state; and a denotes a hardening parameter, defined as the point on p'-axis where the evolving elliptic arcs of the yield surface intersect the critical state line (CSL). For \chi = 1, a = p_{\mathrm{c}}^{\prime }/ 2(p_{\mathrm{c}}^{\prime } denotes pre-consolidation stress).

The projection of the critical state surface on p'–t plane corresponds to a straight line passing through the origin with slope M. Specifically, slope M of the CSL in p'–t plane can be calculated from the internal friction angle \phi obtained from the triaxial test results (Eq. (2b)). Parameter t is a measure of shear stress, as defined in Eq. (2c), and K denotes the ratio of the flow stress in triaxial tension to the flow stress in triaxial compression (0.778 ≤ K≤ 1.00). By setting K = 1, the yield surface becomes independent of the third stress invariant (r) [22,25]. In Eq. (2d), q denotes the deviator stress.

(2)F(\{{\sigma }^{\prime }\}, \{\mathit{k}\})=\frac{1}{{\chi }^2}{\left(\frac{p^{\prime }}{a}-1\right)}^2+{\left(\frac{t}{Ma}\right)}^2=0,

\begin{array}{}\text{(2a)}& p^{\prime }=\frac{{\sigma }_1^{\prime }+{\sigma }_2^{\prime }+{\sigma }_3^{\prime }}{3},\end{array}

\begin{array}{}\text{(2b)}& M=\frac{6\sin {\phi }^{\prime }}{3-\sin {\phi }^{\prime }},\end{array}

\begin{array}{}\text{(2c)}& t=\frac{q}{2}\left[1+\frac{1}{K}-\left(1-\frac{1}{K}\right){\left(\frac{r}{q}\right)}^3\right],\end{array}

\begin{array}{}\text{(2d)}& q=\sqrt{\frac{1}{2}[{({\sigma }_1^{\prime }-{\sigma }_2^{\prime })}^2+{({\sigma }_2^{\prime }-{\sigma }_3^{\prime })}^2+{({\sigma }_1^{\prime }-{\sigma }_3^{\prime })}^2]}.\end{array}

2.2.2 Drucker–Prager/Cap model

The modified Drucker–Prager/Cap plasticity model has been widely used in finite element analysis programs for various geotechnical engineering applications [22,26]. The cap model is appropriate for simulating soil behavior because it is capable of considering the effects of stress history, stress path, and dilatancy, as well as the effect of the intermediate principal stress. The yield surface of the modified Drucker–Prager/Cap plasticity model consists of three parts: a Drucker–Prager shear failure surface; an elliptical cap that intersects the mean effective stress axis at a right angle; and a smooth transition region between the shear failure surface and cap. An elasticity model, in which the bulk elastic stiffness increases as the material undergoes compression, can be used to calculate the elastic strain. The onset of plastic behavior is determined by the Drucker–Prager failure and cap yield surfaces. The Drucker–Prager failure surface can be expressed as follows.

(3)F_{\text{s}} (\{{\sigma }^{\prime }\},\{\mathit{k}\})= t-p\tan \beta -d=0,

where β denotes the soil’s angle of friction and d denotes its cohesion in p'–t plane. The cap yield surface is an ellipse in p'–t plane. The cap yield surface is dependent on the third stress invariant (r). The cap surface hardens (expands) or softens (shrinks) as a function of volumetric plastic strain. When the stress state causes yielding on the cap, volumetric plastic strain (compaction) occurs. This in turn leads to expansion (hardening) of the cap. However, when the stress state causes yielding on the shear failure surface, volumetric plastic dilation occurs, causing the cap to shrink (softening) [22,25]. The cap yield surface is can be expressed as follows.

(4)\begin{array}{rl}& F_{\text{c}}(\{{\sigma }^{\prime }\},\{\mathit{k}\})\\ =& \sqrt{{(p-p_{\text{a}})}^2+\frac{Rt}{1+a-a/\cos \beta }-R(d+p_{\text{a}}\tan \beta )}=0,\end{array}

where R denotes the cap eccentricity parameter, which must be greater than zero (typically 0.0001 ≤ R≤ 1000.0) [22,25], and α denotes the transition surface radius parameter, which is a small number (typically, 0.01 to 0.05) that can be used to define a smooth transition surface between the Drucker–Prager shear failure surface and cap (Eq. (5)) as follows.

(5)\begin{array}{rl}{F}_{\text{t}}(\{{\sigma }^{\prime }\},\{\mathit{k}\})=& \sqrt{{(p-p_{\text{a}})}^2+{\left[t-\left(1-\frac{a}{\cos \beta }\right)(d+p_{\text{a}}\tan \beta )\right]}^2}\\ & -a(d+p_{\text{a}}\tan \beta )=0,\end{array}

where p_a denotes an evolution parameter that controls the hardening–softening behavior as a function of the volumetric plastic strain (Eq. (6)).

(6)p_{\text{a}} =\frac{p_{\text{b}}-Rd}{1+R\tan \beta }.

The hardening–softening behavior can be described by a piecewise linear function that relates the mean effective (yield) stress p_b and volumetric plastic strain p_b = p_b({\epsilon }_{\mathrm{v}\mathrm{o}\mathrm{l}}^{\mathrm{p}\mathrm{l}}). This function can be obtained from the results of the isotropic consolidation test. Using p_b, p_a is calculated.

2.2.3 Mohr–Coulomb model

The Mohr–Coulomb model is based on Mohr’s circle for the states of stress at failure in the plane of maximum and minimum principal stresses [22]. The Mohr–Coulomb failure criterion is adopted as the yield function (Eq. (7)).

(7)F(\{{\sigma }^{\prime }\}, \{\mathit{k}\})=R_{\text{mc}}q-p\tan \phi -c=0.

where R_mc denotes the Mohr–Coulomb deviatoric stress measure, which can be defined as follows.

(8)R_{\text{mc}}( \theta ,\phi )=\frac{1}{\sqrt{3}\cos \phi }\sin \left(\theta +\frac{\text{π}}{3}\right)+\frac{1}{3}\cos \left(\theta +\frac{\text{π}}{3}\right)\tan \phi ,

where θ denotes the deviatoric polar angle [27], defined by Eq. (9).

(9)\cos (3 \theta )= {\left(\frac{r}{q}\right)}^3.

2.3 Model calibration

A representative soil sample of the slope in class GC was selected and tested using a one-dimensional consolidation test [28] (Fig.2) and consolidated drained (CD) triaxial compression test [29] (Fig.3). The parameters required in the constitutive models were determined to simulate the triaxial test using ABAQUS software (Tab.2). Subsequently, by comparing the numerical and laboratory results of the triaxial test, a constitutive model, which can determine the shear strength parameters (c^{\prime } and {\phi }^{\prime }) of the representative soil of the area with a higher degree of accuracy, was selected.

Based on the laboratory tests, k (permeability), e (void ratio), γ_t (wet specific weight), γ_d (dry specific weight), G_s (specific gravity), and E (elastic modulus) were obtained. Specifically, λ (the slope of loading for e–lnp′ curve) and κ (the slope of unloading–reloading) were calculated from the oedometer test (Fig.2(a)). Considering the extended Cam–clay plasticity, the stress ratio M is calculated using Eq. (2b), where a denotes a hardening parameter, which was determined via a trial-and-error procedure. Furthermore, \chi (the shape controller of the yield surface) in the extended Cam–clay, R (the curvature of the cap part of the yield surface) and α (the transition surface parameter) in the modified Drucker–Prager/Cap plasticity, and K (the ratio of the flow stress) in extended Cam–clay and modified Drucker–Prager/Cap plasticity, were calibrated using ABAQUS. For the modified Drucker–Prager/Cap plasticity, cohesion, d, and friction angle, β, were calculated as d={\displaystyle \frac{18c\cos {\phi }^{\prime }}{3-\sin {\phi }^{\prime }}} and \tan \beta ={\displaystyle \frac{6\sin {\phi }^{\prime }}{3-\sin {\phi }^{\prime }}} [20]. The cap hardening curve is obtained from the isotropic consolidation test results as shown in Fig.2(b). The initial cap position {\epsilon }_{\mathrm{v}\mathrm{o}\mathrm{l}}^{\mathrm{p}\mathrm{l}} is equal to zero, which corresponds to the confining pressure (4.5 kg/cm²). The transition surface parameter α should be a small number when compared to unity. Additionally, Mohr–Coulomb shear strength parameters were obtained from the triaxial tests.

A cylindrical specimen of soil with a diameter of 3.8 cm and height of 7.6 cm was modeled using the FE method. Furthermore, CPE8RP (20-nodes, with pore pressure, quadratic displacement, linear pore pressure, and reduced integration) was selected to mesh the soil environment. Based on the results obtained from the triaxial device (Fig.3), the vertical displacement rate of the upper cap was 0.026 mm/min. Fig.4 shows the boundary conditions used in the numerical model. Fig.5 presents the deviator stress obtained from the laboratory test and numerical model for different constitutive models. As shown, the extended Cam–clay determines the shear strength parameters of the area’s representative soil with higher accuracy than the modified Drucker–Prager/Cap or Mohr–Coulomb methods.

3 Numerical modeling

3.1 Use of extended Cam–clay constitutive model in the strength reduction method

According to the elastoplastic behavior of the soil in the extended Cam–clay in ABAQUS, the relationship between the total stress (σ') and strains (ε) can be expressed as follows [22,30].

(10)\frac{\{\mathrm{\Delta }{\sigma }^{\prime }\}}{[K^{\text{ep}}]}=(\{\mathrm{\Delta }\epsilon \}),

where Δε^e denotes the total strain, including Δε^e (elastic strain) and Δε^p (plastic strain), and K^ep denotes an elastoplastic consecutive stiffness matrix. Furthermore, [K^ep] can be expressed as follows.

(11)[K^{\text{ep}}]=[K^{\mathrm{e}} ]-\frac{[K^{\mathrm{e}}]\left\{{\displaystyle \frac{\mathrm{\partial }P(\{{\sigma }^{\prime }\},\{\mathit{m}\})}{\mathrm{\partial }{\sigma }^{\prime }}}\right\}{\left\{{\displaystyle \frac{\mathrm{\partial }F(\{{\sigma }^{\prime }\},\{\mathit{k}\})}{\mathrm{\partial }{\sigma }^{\prime }}}\right\}}^{\text{T}}[K^{\text{e}}]}{{\left\{{\displaystyle \frac{\mathrm{\partial }F(\{{\sigma }^{\prime }\},\{\mathit{k}\})}{\mathrm{\partial }{\sigma }^{\prime }}}\right\}}^{\text{T}}[K^{\text{e}}]\left\{{\displaystyle \frac{\mathrm{\partial }P(\{{\sigma }^{\prime }\},\{\mathit{k}\})}{\mathrm{\partial }{\sigma }^{\prime }}}\right\}+Y},

where [K^e], P({{\sigma }^{\prime }},{m}), and Y denote the elastic stiffness matrix, plastic potential function, and plastic parameters, respectively. The type of plastic behavioral model, with perfect plasticity or strain hardening−softening plasticity, is dependent on parameter Y. In this study, the soil in the extended Cam–clay is considered as perfectly plastic (Y = 0, and K is constant) [20]. The behavior of the soil is nonlinear (elasto–plastic); therefore, the constitutive stiffness matrix is not constant during analysis, but varies with stress. By defining a loop and the UFIELD subroutine (allowing the user to prescribe predefined field variables at the nodes of a model), first, the initial stresses are assigned to the model by defining the initial conditions (geostatic stresses), and then, they are used in Eq. (11). Using Re(n), a scalar reduction function, as a coefficient {Δ{\sigma }^{\prime }} in Eq. (10), the stresses on the nodes can be reduced (Eq. (12)) as follows.

(12)\frac{Re(n).\{\mathrm{\Delta }{\sigma }^{\prime }\}}{[K^{\text{ep}}]}=(\{\mathrm{\Delta }\epsilon \}).

By assembling the global stress, stiffness, and strain matrices for all elements, Eq. (12) is solved simultaneously for all elements of the model. To calculate the SF, when the first significant strain occurs (the first divergence occurs when solving Eq. (12)), the initial shear stress in the non-convergent elements (stabilizing stress) is divided by the maximum shear stress occurring in those elements (destabilizing stress). Fig.6 shows the procedure, as a flowchart, for obtaining the sliding surfaces and SF using the proposed technique.

3.2 Modeling procedure

First, verification was performed to ensure the performance of the proposed method and accuracy of the results. The slope reported in Ref. [2] (Fig.7) is used in this study. Specifically, the shear surface and SFs are determined using the technique described in this study.

To realize acceptable results, three reduction rates are considered: τ < 1; τ = 1; and τ > 1 (Fig.8(a)). Given that the problem must be solved under quasi-static conditions and ratio of kinetic energy to internal energy should be less than 10% (E_k/E_i < 10%) [22], this ratio is checked at the end of the analysis (Fig.8(b)). In this study, after trial and error, the nonlinear reduction function with downward concavity (τ = 1/3) is selected as the best function for stress matrix reduction with minimal dynamic effects. Fig.9 indicates that the results of the numerical technique proposed in this study are in good agreement with the LEM results [2] in terms of forming slide surfaces and SF values. Additionally, the proposed technique leads to conservative values for the SF when compared with that of LEM.

As shown in Fig.9, five shear zones emerge on the terraced slope. The SFs calculated using the FEM are nearly identical to the LEM results. After the verification procedure, the effects of geometric variables and water table elevations on the SF are evaluated using the main model (Fig.10).

In this study, the SF depends on six geometric parameters: SF = f (H₂, H₁, L, W,β₂,β₁), where H₂, H₁, L, W, β₂, and β₁ are the upper slope height, lower slope height, horizontal offset, water table elevation, upper slope angle, and lower slope angle, respectively. The sensitivity of the aforementioned parameters and their effects on the SF are examined. Specifically, W varies in different scenarios, which consequently affects the effective stress in the soil. Moreover, owing to the hydraulic gradient, the seepage force affects the SF of the slope. The geometric boundaries at the sides are fixed in the direction of length (x), whereas at the bottom, they are fixed in the directions of length and width (x,y). By defining infinite elements at the boundaries, the geometric boundaries do not impact the results. The soil is considered as a monolithic elastoplastic material, and an extended Cam–clay behavioral model is implemented to simulate the problem. The soil properties are similar to the characteristics of the area’s representative soil (Tab.2, extended Cam–clay).

After mesh and element sensitivity analyses, CPE4 (4-node bilinear plane strain quadrilateral) element is selected for meshing the model [31–36]. Different scenarios are considered to investigate the impact of different parameters on the SF and slope stability as listed in Tab.3.

Based on these scenarios, the stability of the terraced slope is analyzed, and a sensitivity analysis is performed for each variable. To interpret the results, two terms must be defined: local and global. Global failure occurs when a slide surface forms on the entire slope and the slope collapses. Conversely, local failure occurs when more than one slide surface forms at different locations on the slope. The downward and upward displacements (δ) are considered as positive and negative displacements, respectively. Four benchmark points A, B, C, and D are considered to examine slope displacements (Fig.10). Nodes A and B denote the crown and toe of the upper slope, respectively, whereas nodes C and D denote the crown and toe of the lower slope, respectively. Thus, the stability of the entire slope (AD) and each slope (AB and CD) are checked separately and simultaneously. For the parametric study, a dimensionless number is used to investigate the vertical displacements. For all scenarios, the minimum SF calculated for the slopes is considered as the main SF for the stability of the terraced slope. In the next section, the impacts of the geometric parameters and water table elevation on the SF are investigated. Fig.11 shows the meshing pattern, formation of the sliding surfaces, and displacement vectors in the model for scenario S-18.

4 Results and discussion

4.1 Determining the slope safety factor in different scenarios

According to the results in Fig.12, in case (a), H₁ = 30 m, the upper and lower slopes collapse simultaneously, and the SF is equal to 1.579 for the entire slope (global failure occurred). Nodes A, B, and C recorded continuous downward displacements. Node D initially exhibited downward movement, which was replaced by upward movement upon the formation of a global slide surface and the collapse of the entire slope at SF = 1.579, and thereby, indicating toe uplift. Cases (b), (c), and (d) followed the same pattern as in case (a). In case (e), H₁ = 10 m, node B recorded upward displacements, and the upper slope collapsed prior to global failure. Hence, nodes A and B, located on the upper slope, exhibited the same SF (SF = 1.808). After the collapse of the upper slope, the entire slope failed at SF = 1.891 (nodes C and D). In case (f), H₁ = 5 m, and the upper slope failed prior to the global collapse at SF = 1.9.

As illustrated in Fig.13, a reduction in the upper slope height generates dramatic displacements at all benchmark nodes, with global failure occurring for H₂ = 30, 25, and 15 m, and global slip surfaces emerge at SF = 1.595, 1.635, and 1.77, respectively. In cases (d) and (e), the lower slope failed earlier than the entire slope at SF = 1.819 and SF = 1.823. At lower upper slope heights, the SFs increased because of the overburden load reduction imposed by the upper slope. Furthermore, for H₂ < 10 m, the impact of H₂ reduction was extremely low.

According to Fig.14, as the horizontal offset increases, the SF decreases, and in all cases, global slip surfaces emerged because the critical SFs occurred simultaneously. The SFs were SF = 1.756, 1.668, and 1.602 for L = 25, 15, and 5 m, respectively. Moreover, the influence of the different horizontal offset values on the SFs was minimal. For instance, when L was increased from 5 to 25 m (80% increase), the SF changed by approximately 8%.

With a decrease in the lower slope angle (Fig.15), in case (a), two slide surfaces were formed: first, the lower slope collapsed, and then, the entire slope failed (local failure). This failure pattern was repeated until case (c). After that, with diminishing β₁, the upper slope initially failed at SF = 1.832, and then, the entire slope collapsed at SF = 1.847. As observed in case (d), node B exhibits two types of displacements. When the upper slope failed, because this node was located at the toe, it exhibited upward displacement. Subsequently, when the entire slope collapsed, the displacement of the node suddenly decreased. As shown in Fig.15, the same occurred for case (e) regarding displacements, but the rate of upward displacement was higher than that in case (d).

Considering the changes in β₂ (Fig.16) in cases (a) and (b), the upper slope failed when the entire slope collapsed. Conversely, with a greater decrease in β₂, first, the lower slope failed, and then, the entire slope collapsed. The results indicate that the impact of different angles on the SF of the terraced slope was considerable. For example, by lowering β₂ from 45° to 15° (60% reduction), the SF increased by approximately 22%.

Soil slopes located near rivers exhibit different behaviors in different seasons because of the changes in the elevation of the water table. In dry seasons within the southwestern portion of the United States, rivers do not ordinarily exhibit water flow; hence, the water table elevation is lower than the riverbed (case (e)). During rainy seasons, the water table elevation increases, and water flows from the slope to the river because of the hydraulic gradient (cases (a), (b), and (c)). The flow conditions and rates in a river vary throughout the year. Case (d) represents the SF values for low water surface elevation conditions in a river. In each case, the behavior of the slope differs because of variations in the pore pressure and specific weight of the soil. Therefore, based on the water table elevation, this parameter positively influenced the SF at w = 10, 20, and 30 m; however, it negatively impacted the SF for w > 40 m (Fig.17). In all scenarios, the upper and lower slopes failed at the same time, and one slide surface was formed. Additionally, the SF decreased with a low rate of increase in the water table elevation.

4.2 Sensitivity analysis of the parameters affecting the slope safety factor

To investigate the impact of geometric variables on the SF, the normalized parameters, Ĥ₁, Ĥ₂, {\hat{\beta }}_1, {\hat{\beta }}_2, and \hat{L} were defined (Fig.18). These parameters were obtained by dividing the secondary value by the initial value of each variable. By decreasing H₁, H₂,β₁, and β₂, the main SF of the terraced slope was increased. However, by lowering L, SF was decreased slightly. Owing to the reduction in overburden and passive pressure by decreasing H₂ and H₁, for H_1,2≤ 0.5H_initial1,2, decreasing H₁ and H₂ did not significantly affect the SF. Regarding β₁ and β₂, when β_1,2≤ 0.5β_initial1,2, the SF was increased by reducing the angles. Furthermore, different slope stability behaviors due to the water table elevations are illustrated in Fig.18(f). As shown, depending on the water-table elevation and changes in the soil unit weight and seepage force, the SF exhibited a parabolic pattern. When W was lower than the toe of the terraced slope, the increased pore pressure (uplift forces) at the depth aided the slope’s stability and prevented slide surface formation. When W was higher than the crown of the lower slope, owing to the impact of pore pressure on the slope and soil saturation, an increase in W led to a destructive effect on the SF. Tab.4 lists the results for each scenario.

The results show that in the design of terraced slopes, because of the presence of surcharge loads due to the upper slope and interaction of slopes, special attention must be given to the geometrical conditions and groundwater level. For instance, based on the different scenarios in this research (33 scenarios), depending on the height of the upper and lower slopes, sliding surfaces can be formed initially in the upper slope, lower slope, or simultaneously in upper and lower slopes. Moreover, the horizontal offset between the upper and lower slopes did not affect the stability of each slope. This geometric factor can only affect the stability of the entire slope and formation of the global sliding surface. Another aspect of this research involved investigating the impact of the groundwater level on the stability of terraced slopes. In this regard, the use of a constitutive model, which can represent the correct behavior of the soil, is necessary because the soil in the study area contained significant amounts of clay that is sensitive to water. The selection of the extended Cam–clay constitutive model for the soil in the study area showed that the groundwater level can be a stabilizing or destabilizing factor for terraced slopes.

5 Conclusions

In this study, a numerical analysis was performed on the stability of a terraced slope (consisting of two slopes) to investigate the effects of dimensional parameters and water-table elevation on the slope SF. Based on the extended Cam–clay constitutive model for soil, the SF was calculated using a novel solution based on FEM in ABAQUS. Ultimately, the effects of the upper and lower slope heights, horizontal offsets, angles, and water table elevations on the SF were investigated separately.

The main findings of this study are as follows.

Using the proposed numerical technique, as opposed to the routine and simple Mohr−Coulomb failure criterion, the SF of slopes can be calculated in different complex constitutive models, such as extended Cam–clay. The application of the proposed numerical technique, along with the nonlinear reduction function with downward concavity for force matrix reduction, had minimal dynamic effects on the extended Cam–clay constitutive model in GC soil. An increase in the upper and lower slope heights (50% more than the initial value) decreased the SF. Furthermore, upon increasing the horizontal offset, the SF slightly increased. Additionally, by increasing the upper and lower slope angles (50% more than the initial value), a substantial decrease in the SF was observed. The effect of the slope angle on the stability of the terraced slope was higher than that of other factors such as the slope height and horizontal offset. Depending on the water table elevation, when the level was under the river bed, it had a positive effect on the SF due to uplift forces, but when the level was above the horizontal offset line, it had a negative effect.

References

[1] Camera C, Apuani T, Masetti M. Modeling the stability of terraced slopes: An approach from Valtellina (Northern Italy). Environmental Earth Sciences, 2015, 74(1): 855–868

[2] Cheng Y M, Lansivaara T, Wei W B. Two-dimensional slope stability analysis by limit equilibrium and strength reduction methods. Computers and Geotechnics, 2007, 34(3): 137–150

[3] Bojorque J, de Roeck G, Maertens J. Comments on ‘Two-dimensional slope stability analysis by limit equilibrium and strength reduction methods’ by Y.M. Cheng, T. Länsivaara and W.B. Wei [Computers and Geotechnics 34 (2007) 137–150]. Computers and Geotechnics, 2008, 35(2): 305–308

[4] Xu Q, Yin H, Cao X, Li Z. A temperature-driven strength reduction method for slope stability analysis. Mechanics Research Communications, 2009, 36(2): 224–231

[5] Zhao L, Xia P, Xie R, Li L, Zhang Y, Cheng X. Stability analysis of homogeneous slopes with benches. Geomechanics and Engineering, 2017, 13(3): 517–533

[6] Zhao L, Yang F, Zhang Y, Dan H, Liu W. Effects of shear strength reduction strategies on safety factor of homogeneous slope based on a general nonlinear failure criterion. Computers and Geotechnics, 2015, 63: 215–228

[7] Zhang Y, Xiang C, Yu P, Zhao L, Zhao J X, Fu H. Investigation of permanent displacements of near-fault seismic slopes by a general sliding block model. Landslides, 2022, 19(1): 187–197

[8] Zhang Y, Xiang C, Fu H, Liu J, Cheng Y. Effect of excitation-applied manners on permanent displacements of planar slopes using dynamic sliding blocks analysis. International Journal of Geomechanics, 2022, 22(5): 04022042

[9] Zhang Y, Chen G, Zheng L, Li Y, Zhuang X. Effects of geometries on three-dimensional slope stability. Canadian Geotechnical Journal, 2013, 50(3): 233–249

[10] Zhang Y, Chen G, Wang B, Li L. An analytical method to evaluate the effect of a turning corner on 3D slope stability. Computers and Geotechnics, 2013, 53: 40–45

[11] Michalowski R L. Stability of intact slopes with tensile strength cut-off. Geotechnique, 2017, 67(8): 720–727

[12] Zhao L H, Cheng X, Zhang Y, Li L, Li D J. Stability analysis of seismic slopes with cracks. Computers and Geotechnics, 2016, 77: 77–90

[13] Li Z W, Yang X L, Li T Z. Static and seismic stability assessment of 3D slopes with cracks. Engineering Geology, 2020, 265: 105450

[14] Yang X L, Zhang S. Stability analysis of 3D cracked slope reinforced with piles. Computers and Geotechnics, 2020, 122: 103544

[15] Zhang Y B, Liu Y, Yuan R, He Y. Comparison of seismic stability for slopes with tensile strength cut-off and cracks. Journal of Mountain Science, 2021, 18(12): 3336–3347

[16] BrinkgreveR B. Selection of soil models and parameters for geotechnical engineering application. In: Proceedings of Geo-Frontiers Congress. Austin: ASCE, 2005, 69–98

[17] LadeP V. Overview of constitutive models for soils. In: Proceedings of Geo-Frontiers Congress. Austin: ASCE, 2005, 1–34

[18] Taborda D M, Pedro A M, Pirrone A I. A state parameter-dependent constitutive model for sands based on the Mohr−Coulomb failure criterion. Computers and Geotechnics, 2022, 148: 104811

[19] Basack S, Karami M, Karakouzian M. Pile-soil interaction under cyclic lateral load in loose sand: Experimental and numerical evaluations. Soil Dynamics and Earthquake Engineering, 2022, 162: 107439

[20] PottsD MZdravković LAddenbrookeT IHigginsK GKovačević N. Finite Element Analysis in Geotechnical Engineering: Application. London: Thomas Telford, 2001

[21] Borja R I. Cam−clay plasticity, Part II: Implicit integration of constitutive equation based on a nonlinear elastic stress predictor. Computer Methods in Applied Mechanics and Engineering, 1991, 88(2): 225–240

[22] DassaultSystemes. Theory Manual for ABAQUS Version 6.12-3, 2012

[23] SchofieldA NWrothP. Critical State Soil Mechanics. London: McGraw-Hill, 1968

[24] ParryR H G. Stress−Strain Behaviour of Soils. England: G. T. Foulis & Co., 1972

[25] HelwanyS. Applied Soil Mechanics with ABAQUS Applications. New Jersey: John Wiley & Sons, 2007

[26] Resende L, Martin J B. Formulation of Drucker−Prager cap model. Journal of Engineering Mechanics, 1985, 111(7): 855–881

[27] ChenW FHan D J. Plasticity for Structural Engineers. New York: Springer, 1988

[28] D2435/D2435M-11. Standard Test Methods for One-Dimensional Consolidation Properties of Soils Using Incremental Loading. West Conshohocken: ASTM, 2011

[29] D7181-11. Method for Consolidated Drained Triaxial Compression Test for Soils. West Conshohocken: ASTM, 2011

[30] Potts D M, Zdravkovic L. Accounting for partial material factors in numerical analysis. Geotechnique, 2012, 62(12): 1053–1065

[31] Karami M, Kabiri-Samani A, Nazari-Sharabian M, Karakouzian M. Investigating the effects of transient flow in concrete-lined pressure tunnels, and developing a new analytical formula for pressure wave velocity. Tunnelling and Underground Space Technology, 2019, 91: 102992

[32] Karakouzian M, Karami M, Nazari-Sharabian M, Ahmad S. Flow-induced stresses and displacements in jointed concrete pipes installed by pipe jacking method. Fluids, 2019, 4(1): 34

[33] Karakouzian M, Nazari-Sharabian M, Karami M. Effect of overburden height on hydraulic fracturing of concrete-lined pressure tunnels excavated in intact rock: A numerical study. Fluids, 2019, 4(2): 112

[34] Liu Y, He Z, Li B, Yang Q. Slope stability analysis based on a multigrid method using a nonlinear 3D finite element model. Frontiers of Structural and Civil Engineering, 2013, 7(1): 24–31

[35] Lian J, Wu J. Surficial stability analysis of soil slope under seepage based on a novel failure mode. Frontiers of Structural and Civil Engineering, 2021, 15(3): 712–726

[36] Belandria N, Úcar R, León F M, Hassani F. Stability analysis of slopes with planar failure using variational calculus and numerical methods. Frontiers of Structural and Civil Engineering, 2020, 14(5): 1262–1273

RIGHTS & PERMISSIONS

Higher Education Press