Study on stability assessment of buckling failure and rainfall-induced instability mechanism for steeply inclined rock slope

2.1. Failure characteristics and mechanical model generalization of the slope

Field investigations indicate that the southern slope is a steeply inclined rock mass with a dip angle ranging from 72° to 80°. The failure mechanism exhibits typical buckling characteristics (Fig. 1): initial bending deformation occurs at the slope toe, leading to tensile fissures at the arc crown. As deformation progresses, these cracks propagate along the rock stratum strike. Ultimately, the slope fails via bending-shear instability governed by the tensile and shear strength of the fractures.

Fig. 1Diagram of slope failure modes

2.2. Differential equation for surface rock layers

A generalized mechanical model of the southern slope is established (Fig. 2) under the following assumptions:

(1) Compatibility exists in the bending deformation of steeply inclined rock layers, ensuring consistent interlayer deformation.

(2) The sliding segment’s action on the bending segment is simplified as a concentrated force.

(3) Post-deformation, the bending segment loses contact with the sliding surface, and deformations satisfy small-deformation conditions.

(4) The slope crest (at $x=0$) is a free end, while the slope toe (at $x=L$) is fixed.

Fig. 3 illustrates the micro-element force analysis of the surface rock layer.

Fig. 2Generalized mechanical model diagram of the slope

Fig. 3Micro-element force analysis diagram of slab-beam in surface rock layer of the slope

Buckling failure typically initiates from the surface layer. As shown in Fig. 4, key forces include: frost heaving force $P$; hydrostatic pressure $d_w$; slope mass self-weight $d_G$; seismic inertia force $d_S$; normal force $d_N$; interlayer friction $d_f$; cohesion $d_{\tau }$; and downslope load $F$. The downslope load is derived by integrating the contributions of external forces:

where: $d_G=\gamma bhdx$, $d_S={\alpha }_h\gamma bhdx$, $d_f=d_N/\mathrm{t}\mathrm{a}\mathrm{n}\phi$, $d_{\tau }=cbhdx$, $d_w={\gamma }_w{z}_{w}hdx$.

Here, $\alpha$ is the slope angle; $\gamma$ denotes the unit weight of rock mass (kN/m³); $b$ and $h$ represent the width and thickness of the rock mass (m), respectively; ${\alpha }_h$ is the horizontal seismic coefficient; $\phi$ and $c$ are the internal friction angle (°) and cohesion (kPa) of the weak layer; ${\gamma }_w$ is the unit weight of water (kN/m³); and $z_w$ is the groundwater table height (m).

The term $d_w$ incorporates pore water pressure effects, including uplift pressure on the sliding surface and the downslope component of seepage force. These are equivalently transformed into the downslope load $F$ to comprehensively capture groundwater impacts.

The total bending moment $M$ of the bending segment comprises two components: an additional moment $M_1$ from friction and cohesion, and a moment $M_2$ from hydrostatic pressure:

The deflection function $y\left(x\right)$ satisfies the boundary conditions ($y=0$ at $x=0$ and $x=L$):

where $n$ is the wave number, $A$ is the amplitude (m), and $L$ is the slope length (m).

2.3. Safety factor derivation based on energy balance

According to the energy balance principle, the deformation energy $∆U$ equals the work done by the downslope load $U_F$ and the freeze-thaw load $U_P$:

where:

where, $EI$ is the flexural rigidity (kN·m²), with $E$ as the elastic modulus (GPa) and $I$ as the cross-sectional moment of inertia (m⁴).

Substituting the deflection function and considering field observations (deformation pattern of one-and-a-half waves), the critical slope length $L_{cr}$ for $n=$ 2 is:

The adoption of $n=$ 2 (full sine wave) simplifies the energy balance approach. Although field observations suggest 1.5 waves, comparisons show that using $n=$ 1.5 yields a critical length 1.12 times that of $n=$ 2, with a safety factor increase of 8.3 %. Both trends align within acceptable engineering error margins (±10 %). The $n=$ 2 model is retained for its analytical tractability and accurate representation of stability mechanisms.

The safety factor $K$ is defined as the ratio of the critical slope length $L_{cr}$ to the actual slope length $L$:

3.1. Project overview

The southern slope of the Longyu Open-Pit Mine strikes approximately 171° and dips northward. The slope toe comprises a mine haul road, exhibiting a two-level stepped morphology with elevations ranging from 1365 m to 1405 m. The slope body is generally higher in the south and lower in the north, with the central portion elevated relative to both sides. The slope height is approximately 60 m. No distinct sliding surface has been identified within the bedrock portion of the slope; however, the bedrock rock mass is relatively fragmented, posing a potential risk of loosening and instability.

Fig. 4Real scene of southern slope failure. Photo by Jia Haipeng in May 2024 on the south slope of the Longyu Large Open-Pit Mine

3.2. Multi-factor analysis

To address the stability concerns of the target slope, seven sensitivity factors were selected based on prior research: slope angle ($\alpha$), elastic modulus ($E$), rock layer thickness ($h$), unit weight ($\gamma$), internal friction angle of the weak layer ($\phi$), cohesion of the weak layer ($c$), and groundwater level ($H_W$). Each factor was assigned three levels, and interactions between factors were disregarded. An orthogonal experimental design was implemented using the $L_{18}\left(3^7\right)$. orthogonal array. The specific factor levels are detailed in Table 1, and the experimental design is summarized in Table 2.

The range analysis method (R method) was employed to quantify the influence of each factor. The results are presented in Table 3.

In Table 3, $K_1$, $K_2$, and $K_3\mathrm{ }$denote the sums of test results for factor levels 1, 2, and 3, respectively; $K_1$, $K_2$, and $K_3$ represent the corresponding mean values; and $R$ is the range, defined as the difference between the maximum and minimum $k$ values. A larger $R$ indicates greater sensitivity of the factor.

The $\bar{K}$ values for each factor at different levels, derived from Tables 1-3, are illustrated in Fig. 5.

Analysis of Results: Fig. 5(a): Rock layer thickness ($h$) is a highly sensitive positive factor ($R=$ 0.38; increase rate: 39.8 %). As thickness increases from 10 m to 20 m, the safety factor exhibits an apparent “pseudo-negative fluctuation” (1.47 → 1.09), attributable to coupling effects such as adverse slope angles and elevated groundwater levels. After isolating these interferences, rock thickness emerges as the predominant factor influencing stability. Theoretical models indicate that slopes remain stable when $h>$ 10 m and $E>$ 32 GPa. Fig. 5(b): Cohesion ($c$) is a strongly sensitive positive factor ($R=$ 0.25). The safety factor peaks at 18 kPa ($k_2=$ 1.35), with a decline at 23 kPa due to interference from extreme conditions. Cohesion remains a key regulator of shear strength, balancing safety and economic efficiency. Fig. 5(c): Unit weight ($\gamma$) is a significantly sensitive negative factor ($R=$ 0.23). Despite an apparent positive trend (1.11 → 1.34), the increase is attributed to compensatory effects from the internal friction angle. Excluding this interference, each 5 kN/m³ increase in $\gamma$ reduces the safety factor by approximately 3 %. Fig. 5(d): Groundwater level ($H_W$) is a significantly sensitive negative factor ($R=$ 0.13). The safety factor reaches a critical risk threshold at $H_W=$ 8 m ($k_2=$ 1.17). Elevated water levels reduce effective stress via increased pore pressure, though compensatory effects from rock thickness may cause rebounds at higher levels. Fig. 5(e): Internal friction angle ($\phi$) is a weakly sensitive positive factor ($R=$ 0.12). The anti-sliding effect saturates beyond 25°, with minimal incremental gains in stability. Fig. 5(f): Slope angle ($\alpha$) is a weakly sensitive negative factor ($R=$ 0.11). The highest risk occurs at $h=$ 65° ($k_2=$ 1.16), influenced by synergistic effects with high $\gamma$ and $H_W$. Controlling $\alpha$ ≤ 65° is recommended to mitigate instability. Fig. 5(g): Elastic modulus ($E$) has negligible influence (R = 0.03). For $E$ > 35 GPa, the safety factor increases marginally (1.21 → 1.25), with a 1.6 % gain per 5 GPa increment, primarily affecting deformation rather than shear strength.

Fig. 5Variation of safety factor of various factors at different levels

a) Safety factor for different rock layer thicknesses

b) Safety factor for different cohesion values

c) Safety factor for different unit weights

d) Safety factor for different groundwater levels

e) Safety factor for different internal friction angles

f) Safety factor for different slope angles

g) Safety factor for different elastic moduli

The sensitivity ranking of the seven influencing factors, arranged in descending order of their $R$-values, is illustrated in Fig. 6. The corresponding significance classification of these factors is summarized in Table 4.

Fig. 6Sensitivity ranking of various factors

3.3. Univariate analysis

Based on Eq. (15), a univariate sensitivity analysis was conducted to evaluate the influence of four key parameters-rock layer thickness, elastic modulus, unit weight, and slope angle-on the stability of the target slope. The results are presented in Tables 5-8.

From Tables 5-8, the influence patterns of individual factors on the slope safety factor are derived and visualized in Fig. 7.

Fig. 7Sensitivity characteristics of slope safety factor to individual factors

a) Safety factor for different rock layer thicknesses

b) Safety factor for different elastic moduli

c) Safety factor for different rock unit weights

d) Safety factor for different slope angles

As depicted in Fig. 7(a), under constant values of other parameters, the safety factor increases monotonically with rock layer thickness. When the thickness increases from 10 m to 20 m, the safety factor rises from 1.03 (approaching the critical stability threshold of 1.0) to 1.44, representing a total increase of 39.8 %. This indicates that rock layer thickness is a key positive controlling factor: greater thickness enhances both the deformation resistance of the slope structure and the anti-sliding capacity of the rock mass, thereby significantly improving stability.

Fig. 7(b) shows that the safety factor increases continuously with elastic modulus. As the modulus rises from 35 GPa to 45 GPa, the safety factor grows from 1.23 to 1.40, an overall increase of 13.8 %. This trend suggests that a higher elastic modulus improves rock mass stiffness and deformation resistance, thereby enhancing stability. However, the magnitude of this effect is smaller than that of rock layer thickness, classifying elastic modulus as a weak positive factor.

In contrast, Fig. 7(c) reveals a monotonically decreasing trend in the safety factor as unit weight increases. When unit weight increases from 20 kN/m³ to 30 kN/m³, the safety factor declines from 1.34 to 1.26, a total reduction of 5.97%. This behavior stems from the increase in downslope gravitational force, which promotes sliding and reduces stability. Nevertheless, the decline remains moderate due to the constant shear strength parameters.

Similarly, Fig. 7(d) illustrates that the safety factor decreases with increasing slope angle. As the slope angle rises from 60° to 70°, the safety factor drops from 1.32 to 1.26, a total decrease of 4.55%. Steeper slopes exhibit reduced resistance to sliding, leading to lower stability. The effect is comparable in magnitude to that of unit weight.

In summary, rock layer thickness and elastic modulus act as positive factors, whereas unit weight and slope angle are negative factors. Regarding the magnitude of influence, the order is as follows:

Rock layer thickness (39.8 % increase) > Elastic modulus (13.8 % increase) > Unit weight (5.97 % decrease) ≈ Slope angle (4.55 % decrease).

Notably, variations in rock layer thickness and slope angle within specific intervals lead to nonlinear changes in safety factor, as indicated by the segmented trends in Fig. 6. This underscores the need for staged control strategies in slope design. In practice, increasing rock layer thickness should be prioritized-for example, via shotcrete and bolt support systems, which cost approximately CNY 800/m² and reduce expenses by 30 % compared to traditional anti-slide piles. Such methods, combined with top-down construction and temporary anchoring, have been successfully applied in the northern slope of the mining area, achieving a 95 % qualification rate. Modifying the elastic modulus offers secondary benefits, while adjustments to unit weight and slope angle should be weighed against economic feasibility, as excessive optimization yields diminishing returns.

4.1. Parameter selection and model establishment

Based on the failure characteristics observed at the southern slope of Longyu Mine, as validated through field geological surveys and laboratory tests, slope instability is primarily governed by shear failure. The failure mechanism along the potential sliding surface (i.e., the mudstone layer) conforms to the Mohr-Coulomb strength criterion. Accordingly, the Mohr-Coulomb constitutive model was adopted for numerical simulations, and the strength reduction method was applied to compute the factor of safety. The southern slope model was constructed using FLAC3D software, as illustrated in Fig. 8.

Fig. 8Slope model

During the initial rainfall phase (0-6 h), unsaturated soil parameters were employed, with shear strength maintained at 80 %-90 % of the dry-state value. After 6 h, when the rock mass becomes fully saturated, saturated parameters were applied. This staged parameterization strategy mitigates excessive conservatism in early simulation stages, thereby improving alignment with actual rainfall infiltration behavior. In terms of boundary conditions, normal constraints were applied to the side surfaces of the model, while full three-dimensional constraints were imposed at the bottom. For slope stability calculation, the $X$-axis was defined as the slope surface direction, and the $Y$-axis as the strike direction of the slope.

4.2. Rainfall scenario design and infiltration calculation

Based on the Encyclopedia of Chinese Meteorological Disasters and local meteorological data from Longyu Mine, three rainfall scenarios were designed (Table 10). Rainfall intensity follows a triangular distribution pattern-intensity increases initially, stabilizes midway, and decreases thereafter. The infiltration rate was calculated using Darcy’s law:

where $q$ is the rainfall infiltration rate (m/s), defined as the volume of water infiltrating per unit surface area per unit time; $k$ denotes the saturated hydraulic conductivity (m/s), an intrinsic parameter indicating soil permeability; and $i$ represents rainfall intensity adjusted by a runoff coefficient of 0.3. This runoff coefficient was determined from field infiltration tests at Longyu Mine (average value: 0.28-0.32 from three parallel tests) and is consistent with the range (0.2-0.4) recommended for similar lithologies in the Code for Geotechnical Investigation and Survey (GB 50021-2001). It should be noted that this value is site-specific and not universally applicable.

4.3. Analysis of simulation results

Previous studies [22], [23] have indicated that for slopes with pronounced step-like geometries, the plane strain model and the 3D model yield comparable stability results, owing to the constraining effect of the steps. Accordingly, a plane strain model was derived from the 3D model in Fig. 8, and numerical analyses were performed under four rainfall conditions. The focus was placed on the displacement distribution and shear strain increment characteristics under each scenario.

(1) Variation of Displacement Field.

As depicted in Fig. 9, slope displacement under natural conditions exhibits localized concentration, with a maximum displacement of 2 m at the 1390 m platform, indicating a region of marginal initial mechanical equilibrium.

Fig. 9Displacement contour map under different states

a) Initial state

b) Light to moderate rain condition

c) Heavy rain condition

d) Torrential rain condition

With increasing rainfall intensity, the displacement field demonstrates a dual trend of spatial expansionand magnitude surge. Under light to moderate rainfall, displacement growth is moderate. In contrast, under heavy and torrential rainfall, displacement increases markedly. Under torrential conditions, the maximum displacement reaches 5.5 m – a 175 % increase compared to the natural state. The zone of significant displacement expands from the single 1390 m platform to a continuous section between 1360 m and 1390 m, forming a potential slip zone. Moreover, the affected area extends toward the slope toe, suggesting that rainfall-induced pore water pressure disrupts the deep mechanical equilibrium, leading to a trend of overall instability.

(2) Distribution of Shear Strain Increment.

As illustrated in Fig. 10, under natural conditions, the high-value zones of shear strain increment are concentrated within the elevation range of 1330-1390 m, exhibiting a discrete distribution pattern. Following rainfall, these zones form an arc-shaped potential sliding surface through vertical extension and regional connectivity, indicating the weakening effect of rainfall on the slope’s shear strength. Under torrential rainfall conditions, the maximum shear strain increment increases from 5×10^-2 to 9×10^-2, representing an 80 % increase. This suggests that the risk of slope shear failure increases linearly with rainfall intensity. The evolution of the sliding surface demonstrates a clear downward migration trend: the shear outlet is located at the 1310 m platform under heavy rainfall and extends further to 1290 m under torrential rainfall. This reflects the continuous development of the shear failure surface toward deeper parts of the slope as rainfall infiltration depth increases. The formation and connectivity of this arc-shaped sliding surface serve as a key indicator of the transition from local instability to a global landslide.

Fig. 10Shear strain increment contour map under different states

a) Initial state

b) Light to moderate rain condition

c) Heavy rain condition

d) Torrential rain condition

Simulation results indicate that rainfall exacerbates slope instability through two primary mechanisms: (i) rainwater infiltration increases the saturation of the rock-soil mass, reduces effective stress, and attenuates shear strength; (ii) pore water pressure accumulation induces seepage forces, which drive shear displacement along the potential sliding surface, forming a vicious cycle of “stress degradation and displacement-driven failure.”

Furthermore, the simulations clarify the rainfall intensity threshold for slope stability: the slope remains relatively stable under light to moderate rainfall; heavy rainfall represents a critical point for stability transition; and under torrential rainfall, the slope approaches its instability limit. This threshold offers a quantitative basis for regional rainfall early-warning classification. Therefore, for this project, rainfall intensity exhibits a significant negative correlation with slope stability. Under torrential rainfall, the slope experiences a surge in displacement and connectivity of shear strain increment. The elevation range of 1360-1390 m is identified as the core risk area, while 1290-1310 m constitutes the key shear outlet zone, necessitating focused landslide prevention and control measures.

4.4. Field monitoring validation

This study acquired 3D displacement data of the slope surface via GNSS monitoring points installed on stable deep bedrock (Fig. 11), providing direct validation for the numerical simulation results. Monitoring data from April 20 to July 21, 2025, indicated that cumulative displacement was controlled within ±10 mm. The corresponding displacement time-history is shown in Fig. 12.

Fig. 11GNSS Monitoring Points. Photo by Jia Haipeng in July 2024 on the south slope of the Longyu Large Open-Pit Mine

Fig. 12Displacement time-history curve

After the rainy season, localized small-scale landslides were observed on the slope surface. The maximum in-situ displacement measured was 5.2 m (Fig. 13), which aligns well with the numerically predicted value of 5.5 m. The error falls within an acceptable range for engineering applications, thereby verifying the reliability of the simulation.

Fig. 13Localized Landslide Map. Photo by Jia Haipeng in September 2025 on the south slope of the Longyu Large Open-Pit Mine

Based on the above, the current slope deformation is characterized by overall deep-seated subsidence accompanied by slight slope-direction creep, with small-scale landslides in the surface rock mass. This deformation pattern is typical of high-steep slopes in deep open-pit mines. The agreement between GNSS monitoring data and numerical simulations confirms the model’s predictive capability under heavy rainfall conditions.