#### 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.

#### 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.

**(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. Section 2 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 (Figure 1), is assumed as a vertical and cylindrical excavation filled with caved rocks (Figure 2(a)).

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]: where is the horizontal stress from caved rocks to the surrounding rocks, is the average density of caved rocks, is the gravity of earth, is the radius of the cross section of caved space, is a variable representing the depth measured from the free surface of caved rocks, and 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 (Figure 2(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 (, , ), original due to Krisch: where , , and are the radial, tangential, and axial local normal stress near the rocks around the cylindrical excavation, respectively, , , and are the radial, tangential, and axial local shear stress, respectively, and , , and are the vertical, maximum, and minimum horizontal in situ stress, respectively; the in situ stress linearly varies with the depth from ground surface () and can be expressed as , , and ; , , , , , and 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); 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., , as illustrated in Figure 2(b)); represents the direction of maximum horizontal in situ stress (), and represents the direction of minimum horizontal in situ stress (); is Poisson’s ratio of the surround rocks.

The normal and shear stress at caved space wall () can be accordingly calculated from Equation (2):

Equation (3) reveals that the relationship between the three principal stress (, , and ) depends on not only the constants from in situ or laboratory test (e.g., , , , , , ,, and ) but also the depth () 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., and , respectively). For instance, in shallow depth, the principal stress is likely to satisfy . However, when the depth increases, the relationship of can be satisfied, because will be near to its ultimate value (). 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] where , , and are the maximum, intermediate, and minimum applied principal stress, respectively; is rock strength; can be replaced by long-term strength of rocks (), if this criterion is utilized to analyze rock stability in a long duration of time [28]; the long-term strength of rocks () can be converted from the rock strength without time effect () by ; , for the intact rocks without preexisting fractures [29]; ; is the internal friction angle of rocks.

Assuming the tangential, radial, and axial principal stress (i.e., , , and ) are the maximum, intermediate, and minimum applied principal stress (i.e., , , and ), respectively, by substituting Equation (3) into Equation (4), noticing , , and , we obtain the following relationships, which should be satisfied, if the surrounding rocks are expected to maintain stable in a long duration of time

Because the tangential principal stress () is the maximum applied principal stress (i.e., either or is satisfied), Equation (5) can be reduced to the following equations:

Substituting , , , , and Equation (1) into Equations (6a) and (6b), we obtain the following equations to test the stability of the rocks around the cylindrical excavation: where is a variable that indicates the depth from ground surface and is the depth of the free surface of caved rocks from ground surface.

Equations (7a) and (7b) enable the prediction for the depth () 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 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., 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 () but also the stress due to caved rocks (). Equations (7a) and (7b) provide more details regarding the contribution of the caved rocks. Because of and , the increase of the density (i.e., increases) or height (i.e., 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 () 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].

Figure 3(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 () and the normal to the discontinuities. The slip failure along the discontinuities is most likely to take place, when satisfies 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.

**(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: where is the cohesion of discontinuities, is the coefficient of internal friction on discontinuities, and . When tends to or 90°, the maximum applied principal stress (), 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 , the slip failure along discontinuities due to () is likely to take place before the slip failure due to (), 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 () 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)): where is angle between the stress that accounts for the slip failure () and the normal to the discontinuities on the plane consisted by confining stress () and the stress accounts for the slip failure ().

##### 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., , , and ) are the maximum, intermediate, and minimum applied principal stress (i.e., , , and ), respectively, we substitute Equation (3) into Equation (9), noticing , , and ; we obtain the following relationships, which should be satisfied, if the surrounding rocks are expected to keep stable from the slip failure along discontinuities:

Because the tangential principal stress () is the maximum applied principal stress (i.e., either or is satisfied), Equation (10) can be reduced to the following equations:

Substituting , , , 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:

Equation (12) enables the prediction for orientation () and depth () where slip failure along the inclined discontinuities will take place, and it shows the occurrence of slip failure heavily depends on the property ( 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 is satisfied) is likely to lead to the subsequent failure in other directions (e.g., or is satisfied).

On the other hand, Equations (13a)–(13d) show the impact of caved materials on the stability of surrounding rocks. Because of 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 () 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 ().

#### 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., ) in Equation (1). By regressing 3 sets of laboratory test data, the results demonstrate (), as illustrated in Figure 4.

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:

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.

Table 1 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.

The in situ stress and involved parameters reveal the three principal stress (, , ) satisfies under the free surface of caved rocks (). 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., , , and , 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: where is the angle between discontinuities’ strike and the maximum in situ stress () 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.

As expected, the variation of the critical depth related to rock shear failure (i.e., “shear failure” in Figure 6) is symmetric on either sides of the direction where the maximum effective principal stress () has maximum value (i.e., , N10^{°}W, or the direction of minimum in situ stress). The minimum value of such critical depth appears at , in the direction of minimum in situ stress. However, the depth of implemented undercut locates at , and this means rock shear failure is not responsible for extension of surface subsidence in Xiaowanggou iron mine.

On the other hand, Figure 6 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 , will not take place above the depth in this project. When the 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 (i.e., “slip failure due to ” 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., or ), but the offset towards the direction of minimum in situ stress can be observed. The variation of and 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 and is smaller than the one between and ). The minimum critical depth of undercut appear at when , 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., 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., or “undercut depth” in Figure 6), Figure 6 shows the slip failure due to will take place between and , 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.

The extension process of surface subsidence is illustrated in Figure 7, 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.

#### Data Availability

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

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This study was funded by the Fundamental Research Funds for the Central Universities (grant numbers 2020IVA083 and WUT2019III187), the National Key R&D Plan (grant numbers 2018YFC0808405 and 2018YFC0604401), and the Natural Science Foundation of China (grant number 51774220). The authors would like to thank the Institute of Crustal Dynamics (ICD), China Earthquake Administration, for providing the in situ stress data near Xiaowanggou iron mine.