A new analytical method for stability analysis of rock blocks with basal erosion in sub-horizontal strata by considering the eccentricity effect

The basal cavity of a rock block formed due to differential weathering is an important predisposing factor for rockfall in hard–soft interbedded rocks, which induces an eccentricity situation at the base of the rock block. Rock block falling due to the non-uniform distribution with the failure modes of toppling or sliding is defined as biased rockfall in this study. Taking into account the non-uniform stress distribution due to the eccentricity effect, a new analytical method is proposed for three-dimensional stability analysis of biased rockfall. The development of non-uniform stress distribution stress calculated by this analytical method was verified by numerical simulation. The biased rockfall progresses from partial damage of the soft underlying layer, caused by non-uniform distributed stress, to toppling and sliding of overhanging hard rock block due to overall unbalanced force. Therefore, a set of factors of safety (Fos's) against partial damage (compressive and tensile damage of the soft underlying layer) and overall failure (toppling and sliding of the hard rock block) are used to determine the rockfall susceptibility level. The analytical method is applied and validated using biased rockfalls on the northeastern edge of the Sichuan Basin in southwest China, where a significant number of rockfalls consisting of overhanging thick sandstone and underlying mudstone occur. The evolution process of biased rockfalls is divided into four stages, initial state, basal cavity formation, partially unstable and failure. The proposed method is validated by calculating Fos's of the typical unstable rock blocks in the study area. As the cavity continues to grow, the continuous retreat of mudstone causes stress redistribution between the hard and soft rock layers. This results in damage to the underlying soft rock layer due to the development of the non-uniform distribution, ultimately leading to the failure of the hard rock block. The critical retreat ratio is determined to be 0.33, which is used to classify the low and moderate rockfall susceptibility in the eastern Sichuan Basin. The proposed analytical method provides insights into the evolution of biased rockfall and a means for early identification and susceptibility assessment of rockfall.

1 Introduction

Rockfall is defined as the detachment of a rock block from a steep slope along a surface, on which little or no shear displacement takes place (Cruden and Varnes, 1996). Rockfalls frequently occur in mountainous ranges, on cut slopes and on coastal cliffs, and they may cause significant facility damage and casualties in residential areas and transport corridors (Chau et al., 2003; Volkwein et al., 2011; Corominas et al., 2018). Stability analysis of rock blocks is crucial for risk management and early warning of rockfall (Kromer et al., 2017).

Figure 1 Potentially unstable blocks and basal cavities caused by differential weathering.

Rockfalls are prone to occur in soft–hard rock formations, and the non-uniform stress distribution caused by differential weathering of rock formations is the main reason for the failure of rockfall. In the eastern Sichuan Basin, southwest China, rockfall is widespread and poses a high risk (Chen et al., 2008; Chen and Tang, 2010; Zhang et al., 2016; Zhou et al., 2017, 2018). The rockfall in this area is attributed to the tectonic setting of Jura-type folds and the stratum sequence, which is characterized by the interbedding of hard and soft layers. An alternation of thick sandstone and thin mudstone layers is formed in the wide and gentle-angle synclines (Zhang et al., 2015; Wu et al., 2018). Weathering is known to be one of the main predisposing factors in rockfall (Jaboyedoff et al., 2021; Zhan et al., 2022). The cliff comprised of hard sandstone is the source of rockfall, and the underlying mudstone is more susceptible to weathering. Along with the retreat of basal cavities in the mudstone layer, the gravity centre of the overhanging sandstone block moves outward relative to the mudstone. In this case, the stress distribution in the contact surface of sandstone and mudstone is non-uniform. The mudstone on the outer side bears higher compressive stress than that on the inner side. This phenomenon can be defined as an eccentricity effect, which leads to mudstone damage and failure of the overhanging sandstone by toppling or sliding. This type of rockfall is defined as biased rockfall in this study (Fig. 1). Similar rockfall patterns have been widely reported in other regions, such as Joss Bay in England (Hutchinson, 1972), Okinawa Island in Japan (Kogure et al., 2006) and the Colorado Plateau of the southwestern United States (Ward et al., 2011). Retreat of the basal cavity is a main cause for the failure of the overhanging block. Therefore, it is necessary to establish an analytical method, considering the development of the basal cavity, to analyse the stress distribution and stability of rock blocks, which is fundamental to the susceptibility assessment and risk control of biased rockfall.

Rockfall stability analysis methods include statistical analysis (Frattini et al., 2008; Santi et al., 2009), empirical rating systems (Pierson et al., 1990; Ferrari et al., 2016) and mechanical analysis (Lin and Fairhurst, 1988; Jaboyedoff et al., 2004; Derron et al., 2005; Matasci et al., 2018). The statistical analysis and empirical rating systems are suitable for rockfall hazard assessment at a regional scale. The accuracy of statistical analysis depends on the completeness of rockfall inventories (Chau et al., 2003; Guzzetti et al., 2003; D'Amato et al., 2016). However, its application to rockfall hazards is limited due to the lack of complete inventory data (Budetta and Nappi, 2013; Malamud et al., 2004). Empirical and semi-empirical rating systems are used where site-specific rockfall inventories are either unavailable or unreliable. Therefore, rockfall susceptibility can be assessed by heuristic ranking of selected predisposing factors (Frattini et al., 2008; Budetta, 2004). Mechanical analysis based on static equilibrium theory is the main method to analyse the stability of site-specific rockfall using the factor of safety (Fos). Ashby (1971) conducted stability analysis with a parallelepiped block resting on an inclined plane (Fig. 2a), and the solution was subsequently modified by Bray and Goodman (1981) and Sagaseta (1986). Kogure et al. (2006) utilized a cantilever beam model to determine the critical state of limestone cliffs. Frayssines and Hantz (2009) proposed the limit equilibrium method (LEM) to predict block stability against sliding and toppling in steep limestone cliffs (Fig. 2c). Chen and Tang (2010) established a stability analysis method of three types of unstable rocks in the Three Gorges Reservoir area with the LEM. Alejano et al. (2015) studied the influence of rounding of block corners on the block stability. Zhang et al. (2016) defined Fos based on fracture mechanics and studied the progressive failure process by analysing crack propagation. Alejano et al. (2010) and Pérez-Rey et al. (2021) deduced a formula for the Fos of blocks with more complex geometry.

Figure 2 Traditional force analysis diagrams of the rock block. Panels (a) and (b) are stability analysis diagrams of rock blocks under dynamic conditions, resting on an inclined plane with a dip angle of α. The rock block is generalized as a cuboid with dimensions b×h and weight W (as modified from Ashby, 1971; Bray and Goodman, 1981; and Sagaseta, 1986). (c) Force description of the toppling model proposed by Frayssines and Hantz (2009). In the above assumptions, N, T and W are regarded as forces applied at a point.

The supporting force on the contact surface is assumed to be applied at a point in the current LEMs (i.e. N in Fig. 2b and c). However, the supporting force is actually a distributed force. The external erosion generates an eccentricity situation on the overhanging rock block and results in a non-uniform distribution of the supporting force on the contact surface, which is not considered in the traditional LEM. The presence of non-uniform stress distribution plays a critical role in inducing localized damage within a rock mass. Damage is frequently considered an indicator or a threshold for the onset of accelerated failure in rock masses (Zhang et al., 2016). Therefore, it is imperative to consider the non-uniform stress distribution for the rockfall stability analysis. Furthermore, most studies simplified the three-dimensional geometry of the rock block using one cross-section to represent the critical features of the slope structure. Nevertheless, for natural blocks with basal cavities, the cavities usually present different depths along different directions (Pérez-Rey et al., 2021). Therefore, a three-dimensional model is necessary to accurately calculate the stability. In addition, when a block has multiple free faces and a complex structure, its potential failure is dominated by different modes, including rock mass damage and overall block failure. Therefore, the probable failure modes should be determined prior to the calculation of Fos.

Based on rockfall investigation in the eastern Sichuan Basin, China, the main objective of this study was to propose a new three-dimensional method for the determination of failure modes and the Fos of biased rockfall, considering the non-uniform force distribution on the contact surfaces. Compared with the traditional LEM, this study takes into account the partial damage of the underlying soft rock and the overall instability of the overhanging hard rock blocks and can evaluate the stability of biased rockfall more comprehensively. Fos's of the typical unstable rock blocks in the study area are calculated to validate the proposed method. In addition, the critical mudstone retreat ratio in this area is analysed. This study is an extension of the basic LEM for rockfall, which can promote the accuracy of rockfall stability analysis and facilitate rockfall prevention and risk mitigation.

2 Study area

2.1 Geological setting

The study area is located on the northeastern edge of the Sichuan Basin, China (Fig. 3a). Continuous erosion processes generate moderate–low mountain and valley landforms (Yu et al., 2021). The tectonic structure of this area is characterized by a series of ENE anticlines and synclines (Fig. 3b and c). In the anticline area, the rock layers dip relatively steeply, where translational rockslides are the main mode of slope failure. The syncline area is dominated by gently dipping strata and is prone to rockfall (Zhou et al., 2018). The study area is located in the core of the Matouchang syncline, where the rock layers are sub-horizontal (Fig. 3d and e). In this valley, due to the long-standing fluvial incision, the relative relief is approximately 500 m and the valley flanks are extremely steep (Fig. 3e). In addition, the toes of the hillslopes are reshaped because of the construction of the G318 national road, which is the main traffic line and is always threatened by rockfalls from steep rock slopes (shown in Fig. 3d and Table 1).

Figure 3 (a) Location of the study area in China; (b) geological map of the study area; (c) tectonic sketch profile of A–A^′, whose location is shown in (b); (d) rockfall-prone segment and key investigation areas. The red dots are the positions of historical rockfall events, corresponding to the numbers in Table 1. (e) Geological cross-section of the hillslope in the Jitougou section of G318 national road, which is marked by a red rectangle in (d). Publisher's remark: please note that the above figure contains disputed territories.

Table 1 Historical rockfall events along the G318 national road in the study area.

^* Note: the exact time is unknown.

Figure 4 Characteristics of biased rockfalls in the study area. (a) Multiple layers of rockfall sources, which is consist of thick sandstone. (b) Two sets of sub-vertical joints (F1 and F2) recognized by the UAV photos. (c) Large basal cavity developed in the underlying mudstone. (d) Dense fractures on the mudstone surface generated by weathering and compression. (e) Vertical tension crack in the rear of the block, through which precipitation can infiltrate.

2.2 Rockfall characteristics

The slopes in the study area consist of a sub-horizontally interbedded sandstone and mudstone layer. Therefore, there are multiple layers of potentially unstable rock blocks in the hillslopes (Fig. 4a). The thick sandstone has two sets of sub-vertical joints (Fig. 5), which cut the rock mass into blocks as the potential rockfall source (Fig. 4b). Basal cavities have formed in the underlying mudstone layer (Fig. 4c and d). Joints and bedding planes (BPs) constitute the detachment surfaces between the blocks and steep slope (Fig. 4e). The eccentricity effect produced by the mudstone external erosion plays an important role in the evolution process of rockfall. When the basal mudstone cannot provide adequate supporting force, the blocks detach from the steep slope, and biased rockfall occurs. Sliding and toppling are two possible failure modes of biased rockfall.

Figure 5 Stereonet produced using compass-clinometer survey data, showing the densities and orientations of five clusters. The data were collected in the rockfall-prone area shown in Fig. 3d.

According to the historical rockfall events in this area, precipitation is considered a triggering factor of rock instability. The precipitation mainly infiltrates along the sub-vertical joints or cracks of the sandstone (Fig. 4e). However, the drainage of fissure water is hysteretic due to the obstruction of basal mudstone. Therefore, transient steady flow exists in vertical cracks during heavy rainfall, and the hydrostatic pressure triggers the detachment of rock blocks. Thus, typical scenarios (such as rainfall intensity and earthquake) need to be considered in the stability analysis model.

Figure 6 The unstable blocks were labelled W02, W08, W18, W04 and W21, which are detached by the dominating discontinuities in Fig. 5. Basal cavities can be identified under the bedding planes of sandstone.

Figure 7 The evolution process of rock blocks from a stable state to failure.

3 Calculation method

3.1 Geological models and assumptions

A detailed geological investigation of unstable rock blocks was carried out in the study area (Fig. 6). The geological model of the rock block is mainly composed of the overhanging sandstone and the underlying mudstone. The sandstone block is assumed to be a rigid body, which is divided by two sets of orthogonal vertical smooth joints without friction resistance. According to the relatively persistent sub-vertical fractures observed in the field, the vertical joints are assumed to be fully persistent in the geological model. The sandstone block is assumed to be a complete body without persistent discontinuity, and it will not disintegrate before it falls. Due to the basal cavity in mudstone, the contact surface between sandstone and mudstone exhibits an eccentricity situation where non-uniform stresses are distributed at different positions. Mudstone is mainly loaded by compressive stress and tensile stress. When the compressive stress of mudstone exceeds its strength on the outer side, some initial damage appears. The effective contact surface between mudstone and sandstone is reduced, which aggravates the non-uniform distribution of stress. In this way, the ability of mudstone to resist the sliding and toppling of overhanging sandstone is reduced. In the field, compression deformation of mudstone can be observed, which usually manifests as micro-fractures and cleavages (Fig. 4d). The deformation is very slight and slow in the short term. In addition, the LEM is essentially a force and stress approach that does not take into account the deformation. Therefore, in this study, it is assumed that the mudstone is not subjected to deformation. The rock block remains in the state of static equilibrium prior to the final overall failure. Figure 7 displays the four evolutionary stages of biased rockfall. In the initial stage, the base cavity has not yet formed, and the normal force acting on the contact surface is uniform in different positions. The eccentricity effect leads to a non-uniform supporting force as the basal cavity grows, and partial damage gradually develops when the non-uniform stress exceeds the compressive or tensile strength of the mudstone. Under the triggering effects of rainfall or earthquakes, the rock blocks are separated by sliding or toppling.

Figure 8 represents the mechanical model of the force equilibrium analysis of a rock block with two or three free faces. The rock block (the overhanging sandstone) is generalized as a parallelepiped block. The underlying mudstone is impermeable, so rainfall can fill the joints and transmit horizontal hydrostatic pressure. The shear strength of the underlying mudstone is assumed to obey the Mohr–Coulomb criterion. Rainfall and earthquakes decrease Fos by generating hydrostatic pressure H in the vertical crack and horizontal seismic force E on the block.

Figure 8 Diagram of the force equilibrium analysis of the rock block model. Panels (a) and (b) represent the case of unstable rock blocks with two or three free vertical surfaces, respectively.

A Cartesian coordinate system is established in three-dimensional space for the force analysis. The origin O is located at the centre of the contact surface between sandstone and mudstone. For the case with two free surfaces, the orientation of the free surfaces is set to be the positive direction of the x axis and y axis. For the case with three free surfaces, the negative direction of the x axis is also a free surface. Joint J2 is perpendicular to the x axis, and joint J1 is perpendicular to the y axis.

3.2 Calculation processes

3.2.1 Stress distribution at the block base

The following formulas are used to calculate the apparent dip of α (θ₁ and θ₂):

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }{\mathit{\theta }}_{\mathrm{1}}=\arctan \left(\tan \mathit{\alpha }\cdot \cos {\mathit{\omega }}_{\mathrm{1}}\right),\text{(2)}& {\displaystyle }& {\displaystyle }{\mathit{\theta }}_{\mathrm{2}}=\arctan \left(\tan \mathit{\alpha }\cdot \cos {\mathit{\omega }}_{\mathrm{2}}\right),\end{array}$$

where ω₁ and ω₂ are the angles between the trend of the contact surface and the x direction or y direction, respectively.

As shown in Fig. 8b, with respect to the x axis, gravity, seismic forces and hydrostatic pressure create a non-symmetrical stress distribution on the foundation. The bending moment of gravity with respect to the x axis ($M_{\mathrm{b}{W}_x}$) is

$$\begin{array}{}\text{(3)}& M_{\mathrm{b}{W}_x}=W\cdot {\displaystyle \frac{d_{\mathrm{1}}-d_{\mathrm{3}}}{\mathrm{2}}}\cos {\mathit{\theta }}_{\mathrm{1}}.\end{array}$$

Assuming that the height of the water in the fracture is h_w, the hydrostatic pressure along the x direction (H_x) and its bending moment ($M_{\mathrm{b}{H}_x}$) are, respectively, expressed as

$$\begin{array}{}\text{(4)}& {\displaystyle }{H}_x={\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{w}}{h}_{\mathrm{w}}^{\mathrm{2}}}{\mathrm{2}}}\left(b-d_{\mathrm{2}}\right),\text{(5)}& {\displaystyle }\begin{array}{rl}{M}_{\mathrm{b}{H}_x}=& \underset{-\frac{b-d_{\mathrm{2}}}{\mathrm{2}}}{\overset{\frac{b-d_{\mathrm{2}}}{\mathrm{2}}}{\int }}\underset{\mathrm{0}}{\overset{h_{\mathrm{w}}\cos {\mathit{\theta }}_{\mathrm{1}}}{\int }}{\mathit{\gamma }}_{\mathrm{w}}\left(h_{\mathrm{w}}-{\displaystyle \frac{z}{\cos {\mathit{\theta }}_{\mathrm{1}}}}\right)\\ & \left({\displaystyle \frac{z}{\cos {\mathit{\theta }}_{\mathrm{1}}}}+{\displaystyle \frac{a-d_{\mathrm{1}}-d_{\mathrm{3}}}{\mathrm{2}}}\cdot \sin {\mathit{\theta }}_{\mathrm{1}}\right)\mathrm{d}z\mathrm{d}y.\end{array}\end{array}$$

The horizontal seismic force along the x direction (E_x) and its bending moment ($M_{\mathrm{b}{E}_x}$) are, respectively, expressed as

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }{E}_x=k_{\mathrm{e}}W\text{(7)}& {\displaystyle }& {\displaystyle }{M}_{\mathrm{b}{E}_x}=E_x\left({\displaystyle \frac{h}{\mathrm{2}}}-{\displaystyle \frac{d_{\mathrm{1}}-d_{\mathrm{3}}}{\mathrm{2}}}\sin {\mathit{\theta }}_{\mathrm{1}}\right).\end{array}$$

The total applied vertical load (N_z) and the total bending moments along the x direction ($M_{{\mathrm{b}}_x}$) can be derived as

$$\begin{array}{}\text{(8)}& {\displaystyle }\begin{array}{rl}{N}_z& =W\cos \mathit{\alpha }-\left(H_x\cdot k_{\mathrm{1}}\cdot k_{\mathrm{3}}+E_x\cdot k_{\mathrm{2}}\right)\sin {\mathit{\theta }}_{\mathrm{1}}\\ & -\left(H_y\cdot k_{\mathrm{1}}+E_y\cdot k_{\mathrm{2}}\right)\sin {\mathit{\theta }}_{\mathrm{1}},\end{array}\text{(9)}& {\displaystyle }{M}_{{\mathrm{b}}_x}=M_{\mathrm{b}{W}_x}+M_{\mathrm{b}{H}_x}\cdot k_{\mathrm{1}}\cdot k_{\mathrm{3}}+M_{\mathrm{b}{E}_x}\cdot k_{\mathrm{2}},\end{array}$$

where k₁, k₂ and k₃ are the coefficients set to make Eqs. (8) and (9) compatible with different calculation scenarios. Therefore, Eqs. (8) and (9) and the following formulas can be expressed in a unified form. In the natural scenario, k₁ and k₂ are both equal to 0. In the rainfall scenario, k₁=1. In the earthquake scenario, k₂=1. For the case of two free faces, k₃=1. For the case of three free surfaces, k₃=0.

Based on bending theory (Adrian, 2010), the eccentricity distance along the x direction (e_x) can be expressed as

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& e_x={\displaystyle \frac{M_{{\mathrm{b}}_x}}{N_z}}=\\ & {\displaystyle \frac{M_{\mathrm{b}{W}_x}+M_{\mathrm{b}{H}_x}\cdot k_{\mathrm{1}}\cdot k_{\mathrm{3}}+M_{\mathrm{b}{E}_x}\cdot k_{\mathrm{2}}}{W\cos \mathit{\alpha }-\left(H_x\cdot k_{\mathrm{1}}\cdot k_{\mathrm{3}}+E_x\cdot k_{\mathrm{2}}\right)\sin {\mathit{\theta }}_{\mathrm{1}}-\left(H_y\cdot k_{\mathrm{1}}+E_y\cdot k_{\mathrm{2}}\right)\sin {\mathit{\theta }}_{\mathrm{1}}}}.\end{array}\end{array}$$

The same method can be used to obtain e_y:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& e_y={\displaystyle \frac{M_{{\mathrm{b}}_y}}{N_z}}=\\ & {\displaystyle \frac{M_{\mathrm{b}{W}_y}+M_{\mathrm{b}{H}_y}\cdot k_{\mathrm{1}}+M_{\mathrm{b}{E}_y}\cdot k_{\mathrm{2}}}{W\cos \mathit{\alpha }-\left(H_x\cdot k_{\mathrm{1}}\cdot k_{\mathrm{3}}+E_x\cdot k_{\mathrm{2}}\right)\sin {\mathit{\theta }}_{\mathrm{1}}-\left(H_y\cdot k_{\mathrm{1}}+E_y\cdot k_{\mathrm{2}}\right)\sin {\mathit{\theta }}_{\mathrm{1}}}}.\end{array}\end{array}$$

According to the stress distribution of a rectangular-shaped foundation (Adrian, 2010), the stress in the (x, y) coordinates, p(x,y), is

$$\begin{array}{}\text{(12)}& p(x,y)={\displaystyle \frac{N}{A}}+{\displaystyle \frac{N{e}_x}{I_y}}x+{\displaystyle \frac{N{e}_y}{I_x}}y,\end{array}$$

with the formulas

$$\begin{array}{}\text{(13)}& {\displaystyle }& {\displaystyle }{I}_x={\displaystyle \frac{\left(a-d_{\mathrm{1}}\right){\left(b-d_{\mathrm{2}}\right)}^{\mathrm{3}}}{\mathrm{12}}},\text{(14)}& {\displaystyle }& {\displaystyle }{I}_y={\displaystyle \frac{\left(b-d_{\mathrm{2}}\right){\left(a-d_{\mathrm{1}}\right)}^{\mathrm{3}}}{\mathrm{12}}},\text{(15)}& {\displaystyle }& {\displaystyle }A=\left(a-d_{\mathrm{1}}-d_{\mathrm{3}}\right)\left(b-d_{\mathrm{2}}\right).\end{array}$$

By substituting Eqs. (13)–(15) into Eq. (12), p(x,y) can be derived as

$$\begin{array}{}\text{(16)}& \begin{array}{rl}p(x,y)=& {\displaystyle \frac{N}{A}}\left[\mathrm{1}+{\displaystyle \frac{\mathrm{12}{e}_x}{{\left(a-d_{\mathrm{1}}-d_{\mathrm{3}}\right)}^{\mathrm{2}}}}x+{\displaystyle \frac{\mathrm{12}{e}_y}{{\left(b-d_{\mathrm{2}}\right)}^{\mathrm{2}}}}y\right],\\ & x\in \left[-{\displaystyle \frac{a-d_{\mathrm{1}}-d_{\mathrm{3}}}{\mathrm{2}}},{\displaystyle \frac{a-d_{\mathrm{1}}-d_{\mathrm{3}}}{\mathrm{2}}}\right],\\ & y\in \left[-{\displaystyle \frac{b-d_{\mathrm{2}}}{\mathrm{2}}},{\displaystyle \frac{b-d_{\mathrm{2}}}{\mathrm{2}}}\right].\end{array}\end{array}$$

p_max and p_min can be derived from Eq. (16) as

$$\begin{array}{}\text{(17)}& {\displaystyle }& {\displaystyle }{p}_{\max }=p\left({\displaystyle \frac{a-d_{\mathrm{1}}-d_{\mathrm{3}}}{\mathrm{2}}},{\displaystyle \frac{b-d_{\mathrm{2}}}{\mathrm{2}}}\right),\text{(18)}& {\displaystyle }& {\displaystyle }{p}_{\min }=p\left(-{\displaystyle \frac{a-d_{\mathrm{1}}-d_{\mathrm{3}}}{\mathrm{2}}},-{\displaystyle \frac{b-d_{\mathrm{2}}}{\mathrm{2}}}\right).\end{array}$$

The mudstone foundation has both compressive strength and tensile strength, so the value of p(x,y) is modified to obtain the two piecewise functions

$$\begin{array}{}\text{(19)}& {\displaystyle }& {\displaystyle }{p}_{\mathrm{p}}(x,y)=\left\{\begin{array}{ll}{\mathit{\sigma }}_{\mathrm{cmax}},& p(x,y)\ge {\mathit{\sigma }}_{\mathrm{cmax}},\\ p(x,y),& \mathrm{0}<p(x,y)\le {\mathit{\sigma }}_{\mathrm{cmax}},\\ \mathrm{0},& p(x,y)<\mathrm{0},\end{array}\right.\text{(20)}& {\displaystyle }& {\displaystyle }{p}_{\mathrm{n}}(x,y)=\left\{\begin{array}{ll}\mathrm{0},& p(x,y)<-{\mathit{\sigma }}_{\mathrm{tmax}},\\ p(x,y),& -{\mathit{\sigma }}_{\mathrm{tmax}}\le p(x,y)<\mathrm{0},\\ \mathrm{0},& p(x,y)\ge \mathrm{0}.\end{array}\right.\end{array}$$

Here, p_p(x,y) provides support force for the overhanging sandstone, and p_n(x,y) provides tension force.

3.2.2 Calculation of factors of safety

According to the Mohr–Coulomb criterion, the ultimate shear strength τ_max is

$$\begin{array}{}\text{(21)}& {\mathit{\tau }}_{\max }=\underset{-\frac{a-d_{\mathrm{1}}-d_{\mathrm{3}}}{\mathrm{2}}}{\overset{\frac{a-d_{\mathrm{1}}-d_{\mathrm{3}}}{\mathrm{2}}}{\int }}\underset{-\frac{b-d_{\mathrm{2}}}{\mathrm{2}}}{\overset{\frac{b-d_{\mathrm{2}}}{\mathrm{2}}}{\int }}\left[p_{\mathrm{p}}(x,y)\tan \mathit{\phi }+c\right]\mathrm{d}y\mathrm{d}x.\end{array}$$

Therefore, Fos against sliding, Fos_sl, can be defined as

$$\begin{array}{}\text{(22)}& \begin{array}{rl}& {\mathrm{Fos}}_{\mathrm{sl}}={\displaystyle \frac{S_{\mathrm{stabilizing}}}{S_{\mathrm{sliding}}}}=\\ & {\displaystyle \frac{{\mathit{\tau }}_{\max }}{W\left|\sin {\mathit{\alpha }}_{\mathrm{s}}\right|+H_x\cdot \cos {\mathit{\omega }}_{\mathrm{s}}\cdot \cos {\mathit{\alpha }}_{\mathrm{s}}\cdot k_{\mathrm{1}}\cdot k_{\mathrm{3}}+H_y\cdot \left|\sin {\mathit{\omega }}_{\mathrm{s}}\right|\cdot \cos {\mathit{\alpha }}_{\mathrm{s}}\cdot k_{\mathrm{1}}+E\cdot \cos {\mathit{\alpha }}_{\mathrm{s}}\cdot k_{\mathrm{2}}}}.\end{array}\end{array}$$

When the block can slide freely, α_s=α and ω_s=0; when the block is constrained to slide along a joint plane (e.g. J1), α_s=θ₁ or θ₂ and ω_s=ω₁ or ω₂. For the case of an anaclinal slope, the sliding direction is opposite to the free surface. Therefore, the rock block does not slide, and Fos_sl is not considered in the model.

With regard to stability against toppling, along the x direction, the part of the block above the mudstone base provides the stabilizing moment $M_{W_{{\mathrm{in}}_x}}$, and the part of the block above the basal cavity provides the overturning moment $M_{W_{{\mathrm{out}}_x}}$. When tension exists, there is an additional stabilizing moment. $M_{W_{{\mathrm{in}}_x}}$, $M_{W_{{\mathrm{out}}_x}}$ and $M_{{\mathrm{p}}_x}$ can be derived as

$$\begin{array}{}\text{(23)}& {\displaystyle }& {\displaystyle }{M}_{W_{{\mathrm{in}}_x}}=W{\displaystyle \frac{a-d_{\mathrm{1}}}{a}}\cos {\mathit{\theta }}_{\mathrm{1}}\cdot \left({\displaystyle \frac{a-d_{\mathrm{1}}}{\mathrm{2}}}\right),\text{(24)}& {\displaystyle }& {\displaystyle }{M}_{W_{{\mathrm{out}}_x}}=W{\displaystyle \frac{d_{\mathrm{1}}}{a}}\cos {\mathit{\theta }}_{\mathrm{1}}\cdot {\displaystyle \frac{d_{\mathrm{1}}}{\mathrm{2}}},\end{array}$$

$$\begin{array}{}\text{(25)}& M_{{\mathrm{p}}_x}=-\underset{-\frac{b-d_{\mathrm{2}}}{\mathrm{2}}}{\overset{\frac{b-d_{\mathrm{2}}}{\mathrm{2}}}{\int }}\underset{-\frac{a-d_{\mathrm{1}}-d_{\mathrm{3}}}{\mathrm{2}}}{\overset{\frac{a-d_{\mathrm{1}}-d_{\mathrm{3}}}{\mathrm{2}}}{\int }}{p}_{\mathrm{n}}(x,y)\cdot \left({\displaystyle \frac{a}{\mathrm{2}}}-d_{\mathrm{1}}-x\right)\mathrm{d}x\mathrm{d}y,\end{array}$$

and $M_{H_x}$ and $M_{E_x}$ can be derived as

$$\begin{array}{}\text{(26)}& {\displaystyle }\begin{array}{rl}{M}_{H_x}=& \underset{-\frac{b-d_{\mathrm{2}}}{\mathrm{2}}}{\overset{\frac{b-d_{\mathrm{2}}}{\mathrm{2}}}{\int }}\underset{\mathrm{0}}{\overset{h_{\mathrm{w}}\cos {\mathit{\theta }}_{\mathrm{1}}}{\int }}{\mathit{\gamma }}_{\mathrm{w}}\left(h_{\mathrm{w}}-{\displaystyle \frac{z}{\cos {\mathit{\theta }}_{\mathrm{1}}}}\right)\\ & \left({\displaystyle \frac{z}{\cos {\mathit{\theta }}_{\mathrm{1}}}}+\left(a-d_{\mathrm{1}}\right)\sin {\mathit{\theta }}_{\mathrm{1}}\right)\mathrm{d}z\mathrm{d}y,\end{array}\text{(27)}& {\displaystyle }{M}_{E_x}=E_x\left({\displaystyle \frac{h}{\mathrm{2}}}+\left({\displaystyle \frac{a}{\mathrm{2}}}-d_{\mathrm{1}}\right)\sin {\mathit{\theta }}_{\mathrm{1}}\right).\end{array}$$

Therefore, the Fos against toppling along the x direction, Fos${}_{{\mathrm{to}}_x}$, results in

$$\begin{array}{}\text{(28)}& {\mathrm{Fos}}_{{\mathrm{to}}_x}={\displaystyle \frac{M_{\mathrm{stabilizing}}}{M_{\mathrm{overturning}}}}={\displaystyle \frac{M_{W_{{\mathrm{in}}_x}}+M_{{\mathrm{p}}_x}}{M_{W_{{\mathrm{out}}_x}}+M_{H_x}\cdot k_{\mathrm{1}}\cdot k_{\mathrm{3}}+M_{E_x}\cdot k_{\mathrm{2}}}}.\end{array}$$

Similarly, Fos${}_{{\mathrm{to}}_y}$ can be obtained as

$$\begin{array}{}\text{(29)}& {\mathrm{Fos}}_{{\mathrm{to}}_y}={\displaystyle \frac{M_{\mathrm{stabilizing}}}{M_{\mathrm{overturning}}}}={\displaystyle \frac{M_{W_{{\mathrm{in}}_y}}+M_{{\mathrm{p}}_y}}{M_{W_{{\mathrm{out}}_y}}+M_{H_y}\cdot k_{\mathrm{1}}+M_{E_y}\cdot k_{\mathrm{2}}}}.\end{array}$$

The smaller value is selected as the Fos of the toppling failure mode Fos_to:

$$\begin{array}{}\text{(30)}& {\mathrm{Fos}}_{\mathrm{to}}=\min \left({\mathrm{Fos}}_{{\mathrm{to}}_x},{\mathrm{Fos}}_{{\mathrm{to}}_y}\right).\end{array}$$

When the stress on mudstone exceeds its strength, it causes partial damage and decreases the stability of the rock block.

Therefore, Fos with the consideration of compressive strength (Fos_co) and tensional strength (Fos_te) can be derived as

$$\begin{array}{}\text{(31)}& {\displaystyle }& {\displaystyle }{\mathrm{Fos}}_{\mathrm{co}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{cmax}}}{p_{\max }}},\text{(32)}& {\displaystyle }& {\displaystyle }{\mathrm{Fos}}_{\mathrm{te}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{tmax}}}{-p_{\min }}}.\end{array}$$

Fos_co and Fos_te represent the current damage degree of mudstone due to compressive stress and tensile stress, respectively. When the stress exceeds the ultimate strength, the strength of the mudstone is reduced to the residual value, and the initial deformation appears. The ability of mudstone to provide resistance to the sliding and toppling of sandstone blocks is thus reduced, and Fos_sl and Fos_to subsequently decline. The smaller the value of Fos_co and Fos_te, the greater the damage to the underlying mudstone. The effective contact area between sandstone and mudstone becomes smaller as compressive and tension damage develop, which significantly affects the stability of the overhanging sandstone block.

Summarizing, four Fos types of unstable rock block are obtained. Fos_sl and Fos_to are routine indicators directly representing the stability of sandstone blocks. Fos_co and Fos_te are two indicators proposed in this study for the stability analysis of biased rockfall, which describe the damage state of the underlying mudstone base. It is necessary to simultaneously consider four Fos types to evaluate the stability of unstable biased rockfall. The entire calculation process is shown in Fig. 9.

Figure 9 Calculation process of the Fos of the unstable rock blocks.

4 Validation of analytical methods by numerical simulation

The damage mechanisms at the base of the rock block play an important role in the rockfall evolution process. However, the stress distribution on the contact surface calculated by the proposed analytical methods is difficult to validate with the field data. Therefore, numerical simulation of a biased rockfall was conducted in this study to determine the stress distribution on the contact surface between overhanging sandstone and underlying mudstone. Numerical simulations can take into account material deformation, unlike the analytical methods. FLAC3D, professional software that utilizes the finite-difference method (FDM) for three-dimensional analysis of rocks, soils and other materials, was employed for the 3D numerical simulation. Based on the geological models, a 3D numerical simulation model was conducted with FLAC3D 6.00 to analyse the stress distribution on the contact surface (Fig. 10).

Figure 10 Numerical model built in FLAC3D.

The model is mainly composed of sandstone and mudstone, where “Overhanging sandstone 1” represents an unstable rock block (dimensions $a\times b\times h$ are 6, 8 and 10 m, respectively), and the weathering process of the mudstone is represented by excavating in stages in the underlying mudstone. Considering the deformation of the materials, sandstone was considered an elastic model and mudstone was assigned as a Mohr–Coulomb model. Material properties were determined by referring to published literature and investigation reports in the study area. The sandstone has a unit weight (γ_s) of 25 kN m⁻³, an elastic modulus of 4530 MPa and a Poisson's ratio of 0.183 (Tang et al., 2010), and the mudstone has a unit weight of 22.54 kN m⁻³, an elastic modulus of 587 MPa and a Poisson's ratio of 0.232. The friction angle of the contact surface (φ) is set to 25^∘, and the cohesion (c) is set to 70 kPa (Zhang et al., 2016). Because of the strength degradation of mudstone foundations due to intense weathering, the maximum compressive stress of mudstone (σ_cmax) is replaced by the bearing capacity of mudstone foundations (2300 kPa), which is obtained through plate load tests in adjacent areas (Zhang et al., 2021). In addition, the maximum tensile stress of mudstone (σ_tmax) is valued as $\mathrm{1}/\mathrm{9}$ of σ_cmax. The west, north and bottom boundaries of the model are constrained by roller boundary conditions. The cohesion and internal friction angle of the interface between “Overhanging sandstone 1” and “Overhanging sandstone 2” are set to 0. After reaching the initial force–equilibrium state, the mudstone was excavated to simulate the weathering process, and the vertical stress distribution on the sand–mudstone interface at different basal cavity depths was obtained, as shown in Fig. 11.

Figure 11 Diagram of stress distribution in the vertical direction on the contact interface through different methods: (a) the results of numerical simulation by FLAC3D and (b) the results of the proposed analytical method. SVD stands for singular value decomposition.

When there is no cavity present, represented by d=0 m, the stress distribution is uniform compressive stress (according to the FLAC3D software, compressive stresses are negative). At d=0.5 m, the stress remains entirely compressive, but non-uniform stress distribution occurs on the contact surfaces. At d=1 m, the vertical stress value in the upper-left corner of the contact interface surpasses 0 (Fig. 11), indicating the presence of tensile stress. As d increases to 1.5 or 2 m, the tensile stress in the upper-left corner gradually intensifies, exacerbating the non-uniform stress distribution. The results obtained from the numerical simulation align with those from the analytical method, confirming the existence of tensile stress at the contact interface in the biased rockfall due to external erosion development (Fig. 11). Tensile stress commonly emerges within the contact surface, making it challenging to observe directly in the field.

In the context of the limit equilibrium method, the contact area plays a vital role in stability analysis, as shown in Eqs. (21)–(30) in Sect. 3. The numerical simulation process provides an intuitive understanding of the influence of non-uniform stress distribution on the contact surfaces on the stability of rock blocks. Whether subjected to tension or compression, the rock layer has an ultimate strength. In Fig. 11, when d=1.5 or 2 m, the tensile stress exceeds the ultimate tensile strength, leading to tension failure in the upper-left corner of the stress distribution diagram. The region enclosed by a dashed yellow line represents ineffective contact, where no anti-slip force or overturning moment can be generated due to tension failure at the contact surface. Therefore, this area needs to be subtracted from the total contact area when calculating Fos_sl and Fos_to. Similar situations occur when the compressive stress exceeds the ultimate compressive strength. The current maximum compressive stress has not reached the ultimate compressive strength in Fig. 11. However, as d continues to increase, the area of compression failure will appear in the lower-right corner of diagram in Fig. 11. This occurrence diminishes the area capable of providing anti-slip force or overturning moment, thereby reducing the stability of the rock blocks.

The traditional LEM does not account for distributed forces and fails to consider changes in the contact surface. The method proposed in this study addresses this issue and is applied to the calculation of Fos_sl and Fos_to as presented in Eqs. (21), (25) and (26).

5 Results

A detailed field investigation was carried out in the source area of rockfall (Fig. 3d). The size of the blocks was determined by on-site measurement with tape and a laser rangefinder. The external erosion in mudstone was measured with a steel ruler, and the morphological characteristics of mudstone foundation were mainly described with the average erosion depth of the basal cavity. The attitude of discontinuities was measured by compass. The mechanical parameters are given in Sect. 4. The height of the water level (h_w) is set to be $\mathrm{1}/\mathrm{3}$ of h, and an earthquake contribution coefficient k_e of 0.05 is considered in stability calculations. The data obtained from the field survey were organized according to the coordinate system of the geological model in Sect. 3.1, and Fos was calculated according to the calculation steps in Sect. 3.2. The calculated geometric parameters and Fos results are shown in Table 2.

Table 2 Geometric parameters of rock blocks in the study area and Fos results.

Note: when there is no tensile stress in the mudstone foundation, Fos_te has no value. For the case of an anaclinal slope, blocks do not slide and Fos_sl has no value. Both parameters are replaced by “–”.

6 Discussion

6.1 Characteristics of rock block stability

There are up to 12 results of Fos per potentially unstable block with the consideration of three scenarios and four failure modes (i.e. partial damage and overall failure). Most Fos_te values are less than 1 in all scenarios (yellow points in Fig. 12), except for two blocks (i.e. W17 and W20), whose Fos_te values are also close to 1 under rainfall or earthquake scenarios. Although most of the Fos_co values (green points in Fig. 12) are greater than 1, they are closer to the critical state of Fos=1 than Fos_sl and Fos_to (represented by blue and orange points in Fig. 12, respectively). The compression damage of the exposed mudstone can be investigated in the field survey (Fig. 4d). However, it is difficult to observe the phenomenon of tensile damage inside the mudstone base. In the case of weak tensile strength, the mudstone base suffers from tension failure, and compression failure usually occurs before tension failure. According to the results, the Fos_te and Fos_co of these blocks are less than 1 or close to 1, which means that the underlying mudstone has been partially damaged due to slight compressive or tension failure, and the blocks are potentially unstable with the current depth of the basal cavity. However, most of the blocks do not exhibit overall failure, and they still exist on the slope. Moreover, their Fos_sl and Fos_to values are greater than 1 in different scenarios, which is consistent with this actuality. The results indicate that most of the blocks are close to a critical state, in which they are partially damaged but the whole block is still stable.

Figure 12 Distribution of Fos values in different scenarios. Shapes represent different scenarios, and colours represent different failure modes.

6.2 Relationship between Fos and geometric parameters

Figure 13 presents the relationship between Fos_min and two main geometric parameters, the dip of the contact surface and the retreat ratio. In general, the dip angle of the contact surface (α) is the key factor influencing the sliding failure mode. The horizontal axis in Fig. 13a is α between the rock blocks and underlying mudstone. Most of the points in Fig. 13a are in the interval [0, 8^∘], which is consistent with the features of sub-horizontal strata in the study area. The shade of the points does not change significantly in the x-axis direction, as Fig. 13a shows. Therefore, compared with the maximum retreat ratio (r_max), the dip of the contact surface has less influence on rockfall stability in the study area. There was a significant positive correlation between the retreat ratio (r_max) and Fos_min. In Fig. 13b, as the retreat ratios increase in the positive direction of the x axis and y axis, the rock blocks show a notable tendency to be unstable.

Figure 13 Correlation between Fos values and the dip of contact surface and retreat ratio. Here, α is the dip angle of the contact surface between the rock block and underlaying mudstone; r_x and r_y are the retreat ratio along the x direction and y direction, respectively, equal to $d_{\mathrm{1}}/a$ and $d_{\mathrm{2}}/b$; and r_max is the larger of r_x and r_y.

6.3 Definition of rockfall susceptibility

To explore the variation in Fos values with the progressive erosion process of the mudstone on the blocks, the mudstone retreat velocities in different directions are assumed to be equal (5 mm yr⁻¹; Zhang et al., 2016). Figure 14 shows the variations in Fos values of two specific blocks during the evolution process of the basal erosion in the mudstone. In the initial stage, the basal cavity is small, and the overhanging block is stable; all Fos values are greater than 1.0. The basal cavity expands over time as the mudstone weathers; then, the contact area decreases, and non-uniform distributed stress arises. When the stress exceeds the ultimate strength of mudstone in a partial area, Fos_co and Fos_te decrease significantly, as shown in Fig. 14. The instability of the blocks starts from the failure (or damage) of the foundation. Fos_te and Fos_co reach the critical state much earlier than Fos_sl and Fos_to. For these two specific blocks, when r_max increases to 0.4, Fos_sl and Fos_to are still higher than 1.0. This means that the rock blocks can remain globally stable in this condition.

Figure 14 Variation in Fos values with r_max. Panel (a) and (b) are the results for W01 and W09, respectively, which represent the situation of the blocks with two and three free faces. The dashed black line (CS) approximately represents the current state of the unstable blocks. The dashed red lines correspond to the critical values of r in different scenarios.

Figure 15 Rockfall susceptibility based on the combination of four Fos types. The susceptibility is defined in three levels, represented by red, yellow and blue. Panel (a) shows the progressive failure process of the rock block changing from low susceptibility to moderate susceptibility as the mudstone retreat ratio increases (illustrated by W01 and W09 in the natural scenario. Panels (b) and (c) show the change in susceptibility of biased rock blocks, when the scenario changes from natural conditions to rainfall and earthquake conditions.

These results further elucidate the stability analysis model proposed in this study. Fos_co and Fos_te introduced in this model present the damage state of basal mudstone caused by compressive and tensile stresses, which do not provide global instability of the overhanging block as sliding and toppling. However, Fos_co and Fos_te are important preliminary signs of subsequent global failure of the rock block, as presented through the numerical simulation in Sect. 4. The damage in the basal mudstone can significantly accelerate weathering and prompt expansion of the basal erosion, which will lead to global failure. The lower Fos_co and Fos_te are, the lower the safety margin of the blocks. Therefore, the four Fos types used in this study can provide a more comprehensive quantification of rockfall stability.

This result is consistent with Fig. 12, in which 63.7 % of the yellow and green points (Fos_te and Fos_co) are located between Fos=0.7 and Fos=2.0. This result can be validated by the field phenomena. In the study area, rock damage (e.g. micro-fractures and cleavages) can be observed in the underlying mudstone. However, most overhanging rock blocks are stable at the present time. This means that even if Fos_sl or Fos_to is higher than 1, its foundation has begun to be damaged. In the case of heavy rain or earthquakes, Fos_sl and Fos_to may be reduced to less than 1, and the rockfall occurs.

Based on the meaning of four Fos types, rockfall susceptibility can be divided into three levels. When both Fos_co and Fos_te are greater than 1, the overall rock block is stable and the mudstone base is not damaged, which is defined as “low susceptibility” and represented by the blue area in Fig. 15. With the development of basal erosion, when Fos_co or Fos_te is less than 1 and Fos_sl and Fos_to are higher than 1, the base undergoes damage and the overhanging sandstone blocks remain relatively stable. This state is defined as “moderate susceptibility” and is represented by the yellow area. When Fos_sl or Fos_to is less than 1 in some scenarios, the rock blocks are in a state of “high susceptibility”, which means that rockfalls are highly likely to occur. Figure 15a indicates that along with the increase in the mudstone retreat ratio, the susceptibility of W01 and W09 changes from low susceptibility to moderate susceptibility in the natural scenario. As Fig. 15b and c show, when rainfall or earthquake occurs, the Fos_sl or Fos_to of some blocks is less than 1, which means that some blocks have evolved to the state of high susceptibility and the overall sandstone blocks are unstable.

6.4 Critical retreat ratio in the study area

The basal erosion in the mudstone plays an important role in the progressive failure process of biased rockfall. To analyse the effect of the retreat ratio on the stability of rock blocks, all blocks in the study area were selected to calculate their Fos values and susceptibility level with the increasing r, whose retreat velocities in different directions are assumed to be equal. Figure 16 shows that along with the increase in the retreat ratio, the susceptibility level of rock blocks changes from low to moderate susceptibility. Corresponding to the critical state of $\min \mathit{\{}{\mathrm{Fos}}_{\mathrm{co}},{\mathrm{Fos}}_{\mathrm{te}}\mathit{\}}=\mathrm{1}$ of all blocks, the minimum retreat ratio is 0.26, and the maximum retreat ratio is 0.41, as marked by the vertical dashed black line in Fig. 16. According to the statistical analysis of critical retreat ratios, both the mean and the median are 0.33. Therefore, the critical retreat ratio of the rock blocks in the study area can be determined as 0.33, which is marked by the vertical dashed red line in Fig. 16. The critical retreat ratio calculated by this method can be used for the preliminary identification of potentially unstable rock blocks in a specific area, which can help concentrate limited risk treatment resources on these priorities. It should be emphasized that the mechanical parameters and analysis scenarios significantly affect the critical value. Therefore, the elaborative risk control of a given rockfall should be arranged based on its specific parameters and analysis scenarios.

Figure 16 Effect of the retreat ratio (r) on the Fos of the rock block, which is illustrated by all blocks in the study area.

6.5 Limitations

This study involves the development of an analytical model for the three-dimensional stability of biased rockfall, combining the basic LEM and the consideration of the non-uniform distribution. Due to the complexity of rock structure and force analysis, it is necessary to highlight the limitations of this model.

First, this study uses a three-dimensional coordinate system and bending theory. It is difficult to consider diverse shapes of rock blocks, and the rock block was simplified as a prismatic column. The assumption of fully persistent discontinuities may underestimate the stability of rock blocks and ignores the stress transmission in joints or rock bridges. Then, following the basic framework of the general LEM, this study assumed that the rock is not subjected to deformations. The complete stress–strain behaviour, such as the deformation in the mudstone layer, is not considered in this study. The mode of tension failure is very difficult to observe in the field, and it is currently verified by means of numerical simulation. Furthermore, the block stability is strongly influenced by the uncertainty in mechanical parameters. However, because of the difficulties in sampling strongly weathered mudstone, it is difficult to obtain adequate parameter values for uncertainty statistics. These limitations will be important considerations in future studies.

7 Conclusions

Due to differential weathering in sub-horizontally interbedded hard rock and soft rock, multi-layer biased rockfalls develop on steep slopes. In mountainous ranges, cut slopes and coastal cliffs, rockfall may cause significant facility damage and casualties in residential areas and transport corridors. The aim of this study was to present a new three-dimensional analytical method for the stability of rock blocks with basal cavities. In this method, a non-uniform distributed stress due to the eccentricity effect is applied at the contact surface instead of a point force. The development of non-uniform distributed stress calculated by the proposed analytical methods was validated by numerical simulation, which presents the evolution process of biased rockfall from partial damage of the soft underlying layer, caused by non-uniform distributed stress, to toppling and sliding of overhanging hard rock block due to overall unbalanced force. The method considers four failure modes according to the rockfall evolution process, including partial damage of the soft foundation (Fos_co and Fos_te) and overall failure of the rock block (Fos_sl and Fos_to).

Taking the northeastern edge of the Sichuan Basin in southwest China as the study area, the proposed method is used to calculate the Fos of biased unstable rock blocks. The results show that in the natural scenario, the underlying mudstone of some rock blocks has been partially damaged, and compression failure of the mudstone has been observed in the field. Some rock blocks are expected to fail as a whole in rainfall or earthquake scenarios. The statistical analysis indicates that the retreat ratio is the crucial factor influencing the Fos of biased rockfall. On the basis of different combinations of four Fos types, rockfall susceptibility was classified into three levels. As the retreat rate increases, the rock blocks undergo an evolution process from stability to partial instability and then overall instability. Based on the current mechanical parameters of the eastern Sichuan Basin, the critical retreat ratio from low to moderate rockfall susceptibility is 0.33.

The proposed method improves the three-dimensional mechanical model of a rock block with the basal erosion by considering non-uniform distributed stress at the contact surface, which could promote the accuracy of rockfall stability analysis. Due to the assumptions adopted and the complexity of the failure mechanism of biased rockfall, there are some limitations in this method, mainly including the simplification of boundary conditions and rock deformation. These limitations will be important considerations in future studies.

Appendix A: List of symbols

a	length of the block along the x direction
A	area of contact surfaces
b	width of the block along the y direction
c	cohesive force of the mudstone
di	width of the basal cavity in a certain direction
ex	eccentric distance along the x direction
ey	eccentric distance along the y direction
Ex	horizontal seismic force along the x direction
Fos	factor of safety
h	height of the block
hw	height of the water in the fracture
Hx	water pressure along the x direction
Ix	moment of inertia with respect to the x axis
Iy	moment of inertia with respect to the y axis
ke	earthquake contribution coefficient

k1	rainfall coefficient, taking 1 in the rainfall
	scenario and 0 in the non-rainfall scenario
k2	earthquake coefficient, taking 1 in the seismic
	scenario and 0 in the non-seismic scenario
k3	free-surface coefficient, taking 1 for two free
	surfaces and 0 for three free surfaces
$M_{{\mathrm{b}}_x}$	total bending moments with respect to the
	x axis on the mudstone foundation
$M_{{\mathrm{b}}_y}$	total bending moments with respect to the
	y axis on the mudstone foundation
$M_{\mathrm{b}{E}_x}$	bending moment of Ex with respect to the
	x axis on the mudstone foundation
$M_{\mathrm{b}{H}_x}$	bending moment of Hx with respect to the
	x axis on the mudstone foundation
$M_{\mathrm{b}{W}_x}$	bending moment of W with respect to the
	x axis on the mudstone foundation
$M_{E_x}$	overturning moment provided by Ex
	along the x direction
$M_{H_x}$	overturning moment provided by Hx along
	the x direction
$M_{{\mathrm{p}}_x}$	stabilizing moment of pn along the
	x direction
$M_{W_{{\mathrm{in}}_x}}$	stabilizing moment provided by W along
	the x direction
$M_{W_{{\mathrm{out}}_x}}$	overturning moment provided by W along
	the x direction
Nz	total applied vertical load on the
	mudstone base
O	origin of the (x, y) coordinates
p(x,y)	pressure magnitude at point (x, y)
ri	the basal mudstone retreat ratio equal to
	the ratio of the basal cavity width to the block
	width in a certain direction
W	weight of the block
x	distance to O along the x axis
y	distance to O along the y axis
α	true dip of the contact surface
γs	unit weight of sandstone
γw	unit weight of water
θ1	apparent dip of α on plane J1
θ2	apparent dip of α on plane J2
σcmax	ultimate compressive strength of the
	mudstone
σtmax	ultimate tensile strength of the
	mudstone
τmax	ultimate shear strength of the
	mudstone
φ	friction angle of the mudstone
ω1	angle between the trend of the contact
	surface and the x direction
ω2	angle between the trend of the contact
	surface and the y direction