Abstract
To deeply understand the mechanism of the rock mass caving and associated surface subsidence during sublevel caving mining, the Xiaowanggou iron mine was selected as an engineering project case study. The study area was analyzed by means of an in situ geological investigation and numerical simulation. First, a borehole television (BHTV) system and a GPS monitoring system were used to monitor the caving process of the roof rock mass and the development of the surface subsidence; the monitoring time was thirteen months. Then, a numerical simulation was used to analyze the damage evolution of the rock mass. Research shows the following: in situ geological monitoring results indicate that the caving process of the roof rock mass presents intermittent characteristics, where slow caving and sudden caving are conducted alternatively and the arched-caving trend is more pronounced with continuous mining. The surface subsidence, horizontal displacement, and horizontal deformation of the hanging wall are higher than that of the footwall, and the subsidence center gradually deflects to the hanging wall in the late stage of the +45m sublevel mining. Numerical simulation results indicate that the extension and penetration of the shear and tensile cracks along the joints and intact rock bridges are the main factors causing the rock mass caving and the existence of the stress arch and its evolution process is the fundamental reason for the intermittent caving of the rock mass. The rapid development of damage to the hanging wall (the damage angle reduced) is the main cause of the deflection of the subsidence center affected by joints and the mining size. In the future of mining, large-scale subsidence will occur on the surface of the hanging wall.
1. Introduction
There are a large number of iron ore resources in the northeastern part of China, and the total reserves of these resources account for two-fifths of the national gross proved reserves. Due to the influence of geological conditions, orebody distribution characteristics, and mining technology, most of these underground iron mines are mined by the sublevel caving method [1–3], which is characterized by extraction on a large scale and high production efficiency with low mining cost. However, sublevel caving mining inevitably results in the formation of underground goafs [4], and the existence of underground goafs will cause the imbalance and redistribution of the stress around the opening, which will cause the displacement, deformation, and failure of the roof strata [5, 6]. More seriously, the failure will cause the roof strata caving to the goaf, greatly threatening the mine operations and personnel safety [7, 8]. The caving process propagates upwards to the surface; thus significant surface subsidence appears and even large-scale collapse pits are formed [9]. As a result, this poses a serious threat to buildings, haulage roads, rivers, industrial facilities, and farmland [10–12]. To control the dangerous situation, it is of great importance to study the mechanism of the rock mass caving and associated surface subsidence.
Rock mass caving is the product of a complex rock mass response to excavation [13], and rock failure is the most important factor causing the occurrence of a caving. Because of the complexity of the underground rock distribution, the partial failure of the rock will lead to a series of chain failure [6] and thus induce massive caving of the roof strata. The formation of surface subsidence is the inevitable result of the caving to the surface.
The rock mass failure and surface subsidence caused by underground mining depend on many factors such as mining depth, geological conditions, and rock mass structure [14–16]. Many studies have conducted in-depth analysis of the mining-induced rock mass failure and surface subsidence, which can be divided into three categories. The first category is the physical similarity simulation method. By establishing a model similar to the actual engineering, the displacement, stress, and failure under the influence of mining can be analyzed. Ju et al. [17] studied the surface stepped subsidence mechanism and proposed the important role of primary key stratum and subkey stratum in the development of subsidence. Sun et al. [18] studied the characteristics of roof strata movement and the evolution of water-conducting fractures and proposed SBM technology to effectively control strata movement and deformation. Huang et al. [19] studied the characteristics of the goaf-roof displacement and strain through physical similarity experiments and reasonably determined the key parameters required for mining design. The second category is on-site monitoring method, such as total station surveys, GPS field surveys, differential InSAR techniques, borehole monitoring, and microseismic monitoring. Xia et al. [20] used borehole monitoring technology to analyze the caving process of the goaf-roof and pointed out that the caving evolution was affected by the structural character of the rock mass. Wang et al. [21] used a combination of on-site monitoring and numerical simulation to analyze the fracture development of the rock mass and surface subsidence characteristics near the riverbed and proposed the Four-Element Structure of antiseepage for the riverbed. Zhao et al. [22] used microseismic technology to analyze the stability of goafs and slope during the transition from open-pit to underground mining and to study the rock failure mechanism. The third category is numerical simulation methods, such as the finite element method (FEM), finite difference method (FDM), and discrete element method (DEM), which have developed rapidly in recent years. Zhao et al. [23] analyzed the rock movement mechanism and surface subsidence characteristics under different stress fields. Xu et al. [24] proposed an EJRM-based numerical method to study the strata and surface movement under the influence of a jointed rock mass. He et al. [25] used a discontinuous deformation analysis to study the formation and stabilization of a stress arch in laminated rock mass.
Based on the above analysis results, each method has its own advantages. However, the physical similarity simulation method is mainly focused on underground coal mining, and the rock structure of metal mines is different from that of coal mines [23]. The on-site monitoring method is mainly used to evaluate and analyze the trend of rock caving and surface subsidence, but cracks and stress evolution in the rock mass cannot be represented visually. When we use numerical methods for analysis, the simulation results may not be consistent with actual site conditions [26, 27]. The needs of accurate study of the mechanism of caving and subsidence cannot be met solely using one method, and there are relatively few studies on the caving mechanism of the rock mass in metal mines. Therefore, using a combination of on-site monitoring and numerical simulation, it is possible to realize a comprehensive explanation of the sublevel caving mining from the macroscopic rock mass collapse model and surface subsidence trend of the caving mechanism of roof strata and the microscopic crack propagation law. This connects mining and rock failure in time and space, thus providing support for mining safety and the protection of surface environment and facilities.
In this paper, the Xiaowanggou iron mine was selected as an engineering project for an in-depth study. The rock mass caving mechanism and surface subsidence characteristics during the sublevel caving mining were comprehensively studied using an in situ geological investigation and numerical analysis method (the FEM-based numerical software RFPA 2D [28–31]). The overall understanding of the process from the crack propagation, stress evolution, and damage development to the rock caving and surface subsidence is obtained. The research can provide guidance for other metal mines with similar conditions.
2. Background of Geologic and Mining Conditions
The Xiaowanggou iron mine is a sedimentary metamorphic magnetite deposit and is located in Liaoyang City, Liaoning Province, northeast China. The mining area is approximately 7.0 km long, and the elevation is approximately 250m. It has a temperate continental monsoon climate with an average annual precipitation of approximately 800mm. The ore body is stratoid, and the strike is north by east 50°, with the dip direction of south-east and the dip angle ranging from 30° to 50°. The strike length is approximately 350m and the average thickness is 100m. Thus, it is a typically inclined metal mine, the deep mining area extends from 190 m below the surface to over 340 m, and the surrounding rock is mainly mixed granite.
The Xiaowanggou iron mine is mined by partition. The entire deposit from top to bottom is divided into three mining areas, namely, the upper mining area, the middle mining area, and the deep mining area, and there is a safety pillar between the partitions [32]. The deep mining area was selected as the research object of this study. The sublevel caving method was used in this mining area, and the first mining section was located at the +60m level, with a subsection height of 15m. At present, the major mining activities are carried out at the +30m level, and the exploration system has been arranged to the 0 m level. In the mining process, the overlying strata face continuous fracturing and caving. In the late stage of +45m sublevel mining, an oval collapse pit with a maximum depth of 27m and a diameter of about 80m was formed on the ground in October 2015. Multiple fracture lines were formed around the collapse pit, which poses a serious threat to the environment, surface industrial facilities, and haulage roads (Figure 1).

This is of great significance for grasping the process of rock mass caving and the trend of surface subsidence in advance, which can effectively reduce the potential safety hazard caused by mining on surface roads and facilities.
3. Analysis of the In Situ Geological Investigation
3.1. In Situ Observations of Rock Mass Caving in Boreholes
The mining of the +45 m section started in May 2014. In order to monitor the roof rock caving situation in real-time, the mine had four geological monitoring bore holes arranged along the 3# exploration line (Figure 2), to be used in combination with the geological exploration results. The average spacing of the boreholes was 30 m, covering the main subsidence area affected by the mining activities. The relevant technical parameters are shown in Table 1.

Monitoring was performed beginning with the completion of the geological borehole project. The BHTV system is mainly used to monitor the caving process of the roof rock mass and it is also widely used to obtain statistics about joints and fissures [4]. The monitoring time was thirteen months and the initial monitoring was carried out once per month. In the late stage of monitoring, with the acceleration of the caving of the rock mass, the monitoring time interval was adjusted. When the roof rock mass in the 2# borehole was only approximately 30 m away from the ground surface, the personnel and equipment were evacuated immediately. On October 5, 2015, a massive collapse occurred on the ground surface.
Scanning images in Figure 3 show that a clear fracture zone appears at a certain distance above the bottom of the borehole. It is noteworthy that the fracture position of the two boreholes differs by 7.3 m; this may be related to the size effect of the rock mass structure [33]. Through regular monitoring, accurate data of the roof caving process are obtained, and the monitoring analysis results are shown in Figure 4.

(a)

(b)

(a)

(b)
In the early stage of monitoring, the stability of the roof mixed granite was good. As the mining advanced, the goaf-roof was always in the slow caving stage (five months). As the mining span continued to increase, the first large-scale collapse took place between January and February 2015, and the caving curve was an obvious arch (Figure 4(b)). During the entire monitoring process, a total of three sudden caving period took place, when the caving heights were 15.2 m, 23.6 m, and 48.1 m (Figure 4(a)). The entire caving project could be described as alternating development with slow sudden caving. In the later monitoring period, the roof strata collapsed rapidly and caved to the surface in October 2015. The total height of the caving was approximately 80 m, and the time to collapse to the surface was only 20 days, indicating that the roof strata activity in this stage was more intense.
3.2. Monitoring of Surface Movement and Subsidence
Surface subsidence is the inevitable result of the rock caving to the surface [34]. In order to ascertain the range and degree of surface rock movement and subsidence caused by mining, GPS monitoring is conducted in the mining influenced area.
According to the geological and mining conditions, a surface monitoring system consisting of two monitoring lines was established along the strike and inclination of ore bodies. As shown in Figure 5, fifteen monitoring points (R1, R2, R3, ... R15) were arranged along the strike and fourteen monitoring points (N1, N2, N3, ... N14) were arranged along the inclination and the interval between monitoring points was 20m. These monitoring points were controlled by six monitoring stations (C1, C2, C3, …C6).

The recovery of the +45m sublevel ore body began in May 2014 and was finished by the end of July 2015, with a total recovery time of approximately fourteen months, and surface subsidence monitoring was carried out monthly. The haulage road is located within the range of rock movement affected by mining; therefore, it is crucial to carry out risk assessment of the haulage road. Since the influence of surface rock movement along the inclination is significant for the haulage road, the study will focus on the analysis of the characteristics of the surface deformation and subsidence along the inclination.
The main characteristics of surface subsidence and deformation are summarized as follows:
With continuous mining of the +45m sublevel ore body, the mining-induced surface subsidence gradually increased. The lowest point of subsidence mainly appeared in the center of the ore body, and the subsidence rate increased, as shown in Figure 6(a). After January 2015, the subsidence center slowly migrated to the hanging wall and the maximum subsidence rate decreased. By comparing the far-field subsidence value, the largest cumulative subsidence in the hanging wall was clearly higher than that in the footwall. By the end of July 2015, the maximum subsidence caused by the +45m sublevel mining was of 1013.3 mm, while the maximum subsidence near the haulage road was approximately 32.5 mm.

(a)

(b)

(c)

(d)
The maximum horizontal displacement appeared in the vicinity of the hanging wall and footwall boundaries of the ore body, not the center. On the contrary, the horizontal displacement at the center of the ore body was almost zero, with the largest subsidence. The maximum horizontal displacement in the hanging wall was about 2 times that in the footwall, as shown in Figure 6(b). During the mining process, the zero point of horizontal displacement moved slowly towards the hanging wall. By the end of extraction, the zero point of horizontal displacement moved to the right by approximately 11.4m, and the maximum horizontal displacement was -285.7mm in the hanging wall. The maximum horizontal displacement near the haulage road was about 27.8m, where the negative value indicated the horizontal displacement in the hanging wall and the positive value indicated the horizontal displacement in the footwall.
The surface horizontal deformation mainly showed wave-shaped changes along the inclination. As the mining advanced, the surface horizontal deformation in the far-field mainly alternated with tensile and compressive deformation. Near the mining area, the horizontal deformation was mainly manifested as tensile deformation and increased gradually, as shown in Figure 6(c). The horizontal deformation near the center of the ore body decreased, and it was worth noting that the horizontal deformation in the hanging wall was significantly higher than that in the footwall and the maximum horizontal deformation in the hanging wall was approximately 1.4~1.9 times that in the footwall. Near the haulage road the maximum tensile deformation was 0.95 mm/m and the maximum compressive deformation was -0.53 mm/m, and several cracks had appeared on the pavement recently (Figure 7). Positive values represented horizontal tensile deformation and negative values represented horizontal compression deformation.

Based on the monitoring results of monitoring point N8, the dynamic subsidence of the surface was mainly divided into three stages, the slow subsidence period, the active subsidence period, and the attenuation subsidence period, as shown in Figure 6(d). When the mining operation reached 28 m, it entered the active subsidence period and the subsidence and subsidence rate increased significantly. The severe subsidence occurred in the range of 42~64 m and the maximum subsidence rate reached 10.7 mm/d. The active subsidence period under the mining-influence had a long duration of approximately 193 days and the subsidence reached 91.5%. When the mining had advanced to 80 m, the subsidence and subsidence rate tended to be stable and the attenuation subsidence period was entered. It was noteworthy that the maximum subsidence rate did not represent the largest subsidence.
Based on the distribution characteristics of the curves of subsidence (similar to s-type) and subsidence rate (normal distribution) in Figure 6(d), a dynamic subsidence function expression is established:
where D0 is the final settlement value of the surface at the monitoring point, mm; x is the horizontal distance between the monitoring point and the mining face, m; is the surface settlement value at the monitoring point when the distance is x, mm; l is the distance between the monitoring point and the mining open-cut, m; v is the mining velocity, m/d, m and n are the fitting coefficients. When mining operations are carried out at a constant velocity v, the surface settlement rate at the monitoring point is expressed as follows:
According to (1) and (2), the curves of subsidence and subsidence rate are shown in the solid line in Figure 6(d), where the fitting curve basically agrees with the monitoring result. The above expression can provide guidance for the dynamic prediction of the process of surface subsidence.
4. Modeling Work
Through on-site monitoring, the macroscopic process of rock mass caving and the development trend of surface subsidence can be clearly observed. However, it is difficult to accurately explain the mechanism of rock mass caving and surface subsidence through on-site monitoring. In order to further study the mechanism of rock mass caving and surface subsidence in the Xiaowanggou iron mine, RFPA 2D numerical software is used in this study.
4.1. The Principle of RFPA 2D
RFPA 2D can be used to simulate the nonlinear deformation and discontinuous medium mechanics of quasi-brittle rocks. The material properties of heterogeneous rocks [35, 36], including Young's modulus, Poisson's ratio and intensity properties, are randomly distributed in the entire analysis domain. Assume these elements are subject to the Weibull's distribution [37], defined as follows:
where is the parameter of the element (such as the Young’s modulus, Poisson’s ratio or strength properties); is the average of the element parameter; m is a parameter defined by the shape of the distribution function and represents the degree of material heterogeneity.
In RFPA 2D, the maximum tensile strain and Mohr-Coulomb criterion are used to define the damage threshold, the former is used to determine whether the element is subject to tensile damage, and the latter is used to determine whether the element is subject to shear damage. The elastic modulus of the damaged material is defined as follows:
where D is the damage variable and E and E0 are the elastic modulus of the damaged and undamaged elements, respectively.
Figure 8 presents the elastic-brittle damage constitutive relation with a given specific residual strength. When the tensile stress in the element reaches the tensile strength , i.e., , the damage variable D can be defined as

(a)

(b)
where is the residual tensile strength coefficient, , and is residual tensile strength. is the strain at the elastic limit and is the ultimate tensile strain at which the element completely damaged in tension, as shown in Figure 8(a). Ultimate tensile strain is given as , is the ultimate strain coefficient.
The Mohr-Coulomb criterion is used as a second damage criterion to describe the element damage under compressive stress conditions (Figure 8(b)), and the expression is as follows:
where and are the maximum and minimum principal stress, respectively; is the friction angle; is the uniaxial compressive strength. The damage variable under compression is described as follows:
where is the compressive strain of the elastic limit, is the residual strength coefficient, and is assumed to be true when the element is under compression or tension.
In numerical analysis, when the element is damaged, its stiffness and strength are reduced and the model is then recalculated under the new parameters. The next load increment is added only when there are no more elements strained beyond the strength-threshold corresponding to the equilibrium stress field and a compatible strain field.
When using RFPA2D for numerical simulation, black cracks will appear in the postprocessing graph after the element exceeds the tensile limit. The simulation of the crack is similar to the fuzzy crack model; i.e., the crack is smeared over the whole element, which greatly simplifies the simulations of crack generation, expansion, and propagation.
4.2. Determination of Mechanical Parameters in Numerical Model
According to the geological survey results, the roof strata mainly consist of mixed granite and there is no major fault-broken zone in the study area. The rock mass is obviously influenced by gneiss structure; it is the main factor in the formation of joints (Figure 9). The joints are mainly composed of tectonic and weathered fractures. The existence of joints has a significant influence on the development of the caving and the characteristics of surface rock movement [24, 38–40], which usually leads to a large number of unstable rock blocks in the formation. The detailed survey of the structural plane of the roof rock mass was carried out by manual in underground roadways, the main survey parameters included occurrence, inclination, dips, and spacing of joints, and the measured parameters were then analyzed using DIPS software [4]. Ultimately, we identified the dominant structural planes. Based on geological survey results, there are mainly three joint sets in the rock mass. Figure 10 shows the distribution of joint sets in the study area. The relevant parameters of the joint sets are shown in Table 2.


As shown in Figure 10 and Table 2, the optimal joint sets obtained from the calculation and analysis consist of a gently dipping joint and two almost orthogonal steeply dipping joints; this rock mass structure is the most favorable for collapse [41]. For simplicity and convenience, two-dimensional numerical models are constructed by appropriately simplifying the three joint sets with ∠0°, ∠70°, and ∠80°. The mechanical parameters of the roof mixed granite are mainly obtained through laboratory tests [42, 43], as shown in Table 3. The elastic modulus of the rock mass is calculated using formula (8) [24].
where E is the elastic modulus of the rock mass and RMR is the geomechanics classification index of the rock mass [44, 45].
When the in situ geological model is used for numerical simulation, it is very important to correctly choose the input parameters but also the correct boundary conditions and dimensions of the model, which directly determine the reliability of the simulation results. In combination with the actual surface subsidence range obtained in the field, the model dimensions should be greater than this range, and the boundary conditions can refer to relevant references [28, 29, 31]. Here, we highlight how to determine the physical-mechanical parameters used in numerical simulation through the measured physical-mechanical parameters of the rock mass and the homogeneity coefficient.
In order to solve the practical problems, a numerical model is first established. Its geometrical dimensions are on the same order of magnitude as the actual engineering problems to be studied. The numerical model size is 200 m×100 m and is divided into 20,000 finite element grids. In this paper, the coefficient of homogeneity of the sample is selected as m=3, and the relevant parameters input in the model are selected according to the fitting formulas (9) and (10) [46].
where E0 and fc0 represent the elastic modulus and intensity values calculated in the model, respectively. E and fc represent the macroscopic elastic modulus and intensity values of the model (measured values), respectively.
As shown in Table 3, the compressive strength (fc) of the mixed granite is 55.3 MPa and the elastic modulus (E) is 3.5 GPa. The calculation inputs for the compressive strength and elastic modulus are calculated according to formulas (9) and (10); the values are 68.9 MPa and 11.3 GPa, respectively. Subsequently, the uniaxial compression numerical simulation of the model is carried out, and we observe the peak intensity of the model. Because the numerical calculation has good repeatability and maneuverability, the initial value of the compressive strength and elastic modulus can be adjusted within a narrow range, until the compressive strength of the model obtained is consistent with the rock mass strength of the actual engineering. As shown in Figure 11, the peak strength of the model has reached 54.9 MPa, which is very close to the measured compressive strength of the mixed granite.

The average length and spacing of the joints based on the geological survey results are approximately 25 m and 15 m, respectively. It should be noted that the joints used in the model are assumed to be made of a ‘‘weak material’’ with lower strength and stiffness. Wong et al. [47] reported that the homogeneity index should be greater than 2.0 but then fall within the typical range of 2.0–6.0. Therefore, m=3 is chosen to characterize the heterogeneity of the rock mass. Finally, we determine the input parameters of the numerical model as shown in Table 4.
4.3. Numerical Modeling
In this study, the geological profile of exploratory line No. 3 is selected and modeled for the subsequent numerical simulation (Figure 2: minor surface undulations and contact zones are ignored). The numerical model is shown in Figure 12; the analysis domain is 600×250 m, and the model is divided into a mesh with 150,000 elements.

The top boundary is set as free, and normal displacements are constrained on the right, left, and bottom boundaries; the model is affected by its own gravity. The draw is simulated by gradually removing the rock mass at a depth of 175 m (the 1st mining layer) and, for each step, the excavation size is 15 m in height and 10 m in length. The model is calculated in a quasistatic manner to achieve a balanced state, and the calculation process is assumed to be a plane strain problem. Then, we will describe the progression of rock mass caving and associated surface subsidence through the excavation with an advance increment of 10 m/step.
In the caving process fracture initiation always starts from the joint sets and propagates into the rock bridges [48]. In order to highlight the role of rock bridges in the caving of rock masses, we assigned a persistence <1 to the joint sets and the nonpenetrative area between joints in the numerical model is considered as the rock bridges. The uneven spatial distribution of the structural planes in the rock mass must inevitably lead to uneven distribution of length and position of rock bridges [49, 50]. Affected by the existence of joints, rock bridges are generally randomly distributed in the rock strata. Therefore, in numerical models randomly distributed rock bridges were used for simulation and the rock bridges were not quantified in detail.
5. Numerical Simulation Results
5.1. Rock Mass Caving
Figure 13 shows the numerical results, where the crack failure process, stress distribution, and damage development can be more intuitively demonstrated throughout the numerical simulation process. In the initial stage of excavation, the initial balance within the rock mass is disturbed, causing the stress to be concentrated on both sides of the undercut (the shading intensity indicates the relative magnitude of the shear stress). At the same time, the goaf-roof produces tensile and shear cracks along the joints under self-weight. In the damage development graph, the red cells represent tensile damage and the white cells represent shear damage (Figure 13(a)). As the excavation continues, when the advancing distance reaches the ultimate caving span [51], tensile and shear cracks interconnect along the joints and intact rock bridges, and the roof strata begin to collapse in the goaf until it forms a self-stabilizing structure. Tensile damage occurs at the bottom of the roof strata, and the shear damage continues to expand upward along the joints in this process. However, the rock mass on both side walls shows different caving mechanisms. The sidewall rock mass at the bottom is mainly sliding failure, which is conducive to the development of failures towards the footwall. It gradually turns into a toppling failure with the continuous caving, which is conducive to the development of failures towards the hanging wall (Figure 13(c)). It is worth noting that one or more obvious stress arches [52, 53] are formed above the self-stabilizing roof strata (the zone presented with brighter color) and the stress is mainly concentrated within the stress arch, while the stresses below are diminished. The rock mass caving is mainly due to the tensile failure. The presence of stress arches limits the further collapse of the roof strata, and it also causes the rock mass to exhibit an arched-caving characteristics. The roof strata entered a slowly caving period at this time, as shown in Figure 13(b). The stress arch must be broken for mining to continue; this will force the support points of the stress arch to move upwards and, with the continuous caving of the rock mass below the stress arch, it enters the suddenly caving period. When the caving line approaches the new stress arch, it once again enters the slowly caving period and the arched-caving trend is more obvious. Meanwhile, there has been a small subsidence on the surface, but the roof strata have not completely collapsed and the rocks of the side wall remain in a stable state (Figure 13(d)). In the whole mining process, the collapse rate of roof rock mainly presents periodic characteristics. It can be expressed as alternating development with slowly suddenly collapse and eventually surface subsidence. This caving mechanism is basically consistent with the on-site monitoring results of the Xiaowanggou iron mine.

(a)

(b)

(c)

(d)
5.2. Surface Subsidence
The numerical results of the +45 m sublevel excavation are shown in Figure 14. When the roof rock caves to the surface, the stress arch has completely disappeared and forms a significant subsidence trough on the surface, the tensile and shear damage has fully developed to the surface, and the subsidence is large in the central area and decreases progressively towards the sidewalls. At the same time, we observe that the stress concentration in the hanging wall is more pronounced than that in the footwall and in turn accelerates the development of destruction and caving of the hanging wall rock mass; surface disturbances within this zone are caused by changes in the stress field, and the damage development has fully extended to the surface and continues to extend to the sidewall (Figure 14(a)). With further excavation, the rock mass in the nonsteady area of the hanging wall is separated and destroyed, which resulted in a corresponding increase of the surface subsidence range and internal damage development (Figure 14(b)). When the excavation distance reached 80 m, relatively large toppling failure of the rock mass on both sidewalls has occurred along the joints (Figure 14(c)) and these joints are more conducive to the damage development of the hanging wall rock mass. After +45 m sublevel excavation was completed, the rock mass on both sidewalls has been destroyed and subsided completely and the damage angle of each sidewall has further decreased (Figure 14(d)). In this study, the damage angle is obtained based on the development of tensile damage.

(a)

(b)

(c)

(d)
The numerical results of surface subsidence for the +45 m sublevel at different excavation distances are shown in Figure 15. The depth of the subsidence trough increases as the excavation distance increases. The maximum subsidence value is produced in the mining center area; subsidence decreases as the distance from the central area increases. Eventually, the center of the subsidence trough deflects towards the hanging wall and the asymmetrical feature is remarkable. The concept of the subsidence trough was gradually used to assess the degree of surface subsidence caused by mining; this concept takes into account one of the most important observed facts about subsidence, namely, that the surface area affected is larger than the mined area [54]. There is a certain deviation between the maximum subsidence values obtained by numerical and on-site monitoring, with the main reason being that the bulking effect of the collapsed rock mass is not taken into account in the numerical analysis, which may lead to a large value in the numerical monitoring. However, the numerically obtained characteristics of the surface subsidence are basically the same as the on-site monitoring as already mentioned, which can better explain the on-site surface subsidence mechanism.

6. Analysis and Discussion
In this paper, the Xiaowanggou iron mine is taken as an engineering case study, and the combination of on-site monitoring and numerical simulation is used to study the mechanism of the rock mass caving and related surface subsidence in sublevel caving. The process of rock mass caving and the development trend of surface subsidence are obtained through an in situ geological survey. The processes of crack, stress, and damage evolution in the rock mass during mining are obtained through a numerical inversion. This can be used to recognize the whole failure process from the beginning of the caving to the formation of surface subsidence, allowing for an intuitive understanding of the rock mass caving and related surface subsidence mechanism. The research results can provide guidance for other mines with similar mining conditions.
According to the in situ geological survey, the arched-caving development of the rock mass has been verified. The entire caving process can be described as alternating development with slow sudden caving; it is concluded that this type of development is related to the arch effect of rock mass caving. The effect of mining the +45 m sublevel ore body on surface subsidence and horizontal deformation is significant. Although the subsidence rate of the hanging wall is larger than that of the footwall, the horizontal deformation near the haulage road reached 0.95 mm/m by the end of May 2015 and some cracks had been observed in the field. There is no doubt that the mining of the +30 m and +15 m sublevel ore bodies will increase the risk of road breakage and collapse.
In the analysis of rock mass caving, the crack failure occurs mainly along the randomly distributed joints (weakly constructed) and extends into the intact rock bridges, as the failure path of rock mass is not straight but flexural. The roof strata will not collapse to the goaf until the shear and tensile cracks have completely penetrated. Once the intensity of the stress concentration in a caved zone of the roof strata is released, a stress shadow area (the zone presented with greyer color) and an apparent stress arch (the zone presented with brighter color) above the caved zone will be produced, as shown in Figure 13 (the stress distribution). When one or more key rock bridges are damaged, the stress arch will break immediately. With the increase in advancing distance, the stress arch is repeatedly formed and dissipated and continues to develop upward. Finally, the roof strata collapse to the surface to form a significant subsidence trough with the destruction of the last stress arch. Understanding the caving mechanism of the rock mass during the mining is an essential aspect for accurate prediction of the rock fracture propagation and the surface subsidence development; it can aid in the timely adjustment of the mining sequence, put forward preventive measures, and provide an effective guarantee for mine safety exploitation and surface environment and structural protection.
In the analysis of ground subsidence, the concept of the damage angle is introduced. The curve of the damage angle with different excavation distances is shown in Figure 16; all data are obtained from the numerical models and each point corresponds to a different model, where damage angles decrease with an increasing excavation distance. It is worth noting that, in the 1st mining layer (+60 m sublevel), the damage angle of the hanging wall is always larger than that of the footwall, and the angle difference decreases as the advancing distance increases, indicating that the rock mass in the footwall is affected greatly in this stage. With the increase in excavation distance, this effect is gradually reduced. Because of the existence of the stress arch, there is no significant subsidence on the surface. After the ore body of the +45 m sublevel is excavated, the damage angle has been deflected and the destruction and caving of the hanging wall rock mass gradually play a major role. In the excavation process of the first 60m, the subsidence trough is basically located at the center of the mining area. This is the result of the interaction of sliding and toppling damage under the influence of joints. In subsequent excavations (80~100 m), the damage development rate of the hanging wall is obviously higher than that of the footwall affected by joints. The subsidence trough is gradually deflected to the hanging wall. The comprehensive analysis has shown that the joints and mining size are two important factors that cause different subsidence mechanisms at different mining stages. It can also be observed that failure potentially occurs at the hanging wall zone with continuous mining. With the decrease in the bearing capacity of the sidewall rock mass induced by mining, the subsidence range will be further expanded. After the mining of the +30 m, +15 m, and 0 m sublevels is finished, the haulage road is bound to be destroyed and large-scale subsidence will occur at the surface of the hanging wall.

7. Conclusions
In situ geological surveys and numerical analyses were conducted to research the mechanism of rock mass caving and surface subsidence in the Xiaowanggou iron mine. The main conclusions were the following:
Based on an in situ geological survey, the entire caving project can be described as alternating development with slow sudden caving, and the roof strata activity in the later monitoring period is more intense. With the continuous extraction, the subsidence center slowly migrated to the hanging wall and the maximum subsidence rate decreased. The maximum horizontal displacement was -285.7mm in the hanging wall, and the horizontal deformation of the hanging wall was significantly higher than that of the footwall. The mining of the +30 m and +15 m sublevel ore bodies will increase the risk of road breakage and collapse.
Based on the numerical simulation, the crack failure occurs mainly along the randomly distributed joints and extends into the intact rock bridges; when one or more key rock bridges are damaged, the stress arch will break immediately. When the rock mass is not completely collapsed to the surface, the damage angle of the hanging wall is always larger than that of the footwall. After the rock mass collapse to the surface, the damage development rate of the hanging wall is obviously higher than that of the footwall. Ultimately, the subsidence trough is gradually deflected to the hanging wall and large-scale subsidence will occur at the surface of the hanging wall.
The arched-caving development of the rock mass is verified in the study; the root cause of the rock mass caving is that the shear and tensile cracks are interlinked along joints and intact rock bridges. The intrinsic factor of the slow sudden caving being periodic is the result of the stress arch evolution. The rock mass will not collapse to the surface until the last stress arch breaks.
The development of surface subsidence is seriously affected by joints and mining size. The numerical result is basically consistent with the result of the on-site monitoring; it is feasible and reasonable to use the damage angle to analyze the characteristics of surface subsidence. With the continued mining of the ore bodies, the subsidence area will be further expanded, and the haulage road on the footwall is bound to be destroyed. It is necessary to take protective measures as soon as possible.
Combined with the in situ geological investigation and numerical simulation, the rock mass caving mechanism and ground subsidence characteristics in a jointed rock mass are clarified. This research can be significant for mining safety, allowing mining technicians to better understand the internal causes of rock caving and surface subsidence; this can prevent the hazards caused by the caving and subsidence in advance.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Authors’ Contributions
Fengyu Ren, Dongjie Zhang, and Miao Yu contributed to the formulation of overarching research goals and aims; Dongjie Zhang and Jianli Cao analyzed the data; Dongjie Zhang wrote the paper; Dongjie Zhang and Shaohua Li contributed to the numerical calculation and simulation.
Acknowledgments
The study is jointly supported by grants from the Key Program of National Natural Science Foundation of China (Grant no. 51534003); National Key Research and Development Program of China (Grant no. 2016YFC0801601); National Key Research and Development Program of China (Grant no. 2016YFC0801606). The authors are grateful for these supports.