Stability analysis for two-layered slopes by using the strength reduction method

The aim of this article is to present the slope stability charts for two layered soil slopes by using the strength reduction method (SRM). The primary focus is to provide a quantitative estimation of the improvement of slope stability when a stronger layer is placed over the weaker layer. The SRM carried in this work comprises a series of finite element lower bound (LB) and upper bound (UB) limit analysis in conjunction with nonlinear optimization. Unlike the limit equilibrium method (LEM), there is no need to consider any prior assumptions regarding the failure surface in SRM. The study is carried out for different combinations of (i) slope angles (β), (ii) strength properties of the top and the bottom layer (c, ϕ) and (iii) different thickness of the top layer. The failure patterns are shown for a few cases.

slopes by applying the variational method and the strength reduction technique considering pseudostatic analysis.
The stability of layered slopes was also studied rigorously especially in the past few years [4, 5, 7, 17, 18, 24, 27-29, 37, 42]. By assuming the log-spiral failure mechanism, Chen et al. [7] evaluated the upper bound stability of non-homogeneous slopes corresponding to varying cohesion with depth. Kumar and Samui [24] evaluated the stability of layered soil slopes subjected to pore-water pressure and horizontal seismic force by using the rigid block upper bound limit analysis. Hammouri et al. [17] carried out the stability analysis of the layered slopes by considering the effects of rapid drawdown and tension cracks with the aid of PLAXIS 8.0 and SAS-MCT 4.0 software. By using the finite element limit analysis, [27,28], Shiau et al. [42], and Qian et al. [37] analyzed the undrained stability of the non-homogenous cohesive soil slopes. Chang-yu et al. [18] considered rotational mechanism with logarithm helicoids surface and assessed the 3D stability of non-homogeneous slopes. Lim et al. [29] proposed stability charts for frictional fill material placed on purely cohesive soil by using finite element LB limit analysis method. By using SRM, Chatterjee and Krishna [5] analyzed two and three-layered soil slopes considering a fixed slope angle (26.57°) and different combinations of three different chosen soils. The literature review clearly indicates that the rigorous analysis for the two-layered cohesive-frictional soil slope is quite limited. Although a few stability studies [24,38,40,41] were previously carried out by considering weaker layer overlying on strong layer, however, as per the authors' findings except the work of Sazzad et al. [40] hardly any study seems to be available for the case where strong layer is considered to be placed over weak stratum. The work of Sazzad et al. [40] also pertains to a specific combination of layered system (Top Soil: c 1 = 10 kPa, φ 1 = 18° and Bottom Soil: c 2 = 6 kPa, φ 2 = 10°). Hence, there is a requirement to carry out an extensive and rigorous investigation to estimate the improvement in stability by placing a stronger layer over a weaker layer. This is the prime motivation to carry out the present work. In this present article, strength reduction method is employed to analyze the two-layered soil slopes and to determine the factor of safety. The factor of safety is obtained for different combinations of (i) slope geometry (i.e. slope angle, β), (ii) strength properties of the top (c 1 , φ 1 ) and bottom layer (c 2 , φ 2 ) (iii) and the thickness of the top layer (t). The effect of placing different stronger layers over the weaker bottom layers is thoroughly investigated.

Strength reduction method (SRM)
The work of Zienkiewicz et al. [52] appears to be the first where SRM was used to solve the slope stability problem. Following his work, many further studies [5,9,10,13,15,16,20,26,30,31,34,35,40,45,47,50,51] were performed for slope stability problems by using SRM. This method is mainly based on the frame work of finite element method (FEM) and hence, all the advantages of FEM are retained in this method. The significant advantages of this method are the following: (i) the method is suitable to apply for complex geometries, complicated boundary and loading conditions, and (ii) there is no need to consider any assumptions regarding interslice shear forces and the critical failure surface. The critical failure surface is obtained automatically from the shear strength reduction. The additional information regarding stresses, strains, and pore pressures can also be obtained from this method. SRM is also applied on the basis of finite difference method [10]. Generally, SRM is applied to determine the factor of safety by successively reducing or increasing the shear strength of the material until the slope reaches the limiting equilibrium state.
The SRM is commonly used with the linear Mohr-Coulomb criterion where the failure strength is characterized by cohesion (c) and the internal friction angle (φ). The Mohr-Coulomb model is expressed in Eq. (1).
where τ is the maximum amount of shear stress the soil can resist for a certain applied normal stress (σ n ). The analysis is carried out by reducing the strength parameters (c,φ) progressively until the slope becomes unstable. In conventional SRM, both parameters are reduced by the same factor, or in other words, the reduction path of the cohesion and the friction are identical. The reduced cohesion and the friction angle are computed from Eq. (2).
where (i) c r and φ r are the reduced strength parameters and (ii) F s is the strength reduction factor. These reduced parameters are then reinserted into the model till the failure occurs. The main objective of SRM approach is to compute the strength reduction factor and the reduced material parameters that lead to collapse state.
In the present analysis, Optum G2 [36] is used for estimating the factor of safety of the slope through strength reduction method. OptumG2 is a finite element limit analysis (FELA) based software developed by OptumCE. For obtaining the limiting solutions, OptumG2 uses second order cone programming to solve the plane-strain stability problems. This scheme works by infeasibility detection in a very controllable way [22,43]. Following steps are adopted for the formulation: Step 1: Assuming F min and F max ; where, F min and F max are the minimum and maximum value of factor of safety. Generally, F min is chosen to be zero and F max is taken to be a large number within the range of machine precision.
Step 2: Initializing the value of F s and computing reduced strength parameters with the help of Eq. (2).
Step 3: Checking feasibility through the interior point method by using the reduced strength parameters.
Step 4: If the problem is feasible, assign F min = F S and evaluate a new factor of safety by using the harmonic mean as depicted in Eq. (3).
Otherwise, if the problem is infeasible, assign F max = F S and evaluate a new factor of safety by following the arithmetic mean as expressed in Eq. (4).
Step 5: Continue the iterative process (Step 1-Step 4) until the following convergence condition as mentioned in Eq. (5) is fulfilled: here, the tolerance limit, T L , is kept as 0·01. It is worth mentioning that based on the solution process, F min and F max provides the limiting extremities of the bound theorem. F min and F max represent rigorous lower and upper bound on the factor of safety corresponding to the statically admissible stress field domain and kinematically admissible velocity field domain, respectively. The numerical values presented herein are the average of both the limiting values. Figure 1 shows a two-layered soil slope having an angle, β. The strength parameters of the top and the bottom layer are represented by c 1 , φ 1 and c 2 , φ 2 , respectively. Slope height (H) in the present analysis is taken as 20 m for representing the high cut slopes. With the aid of strength reduction method, it is intended to analyze and subsequently compute the factor of safety for three different two-layered slopes (namely, 25°, 35° and 45°) consisting of different soil materials.

Problem statement and methodology
For performing the analysis, the size of the domain is considered adequately high so that the failure surface remains contained well within the domain. Based on trials, the height (D) and length (L) of the domain are kept as 2H and 9H, respectively. The boundary conditions are mentioned in Fig. 1. Vertical and horizontal displacements are restrained along the base of the considered domain. Along the left and right boundaries, horizontal displacement is not allowed to occur. The soil mass is discretized by using three nodded linear triangular elements. The soil plasticity is governed by the Mohr-Coulomb failure criterion and associated flow rule. Adaptive mesh refinement based on plastic shear dissipation has been used. Three iterations of adaptive meshing with 10,000 elements have been considered for all analyses. A nonlinear optimizer named sonic is used in Optum G2 software for optimization.

Results and discussions
In the present work, solutions are presented in terms of factor of safety for different combinations of the (a) slope angles (β), (b) soil strength properties of the top and the bottom layers (c 1 , φ 1 and c 2 , φ 2 ) and (c) top layer thickness (t). Stability charts are presented in tabular form. Tables 1, 2, 3 show the factor of safety values computed for three different two-layered slopes (β = 25°, 35° and 45°). In this article, nine different stronger soil layers [(35,0), (35,5), (35,8), (40,0), (40,5), (40,8), (45,0), (45,5), (45,8)] were considered to be placed over twelve different weaker bottom layer [ (20,5), (20,10), (20,15), (20,20), (25,5), (25,10), (25,15), (25,20), (30,5), (30,10), (30,15), (30,20)]; the first and second part within the parenthesis indicate the frictional (in degrees) and cohesive strength (in kPa), respectively. The thickness of the stronger layer was varied between 20%-80% of the domain height. Both the limiting values (lower and upper) were obtained. A total number of 3240 computations were performed. Following observations are made from the numerical results: (a) As expected, the factor of safety (F s ) decreases with an increase in slope angle (β). However, this decrement depends on the strength of the soil layers. For an example, as β varies from 25° to 45°, F s reduces markedly. This reduction differs by 14% (from 53 to 39%) as the cohesive strength of the top layer, (φ 1 = 35°, and, t/D = 0.4) which is rested upon a certain bottom layer (c 2 = 20 kPa, φ 2 = 25°), increases up to 8 kPa from 0 kPa. (b) Placing a stronger layer over a weaker stratum undoubtedly improves the stability of the slope and this improvement is more significant for steep slopes. For the previous example, (i) if the cohesive strength of the top layer rises from 0 to 8 kPa (keeping φ 1 equals to be 35°), the improvement in F s for 25° and 45° slope occurs by 19% and 56%, respectively; and, (ii) if the frictional strength increases from 35° to 45°, F s improves by 26% and 42% for 25° and 45° slope, respectively. (c) When the thickness of the top layer is within a certain limit, the strength of the bottom layer also influences the stability of the slope. There is almost a linear relationship between the improvement of F s with the increase in cohesive strength of the bottom layer. However, the relation between the improvement of F s with the increase in frictional strength of the bottom layer is highly nonlinear. (d) The tabulated data clearly reveals that the improvement in F s is quite significant as the top layer thickness changes from 0.2D to 0.4D. On the contrary, when t/D ratio varies from 0.6 to 0.8, the improvement in stability is almost negligible. The effect of the thickness of the top layer is further studied graphically. It is quite evident that higher the strength of the top layer higher would be the safety factor, F s . These figures depict that there is a certain t/D beyond which there is hardly any improvement in stability of the slopes. This particular top layer thickness is termed as optimum thickness Table 1 Proposed stability chart (indicating the factor of safety) corresponding to φ      Table 3 Proposed stability chart (indicating the factor of safety) corresponding to φ 2 = 30°φ   and is referred here as dimensionless parameter, 't opt /D' . The value of t opt /D increases with the increase in the strength of the top layer. The figure shows that the dependence of t opt /D on the cohesive strength of the bottom layer is further influenced by the frictional strength of the top layer; when φ 1 = 40°, t opt /D decreases with increase in c 2 however, when φ 1 = 45° there is no impact of c 2 on the computed value of t opt /D. Figure 3 shows the variation of F s with t/D for β = 25° and 45°, corresponding to two different φ 2 , namely, 20° and 30°. The properties of the top layer are kept to be constant and the cohesive strength of the bottom layer is varied within the range of 5-20 kPa. The figures clearly reveal that for the same soil properties, optimum thickness of the top layer is significantly smaller for the steeper slopes. As the frictional strength of the bottom Fig. 2 The variation of F s with t/D for a two-layered slope (β = 25°) corresponding to varying c 2 and (a) layer increases, the magnitude of t opt /D further reduces. The numerical solutions give an impression that the impact of the strength of the top layer on the stability is much higher than the strength of the bottom layer. Figure 4 shows the mesh pattern at the collapse state for three different slope angles, namely, 25°, 35° and 45°. The soil profiles for these three cases are kept to be the same. It is to be noted that adaptive mesh refinement technique continuously updates the sizes of all the elements in an optimal fashion by computing the variations of stresses and velocities. Finer elements were automatically placed in the shear failure zone. Hence, these meshes indirectly depict the failure patterns. The figure shows that the size of the failure zone decreases with the increase in slope angle. Moreover, as the slope angle increases the failure is likely to become toe failure. This observation is in accordance with the studies of Lim et al. [29] and Sazzad, et al. [40] who had earlier reported that if the top layer is considered to be stronger than the bottom layer and the slope angle is considered to be less than equal to 45° the incipient state of collapse in the soil slope will be triggered by developing base failure. As the steepness of the slope increases, the extent of the failure zone seems to be restricted closer to the slope surface. Figure 5 illustrates failure state corresponding to three different thickness of the top layer. The soil properties of the top layer as well as the bottom layer are the same for all the three cases. The figure demonstrates that as the thickness of the top layer increases the type of failure surface turns from toe to base. However, beyond a certain thickness, the slope collapses by developing the toe failure surface and the shear zone seems to be confined within the top layer. This observation substantiates the existence of t opt /D. Figure 6 depicts the mesh pattern at the collapse state corresponding to three different frictional angle of the top layer. All other geometrical and material strength parameters are kept to be the same. The figure illustrates that as the frictional strength of the top layer increases the failure zone grows in size. However, the extent at which the finer elements are laid at the collapse state goes thinner with the increase in φ 1 . It gives an impression that the thickness of the shearing zone (i.e. shear band) becomes smaller with the placement of stronger layer over a weaker stratum.

Comparison of results
Comparisons of both the limiting solutions, for the homogenous and layered slopes, are presented in Tables 4 and 5, respectively. In most of the cases the difference seems to be in the second decimal place. Closeness of the lower and upper bound solutions further depicts the accuracy in the computed solutions. Limit theorems suggest that the true solution will lie somewhere between these bounding values. It should be recalled that the safety factors charts presented in Tables 1, 2, 3, are the average value of the two extremities. Table 6 shows the comparison of the present solutions computed with the numerical results provided by Dawson et al. [10] for a homogenous slope of 10 m height having unit weight of soil, γ = 20 kN/m 3 and cohesion, c = 12.38 kPa. Dawson et al. [10] had employed strength reduction method by using the explicit finite difference code, FLAC. As the frictional strength of the soil increases the present method provides higher stability value. For the same soil, the difference between these two solutions reduces as the steepness of the slope increases. The reason may be attributed  [10] had discretized the chosen domain with four nodded rectangular elements, whereas, in the present work three nodded linear triangular elements are used. The same trend is also observed while the solutions of Dawson et al. [10] are compared with the upper bound solutions (assuming log spiral mechanism) obtained by Chen [6]. It is to be noted that the present finite element limit solutions are quite smaller than those rigid block upper bound solutions provided by Chen [6]. It shows the improvement of the solutions when finite element limit analysis is employed for the analysis. Table 7 shows the comparison of the present solutions with the results provided by Kumar and Samui [24] by using the rigid block upper bound method considering logspiral failure mechanism. The comparison is carried out for 45° slope corresponding to different soil layer properties and varying top layer thickness. Similar to the previous observation, it is well noted that the present solutions become quite smaller than the reported solutions of Kumar and Samui [24] as the strength of the soil layer increases. Table 8 illustrates the comparison of current solutions with the solutions provided by Chatterjee and Krishna [5] for non-homogeneous slopes. Chatterjee and Krishna [5] Table 9 depicts the comparison of present solutions with limit equilibrium solutions presented by Sazzad et al. [41] for layered soil slopes. Sazzad et al. [41] used Bishop Method [3] for LEM analysis. The present solutions appear to be smaller than the LEM solutions. This can be attributed to the fact that LEM solutions generally overestimate the factor of safety due to the usage of statical and kinematical assumptions. This is also observed in earlier studies as well [48].

Conclusions
In the present article, the stability of two-layered soil slopes is analyzed by using strength reduction method. A series of upper and lower bound limit analyses are carried out in Optum G2 software by placing different stronger layers of varied thickness over the weaker stratum. Stability charts are prepared in the form of the factor of safety for different soil properties, slope geometries and top layer thickness. The amount of improvement in the stability by placing a layer of stronger soil over the weaker stratum is numerically investigated. The optimum thickness of the top layer is reported for a wide range of slopes. The extent and the type of the failure zones are presented for several cases. The obtained solutions compared quite well with the available solutions. The proposed design charts may be useful to the practicing engineers.