A Gas Seepage Modeling Study for Mitigating Gas Accumulation Risk in Upper Protective Coal Seam Mining Process

Beijing Key Laboratory for Precise Mining of Intergrown Energy and Resources, China University of Mining & Technology (Beijing), Beijing 100083, China State Key Laboratory of Coal Resources and Safe Mining (CUMT), Xuzhou, Jiangsu 221116, China Key Laboratory of Deep Coal Resource Ministry of Education of China, School of Mines, China University of Mining & Technology, Xuzhou, Jiangsu 221116, China School of Physics, China University of Mining & Technology, Xuzhou, Jiangsu 221116, China Mining College, Guizhou University, Guiyang 550025, China


Introduction
Coal and gas outburst accident is a complex dynamic phenomenon, which occurs frequently in coal mines in China.Because of its sudden occurrence and intense destructiveness, the accident endangers both coal mine production and personnel safety.It also causes immense economic losses to mining companies [1][2][3].To address this issue, many research studies have been conducted nationally.These studies have elucidated the outburst mechanism in coal and gas; moreover, outburst risk indexes and mitigation measures have been proposed in these studies [4][5][6].Protective coal seam mining (PCSM) is preferentially used to release coal seam pressure and enhance permeability.It is one of the main outburst mitigation measures in China, mainly because gas content and its associated outburst risk are relatively low.By extracting the protective seam, the associated pressure is released and it substantially increases the permeability of the protected seam.Pressure relief gas extraction is enhanced consequently [7][8][9].During protective coal seam mining (PCSM), the phenomenon of methane concentration exceeding the limit usually occurs in the protective seam working face (PSWF).Because the protected seam offers gas pressure relief, a large amount of gas is desorbed and flows into PSWF.This is the main reason why methane concentration exceeds the limit in the PCSM process.To guarantee safe mining and to improve methane mining, pressure relief gas extraction technology is used currently [10][11][12].Presently, more attention is paid to pressure relief effect of the protected seam, but factors causing methane concentration to exceed the limit are ignored during PCSM [13,14].A qualitative method is mainly used to estimate pressure relief gas seepage, which is not applicable for gas drainage and control [15][16][17].Besides, there are few studies on the main influencing factors of pressure relief gas seepage.It is also not conducive to the extraction and control of pressure relief gas.Thus, a quantitative analysis of seepage characteristics of pressure relief gas is of great significance to the layout of gas drainage boreholes and the mining parameters of protective coal seam.
In this paper, previous studies on coal and rock permeability have been summarized.Stress-strain relationship, stress changes, and failure depth of the floor have been determined in the process of PCSM.Based on review results, GSM is established and its equation is used for calculating the quantity of gas desorbed from the protected seam into the protective seam during PCSM.The influencing factors were also identified by analyzing the seepage process.Finally, the model is validated with the data from an in situ ventilation air methane of the PSWF in a coal mine, which is located in Hancheng, China.The results indicate that the established GSM completely agreed with obtained field measurements.

Previous Studies on the Evolution Law of Permeability
Ever since the proposition of Darcy's law, rock permeability research has always been an attractive and ongoing research field.After years of development and research, two main types of the permeability model were established: porositypermeability model and stress-permeability model.These models are represented by Palmer-Mansoori's (P-M) and Gray's model, respectively [18].In the coal permeability model proposed by Gray [19], ground and expansion/shrinkage stress caused by adsorption/desorption gases was investigated comprehensively.
Here, σ e h is the effective horizontal stress, which is expressed in MPa; σ e h0 is the initial effective horizontal stress, which is expressed in MPa; Δp s is the change in adsorption equilibrium pressure, which is expressed in MPa; Δε s is the strain change in unit, which is caused by change in adsorption equilibrium pressure; v is the Poisson's ratio; p is gas pressure and p 0 is initial gas pressure, which are expressed in MPa; and E is the elastic modulus in GPa.
Herein, c f is the cleat volume compressibility with respect to changes in effective stress normal to the cleat, k 0 is initial cleat permeability at initial effective stress σ 0 , and k is cleat permeability at effective stress σ.
Equation ( 2) has been used in several laboratory experiments, which have been conducted to develop fitting regression equations.As shown in (3), Ren and Edwards [21] proposed that permeability and stress relation were associated with intact rock.k = k 0 e −0 25 σ−σ 0 3 Whittles et al. [22] summarized the relationship between the permeability of rock fracture surface and principal stress.Equation ( 4) is used for calculating fracture surface permeability.
Here, σ 1 and σ 3 are minimum and maximum confining stresses, respectively, which are expressed in MPa; k 1 is intrinsic permeability expressed in m 2 when σ 1 + σ 3 /2 = 1 MPa; and m is a material parameter found by regression analysis.In a study conducted by Whittles et al., the values of −0.8616 and 2.613 × 10 −13 m 2 were determined for m and k 1 , respectively.[22].
The most representative permeability model proposed by Palmer and Mansoori [23] and Palmer et al. [24] was expressed as where φ is porosity; φ 0 is initial porosity; K is the bulk modulus in GPa; M is the constrained axial modulus in GPa; γ is the solid compressibility, which is also known as grain compressibility; p s is atmospheric pressure under standard condition, which is expressed in MPa; and f is a fraction whose value varies between 0 and 1.A value of 0 and 1 is recommended for vertical and horizontal cleats, respectively.In principle, the porosity-permeability or stresspermeability model is defined as the change in effective stress, which changes the internal structure of coal and rock fracture and deforms the matrix or even leads to its failure.The permeability of the rock mass is affected consequently.This is also the main theoretical principle for evaluating permeability enhancement and pressure relief of protective seam mining.Cheng et al. [25] proposed a theoretical model for permeability evolution and pressure relief of coal, which contains gas in deep PCSM.The existing ground stress in the 2 Geofluids deep seam, which changes effective stress, was considered, and it controlled permeability directly or indirectly.Furthermore, the existing fractures as well as new fractures were further extended and governed by the PCSM process.The influence of coal matrix deformation on permeability alteration can be ignored.By assuming effective stress, gas adsorption/desorption, and strain as variables, a theoretical model for enhancing coal permeability was established [25]: Here, ε is the strain, the subscript "0" refers to the initial point, the subscript "m" refers to the coal matrix, the subscript "V" refers to the volume, and f m is the influence factor, which ranges between 0 and 1, and the desorption and shrinkage of the coal matrix affected the strain of fracture.
Because the relationship between permeability and stress is different for fractured and intact coal and rock masses, the failure of the rock mass between protective and protected seams was considered in any study, which was associated with GSM during PCSM.

Gas Seepage Model of Upper Protective Seam Mining
The gas seepage model (GSM) was established with upper PCSM, which mainly depends on the gas source and seepage media.The gas source is mostly dependent on the protected seam, and seepage media is largely associated with interbedded rock mass.As shown in Figure 1, King [26] determined various types of coal seam gas migration by using numerical simulation.In the protective seam mining process, a large amount of gas was desorbed from the bottom of the protective seam, providing pressure relief as well as coal and rock expansion.This was the main source of gas, which flows into PSWF through the interbedded rock mass.Different failure mechanisms are associated with the interbedded rock mass, which depends on its lithology and its distance to the floor of the protective seam.The corresponding gas permeability also varies in this experiment.There are two types of interbedded rock masses, namely, elastic intact and plastic fractured rock masses.Flow process could be divided into three stages from the protected seam to PSWF: (1) gas desorption occurs by providing pressure relief in the protected seam, (2) the gas seepage stage occurs through the intact rock mass, and 3) the gas seepage stage occurs through the fractured rock mass.To simplify the developed model, the following assumptions were made for analysis: (1) The entire model system has the same constant temperature (2) The interaction between water and other gases was ignored, i.e., only gas is modeled (3) Gas flows along the shortest path in the fractured rock (4) Gas is regarded as an ideal gas; the dynamic viscosity and the diffusion coefficient are considered as constant values Several studies show that there are three main states for in situ coal seam gas, namely, absorbed gas, free gas, and dissolved gas.In the process of pressure relief and expansion of the coal seam, the released gas was mainly absorbed with the formation of the coal seam.The desorption process works as the reverse mechanism of the absorption process, which is governed by the Langmuir equation [27].In gas desorption stage, gas diffuses into elastic rock as pores of the intact rock mass are relatively small.The process of diffusion obeys an unsteady and nonequilibrium desorption model, which is known as Fick's second law.
Herein, Q is the quantity of gas diffusing in unit time, which is expressed in m 3 /s; D is the diffusion coefficient; c − c 0 is the concentration difference; H is the distance of the cleat, which is expressed in meters; and A is the pore area, which is expressed in m 2 .The normal gas is assumed to obey the state equation of ideal gas: Herein, P is the ideal gas pressure, which is expressed in MPa; V is the ideal gas volume, which is expressed in m 3 ; N is the amount of substance of gas, which is expressed in mol; T is thermodynamic temperature of ideal gas, which is expressed in Kelvin; R is the ideal gas constant.
Equation 9 illustrates the relationship between the concentration of gas and the amount of gas substrate: Gas pressure of the protective seam was denoted as P 1 , and the distance was denoted as H. Gas pressure of intact rock mass was assumed to be equal to atmospheric pressure (P 0 ).Equation (10) was used for calculating gas quantity (Q 1 ), which diffuses into intact rock mass: When gas diffuses into the intact rock mass, it enters the seepage stage of intact and fractured rock masses.3 Geofluids Equation (11) shows that gas seepage in the rock mass follows Darcy's Law.
Here, Q is the quantity of seepage gas, which is expressed in m 3 /s; K is permeability of rock, which is expressed in m 2− ; P − P 0 is the pressure difference, which is expressed in Pa; l is the distance of seepage, which is expressed in meters; A is seepage area, which is expressed in m 2 ; and μ is the dynamic viscosity coefficient, which is expressed in MPa•s.Equation (12) was used to calculate gas pressure after the stage of intact rock mass.As shown in (12), the vertical length of the intact rock mass is known as l 1 .According to flow conservation principle, the permeability of the intact rock mass is known as k 2 .
Gas flows (Q 2 ) from the protected seam to PSWF in fractured rock mass, which was calculated by utilizing (13).Rock permeability is called k 3 , l 2 is called vertical length, and ventilation pressure of PSWF is p 2 .
After PCSM, ( 14) was used for calculating maximum and minimum principle stresses at any point for the fractured rock mass [28].
Herein, q 1 is the vertical stress that acts on the floor, which is expressed in MPa; q 2 is the horizontal stress, which is expressed in MPa; and α is the angle of the coal seam, which is expressed in degrees.As shown in Figure 2, θ = θ 1 − θ 2 , where θ 1 and θ 2 are angles of the connection between the point and the floor, and they are also expressed in degrees; l is the vertical distance between the point M and the floor of the protective seam, which is expressed in meters; and γ is the rock density, which is expressed in MN/m 3 .
The permeability-stress model was chosen for determining permeability of intact and fractured rock masses.By using laboratory data fitted with (3) and ( 4), the shortest path was chosen, i.e., θ = 0. Suppose that stress of the intact rock mass is γl, the original rock mass stress is γh, and h is the buried depth of the protected seam.This assumption is based on the fact that in the actual production process, most gases released from the protected seam do not diffuse into PSWF.By continuously advancing the working face, gas pressure of the protected seam is reduced gradually.Equation ( 15) must be multiplied with a loss of coefficient δ.By substituting parameters into (13), ( 15) is obtained as follows:  The coal mining process transfers the overburden load to the coal body ahead of the working face and compacts falling rock in the goaf by forming a peak region for bearing pressure.In front of the working face, the floor rock was in a state of compression, experiencing compressive displacement.In the goaf behind the working face, the floor rock was in a stress relief zone, and it experienced expansive displacement.Shear failure normally occurs in floor rocks because of shear deformation, which are located at the border of compression and expansion regions.In a conventional advancing process of the working face, there are three stages of stress in the floor: stress escalation before mining, stress drop after mining, and stress recovery.These stages of stress repeatedly occur with advancing of the working face.In surrounding floor areas of the working face, stress distribution on the floor rock mass was affected by mining.It was divided into four horizontal regions: (1) original rock stress region, (2) compression region, (3) expansion region, and (4) stress recovery region.The distribution of floor strata was divided into the following four vertical zones: (i) broken zone, (ii) newly damaged zone, (iii) originally damaged zone, and (iv) fractured zone [29].Figure 3 shows the areas of the floor, which were influenced by bearing pressure.
According to elastic-plastic mechanics and Mohr-Coulomb strength theory, a plastic zone was formed in the floor rock with plastic deformation.It was within a certain range of the working face floor when bearing pressure reached or exceeded the pressure limit value.
When the bearing pressure reached maximum capacity, part of the rock failed completely by joining plastic zones of rock mass in the surrounding region of bearing pressure.With plastic deformation of the rock mass, a continuous slip surface is formed in the goaf direction.
The maximum failure depth of the floor, which is caused by stress concentration under ultimate bearing pressure conditions of the floor rock, can be obtained by using (16).The ultimate bearing capacity of strata was calculated using Mohr-Coulomb failure criterion [30].
The width of the plastic area was calculated by using (16).Equation ( 17) was used for determining the maximum failure depth of the floor.
Here, a, which is expressed in meters, is twice of roofcontrol distance in the working face; n m is the maximum coefficient of stress concentration; σ M is the ultimate stress of the coal body edge, which is expressed in MPa; H is the vertical distance from the working face to the ground surface, which is expressed in meters; and β is the internal friction angle of floor strata, which is expressed in degrees.5 Geofluids Equation ( 18) is used for theoretical calculation of failure depth of the floor [31].
Here, h is the depth of the broken zone and L is the length of the working face, both of which are expressed in meters, and R c is the uniaxial compressive strength of the rock, which is expressed in MPa.
In previous studies, statistical analysis was conducted to determine failure depth of the floor during the coal mining process in various mines.Corresponding linear regression equation was formed by using multivariate linear and nonlinear statistical analysis methods.Some of the typical equations are as follows [31][32][33]: Here, M h is the mining height of the coal seam, which is expressed in meters; D is the failure resistance capability of the floor; and I is the factor related to existence of cutting fault or failure situation of the broken zone.By utilizing the aforementioned equations for estimating floor failure depth, the length of intact and fractured rock masses in GSM was obtained.
For a case study in a certain coal mine in Hancheng city, Shanxi province, China, the verification procedure and associated calculations are described in the following subsections.
The actual two PSWFs were used to verify the model.The distance along the strike was 1142 m, with 180 m of inclined width for PSWF4216.For PSWF4215, the distance along the strike is 1158 m with 200 m of inclined width.These two working faces are located in the No.2 coal seam; the thickness of the coal seam was about 0-2.5 m and the average thickness was 1 m.The depth of the No.2 coal seam was 430 m.The spacing of the No.2 coal seam and the underlying No.3 coal seam (the protected seam) was in the range of 9.8-20.7 m, with the average spacing being 15.5 m.The gas content of this protective coal seam is not high, and the pressure was 0.55 MPa.The lower No.3 coal seam had outburst risk.When the relative gas emission of the No.3 coal seam rose to 27.82 m 3 /t, gas pressure became 2.03 MPa and coal seam was of low permeability.Table 1 shows the surrounding rock of the working face and mechanical parameters of the rock mass.Figure 4 shows the borehole diagram of the working face.
In order to verify the accuracy of the developed GSM in the protected seam, floor failure depth during excavation process of the working face was determined first.Using theoretical equations ( 16)- (18) and the data obtained from mining (H is 430 m, β is 40 °, L is 180 m, γ is 0.025 kg/cm 3 , and R c is 21.59 MPa), the failure floor depths were estimated to be 15.2 m and 17.52 m.Using regression equations 19-21 and the obtained values (M h is 1.2 m, α is 0 °, D is 0.2, I is 0, and f is 3), floor failure depths were estimated to be 13.73 m, 15.80 m, and 16.44 m.By considering the results of both theoretical and regression equations, the calculated floor failure depth was found to be in the range of 13.73 to 17.52 m, with an average failure depth of 15.625 m.This value was found to be very close to the actual average interlayer spacing in the mine.To further determine failure depth of the floor,   7 Geofluids is determined by normal (K n ) and shear (K s ) stiffness of the contact.The contact properties, cohesion, tensile strength, and friction angle define the contact strength.Fractures are created when the stress level at the interface exceeds a threshold value in either the tension or shear.The constitutive contact model is shown in Figure 6.
The PSWF starts excavation 50 m away from the left side edge of the model, excavating 5 m at each step for a 100 m distance.Figure 7 illustrates that the working face of advancing floor fractures develops gradually in depth while the range of fracture development enlarges continuously.When PSWF achieves an advancing distance of 30 m, the rock mass between protective and protected seams is connected by fractures but the development of fracture is sparse.The development gets more intensive when the advancing distance is 35 m.At an advancing distance of 40 m, the protected seam is covered by fractures.This causes the gas to desorb from the protected seam and flow into PSWF.
Based on the results of the numerical simulation and theoretical calculations, it was inferred that floor fracture develops on the protected seam during the PCSM process.Consequently, it can be considered that there is no intact rock mass between protective and protected seams.This assumption was deliberated in the following calculations.

Gas Seepage Model Verification
It is challenging to physically measure most of the parameters in GSM, while some parameters remain unchanged under same geological conditions.However, it can be concluded 8 Geofluids from the aforementioned analysis that interbedded rock masses are all fractured rock masses.Equation ( 15) can be simplified into Based on the actual data for the considered PSWF4216 in Hancheng, China, the parameters associated with gas are as follows: P 1 = 2 03 MPa (the gas pressure of the protected coal seam), P 2 = 0 1 MPa (the gas pressure of PSWF), γ = 0 025 MN/m 3 , and the floor area (inclined width of PSWF × floor width) of the working face is 180 m 2 and 200 m 2 for PSWF4216 and PSWF4215, respectively.The value of μ for gas is 1.75 × 10 −11 MPa•s.The m and k 1 are assigned as −0.8616 and 2.613 × 10 −13 m 2 , respectively, which were introduced in the literature [22].By substituting the above data into (22), the gas emission from the protected face was calculated by using Q g m 3 /min = 60Q 2 = 53986 01 × δl −0 8616  23 Between protective and protected seams, the altering values associated with interlayer spacing were obtained from a gas drainage borehole.This was placed in the floor of the air return gateway.After the working face advances to a distance of 30 m, floor fractures of PSWF can develop into the protected seam.The data from ventilation air methane in the working face were also used after it advances to a distance of 30 m.As shown in Figure 8, the scatter diagram and fitting curve of ventilation air methane were obtained with interlayer spacing of PSWF4216 along the corresponding interlayer spacing data.
The loss coefficient (δ) was calculated by using the scatter diagram of ventilation air methane (Figure 8).The correlation coefficient (R 2 ) of the fitting equation was 0.9022, which shows an acceptable fitting situation for engineering measured data.These findings confirm that (22) was used for calculating gas emission quantity in the protective seam, which corresponds precisely with field measurement.In addition, a loss coefficient (δ) of 1.012 × 10 −3 was obtained.In order to further verify the correctness of GSM, the loss coefficient (δ) of 1.012 × 10 −3 and parameters of PSWF4215 (only the floor area of PSWF4215 was different from PSWF4216.)were substituted into (22).Figure 9 shows the theoretical curve, which was drawn according to (22) of PSWF4215.Figure 9 shows that the measured data can match well with the theoretical curve, which further validates the feasibility and effectiveness of the proposed GSM.

Discussions
In the process of deep seam mining, PCSM has become one of the most significant measures of regional outburst mitigation as the inherent outburst risk is high.During the PCSM process, the phenomenon of methane concentration exceeding the limit usually occurs in the working face.It is very important to understand the factors affecting gas emission of PSWF and to develop an equation for calculating its quantity.In this study, a GSM was developed for the process of the upper PCSM.Using the equation of GSM, the quantity of gas desorbed from the protected seam into the protective seam was estimated.Fick's second law of diffusion and Darcy's flow law were used for this estimation.Moreover, the relationship between permeability and stress was investigated in an elastic-plastic state.The mechanics of the surrounding rock was also investigated in this paper.There are four assumptions and several simplifications used in the establishment of GSM; it greatly simplifies computational complexity and reduces the selection of parameters in engineering calculations.As shown in (15), it is simpler and easier to use GSM.This model does not need to spend a lot of time and money for additional field tests, and there are no additional uncertainty factors in GSM.For engineering application, this model has distinctive advantages and high efficiency.However, temperature, gas component, and moisture content are factors governing assumptions and simplifications.They also affect the permeability of coal and rock, reducing the Gas seepage model (eq.( 15)) Gas quantity diffusing into the intact rock (eq.( 10)) Gas pressure after the stage of intact rock mass seepage (eq.( 12)) Gas quantity f rom the protected seam to the working face (eq.( 13)) Failure depth of the floor Theoretical model (eqs.( 16)-( 18)) Typical model (eqs.( 19)-( 21)) Numerical simulation (UDEC) Permeability of the intact rock (eq.( 3)) Permeability of the fracture rock (eq.( 4) and eq.( 14)) Gas desorption from the protected coal seam (eq.( 7

Conclusions
In the process of PCSM, surrounding rock permeability was divided into two types, namely, failure plastic and elastic rock masses.The relationships between rock permeability and stress were also summarized under these two evolution types, which provides theoretical basis for GSM between protective and protected seams.Eventually, gas emission from the protected seam was classified in three stages, namely, (1) the gas desorption stage from the protected seam with pressure relief, (2) the gas seepage stage in the intact rock mass, and (3) the gas seepage stage in fractured rock mass (Figure 10).From these equations, it can be concluded that the initial gas pressure of the protected seam, the characteristics of interlayer rocks, and the ventilation pressure of the working face were the main influencing factors of gas emission quantity, which flowed from the protected coal seam to PSWF.
The actual geological conditions of a protective seam were investigated in a coal mine, which was located in Hancheng, China.The GSM was used to calculate gas emission quantity from the protected seam.The modeling results correspond precisely with field measurement by considering the fitting curve of ventilation air methane quantity for PSWF4216 and PSWF4215 with an interlayer spacing and loss coefficient (δ) of 1.012 × 10 −3 .It is worth noting that actual conditions are simplified by the model proposed in this paper.These conditions were based on some aforementioned assumptions.In future study, gas emission from the goaf and its influencing factors were identified.Furthermore, some of the parameters in these equations were adjusted based on the analysis of additional field data measurements from the working face of the PCSM process.Future study will help us to further mitigate gas accumulation risk in the PCSM process.

Figure 2 4 .
Figure2summarizes the entire GSM of the upper PCSM.Equation(15) shows that there are three main factors associated with gas quantity, which can cause gas to flow into PSWF from the protected seam, namely (1) the initial gas pressure of the protected seam (2) interlayer spacing of interbedded rock mass, the percentage of fractured and intact rock mass, and the permeability coefficient of fractured and intact rock mass (3) the ventilation pressure of PSWF and seepage area of the floor 4. Analysis of Floor Failure Depth 4.1.Theoretical Calculation Model for Floor Fracture Development Depth.The most significant part of the analysis is to determine floor failure depth according to GSM of the protected seam, which was introduced in Section 3. The methods of calculating floor failure depth consisted of theoretical calculation and empirical regression equation.The coal mining process transfers the overburden load to the coal body ahead of the working face and compacts falling rock in the goaf by forming a peak region for bearing pressure.In front of the working face, the floor rock was in a state of compression, experiencing compressive displacement.In the goaf behind the working face, the floor rock was in a stress relief zone, and it experienced expansive displacement.Shear failure normally occurs in floor rocks because of shear deformation, which are located at the border of compression and expansion regions.In a conventional advancing process of the working face, there are three stages of stress in the floor: stress escalation before mining, stress drop after mining, and stress recovery.These stages of stress repeatedly occur with advancing of the working face.In surrounding floor areas of the working face, stress distribution on the floor rock mass was affected by mining.It was divided into four horizontal regions: (1) original rock stress region, (2) compression region, (3) expansion region, and (4) stress recovery region.The distribution of floor strata was divided into the

Figure 3 :
Figure 3: Floor failure depth caused by the bearing pressure.

Figure 5 :
Figure 5: Schematic diagram of the numerical model.

Figure 10 :
Figure 10: Calculation process of the gas seepage model.

Table 1 :
Physical and mechanical parameters of the rock and coal.
4.2.Numerical Simulation for Floor Fracture Development Depth.The Mohr-Coulomb constitutive model was used for simulation, and rock mechanics parameters are presented in Table 1.A schematic diagram of the model is presented in