A numerical investigation of hydraulic fracturing on coal seam permeability based on PFC-COMSOL coupling method

Hydraulic fracturing and permeability enhancement are effective methods to improve low-permeability coal seams. However, few studies focused on methods to increase permeability, and there are no suitable prediction methods for engineering applications. In this work, PFC2D software was used to simulate coal seam hydraulic fracturing. The results were used in a coupled mathematical model of the interaction between coal seam deformation and gas flow. The results show that the displacement and velocity of particles increase in the direction of minimum principal stress, and the cracks propagate in the direction of maximum principal stress. The gas pressure drop rate and permeability increase rate of the fracture model are higher than that of the non-fracture model. Both parameters decrease rapidly with an increase in the drainage time and approach 0. The longer the hydraulic fracturing time, the more complex the fracture network is, and the faster the gas pressure drops. However, the impact of fracturing on the gas drainage effect declines over time. As the fracturing time increases, the difference between the horizontal and vertical permeability increases. However, this difference decreases as the gas drainage time increases. The higher the initial void pressure, the faster the gas pressure drops, and the greater the permeability increase is. However, the influence of the initial void pressure on the permeability declines over time. The research results provide guidance for predicting the anti-reflection effect of hydraulic fracturing in underground coal mines.


Introduction
Gas disasters significantly affect mine safety. Due to increased mining depth and intensity in China, many lowpermeability gas mines have been transformed into highpermeability gas or gas outburst mines, adversely affecting safety (Al-Rubaie and Ben Mahmud 2020; Du et al. 2020;Liu et al. 2011;Wang et al. 2013). Coal bed methane exploitation can prevent coal and gas outbursts and provide clean energy. An increase in the mining depth has intensified the in-situ stress, decreased the permeability of coal seams, and increased the difficulty of coalbed methane exploitation (Huang and Liu 2018;Liu et al. 2019;Wang et al. 2018). Therefore, improving the permeability and preventing gas disasters are crucial in coalbed methane development (He et al. 2012;Huang et al. 2015;Li et al. 2016;Wang et al. 2017b;Zhao et al. 2020;Zhou et al. 2021). Research has shown that hydraulic fracturing increases coal seam permeability, which is conducive to the exploitation of coalbed methane and the prevention of coal and gas outburst (Cheng et al. 2018;Wang et al. 2014;Zhang 2014).
Hydraulic fracturing of coal seams refers to injecting high-pressure liquid into the coal seam through the borehole using a high-pressure pump. The pressure overcomes the in-situ stress, tensile strength, and cohesion of coal seams, forming a fracture network that improves gas flow. Wang et al. (2014) used PFC2D software to analyze the effect of the macro mechanical properties on crack initiation and size during hydraulic fracturing. The results showed that the main fracture extended into the direction of maximum principal stress. Zhou et al. (2017) used a two-dimensional particle flow simulation program and a smooth joint model to simulate hydraulic fracturing. The authors studied the interaction between hydraulic fractures and natural fractures and compared the results with indoor experimental results and analysis model results. Wang et al. (2018) used the finite element method to study the fracture propagation mechanism of hydraulic fracturing in a coal seam with discontinuous natural fractures. The results showed that the hydraulic fracture network was spindle-shaped with a multilevel branch structure at a high stress difference. Secondary fractures accounted for a large proportion, representing an important part of the hydraulic fracture network in the coal seam. Wang et al. (2017a) used the particle flow code (PFC) to simulate the evolution of the overburden fissure field; the research results provided valuable guidance for gas drainage in coal mine fields. Lyu et al. (2020) analyzed the influence of natural fractures in coal seams on hydraulic fracturing performances. It was observed that the propagation of hydraulic fractures depended on the direction of natural fractures and the principal stress. Yuan et al. (2012) established a mathematical model of hydraulic fracturing in a low-permeability coal seam and analyzed the influence of the water injection pressure and other factors on the characteristics of fracture expansion. The results demonstrated that the fracture length increased linearly, and the fracture width increased exponentially with an increase in injection pressure. The fracture width was larger than the fracture length in the later fracturing stage. Zhao et al. (2019) studied the influence of the coal rock type (bright coal, semi-bright coal, semi-dark coal, and dark coal) and perforation location on hydraulic fracture propagation. It was found that the fracture morphology and proppant distribution were better for bright coal than dark coal.
The literature indicates that many scholars have investigated the fracture propagation type, fracturing influence range, and the influence of existing fractures using experiments and simulations (Adachi et al. 2007;Li et al. 2014;Liang et al. 2017). Many studies were conducted on the influence of coal seam hydraulic fracturing on the gas drainage effect. Fan et al. (2019) proposed a fully coupled mathematical model of hydraulic stress damage considering gas-water two-phase flow to simulate hydraulic fracturing to enhance underground gas drainage. Zhang (2014) established a twophase three-dimensional seepage hydraulic fracturing model for dual porous media and conducted simulations using well test data from a basin in Western China. The results showed that hydraulic fracturing promoted the desorption and diffusion of coalbed methane and significantly improved its production. Zhang et al. (2018) conducted a hydraulic fracturing study of the Nantong Mine in the southeast of the Sichuan Basin. Field investigations showed that hydraulic fracturing significantly increased the methane extraction rate; it was more than 10 times higher than that of conventional boreholes. The field research of Huang et al. (2012) showed that hydraulic fracturing improved the downhole gas permeability, and the gas drainage capacity increased by 15 times. The field application results of pulse hydraulic fracturing (Xu et al. 2017) showed that the proportion of micropores decreased by 7.7%, the proportion of mesopores increased by 23.1%, and the proportion of macropores increased by 2.9%, significantly improving the permeability of the coalbed methane reservoir. Li et al. (2015) observed that the gas desorption index k 1 of the driving face fell below the critical value after hydraulic fracturing. The gas drainage volumes of the fracturing boreholes and pilot boreholes were 3.32 times and 3.07 times higher, respectively, than that of the normal boreholes.
Scholars have investigated coalbed methane development and the impact of hydraulic fracturing on coal seams. However, few studies analyzed the anti-reflection effect of hydraulic fracturing, and no suitable method exists to predict the antireflection effect of coal seams. In this paper, PFC discrete element software is used for hydraulic fracturing simulation. The connected fractures are extracted and imported into COMSOL Multiphysics numerical analysis software for gas drainage simulation to evaluate the permeability enhancement caused by hydraulic fracturing of the coal seam. The research results can be used to predict the anti-reflection effect of hydraulic fracturing in underground coal mines.

Coupled model for hydraulic fracturing
The PFC5.0 software was used to establish a discrete element coupled model to simulate hydraulic fracturing. The fluid was stored in the pore grid, as shown in Fig. 1. Fluid exchange occurs in the adjacent pore grid due to the fluid pressure difference. The coupling of the fluid and solid is achieved by altering the contact force to change the channel pore pressure. The pressure is changed by modifying the mechanical characteristics of the study area. The pore pressure exerts a force on the particles inside the pores.
The flow rate of the fluid exchange can be expressed by the Hagen-Poiseuille equation (Hubbert and Willis 1957;Shimizu et al. 2011): where q is the flow rate, a is the opening of the fluid channel, which is related to the normal force of the two particles, k is the permeability coefficient, ▵ p is the pressure difference between the two pore basins, and L is the length of the fluid channel.
When the bond between the two particles is broken, the opening a is (Shimizu et al. 2011): where d is the distance between two particles, R 1 and R 2 are the radii of two particles, and is the dimensionless multiplier.
At time ▵ t , the change in pore fluid pressure due to fluid flow is (Hubbert and Willis 1957;Shimizu et al. 2011): where K f is the compressive modulus of the fluid, V d is the "domain", i.e., the pore volume, ▵ V d is the pore volume change, and ∑ q is the total flow of the fluid.

Parameter calibration
Parameter calibration is performed to ensure that the macromechanical properties measured in the laboratory correspond to the meso-mechanical properties in the simulation experiments. The PFC contact model used for rock and soil mechanics analysis is the linear parallel bonding model (Pb model). The microscopic parameters of the Pb model correspond to the macro-mechanical parameters, and the microscopic parameters are calibrated by a uniaxial compression test and a uniaxial tensile test. As shown in Table 1, a model with a width of 1 m and a height of 2 m was adopted for the simulation. The interval of the particle radius between the maximum and minimum values is uniform. The minimum radius is 0.01 m, the ratio of the maximum radius to the minimum radius is 1.6, and the porosity is 0.06. Uniaxial compression and uniaxial tension numerical experiments are used to calibrate the parameters. The calibration result is shown in Fig. 2, and the meso-mechanical parameters after calibration are listed in Table 2. These data are used to perform the uniaxial compression and uniaxial tension simulations to obtain the macroscopic mechanical properties corresponding to the meso-mechanical parameters.   Table 3 indicates that the measured and simulated macroscopic mechanical properties are very close.

Analysis of fracturing effect
The numerical model is a two-dimensional model, as shown in Fig. 3. The length and width of the model are 5 m, the minimum and maximum radii of the particles are 0.010 m and 0.016 m, respectively, and the interval between the maximum radius and the minimum radius is uniform. The number of elements in the model is 43,457. The initial stress is simulated by adjusting the wall speed so that the horizontal stress of the model reaches 5 MPa and the vertical stress is 8 MPa. The water injection hole is located in the middle of the model and has a radius of 0.1 m. The fluid is injected at a rated pressure of 15 MPa. The mechanical model of the coal seam is a linear Pb model, and the parameter values are the micro-mechanical parameters obtained from the calibration. Due to the long-term triaxial stress of the overlying strata, the coal seam has a certain strength. According to the stress conditions around the fracturing hole and the classical fracturing theory, many scholars have found that the fracture pressure is related to the horizontal effective stress, tensile strength, and pore pressure (Li et al. 2015). Therefore, the fracture pressure of the coal seam can be calculated as follows: where σ t is the tensile strength of the coal seam (MPa). σ 1 and σ 3 are the horizontal maximum and minimum principal stress (MPa), respectively. P pore is the pore pressure.
The pressure near the injection hole is measured in a circle with a radius of 0.12 m. The relationship between the injection pressure and the number of cracks over time is shown in Fig. 4. The peak pressure of 7.42 MPa is called the burst pressure. The rupture pressure of the  (4) is 9 MPa, which is 1.58 MPa higher than the simulated rupture pressure of 7.42 MPa. The difference can be attributed to the excessively large initial pressure of the model and the short pressure accumulation time. Since the liquid is continuously injected, it accumulates at the injection orifice and the previously formed fracture. The pressure at the crack tip rises over time until it reaches the fracture pressure of the coal seam. New cracks appear in the coal seam, causing a pressure drop at the injection orifice. Therefore, the injection pressure curve shows fluctuations as the fracture expands. It is observed in Fig. 4 that the rate of fracture propagation changes from fast to slow over time. Moreover, the crack propagation is discontinuous, and the pressure at the tip of the crack is not sufficient to damage the coal body initially. The crack continues to expand after the pressure has reached the failure pressure.
The fracturing effect at different time steps is shown in Fig. 5. The fractures expand in the horizontal and vertical directions but dominantly in the maximum principal stress direction. The propagation distance of the fracture is 2.38 m at 100,000 steps. As shown in Fig. 6a, the pressure of the fracturing fluid continues to increase at the fracture tip. Under the combined action of the fluid's pressure and ground stress, the horizontal force of the particles at the fracture tip (the direction of the minimum principal stress) increases. As a result, the speed of horizontal crack propagation increases. Eventually, the particles move in the direction of the minimum principal stress, and the cracks continue to expand in this direction. As shown in Fig. 6b, the particles have all moved in the direction of the minimum principal stress, and the closer they are to the injection hole, the greater the displacement is, the more pronounced the crack propagation is.

Control equation of coal seam deformation
According to the total strain equation of coal, the Langmuir volume strain equation ε s = ε L p p+p L , the balance equation σ ij,j + f i = 0 , and the Cauchy equation σ ij,j + f i = 0 , the governing equation for coal deformation is obtained as: where ε ij is the component of the total strain tensor, G is the shear modulus of coal, K is the bulk modulus of coal, E is Young's modulus of coal, is the Poisson's ratio of coal, and p is the gas pressure in the pores. The effective stress component is defined as � ij = ij + p ij , = 1 − K∕K s is the Biot coefficient, K s is the bulk modulus of the coal particles, ij is the Kronecker number, u i is the component of displacement, σ ij is the component of the total stress tensor, f i is the component of the force on the object, ε L is the Langmuir pressure constant, and p L is the pore pressure. At this time, the measured volumetric strain is equal to 0.5ε L (Zhang et al. 2008).

Gas flow control equation
According to the ideal gas law, mass balance equation, and Darcy's law, the gas flow control equation is: where p a is then atmospheric pressure (101.325 kPa), k is the permeability, is the dynamic viscosity of the gas, ϕ is the porosity, ga is the gas density under standard conditions, c is the coal density, g is the gas density, g is the Darcy velocity vector, Q S is the gas source, t is the time, s and m are the gas contents of the free state gas and adsorbed gas, respectively (Mi et al. 2018).

Porosity and permeability models
It is assumed that the adsorption strain of coal is the same as the adsorption strain of the pore space. When the initial pressure is p 0 , the initial porosity is ϕ 0 , and the initial volumetric strain is zero. According to the site conditions, it is assumed that ε 33 is the uniaxial strain direction, and the overburden load direction is unchanged. ε 11 and ε 22 are lateral strains and both are zero (Kumar et al. 2016;Mi et al. 2018;Pan and Connell 2012;Zhang et al. 2008). The permeability expression for the given uniaxial strain and a constant overburden load is: where M is the limiting axial elastic modulus, The effect of grain compression is considered for the rebound pressure of the fully coupled pore model. If the uniaxial strain and the overburden load remain unchanged, the expression of the rebound pressure is (Mi et al. 2018): The partial derivative of the porosity ϕ with respect to time t is substituted into Eq. (6) to obtain the fully coupled control equation of gas flow under the influence of coal seam deformation: Therefore, Eqs. (5), (6), (7), (8), and (9) define the fully coupled model of coal seam deformation and gas flow. The coupling relationship between coal seam deformation and gas flow is shown in Fig. 7. (7)

Simulation scheme and parameters
The above-mentioned control equations constitute a mathematical model describing coal deformation and gas flow. COMSOL Multiphysics software is used to solve the numerical model using the finite element method. Fractures are important flow channels during gas drainage; thus, different fracturing times result in different fracture network conditions and gas drainage effects. Different initial pore pressures also affect the gas drainage effects. Therefore, in this work, we analyze the effects of different fracturing times and different initial pore pressures on the gas drainage effect. The specific simulation schemes are shown in Table 4. First, the fracture propagation diagram for different fracturing time steps is imported into a CAD program, and the connected fractures are extracted, as shown in Fig. 8. The results are imported into COMSOL Multiphysics software. The calculation grid is densified to ensure the integrity of the fracture boundary. The final physical model of gas drainage is shown in Fig. 9. The size of the model is all 5 m × 5 m, the suction pressure is one atmosphere, and the boundary of the model is rolling support. The main parameters used in the numerical simulation are listed in Table 5.

Fracture and non-fracture models
The established mathematical model is numerically solved using the parameters of the physical model. Figure 10a and b show the gas pressure distribution of the non-fracture   The results show that the gas pressure decreased rapidly in the first few days, but the decrease rate decreased slowly with an increase in the drainage time. The gas pressure of the fracture model decreased faster than that of the non-fracture model. The former shows differences in the distribution of the gas pressure in the horizontal and vertical directions. The gas pressure was lower in the vertical direction (main direction) of fracture propagation than in the horizontal direction. For the quantitative analysis, we used a detection point (1.5, 0) in the non-fracture model (because the model is symmetric in the X and Y directions, only one monitoring point is required) and two in the fracture model, point 1 (1.5, 0) and point 2 (0, 1.5), as shown in Fig. 11. The results show that the gas pressure trend of  the two models was the same. The gas pressure dropped quickly at the beginning of the drainage, and the gas pressure decline rate rapidly decayed and approached zero as the drainage time increased. The gas pressure of the fracture model decreased faster than that of the non-fracture model. In the fracture model, the gas pressure decreased faster at  Jia et al. (2020) at different distances from the fracturing hole is shown in Fig. 12. The results show that the gas pressure decreased rapidly in the initial stage of drainage and then gradually approached 0 with an increase in the drainage time. The further the distance from the fracturing hole, the slower the gas pressure drop was. This result shows that the anti-reflection effect of hydraulic fracturing becomes worse with the increasing distance from the fracturing hole, which is consistent with the simulation results in this paper. Figure 13 shows the permeability increase rate (k − k 0 )∕k 0 of the non-fracture model and the fracture model at different times of drainage. The results show that as the extraction time increases, the permeability increase rate is higher for the fracture model than the non-fracture model. However, the increase rate of permeability in the fracture model is different in the horizontal and the vertical directions, and the increase rate of the permeability is slightly higher in the main direction of fracture propagation. A detection line was created between points (0, 0) and (2.5, 0) of the nonfracture model and the fracture model to quantify the rate of change of permeability over time, as shown in Figs. 14 and 15. The results show that the further from the center point, the smaller the increase rate of permeability was. The decrease rate gradually approached 0, changing from fast to slow, and the permeability increase rate was higher for the fracture model than the non-fracture model. As the drainage time increased, the increase rate of permeability rose significantly. After 20 d of drainage, the permeability at the center and edge of the non-fracture (fracture) model increased by 2.1% and 1.1% (2.27% and 1.7%), respectively. Figure 16 shows the increase rate of gas pressure and permeability at different time steps (S2-S5) of fracturing for different gas drainage times. As the drainage time increased, the gas pressure declined, and the permeability increase rate increased. The gas pressure decrease rate and permeability increase rate were higher in the direction of maximum principal stress (main direction of fracture propagation) than the direction of minimum principal stress for all fracturing time steps. As the fracturing time increased, this difference also increased. Two detection points were used in the S2-S5 model: monitoring point 1 (1.5, 0) and detection point 2 (0, 1.5). The gas pressure at different fracturing time steps is shown in Fig. 17. The trends of the four models were the same. The gas pressure dropped rapidly in the early stage of drainage. As the drainage time increased, the gas pressure decline rate decreased sharply and tended to zero. The longer the fracturing time, the more complex the fracture network was, and the faster the gas pressure dropped. Moreover, as the fracturing time increased, the difference between the horizontal and vertical gas pressure decrease rates increased, and the rate was significantly higher in the direction of the maximum principal stress than the direction of the shows that the influence of fracturing on the gas drainage effect decreases over time. The permeability increase rate at different fracturing time steps is shown in Fig. 18. As the drainage time increased, the rate of increase in permeability changed from fast to slow. Since test point 2 was located in the main fracture propagation direction, the increase in permeability of test point 2 was always greater than that of test point 1. The difference because more pronounced as the fracturing time increased. However, the difference decreased as the extraction time increased. The fracture radius and fracture width near the fracturing orifice during the fracturing process measured by Bao et al. (2021) are shown in Fig. 19. The fracture radius and width increased rapidly at the initial fracturing stage, but the rate of increase decreased sharply over time. This result shows that the influence of fracturing on the gas drainage effect decreased over time, which is consistent with the simulation results in this paper.

Different initial pore pressures
The same detection points were used in the S5-S8 models. Figure 20 shows the gas pressure of the models with different initial pore pressures (S5-S8). The results show that the gas pressure decline rate changes from fast to slow, and the gas pressure difference between point 1 and point 2 decreased over time. The higher the initial pore pressure, the faster the gas pressure dropped. It approached the orifice pressure as the extraction time increased. It required 8 d (point 1) and 6 d (point 2) for the gas pressure to drop below 0.74 MPa for model S5, 10 d and 7 d for model S6, 10 d and 8 d for model S7, and 11 d and 8 d for model S8. The initial pore pressure of model S8 was 1.5 times higher than that of model S5. However, after fracturing, the gas pressure dropped below the standard due to the drainage effect, and it only took 3 more days for extraction. Figure 21 shows the permeability increase rate of the models with different initial pore pressures. The trends of the different models were consistent. The permeability increased rapidly at the beginning of the drainage, but the increase rate decayed rapidly and approached zero as the drainage time increased. The greater the initial pore pressure, the greater the increase in permeability was. However, the influence of the initial pore pressure on the increase in permeability gradually decreased. For a quantitative characterization of the decline rate of the coal seam gas pressure, the gas pressure decline rate D p was defined by the ratio of the gas pressure decline (difference between adsorbed gas pressure p 0 and residual gas pressure p t ) and the adsorbed gas pressure at the time of gas extraction t: Zhang et al. (2014) obtained the gas pressure and gas pressure decline rate over time from models with different initial pore pressures, as shown in Fig. 22. The experimental results show that the gas pressure decreased from fast to slow as the drainage time increased. The higher the initial pore pressure, the faster the gas pressure decreased, which is consistent with the simulation results in this paper.

Application prospects
This study analyzed the influence of hydraulic fracturing on coal seam permeability for different fracturing times and different initial pore pressures. The research results showed  Variation curves of gas pressure and gas pressure drop rate of different initial pore pressure models accordingly. The higher the initial gas pressure, the faster the gas pressure dropped, and the greater the increase in permeability was. The results can be used to guide hydraulic fracturing to enhance permeability in mines with different gas contents. Due to in-situ stress, the main direction of fracture propagation is in the direction of maximum principal stress. Therefore, a difference exists in the permeability change between the horizontal and vertical directions, and this difference increases with the fracturing time. We will focus on eliminating or exploiting this difference in our next study.

Conclusions
In this work, a coupled mathematical model of the interaction between coal seam deformation and gas flow was established, and simulations were conducted to determine the effects of hydraulic fracturing on gas drainage. The research conclusions can be summarized as follows: (1) The fracture growth rate changed from fast to slow as the fracturing time increased and crack propagation was discontinuous. When the pressure at the crack tip was insufficient to damage the coal, no new cracks occurred, and when the pressure was equal to the breaking pressure, the cracks continued to expand. The particles moved in the direction of the minimum principal stress, and the main fractures expanded in the direction of the maximum principal stress. The closer to the injection hole, the greater the displacement was, and the more adequate the fracture expansion was. (2) The gas pressure trends were the same for the fracture model and the non-fracture model. As the drainage time increased, the decline rate of gas pressure changed from fast to slow. The gas pressure decrease rate and permeability increase rate of the fracture model were significantly higher than those of the non-fracture model. However, the change in the gas pressure and permeability was different in the horizontal and vertical directions.
(3) The longer the fracturing time, the more complex the fracture network was, the faster the gas pressure decreased, and the permeability increased, but the greater the difference was between the horizontal and vertical directions. However, this difference gradually decreased with an increase in the extraction time. The fracturing effect on the gas drainage decreased over time.
(4) The greater the initial pore pressure, the faster the gas pressure dropped. It approached the orifice pressure as the extraction time increased. The increasing trend of the permeability was the same for the S5-S8 models. The increase rate of permeability decreased exponentially with the increase in the extraction time. The greater the initial pore pressure, the larger the permeability increase was. However, the influence of the initial pore pressure on the rise in permeability gradually declined over time.