Limit analysis can be classified into two categories: upper bound analysis and lower bound analysis. The former is particularly effective in analyzing the failure mechanism of slopes (Sloan, 2013). According to the upper bound theorem, the limit load, which is the minimum value among all the loads corresponding to admissible velocity fields, can be mathematically expressed as Eq. 1: ∫ S T i v i d S + ∫ V f i v i d V ≥ ∫ V σ i j * ε i j * d V (1) where S and V are the surface area and volume region, respectively; T _i and fi are surface forces and body forces, respectively; v _i is the admissible velocity corresponding to the failure mechanism; and σ _ij and ε _ij are stress and strain, respectively.

The yield criterion can be described using the Mohr-Coulomb model and considering non-flow criterion for soils, which can be expressed as Eqs 2, 3: Y = σ 1 − σ 3 + σ 1 + σ 3 sin ⁡ ϕ * − 2 c * ⁡ cos ⁡ ϕ * (2) tan ⁡ ϕ * = cos ⁡ ϕ ⁡ cos ⁡ ψ 1 − ⁡ sin ⁡ ϕ ⁡ sin ⁡ ψ tan ⁡ ϕ c * = cos ⁡ ϕ ⁡ cos ⁡ ψ 1 − ⁡ sin ⁡ ϕ ⁡ sin ⁡ ψ c (3) where Y is the yield function; σ ₁ and σ ₃ are the maximum and minimum stresses, respectively; ψ is the dilation angle; c and ϕ are the cohesion strength and friction angle, respectively; and c ^* and ϕ ^* are the cohesion strength and friction angle influenced by the dilation angle, respectively.

The factor of safety, F is expressed as Eq. 4: F = c c , = tan ⁡ ϕ tan ⁡ ϕ , (4) where c, and ϕ, are reduced cohesion strength and friction angle, respectively.

Although limit analysis method is a powerful method in geotechnical stability analysis as explained above, they still have limitations due to the need for assuming the shape of slip planes before calculation. To address this limitation, Sloan (2013) proposed the finite element limit analysis method, which combines the benefits of finite element method and limit analysis method to determine the safety factor through an automated search process. Furthermore, by incorporating mesh adaptive techniques into the calculation, it enables to identify slip plane directly. In the finite element limit analysis method, the object domain is discretized using nonlinear elements. Subsequently, the solution is transformed into an optimization problem, represented by the following (Eq. 5): min ⁡ imize  σ T B ¯ u − ∫ A t t T u d A − ∫ V g T u d V  power dissipation rate of work done done by fixed external forces subject to  B ¯ e u e = λ ˙ e ∇ f σ e  flow rule conditions λ ˙ e ≥ 0  plastic multiplier λ ˙ e f σ e = 0  consistency conditions A u = b  velocity conditions f σ e ≤ 0  yield conditions (5) where σ is a global vector of stresses; B ¯ is a global strain–displacement matric; u is a global vector of nodal velocities; t is a vector of known tractions; g is a vector of fixed body forces; B ¯ e is the strain-displacement matrix of element e; u e is the vector of nodal velocities of element e; λ ˙ e is the vector of plastic multiplier of element e; A is a metric of constants, b is a vector of constants; and f σ e is the yield condition of element e, respectively.

In this study, the computational analysis was conducted using the software OptumG2 (Krabbenhoft et al., 2015), where the crack was modeled by the Shear Joint element with smaller strength than adjacent soil. OptumG2 is a finite element program that offers various analysis types, including limit analysis and elastoplastic deformation analysis. The unique feature of OptumG2 lies in its implementation of the former analysis, which is based on the limit analysis theorem. This method provides a more efficient and convenient way for directly analyzing geotechnical stability, in contrast to the traditional elastoplastic methods and limit equilibrium methods. Moreover, OptumG2 enables the acquisition of rigorous lower and upper bounds, which help to bracket the true solution from below and above. In terms of modeling capabilities, the software facilitates easy and intuitive geometry modeling by automatically recognizing intersections, closed surfaces, and other relevant features. Additionally, OptumG2 offers a range of constitutive models and special elements specifically designed for modeling interfaces, reinforcements, and joints. For further information regarding OptumG2, more details can be referred to Krabbenhoft et al. (2015). For detailed calculation principles and software features, readers can refer to references Sloan (2013) and Xu et al. (2018).

In the following sections, the factor of safety, F, and the corresponding slip planes were determined using the flowchart illustrated in Figure 1.

Gao et al. (2005) reported centrifuge model tests results on a slope with three layers. The dimension of the slope and shear strength of soils are shown in Figure 2A. The slope model in the centrifuge test was accelerated to 60 g, and vertical cracks were observed at the top of the slop. The slip planes from the centrifuge model test and the proposed method in this study are shown in Figure 2A. Results in Figure 2A indicate that the shape of the slip planes from the proposed method are consistent with those from the model tests.

Chen et al. (2014) performed a numerical analysis on a slope with cracks using different limit equilibrium methods implemented in the Sliding software. The slope, as depicted in Figure 2B, has a height of 30 m and consists of three layers with parameters listed in Table 1. The maximum length of the vertical crack is 5 m. The cohesive strength and friction angle of the cracks were taken as 14 kPa and 4° in numerical analysis (Chen et al., 2014). The vertical crack located at the top of the slope with a depth of 5 m and a horizontal spacing of 1 m.

Computed results of safety factor, F from different methods for slopes with and without crack are listed in Table 2. Results in Table 2 indicate that the safety factor, F obtained from the finite element limit analysis method in this study are close to those from different limit equilibrium methods reported by Chen et al. (2014). Additionally, the safety factor obtained from this study falls within the range reported by Chen et al. (2014). The results presented in Table 2 also indicate that taking into account the vertical cracks within the slope results in lower safety factors as compared to that without cracks.

Figure 4 shows the variations of η with crack lengths, l _c . Results in Figure 4 indicate that: ① When the slope ratio, α = 1:1, η decreases linearly with normalized crack length, l _c /H. However, it is important to note that the F remains larger than 1.0 within the parameters analyzed in this study, indicating that the slope is in a stable state. ② Compared with the values of η obtained when α = 1:1, larger values of η are observed when α = 1:1.5 and 1:2. Additionally, unlike the linear trend observed for α = 1:1, the calculated values of η for α = 1:1.5 and 1:2 initially decrease slowly with l _c /H until reaching a critical value of l _c /H =25%, beyond which a significant decreasing trend is observed. ③ The use of reinforcements enhances the stability of slopes with cracks, especially when l _c /H > 25%. Furthermore, the effectiveness of the reinforcement is greater when the crack length is larger. For example, for slopes with l _c = 2 m and 4 m, the safety factor can be increased by approximately 15% and 18% respectively as compared to the those without reinforcement. ④ Overall, the results in Figure 4 indicate that crack length have a negative impact on slope stability, especially when the slope ratio α is large.

To further analyze the reasons for the variation of the safety factor, η with crack length, l _c shown in Figure 4, Figure 5 shows the shape of slip plane for the slope with an angle ratio α = 1:1.5 but different crack lengths, l _c . Numerous previous studies have indicated that slip planes typically originate from the slope crest and extend downward to the slope toe (Griffiths and Marquez, 2007; Li et al., 2009; Tschuchnigg et al., 2015), which is consistent with the results shown in Figure 5. Consequently, the cracks located to the right of the slope toe (Figure 3A) have minimal impact on the stability.

Results in Figure 5 indicate than: ① The slip plane is arc-shaped when l _c = 0 (Figure 5A), and the exit of the slip plane is located at the toe of the slope. ② The shape of the slip plane in case of l _c = 1 m (Figure 5B) is almost the same as that in Figure 5A (l _c = 0), and there is almost no intersection between the slip plane and the crack. As a result, the corresponding values of η in Figure 2 are close to each other. ③ When the crack length, l _c increases to 2 m and 3 m (Figures 5C, D), the slip planes intersect with several cracks. Due to the lower strength of the cracks compared to the adjacent soil, the values of η in Figure 4 decrease significantly at these points. ④ Compared with the unreinforced case (Figure 5C), the slip plane is located farther away from the slope surface when reinforcements are used (Figures 5E, F), especially the reinforcement length is larger, resulting in less steep slip plane. Therefore, the corresponding values of η are relatively larger in these cases.

The safety factors, F for slopes with different crack inclinations β, crack lengths l _c , and reinforcement lengths l _g are listed in Table 4. An inclination of β = 0° denotes a vertical crack, while β = 45° and −45° represent cracks rotating counterclockwise and clockwise, respectively, at an angle of 45° from the downward vertical line.

Results shown in Table 4 for slopes without reinforcement indicate that: ① When the crack length, l _c is 1 m, the inclination of the crack, β has almost no influence on the safety factor of slopes. This phenomenon may be attributed to the short length of the crack, leading to no crack intersecting with the slip plane. ② The minimum and maximum values of the safety factor, F are obtained when the crack inclination, β is 0°and −45°, respectively in the case of l _c = 2 m. ③ The influence of the crack inclination, β on the safety factor, F is more significant when l _c increases to 3 m as compared to those in cases of l _c = 1 m and 2 m.

Results shown in Table 4 also indicate that the utilization of reinforcement leads to an increase in values of the safety factor, F. Additionally, this effect is more significant for slopes with longer lengths and smaller crack inclinations.

The reasons for the variation of the safety factor, F with crack inclination, β shown in Table 4 can be explained by examining the shape of slip planes illustrated in Figure 6 for slopes with crack length l _c = 3 m. Results in Figure 6 indicate that: ① Cracks intersecting with the slip planes display noticeable deformation when β = −45° (Figure 6A). Although there are several inclined cracks in the slope, only relative small deformation is observed between two adjacent sliding soil masses, especially those in the middle of the slope. ② The failure mechanism of the slope is toppling collapse when the cracks are vertical (β = 0° shown in Figure 6B). Larger relative movement along the cracks is observed in this case, leading to a lower value of safety factor in Table 4. ③ The failure mechanism shown in Figure 6C is similar to that in Figure 6A, where only minimal relative deformation is observed between two adjacent sliding soil masses at the top of the slope. Besides, the slip planes overall align with the direction of the cracks, leading to a smaller value of F as compared to that in Figure 6A. ④ Compared with the failure mechanisms without reinforcement shown in Figures 6A–C, a global failure mechanism is observed in Figures 6D–F due to the reinforcement, which enhances slope integrity. The slip planes and deep cracks intersect with the reinforcement, leading to the soil masses moving together, and meanwhile, a greater value of safety factor is obtained. Furthermore, theses slip planes are located farther away from the slope surface as compared to those shown in Figures 6A–C.

Figure 7 shows the relationship between η and the crack friction angle, ϕ _c . Results in Figure 7 indicate that: ① Values of η for slopes without reinforcement are almost independent of the crack friction angle, ϕ _c in the case of slope with shallow cracks (i.e., ϕ _c = 1 m). This could be attributed the deep slip plane, leading to no intersection with cracks. ② As the crack length, ϕ _c increases to 2 m and 3 m, η remains constant regardless of the crack strength, ϕ _c when ϕ _c is small. However, these η values show a notable increase when ϕ _c exceeds 6°. Furthermore, the influence of ϕ _c on η becomes more significant with the crack length, l _c . ③ The values of η are greater for reinforced slopes, especially when the crack strength is lower. For example, for a slope with a crack length of 2 m and crack friction angle ranging from 0 to 10°, the use of reinforcements can increase the safety factor by approximately 13%.

Figure 8 shows the shape of slip planes for slopes with different soil friction angles when the crack length is 3 m. Results in Figure 8 indicate that the soil friction angle affects the shape of the slip planes. When the soil friction angle is relatively large, the slip planes are more close to the slope surface, resulting in a smaller failure zone and a larger safety factor.