Subsurface strata failure and movement based for improving gas emission control: Model development and application

In underground mining operations, the overlying strata move, deform, and crack as the underground working face advances due to strata breakage, failure, and subsidence, forming a “vertical three zone” over the mined out area. The permeability of strata under mining changes differently in different zones. Therefore, studying the distribution of strata permeability is very important to obtain knowledge of gas migration pattern in the overlying strata in order to increase the gas extraction efficiency. In this study, a series of numerical models are developed to describe the overlying strata movement for three different zones. In detail, models for calculating the strata subsidence for fracture and bending zones and analyzing strata collapsing development for the caving zone are proposed based on the influence function method and mechanical analyses. In addition, on the basis of such models, a new comprehensive model to describe the strain‐porosity‐permeability of strata is developed. In the case study of applying such models in an underground mine, the reliability of the strata movement model is verified and the permeability change is also calculated, which gives a good guidance to design the location of gas drainage roadway practices in the future. The field engineering practice shows its good application.


| Strata Movement Model in fracture zone and bending zone
The Knothe's influence function is the most successful one that has been adapted by many major coal-producing countries. It states that the distribution of the subsidence caused by the extraction of one unity area can be expressed by a modified normal probability distribution function. The influence functions for subsidence and horizontal displacement in two-dimensional case (good for predicting final subsidence horizontal displacement along major cross-sections) are shown in Equation (1) and (2). 37 where m is the mining height, a is the subsidence factor, R is the radius of major influence, x' is the horizontal distance between the extracted element and the surface point where final subsidence is to be calculated and H is the mining depth.
The derived mathematical expressions for the final subsidence and horizontal displacement at the surface point x are where d is the offset of inflection point and W is the longwall panel width.
To broadly apply the models into subsurface strata for predicting the subsurface strata movement, significant modifications should be made. The proposed mathematical model for predicting subsurface strata movement would address the mentioned problems. 38 The subsurface strata movement can be determined as the following procedure: Collecting the geologic columnar profile over the underground longwall panel and naming the stratum from immediate roof to the surface as 1st, 2nd, 3rd,…, nth rock layer as Figure 1.
Defining the subsurface influence function. It should be noted that the influence functions defined here for the vertical and horizontal movements, which are applied in studying subsurface strata, would generate different results, such as inverted parabolas profiles that represent the subsurface strata movement and horizontal movement, at different rock layers. The following Equation 39 is used for subsurface subsidence along a major cross-section: In this equation, the coordinate x ′ is the horizontal distance between the point of "influence" to cause subsidence and the prediction point on the top surface of the layer. The term S max (z i−1 ) is the predicted maximum final subsidence of the past layer z i−1 , which can be obtained from the calculation results for the last layer. a i and R i are the subsurface strata movement factor and radius of major influence for the ith layer. For the first layer immediately above the mined coal seam, the mining height should be used in the place of S max (z i−1 ) in the influence function. 40 The influence function for horizontal displacement along a major cross-section is derived from the influence function for subsidence (Equation 6) as shown in the following equation.
Integrating the influence function within a proper horizontal interval for the final subsurface subsidence and horizontal displacement. 41 Similar to calculate the surface deformation using the influence method, the final subsurface subsidence and horizontal displacement at a prediction point (x,z i ) are obtained by integrating influence functions between left and right inflection points at the corresponding rock layer as shown in Equations (7) and (8). In equations, W is the width of the mined longwall panel and d i is the offset of inflection point.

| Calculation of extreme spans for initial breakage and periodic breakage of the strata
Two parameters, the extreme span (L j ) of the initial caving of the strata and the extreme span (L zj ) of the periodic caving of the strata, are required to determine the initial or the periodic caving, which occurs once the length of the suspending strata is longer than the corresponding extreme span. 42

Extreme span of initial caving of the strata
The strata break down when its suspension length reaches the extreme span and then accumulate in the gob. Then, the extreme span of the fracture failure can be inferred by analyzing the supporting mechanical model of strata. In practical situations, the stress state of the strata is very complicated. In order to simplify the analysis process, the strata can be simplified to clamped beam for calculation and be simplified to the cantilever beam after initial breakage for continuing the calculation related to the periodic breakage.
According to the rock mechanics theory, the Mohr-Coulomb yield criterion is usually used in the engineering calculation as the criterion for distinguishing shearing and breakage. The Griffith criterion is used as the criterion for the rock strength of the tension stress damage. 43 These theories can be similarly used as the criterion for distinguishing the breakage of the overlying strata breakage under mining. Before the initial breakage of the overlying strata under mining, the strata can be considered as the clamped beam by means of elasticity and the stress state of the clamped beam can be obtained as shown in Figure 2. The horizontal stress of clamped beam in coordinates is calculated according to Equation (9) 44 : In the Eq.: l-suspension length of the strata, m; h-thickness of the strata, m; -rock Poisson's ratio; q-overlying load of the strata, Pa.
When x = ±1/2 and y = h/2, the horizontal stress x can be calculated according to Equation (10): When x exceeds the tensile strength [ t ] of the strata, the strata damage. The extreme span of the initial breakage of the strata is calculated according to Equation (11): In the Equation (11), L j indicates the extreme span of the strata if the breakage is in the form of the clamped beam. Once the suspension length of the strata reaches or exceeds such value, the overlying strata suffer the initial breakage.

Extreme span for initial breakage of the strata
When the strata suffer the initial breakage, as the underground working face continues to advance, the strata on one side of working face can be considered as a cantilever beam. The strata similar to the cantilever beam in the ground pressure control may suffer periodic breakage, thus to form the periodic pressure. When the cantilever beam bears deadweight and uniform load, horizontal stress x , vertical stress y , and total stress xy are calculated according to Equation (12), (13), and (14), respectively: x exceeds the tensile strength [ t ] of the strata, the strata damage and the extreme span of cantilever beam strata is shown in Equation (15): In the Equation (15), L xj indicates the extreme span of the strata when the breakage is in the form of the cantilever beam. Once the mining length of working face reaches or exceeds L j + L xj , the overlying strata start to suffer the periodic breakage.

Movement of Breakage Strata
After the initial breakage of the strata, as the mining continues to advance, the strata on one side of the working face suffer the continuous periodic breakage in the form of the cantilever beam and the strata on the side of open-off cut are considered as the cantilever beam. Therefore, the strata above it are suffering the periodic breakage, and the state of the cantilever beam strata on the side of open-off cut keeps relatively stable.
The bended and deformed curve of the strata on one side of the working face can be calculated by employing the deflection curve equation of the cantilever beam under the uniform load in elasticity and material mechanics, as shown in Equation (16): In the Equation (16), the value range of l is (0, L xj ). As the working face advances and the suspension length l of the cantilever beam strata increases, once l is equal to L xj , it shows the strata complete a periodic breakage and start a new periodic movement. Therefore, the horizontal length of subsidence curve of certain cantilever beam strata always varies as the working face advances, as immediate roof, main roof, and strata region in the third layer of overlying strata in dashed line frame one side of the working face as shown in Figure 3. The strata on one side of open-off cut are considered as the cantilever beam, as strata region in left dashed line frame. However, suspension lengths of immediate roof, main roof, and the third layer of strata are fixed. Therefore, it can be considered that the cantilever beam strata on one side of the open-off cut are consecutive after the strata suffer breakage and the calculation results by strata movement model established based on influence function method are followed. Because the state of the cantilever beam strata on one side of open-off cut is relatively stable, the horizontal distance from the right end of the cantilever beam on the left on the n th layer to the open-off cut can be calculated by the following Equation (17): For broken rocks filling most space of the gob, the rock may suffer the volume expansion compared with the in situ stress condition, namely the embodiment of the broken expansion. This characteristic is one of the important causes for preventing the continuous breakage of the strata in the fractured zone. Therefore, after breakage and accumulation of the strata, the residual broken expansion coefficient K is employed to calculate subsidence and horizontal displacement of the strata. Without considering the expansion of the strata in the tendency, subsidence and horizontal displacement of caving and breakage part of the strata on the ith layer can be calculated by the following Equations (18) and (19): In the Equations: m-mining height, m; K -residual broken expansion coefficient of rock, dimensionless; h-thickness of the strata, m.

| Calculation process of Strata movement model
In the advancing process of the working face, the movement of the overlying strata is a dynamic process. The roof strata suffer the initial breakage firstly, and the roof strata on the mining side of the working face form the cantilever beam as the advancing distance W of the working face increases. After the immediate roof (the first layer of the strata) suffers the breakage, the upper layer of the strata (the second layer of the strata) still bears a certain load. Therefore, the second layer of the strata can be considered as the clamped beam with the continuous initial breakage of the strata until the end of the mining. As the working face advances, the cantilever beam of the strata on one side of the working face continuously suffers the breakage. The upper strata may suffer the breakage after the suspension length of the clamped beam or the cantilever beam exceeds the extreme span.
The continuity of the numerical model established based on influence function method is not suitable for calculating the movement of the broken strata. Therefore, conditions for the breakage of the strata should be discussed. Figure 4 shows the judgment of a point (x, z i ) on the ith layer of the strata and the chart of calculation process.

| The Study of permeability distribution in strata of overlying vertical three zones
The rock permeability is the main parameter characterizing the rock permeability. 45 When the overlying strata are affected by mining, factors such as in situ stress and temperature change to different degrees, while the permeability of the strata may vary with these factors. Generally speaking, when the coal seam is mined, the overlying strata are disturbed and broken and bended to different degrees forming "vertical three zone" from underground to surface, as shown in Figure 5: In order to more scientifically reveal the movement of the overlying strata with the mining coal seam as well as the produced fracture distribution and development and to finally provide theoretical and technical basis for underground gas control, the established strata movement model is combined in this paper to study on the change law of the permeability with the rock stress after the rock strata subside and move.

| Model of Strain-Porosity-Permeability change in Overlying Strata
In the mining process of coal seam, the emission of gas causes extreme threat to the mine safety production. The coal rib, the broken coal pieces, and the gob are three sources of gas emission, in which gob gas contributes a relatively high proportion of gas emission. The caving zone mainly focuses on the pore permeability, fractured zone, and bending zone focus on fracture permeability.
The fractured in the overlying strata under mining can be studied according to characteristics of overlying strata zoning: (a) caving zone: It is the porosity-fractured zone; (b) fractured zone and bending zone: They are main development regions of fractures. The value solution method and influence parameters of the permeability of the fractured field in each zone are different.

| Determination of permeability of porous coal of caving zone
The caving zone is mainly formed by accumulation of residual coal and broken rocks, and the diffusion method of gas is mainly the movement in the form of laminar and diffusion. The permeability of caving zone can be expressed as Equation (20): 46 where K k -permeability of caving zone, m 2 ; -porosity, %; K 0 -original permeability of the strata, m 2 .
The porosity of the strata can be calculated according to the following Equation (21): where 0 -strata original porosity of the strata, %; t -total strain of strata, mm/m.
The total strain t of the strata can be calculated by the horizontal stress x and the vertical stress z produced by the strata, and Equation (22) and Equation (23) are obtained by taking the derivative of the Eq.: The total strain t is shown in Equation (24): The strata may suffer breakage when the suspension length reaches the extreme span, and the broken rocks may accumulate in the gob, mainly focusing on the porosity development. Equation (24) is substituted into Equation (20) to obtain Equation (25) for calculating the permeability of the caving zone:

| Determination of fracture permeability in fractured zone and bending zone
The fractures with different widths are distributed in the fractured zone of the overlying strata. For fracture and length in the rock, the permeability can be expressed as Equation (26) by using the parallel plate model, 47 in which the fracture width b is calculated by the tension fracture, as shown in Equation (27), 48 and the fracture depth Z involved in the calculation process can be calculated according to Equation (28): where K ls -permeability of fracture per unit area and length, m 2 ; f -fluid density, kg/m 3 ; g-gravitational acceleration, m/ s 2 ; b-fracture width, m; C s -proportional coefficient, usually 0.003~0.015; Z-fracture depth, m; E-elastic modulus, MPa; x -horizontal strain, mm/m, which can be calculated according to Eq. (22); r -bulk density of the strata, kN/m 3 ; v-Poisson's ratio.
Equation (22), Equation (27), and Equation (28) are substituted into Equation (26) to obtain the Equation (29) for calculating the permeability of the unbroken strata in fractured zone and bending zone:

| Calculation process of "verticalthree-zone" permeability model
In continuous mining, the overlying strata of the working face suffer initial breakage and periodic breakage to form caving zone, above which are fractured zone and bending zone in turn. The calculation process of strain-porosity-permeability model established aiming at "vertical three zone" is shown in Figure 6, and a point (x, z i ) on the ith layer of the strata is taken for an example. The involved parameters, such as the horizontal strain x of strata and the total strain t of strata, should be obtained by taking the derivative of Eqs for calculating horizontal and vertical movement of the strata. Therefore, the permeability of overlying strata under mining should be associated with the displacement, strain, and porosity of the overlying strata under mining.

| Background of Working Face
For mining of #15 coal seam in a mine, the width of the working face is 230 m and the length is 1584 m. The average thickness of coal seam is 6.6 m. The maximum mining depth is 526 m, and the average mining depth is 468 m. Figure 7 shows part of physical parameters of comprehensive geological histogram #15 coal seam and the corresponding strata, and all studied overlying strata lithology and relevant physical parameters are shown in Table 1.

| Comparative analysis of calculation
results of strata movement model and 3DEC numerical simulation results

3DEC modeling
The 3DEC simulation software is a simulation software based on discrete element method, which can be used to study the movement of the strata bearing static load or dynamic load under noncontinuum. 49 Before establishing the overlying strata numerical model, main steps include the following: Selecting correct and appropriate constitutive model, which can ideally reflect the mechanical behavior of strata under mining. 33 3DEC simulation software provides three constitutive models aiming at movement of blocks, and Mohr-Coulomb constitutive model is usually selected as the constitutive model for simulation of strata movement. Because the rock is considered as block in the simulation process, Mohr-Coulomb constitutive model is selected as constitutive model complying with overlying strata of working face for resolution.
Preliminary research and data collection. The main purpose is to select reasonable physical-mechanical parameters. Mohr-Coulomb constitutive model involves massive basic physical parameters, which should be set and perfected in advance, such as rock density, elastic modulus, Poisson's ratio, internal friction angle, cohesion, and tensile strength. Parameters of the strata are determined and shown in Table 2. Once determined the specific value of physical parameters of strata, the physical parameters of contacts could be obtained from the empirical value based on whole strata's in-stope physical parameters, and parameters of the interface between blocks are shown in Table 3 after the division of blocks.
Setting of boundary condition. According to the stress environment of the strata, appropriate stress and speed boundary conditions are selected. The height of the strata involved when simulating the movement of overlying strata of the working face is 74.32 m, and 16 MPa of pressure is loaded on the top layer of the strata surface.

| Comparative analysis of calculation results of strata movement model and simulation results
In order to study the development of porosity and fracture in the overlying strata over the working face, the movement of 12 layers of strata in 76.14 m range over the working face is calculated. It is very important to judge whether the strata suffer the initial breakage as well as the periodic breakage. The extreme span for the initial strata breakage, the periodic breakage span, and the residual broken expansion coefficient K of the corresponding strata blocks after the breakage are shown in Table 4. One of the very important factors for the suspended subsidence and the movement of the overlying strata in the gob is the broken expansion of blocks in the caving zone. When the rock suffers the breakage after liberated under the in situ stress, the volume of broken rocks may suffer a certain expansion. When the volume increment after the expansion of broken and accumulated rock completely fills the volume of the mined out space, the strata above the caving zone at this time will not continue the breakage. Figures 8 and 9 show breakage times, predicted movement curve, and simulated movement of strata of overlying strata when the working face advances to 50 and 150 m. The schematic of the breakage of the overlying strata shows the initial breakage and the periodic breakage of the strata, and the 3DEC numerical simulation shows the trend of the strata with the movement of the working face in the mining process. It can be Therefore, it can be considered that the volume of breakage blocks in #1~#8 stratum after expansion has filled the gob volume caused by mining. The strata above #8 strata will no longer has fierce breakage phenomenon. In the whole advancing process of the working face (advancing length is between 150 m and 1584 m), the height of overlying strata of the working face in gob in caving zone is #1~#8 stratum, 55.37 m of height from the coal seam.
(3) Figure 8 shows breakage law of overlying strata, movement curve, and simulated strata when the working face advances to 50 m. Upon calculation, the length of initial breakage of coal seam immediate roof ( #1 stratum) is 53.86 m; therefore, the strata are only bended and deformed but not broken, and the movement of the whole overlying strata is only limited to the bended and deformed movement. Therefore, fractured zone and bending zone are only formed on the whole overlying strata. The movement prediction curve of the strata shows that the subsidence and movement curve of coal seam immediate roof have obvious flat-bottom phenomenon. cut has almost breakage sign, which reflects 1# strata will suffer initial breakage. The calculated subsidence quantity of 1#~12# strata of the whole overlying strata reaches 1.2~6.6 m, and the radius of the subsidence basin of the strata increases as the vertical distance of coal and rock increases.
When the working face advances to 150 m, 1#~8# strata suffer the breakage, and the height of broken rocks can balance the coal seam thickness after the broken expansion. Therefore, the maximum height of the development of the caving zone is 55.37 m, while simulation results show the development height of the caving zone reaches 57.3 m, as shown in Figure 9(3). Three obvious separation fractures occur in the overlying strata on one side near the working face. Figure 9(1) shows that 8# strata suffer two times of periodic breakage; at the same time, 3# strata suffer up to 10 times of periodic breakage, while 6# strata suffer one time of periodic breakage, and other strata, respectively, suffer 5~7 times of periodic breakage. The cantilever beam strata on one side of the working face after the breakage of the strata are shown in broken circle in Figure 9(2), and lengths of the cantilever beam are, respectively, 0.5 m (1#), 0.03 m (2#), 6.54 m (3#), 9.08 m (4#), 2.38 m (6#), and 1.96 m (8#). Figure 9(3) shows the result of simulated strata movement. After 8# strata suffer breakage, the height of the caving zone continues to develop compared with the working face advancing to 100 m, the height (55.42 m) of subsidence prediction curve of 8# strata in the gob is almost near the original height (55.37 m), namely, rocks accumulate after the breakage of 1#~8# strata, because the height of the rock after the volume expansion reaches the height of 8# strata, the strata above 8# strata will T A B L E 5 The comparative analysis of strata subsidence using 3DEC simulation and mathematical model calculation when the working face advancing distance varying from 0 meter to 200 meters not have obvious breakage phenomenon because of the support of broken rocks in the process of bending and sinking. Four monitoring lines are taken along the trend in the operational process of 3DEC simulation and numerical model, respectively, 25, 75, 125, and 175 m from the open-off cut. Nine subsidence displacement monitoring points are taken according to strata stratification of each monitoring line vertically, located in the middle of the strata. As the working face advances, the simulation and comparison of simulation calculation results are shown in Table 5. The comparative analysis of the difference in both is shown in Figure 10, and the basic trend of 3DEC simulation results and model calculation results is consistent in the subsidence process of overlying strata. When the working face advances to 50m, the error between simulation results and model calculation results is large, which corresponds to 50m (No.1) in Figure 10 Combining the existing subsidence model with strain-porosity-permeability relation employed, the permeability distribution model of overlying strata in gob is obtained. To obtain distribution prediction of the permeability K of the overlying strata under advancing length condition of different working faces as shown in Figures 11 and 12, which involves part of strata, parameters (bulk density r , Poisson's ratio v, elastic modulus E) can refer to data of 1#~12# strata in Table 2. The original porosity 0 of the strata and original porosity K 0 is shown in Table 6.
According to relevant research findings of petrophysical phase of gas reservoir 50 and the permeability of coal reservoir described by well log, 51 the strata are divided into three different levels and types according to the permeability of the strata, namely I-level highly permeable area, the permeability is more than 1 md; II-level moderately permeable area, the permeability is between 0.1 and 1 md; and III-level lowly permeable area, the permeability is less than 0.1 md. When the advancing distance of the working face is 50 m, the mining effect has the small disturbance to the overlying strata. Therefore, all the overlying strata are considered being in fractured zone or bending zone and the variation range of the permeability is small, approx. 0.6~1.7 md. Only the permeability of part strata near the coal seam significantly changes, as the strata area with distances from A, B, C, and D less than 32.1 m from coal seam shown in Figure 11. The strata below 5# strata, in which the area (B, C) of the blue isoline is the area where the permeability decreases compared with the original rock permeability, which is usually less than the original permeability of the corresponding strata and the range is 0.6~1 md belonging to II-level moderately permeable area. The strata area (A, D) of red isoline far away from the middle of gob is the area where the rock increases. The range is 1~1.7 md, belonging to I-level highly permeable area. The permeability of other strata above 5# strata increases or decreases; however, the variation range of the permeability and the affected rock range are small. A', B', and C' in Figure 11 are, respectively, areas where 6#~12# strata are disturbed, and the permeability of the strata in this area slightly increases or decreases. When the working face advances to 150m, the change distribution of permeability K of overlying strata in gob is shown in Figure 12. At this time, the advancing length of such working face just reaches the length (149.46 m) of completed development of overlying strata in caving zone, and the final development height in caving zone is 55.37 m, namely the accumulated thickness of 1#~8# strata. From Figure 12, the maximum area of permeability K is 1#~4# strata (B) near coal seam, and the next is 5#~7# strata (B'), which increases by 12.45 m in the height of area B' in Figure 11, and the permeability of 8# and 9# strata increases to a certain extent. The maximum permeability of 1#~3# strata is 5.4 md, and B belongs to I-level highly permeable area; the maximum permeability of 4#~7# strata is approx. 3.2 md, and B' belongs to II-level moderately permeable area, the permeability is between 0.4 and 1 md.

| Design of high-level drainage roadway of working face and assessment of effects
The permeability of coal and surrounding rock of the working face is small, and emission quantity of gas is nearly from the coal seam before breakage of the strata in gob and the relative gas emission quantity is below 10 m 3 /t. However, after strata breakage, gas of surrounding rock emits largely to the working face of coal mining with the mining fracture. The absolute gas emission quantity of the mine reaches 190 m 3 / min, in which over 90% of gas comes from the adjacent coal seam and surrounding strata.
The average thickness of #15 coal seam of working face is 6.6 m, and the height of fractured zone locates between 39 and 230 m over #15 coal seam and the high-level drainage roadway of working face is finally designed at the position at 52 m from #15 coal seam over the working face. Figure 13 shows a plan view for arrangement of high-level drainage roadway of working face. The gas tube pipeline arranged in high-level drainage roadway is connected to ground gas extraction pump station. The horizontal distance from high-level drainage roadway to return airway of working face is 20 m, and the vertical distance from the working face is 52 m. Fractures develop in the "O" form range over the longwall panel. Such areas are good for doing gas extraction. The forming negative pressure by the gas pump at the upper corner of the working face could effectively extract and gather mine gas in strata and around working face.
According to the observation daily report of high-level drainage roadway of the working face, as shown in Table 7, it shows that gas flowrate drained by the return airway and tail roadway is small. The maximum gas flowrate in the return airway is 1.14m 3 /min, and the gas flowrate of tail roadway is 0.08~0.1 m 3 /min; the gas flowrate of high-level drainage roadway takes a large proportion in total emission of gas and extraction flowrate reaches 39.2~45.85 m 3 /min, with the proportion up to above 96.58%. Thus, it can be seen that highlevel drainage roadway of working face is the most important gas extraction way of the working face, which can effectively control the gas concentration of the working face.

| CONCLUSIONS
Because the overlying strata over the mine gob may form different movement states of caving zone, fractured zone, and bending zone after affected by the mining, different calculation models for strata movement are established. Reliability validation is carried out by using numerical simulation method. The strain-porosity-permeability relation of the strata is studied, and different basic equations of the calculation permeability are employed aiming at the "vertical three zone" of overlying strata under mining. The prediction analysis model of strainporosity-permeability of the strata is developed. It is applied to an actual working face, which shows the variation situation of the permeability after deformation and destruction of upper overlying strata after mining. Through the engineering practice of designing the high-level drainage roadway of a working face, the gas extraction efficiency is improved, which also show developed models with the good application value.