Numerical Simulation of Crack Evolution Mechanism and Subsidence Characteristics Effected by Rock Mass Structure in Block Caving Mining

In the block caving mining, the significant rock mass deformation and surface subsidence will be formed with the continuous extraction of ore. However, the internal crack evolution mechanisms in rock mass and associated subsidence characteristics present one of the key issues in rock mining engineering. Although block caving method has been used for many years, current knowledge of the crack evolution mechanisms, the subsidence characteristics under the influence of rock mass structure and subsidence prediction capabilities are limited. Based on the rock mechanics model provided by CEMI (Centre for Excellence in Mining Innovation), crack evolution mechanisms and subsidence characteristics effected by the rock mass structure in block caving are numerically investigated using RFPA 2D, a numerical code based on FEM (Finite Element Method). Crack formation, propagation and coalescence in the overlying strata and the stress-balancing arch evolution in the stress field are represented visually during the whole process of extraction. The numerically obtained crack evolution shows that the stress-balancing arch has a significant influence on the fracture development of rock mass, and directly determines the slump form and rate of the rock mass. After understanding of the crack evolution mechanism in rock mass, the characteristics of surface subsidence are analyzed. Numerical experiments emphasize the geometrical configuration of joints and faults about mechanisms of subsidence development, including joints orientation, faults location and inclination, which can provide significantly meaningful guides for investigation of subsidence mechanisms and implementation of remedial measures.


Introduction
Resources and environmental issues are the two major challenges of world development. With the continuous exploitation of mineral resources, mining depth increased from hundred meters to several kilometers, and the rock movement and surface subsidence caused by mining became one of the most significant problems in geotechnical engineering.
For large low-grade metal mines, block caving is widely used underground mining technology (Brown Abstract In the block caving mining, the significant rock mass deformation and surface subsidence will be formed with the continuous extraction of ore. However, the internal crack evolution mechanisms in rock mass and associated subsidence characteristics present one of the key issues in rock mining engineering. Although block caving method has been used for many years, current knowledge of the crack evolution mechanisms, the subsidence characteristics under the influence of rock mass structure and subsidence prediction capabilities are limited. Based on the rock mechanics model provided by CEMI (Centre for Excellence in Mining Innovation), crack evolution mechanisms and subsidence characteristics effected by the rock mass structure in block caving are numerically investigated using RFPA 2D, a numerical code based on FEM (Finite Element Method). Crack formation, propagation and coalescence in the overlying strata and the stress-balancing arch evolution in the 2003), which is characterized by caving and extraction on a large scale, high production efficiency with low mining cost. The rock movement and surface subsidence caused by mining is closely related to the rock mass structure, which involves a series of rock stratum movement including the generation of tension cracks, spread of scarps, fracturing and caving towards the surface, eventually formation of a large crater on the surface (Beck et al. 2006;Vyazmensky et al. 2010). Therefore, it is necessary to understand the characteristics of rock crack propagation, stress distribution rules and the settlement mechanism under the influence of rock mass structure, which play an important role in mining operational hazards and environmental impact assessments.
For rock mass, the in-situ stress is a natural stress before the occurrence of mining activity. Compared with rock excavation, it can be seen as a constant stress field. In the rock mass, the stress balance will be destroyed when the stress near the excavation is redistributed to produce the secondary stress field, which leads to the deformation and failure of the roof strata (Zhang et al. 2021a, b). The formation of stress arch makes the overlying strata stabilize itself. Then the weight of the overlying strata is mainly converted to compressive stress and transferred to the bottom through the stress arch (Li et al. 2006a, b;He and Zhang 2015) , which leads to stress redistribution and controls the further development of collapse.
However, with the increase of the undercut size, the stress arch is unable to carry the weight of the overburden and collapse again. Thus, the overlying strata are continuously broken and fall into the goaf, causing the collapse process that would eventually spread to the surface. Moreover, this collapse can result in continuous or non-continuous surface subsidence, and deteriorate the surface environment and infrastructure. At present, this development process in the actual mining operations can't be directly observed, and there is an urgent need for a research methodology to elaborate on the above mentioned development process.
The rock movement and surface subsidence result from the response of the complex rock excavation, involving the complicated kinematic mechanism. At present, the study of rock movement and surface subsidence caused by mining is mainly focused on underground coal mining. The main methods include probabilistic integral, profile function, physical similarity test, fuzzy mathematics and numerical simulation (Cao et al. 2020;Zhang et al. 2020;Sui et al. 2021;Hekmatnejad et al. 2021). The collapse is the result of the response of a complex rock excavation. This response includes comprehensive failure of the rock mass in tensile and shear, along the weak fractures and through intact rock bridges, which lead to the limitation of using existing theory to study the mechanism of rock movement and ground subsidence (Zhang et al. 2021a, b). Meanwhile, the technologies of surface subsidence monitoring (Li et al. 2019;Brencher et al. 2021;Mei et al. 2021), including leveling, total station surveys, GPS field surveys and differential InSAR techniques, have been used for mining-induced subsidence monitoring. These methods are mainly used to monitor surface deformation, while crack evolution and failure-induced stress redistribution in the rock mass are still unable to observe. Therefore, it is very important to study the mechanism of internal crack propagation in rock mass and surface subsidence vis-a-vis its influence on rock mass fabric, which is critical for protecting both surface environment and ground infrastructure.
Although some studies related to rock movement and surface subsidence have been carried out (Chen et al. 2016;Chang, et al. 2021;He, et al. 2021), it is still a challenge for the quantitative assessment of the crack evolution and subsidence mechanisms caused by block caving mining. One of the main obstacles is that the process of crack evolution is not visible in the overlying strata. Through using advanced numerical tools, crack evolution rules and surface subsidence characteristics under the influence of rock mass fabric can be systematically studied, which facilitates a more detailed understanding of the surface subsidence mechanism.
In this paper, we study the characteristics of surface rock movement and subsidence affected by joints and faults through a series of numerical experiments. The mechanism of crack evolution and surface subsidence is clearly investigated.

The Principle of RFPA 2D
Block caving method is suitable for both loose fragmented and hard rock conditions, the main failure in the mining is a brittle fracture. The surface subsidence is a result of the excavation response to complex rock mass, involving comprehensive failure under tensile and shear stress. Therefore, a numerical method, which is based on elastic-brittle fracture mechanics and combined with the coupling compression and tensile strength standard, are used in the study.
RFPA 2D is a two-dimensional numerical tool used for the progressive damage to overall instability in quasi-brittle materials. This software was developed by Tang (1997) and improved at Mechsoft, China. By leading the non-uniformity of the rock material into the numerical model, the nonlinear deformation and discontinuous medium mechanics of quasi-brittle rocks can be simulated. The material properties of heterogeneous rocks, including Young's modulus, Poisson's ratio and intensity properties, are randomly distributed in the entire analysis domain. The rock material used in the numerical simulation consists of many elements that are assumed to be subject to the Weibull's distribution (Weibull 1951) and defined as follows: where is the parameter of the element (such as the Young's modulus, Poisson's ratio or strength properties), 0 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. m is homogeneous material.
In RFPA 2D, the maximum tensile strain and Mohr-Coulomb criterion are used to define the damage threshold, where the former is for tensile damage and the latter is for shear damage. The elastic modulus in the elastic damage mechanics, as expressed in Eq. 2, it will decreases with the progress of the damage. ( where D is the damage variable, E and E 0 are the elastic modulus of the damaged and undamaged elements, respectively. Figure 1 shows the relationship between elasticbrittle damage constitutive and specific residual strength. When the tensile stress in the element reaches the f t0 (i.e.,) 3 ≤ −f t0 ), the damage variable D can be defined as: where is the residual tensile strength coefficient, f tr = ⋅ f t0 , f tr is a residual tensile strength, t0 is the strain at the elastic limit, and tu is the ultimate tensile strain at which the element completely damaged in tension (Fig. 1a). The ultimate tensile strain is given as tu = ⋅ t0 , and is the ultimate strain coefficient.
The Mohr-Coulomb criterion is utilized as a second damage criterion to describe the element damage under compressive stress conditions (Fig. 1b). The expression is as follows: where 1 and 3 are the maximum and minimum principal stress, respectively. is the friction angle, and f c0 is the uniaxial compressive strength. The damage variable under uniaxial compression is described as follows: where c0 is the compressive strain at the elastic limit, is the residual strength coefficient, and f cr ∕f c0 = f tr ∕f t0 = is assumed to be true when the element is under uniaxial compression or tension. According to the derivation of the above mentioned formula, the elastic modulus value of the damaged element under different stress field can be obtained.
In numerical analysis, when the element is damaged, its stiffness and strength will be reduced and the model will be subsequently 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, the black cracks will appear in the post-processing graph after the element exceeds the tensile limit. The simulation of the crack is similar to the fuzzy crack model; the crack is smeared over the whole element, which greatly simplifies simulations of crack generation, expansion, and propagation. With the development of computer technology, more researchers have attempted to use a similar principle in solving discontinuous problems through continuum mechanics (Fang and Harrison 2002;Ma et al. 2011). The main feature of RFPA2D is that simulating the crack formation and propagation without a pre-existing crack. The advantage is that the mesh topology is left untouched (Pearce et al. 2000). Similar methods are also applied to numerical modeling (Pietruszczak and Xu 1995;Pan et al. 2014).

Modeling Work of Numerical Simulation
As shown in Fig. 2, an obvious collapse pit will be formed above the center of the mining area, and the surrounding rock cracks on the sidewall of the collapse pit will develop and move toward the collapse pit. Surface rock movement can be divided into three Pre-mining topography characteristic areas: caving area, fractured area and subsidence area. The caving process is the macroscopic mechanical response of the rock mass in the two failure states of tension and shear failure, this response is accumulated by local failure (fracture initiation, expansion along the structural plane, rock bridge failure, etc.). The rock mass stress evolution characteristics are important reference for analyzing the development law of rock movement. Using numerical analysis methods, this feature can be displayed intuitively, which helps us better understand the caving characteristics of rock mass and the mechanism of surface rock movement in caving mining.
Therefore, Centre for Excellence in Mining Innovation (2010) provides a rock mechanics model as shown in Fig. 3, requiring a numerical analysis method to study the variation of the caving process and the surface deformation along with the mining. This study is based on the numerical modeling of the rock mechanics model to clarify the process of crack evolution in the rock mass during the block caving mining and confirm the existence of the stress-balancing arch. At the same time, the mechanism of rock migration and surface subsidence under the influence of rock mass structure (joints and faults) is clarified.
In the model of Fig. 3, the deposit is composed of four rocks, including Rhyolite, Quartz-Monzodiorite, Sandstone and Siltstone, and Biotite granodiorite. The ore body boundary is not clearly defined. There are three faults in the deposit, and the existence of faults has a significant effect on the failure pattern of the rock mass. When a significant fault exists in the mining area, a moderate to steep subsidence movement will occur along the fault, regardless of the cave angle through intact rock (Song et al. 2010;Villegas et al. 2011). However, the main focus of modeling analysis in this study is to clarify the process of crack evolution in rock mass and the mechanism of surface subsidence. This paper adopts the following three research methods: (1) Under the excavation initiation depth, the crack propagation and stress evolution of the rock mass are studied, and joints are randomly distributed from the depth of 160 to 320 m with the 10 m increments. (2) Under the 160 m undercut size, the influence of joint sets orientation on the development of surface is investigated, and the whole surface subsidence profile is defined by using the subsidence angle (the subsidence limit > 0.05 m). The effect of the faults is ignored, and the faults are considered to be a perfectly strong interface. (3) At a certain undercut size (160 m), the influence of fault location and inclination to the development of surface subsidence is further discussed, and random distribution joints are adopted.
In this paper, RFPA 2D is used to simulate the rock failure process caused by block caving, and the reliability of this software has been validated in previous studies (Li et al. 2006a(Li et al. , b, 2011a. Flores and Karzulovic (2002) studied some block caving mines and reported average caved ore block heights of around 200 m. The numerical model is shown in Fig. 4, the analysis domain is 600 × 200 m. The model is divided into a mesh with 120,000 elements, and the caving depth is 160 m. We arrange seven measuring points above the left side of the bottom cutting, the distance between P1 ~ P4 measuring points is 10 m, and the distance between P4 ~ P7 measuring points is 30 m.
The top boundary is set as free, and normal displacements are constrained on the right, left and bottom boundaries. The draw is simulated by gradually removing the rock mass at a depth of 160 m until the The model is calculated in a quasi-static manner to achieve a balanced state, and the calculation process is assumed to be a plane strain problem. In order to investigate the crack propagation and evolution characteristics of the rock mass, randomly generated joints are applied in the model. The average length and spacing of the joints are approximately 30 m and 20 m, respectively. Wong et al. reported that the Weibull parameters should be greater than 2.0 (Wong et al. 2006), but then fall within the typical range of m = 2.0 ~ 6.0. In this paper, combined with the previous research results, when using RFPA software for analysis, the heterogeneity coefficient of different rocks is generally between 3 and 4 (Zhang et al. 2009Li et al. 2011a, b). Therefore, Weibull parameters (m = 3.5) are chosen to characterize the heterogeneity of rock mass. It is assumed that Young's modulus and other mechanical parameters satisfy the respective homogeneity index and vary independently. The physical and mechanical parameters of rock used in this model are listed in Table 1.
A series of parametric numerical experiments were performed to clarify the crack propagation and stress evolution rules of the rock mass in block caving, and the relative significance of rock mass structure (joints and faults) in the surface subsidence development was investigated the following sections.  the figure represent subtle fissures, which are connected to form fractures. During the excavation process, the rock mass continues to destroy until collapsing to the surface for the formation of a clear subsidence trough. As shown in Fig. 5d, the subsidence in the central area of the ore block is large and decreases progressively towards the sides. The numerical simulation results in Fig. 5 can clearly exhibit the entire development of the rock from fracture to collapse, which is also a prominent advantage of the RFPA2D program. It can be seen that the fracture and failure are located along existing structural plane and through intact rock bridges in the action of tensile and shear stress.

Crack Evolution in Rock Mass
During the excavation process, the stress evolution in the rock mass is shown in Fig. 6. The shadow intensity represents the magnitude of the maximum shear stress in the element. With the increase of undercut size, the macro-cracks cause the stress-concentrated area extends upward to the shape of the arch. Then a plurality of stress-balancing arches are formed above the excavation area (Fig. 6a), namely stress-balancing arch, which can limit the further collapse of the roof rock stratum and is basically consistent with the voussoir beam theory. The notion of the voussoir was first established by Evans to explain the stability of a jointed or cracked beam, which now has been applied to explain the roof collapse behavior and assess the roof stability (Hatzor and Benary 1998;Bakun-Mazor et al. 2009). Classical voussoir beam theory states that joints might extend upwards until a stable arch is formed inside the rock mass under gravity. Moreover, the voussoir beam could stabilize itself through forming an arch, which is the stress-balancing arch. Through numerical simulation, it can be found that the failure path of the rock mass is not straight but flexural. Generally, there are two kinds of damage modes for different rock materials, namely high-stress damage and low strength damage. High-stress damage occurs mostly in homogeneous materials while low-strength damage occurs mostly in heterogeneous materials, which can be attributed to fact that micro-cracks, joints and faults could generate the development of fissure at the weaker position in the heterogeneous material (e.g., rock). Therefore, the heterogeneity of the material is the key factor for occurring damage in which the stress is not the largest. Once the stress concentration strength of the roof rock above the excavation area is released, an apparent stress-balancing arch will be produced (Fig. 6b). , the displacements of measurement points 1 and 2 are significantly increased to 2.81 m, and other monitoring points did not produce a large displacement. When the excavation reaches 110 m (point C), the measurement point 5, where the roof began to collapse, and the maximum displacement is 3.58 m. The top two measuring points 6 and 7 collapses after 130 m excavation, resulting in a maximum vertical displacement of 3.72 m. Based on the observation of the measuring point displacement, it is suggesting that there is an significant delay and jump in the displacement of the roof strata with the increase of the undercut size. The collapsed form is usually not a large-scale disposable but an intermittent failure. The basic collapse process can be described as follows: slowly development, suddenly collapsed, and eventually surface subsidence. The collapse process clearly illustrates the dynamic equilibrium process. Furthermore, the stress-balancing arch dose not exists, the rock mass will not completely cave to the surface until all the stress-balancing arches are broken. In the measuring point displacement at the end of excavation, the reason for the slight change may be related to the sinking effect of the bottom broken rock mass.
According to the simulation results, it is indicated that the crack propagation occurs mainly along the randomly distributed joints (weakly constructed) and extends into the intact rock bridge, thus forming a stress-balancing arch. Once one or more key rock bridges are damaged, the stress-balancing arch would be broken immediately. With the increase of undercut size, the stress-balancing arch is repeatedly formed and dissipated, and continues to develop upward until the last stress-balancing arch breaks. The roof rock mass collapses to the surface and form a significant subsidence area. Mahtab et al. 1973 stated that the rock mass structure, which is most favorable for collapse, consists of a gently and two almost orthogonal steeply dipping joints. This study will focus on the influence of joints orientation on surface subsidence, and presents the following five modeling schemes (Table 2). Among the above mentioned five models, A1, A2, and A3 are used to study how the orientation of joint sets effect the surface subsidence development. A4 and A5 based on the A3 are applied to assess the surface subsidence development significance of changing the orientation for gently dipping joints and the presence of additional vertical joints. Figure 8 shows the mechanism of surface subsidence development under five models with 160 m undercut size. With the increase of excavation distance, the crack inside the rock mass is continuously expanded and unloaded, which has favorable movement conditions for the separation of the overlying strata adjacent to the goaf. Eventually, the surface forms a clear subsidence crater. After the rock mass caves to the surface, the increase of the mining span leads to the continuous decline of the fragmented rock in the crater, decreasing the lateral support to the crater walls, and it causes the sliding and/or topping failure along the crater wall, which further increases the subsidence deformation range. The ground subsidence mechanism is consistent with that proposed by Abel et al. (Abel and Lee 1980) .

Impact of Joints
It can be inferred that the expansion of cracks in the rock mass, the propagation of the caving to the surface, and the destruction mechanism of the nearsurface rock mass are closely related to the joint pattern. The joints pattern controls not only the crack propagation mode but also the surface subsidence range. In order to quantify the influence of joints pattern on the deformation characteristics of surface subsidence, this study adopts 5 cm displacement threshold to express the main surface deformation extent. The contours of 5 cm displacements at 160 m undercut size for A1, A2, A3, A4 and A5 models are used to define the Rock mass Deformation Range (RDR) (Fig. 9). Then, subsidence angle is adopted to describe the range of the main surface subsidence caused by excavation, as shown in Fig. 8.
The subsidence angles with different joints orientation that controls major surface deformation areas are compared. The subsidence angles range from 67° to 82° on the right side of the undercut, and it generally does not exceed the inclination of the steeply sloping joint sets. It indicates that the failure mode of the rock mass is mainly controlled by the inclination of the sub-vertical/steeply sloping joint sets, mainly the sliding failure along the steeply sloping joint sets, limiting the extent of the rock mass deformation and surface subsidence range. On the left side of the undercut, the subsidence angles range from 72° to 83°, it is mainly the topping failure along the joint sets. It is found that the subsidence angles are generally higher than the right side (except for A2 model). This variation tendency can tentatively be explained that the sliding failure of the rock mass along the steeply sloping joints is greater than the extent of topping failure.
In order to further study the surface subsidence deformation characteristics in block caving, the major surface disturbance area can be characterized by the overall subsidence range, where the vertical and horizontal surface displacement (≥ 5 cm) are used for analysis (Figs. 10 and 11).
Taking A1 model as a reference, with a 10° joint rotation(A2 model), the surface subsidence range is increased by 11%, and the deformation volume of the rock mass increases by about 4% (Fig. 9), rotation of the joints by 20° (A3 model) leads to an increase in the range of surface subsidence by about 45%, and the RDR by about 15%, it is pointed out that the influence of the joints rotation angle on the surface subsidence range is obvious, and this changing trend mainly depends on the toppling and sliding failure mechanism of rock mass controlled by steeply A5 inclined joint sets. In the A4 model, the extent of surface subsidence is reduced by about 22%, and the RDR by about 15% (relative to the A3 model), indicating that by reducing the inclination of the gently dipping joints under certain conditions of the steeply sloping joints can effectively limit the range of surface subsidence. In model A5, the additional vertical joints increase the subsidence angles on the left and right undercut by 5°and 8°, respectively, the rock deformation is reduced by about 4% (Fig. 9), and the range of final subsidence is reduced by about 12%. It can be concluded that the existence of the additional vertical joints reduces the importance of the steeply and gently dipping joints setting, which plays a leading role in the development of the surface subsidence.
In order to study the symmetry relationship of surface subsidence caused by block caving, symmetry index (SI) is introduced. The index is defined as the ratio of the minimum to maximum subsidence angles on both sides of the undercut. Absolute symmetry corresponds to SI of 1. The symmetry index of the five-group models is shown in Fig. 8, and the highest symmetry index is 0.987 in A1, followed by 0.974 in A5 and lowest at 0.930 in A3. It can be found the symmetry index decreases gradually in the three groups from A1 to A3, suggesting that the rotation pattern of the joints has an important influence on the asymmetry development of surface subsidence. For instance, the extent of asymmetry increases with the increase of the joints rotation angle. The difference is that the symmetry index in A4 and A5 is higher than A3, indicating that reducing the inclination of the gently dipping joints can limit the development rate of subsidence on the right undercut and the additional vertical joints are more favorable for the symmetrical development of surface subsidence on both sides. As shown in Figs. 10 and 11, surface deformation is obvious with the change of joints inclination, the overall deformation of the major surface vertical and horizontal displacements in A3 reaches a maximum of 272 m and 296 m, respectively. For all models, the total extent of the major surface horizontal displacements is mainly larger than of the vertical displacements, indicating the collapse of the near-surface fractured rock mass is mainly in a lateral direction, which can be attributed to the topping failure mechanism of rock mass along the joints.

Impact of Faults
Combining the rock mechanics model provided by CEMI (Fig. 3), the faults influence on the surface subsidence is systematically studied. Two major faults are provided in the mechanical model, namely the North fault (79°) and the South fault (90°) along the A-B section. In the two-dimensional numerical model, faults are defined as severely weak interfaces.
This study investigates the effect of faults on surface subsidence deformation through multiple sets of models, which involves north fault with different locations and inclination. Randomly generated joints are included in the model, and the geometric model under different fault conditions is shown in Fig. 12. B1 and B2 models are used to illustrate the faults effect on the surface subsidence. B3 and B4 models based on the B2 model are utilized to evaluate the effect of location of the north fault on the surface subsidence. The north fault is moved forward by 50 m in B3 and moved back by 50 m in B4. B5 and B6 models based on the B4 are applied model to evaluate the effect of the north fault inclination on the surface subsidence, in which the north fault inclination is 59° (reduced by 20°) in B5 and 99° (increased by 20°) in B6.
The subsidence angle and RDR in different fault conditions are shown in Fig. 13. By comparing B1 and B2 models, the main fracture mechanism of the rock mass is shear and tensile failure along the preexisting joints and the complete rock bridge without faults (B1). The development of broken rock mass is mainly distributed between the north and south faults with regard to faults (B2), and the tensile and shear stress in the rock mass cannot be effectively propagated to the surrounding rock mass beyond the faults. Only the main surface subsidence occurs between the two faults, no obvious fissure outside the fault was found. Compared with the B1 model, the subsidence angle in B2 model increases by 6° on the left side and 4° on the right side, and the RDR is reduced by 4% (Fig. 14). As a result, the fault has an important effect on the fracture mechanism and the development of the surface subsidence. In addition, the influence is greater than joints. Based on subsidence angles under different fault conditions, it can conclude that the subsidence deformation is basically sliding along the south fault at the right side and the subsidence angles are consistent at 78° (4° higher than B1). For B3 and B4 models, significant subsidence development mechanisms were observed by changing the position of the north fault. When the north fault is moved forward by 50 m (B3), namely the north fault is located inside the excavation boundary, the footwall rock mass with caving that induced unloading could triggers translational failure along the fault interface. Eventually, it leads to the hanging wall collapse to the goaf, and only minor topping damage occurs near the fault interface. The subsidence angle is 95°, which is 11° higher than B2 model, and the symmetry index is only 0.821 indicating the surface subsidence is obviously asymmetric in the direction towards the excavation. The existence of fault effectively limits the destroying development of the hanging wall rock mass. When the north fault is moved back by 50 m (B4), namely the north fault is located outside the excavation boundary, there is no significant effect on the surface subsidence, as shown in Fig. 13. At this time, the subsidence angle is basically the same as that without fault (B1). Compared with the B2 model, the subsidence angle is reduced by 7° and RDR is increased by 3%, indicating the extent of impact on the surface subsidence area is minimal as the fault is far from the excavation area to a certain distance. B5 and B6 models are based on the B3 model, the fault is located inside the excavation boundary, by changing the north fault inclination, different subsidence development mechanisms are also observed. For 59° fault model (B5), the rock mass is mainly sliding failure along the footwall, the existing fault does not play a key role in limiting the rock mass movement in the hanging wall. The subsidence angle decreased from 95° to 86°, and RDR increased by 5%. Tensile fracturing can be noted in the hanging wall during the vicinity of the caving boundary, indicating the weakening effect of the fault on the hanging wall rock mass is enhanced with the reduction of fault inclination, and also resulting in large-scale topping failure. This result is consistent with the hanging wall failures proposed by Brown and Ferguson (1979) . For the 99° fault model (B6), the hanging wall and footwall are changed. It is found that the hanging wall rock mass is mainly sliding failure along the fault interface indicating the failure is mainly controlled by the plane of weakness provided by the fault, continuous ore extraction leads to the hanging wall rock mass breaking into segments and sliding to the cave. The footwall rock mass is only partial failure, and this failure is dominated by randomly distributed joints. The subsidence angle is only 3° higher than B5 model and RDR decreased by 2% with increasing the symmetry of surface subsidence (relative to B3 model). Figures 15 and 16 show the results of horizontal and vertical surface displacement (≥ 5 cm) with different fault conditions. Taking B2 model as a reference without fault (B1 model), the main horizontal and vertical displacement of the overall range increased by about 20% and 14%, respectively. It is suggesting that the subsidence range is mainly controlled by faults when faults located in the primary impact area. For the B3 model, the subsidence area is reduced by about 15%, and the minimum displacement is obtained. It is interesting to note that the surface horizontal displacement is basically the same as the vertical displacement. For the B4 model, the main horizontal displacement and vertical displacement of the overall range increased by about 11% and 5%, respectively, and the subsidence angle is basically consistent with the B1 model, which indicates the fault away from the main affected excavation area at a certain distance will reduce the impact on the surface subsidence displacement.

B1
Taking B3 model as a reference, B5 and B6 models illustrate the effect of fault inclination on the range of surface subsidence. When fault inclination is 59° (B5), the horizontal and vertical displacement increased by 7% and 11%, respectively. It indicates the weakening effect of the hanging wall rock mass is increased as the fault inclination is slowed (reduced by 20°), providing the conditions for the unstable development of the rock mass along the fault plane and causing large-scale topping failure. This is the main reason for the overall surface subsidence increasing. When the fault inclination is 99° (B6), the main horizontal and vertical displacement increased by 5% and 9%, and the overall subsidence range is reduced relative to the B5 model. This demonstrates the footwall movement was more limited in exposure as the angle of the fault inclination increased by 20°. The passive support provided by the broken rock section prevent the development of major internal instability, which lead to the overall sliding failure of the hanging wall rock mass along the fault interface. Nevertheless the footwall rock mass is only damaged and form a smaller subsidence range. It is worthy of note that removal of this support rock mass will cause further failure, which is the main reason for surface collapse formation and expansion.

Conclusions
(1) Crack formation, propagation, coalescence and the stress-balancing arch evolution can be clearly demonstrated during the whole caving process.
The numerical results show that the stress-balancing arch plays a key role in the development of rock mass fracture. The rock mass fails to cave to the surface until all the stress-balancing arches disappear, and the surrounding collapse is usually characterized by intermittent failure. (2) Without the faults, the steeply (vertical) joints control the main development mechanism of surface subsidence. The shallower inclination corresponds to the larger subsidence range, and the worse symmetry of the surface deformation. When the steeply inclined joints are invariant, reducing the inclination of the gently dipping joints can limit the surface subsidence in a certain range. When the steeply joints possess different inclinations coexist, the steepest joints have the dominant effect on the development of subsidence, and the collapse of fractured rock mass near the surface is mainly in the lateral direction.
(3) When the fault located in the main influence area, it can weaken the influence of joints on the surface subsidence. The fracture development of the rock mass is mainly concentrated on area between the faults. The fault partially intersecting the caving area controls the extent of surface subsidence deformations. Gently dipping fault will increase while steeply dipping fault will reduce the failure extent of the hanging wall rock mass. The gently inclined fault will result in larger failure extent of the overlying rock mass. (4) Understanding the crack evolution process and subsidence mechanism during the block caving is an essential aspect for accurate prediction of the rock fracture development and the surface subsidence range, which can timely adjust the mining order and put forward preventive measures. Moreover, it provides the effective guarantee for mine safety exploitation, surface environment and structural protection.