#### Abstract

The evolution of the fracture network induced by mining has an important influence on the mechanical behavior of the rock and the safety of the mine. In this study, a new 3D modeling approach based on discrete fracture network (DFN) and fractal theory was developed and applied to the Tashan coal mine. The model results of the evolution of fracture in the coal seam are consistent with field observations. The evolution law of fracture network in the overburden strata was studied, and the results show that the fractal dimension of the fractures in the lower strata increased linearly and stabilized quickly within 50 m behind the mining face. In the higher strata, most of the fractures were generated behind the mining face and continued to develop farther than 100 m. More fractures were generated in the lower strata than in the higher strata, and the fractures more easily developed and expanded in the soft rock than in the hard rock. The evolution of the fractures of the main thick hard roof in the lower strata had a great impact on the generation of fractures in the higher strata.

#### 1. Introduction

The longwall top coal caving (LTCC) method is a new type of mining technique that is popular for thick coal seams in China due to its high efficiency and output [1]. However, these large-scale excavations seriously disturb the original equilibrium state underground, and the mining-induced redistribution of stress leads to the initiation and propagation of fractures, creating fractured zones with high connectivity around the mining excavations. These fractured zones weaken the strength of the rock mass and provide a pathway for confined water and gas inflow, leading to the movement of overburden strata and water and gas outbursts, which severely affect mine safety [2–6]. It is an essential issue in mining science to recognize the evolution of fracture networks during mining processes. Many researchers have studied this topic and achieved fruitful results.

Xu and Qian [7, 8] studied mining-induced fractures in the overlying rock layer of a coal mine and revealed their “o”-shaped distribution in the longwall face. Li et al. [9] discussed the effect of key strata on the distribution of fractures and concluded that an elliptic zone would be formed in the overlying strata. Kang et al. investigated the mechanism of the occurrence of mining-induced fractures and analyzed the morphological characteristics of a fracture zone [10–12].

Nevertheless, it is difficult to develop potential prevention mechanisms and the associated spatiotemporal effects throughout the mining process from a quantitative perspective due to the extreme complexity and disorder of mining-induced fractures [13]. The fractal geometry proposed by Mandelbrot [14] can quantitatively describe complex objects and irregular phenomena. Xie et al. [15] introduced the fractal geometry of fractures and discovered the self-similarity behavior of the spatial distribution of fracture networks. Since then, the evolution of fracture networks has been studied with fractal theory based on similar material models [16, 17], on-site field observations [18, 19], and 2D discrete element methods [20, 21]. However, field observations are limited by local conditions and observation tools. In similar material methods, the mining-induced fractures need to be determined manually from the model with subjectivity. 2D discrete element methods present the fracture evolution in only a single plane.

In this study, a new 3D modeling approach for the calculation of the fractal dimension of a fracture network was developed. This approach makes use of the built-in discrete fracture network (DFN) module in the commercial software FLAC3D to create stochastically distributed fractures. After each longwall face excavation, deterministic global stress analysis is first performed, followed by an evaluation of all the discrete fractures to identify fracture slip. A fracture is deemed to “generate” if it slips. Finally, the fractal dimension of the generated fracture network was calculated based on specified section planes. With the advance of mining, the evolution of the mining-induced fracture network can be determined.

This modeling approach was applied to the Tashan coal mine in Shanxi Province, China, where very thick coal seams were extracted using the LTCC method. Section 2 presents the detailed methodology of the numerical model. The establishment and running of the model based on the 8212 working panel are introduced in Section 3. The evolution of the mining-induced fracture network in the coal seam and overburden strata is discussed in Section 4.

#### 2. Description of the Modeling Procedure and Method

The modeling procedure involves three main steps. First, a 3D model representing the working longwall panel and the strata is built with a group of discrete fractures (a DFN), and then mechanical properties and initial conditions are applied to zone elements to reach an initial equilibrium state. Next, progressive longwall mining is simulated by assigning null mechanical properties to zone elements. At the next excavation step, together with the simulation of coal extraction, the previously “removed” elements are assigned goaf material to simulate the recompaction of the goaf. In the final modeling step, a fracture stress calculation is performed to determine whether fracture slip occurs based on the prevailing stress state. Specific section planes are set up, and the fractal dimensions of the slipped fracture network on the section planes are calculated. Steps 2 and 3 are repeated in each coal excavation step until the completion of the longwall mining simulation in the study area. Notably, the model progress from Step 2 to Step 3 is a one-way process, which means that the state of the fractures is determined by the zone element stress states, while it is assumed that the stress states of the rock mass are not affected by fracture slippage.

##### 2.1. DFN in FLAC3D

A randomly distributed network of disc-shaped discrete planar fractures can be created in the FLAC3D model by setting a series of attributes. These attributes include the fracture stopping criterion, fracture size distribution, and fracture orientation.

The total number of fractures was used as the stopping criterion to generate the fracture network. In consideration of the model size and computing ability, the total number of fractures was set to twice the total number of zone elements. The power law relationship [22, 23] has been adopted to describe the distribution of fracture sizes:where is the number of fractures with size per unit volume, is the density term, and *a* (>1) is the power law scaling exponent and usually ranges between 2 and 4.

The total number of fractures within the range [*l*_{min}, *l*_{max}] can be obtained by integrating equation (1):

The cumulative distribution function is then given by

In view of the meshed element size and the scale of research, the minimum fracture size was set to the shortest zone element length, and the maximum fracture size was set to 15 m.

The fracture position was assumed to follow a uniform spatial distribution throughout the model domain. Due to the lack of geological measurements, the dip angles and dip directions that determine the fracture orientations also follow a uniform distribution.

##### 2.2. Fracture Slip Criterion

The Mohr–Coulomb slip criterion is widely used to evaluate fracture slippage. The fracture is considered to have slipped if the shear stress along the fracture plane exceeds its resistance :where is determined by the properties of the fracture surfaces aswhere *σ* is the normal stress on the fracture surface and and are the static friction coefficient and the cohesion along the fracture, respectively.

Alternatively, the fracture shear strength in equation (5) can be estimated by applying a strength reduction factor to the intact rock strength given by the following equation:

The shear and normal stresses used in equations (5) and (6) can be calculated using the stresses on the zone element closest to the centroid of a fracture. Salamon [24] used this method for fracture slip evaluation. This approach is valid if the fracture is small compared to a zone element. As mentioned above, the fracture is larger than a zone element, so following Board [25] and Cao et al. [22], this study adopted nine test points for slip, one testing point at the center, and eight testing points uniformly distributed around the periphery of the fracture surface. The testing point is deemed to slip if the zone element closest to it is in the failure state or if equation (6) is satisfied, and the fracture is thought to slip when all nine testing points have slipped. The workflow of the evaluation of fracture slip is presented in Figure 1.

##### 2.3. Calculation of the Normal and Shear Stresses at a Test Point

The test points on a fracture surface are specified in local coordinates with their origin at (*x*_{0}, *y*_{0}, and *z*_{0}) in the global coordinate system (Figure 2). The direction cosines of the axis OX′ are *l*1, *m*_{1}, and *n*_{1}, and those of axes OY′ and OZ′ are *l*2, *m*2, and *n*2, and *l*3, *m*3, and *n*3, respectively. The transformation of the *i*th test point from its local coordinates to the global coordinates is given by

The normal stresses along the local coordinate axis directions at the test point can be written aswhere , , and and , , and are the normal and shear stresses of the zone element closest to the test point, respectively. The corresponding normal and shear stresses on the fracture surface associated with the test point are then given by

The computed normal and shear stresses are used in equations (5) and (6) to evaluate fracture slippage.

##### 2.4. Calculation of the Fractal Dimension of the Fracture Network

The calculation of fractal dimension on the section plane is done by the “counting the boxes” method [26], and the diagram of the calculation is shown in Figure 3. A square mesh of scale *R* is used to cover the target area, and the number of grid cells in which the fracture length is equal to or greater than the corresponding grid size is counted and denoted by *N* (*R*). An array of corresponding *N* (*R*) values can be obtained with a series of *R* values to study the changes in log (*R*) and log (1/*R*). The fractal dimension *D* is defined as follows:

**(a)**

**(b)**

**(c)**

#### 3. Establishment of the Numerical Model

##### 3.1. Introduction of 8212 Working Facing

The Tashan coal mine is one of the largest underground coal mines in the world and produced more than 20 Mt of coal in 2019. The 8212 working face is mined using the LTCC method and is located in the northwest region of the second mining district 2614 m along strike and 230 m along dip. The working face mines the very thick #3–5 coal seam, and the average thickness of this coal seam is 11.6 m. The dip angle of the coal seam is 4°; the coal seam can be viewed as an approximately horizontal structure. The average mining height and caving height are 3.5 m and 8.1 m, respectively. The main roof strata in 8212 working face are sandstone which would accumulate large elastic energy before failure [27, 28], and the breaking of roof strata can cause dynamic disasters such as rockburst [29–31] and coal-gas outburst [32, 33]. The detailed rock types and thicknesses of the 8212 working panel strata are presented in Table 1.

##### 3.2. Model Setup and DFN Generation in FLAC3D

A 3D numerical model based on the 8212 working face was constructed in FLAC3D. The size of the model was 600 m × 300 m × 500 m in the *x*, *y*, and *z* directions, respectively, and the mining area was 300 m × 230 m. The model comprised 620400 zone elements and 643720 grid-points. As discussed in Section 2.1, the number of fractures was set to twice the number of zone elements, i.e., 1240800. As the element size of the coal and roof strata was approximately 5 m × 5 m × 3 m, the fracture sizes, following a power law distribution, ranged from 3 m to 15 m, and the numerical model and the size of the generated fracture are shown in Figure 4. Notably, this paper focuses on the evolution of the fractures in only the coal seam and overburden strata, so discrete fractures were not generated in the floor strata.

**(a)**

**(b)**

**(c)**

##### 3.3. Property and Initial Conditions

The Mohr-Coulomb model was used to represent shear failure in the coal and strata. The detailed physical and mechanical properties of the strata at the 8212 working panel are presented in Table 2 [34]. The falling waste in the gob was considered an elastic material, and as the working face advances, the waste was progressively compressed by the overburden. Experimental results show that the density, Young’s modulus, and Poisson’s ratio can be characterized as follows [1]:

In the model, the horizontal displacement of the two sides was limited, and the vertical displacement of the bottom was limited. A vertical load () was applied on the top to simulate the overburden weight.

##### 3.4. Mining Excavation Procedure

After the initial equilibrium state, model mining started with the excavation of the intake airway and return airway, and thirty excavation steps were conducted to simulate progressive longwall mining. The excavation step is 10 m in the mining direction, representing two days of face advance. Each excavation also included the “backfill” of the previously removed elements by assigning waste material properties. A special subroutine was invoked after each step to determine and record the slippage of the fractures, and a fractal dimension calculation was implemented.

#### 4. Model Results

##### 4.1. Section Plane Setup

To describe the evolution of the fracture network with the mining process, eleven section planes are considered in the model, as shown in Figure 5. The first four vertical section planes are along the *X* direction (face advance direction) with coordinates of 250, 300, 350, and 400. Note that the coordinate in the *X* direction for the mining start line is 150, so as the mining face advances, and these vertical sections can record the evolution of the fracture network in front of and behind the mining face. Another seven horizontal section planes are placed at the center of each stratum to present the evolution of the fracture network in each stratum during the mining process.

##### 4.2. Evolution of the Coal Seam Fracture Network

As the mining face advances, the redistribution of the stress in the coal seam leads to the initiation and propagation of fractures. Figure 6 shows the evolution of fractures in the coal seam in section plane 3. Fractures continue to grow with the approach of the mining face. To quantitatively describe the evolution of the fractures in the coal seam, the fractal dimensions of the fractures in the section plane were calculated based on the method mentioned above, as shown in Figure 7. The fractures exhibited a two-stage evolution as the mining face progresses: a stable growth stage (at a distance of 170∼60 m from the mining face) and a rapid growth stage (at a distance of 0∼60 m from the mining face). In the first stage, fractures were generated under the influence of dynamic mining, and the number of fractures gradually increased as the mining face progressed. In the second stage, as the mining face continued to approach, the impact of mining became stronger, and the number of fractures increased rapidly. At approximately 10 m in front of the mining face, most of the fractures developed and interconnected due to failure of the coal body.

**(a)**

**(b)**

**(c)**

**(d)**

Field experiments, including borehole observation and abutment pressure monitoring, were conducted by Gao et al. [13] at the 8212 working panel, and the results are shown in Figure 8. The on-site field data present great fluctuations due to the changing geological environment. However, these data reveal a two-stage fracture evolution, and the abutment pressure also shows that the coal rock would break and fail at approximately 10 m in front of the mining face, leading to a decrease in abutment pressure. In general, the model results of the fracture evolution in the coal seam are consistent with the field observations.

**(a)**

**(b)**

##### 4.3. Evolution of the Fracture Network in the Overburden Strata during the Mining Process

Figure 9 shows the relationship between the fracture network and distance from the mining face in section plane 2. The different colors represent different strata. The fractures continue to develop and grow steadily with the advance of the mining face. The fractal dimension of the fractures in the vertical section plane with distance from the mining face is obtained to quantitatively analyze the fracture evolution, as shown in Figure 10. The development of the fractal dimension of fractures with distance from the mining face can be divided into two stages: a linear growth stage (in front of the mining face) and a slow growth stage (behind the mining face). In the first stage, as the mining face moved forward, the number of fractures increased steadily due to the influence of mining, and the fractal dimension increased linearly with the mining advance. Notably, these fractures were generated mainly in the coal seam and lower strata. In the second stage, the fractures in the lower strata continued to develop, leading to the deformation and movement of the lower strata. Meanwhile, fractures in the higher strata continued to develop, presenting a saddle shape in the vertical section plane. Therefore, the overall fractal dimension of the fractures in the second stage slowly increases and gradually stabilizes behind the mining face.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

##### 4.4. Evolution of the Fracture Network in the Overburden Strata in the Vertical Section Planes

Section 4.3 discusses the overall fracture network evolution with distance from the mining face. In this section, a detailed study of each overburden stratum is conducted. Taken from the vertical section planes, the development of the fracture network in each stratum with distance from the mining face can be obtained, and the fractal dimensions were calculated, as shown in Figure 11. The evolution of the fractures in each stratum is related to its location.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

In the lower strata such as the pebbly sandstone and siltstone, the evolution of the fractures is directly influenced by the mining excavation. The fractal dimensions of the fractures in the different section planes present similar characteristics. In front of the mining face, the number of fractures increased as the mining face advanced, and the fractal dimension of the fractures exhibited a corresponding linear increase. Behind the mining face, the strata collapsed rapidly and fill the goaf. Most of the fractures in the strata were well developed around the mining face, and the fractal dimension of the fractures stabilized quickly within 50 m behind the mining face.

The higher strata were not affected at the beginning of the mining process until the impact of the mining extended upward. Therefore, the evolutions of the fractures in each section plane are different. In addition, in contrast to the lower strata, in the higher strata, only a few fractures were generated in front of the mining face, and most of the fractures continued to develop more than 100 m behind the mining face.

##### 4.5. Evolution of the Fracture Network in the Overburden Strata in the Horizontal Section Planes

Figures 12 and 13 show the evolution of the fracture network in the typical lower rock strata (pebbly sandstone) and higher rock strata (medium sandstone) with mining advance based on the horizontal section planes. Fractures were generated around the mining excavations as the mining advanced, exhibiting an elliptical distribution in the strata. To quantitatively compare the evolution of the fracture network in each stratum, the fractal dimension of the fracture network was calculated based on the horizontal section planes, as shown in Figure 14. It should be note that the height “*h*” in the figure means the distance from the section planes to the coal seam. Two distinct features were presented from the figure:(1)In terms of the number of fractures, more fractures were generated in the lower strata than in the higher strata. In addition, a fracture could more easily develop and expand in the soft rock (mudstone) than in the hard rock (fine sandstone).(2)The higher strata (medium sandstone) were almost unaffected in the first 100 m of mining until massive fractures were generated in the main thick hard roof (siltstone), which shows that the evolution of fractures of the main thick hard roof in the lower strata had an important influence on the generation of fractures in the higher strata.

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

#### 5. Conclusions

It is an essential issue in mining science to recognize the evolution of fracture networks during mining. In this paper, a new 3D modeling approach based on a DFN in FLAC3D and fractal dimension was developed and applied to the Tashan coal mine. The model results of the fracture evolution in the coal seam are consistent with field experiments. Finally, the evolution of a fracture network in overburden strata was studied, and the following conclusions were drawn:(1)The fractures in the coal seam exhibited a two-stage evolution as the mining face progresses: a stable growth stage and a rapid growth stage.(2)In the lower strata, the fractal dimension of the fractures linearly increased as the mining face moved forward and stabilized quickly within 50 m behind the mining face.(3)In the higher strata, only a few fractures were generated in front of the mining face, and most of the fractures in the higher strata continued to develop further than 100 m behind the mining face.(4)More fractures were generated in the lower strata than in the higher strata, and the fractures were more easily developed and expanded in the soft rock than in the hard rock. In addition, the evolution of the fractures of the main thick hard roof in the lower strata had an important influence on the generation of fractures in the higher strata.

#### 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 they have no conflicts of interest.

#### Acknowledgments

The authors would like to acknowledge the support by the National Natural Science Foundation of China Youth Fund (grant no. 51822403).