# Cylindrical Caved Space Stability Analysis for Extension Prediction of Mining-Induced Surface Subsidence

Authors: Yang Liu; Yongxiang Ge; Congrui Zhang; Fengyu Ren; Junsheng Ma; Gaofeng Ren
Journal: Geofluids
 (2021)

---

## Abstract

Subsequent extension of surface subsidence after vertical caving leads to large-scale surface destruction, as well as associated geological hazards. The extension prediction for cylindrical caved space, which appears circular surface subsidence, is still an intractable issue, due to the absence of robust models. To fill such a research gap, this paper provides an analytical model for the depth and orientation where the shear failure of isotropic rocks around the caved space is firstly observed. The anisotropy of surrounding rocks is further involved to enable this model to analyze the slip failure along discontinuities in anisotropic stress state. The prediction for the extension of the surface subsidence in Xiaowanggou iron mine is conducted, and the comparison between the prediction and the observation in satellite images demonstrates the validity of the proposed model. Even though this model cannot provide a definite boundary after extension, the prediction for the orientation surface subsidence extends to contribute to mitigating the effect of geological hazards. Another contribution of this work is to provide guidance to mitigate the impact of surface subsidence on safety and environment, such as filling the interspace between large-sized caved rocks by dumping small-sized waste rocks or backfilling the caved space with waste rocks.


---

## Body

## 1. Introduction

The geological hazards due to mining-induced surface subsidence raise public concerns [1, 2]. The employment of caving mining is restricted because of the surface subsidence [3, 4], even though it is being favored in underground mining projects due to high cost efficiency [5]. On the other hand, in situ observations have demonstrated that the surface subsidence is likely to take place in open-stope-based mines [6], even in the filling-based ones [7, 8]. Such mining-induced surface subsidence can be categorized into 6 types, such as crown hole, chimney caving, plug subsidence, solution cavities, block caving, and progressive caving [9]. The progressive caving indicates the process that surrounding rocks progressively cave into the original caved space formed after chimney caving, plug subsidence, or block caving; meanwhile, the caved space will be extended accordingly during such process. In situ or remote monitoring from many mining projects shows the progressive caving is a primary contributor to the large-scale extension of surface subsidence [10–13].To predict the progressive caving, as well as the extension of surface subsidence, some analytical solutions have been proposed. The most commonly cited is the limiting equilibrium analysis by Hoek [14]. This model relates the stability of the hanging wall of caved space to the effective cohesion of rock mass, thrust on shear plane due to caved rocks, and the base area and weight of wedge of sliding rock mass. Brown and Ferguson [15] extended Hoek’s model by involving the impact of water pressure and sloping ground surface on shear failure plane. Lupo [16] modified Hoek’s model to enable it to predict the rock failure on both hanging wall and footwall. Yet despite having been successfully applied in many mining projects, such as Kiruna iron mine in Sweden, El Teniente copper mine in Chile, Río Blanco copper mine in Peru, and Gaths asbestos mine in Zimbabwe, the employment of both Hoek’s and extended models is still restricted because of some basic assumptions: (1) the rock failure and associated surface subsidence occur for a long distance compared with the cross section normal to the strike of ore body, and the analysis can be reduced to a two-dimension problem; (2) the rock mass has homogeneous and isotropic mechanical property, which means the discontinuities (e.g., faults, joints, or beddings) are not taken into account. Such assumptions mean Hoek’s and extended models are invalid to analyze progressive caving either around the cylindrical caved space or with the consideration of the rock anisotropy due to discontinuities.Both in situ monitoring [17, 18] and numerical simulations [7, 19, 20] show the evident impact of discontinuities on the progressive caving. On the other hand, the cylindrical surface subsidence formed after chimney caving, plug subsidence, or block caving has been observed not only in caving-based projects (e.g. Northparkes copper-gold mine, Kimberley diamond mine, and Xiaowanggou iron mine; Figure 1) but also in filling-based one (e.g., Jinchuan nickel mine) [9]. However, the analytical prediction for the progressive caving, as well as the associated extension of surface subsidence in these projects, is still an intractable issue, due to the absence of robust model, especially when the rock anisotropy due to discontinuities is taken into consideration.Figure 1
Circular surface subsidence observed visually in Google Earth satellite images: (a) Northparkes copper-gold mine, Australia; (b) Kimberley diamond mine, South Africa; (c) Xiaowanggou iron mine, China.
(a)(b)(c)Therefore, to fill such a research gap, we introduce an analytical model to relate the rock shear failure around the cylindrical caved space, which accounts for extension of surface subsidence, to in situ stress, and the property of surrounding rocks and caved rocks. Additionally, the rock anisotropy due to inherent discontinuities is involved to enable the proposed model to analyze the extension of surface subsidence because of the slip failure along discontinuities. The rest of this paper is organized as follows. Section2 introduces the solutions to test stability of the isotropic rocks around cylindrical caved space. Section 3 provides the solution to predict the slip failure along discontinuities around the caved space. In Section 4, employing the proposed model, we calculated the depth and orientation where the shear or slip failure take place at the intact rocks around caved space in Xiaowanggou iron mine (Figure 1(c)), and the associated extension of surface subsidence is predicted. Additionally, some implications regarding the proposed model and case study are discussed in this section. Finally, the conclusions and contributions are summarized in Section 5.


## 2. Stability Analysis for Surrounding Isotropic Rocks of Cylindrical Caved Space

The stability analysis for the surrounding isotropic rocks of cylindrical caved space is primarily to establish a relationship to test whether rock failure will take place, by comparing the redistributed stress after vertical caving (e.g., chimney caving, plug subsidence, or block caving) with rock strength by an appropriate failure criterion. To conduct such analysis, the caved space, observed as circular subsidence area on ground surface (Figure1), is assumed as a vertical and cylindrical excavation filled with caved rocks (Figure 2(a)).Figure 2
Coordinates’ transformation from in situ stress (σH, σV, σh) into the local stress at surrounding intact rocks (σθ, σz, σr, τθz, τrθ, τrz) in polar system (θ, z, r). (a) Three-dimension (3D) view of the caved space. (b) Local stress at the surrounding rocks in the cross section of caved space. (c) An element illustrating the local stress at surrounding rocks.When the rocks situated in the in situ stress state caves in the void after deposit excavation, the stress redistribution takes place near the caved space. To analyze the rock shear failure surrounding rocks of caved space, the following inputs are involved, such as rock strength, in situ stress, and the impact of caved rocks. Conventionally, rock strength, in situ stress, and associated orientations can be obtained by some robust methods [21–24]. However, the analysis for the impact of caved rocks is rarely reported.


### 2.1. Impact of Caved Rocks on Stability of Cylindrical Caved Space

In the analytical models proposed in existing literature [14–16], the caved rocks provide thrust to the sliding rocks, which contribute to enhancing the stability of surrounding rocks of caved space. Under the assumption that the subsidence occurs for a long distance compared with the cross section normal to the strike of ore body, this thrust can be expressed as a function of the density and height of caved rocks. Hence, numerous literatures have demonstrated such contribution; Laubscher [13] proposed an empirical solution to describe the effect of the density and height of caved rocks. To describe the stress from caved rocks to the surrounding rocks in a confined space, Ren et al. [25] conduct the scaled laboratory tests, and the results reveal this stress matches the solution for granular materials by Janssen [26]:
(1)σb=Cρbgrb1−e−zb/4rb,where σb is the horizontal stress from caved rocks to the surrounding rocks, ρb is the average density of caved rocks, g is the gravity of earth, rb is the radius of the cross section of caved space, zb is a variable representing the depth measured from the free surface of caved rocks, and C is a constant that can be obtained by in situ or laboratory tests.


### 2.2. Shear Failure at Isotropic Surrounding Rocks

For the cylindrical excavation (Figure2(a)), the in situ stress and redistributed local stress are illustrated in Figures 2(b) and 2(c), respectively. Such redistributed local stress near the rocks around the cylindrical caved space can be expressed in polar system (θ, z, r), original due to Krisch:
(2)σr=σH+σh21−rb2r2+σH−σh21−4rb2r2+3rb4r4cos2θ+σbrb2r2,σθ=σH+σh21+rb2r2−σH−σh21+3rb4r4cos2θ−σbrb2r2,σz=σV−2vσH−σhrb2r2cos2θ,τrz=τrθ=τθz=0,where σr, σθ, and σz are the radial, tangential, and axial local normal stress near the rocks around the cylindrical excavation, respectively, τrθ, τrz, and τθz are the radial, tangential, and axial local shear stress, respectively, and σV, σH, and σh are the vertical, maximum, and minimum horizontal in situ stress, respectively; the in situ stress linearly varies with the depth from ground surface (z) and can be expressed as σV=aVz+bV, σH=aHz+bH, and σh=ahz+bh; aV, bV, aH, bH, ah, and bh are constants which can be obtained by regressing the results from in situ test (e.g., overcoring and hydraulic fracturing) or core-based test (e.g., anelastic strain recovery, differential strain curve analysis, and acoustic emission); r is a variable that indicates the distance between the position stress redistribution that occurs and the axis of the cylindrical excavation; θ is the orientation around the cylindrical excavation measured from the direction of maximum horizontal in situ stress (i.e., σH, as illustrated in Figure 2(b)); θ=0° represents the direction of maximum horizontal in situ stress (σH), and θ=90° represents the direction of minimum horizontal in situ stress (σh); v is Poisson’s ratio of the surround rocks.The normal and shear stress at caved space wall (r=rb) can be accordingly calculated from Equation (2):
(3)σr=σb,σθ=σH+σh−2σH−σhcos2θ−σb,σz=σV−2vσH−σhcos2θ,τrz=τrθ=τθz=0.Equation (3) reveals that the relationship between the three principal stress (σθ, σz, and σr) depends on not only the constants from in situ or laboratory test (e.g., aV, bV, aH, bH, ah, bh,ρb, and rb) but also the depth (z) and direction (θ) at the surrounding rocks. Even though the tangential principal stress (σθ) is the maximum principal stress, it is intractable to determine the relationship between the radial and axial principal stress (i.e., σr and σz, respectively). For instance, in shallow depth, the principal stress is likely to satisfy σθ>σr>σz. However, when the depth increases, the relationship of σθ>σz>σr can be satisfied, because σb will be near to its ultimate value (limzb→+∞σb=Cρbgrb). Therefore, to analyze the rock failure around the cylindrical caved space in varying depths, we employ the extended Mohr-Coulomb failure criterion for the problems in 3D space [27]
(4)σc+qσ3−σ1≥0,σc+qσ1−σ3≥0,σc+qσ2−σ1≥0,σc+qσ1−σ2≥0,σc+qσ3−σ2≥0,σc+qσ2−σ3≥0,where σ1, σ2, and σ3 are the maximum, intermediate, and minimum applied principal stress, respectively; σc is rock strength; σc can be replaced by long-term strength of rocks (σcd), if this criterion is utilized to analyze rock stability in a long duration of time [28]; the long-term strength of rocks (σcd) can be converted from the rock strength without time effect (σc′) by σcd=ασc′; α=0.4~0.6, for the intact rocks without preexisting fractures [29]; q=1+sinϕ/1−sinϕ; ϕ is the internal friction angle of rocks.Assuming the tangential, radial, and axial principal stress (i.e.,σθ, σr, and σz) are the maximum, intermediate, and minimum applied principal stress (i.e., σ1, σ2, and σ3), respectively, by substituting Equation (3) into Equation (4), noticing σ1=σθ, σ2=σr, and σ3=σz, we obtain the following relationships, which should be satisfied, if the surrounding rocks are expected to maintain stable in a long duration of time
(5)σcd+qσV−2vσH−σhcos2θ−σH+σh−2σH−σhcos2θ−σb≥0,σcd+qσH+σh−2σH−σhcos2θ−σb−σV−2vσH−σhcos2θ≥0,σcd+qσb−σH+σh−2σH−σhcos2θ−σb≥0,σcd+qσH+σh−2σH−σhcos2θ−σb−σb≥0,σcd+qσV−2vσH−σhcos2θ−σb≥0,σcd+qσb−σV−2vσH−σhcos2θ≥0.Because the tangential principal stress (σθ) is the maximum applied principal stress (i.e., either σθ>σr>σz or σθ>σz>σr is satisfied), Equation (5) can be reduced to the following equations:
(6a)σcd+qσV−2vσH−σhcos2θ−σH+σh−2σH−σhcos2θ−σb≥0,(6b)σcd+qσb−σH+σh−2σH−σhcos2θ−σb≥0.SubstitutingσH=aHz+bH, σh=ahz+bh, σV=aVz+bV, σcd=ασc′, and Equation (1) into Equations (6a) and (6b), we obtain the following equations to test the stability of the rocks around the cylindrical excavation:
(7a)σcd+Cρbgrb1−e−z−z0/4rb+qaV−aH+ah+21−qvaH−ahcos2θz+qbV−bH+bh+21−qvbH−bhcos2θ≥0,(7b)σcd+1+qCρbgrb1−e−z−z0/4rb−aH+ah−2aH−ahcos2θz+2bH−bhcos2θ−bH+bh≥0,where z is a variable that indicates the depth from ground surface and z0 is the depth of the free surface of caved rocks from ground surface.Equations (7a) and (7b) enable the prediction for the depth (z) where rock failure will take place around the cylindrical excavation (θ). If the depth of implemented undercut exceeds the results by Equations (7a) and (7b), rock failure and associated extension of surface subsidence can be expected to take place. On the other hand, Equations (7a) and (7b) also reveal the predicted depth for rock failure varies in different orientations. For instance, when σθ>σz>σr is satisfied (i.e., Equations (6b) and (7b) are valid to conduct such prediction), rock failure and associated extension of surface subsidence are most likely to take place in the direction of minimum in situ stress, i.e., θ=90° is satisfied. This means the extension of surface subsidence due to the rock shear failure can be expected to be anisotropic.On the other hand, Equations (6a) and (6b) reveal that the stability of surrounding rocks heavily depends on not only the rock strength (σcd) but also the stress due to caved rocks (σb). Equations (7a) and (7b) provide more details regarding the contribution of the caved rocks. Because of q>0 and C>0, the increase of the density (i.e., ρb increases) or height (i.e., z0 decrease) of caved rocks facilitates the stability of surrounding rocks. This provides some measures to prevent the rock failure and associated extension of surface subsidence: (1) filling the interspace between large-sized caved rocks by dumping small-sized waste rocks to increase the average density of the materials in the caved space or (2) backfilling the caved space with waste rocks to increase the height of caved rocks. Some discussions for such measures should be noted. Equation (1) shows the ultimate value of stress due to caved rocks (limzb→+∞σb=Cρbgrb) will increase, if measure (1) is implemented. This means that this measure contributes to the stability of surrounding rocks in all depth. However, the effect of measure (2) only works in shallow depth, because the ultimate value of the horizontal stress due to caved rocks maintains constant when the height of caved rocks changes.


## 2.1. Impact of Caved Rocks on Stability of Cylindrical Caved Space

In the analytical models proposed in existing literature [14–16], the caved rocks provide thrust to the sliding rocks, which contribute to enhancing the stability of surrounding rocks of caved space. Under the assumption that the subsidence occurs for a long distance compared with the cross section normal to the strike of ore body, this thrust can be expressed as a function of the density and height of caved rocks. Hence, numerous literatures have demonstrated such contribution; Laubscher [13] proposed an empirical solution to describe the effect of the density and height of caved rocks. To describe the stress from caved rocks to the surrounding rocks in a confined space, Ren et al. [25] conduct the scaled laboratory tests, and the results reveal this stress matches the solution for granular materials by Janssen [26]:
(1)σb=Cρbgrb1−e−zb/4rb,where σb is the horizontal stress from caved rocks to the surrounding rocks, ρb is the average density of caved rocks, g is the gravity of earth, rb is the radius of the cross section of caved space, zb is a variable representing the depth measured from the free surface of caved rocks, and C is a constant that can be obtained by in situ or laboratory tests.


## 2.2. Shear Failure at Isotropic Surrounding Rocks

For the cylindrical excavation (Figure2(a)), the in situ stress and redistributed local stress are illustrated in Figures 2(b) and 2(c), respectively. Such redistributed local stress near the rocks around the cylindrical caved space can be expressed in polar system (θ, z, r), original due to Krisch:
(2)σr=σH+σh21−rb2r2+σH−σh21−4rb2r2+3rb4r4cos2θ+σbrb2r2,σθ=σH+σh21+rb2r2−σH−σh21+3rb4r4cos2θ−σbrb2r2,σz=σV−2vσH−σhrb2r2cos2θ,τrz=τrθ=τθz=0,where σr, σθ, and σz are the radial, tangential, and axial local normal stress near the rocks around the cylindrical excavation, respectively, τrθ, τrz, and τθz are the radial, tangential, and axial local shear stress, respectively, and σV, σH, and σh are the vertical, maximum, and minimum horizontal in situ stress, respectively; the in situ stress linearly varies with the depth from ground surface (z) and can be expressed as σV=aVz+bV, σH=aHz+bH, and σh=ahz+bh; aV, bV, aH, bH, ah, and bh are constants which can be obtained by regressing the results from in situ test (e.g., overcoring and hydraulic fracturing) or core-based test (e.g., anelastic strain recovery, differential strain curve analysis, and acoustic emission); r is a variable that indicates the distance between the position stress redistribution that occurs and the axis of the cylindrical excavation; θ is the orientation around the cylindrical excavation measured from the direction of maximum horizontal in situ stress (i.e., σH, as illustrated in Figure 2(b)); θ=0° represents the direction of maximum horizontal in situ stress (σH), and θ=90° represents the direction of minimum horizontal in situ stress (σh); v is Poisson’s ratio of the surround rocks.The normal and shear stress at caved space wall (r=rb) can be accordingly calculated from Equation (2):
(3)σr=σb,σθ=σH+σh−2σH−σhcos2θ−σb,σz=σV−2vσH−σhcos2θ,τrz=τrθ=τθz=0.Equation (3) reveals that the relationship between the three principal stress (σθ, σz, and σr) depends on not only the constants from in situ or laboratory test (e.g., aV, bV, aH, bH, ah, bh,ρb, and rb) but also the depth (z) and direction (θ) at the surrounding rocks. Even though the tangential principal stress (σθ) is the maximum principal stress, it is intractable to determine the relationship between the radial and axial principal stress (i.e., σr and σz, respectively). For instance, in shallow depth, the principal stress is likely to satisfy σθ>σr>σz. However, when the depth increases, the relationship of σθ>σz>σr can be satisfied, because σb will be near to its ultimate value (limzb→+∞σb=Cρbgrb). Therefore, to analyze the rock failure around the cylindrical caved space in varying depths, we employ the extended Mohr-Coulomb failure criterion for the problems in 3D space [27]
(4)σc+qσ3−σ1≥0,σc+qσ1−σ3≥0,σc+qσ2−σ1≥0,σc+qσ1−σ2≥0,σc+qσ3−σ2≥0,σc+qσ2−σ3≥0,where σ1, σ2, and σ3 are the maximum, intermediate, and minimum applied principal stress, respectively; σc is rock strength; σc can be replaced by long-term strength of rocks (σcd), if this criterion is utilized to analyze rock stability in a long duration of time [28]; the long-term strength of rocks (σcd) can be converted from the rock strength without time effect (σc′) by σcd=ασc′; α=0.4~0.6, for the intact rocks without preexisting fractures [29]; q=1+sinϕ/1−sinϕ; ϕ is the internal friction angle of rocks.Assuming the tangential, radial, and axial principal stress (i.e.,σθ, σr, and σz) are the maximum, intermediate, and minimum applied principal stress (i.e., σ1, σ2, and σ3), respectively, by substituting Equation (3) into Equation (4), noticing σ1=σθ, σ2=σr, and σ3=σz, we obtain the following relationships, which should be satisfied, if the surrounding rocks are expected to maintain stable in a long duration of time
(5)σcd+qσV−2vσH−σhcos2θ−σH+σh−2σH−σhcos2θ−σb≥0,σcd+qσH+σh−2σH−σhcos2θ−σb−σV−2vσH−σhcos2θ≥0,σcd+qσb−σH+σh−2σH−σhcos2θ−σb≥0,σcd+qσH+σh−2σH−σhcos2θ−σb−σb≥0,σcd+qσV−2vσH−σhcos2θ−σb≥0,σcd+qσb−σV−2vσH−σhcos2θ≥0.Because the tangential principal stress (σθ) is the maximum applied principal stress (i.e., either σθ>σr>σz or σθ>σz>σr is satisfied), Equation (5) can be reduced to the following equations:
(6a)σcd+qσV−2vσH−σhcos2θ−σH+σh−2σH−σhcos2θ−σb≥0,(6b)σcd+qσb−σH+σh−2σH−σhcos2θ−σb≥0.SubstitutingσH=aHz+bH, σh=ahz+bh, σV=aVz+bV, σcd=ασc′, and Equation (1) into Equations (6a) and (6b), we obtain the following equations to test the stability of the rocks around the cylindrical excavation:
(7a)σcd+Cρbgrb1−e−z−z0/4rb+qaV−aH+ah+21−qvaH−ahcos2θz+qbV−bH+bh+21−qvbH−bhcos2θ≥0,(7b)σcd+1+qCρbgrb1−e−z−z0/4rb−aH+ah−2aH−ahcos2θz+2bH−bhcos2θ−bH+bh≥0,where z is a variable that indicates the depth from ground surface and z0 is the depth of the free surface of caved rocks from ground surface.Equations (7a) and (7b) enable the prediction for the depth (z) where rock failure will take place around the cylindrical excavation (θ). If the depth of implemented undercut exceeds the results by Equations (7a) and (7b), rock failure and associated extension of surface subsidence can be expected to take place. On the other hand, Equations (7a) and (7b) also reveal the predicted depth for rock failure varies in different orientations. For instance, when σθ>σz>σr is satisfied (i.e., Equations (6b) and (7b) are valid to conduct such prediction), rock failure and associated extension of surface subsidence are most likely to take place in the direction of minimum in situ stress, i.e., θ=90° is satisfied. This means the extension of surface subsidence due to the rock shear failure can be expected to be anisotropic.On the other hand, Equations (6a) and (6b) reveal that the stability of surrounding rocks heavily depends on not only the rock strength (σcd) but also the stress due to caved rocks (σb). Equations (7a) and (7b) provide more details regarding the contribution of the caved rocks. Because of q>0 and C>0, the increase of the density (i.e., ρb increases) or height (i.e., z0 decrease) of caved rocks facilitates the stability of surrounding rocks. This provides some measures to prevent the rock failure and associated extension of surface subsidence: (1) filling the interspace between large-sized caved rocks by dumping small-sized waste rocks to increase the average density of the materials in the caved space or (2) backfilling the caved space with waste rocks to increase the height of caved rocks. Some discussions for such measures should be noted. Equation (1) shows the ultimate value of stress due to caved rocks (limzb→+∞σb=Cρbgrb) will increase, if measure (1) is implemented. This means that this measure contributes to the stability of surrounding rocks in all depth. However, the effect of measure (2) only works in shallow depth, because the ultimate value of the horizontal stress due to caved rocks maintains constant when the height of caved rocks changes.


## 3. Stability Analysis for Caved Space with Anisotropic Rocks due to Discontinuities

### 3.1. The Impact of Discontinuities on Rock Strength

Inherent discontinuities (e.g., faults, joints, and beddings) are the primary cause for rock strength weakening. Laboratory tests of the failure behavior in anisotropic rocks with polyaxial or triaxial compression demonstrate that the strength of anisotropic rocks varies significantly with the orientation of discontinuities [30–33].Figure3(a) shows such variation of the peak strength of anisotropic rocks, which can be obtained by the experimental method illustrated in Figure 3(b), at a constant confining pressure with the angle (i.e., β illustrated in Figure 3(b)) between the maximum applied stress (σ1) and the normal to the discontinuities. The slip failure along the discontinuities is most likely to take place, when β satisfies
(8)β=π4+ϕ′2,where β is the angle between the maximum applied stress and the normal to the discontinuities, as illustrated in Figure 3(b), and ϕ′ is the internal friction angle on the discontinuities.Figure 3
(a) Variation of the peak rock strength withβ, at a constant confining pressure. (b) Transversely isotropic specimen in triaxial compression. (c) 3D view of a specimen with inclined discontinuities in anisotropic stress state.
(a)(b)(c)To analyze this slip failure along the inherent discontinuities, Jaeger et al. [34] extended Coulomb’s failure criterion and proposed a new relationship for such analysis. The slip failure along discontinuities is expected to be prevented, if the following relationship is satisfied:
(9)σ1≤2c′+μ′σ31−μ′cotβsin2β+σ3,where c′ is the cohesion of discontinuities, μ′ is the coefficient of internal friction on discontinuities, and μ′=tanϕ′. When β tends to ϕ′ or 90°, the maximum applied principal stress (σ1), which accounts for the slip failure, tends to infinity.Jaeger’s model is valid to analyze the slip failure of surrounding rocks, when the discontinuities are parallel to the cylindrical excavation. This is because the analysis for the slip failure can be reduced to a 2D problem on a plane which is perpendicular to the axis of the excavation. The validity of Jaeger’s model has been demonstrated in a lot of underground engineering, such as borehole drilling [35] and reinforcement for underground tunnels [36]. However, when the discontinuities distribute inclining towards excavation’s axis, Jaeger’s model is no longer valid to conduct such analysis, because the stress accounting for the slip failure varies in the 3D space. For instance, as illustrated in Figure 3(c), with the confining stress σ3, the slip failure along discontinuities due to σ2 (βσ2−σ3) is likely to take place before the slip failure due to σ1 (βσ1−σ3), because of the anisotropy of rock strength in different orientations.Therefore, we extend Jaeger’s model to enable it to analyze the slip failure along inclined discontinuities with anisotropic applied principal stress (σ1>σ2>σ3) in 3D space, as illustrated in Figure 3(c). Assuming the redistributed stress is homogenous after vertical caving, an extended model after Jaeger’s is proposed for the slip failure along inclined discontinuities, based on the extended Mohr-Coulomb failure criterion in 3D space (i.e., Equation (4)):
(10)2c′+μ′σ31−μ′cotβσ1−σ3sin2βσ1−σ3+σ3−σ1≥0,2c′+μ′σ11−μ′cotβσ3−σ1sin2βσ3−σ1+σ1−σ3≥0,2c′+μ′σ21−μ′cotβσ1−σ2sin2βσ1−σ2+σ2−σ1≥0,2c′+μ′σ11−μ′cotβσ2−σ1sin2βσ2−σ1+σ1−σ2≥0,2c′+μ′σ31−μ′cotβσ2−σ3sin2βσ2−σ3+σ3−σ2≥0,2c′+μ′σ21−μ′cotβσ3−σ2sin2βσ3−σ2+σ2−σ3≥0,where βσi−σj is angle between the stress that accounts for the slip failure (σi) and the normal to the discontinuities on the plane consisted by confining stress (σj) and the stress accounts for the slip failure (σi).


### 3.2. Slip Failure along Discontinuities at Caved Space Wall

The extended Jaeger model (i.e., Equation (10)) is employed to test the stability of the anisotropic rocks around the cylindrical caved space. Assuming the tangential, radial, and axial principal stress (i.e., σθ, σr, and σz) are the maximum, intermediate, and minimum applied principal stress (i.e., σ1, σ2, and σ3), respectively, we substitute Equation (3) into Equation (9), noticing σ1=σθ, σ2=σr, and σ3=σz; we obtain the following relationships, which should be satisfied, if the surrounding rocks are expected to keep stable from the slip failure along discontinuities:
(11)2c′+2μ′σV−2vσH−σhcos2θ1−μ′cotβσθ−σzsin2βσθ−σz+σV−2vσH−σhcos2θ−σH+σh−2σH−σhcos2θ−σb≥0,2c′+2μ′σH+σh−2σH−σhcos2θ−σb1−μ′cotβσz−σθsin2βσz−σθ+σH+σh−2σH−σhcos2θ−σb−σV−2vσH−σhcos2θ≥0,2c′+2μ′σb1−μ′cotβσθ−σrsin2βσθ−σr+σb−σH+σh−2σH−σhcos2θ−σb≥0,2c′+2μ′σH+σh−2σH−σhcos2θ−σb1−μ′cotβσr−σθsin2βσr−σθ+σH+σh−2σH−σhcos2θ−σb−σb≥0,2c′+2μ′σV−2vσH−σhcos2θ1−μ′cotβσr−σzsin2βσr−σz+σV−2vσH−σhcos2θ−σb≥0,2c′+2μ′σb1−μ′cotβσz−σrsin2βσz−σr+σb−σV−2vσH−σhcos2θ≥0.Because the tangential principal stress (σθ) is the maximum applied principal stress (i.e., either σθ>σr>σz or σθ>σz>σr is satisfied), Equation (10) can be reduced to the following equations:
(12)2c′+2μ′σV−2vσH−σhcos2θ1−μ′cotβσθ−σzsin2βσθ−σz+σV−2vσH−σhcos2θ−σH+σh−2σH−σhcos2θ−σb≥0,2c′+2μ′σb1−μ′cotβσθ−σrsin2βσθ−σr+σb−σH+σh−2σH−σhcos2θ−σb≥0,2c′+2μ′σV−2vσH−σhcos2θ1−μ′cotβσr−σzsin2βσr−σz+σV−2vσH−σhcos2θ−σb≥0,2c′+2μ′σb1−μ′cotβσz−σrsin2βσz−σr+σb−σV−2vσH−σhcos2θ≥0.SubstitutingσH=aHz+bH, σh=ahz+bh, σV=aVz+bV, and Equation (1) into Equation (12), we obtain the following equations to test whether the slip failure along discontinuities will take place around the cylindrical caved space:
(13a)2c′+1−μ′cotβσθ−σzsin2βσθ−σzCρbgr1−e−z−z0/4rb+2μ′+1−μ′cotβσθ−σzsin2βσθ−σzaVz+bV+21−v1−μ′cotβσθ−σzsin2βσθ−σz−4vμ′aH−ahz+bH−bhcos2θ−1−μ′cotβσθ−σzsin2βσθ−σzaH+ahz+bH+bh≥0,(13b)2c′+2μ′+1−μ′cotβσθ−σrsin2βσθ−σrCρbgr1−e−z−z0/4rb−1−μ′cotβσθ−σrsin2βσθ−σraH+ahz−2aH−ahzcos2θ+bH+bh−2bH−bhcos2θ≥0,(13c)2c′−1−μ′cotβσr−σzsin2βσr−σzCρbgrb1−e−zb/4rb+2μ′+1−μ′cotβσr−σzsin2βσr−σzaVz+bV−2vaH−ahzcos2θ−2vbH−bhcos2θ≥0,(13d)2c′+2μ′+1−μ′cotβσz−σrsin2βσz−σrCρbgrb1−e−zb/4rb−1−μ′cotβσz−σrsin2βσz−σraVz+bV−2vaH−ahzcos2θ−2vbH−bhcos2θ≥0.Equation (12) enables the prediction for orientation (θ) and depth (z) where slip failure along the inclined discontinuities will take place, and it shows the occurrence of slip failure heavily depends on the property (c′ and μ′) of discontinuities. This means the slip failure along discontinuities and associated extension of surface subsidence may suddenly occur due to the dramatic decrease of discontinuities’ cohesion, such as when the infill in discontinuities is saturated because of heavy precipitation or underground water [37]. Additionally, laboratory tests have demonstrated the strength of discontinuities will reduce after the slip failure [38]. This reduction will deteriorate the stability of surrounding rocks, because the initial slip failure in one dimension (e.g., when βσθ−σr=π/4+ϕ′/2 is satisfied) is likely to lead to the subsequent failure in other directions (e.g., βσθ−σz=π/4+ϕ′/2 or βσz−σr=π/4+ϕ′/2 is satisfied).On the other hand, Equations (13a)–(13d) show the impact of caved materials on the stability of surrounding rocks. Because of ϕ′<βσi−σj<90°, the caved rocks contribute to keeping the surrounding rocks from slip failure (i.e., Equations (13a), (13b), and (13d)), apart from the condition that the radial principal stress (σb) accounts for the slip failure (i.e., Equation (13c)). This means the measures to prevent the shear failure of isotropic surrounding rocks (i.e., (1) filling the interspace between large-sized caved rocks by dumping small-sized waste rocks or (2) backfilling the caved space with waste rocks) are also valid to enhance the stability of anisotropic surrounding rocks, as well as to mitigate the extension of surface subsidence, if the rock failure is due to tangential principal stress (σθ) or axial principal stress (σz).


## 3.1. The Impact of Discontinuities on Rock Strength

Inherent discontinuities (e.g., faults, joints, and beddings) are the primary cause for rock strength weakening. Laboratory tests of the failure behavior in anisotropic rocks with polyaxial or triaxial compression demonstrate that the strength of anisotropic rocks varies significantly with the orientation of discontinuities [30–33].Figure3(a) shows such variation of the peak strength of anisotropic rocks, which can be obtained by the experimental method illustrated in Figure 3(b), at a constant confining pressure with the angle (i.e., β illustrated in Figure 3(b)) between the maximum applied stress (σ1) and the normal to the discontinuities. The slip failure along the discontinuities is most likely to take place, when β satisfies
(8)β=π4+ϕ′2,where β is the angle between the maximum applied stress and the normal to the discontinuities, as illustrated in Figure 3(b), and ϕ′ is the internal friction angle on the discontinuities.Figure 3
(a) Variation of the peak rock strength withβ, at a constant confining pressure. (b) Transversely isotropic specimen in triaxial compression. (c) 3D view of a specimen with inclined discontinuities in anisotropic stress state.
(a)(b)(c)To analyze this slip failure along the inherent discontinuities, Jaeger et al. [34] extended Coulomb’s failure criterion and proposed a new relationship for such analysis. The slip failure along discontinuities is expected to be prevented, if the following relationship is satisfied:
(9)σ1≤2c′+μ′σ31−μ′cotβsin2β+σ3,where c′ is the cohesion of discontinuities, μ′ is the coefficient of internal friction on discontinuities, and μ′=tanϕ′. When β tends to ϕ′ or 90°, the maximum applied principal stress (σ1), which accounts for the slip failure, tends to infinity.Jaeger’s model is valid to analyze the slip failure of surrounding rocks, when the discontinuities are parallel to the cylindrical excavation. This is because the analysis for the slip failure can be reduced to a 2D problem on a plane which is perpendicular to the axis of the excavation. The validity of Jaeger’s model has been demonstrated in a lot of underground engineering, such as borehole drilling [35] and reinforcement for underground tunnels [36]. However, when the discontinuities distribute inclining towards excavation’s axis, Jaeger’s model is no longer valid to conduct such analysis, because the stress accounting for the slip failure varies in the 3D space. For instance, as illustrated in Figure 3(c), with the confining stress σ3, the slip failure along discontinuities due to σ2 (βσ2−σ3) is likely to take place before the slip failure due to σ1 (βσ1−σ3), because of the anisotropy of rock strength in different orientations.Therefore, we extend Jaeger’s model to enable it to analyze the slip failure along inclined discontinuities with anisotropic applied principal stress (σ1>σ2>σ3) in 3D space, as illustrated in Figure 3(c). Assuming the redistributed stress is homogenous after vertical caving, an extended model after Jaeger’s is proposed for the slip failure along inclined discontinuities, based on the extended Mohr-Coulomb failure criterion in 3D space (i.e., Equation (4)):
(10)2c′+μ′σ31−μ′cotβσ1−σ3sin2βσ1−σ3+σ3−σ1≥0,2c′+μ′σ11−μ′cotβσ3−σ1sin2βσ3−σ1+σ1−σ3≥0,2c′+μ′σ21−μ′cotβσ1−σ2sin2βσ1−σ2+σ2−σ1≥0,2c′+μ′σ11−μ′cotβσ2−σ1sin2βσ2−σ1+σ1−σ2≥0,2c′+μ′σ31−μ′cotβσ2−σ3sin2βσ2−σ3+σ3−σ2≥0,2c′+μ′σ21−μ′cotβσ3−σ2sin2βσ3−σ2+σ2−σ3≥0,where βσi−σj is angle between the stress that accounts for the slip failure (σi) and the normal to the discontinuities on the plane consisted by confining stress (σj) and the stress accounts for the slip failure (σi).


## 3.2. Slip Failure along Discontinuities at Caved Space Wall

The extended Jaeger model (i.e., Equation (10)) is employed to test the stability of the anisotropic rocks around the cylindrical caved space. Assuming the tangential, radial, and axial principal stress (i.e., σθ, σr, and σz) are the maximum, intermediate, and minimum applied principal stress (i.e., σ1, σ2, and σ3), respectively, we substitute Equation (3) into Equation (9), noticing σ1=σθ, σ2=σr, and σ3=σz; we obtain the following relationships, which should be satisfied, if the surrounding rocks are expected to keep stable from the slip failure along discontinuities:
(11)2c′+2μ′σV−2vσH−σhcos2θ1−μ′cotβσθ−σzsin2βσθ−σz+σV−2vσH−σhcos2θ−σH+σh−2σH−σhcos2θ−σb≥0,2c′+2μ′σH+σh−2σH−σhcos2θ−σb1−μ′cotβσz−σθsin2βσz−σθ+σH+σh−2σH−σhcos2θ−σb−σV−2vσH−σhcos2θ≥0,2c′+2μ′σb1−μ′cotβσθ−σrsin2βσθ−σr+σb−σH+σh−2σH−σhcos2θ−σb≥0,2c′+2μ′σH+σh−2σH−σhcos2θ−σb1−μ′cotβσr−σθsin2βσr−σθ+σH+σh−2σH−σhcos2θ−σb−σb≥0,2c′+2μ′σV−2vσH−σhcos2θ1−μ′cotβσr−σzsin2βσr−σz+σV−2vσH−σhcos2θ−σb≥0,2c′+2μ′σb1−μ′cotβσz−σrsin2βσz−σr+σb−σV−2vσH−σhcos2θ≥0.Because the tangential principal stress (σθ) is the maximum applied principal stress (i.e., either σθ>σr>σz or σθ>σz>σr is satisfied), Equation (10) can be reduced to the following equations:
(12)2c′+2μ′σV−2vσH−σhcos2θ1−μ′cotβσθ−σzsin2βσθ−σz+σV−2vσH−σhcos2θ−σH+σh−2σH−σhcos2θ−σb≥0,2c′+2μ′σb1−μ′cotβσθ−σrsin2βσθ−σr+σb−σH+σh−2σH−σhcos2θ−σb≥0,2c′+2μ′σV−2vσH−σhcos2θ1−μ′cotβσr−σzsin2βσr−σz+σV−2vσH−σhcos2θ−σb≥0,2c′+2μ′σb1−μ′cotβσz−σrsin2βσz−σr+σb−σV−2vσH−σhcos2θ≥0.SubstitutingσH=aHz+bH, σh=ahz+bh, σV=aVz+bV, and Equation (1) into Equation (12), we obtain the following equations to test whether the slip failure along discontinuities will take place around the cylindrical caved space:
(13a)2c′+1−μ′cotβσθ−σzsin2βσθ−σzCρbgr1−e−z−z0/4rb+2μ′+1−μ′cotβσθ−σzsin2βσθ−σzaVz+bV+21−v1−μ′cotβσθ−σzsin2βσθ−σz−4vμ′aH−ahz+bH−bhcos2θ−1−μ′cotβσθ−σzsin2βσθ−σzaH+ahz+bH+bh≥0,(13b)2c′+2μ′+1−μ′cotβσθ−σrsin2βσθ−σrCρbgr1−e−z−z0/4rb−1−μ′cotβσθ−σrsin2βσθ−σraH+ahz−2aH−ahzcos2θ+bH+bh−2bH−bhcos2θ≥0,(13c)2c′−1−μ′cotβσr−σzsin2βσr−σzCρbgrb1−e−zb/4rb+2μ′+1−μ′cotβσr−σzsin2βσr−σzaVz+bV−2vaH−ahzcos2θ−2vbH−bhcos2θ≥0,(13d)2c′+2μ′+1−μ′cotβσz−σrsin2βσz−σrCρbgrb1−e−zb/4rb−1−μ′cotβσz−σrsin2βσz−σraVz+bV−2vaH−ahzcos2θ−2vbH−bhcos2θ≥0.Equation (12) enables the prediction for orientation (θ) and depth (z) where slip failure along the inclined discontinuities will take place, and it shows the occurrence of slip failure heavily depends on the property (c′ and μ′) of discontinuities. This means the slip failure along discontinuities and associated extension of surface subsidence may suddenly occur due to the dramatic decrease of discontinuities’ cohesion, such as when the infill in discontinuities is saturated because of heavy precipitation or underground water [37]. Additionally, laboratory tests have demonstrated the strength of discontinuities will reduce after the slip failure [38]. This reduction will deteriorate the stability of surrounding rocks, because the initial slip failure in one dimension (e.g., when βσθ−σr=π/4+ϕ′/2 is satisfied) is likely to lead to the subsequent failure in other directions (e.g., βσθ−σz=π/4+ϕ′/2 or βσz−σr=π/4+ϕ′/2 is satisfied).On the other hand, Equations (13a)–(13d) show the impact of caved materials on the stability of surrounding rocks. Because of ϕ′<βσi−σj<90°, the caved rocks contribute to keeping the surrounding rocks from slip failure (i.e., Equations (13a), (13b), and (13d)), apart from the condition that the radial principal stress (σb) accounts for the slip failure (i.e., Equation (13c)). This means the measures to prevent the shear failure of isotropic surrounding rocks (i.e., (1) filling the interspace between large-sized caved rocks by dumping small-sized waste rocks or (2) backfilling the caved space with waste rocks) are also valid to enhance the stability of anisotropic surrounding rocks, as well as to mitigate the extension of surface subsidence, if the rock failure is due to tangential principal stress (σθ) or axial principal stress (σz).


## 4. Case Study for Extension of Surface Subsidence in Xiaowanggou Iron Mine

Herein is presented a case study for Xiaowanggou iron mine. Some implications from the proposed model and case study are discussed.


### 4.1. Prediction for Subsequent Extension of Surface Subsidence after Vertical Caving

Xiaowanggou iron mine is located in Liaoyang, Liaoning of China, whose deposit is extracted by sublevel caving method. The surface subsidence appears circular after the vertical caving. To predict the extension of surface subsidence, the proposed model is employed. The inputs should be noted firstly, such as in situ stress, rock strength, property of caved materials, and the distribution of discontinuities in 3D space.We conducted the point load test for the strength of surrounding rocks, and it is converted to the long-term strength to analyze the stability of caved space in a long duration of time [22]. On the other hand, we also collected the caved rocks in situ, and scaled laboratory tests are thus conducted for the constant (i.e., C) in Equation (1). By regressing 3 sets of laboratory test data, the results demonstrate C=0.991 (R2=0.993), as illustrated in Figure 4.Figure 4
The variation ofσb with the depth from the free surface of caved rocks by laboratory tests.The in situ stress is obtained by regressing 12 sets of data by anelastic strain recovery and hydraulic fracturing methods, which is provided by the Fundamental Database of Crustal Stress Environment in continental China. The direction of the maximum in situ stress is N80°E [39]. The in situ stress in Xiaowanggou iron mine can be expressed as follows:
(14)σH=0.0304z+2.9033,σV=0.0244z+1.1593,σh=0.0149z+2.6795.The in situ investigation shows the strike and inclination of the discontinuities are N15°E and 76°, respectively. Figure 5 shows the distribution of caved space, in situ stress, and discontinuities in 3D space.Figure 5
(a) The distribution of caved space, in situ stress, and discontinuities in 3D space. (b) In situ observation for the discontinuities in the undercut in Xiaowanggou iron mine. (c) The orientation of mean maximum horizontal in situ stress (represented by black bar) around Xiaowanggou iron mine.Table1 lists the inputs to calculate the critical depth of undercut, above which either rock shear failure or slip failure along discontinuities can be prevented, if the implemented undercut locates.Table 1
Parameters of the surrounding intact rocks, caved rocks, and discontinuities.
Surrounding intact rocksDiscontinuitiesCaved rocksIs50 (MPa)σcd (MPa)ϕ (°)vc′ (MPa)ϕ′ (°)ρB (t/m3)CrB (m)z0 (m)4.09332.744350.253.25201.770.99177.0545The in situ stress and involved parameters reveal the three principal stress (σθ, σz, σr) satisfies σθ>σz>σr under the free surface of caved rocks (z0=45m). This means Equations (13a), (13b), (13d), and (7a) are valid to predict the slip failure along discontinuities and the shear failure occurring at the surrounding rock of caved space, respectively. Additionally, the calculation for the angles related to slip failure, i.e., βσθ−σr, βσθ−σz, and βσz−σr, should be noted. The angle δ is defined to describe the distribution of in situ stress and the strike of the discontinuities on the plane constituted by the maximum and minimum in situ stress (Figure 5(a)), and the angles related to the slip failure in Xiaowanggou iron mine can be calculated by the following equations:
(15)βσθ−σr=θ−δ,βσθ−σz=arctan1tanη·cosθ,βσz−σr=arctantanη·sinθ,where δ is the angle between discontinuities’ strike and the maximum in situ stress (σH) on the plane constituted by the maximum and minimum in situ stress and η is the inclination of the discontinuities.Therefore, we calculate the critical depth of undercut around the caved space, above which the surrounding wall maintains stable from either rock shear failure or slip failure along discontinuities. Because of symmetry, we only present the results in a half circumference, i.e.,θ varies from 0° to 180°, as shown in Figure 6.Figure 6
The variation of the critical depth of undercut around the cylindrical caved space, above which the surrounding rocks maintain stable from either rock shear failure (i.e., “shear failure”), or the slip failure along discontinuities (i.e., “slip failure due toσθ” or “slip failure due to σz”). The undercut locates at the depth of 168 m (i.e., “undercut depth”), and the slip failure of surrounding rocks along discontinuities will take place at the direction between N18°W and N55°W, as well as between S18°E and S55°E due to symmetry, because the critical depth of undercut to prevent slip failure due to σθ (i.e., magenta line) locates above the depth of implemented undercut in these ranges of direction.As expected, the variation of the critical depth related to rock shear failure (i.e., “shear failure” in Figure6) is symmetric on either sides of the direction where the maximum effective principal stress (σθ) has maximum value (i.e., θ=90°, N10°W, or the direction of minimum in situ stress). The minimum value of such critical depth appears at z=406m, in the direction of minimum in situ stress. However, the depth of implemented undercut locates at z=168m, and this means rock shear failure is not responsible for extension of surface subsidence in Xiaowanggou iron mine.On the other hand, Figure6 provides the results of the critical depth of undercut related to the slip failure along discontinuities. The results indicate the slip failure due to σθ, with the confining stress σz, will not take place above the depth z=2500m in this project. When the σr serves as the confining stress, the critical depth related to the slip failure due to σθ (i.e., “slip failure due to σθ” in Figure 6) or σz (i.e., “slip failure due to σz” in Figure 6) presents an asymmetrical distribution. Their ultimate value of the critical depth appears near the direction where Equation (8) is satisfied (i.e., βσθ−σr=55° or βσz−σr=55°), but the offset towards the direction of minimum in situ stress can be observed. The variation of σθ and σz around the caved space is responsible to this offset, as well as the difference of the ultimate value of critical undercut depth in each distribution interval (e.g., the ultimate value of critical undercut depth due to σθ between θ=0° and θ=40° is smaller than the one between θ=90° and θ=150°). The minimum critical depth of undercut appear at z=55m when θ=115°, i.e., N35°W, and this value is smaller than the depth of implemented undercut. This means the slip failure along discontinuities and the associated extension of surface subsidence are likely to be firstly observed in N35°W, as well as in S35°E due to symmetry. Therefore, the extension of surface subsidence in Xiaowanggou iron mine is predicted. In the depth of the free surface of caved rocks (i.e., z=45m or “caved rocks surface” in Figure 6), the surrounding rocks of the caved space maintain stable from either rock shear failure or slip failure along discontinuities. In the depth of implemented undercut (i.e., z=168m or “undercut depth” in Figure 6), Figure 6 shows the slip failure due to σθ will take place between θ=98° and θ=135°, i.e., between N18°W and N55°W, as well as between S18°E and S55°E due to symmetry. This means the surface subsidence is likely to extend to such orientations. To test the validity of this prediction, we compared it with the satellite images by Google Earth, as illustrated in Figure 7.Figure 7
The comparison between the prediction by the proposed model and the observation in Google Earth satellite images for the surface subsidence extension in Xiaowanggou iron mine (41.39°N, 123.63°E). (a) The circular surface subsidence observed after vertical caving on October 1, 2012. (b) The extended surface subsidence area observed on April 25, 2014. (c) The extended surface subsidence area on May 8, 2017. The blue shadow illustrated indicates direction that the surface subsidence extends to, which is obtained from Figure 6, by comparing the depth of implemented undercut to critical depth of undercut to maintain the surrounding rocks from either rock shear failure or slip failure along the discontinuities.The extension process of surface subsidence is illustrated in Figure7, and the blue shadow indicates the predicted orientation that the surface subsidence extends to by the proposed model. Figure 7 shows the observed extension of surface subsidence generally matches the prediction, even though the observed range of orientation that surface subsidence extends to is larger than the predicted. Such results demonstrate the validity of the proposed model for the prediction of subsequent extension of surface subsidence after vertical caving.


### 4.2. Implications from the Proposed Model and Illustrative Example

Herein is presented some implications regarding the proposed model and illustrative example. The difference of orientation range, where surface subsidence extends to, between the observation in satellite image and the prediction by proposed model should be discussed first. The long-term strength is involved to test the stability of surrounding isotropic rock in a long duration of time because of the consideration of time effect, such as rock creep or the saturation by surface or underground water. However, rare analytical solutions regarding long-term strength of discontinuities have been reported, even though the impact of time on discontinuities’ strength has been commonly observed [40, 41]. Thus, the time effect is not involved in the proposed model for the slip failure (Equations (11), (12), and (13a)–(13d)), and this leads to observed orientation range where surface subsidence extension that takes place is larger than the predicted. Meanwhile, the absence of the consideration for the impact of time on discontinuities’ strength is a limitation of the proposed model.The comparison for surface subsidence extension between the observation in satellite images and the analytical prediction demonstrates the validity of the proposed model. Compared with the existing models [13–16], the proposed targets predict the subsequent extension of surface subsidence after vertical caving, which is rarely analyzed in previous literatures. Both rock shear failure and slip failure along discontinuities are involved in this model, and the case study shows such two failures are likely to take place either in the direction of minimum in situ stress or near the direction where the anisotropic rocks have minimum strength, respectively. Moreover, the impact of caved rocks on the stability of surrounding rocks of caved space is addressed, and some approaches to prevent the extension of surface subsidence are thus suggested, e.g., filling the interspace between large-sized caved rocks by dumping small-sized waste rocks or backfilling the caved space with waste rocks.


## 4.1. Prediction for Subsequent Extension of Surface Subsidence after Vertical Caving

Xiaowanggou iron mine is located in Liaoyang, Liaoning of China, whose deposit is extracted by sublevel caving method. The surface subsidence appears circular after the vertical caving. To predict the extension of surface subsidence, the proposed model is employed. The inputs should be noted firstly, such as in situ stress, rock strength, property of caved materials, and the distribution of discontinuities in 3D space.We conducted the point load test for the strength of surrounding rocks, and it is converted to the long-term strength to analyze the stability of caved space in a long duration of time [22]. On the other hand, we also collected the caved rocks in situ, and scaled laboratory tests are thus conducted for the constant (i.e., C) in Equation (1). By regressing 3 sets of laboratory test data, the results demonstrate C=0.991 (R2=0.993), as illustrated in Figure 4.Figure 4
The variation ofσb with the depth from the free surface of caved rocks by laboratory tests.The in situ stress is obtained by regressing 12 sets of data by anelastic strain recovery and hydraulic fracturing methods, which is provided by the Fundamental Database of Crustal Stress Environment in continental China. The direction of the maximum in situ stress is N80°E [39]. The in situ stress in Xiaowanggou iron mine can be expressed as follows:
(14)σH=0.0304z+2.9033,σV=0.0244z+1.1593,σh=0.0149z+2.6795.The in situ investigation shows the strike and inclination of the discontinuities are N15°E and 76°, respectively. Figure 5 shows the distribution of caved space, in situ stress, and discontinuities in 3D space.Figure 5
(a) The distribution of caved space, in situ stress, and discontinuities in 3D space. (b) In situ observation for the discontinuities in the undercut in Xiaowanggou iron mine. (c) The orientation of mean maximum horizontal in situ stress (represented by black bar) around Xiaowanggou iron mine.Table1 lists the inputs to calculate the critical depth of undercut, above which either rock shear failure or slip failure along discontinuities can be prevented, if the implemented undercut locates.Table 1
Parameters of the surrounding intact rocks, caved rocks, and discontinuities.
Surrounding intact rocksDiscontinuitiesCaved rocksIs50 (MPa)σcd (MPa)ϕ (°)vc′ (MPa)ϕ′ (°)ρB (t/m3)CrB (m)z0 (m)4.09332.744350.253.25201.770.99177.0545The in situ stress and involved parameters reveal the three principal stress (σθ, σz, σr) satisfies σθ>σz>σr under the free surface of caved rocks (z0=45m). This means Equations (13a), (13b), (13d), and (7a) are valid to predict the slip failure along discontinuities and the shear failure occurring at the surrounding rock of caved space, respectively. Additionally, the calculation for the angles related to slip failure, i.e., βσθ−σr, βσθ−σz, and βσz−σr, should be noted. The angle δ is defined to describe the distribution of in situ stress and the strike of the discontinuities on the plane constituted by the maximum and minimum in situ stress (Figure 5(a)), and the angles related to the slip failure in Xiaowanggou iron mine can be calculated by the following equations:
(15)βσθ−σr=θ−δ,βσθ−σz=arctan1tanη·cosθ,βσz−σr=arctantanη·sinθ,where δ is the angle between discontinuities’ strike and the maximum in situ stress (σH) on the plane constituted by the maximum and minimum in situ stress and η is the inclination of the discontinuities.Therefore, we calculate the critical depth of undercut around the caved space, above which the surrounding wall maintains stable from either rock shear failure or slip failure along discontinuities. Because of symmetry, we only present the results in a half circumference, i.e.,θ varies from 0° to 180°, as shown in Figure 6.Figure 6
The variation of the critical depth of undercut around the cylindrical caved space, above which the surrounding rocks maintain stable from either rock shear failure (i.e., “shear failure”), or the slip failure along discontinuities (i.e., “slip failure due toσθ” or “slip failure due to σz”). The undercut locates at the depth of 168 m (i.e., “undercut depth”), and the slip failure of surrounding rocks along discontinuities will take place at the direction between N18°W and N55°W, as well as between S18°E and S55°E due to symmetry, because the critical depth of undercut to prevent slip failure due to σθ (i.e., magenta line) locates above the depth of implemented undercut in these ranges of direction.As expected, the variation of the critical depth related to rock shear failure (i.e., “shear failure” in Figure6) is symmetric on either sides of the direction where the maximum effective principal stress (σθ) has maximum value (i.e., θ=90°, N10°W, or the direction of minimum in situ stress). The minimum value of such critical depth appears at z=406m, in the direction of minimum in situ stress. However, the depth of implemented undercut locates at z=168m, and this means rock shear failure is not responsible for extension of surface subsidence in Xiaowanggou iron mine.On the other hand, Figure6 provides the results of the critical depth of undercut related to the slip failure along discontinuities. The results indicate the slip failure due to σθ, with the confining stress σz, will not take place above the depth z=2500m in this project. When the σr serves as the confining stress, the critical depth related to the slip failure due to σθ (i.e., “slip failure due to σθ” in Figure 6) or σz (i.e., “slip failure due to σz” in Figure 6) presents an asymmetrical distribution. Their ultimate value of the critical depth appears near the direction where Equation (8) is satisfied (i.e., βσθ−σr=55° or βσz−σr=55°), but the offset towards the direction of minimum in situ stress can be observed. The variation of σθ and σz around the caved space is responsible to this offset, as well as the difference of the ultimate value of critical undercut depth in each distribution interval (e.g., the ultimate value of critical undercut depth due to σθ between θ=0° and θ=40° is smaller than the one between θ=90° and θ=150°). The minimum critical depth of undercut appear at z=55m when θ=115°, i.e., N35°W, and this value is smaller than the depth of implemented undercut. This means the slip failure along discontinuities and the associated extension of surface subsidence are likely to be firstly observed in N35°W, as well as in S35°E due to symmetry. Therefore, the extension of surface subsidence in Xiaowanggou iron mine is predicted. In the depth of the free surface of caved rocks (i.e., z=45m or “caved rocks surface” in Figure 6), the surrounding rocks of the caved space maintain stable from either rock shear failure or slip failure along discontinuities. In the depth of implemented undercut (i.e., z=168m or “undercut depth” in Figure 6), Figure 6 shows the slip failure due to σθ will take place between θ=98° and θ=135°, i.e., between N18°W and N55°W, as well as between S18°E and S55°E due to symmetry. This means the surface subsidence is likely to extend to such orientations. To test the validity of this prediction, we compared it with the satellite images by Google Earth, as illustrated in Figure 7.Figure 7
The comparison between the prediction by the proposed model and the observation in Google Earth satellite images for the surface subsidence extension in Xiaowanggou iron mine (41.39°N, 123.63°E). (a) The circular surface subsidence observed after vertical caving on October 1, 2012. (b) The extended surface subsidence area observed on April 25, 2014. (c) The extended surface subsidence area on May 8, 2017. The blue shadow illustrated indicates direction that the surface subsidence extends to, which is obtained from Figure 6, by comparing the depth of implemented undercut to critical depth of undercut to maintain the surrounding rocks from either rock shear failure or slip failure along the discontinuities.The extension process of surface subsidence is illustrated in Figure7, and the blue shadow indicates the predicted orientation that the surface subsidence extends to by the proposed model. Figure 7 shows the observed extension of surface subsidence generally matches the prediction, even though the observed range of orientation that surface subsidence extends to is larger than the predicted. Such results demonstrate the validity of the proposed model for the prediction of subsequent extension of surface subsidence after vertical caving.


## 4.2. Implications from the Proposed Model and Illustrative Example

Herein is presented some implications regarding the proposed model and illustrative example. The difference of orientation range, where surface subsidence extends to, between the observation in satellite image and the prediction by proposed model should be discussed first. The long-term strength is involved to test the stability of surrounding isotropic rock in a long duration of time because of the consideration of time effect, such as rock creep or the saturation by surface or underground water. However, rare analytical solutions regarding long-term strength of discontinuities have been reported, even though the impact of time on discontinuities’ strength has been commonly observed [40, 41]. Thus, the time effect is not involved in the proposed model for the slip failure (Equations (11), (12), and (13a)–(13d)), and this leads to observed orientation range where surface subsidence extension that takes place is larger than the predicted. Meanwhile, the absence of the consideration for the impact of time on discontinuities’ strength is a limitation of the proposed model.The comparison for surface subsidence extension between the observation in satellite images and the analytical prediction demonstrates the validity of the proposed model. Compared with the existing models [13–16], the proposed targets predict the subsequent extension of surface subsidence after vertical caving, which is rarely analyzed in previous literatures. Both rock shear failure and slip failure along discontinuities are involved in this model, and the case study shows such two failures are likely to take place either in the direction of minimum in situ stress or near the direction where the anisotropic rocks have minimum strength, respectively. Moreover, the impact of caved rocks on the stability of surrounding rocks of caved space is addressed, and some approaches to prevent the extension of surface subsidence are thus suggested, e.g., filling the interspace between large-sized caved rocks by dumping small-sized waste rocks or backfilling the caved space with waste rocks.


## 5. Conclusion

Surface subsidence due to underground mining leads to severe geological hazards. The extension of surface subsidence after vertical caving is a primary contributor to the large-scale ground destruction, as well as the associated geological hazards. However, it is still intractable to analytically predict the extension of circular surface subsidence, because of the absence of robust model.To fill such a research gap, employing in situ stress, property of surrounding intact rock, inherent discontinuities, and caved rocks, an analytical model is proposed to calculate the depth around the cylindrical caved space, above which the surrounding rocks maintain stable from either rock shear failure or the slip failure along discontinuities. The comparison between the predicted extension of circular surface subsidence in Xiaowanggou iron mine and the in situ observation demonstrates the validity of the proposed model. The proposed model reveals the surface subsidence tends to extend to the direction of minimum in situ stress, when the rock shear failure accounts for such extension. Meanwhile, the extension of surface subsidence due to slip failure along discontinuities is most likely to take place near the direction where the anisotropic rocks have minimum strength. On the other hand, the impact of caved rocks clarified by the proposed model shows the increase of either density or height of caved rocks facilitates the stability of the rocks around the caved space. Some measures are accordingly suggested to prevent the extension of surface subsidence, including filling the interspace between large-sized caved rocks by dumping small-sized waste rocks or backfilling the caved space with waste rocks.Finally, the contribution of this study should be addressed. Theoretically, this paper proposed an analytical solution to predict the subsequent surface subsidence after vertical caving, by analyzing either rock shear failure or slip failure along discontinuities that occurred at surrounding rock. Practically, the validity of this solution has been validated, which means it has the potential to guide the underground metal mines to mitigate the impact of surface subsidence on safety and environment.



---

*Source: 4309530-2021-05-27.xml*