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

School of Resource and Environmental Engineering, Wuhan University of Technology, Wuhan, Hubei 430070, China Key Laboratory of Mineral Resources Processing and Environment of Hubei Province, Wuhan University of Technology, Wuhan, Hubei 430070, China School of Resources and Civil Engineering, Northeastern University, Shenyang, Liaoning 110819, China China Enfi Engineering Corporation, Beijing 100038, China


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 progres-sively 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][11][12][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 cavingbased 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. 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.

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][22][23][24]. However, the analysis for the impact of caved rocks is rarely reported.

Impact of Caved Rocks on Stability of Cylindrical Caved
Space. In the analytical models proposed in existing literature [14][15][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 σ 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, r b is the radius of the cross section of caved space, z b 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.

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 (θ, z, r), original due to Krisch: where σ r , σ θ , and σ z are the radial, tangential, and axial local normal stress near the rocks around the cylindrical 2 Geofluids 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 h , and b h 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 θ =    3 Geofluids 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 = r b ) can be accordingly calculated from Equation (2): 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., and r b ) 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 ( lim 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 σ 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 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: (1) into Equations (6a) and (6b), we obtain the following equations to test the stability of the rocks around the cylindrical excavation: where z is a variable that indicates the depth from ground surface and z 0 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 4 Geofluids 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., z 0 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 ( lim mented. 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.

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][31][32][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 (σ 1 ) 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. 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 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)): 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 ).

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 5 Geofluids (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: 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), we obtain the following equations to test whether the slip failure along discontinuities will take place around the cylindrical caved space:

and Equation (1) into Equation
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 7 Geofluids 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 ).

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.

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 (R 2 = 0:993), 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: σ H = 0:0304z + 2:9033, σ V = 0:0244z + 1:1593, σ h = 0:0149z + 2:6795: ð14Þ 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 (σ θ , σ z , σ r ) satisfies σ θ > σ z > σ r under the free surface of caved rocks (z 0 = 45 m). 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: 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.
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., θ = 90°, N10°W, or the direction of minimum in situ stress). The minimum value of such critical depth appears at z = 406 m, in the direction of minimum in situ stress. However, the depth of implemented undercut locates at z = 168 m, 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 σ z , will not take place above the depth z = 2500 m 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)  8 Geofluids (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 = 55 m 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 = 45 m 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 = 168 m 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.
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 lip failure due to 휎 휃

Shear failure
Slip failure due to 휎 z 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. 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.
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.

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][14][15][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 largesized caved rocks by dumping small-sized waste rocks or backfilling the caved space with waste rocks.

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 compari-son 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.