#### Abstract

Axisymmetric concave slopes, one special type of three-dimensional (3D) slopes, may be encountered in mining and civil engineering practice. Analysis of 3D slopes is generally complex and mostly relies on complicated numerical simulations. This paper proposes an elastoplastic solution for determining the additional shear resistances due to spatial effects of axisymmetric concave slopes. By incorporating the extra antislide forces, this paper proposes a simplified two-dimensional (2D) limit equilibrium procedure for the stability analysis of axisymmetric concave slopes. Combined with an iteration algorithm, the procedure can obtain the factors of safety for axisymmetric concave slopes in a simple and efficient way. Comparisons of the results from the proposed method and the numerical software FLAC^{3D} are performed to demonstrate the validity of the proposed method for practical applications. Finally, the effects of several key parameters on the stability of axisymmetric concave slopes are investigated through a parametric study.

#### 1. Introduction

Axisymmetric concave slopes, one special type of three-dimensional (3D) slopes, may be encountered in mining and civil engineering projects, such as open-pit excavations and earth-fill cofferdams. These slopes have a round or approximately round shape in a plan view but a 3D geometry in nature. Because of spatial 3D effects of stresses in the slope body, axisymmetric concave slopes cannot be directly treated as infinite long straight slopes with a uniform cross section; hence, typical two-dimensional (2D) limit equilibrium procedures are no longer applicable for axisymmetric concave slopes in a straightforward way. For this reason, further research is needed to investigate the stability of axisymmetric concave slopes.

In the past, 2D limit equilibrium methods (LEM) have been commonly used for stability analyses of slopes in geotechnical engineering (e.g., [1–8]). In a typical 2D limit equilibrium method, a slip surface is often assumed *a priori* and the failure body, encompassed by the slip surface and the slope surface, is then divided into slices of soil mass. After that, force equilibrium and/or moment equilibrium conditions are established to solve the factor of safety (FoS), which is typically defined as the ratio of the shear strength (or resisting moment) to the shear stress (or driving moment).

Considering that the actual failure surface of a slope is generally 3D in nature, a few researchers carried out stability analyses of slopes based on limit equilibrium approaches with 3D failure surfaces, e.g., Leshchinsky et al. [9], Lam and Fredlund [10], Chen et al. [11], and, more recently, Jiang and Zhou [12]. These studies assumed straight slopes with 3D failure surfaces (cylindrical, spherical, or others). However, as compared to straight slopes, axisymmetric concave slopes have a round or approximately round shape in the plan view, which has apparent spatial effects of stresses in the slope body. Ignoring these effects may lead to inaccurate or too conservative calculation of factors of safety for axisymmetric concave slopes.

Numerical simulations are considered to be a powerful tool for analyzing the spatial (or 3D) effects developed in a concave slope. For example, Lorig [13] evaluated the stability of concave slopes using the numerical software FLAC^{3D}, and Sun et al. [14] employed the displacement finite element software ABAQUS to develop stability charts for convex and concave slopes. In addition to numerical simulations, the method of characteristics has been used to analyze 3D stability of concave slopes, e.g., Jenike and Yen [15] and Jahanandish and Keshavarz [16]. However, these studies assumed homogeneous concave slopes.

Although 3D limit equilibrium approaches or 3D numerical simulations may provide more reasonable results for stability analyses of slopes than 2D limit equilibrium approaches, building 3D models is generally time-consuming and usually requires special expertise of designers. In contrast, 2D limit equilibrium approaches are preferable for a preliminary analysis of slope stability in engineering practice because of their simplicity and long history of use. Generally, the FoS obtained by the 2D limit equilibrium approach is smaller than that based on 3D models. However, if additional shear resistance due to spatial effects of 3D slopes can be considered in the analysis, the results of the FoS based on the 2D limit equilibrium approach will be close to the actual value.

Based on a 3D limit equilibrium approach, Zhang [17] proposed a practical method for 3D stability analysis of concave slopes in the plan view. In his method, the failure mass was divided into *n* vertical columns and the force and moment equilibrium conditions of the failure mass were established to calculate the FoS. The most important contribution of this work is that it considered the so-called “end force” *P* caused by the lateral pressure of soil acting on individual columns, i.e., accounted for the spatial (3D) effects of a concave slope. Because the end force *P* (normal to the radial line) is not perpendicular to the sliding direction, it provides additional shear resistance on the failure mass. Zhang [17] calculated the end force *P* empirically using the following equation:where *γ* is the unit weight of soil, *h* is the height of a column, *K*_{a} is the active earth pressure coefficient, and *V*_{xz} is the projection of the base of each column *V* on the *X*-*Z* plane.

In the authors’ opinion, the use of active earth pressure for the approximation of the end force *P* may be too conservative in practical applications. This is because, when the concave slope fails, the sliding columns squeeze each other in the circumferential direction and the soil pressure in the circumferential direction will be much greater than active and static earth pressures. Therefore, a better method concerning the calculation of the soil pressure in the circumferential direction needs to be proposed for the stability analysis of concave slopes, which constitutes a major objective of this paper.

Based on the above discussions, an elastoplastic solution is proposed in this study to give more reasonable estimations of the additional shear resistance due to the spatial effects of an axisymmetric concave slope. For this purpose, the concave slope is divided into a series of arch-shaped slices in the plan view (or thick cylinders in the space view). Then, the earth pressures acting on inner and outer faces of an individual cylinder are simplified as axisymmetric linearly distributed loads with depth, and the surcharge load on the slope surface is simplified as a uniformly distributed load acting on the upper face of the cylinder. Based on these simplified boundary conditions, an elastoplastic solution is deduced for the distributions of stresses in the cylinder. Finally, the solution of the circumferential stress is used to calculate the additional shear resistance acting on individual column elements for a 2D limit equilibrium analysis for concave slopes. With the aforementioned elastoplastic solution, a simplified 2D limit equilibrium procedure for the stability analysis of concave slopes is proposed later in this paper. The proposed method is verified by comparing the calculated results with those from the numerical software FLAC^{3D}.

#### 2. Problem Descriptions

##### 2.1. Geometry of a Concave Slope

Figure 1 shows the geometry parameters of a typical axisymmetric concave slope with horizontal layers including (1) the slope height (*H*), (2) the slope angle (*β*_{c}), and (3) the radius at the toe of the slope (*R*_{c}). In order to conduct a 3D limit equilibrium analysis, the concave slope can be divided into *m* × *n* individual vertical columns. This can be easily done by equally dividing the concave slope into *m* cylinder slices (in the plan view) along the radial direction and *n* wedge-shaped blocks along the circumferential direction. In the plan view, the overlapping area of the *i*-th cylinder slice and the *j*-th wedge-shaped block is designated as column (*i*, *j*), which is a fundamental element for the 3D limit equilibrium analysis of an axisymmetric concave slope.

##### 2.2. Definition of Additional Shear Resistance

Figure 2(a) shows a typical column element (*i*, *j*) taken from an axisymmetric concave slope. Prior to defining the additional shear resistance, it is necessary to analyze all the internal and external forces acting on this element. The 3D limit equilibrium analysis of a general axisymmetric slope includes the following forces: (1) the normal force *N*_{i} and the shear force *S*_{i} at the base of the column; (2) the normal force *P*_{ui} due to pore-water pressure at the base of the column; (3) the external vertical forces *W*_{i} and *Q*_{i} due to soil weight and surcharge load of the slope, respectively; (4) the intercolumn normal forces *E*_{ri−1} and *E*_{ri} in the *r*-direction and *E*_{θi} in the *θ*-direction; and (5) the vertical intercolumn shear forces *X*_{ri−1} and *X*_{ri} in the *r*-direction. Note that the subscript *j* for all these quantities has been omitted considering the inherent axisymmetry of the studied problem.

**(a)**

**(b)**

The forces *N*_{i}, *S*_{i}, *P*_{ui}, *W*_{i}, *Q*_{i}, *X*_{ri−1}, *X*_{ri}, and *E*_{ri} play the same roles in a limit equilibrium analysis for both a straight slope and a concave slope. Therefore, the only difference between the 3D limit equilibrium analyses of a straight slope and a concave slope is that the force *E*_{θ} is not perpendicular to the sliding direction and can provide additional shear resistance on the column element (*i*, *j*). The additional shear resistances account for the spatial (3D) effects of the stresses in a concave slope and hence should be considered in practical applications.

Figure 2(b) shows the column element (*i*, *j*) and the corresponding intercolumn forces in the plan view. Figure 2(b) shows that the intercolumn forces *E*_{θi} can be decomposed into the component forces *E*_{θri} in the direction of sliding and the component forces *E*_{θθi} perpendicular to the direction of sliding. If the central angle *dθ* of the column element (*i*, *j*) is sufficiently small, the following relation can be easily obtained:and the resultant force *P*_{ri} of the component forces *E*_{θi} acting on two sides of the column element (*i*, *j*) can be calculated as

Based on the above discussions, the resultant force *P*_{ri} acting in the opposite direction of movement is defined as the additional shear resistance for the specific column element (*i*, *j*). Because of the presence of the additional shear resistances, a concave slope is always more stable than a straight slope.

#### 3. Determination of Additional Shear Resistances

##### 3.1. Assumptions for Simplification

As discussed previously, a key issue for the stability analysis of an axisymmetric concave slope is the determination of the additional shear resistances. Since the shear resistances are *in fact* generated by the lateral pressures (stresses) in the circumferential directions of an axisymmetric concave slope, it is necessary to determine the inner circumferential stresses developing in the concave slope so as to more reasonably estimate the additional shear resistances. For this purpose, the axisymmetric concave slope is simplified as a set of thick cylinders (or cylinder slices in the plan view) of different inner and outer diameters. Based on the boundary conditions, the distributions of stresses in these cylinders can be derived via the elastoplastic theory.

Figure 3 shows the cross-sectional profile of the *i*-th thick cylinder element in an axisymmetric concave slope, which is also referred to as the *i*-th cylinder slice. To be clear, the following symbols are given: the height of the *i*-th cylinder is *h*_{i}; the inner and outer diameters of the cylinder slice are *r*_{i−1} and *r*_{i}, respectively; the surcharge load applied on the top of the *i*-th cylinder slice is *q*_{i}; the average unit weight of the cylinder is *γ*_{i}; the soil pressures acting on the inner and outer faces of the cylinder are *p*_{i−1} and *p*_{i}, respectively; and the vertical shear stresses acting on these two faces are *τ*_{i−1} and *τ*_{i}, respectively. Before an elastoplastic solution is deduced for the distributions of stresses in the cylinder, some assumptions for simplification of the boundary conditions are proposed:(1)The soil pressures (normal stresses), *p*_{i−1} and *p*_{i}, acting on the inner and outer faces of the cylinder are axisymmetric and linearly distributed with depth.(2)The vertical shear stresses, *τ*_{i−1} and *τ*_{i}, acting on the inner and outer faces of the cylinder, are irrelevant to the calculation of the distribution of the circumferential stress *σ*_{r}, which will be demonstrated later in this paper.(3)The thickness of the cylinder element is assumed to be sufficiently small; hence, the surcharge load *q*_{i} on the top of the cylinder can be treated as a uniform vertical load.

##### 3.2. Elastoplastic Solution for Additional Shear Resistances

Based on the simplification assumptions, an elastoplastic solution is developed for the distributions of stresses in the cylinder element in this section. Then, the solution for the circumferential stress *σ*_{r} is used to calculate the additional shear resistances *P*_{ri} for a specified vertical column element (*i*, *j*) in the limit equilibrium analysis.

The stress equilibrium equation for the *i*-th cylinder can be written asand the simplified boundary conditions are expressed aswhere *k*_{i−1} and *k*_{i} are the gradients of the soil pressures acting on the inner and outer faces of the cylinder, respectively; *τ*_{i−1} and *τ*_{i} are the vertical shear stresses acting on the inner and outer faces of the cylinder, respectively; and *q*_{i} is the vertical surcharge load acting on top of the cylinder. For a specific column element (*i*, *j*), these quantities can be calculated using the following relations:where and are the normal forces acting on the inner and outer faces of the element column (*i*, *j*); and are the vertical shear forces acting on the inner and outer faces of the element column (*i*, *j*); is the vertical surcharge load acting on the upper surface of the column element (*i*, *j*); *b*_{i} is the thickness of the column element (*i*, *j*); and *dθ* is the central angle of the column element (*i*, *j*).

According to the elastoplastic theory, the solution of equation (4) can be decomposed into a characteristic solution without a gravity force and a special solution with a gravity force. Therefore, a general solution is derived by the following equation:

Based on Love’s method [18], the stress components in equation (11) can be expressed aswhere are the stress components; is a potential function; and is Poisson’s ratio. The potential function must satisfy the following biharmonic condition:where is the Laplacian operator.

For a spatial axisymmetric problem, the potential function (*ϕ*) can be defined aswhere *A*_{1} to *A*_{7} are unknown constants. Combining equations (12)–(17) with the stress boundary condition (5), the constants *A*_{1} to *A*_{7} can be derived as

By substituting the potential function *ϕ* in equation (17) and the constants in equations (18)–(24) into equations (12)–(15), the stress components can be further expressed as

From equations (25) and (26), it can be easily seen that the shear stresses *τ*_{i−1} and *τ*_{i} are irrelevant to the stresses and , which is consistent with Assumption (3) in Section 3.1.

Consider the field of gravity as a special solution of equation (4), the vertical component of which isand the general solution for the stress *σ*_{z} can be expressed as

Considering the soil satisfies the Mohr–Coulomb failure criterion, the failure condition can be expressed aswhere and denote the effective cohesion and effective friction angle of the soil at depth *z*, respectively.

From equations (25) and (26), it can be seen that a smaller value of *r* results in a lower value of and a higher value of . Based on this fact and equation (31), the plastic zone first occurs at the radial distance of *r* = *r*_{i−1}. Therefore, on the incipient failure of the concave slope, the plastic stresses and can be deduced from equations (25) and (26) at the radial distance of *r* = *r*_{i−1} and be further simplified into

Substituting equations (32) and (33) into equation (31) leads to

After some necessary manipulations, equation (34) can be transformed into the following relation:

Finally, substituting equation (35) into equation (32) leads to

Note that equation (36) is a linear expression with respect to depth *z*. Therefore, the lateral forces *E*_{θi} acting on two sides of the column element (*i*, *j*) can be calculated after combining equation (36) with equation (6):and the force acts on 1/3 the height of the column element (*i*, *j*).

Finally, considering equation (3), the previously defined additional shear resistance *P*_{ri} can be calculated as

#### 4. Simplified 2D Limit Equilibrium Procedure for Stability Analysis of Concave Slopes considering Additional Shear Resistance

##### 4.1. Equilibrium Conditions

A variety of 2D LEMs can be used to estimate the slope stability, such as Fellenius’ method, Bishop’s simplified method, Janbu’s simplified method, Spencer’s method, and Morgenstern–Price’s method. By incorporating additional shear resistance measured by equation (38) with any of these 2D LEMs, the stability of concave slopes can be assessed. In this section, Bishop’s method is taken as an example to illustrate the simplified 2D limit equilibrium procedure.

Figure 2 illustrates that a typical concave slope can be equally divided into *n* wedge-shaped blocks along the circumferential direction. Figure 4 shows the *j*-th wedge-shaped block extracted from the concave slope, which can be subsequently divided into *m* column elements.

With this configuration, the equilibrium conditions of forces and moments for the wedge-shaped block are established as follows.

After a failure surface is assumed, the weight of the column element (*i*, *j*) can be expressed aswhere is the inclination angle of the bottom face of the column element (*i*, *j*) to the horizontal plane and is the length of the bottom face of the column element (*i*, *j*).

The normal force caused by pore-water pressure at the base of the column element (*i*, *j*) is calculated aswhere *U*_{i} is the average pore water pressure.

The shear resistance due to soil cohesion and friction at the base of the column element (*i*, *j*) is calculated aswhere *F* is the factor of safety, denotes the effective friction angle of the soil in the failure surface of the column element (*i*, *j*), and is the shear resistance due to soil cohesion, which can be calculated aswhere denotes the effective cohesion of the soil in the failure surface of the column element (*i*, *j*). In equation (41), is the normal force acting at the base of the column element (*i*, *j*). The quantity can be determined by the force equilibrium in the *z*-direction as follows:

As in the Bishop algorithm, it is assumed that

Hence, combining equations (42)–(44) leads to the expression of *N*_{i} as follows:where the parameter *m*_{αi} can be calculated by

Substituting equation (45) into equation (41) leads to

Next, the force equilibrium in the radial (*r*) direction leads towhere and are the lateral forces acting on the inner and outer faces of the column element (*i*, *j*) and *P*_{ri} is the previously defined additional shear resistance, which acts at the height of *h*_{i}/3 from the bottom face of the column element (*i*, *j*). Note the innermost lateral force *E*_{r0} at the toe of the concave slope is zero.

Finally, the overall moment equilibrium of the failure body leads towhere *R* is the distance from the bottom face of the column element (*i*, *j*) to the rotation centre *O* and *R*_{ri} is the distance from the action line of the additional shear resistances *P*_{ri} to the rotation centre *O*, which is calculated as

##### 4.2. Iteration Algorithm for Solving Critical FoS

To solve the critical value of the FoS efficiently, an iteration algorithm is proposed in this section as follows:*Step 1.*A potential failure surface is assumed and the *j*-th wedge-shaped block is subdivided into a series of column elements, as shown in Figure 4.*Step 2.*The initial values of the additional shear resistances are set to zero, i.e., *P*_{ri} = 0, and an initial estimate of *F* is given, which is typically set to 1.*Step 3.*As in Bishop’s algorithm, the value of *F* is updated repeatedly using equations (39) to (40), (42), and (45)–(50) until the difference between the values of *F* obtained by the last two iterations is smaller than a certain tolerance (e.g., ).*Step 4.*The values of are extracted from the analysis in Step 3, and then these values are substituted into equation (38) to update the values of the additional shear resistance *P*_{ri,j}.*Step 5.*Based on the updated values of *P*_{ri,j}, repeat Steps 3 and 4 until the values of the additional shear resistances *P*_{ri,j} converge to fixed values as the Euclid norm ||*P*_{ri,j}|| is less than a specific tolerance (e.g., ).*Step 6.*The minimum value of *F* is found by repeating Steps 1 to 5 for all possible failure surfaces.

The minimum value of *F* determined based on the above-mentioned steps is taken as the critical FoS for a concave slope.

It is worth mentioning that the proposed method can be seen as a generalized traditional 2D limit equilibrium approach with the additional shear resistances considered in the analysis. The following section will demonstrate that this approach is simple and efficient for practical uses.

##### 4.3. Verification of the Proposed Method

Based on the strength reduction technique incorporated in the numerical software FLAC^{3D}, Zettler et al. [19] analyzed two 25-meter high slopes: (1) an axisymmetric concave slope with a 12-meter radius of curvature at the toe of the slope and (2) a straight slope. Table 1 lists the material and geometry parameters of the concave slope. For comparison, Table 2 lists the calculated results of the FoS obtained by Zettler et al. [19] and the proposed method, together with those from the Bishop algorithm for the straight slope.

Table 2 shows that the results of the proposed method and Bishop’s method agree well for the straight slope. This is because no additional shear resistances are actually considered in the proposed method for a straight slope. In this case, the proposed method is basically the same as Bishop’s method. For the straight slope, the FoS calculated by FLAC^{3D} using the strength reduction technique is 1.37. The difference in the results (FoS) between the proposed method and FLAC^{3D} is less than 8%, which is small considering that these two methods are quite different. For the concave slope, the results of the FoS from the proposed method and FLAC^{3D} are 1.94 and 1.83, respectively. Their difference is also less than 8%. The above discussions indicate that the proposed method is applicable for analyzing concave slopes in practical applications.

#### 5. Parametric Study

This section presents a parametric study performed to investigate the influence of some important parameters on the stability of a concave slope, including the ratio of the radius of curvature at the toe of a concave slope to the slope height *R*_{c}/*H*, the slope angle *β*_{c}, and the strength parameters soil cohesion and soil friction angle .

Figure 5 shows the variations of the calculated FoS with respect to the ratio *R*_{c}/*H*. For a comprehensive comparison, the results of the FoS obtained from FLAC^{3D} are also shown in Figure 5. It needs to be mentioned that other necessary parameters used in FLAC^{3D} are the same as those in Table 1. Figure 5 shows that the maximum error between the calculated FoS using the proposed method and that obtained by FLAC^{3D} is 9.7%, and the error decreases with an increase of *R*_{c}/*H*. Considering the discrepancy between the two approaches, the difference in results is acceptable.

Figure 5 also shows that when *R*_{c}/*H* is greater than 10, the FoS calculated by the proposed method gradually decreases to the FoS given by Bishop’s method (1.274), indicating that a concave slope can be treated as an infinite long straight slope as long as the ratio *R*_{c}/*H* is sufficiently large.

Figure 6 shows the critical failure planes predicted by the proposed method with different *R*_{c}/*H*. As the *R*_{c}/*H* ratio decreases, the position of the critical failure planes moves upwards and their corresponding FoS increases. This phenomenon may be explained by the fact that the hoop stress effect is more apparent for a concave slope with a higher *R*_{c}/*H*, and hence, the stability of a concave slope is enhanced.

To investigate the effects of the slope angle *β*_{c} and the strength parameters and , a conversion factor *C*_{f} is defined as follows:where FoS_{c} and FoS_{s} correspond to the factors of safety for a concave slope and a straight slope with the same characteristics, except for the analyzed parameter.

Figure 7 shows the variation of the conversion factor *C*_{f} with respect to *R*_{c}/*H* for different slope angles *β*_{c}. It is clearly shown that a concave slope is more stable than a straight slope because *C*_{f} is always greater than 1.0. However, the contrast of stability between a concave slope and a straight slope becomes more apparent for steeper slopes (with higher values of tan *β*_{c}), especially when *R*_{c}/*H* is less than 0.3. Also, it is interesting to note that the slope angle *β*_{c} has an opposite impact on the factor *C*_{f} for the cases with *R*_{c}/*H* smaller or greater than 0.6, indicating that an obvious combined effect exists between the parameters *β*_{c} and *R*_{c}/*H*. However, the reason for this phenomenon remains to be further discussed.

Figures 8 and 9 show the variations of *C*_{f} with respect to *R*_{c}/*H* for different soil cohesions and friction angles, respectively. Figure 8 shows that the factor *C*_{f} increases as the soil cohesion increases, implying that the spatial effect in a concave slope is more apparent when the soil cohesion has a higher value. In contrast, Figure 9 shows that the factor *C*_{f} decreases as the soil friction angle increases, implying that the spatial effect in a concave slope weakens as the soil friction angle increases.

#### 6. Conclusions

Because of its simplicity and efficiency, the 2D limit equilibrium method is still the most widely accepted approach for slope stability analysis in practice. However, the use of a 2D limit equilibrium method for stability analysis of concave slopes may lead to an inaccurate estimate of the FoS unless the additional shear resistances are considered in the analysis.

This paper proposes an elastoplastic solution for calculating the additional shear resistances due to the spatial effects of stresses in an axisymmetric concave slope. Considering the additional shear resistances, a 2D limit equilibrium procedure combined with an efficient iteration algorithm is developed to calculate the FoS of concave slopes.

The comparison of the calculated FoS from the proposed method and the numerical software FLAC^{3D} demonstrated that the proposed method is accurate and efficient for stability analyses of concave slopes. Based on the proposed method, a parametric study was performed to study different effects of the parameters on the stability of concave slopes, indicating that the 3D effects of stresses must be considered in the 2D analysis to give a reasonable estimate of the FoS of axisymmetric concave slopes.

#### Data Availability

The data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

The authors would like to thank the financial support of the National Natural Science Foundation of China (NSFC) (No. 51578230) for this work. Also, the authors appreciate the generous help from Professor Changfu Chen of Hunan University during the preparation of this research manuscript.