Experimental and Numerical Simulation Study of Hydraulic Fracture Propagation during Coalbed Methane Development

The extraction of low-permeability coalbed methane (CBM) has the dual signi ﬁ cance of energy utilization and safe mining. Understanding hydraulic fracturing mechanism is vital to successful development of CBM. Therefore, it is important to improve the law of hydraulic fracture propagation in coal and rigorously study the in ﬂ uencing factors. In this paper, laboratory experiments and numerical simulation methods were used to investigate the hydraulic fracture propagation law of coal in coalbed methane reservoir with natural fractures. The results show that the maximum and minimum horizontal in situ stress and the di ﬀ erence in stress signi ﬁ cantly a ﬀ ect the direction of crack propagation. The elastic modulus of coal, the mechanical properties of natural fractures, and the injection rate can a ﬀ ect the fracture length, fracture width, and the amount of fracturing ﬂ uid injected. To ensure the e ﬀ ectiveness of hydraulic fracturing, a reservoir environment with a certain horizontal stress di ﬀ erence under speci ﬁ c reservoir conditions can ensure the increase of fractured reservoir and the controllability of fracture expansion direction. In order to increase the volume of fractured reservoir and fracture length, the pumping speed of fracturing ﬂ uid should not be too high. The existence of stress shadow e ﬀ ect causes the hydraulic fracture to propagate along the main fracture track, where the branch fracture cannot extend too far. Complex fractures are the main hydraulic fracture typology in coalbed methane reservoir with natural fractures. The results can provide a benchmark for optimal design of hydraulic fracturing in coalbed methane reservoirs.


Introduction
Coalbed methane (CBM) is a valuable resource, and meanwhile, it can affect the safe production of coal. China's CBM reserves are very large, and the amount of such resources at the depth of 0-2000 meters has reached 3:68 × 10 12 m 3 . However, one of the key issues restricting the development of CBM in China is that the permeability is generally very low. The coal reservoir fracture system is the main channel for gas flow that controls the permeability characteristics of the coal reservoir [1]. Through hydraulic fracturing to increase permeability, hydraulic fracture can conduct the original natural fractures in coal to a certain extent, which can greatly increase the production of coalbed methane. The mutual coupling relationship between hydraulic frac-tures and natural fractures has become a research hotspot in recent years [2,3].
Scholars have conducted extensive research on the relationship between natural fractures and hydraulic fractures in shale and tight sandstone reservoirs [4][5][6][7]. The fine description of natural fractures in coalbed methane reservoirs and the quantitative analysis of fractures are particularly important for describing the growth of hydraulic fractures in coal. At present, microseismic and geoelectric methods are widely used in industrial field to monitor the fracture propagation. These methods enable monitoring the range of stress and fluid propagation. The monitoring results are much longer than the real effective support fracture section and cannot distinguish the specific fracture shape. In addition, the experimental method is an effective way to study the problem, but the size of the laboratory test blocks is very small. In addition, the hydraulic fracture in the experiment starts and expands inside the closed test block, where we cannot directly observe the fracture initiation and propagation process. The monitoring methods of the hydraulic fracturing experiment in a laboratory include computer tomography and acoustic emission monitoring, However, these technologies do not allow visual observation of the dynamic propagation process of the fractures [8][9][10].
With the development of numerical simulation technology, many scholars have adopted the method to study the change in morphology and expansion of hydraulic fractures in different reservoirs. The numerical simulation methods include finite element method (FEM) [11][12][13], twodimensional particle flow code (PFC) [14][15][16][17], continuousdiscontinuous coupling method [18][19][20], and discrete element method (DEM) [21][22][23], which can visualize and quantify the fracture morphology and fracturing process. Among these methods, the finite element method is the most widely used one. The cohesive unit method model based on the ABAQUS platform is the most commonly used model for analyzing hydraulic fracturing of reservoirs.
However, the cohesive element method based on the ABAQUS platform cannot simulate the random propagation of hydraulic fractures in fractured reservoirs. Therefore, a new method of random propagation of hydraulic fracture based on a zero thickness cohesive element embedded in finite element mesh is established by using a grid node splitting method. This has been carried out based on the characteristics of natural fracture growth in coalbed methane reservoir and the topological data structure of element nodes.
This method identifies the random expansion process of hydraulic fractures through secondary development, which can make up for the insufficiency of the built-in cohesion unit of the ABAQUS platform to effectively simulate the random propagation of hydraulic fractures. Numerical examples are used to study the influence of horizontal stress difference and natural fractures in the reservoir on the propagation process of hydraulic fractures, which can accurately describe the random propagation behavior of complex hydraulic fractures and provide references for numerical simulation of fractured reservoirs.
In this paper, the effects of stress and fluid injection rate on the hydraulic fracture propagation in the similar material block with fractures were experimentally studied. Also, a global cohesive zone model was established to simulate the interaction between hydraulic fractures and natural fractures, and the effects of four factors on hydraulic fracture propagation were discussed. Furthermore, a comprehensive discussion on the properties and causes and control factors of hydraulic fractures in the coal containing natural fractures was carried out.

Specimen Preparation.
Coal strata often contain fractures of different scales that are usually produced by geological tectonic movement. Under the same tectonic background, the development law of surface rock joints can predict the direction and density of coal reservoir fractures. Coal mining often causes fracture development, where these fractures extend from the surface to the deep part of the stratum (Figure 1(b)). As shown in Figure 1(a), a statistical analysis of the fracture development of gas reservoirs in the Qinshui Basin is concentrated in two directions; the dominant joint development direction is in 40°~80°and 320°~350° [24].
In order to simulate the fracture growth morphology of hydraulic fractures in the coal reservoir of interest, containing natural fractures. The material of cement : gypsum : coal dust : water = 3 : 1 : 1 : 3 was mixed together by experimental measurements, and the mechanical parameters (Table 1) were similar to those of coal seam #15 in Qinshui Basin after it was solidified and maintained for one month. As shown in Figure 2(a), the mentioned proportion of the materials was mixed evenly. First, pour it to half the height of the 100 mm × 100 mm × 100 mm cubic mold. A thin cardboard of 25 mm in length and width (as shown in Figure 2(b)) was inserted into the middle of the test block on the bottom surface. After the material in the mold lost fluidity, fill the mold slowly and vibrate evenly to eliminate the weak contact surface of twice pouring slurry. After the test  2 Geofluids blocks were demolded and maintained at constant temperature and humidity for 28 days, a hole with a diameter of 8 mm and a depth of 55 mm was drilled at the top center of each cube, and the wellbore were glued with high-strength epoxy resin AB adhesive as shown in Figure 2(c) and left to stand for 24 hours before conducting the experiment.

Experimental Setup and Method
2.2.1. Experimental Setup. Figure 3 shows a true triaxial hydraulic fracturing experimental system. The experimental system consists of a high-pressure water pump system, a true three-axis servo loading system, and a monitoring and control system. The triaxial servo loading system applies the true triaxial stress by hydraulic jacks in three directions on the cubic test block to simulate the suit stress. A high-pressure pumping system injects high-pressure water into the test block. The monitoring and control system is used to accurately control the pressure and flow rate of high-pressure water.   3 Geofluids multistep loading mode was used to avoid sample shear damage caused by unbalanced loading of three-dimensional stresses. First, the three-dimensional stress was applied at the value of the minimum horizontal principal stress at the same time and remains stable. Then, the vertical stress and the maximum horizontal stress slowly increased to the value of the maximum horizontal principal stress. Finally, the vertical stress was slowly increased to the design value, and after the triaxial stress loading was completed and kept stable for 30 minutes, the high-pressure hydraulic fracturing operation can be carried out. The fracturing pump was controlled by a constant displacement mode, while the computer collects real-time information on the high-pressure pump pressure, piston displacement, and three-way stress.

Experimental Results and Analysis
2.3.1. Injection Pressure Characteristics. Figure 4 shows the detailed pump pressure curves for samples #2 and #4. Along with the gradual increase of injection pressure, both samples can reach the maximum pump pressure at the initial fracture of the specimens, which corresponded to the opening of the hydraulic fracture, and the subsequent rapid decrease in injection pressure means that the fracturing fluid entered into the hydraulically induced fractures. Subsequently, the injection pressure fluctuated frequently, which means that the hydraulic fractures propagated during the fracturing process, forming more secondary fractures, and the fracturing fluid filled these new fractures. As fluid seeped into the new fractures, the pressure in the old fractures decreased rapidly, which also resulted in a rapid decrease in pump pressure. In this experiment, a constant injection rate was applied, and the fluid buildup in the new fracture space rapidly increased the injection pressure until the opening of the next fractures. This process was repeated until the experiment was completed, and the increase and decrease of injection pressure fluctuation indicate the continuous creation of fractures [26,27]. Figure 5, the experimental results show that fracture propagation is more complicated for the case of small horizontal deviatoric stress and accompanied by the generation of branch cracks. The major fracture growth direction has a large intersection angle  In the hydraulic fracturing process, a small injection flow rate (10 ml/min) formed a more obvious major fracture. Also, the fracture direction was controlled by the stress. However, when the injection rate was increased to 30 ml/min, the fracture plane did not extend along the maximum principal stress direction. It first extended perpendicular to the wellbore axis with a portion deflected and extended toward the maximum stress direction and then through the sample ( Figure 6). This shows that the high rate of injection can cause the sample to not expand along the direction of maximum horizontal stress at the fracture initiation. This is consistent with the findings in the literature [28,29].

Numerical Simulation
The reservoir hydraulic fractures extend far away, and the fracture extension morphology is often derived indirectly by means of microseismic monitoring or geodetic electrical methods. The small size of experimental samples for indoor experiments also makes it difficult to accurately simulate the real hydraulic fracture extension morphology. Numerical simulation can clearly show the hydraulic fracture propagation under the influence of various factors such as ground stress and natural cracks, which are an effective means to invest in the hydraulic fracture propagation. In this paper, ABAQUS software was used to establish a global cohesive zone model which can simulate the hydraulic fracture extension process in reservoirs containing conjugate fractures. Taking Qinshui coalfield as the background, a twodimensional numerical model was established.
3.1. Coupling Equations. In this simulation, the tensileseparation criterion with degraded cell stiffness in ABAQUS was used to simulate the fracture initiation and propagation. As the stress on the cohesive unit became larger, the unit stiffness degraded gradually; and when the unit stiffness dropped to 0, the cohesive unit started to fracture and expand. This simulation used the maximum positive stress criterion that is when the unit in either direction reaches its critical stress, the unit begins to fracture.
where σ max n refers to the most tensile stress that the unit can withstand in the vertical direction and τ max s and τ max t refer to the maximum shear stress that the unit can withstand in both directions.
A dimensionless damage factor (D) is introduced, taking a value range between 0 and 1. When D = 0, no damage occurs to the material; when D = 1, the material is completely damaged, and fractures form and continue to propagate; and when 0 < D < 1, the material is undergoing damage; the following equations are applicable: where σ n ′ ≥ 0 represents a cohesive unit subjected to tensile stress and σ n ′ < 0 represents the cohesive unit subjected to compressive stress. σ n , σ s , and σ t represent the normal stress component and the two tangential stress components of the cell under linear elastic condition before the material is damaged.
The fluid flow direction in the cohesive unit is divided into a tangential flow along the cohesive unit and a normal flow perpendicular to the cohesive unit. The tangential flow follows the Newtonian flow formula: where q denotes the fracturing fluid discharge volume, m 3 /min; d denotes the fracture opening width, m; μ denotes the dynamic viscosity of fracturing fluid, mPa·s; and p is the flow pressure, MPa. The fluid in the fracture is mainly tangential flow, and a small amount of fluid will pass through the normal flow of the fracture to the formation on the upper and lower surface of the cohesive unit. The normal seepage formula of fracturing fluid in the cohesive unit is as follows:   Figure 7 shows the 2D reservoir model, which is a square with a side length of 50 meters and two sets of natural fractures. The minimum horizontal stress direction was assumed to be 30°and 150°with respect to the two natural fracture sets in this model. The center of the simulation model is the injection point, and the injection hole is oriented in the same direction as that of maximum horizontal stress. X is the direction of minimum horizontal principal stress, and the Y-axis is the direction of maximum horizontal principal stress. The boundary condition with full displacement constraints was used around the model to limit its horizontal displacement or deformation. Table 3 lists the parameters used in the model.

Simulation Results
3.3.1. Changes in Stress Field during Fracture Propagation. As the hydraulic fracture opens and expands, the stress field around the hydraulic fracture changes. The hydraulic pressure in the hydraulic fracture will mainly produce induced compressive stress on both sides of the crack, called stress shadow effect. The process of propagation of the dominant hydraulic fracture also forms several small branches of hydraulic fracture, but only one branch can become the main hydraulic fracture. The other fractures close immediately, which is caused by the stress shadow effect. We also observed the change of stress during hydraulic fracture generation and expansion by simulating hydraulic fracturing with horizontal stress difference HSD = 9 MPa.
(1) Variation of the σ H and Shear Stress τ. From Figure 8(a), the effect of Hydraulic fracture expansion on σ H was reflected in a significant increase in a small area of the fracture tip region and a smaller increase in other regions. σ H decreased slightly on both sides of hydraulic fracture. From Figure 8(b), it can be seen that the variation of the shear stress τ is almost 0 at the fracture tip along the fracture length and in the regions on both sides of the fracture. There were obvious changes in shear stress on both tips of the fracture, and the shear stress on both sides of the fracture end point increased on one side and decreased on the other side, which made the fracture end more susceptible to shear damage.
(2) Variation of the Minimum Horizontal Principal Stress σ h . From Figure 8(c), it can be seen that σ h increased to various degrees in the whole model area after the creation of hydraulic fracture, especially on both sides perpendicular to the fracture extension direction. The increase in σ h increased the breaking pressure required to form other hydraulic fractures around the initial fracture and narrowed the width of the generated fractures, which was not conducive to the development of other branch fractures in the same direction. Figure 9 shows the morphology of hydraulic fracture when five stress differences HSD = 1 MPa, 3 MPa, 5 MPa, 7 MPa, and 9Mpa were used in the numerical simulation in this paper. Figures 9(a) and 9(b) show that HDS is small, the expansion of hydraulic fractures was obviously affected by natural fractures, and the expansion path was relatively curved. The length of hydraulic fractures on one wing was obviously greater than that of the other wing. The length of the fracture has obvious characteristics of asymmetric expansion. Figures 9(c)-9(e) show the expansion path of the hydraulic fracture under the condition of increasing HSD which was straighter than that of Figures 9(a) and 9(b), and the expansion of the hydraulic fracture was more balanced in both wings.

Main Fracture Morphology and Controlling Factors.
The stress shadow effect of hydraulic fractures directly affects fracture width, morphology, extension direction, and fracture initiation pressure, which in turn affects the production, recovery, and economic efficiency of horizontal well fracturing design. Wang [12] pointed out that the stress shadow effect and resistance-related fluid distribution were   As shown in Figure 9(f), the hydraulic fracture morphology of this simulation was in good agreement with the complex fracture morphology pointed out by Wang, and the hydraulic  Comparing the experimental results ( Figure 5 (#2)) and numerical simulation results (Figure 9(b)), when the HSD = 3 MPa, there is a more obvious angle deviation between the fracture and σ H . By comparing Figure 5 (#4) and Figure 9(d), when HSD is 7 MPa, the fracture expansion direction tends to the direction of σ H . The experimental and simulation results verify that the in situ stress has an obvious control effect on the horizontal principal stress. In this simulation, when the local stress difference HSD > 5 MPa, that is, when K > 0:625, the fracture direction and σ H were relatively uniform. The results shown in Figure 9 also demonstrate that the propagation of fractures in fractured coal is obviously controlled by the in situ stress. Even when a large number of natural fractures exist, the main propagation path has not completely changed. Especially when HSD is large, the fractures tend to be single-plane fractures. In the case of high stress difference, the control effect of fracture propagation along the main stress direction is very strong, while in the case of small stress difference, the extension direction of the main fracture is more obviously affected by natural fractures.

Expansion Pattern of Hydraulic Fractures Encountering Natural
Fractures. The numerical simulation shows that there are three forms of expansion when hydraulic fractures encounter natural fractures: (1) stagnation of hydraulic fractures after propagation along natural fractures; (2) the hydraulic fracture extends in the natural fracture for a certain distance, then extends out of the natural fracture and continues to expand; and (3) hydraulic fractures propagate directly through natural fractures. Figure 9(a) shows the fracture extension process when HSD = 1 MPa, the hydraulic fracture in the upper part stagnates at the fracture along the natural fracture surface while the hydraulic fracture in the lower part continues to extend forward; Figure 9(b) shows the fracture extension process for the case of HSD = 3 MPa. The hydraulic fracture in the upper part first extends within the fracture and then continues to extend forward from the right boundary of the natural fracture. As shown in Figure 9(c), the hydraulic fracture directly passes through the natural fracture and expanded forward when the HSD = 5 MPa, where the shear effect of the hydraulic fracture is obvious. The HSD determines the form of hydraulic fracture through the natural fracture when the approach angle is large. Similarly, Figures 9(d) and 9(e) similarly show the shearing effect of hydraulic fractures on natural fractures under high stress conditions. The higher shearing effect caused friction particles at the interface to plug some of the fractures and thereby decrease the fracture permeability. This situation is not conducive to the development of coalbed methane.

Parametric Analysis on Fracture Propagation
(1) Effect of Elasticity Modulus. To investigate the effect of elastic modulus on hydraulic fracturing, the elastic modulus of coal was set to 2.6, 4.6, and 6.6 GPa, while the elastic moduli of natural fracture were taken as 0.80 GPa, and the model with HSD = 5 MPa was selected; the numerical results are shown in Figure 10. As elastic modulus of coal increases, the fracture fluid injection volume and maximum fracture width of the main fracture decreased, and the length and fracture pressure of the main fracture increased. The analysis shows that the maximum fracture width of the main crack is strongly correlated with the variation of the elastic modulus. When the elastic modulus increased from 2.6 GPa to 4.6 GPa, the maximum fracture width of the main fracture decreased by about 37.0%, the fracturing fluid injection volume  Geofluids decreased by about 23.9%, the main fracture length increased by 13.5%, and the fracture pressure increased by 34.5%. This is because the increase of the elastic modulus of the matrix makes the fluid carry more energy into the fracture [30].  Table 3. The calculated results ( Figure 11) show that the breaking pressure tends to increase with the increase of natural fracture tensile strength. The length and maximum width of the primary fracture show a decreasing trend, while the fracture volume also shows a decreasing trend. The total volume of hydraulic fractures decreased with increasing tensile strength of natural fractures, because the volume of fractures within the coal matrix produced per unit energy is smaller than that of hydraulic fractures extending within the natural fractures [31][32][33][34][35].
(3) Effect of HSD on Fracture Propagation. To study the effect of in situ stress difference on fracture extension, the HSD were set to 1, 3, 5, 7, and 9 MPa, while keeping the other parameters constant. The results of simulation are shown in  9 Geofluids Figure 12. With the increase of HSD, the length of the main fracture tends to decrease, the width of the main fracture and the volume of fracturing fluid injection increase, and the change in the fracture pressure was not evident. The length of the fracture tends to decrease as the stress difference increases, which is caused by the flatter extension of the fracture when the in situ stress difference increases. When the HSD was large, especially greater than 5 MPa, the hydraulic fracture became more flat and the complexity of the fracture was reduced. The increase of the HSD has a great effect on the fracture fluid injection volume and maximum fracture width of the fracture, while other factors have little influence. When the HSD increases from 1 MPa to 9 MPa, the injection volume of fracturing fluid increases by 33.02%, the maximum fracture width increases by 34.31%, and the fracture length of the main fracture decreases by 7.27%. This shows that HSD controls the geometry of the fractures, and when HSD is small, the angle between the fractures' morphological expansion and the maximum principal stress will increase significantly. This is because the higher HSD increases the initial shear stress in the coal seam, especially in locations with natural fractures, where hydraulic fractures are more likely to expand [36][37][38].
(4) Effect of Injection Rate on Fracture Propagation. In order to explore the influence of injection rate of fracturing fluid on hydraulic fracture propagation, the injection rate was set to 0.1, 0.2, 0.3, and 0.4 m 3 /min. Figure 13 shows the main geometric parameters of hydraulic fractures in different injection rates when the HSD was 5 MPa. The other parameters are presented in Table 3. According to the data in Table 3, the fracture pressure of the fracture has a tendency to increase with the increase of the fracturing fluid injection rate. The fracture volume and length have an obvious tendency to decrease. Meanwhile, from the fracture morphology, the hydraulic fracture is more developed in both wings when the pumping rate is 0.1 m 3 /min. Also, the phenomenon that one wing is fully expanded and the other is not exists when the pumping rate was 0.2, 0.3, and 0.4 m 3 /min. Therefore, the pumping speed of hydraulic fracturing was not conducive to the formation of complex fracture network and the high production of CBM when the pumping speed was too large, which is consistent with the findings by the references in the literature [39,40].
In order to obtain better fracturing effect in industrial field, appropriate stress difference should be utilized during the fracturing process to ensure a larger volume of hydraulic fracture. The stress difference should not be too low in order to make the hydraulic fracture extend farther, pass through more natural fractures, and obtain more complex fractures. The injection rate should also not to be too high in order to ensure that the volume of the hydraulic fracture is larger and the hydraulic fracture extends farther [41,42].

Conclusions
Understanding the law of hydraulic fracture propagation in low-permeability CBM reservoirs with natural fractures is essential for evaluating the complexity of hydraulic fractures and improving the efficiency of CBM development. This study uses a combination of experiments and numerical simulations to study the expansion pattern of hydraulic fractures in fracture-bearing coal reservoirs. The following conclusions can be drawn: (1) The existence of stress shadow effect makes the hydraulic fractures expand along the main fracture trace, and the branch fractures close quickly without extending too far. Complex fractures are the main fracture morphology in naturally fractured CBM reservoirs. It is difficult for hydraulic fractures to form a complex tree-like seam network in coal reservoirs. The stress shadow effect is caused by the increase of 10 Geofluids the minimum principal stress in the surrounding coal of the hydraulic fracture (2) The horizontal stress difference and natural fractures jointly determine the overall shape of the hydraulic fracture network. When HSD is small, the hydraulic fractures are more complex and the morphology of the hydraulic fractures is significantly influenced by the natural fractures (3) To ensure the efficiency of hydraulic fracturing, a reservoir environment with a certain horizontal stress difference under specific reservoir conditions can ensure the increase of fractured reservoir volume and the controllability of fracture expansion direction. In order to increase the fractured reservoir volume and fracture length, the pumping rate of fracturing fluid should not be too large

Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest
No potential conflicts of interest were reported by the authors.