Numerical simulation of fracture propagation in deep coal seam reservoirs

As an unconventional oil and gas resource, coalbed methane has become a key area of unconventional oil and gas exploration worldwide, and hydraulic fracturing is the critical technology for the efficient stimulation of coalbed methane. However, due to the strong heterogeneity of deep coal seam reservoirs and the development of natural cleat fractures, the fracture propagation path is quite complex and difficult to capture. In this study, we propose a hydraulic fracture propagation simulation method using the ABAQUS Software for natural cleat fracture networks and a fluid–solid fully coupled two‐dimensional finite element model with the global embedded cohesive zone model (CZM). CZM is established to simulate the propagation of hydraulic fractures in the heterogeneous combination of coal rock and mudstone. The accuracy of this model is verified by using Blanton's experiment criterion. Finally, the model is used to simulate the influence of cleat angle, cleat density, length ratio of coal rock and mudstone combination, and construction displacement on the propagation behavior of hydraulic fractures. This study aims to deepen the understanding of the mechanism of hydraulic fracture propagation in deep coal seams and provides guidance for deep coal seam fracturing and hydraulic fracturing design.


| INTRODUCTION
Coalbed methane resources are abundant in China and have initially achieved breakthroughs in the exploration and development from shallow to medium-deep layers nowadays.Deep coalbed methane will become a strategic replacement resource for China's natural gas reserves and production, and contribute to the "double carbon" goal. 1,2However, due to the comprehensive influence of coalification, geological structure and depth, and other factors, deep coal reservoirs often have the characteristics of low porosity, low permeability, and strong heterogeneity.In the mining process, they are often accompanied by mixed layers, mainly coal rock and mudstone interaction, which posed a challenge to the exploration and development of deep coal seams.
Hydraulic fracturing can be defined as the process by which a fracture is initiated and propagates because of hydraulic loading applied by fluid inside the fracture. 3As a technique to fracture underground rock formations using pressurized fluid, hydraulic fracturing has been applied extensively in such diverse areas as reservoir stimulation, in situ stress estimation, caving and fault reaction in mining, and environmental subsurface remediation. 4Since the first field test was performed on gas well at the Hugoton field in 1947, hydraulic fracturing has been extensively researched by both academia and industries using experimental tests, field trials, and numerical simulations.With the continuous consumption of fossil energy, unconventional oil and gas resource has gradually become an important energy alternative. 5Hydraulic fracturing has developed greatly over the past few decades and has become one of the most widely utilized methods to increase the production of unconventional reservoirs by creating a new fracture network to improve the reservoirs' permeability.The fundamental principle of hydraulic fracturing in a coal seam is the high-pressure injection of fracturing fluid into cracks, including pre-existing cracks and artificially induced cracks.During the fracturing period, breakdown pressure is achieved, and the cracks are broadened, extended, and combined. 3ompared with conventional natural gas reservoirs, the biggest difference between coal reservoirs is that as a double-porosity rock formation, it has its unique cleat system in addition to the matrix pores of coal rocks. 6,7The face cleats and butt cleats in coal rock are intertwined, and they are roughly perpendicular to each other, forming effective seepage channels.The schematic diagram of the coal rock cleat system is shown in Figure 1.During the process of hydraulic fracturing, high-pressure fluid penetrates into both the matrix and the cleats along many artificially induced cracks, which not only extrude and release gas but also provide transportable passage for gas. 8n the process of hydraulic fracturing, the fracture morphology and propagation law of coal seam hydraulic fracturing are the key to affecting the effect of hydraulic fracturing and evaluating the efficiency of coalbed methane production.Domestic and foreign scholars' research on fracture propagation in deep rock mass mainly includes two aspects: indoor experiments and numerical simulations. 9In terms of indoor experiments, Lamont and Jessen 10 were the first to conduct fracture propagation experiments.Through field experiments, Warpinski 11 found that hydraulic fractures have three propagation modes: through natural fractures (NFs), captured by open NFs, and captured by shear-fractured NFs.Olson et al. 12 embedded picks of different sizes in gypsum rock samples as nonpermeable prefractures through experiments, and observed three propagation modes: bypass, pass, and divert.Chitrala et al. 13 measured the fracture propagation morphology under different horizontal stresses by acoustic emission device.Jin and colleagues 14,15 systematically studied the influence of NFs on the propagation of hydraulic fractures, and creatively considered the influence of NF dip angles on the intersection mode.Zhang et al. 16 designed an experimental model of tight sandstone with a closed cemented pre-existing fracture network to explore the influence of closed cemented NFs on the propagation behavior of hydraulic fracture in tight sandstone formations and reproduced the complex propagation patterns of hydraulic fractures in fractured tight sandstone formations.Li et al. 17 carried out a simulation experiment of intrafracture temporary plugging and diversion fracturing on natural carbonate rock outcrops and studied the influence of different factors on fracture morphology.In terms of numerical simulation, traditional methods include two-dimensional fracture propagation models such as the Kristianovich-Geertsma-de Klerk (KGD) model, the Perkins-Kern-Nordgren) model, the pseudo-3D (P3D) model, and the planar 3D model.However, if NFs are very developed in the reservoir, the propagation trajectory of hydraulic fractures will become very complicated, and conventional analytical models, P3D, planar, or multiplanar hydraulic fracture models will no longer be applicable.In recent years, to simulate the hydraulic fracture (HF) propagation behavior in complex unconventional oil and gas reservoirs, a variety of methods have been proposed by industry and academia, which can be mainly divided into two categories: (1) simulation method based on mechanics of continuous medium: finite element method (FEM), extended FEM, and boundary element method, and so on; (2) simulation methods based on mechanics of discontinuous medium: discrete element method (DEM), discontinuous deformation analysis), combined finite DEM, and so on.Because there are still many difficulties in realizing the complex formation environment of deep rock mass through experimental methods, while the laboratory test methods are expensive and have more influencing factors, it is difficult to effectively simulate the influence of NFs on hydraulic fracture propagation behavior at the reservoir scale. 18Therefore, the numerical simulation of hydraulic fracture propagation has been more widely used.
The cohesive zone model (CZM) based on the ABAQUS platform can well simulate fracture initiation and propagation and has been widely used in hydraulic fracturing simulation 19 CZM is usually defined in terms of the traction and separation of the fracture interface using a special traction-separation (T-S) constitutive law, which can effectively avoid the stress singularity at the fracture tip in traditional linear elastic fracture mechanics.Therefore, it provides an alternative method for describing the fracturing process.However, the conventional CZM needs to set the fracture propagation path in advance during the modeling process, which cannot realize the fracture's random propagation.Chen 20 proposed two methods of embedding cohesive elements: shared nodes and node binding and the results of hydraulic fractures calculated by the CZM are compared with the analytical solutions, which verifies the good accuracy of the CZM in solving hydraulic fractures.Guo et al. 21and Zhang et al. 22 both compared the injection pressure with the field injection pressure through CZM simulation and proved that the numerical solution is in good agreement with the field pressure data.Sobhaniaragh et al. 23 compared the simulation results based on the classical KGD analytical solution to verify the reliability of CZM for evaluating fracture geometry.Wang et al. 24 proposed a merging method of pore pressure nodes to realize the random propagation process of hydraulic fractures by globally embedding cohesive elements.
Therefore, aiming at the characteristics of natural cleat fractures development in coal reservoirs and strong heterogeneity of deep coal seams, this study proposes a hydraulic fracture propagation simulation method using the ABAQUS Software for natural cleat fracture networks and establishes a two-dimensional hydraulic fracture propagation model of heterogeneous coal-mud combination.An innovative finite element meshing scheme was used to model hydraulic fracturing based on the cohesive element method, featuring global cohesive element embedding and global pore pressure node sharing techniques.The natural cleat fracture network characteristics, the influence of cleat angle, cleat density, length ratio of coal rock and mudstone combination, and construction displacement on hydraulic fracture propagation behavior of deep coal seams were investigated.This study aims to deepen the understanding of the mechanism of hydraulic fracture propagation in deep coal seams.The results of our research have important practical guiding significance for accurately describing the random propagation behavior of complex hydraulic fractures in deep heterogeneous coal-mud combinations and guiding deep coal seam fracturing and hydraulic fracturing design.

| MATHEMATIC MODEL 2.1 | Fluid-solid coupling equation
In the model proposed in this work, a coupled stress-seepage-damage field approach is used to simulate the elastic deformation of the rock matrix and fluid seepage process.To simplify the calculation of the mathematical model, the composition of the rock is divided into two parts: the solid matrix and the rock pores. 25In this paper, the reservoir is considered the porous medium, and the stress balance equation is expressed in the form of a volumetric virtual work principle.The rock stress balance equation is as follows 26 : where V is the volume (m 3 ), S is the area (m 2 ), σ is the Cauchy effective stress (Pa), p w is the rock matrix pore pressure (Pa), I is the unit matrix, δε is the virtual strain rate vector (s −1 ), t is the surface force vector (N/m 2 ), δν is the virtual velocity vector (m/s), and f is the volumetric force vector (N/m 3 ).
According to the principle of mass conservation, the mass of fluid flowing into the porous medium rock within a certain period is equal to the sum of the increase of internal fluid and the outflow of fluid.The continuity equation describing fluid seepage is as follows 26 : where J is the formation volume change rate (dimensionless), ρ w is the fluid density (kg/m 3 ), n w is poros- ity (dimensionless), ν w is fluid seepage velocity (m/s), x is space vector (m), and t is time (s).
TIAN ET AL.
| 3561 Assuming that the fluid flow in the rock satisfies Darcy's law 26 where k is the permeability vector (m/s) and g is the acceleration due to gravity (m/s 2 ).

| Criteria for hydraulic fracture initiation and propagation
CZM simulates the initiation and propagation of hydraulic fractures through the T-S criterion of stiffness degradation. 27Figure 2 shows the bilinear constitutive relationship curve obeyed by the damage and failure processes of the cohesive element.As illustrated in Figure 2, where point A represents the initial damage location, which corresponds to the displacement δ 0 and maximum traction T max at the initial damage.Where point B represents the complete damage location, which corresponds to the effective displacement δ f at complete damage.
Before the damage initiation of the cohesive element, the traction force increases with the increase of traction displacement.The relationship between the nominal stress and nominal strain of the cohesive element is expressed as follows 28 : where t is the stress vector (Pa), K is the stiffness matrix, and ε is the strain.
When the damage initiation point is reached, the material begins to be damaged.In this paper, the maximum nominal stress criterion is used to judge whether the initial damage of the cohesive element occurs; it is assumed that when the stress generated in any direction reaches the critical stress set by the cohesive element, the cohesive element begins to form damage and the HF initiation begins.This criterion can be expressed as 28 where t n represents the nominal normal traction force applied on the cohesive element (Pa), t s and t t are, respectively, the tangential traction stress applied in two directions of the cohesive element, t n 0 is the critical normal traction force for the cohesive element to fail (Pa), t s 0 and t t 0 are, respectively, the critical tangential traction force of the cohesive element in two directions (Pa), and symbol  indicates that the cohesive element can only bear tensile stress and will not be damaged when subjected to pure compressive stress compression deformation.When the fracture standard f reaches   f f 1.0 1.0 + tol in the tolerance range (in this paper, f = 0.05 tol ), the damage starts to fracture.
The fracture propagation process is described by the stiffness decay rate D of the cohesive element, it can be expressed as 28 where t n , t s , and t t , respectively, are the cohesive element normal, first tangential, and second tangential stresses predicted under the current strain according to the undamaged frontal elastic T-S criterion (Pa).t n , t s , and t t , respectively, correspond to the stress in the three directions (Pa).
In the above formula, the dimensionless damage factor D represents the overall damage of the cohesive element, and the initial value of D is 0, which changes with the fracturing time.When the damage starts, the D value increases linearly from 0 to 1, and when D is 1, it means a complete failure and the fracture completely expands.In the model, the overall damage variable D is determined by the T-S law and is calculated as 28 where δ m is the current displacement of the cohesive element (m).

| (Flow equation of fracturing fluid in fractures
When the cohesive element is completely damaged, a fracture is formed.The fluid flow in the fracture includes tangential and normal flow. where q is the volume flow vector in the cohesive element (m 3 /s), w is the fracture opening (m); μ is the viscosity of the fracturing fluid (Pa s), p is the fluid pressure gradient within the fracture (Pa/m).
The normal flow can be described by the leak-off coefficient of the cohesive element 29 q c p p q c p p where q t and q b are the normal flow velocity on the top and bottom surfaces of the fracture, respectively (m/ s), c t and c b are the leak-off coefficients of fluid on the top and bottom surfaces of the fracture, respectively (m/Pa s), p t and p b are the pore pressure of the fluid on the top and bottom surfaces of the fracture, respectively (Pa); p i is the pore pressure of the fluid on the middle surface of the fracture.Among them, the fluid in the fracture also needs to satisfy the mass conservation equation 29 where Q t ( ) is the injection rate of fracturing fluid (m 3 / s) and δ x y ( , ) is the Dirac delta function.Because the simulation size is small and the fracturing time is very short, the fracturing fluid leakoff can often be ignored.

| NUMERICAL MODEL
In this study, the finite element analysis software ABAQUS was used to establish a two-dimensional heterogeneous combination of coal rock and mudstone numerical model.The geometry of the model is shown in Figure 4.All the basic parameters were set according to the settings listed in Table 1.The model consists of two parts, the upper layer is the coal seam and the lower layer is the mudstone layer.The liquid injection point of the model is at the center of the model, and a 0.5-m-long perforation is preset above and below the injection point.Among them, cohesive elements are used to prefabricate natural cleat fractures in the upper coal seam.By establishing parallel and equally spaced face cleat lines in the model, the adjacent face cleat lines are connected by the butt cleat lines to finally form an orthogonal coal seam cleat fracture network.The lower mudstone layer maintains a compact homogeneous state.Constant pore pressure and constant displacement boundary conditions, namely, zero displacement and hydrostatic pressure boundaries, are imposed on the model boundary.The initial stress is set to hydrostatic pressure, and there is no stress difference. 27The whole simulation process consists of two steps and the calculation time is 9 s.They are the in situ stress balance stage and the fracturing fluid injection pressure holding and fracture propagation stage. 18Among them, the first second is the in situ stress balance stage and the injection time is 8 s.After the fracturing fluid is injected into the formation with a constant displacement, the pressure will be suppressed and the HFs will continue to expand outward.
To improve the accuracy and convergence of model calculation, the matrix element of reservoir rock is characterized by CPE4P (four-node quadrilateral element considering pore pressure).The global average mesh size of the model is 0.25 m, and the quadrilateral grid and free mapping method are used to divide the reservoir model.To simulate the random propagation of fractures, cohesive elements are embedded in all grids and natural cleat fractures so that all cohesive elements at the intersection share seepage nodes, and the fluid continuity at the intersection is ensured by sharing the nodes to realize the random propagation process of HFs, whose element type is COH2D4P (a four-node 2D cohesive pore pressure element considering seepage).The schematic diagram of the element grid embedded in the global cohesive element is shown in Figure 5.A solid element and its neighboring cohesive element are connected by sharing two common nodes, and two adjacent cohesive elements connect to each other by sharing one node.After the cohesive elements are embedded globally, the model has a total of 7291 CPE4P elements and 14,426 COH2D4P elements.The effectiveness of cohesive elements can be verified by simulating the interaction of HFs and NFs. 30,31To verify the numerical accuracy of the simulation with globally embedded cohesive elements, Blanton's 32 experiments were used for verification.Based on the experimental results of hydraulic fracturing under large triaxial stress in the laboratory, Blanton believes that the extension mode of HF and NF after interaction and intersection under different horizontal stress differences (HSDs) and different angles is that HF stops and passes through and proposed a criterion to predict the interaction between hydraulic fracture and NF.According to the criterion, a curve can be constructed by plotting the approaching angle (between HF and NF) and HSD, as shown in Figure 9.
To verify the presented model, referring to Blanton's experiment, this study established a finite element model of the interaction between HF and prefabricated NF under different HSDs and different angles based on cohesive elements, and merged seepage nodes at the intersection of HF and NF, so that all cohesive elements shared seepage nodes.The basic model is shown in Figure 6 and the main calculation parameters in the model are shown in Table 2.The parameters used are mainly from the research of Wang 33 and Xiong et al. 34 Figure 7 shows the intersection results of HF and NF under the condition that the HSD is 4 MPa and the approximate angle is 30°and 60°, respectively.Figure 8 shows two examples of the simulation results at different HSD.
Table 3 lists the results of the numerical and Blanton's experiments.The numerical simulation results are consistent with Blanton's experimental results.
To fully demonstrate the correctness of the global CZM, we conducted multiple sets of numerical simulations.The comparison between the simulation results and the experimental results of Blanton is shown in Figure 9.As can be seen, the simulation results agree very well with what Blanton's experimental criterion predicts, which verifies the accuracy of the numerical simulation method in this study.| 3565 ANALYSIS This section is divided into four cases to discuss and analyze the influence of coal seam cleat angle, cleat density, length ratio of coal rock and mudstone combination, and construction displacement on hydraulic fracture propagation behavior.

| Influence of cleat angle
To analyze the influence of different cleat angles on the coal seam fracture propagation, two sets of the most typical orthogonal cleat fracture network models were established in this study to characterize the distribution of coal seam cleat fractures.Among them, one set is the oblique cleat model, the angle between the face cleat and the horizontal x-axis direction is 45°, and the angle between the butt cleat and the horizontal x-axis direction is 135°, both of which are perpendicular to each other.The other set is the straight cleat model, the face cleats are distributed along the horizontal direction, and the butt cleats are distributed along the vertical direction, both of which are perpendicular to each other.The distribution cloud images of the fracture openings of the two models at different moments (the fracture openings are enlarged by 100 times in the figure) are depicted in Figure 10.In the oblique cleat model of the coal seam, the coal seam fractures spread evenly along the face cleats and butt cleats on the basis of the main fractures formed at the perforation and develop a tortuous branch fracture at the intersection of the face cleats and butt cleats.Finally, the propagation mode, with main fractures as the main and micro-fractures as supplementary forms, is formed until it passes through the coal rock layer and enters the mudstone layer, and expands vertically outward along the direction of the main fractures.In the straight cleat model, coal seam fractures pass through the face cleats and extend outward along the butt cleats until they pass through the interface between coal rock and mudstone and enter the mudstone layer, finally forming a long straight main fracture.From the comparison of the two, it can be seen that the HFs in tight mudstone formations have a single form, which is one main fracture, and no obvious deflection occurs.However, different orthogonal cleat angles in coal seams have a greater impact on the fracture propagation form.Compared with the orthogonal cleat fractures in the straight cleat model, the orthogonal cleat fractures in the oblique cleat model greatly increase the complexity of fractures after fracturing.The cleat fractures open rapidly under the impact of fracturing fluid, and then continuously connect new natural cleat fractures.It is easier to form a complex fracture network and enhance the regional connectivity of the reservoir.
As shown in Figure 11, the injection pressure rises rapidly before the coal rock breakdown at the initial stage of fracture propagation.At about 2 s, the coal rock in the two models ruptures almost at the same time, but the rock fracture pressure of the oblique cleat model is significantly higher than that of the straight cleat model.In the later stage, when HFs pass through the coal seam and enter the mudstone layer, the pore pressure at the two injection points tends to be gentle and maintains the same variation.Therefore, compared with the orthogonal cleat fractures in the straight cleat model, the construction of coal seams with oblique cleat model cleat fractures requires greater pump pressure.
The variation curve of total fracture length and maximum fracture width under different cleat angles is shown in Figure 12.When the same injection process is used, the maximum fracture opening at the injection point and the total fracture length of the oblique cleat model are larger than those of the straight cleat model.Under this calculation time, the total fracture length and maximum fracture width of the oblique cleat model are 29.09m and 9.77 mm, respectively; the total fracture length and maximum fracture width of the straight cleat model are 15.44 m and 8.74 mm, respectively.It fully indicates that compared with the straight cleat model, the orthogonal cleat fractures in the oblique cleat model greatly increase the complexity of fractures after fracturing, and the synergistic effect on HFs is also more significant.

| Influence of cleat density
The development degree of coal rock cleats affects the fracture propagation morphology and is an important factor in controlling the extension scale of hydraulic fracturing fractures.To analyze the influence of orthogonal cleat fracture systems with different development degrees on coal seam fracture propagation, the cleat density is defined as the number of cleats within a certain measurement distance.Under the condition of ensuring the connectivity of cleats, the cleat density is characterized by changing the number of face cleats in the geometric model, which in turn reflects the development degree of cleats.In this study, two sets of basic models were established, namely, the high-and the low-density model.The high-density model sets the distance between adjacent face cleats to be 1.5 m, and the number of face cleats is 16.In the low-density model, the distance between adjacent face cleats is set to 3 m, corresponding to eight face cleats in the model.Figure 13 depicts the simulated hydraulic fracture propagation morphology under different cleat densities.
In the high-density model, the hydraulic fracture expands uniformly along the face cleats and butt cleats of the coal seam and develops a tortuous branch fracture at the intersection of the face cleats and butt cleats.Finally, a propagation mode combining the main fracture and the associated secondary fracture is formed until the hydraulic fracture passes through the interface along the main fracture direction into the mudstone layer.In the low-density model, most of the HFs propagate forward along the face cleats and a few extend along the butt cleats, and bifurcate at the junction of the face cleats and the butt cleats to form two obvious branch fractures.Finally, it extends vertically into the mudstone layer in the form of double fractures.The comparison indicates that the denser the coal rock cleat system, the higher its development degree, and the more uniform the propagation of HFs along the face cleats and butt cleats.
As shown in Figure 14, the fracture pressure of the high-density model is significantly greater than that of the low-density model at the initial stage of fracture propagation.As the fluid injection pressure reaches the fracture pressure of the coal seam, the formation begins to rupture and the fracture gradually extends forward.When the injection time is about 3 s, a branch fracture is formed in the lower part of the low-density model of the The variation curve of total fracture length (left) and maximum fracture width (right) under different cleat angles.coal seam, which is larger and more obvious than the high-density model.For branch fractures, more energy needs to be accumulated for rupture, and higher injection pressure is required.Therefore, the injection pressure of the low-density model begins to exceed that of the high-density model and remains until the later stage of fracture propagation.
The variation curve of total fracture length and maximum fracture width under different cleat densities is shown in Figure 15.In coal seams with different cleat densities, using the same injection process, the maximum fracture width at the injection point of the highdensity model is larger than that of the low-density model at the initial stage of fracture propagation, but with the formation of branch fractures in the low-density model, the required injection pressure is higher, and the pore pressure of the liquid injection point is increasing continuously, so the maximum fracture width of the liquid injection point is also increasing, gradually exceeding the maximum fracture width of the highdensity model.The changing trend is consistent with the change in the injection pressure curve.With the continuous injection of fluid, the high-density model only develops one branch fracture, while the low-density model develops two branch fractures.Therefore, in the later stage of fracture propagation, the total fracture length of the low-density model exceeds that of the high-density model.The model calculation results show that under the same calculation conditions, the total fracture length and maximum fracture width of the highdensity model are 29.09m and 9.77 mm, respectively, and the total fracture length and maximum fracture width of the low-density model are 30.12m and 10.66 mm, respectively.

| Influence of length ratio of coal rock and mudstone combination
To analyze the influence of the length ratio of the coal rock and mudstone combination on the crack propagation, the length ratio coefficient is defined as N.By changing the length of the coal rock in the geometric model to 5 m, 10 m, and 15 m, the N values correspond to 0.33, 1, 3. Figure 16 depicts the simulated hydraulic fracture propagation morphology under different length ratios of coal rock and mudstone combination.When N = 0.33 and the length ratio of the coal rock and mudstone combination is 1:3, the hydraulic fracture mainly propagates in the homogeneous mudstone layer to form a long straight main fracture.After entering the coal seam, the hydraulic fracture began to expand along the face cleat and the butt cleat and finally formed two tortuous branch fractures in the coal seam.When N = 1 and the length ratio of the coal rock and mudstone combination is 1:1, the hydraulic fracture expands synchronously at the interface of coal rock and mudstone.In mudstone, the fracture continues to expand vertically downward along the direction of the main fracture.In coal rock, the hydraulic fracture bifurcates obviously and three branch fractures are formed at the junction of face cleats and butt cleats.When N = 3 and the length ratio of the coal rock and mudstone combination is 3:1, the hydraulic fracture mainly expands evenly along the face cleats and butt cleats in the coal seams containing orthogonal cleat fractures and forms a branch fracture.Finally, it expands into the mudstone and expands vertically downward along the main seam direction.It can be seen from the comparison that the length ratio of the coal rock and mudstone combination has a greater influence on fracture propagation when the same injection process is used.Therefore, in the actual fracturing construction of the coal-mud interaction reservoir, the coal-mud combination with different length ratios needs to be matched with  different fracturing conditions to meet the field requirements.
Figure 17 shows the injection pressure curve under different length ratios of coal rock and mudstone combination.Due to the influence of formation heterogeneity, the tensile strength of coal rock is much lower than that of mudstone.When N = 0.33, the HFs are mainly concentrated in mudstone, and the fracture pressure is the highest; when N = 3, the HFs mainly propagate in coal rocks, and the fracture pressure is in the middle; when N = 1, the HFs propagate synchronously at the interface between coal rocks and mudstones, and the fracture pressure is the lowest under the joint influence of rock heterogeneity in different formations.With the continuous injection of fracturing fluid, it can be seen that in the later stage of fracture propagation, when N = 0.33, the pore pressure at the injection point is the highest, and when N = 3, the pressure at the injection point is the lowest.This shows that due to the existence of coal rock cleat fractures, the mechanical properties of the rock formation are relatively weaker than those of the mudstone formation, and the injection pressure required for rock fracture in the later stage of fracture propagation is lower.
The variation curve of total fracture length and maximum fracture width under different length ratios of coal rock and mudstone combination is shown in Figure 18.With the same injection process, due to the difference in the mechanical properties of different rocks, the maximum fracture width at the injection point is significantly different under the condition of different length ratios.Fractures are easier to expand in coal rock than in mudstone, so the maximum fracture width at the fracture injection point when expanding in the coal seam is larger than that when expanding in the mudstone.As the fluid is injected, the fracture continues to expand.The maximum fracture width is the smallest when N = 0.33 and the maximum fracture width is the largest when N = 3.As shown in Figure 18, the orthogonal cleat fractures of the coal seam have a significant synergistic effect on the HFs.When N = 0.33, the HFs mainly form long and straight main fractures in the mudstone, so the total length of the fractures is the smallest; when N = 1, the HFs not only expand vertically in the mudstone but also form three obvious branch fractures in the coal seam cleat fracture system; when N = 3, the HFs only forms a tortuous branch fracture in the coal seam cleat fracture system, so the total length of the fracture at N = 1 is slightly larger than that at N = 3.The model calculation results show that under the same calculation conditions, the total fracture length and maximum fracture width of N = 0.

| Influence of construction displacement
In actual fracturing construction, the construction displacement of fracturing fluid has a significant impact on the propagation pattern of the fracture network.Different construction displacements can provide different net pressures for the propagation and extension of fractures, thus affecting the initiation and expansion of fractures.To analyze the influence of construction displacement on fracture propagation and to ensure that the total displacement injected into the formation remains unchanged at 0.06 m 3 , the fracturing fluid is injected into the perforation at peak displacements of 0.002, 0.006, and 0.01 m 3 /s, respectively.Figure 19 depicts the simulated hydraulic fracture propagation morphology under different construction displacements.
With the increase of injection displacement, hydraulic fractures continuously activate NFs opening and expand more uniformly along coal seam face cleats and butt cleats, and the number of communicating natural cleat fractures also gradually increases.At the same time, a higher displacement is conducive to the formation of secondary fractures and promotes the bifurcation of fractures.
Figure 20 shows the curve of injection pressure under different construction displacements.As the injection displacement increases, the primary fracture pressure of the coal seam increases significantly, and the larger the injection displacement, the faster the formation fracture pressure decreases.This shows that increasing the displacement can significantly improve the ability to disturb the formation.With the increase of the displacement, the disturbance effect of the fracturing fluid on the coal rock near the injection point is continuously enhanced.In addition, the large displacement accelerates the flow rate of fracturing fluid in the HF, and the pressure in the fracture increases rapidly until it reaches the formation fracture pressure and then decreases rapidly.
The variation curve of total fracture length and maximum fracture width under different construction displacements is shown in  and the total length of the HF increase, indicating that the high displacement increases the fluid velocity and fluid energy in the HF, which is conducive to promoting hydraulic fractures continue to expand outward.When the liquid injection volume is the same, ensure that the total displacement of the injected formation fluid remains constant at 0.06 m 3 .When the injection rate is 0.002 m 3 /s, the total injection time is 30 s, and the total fracture length and maximum fracture width are calculated to be 22.83 and 7.66 mm, respectively.When the injection rate is 0.006 m 3 /s, the total injection time is 10 s, and the calculated total fracture length and maximum fracture width are 24.37 and 7.92 mm, respectively.When the injection rate is 0.01 m 3 /s, the total injection time is 6 s, and the total fracture length and maximum fracture width are calculated to be 25.24 and 9.21 mm, respectively.It shows that under the condition of sufficient liquid injection, the larger the displacement, the easier the NF to open, and the easier the cleat fracture in the coal seam to penetrate each other.Therefore, the large displacement is beneficial to improve the effect of hydraulic fracturing and achieve the purpose of creating long, wide, and complex fractures.

| CONCLUSIONS
Aiming at the characteristics of strong heterogeneity and natural cleat fracture development in deep coal reservoirs, this paper established a finite element model for hydraulic fracture propagation of a two-dimensional heterogeneous combination of coal rock and mudstone based on the global embedded CZM, and Blanton's 32 experiment criteria verified the accuracy of the global embedded CZM, the simulation results are in good agreement with the experimental results.Finally, four cases were used to discuss the influence of different coal seam cleat angles, cleat density, length ratio of coal rock and mudstone combination, and construction displacement on hydraulic fracture propagation behavior.The results show the following: 1. Compared with the cleat fractures in the straight cleat model, the orthogonal cleats in the oblique cleat model greatly increase the complexity of the fractures after fracturing, making it easier to form complex fracture networks, and the synergistic effect on hydraulic fractures is also more significant.2. The higher the cleat fracture development degree of coal rock is and the denser the cleat density is, the more uniform the HFs are along the face cleats and butt cleats.3. Affected by different rock heterogeneity and coal seam cleat fractures, the length ratio of coal rock and mudstone combination has a greater influence on fracture propagation.Therefore, in the actual fracturing construction, coal rock and mudstone combination with different length ratios need to be matched with different fracturing techniques to meet the field requirements.4. Increasing the displacement can significantly improve the disturbance ability of the formation.Under the condition of sufficient liquid injection, a large displacement is conducive to improving the effect of hydraulic fracturing and achieving the purpose of creating long, wide, and complex fractures.
Figure 3 illustrates the fluid flow within cohesive elements.In the fracture zone, the fluid flows in the way of tangential flow, while in the cohesive zone, the fluid flows in the way of normal flow.Tangential flow is the reason for the forward expansion of HFs, and normal flow causes the filtration of fluid in the HF.This study describes the tangential flow of fracturing fluid in HFs as continuous incompressible Newtonian fluid flow, controlled by the lubrication equation derived from Poise Tuille's law 29 :

T A B L E 1 14 F
Parameters of the numerical model.I G U R E 5 Schematic diagram of the element grid embedded in the global cohesive element.(A) 6-node two-dimensional displacement and pore pressure cohesive elements of 0 thickness and (B) Global embedding cohesive element node connectivity diagram.

F I G U R E 7
Abbreviation: HSD, horizontal stress difference.

10
Fracture propagation morphology under different cleat angles.F I G U R E 11 The injection pressure curve under different cleat angles.F I G U R E 9 Comparison between simulation results and Blanton's experiment criterion.HSD, horizontal stress difference.

F
I G U R E 13 Fracture propagation morphology under different cleat densities.F I G U R E 14 The injection pressure curve under different cleat densities.

F
I G U R E 16 Fracture propagation morphology under different length ratios of coal rock and mudstone combination.

F
I G U R E 15 The variation curve of total fracture length (left) and maximum fracture width (right) under different cleat densities.
33 are 19.05 m and 7.13 mm, respectively; the total length and maximum fracture width of N = 1 are 26.31m and 8.57 mm, respectively; the total length and maximum width of N = 3 are 25.24 m and 9.21 mm, respectively.F I G U R E 17 The injection pressure curve under different length ratios of coal rock and mudstone combination.F I G U R E 18 The variation curve of total fracture length (left) and maximum fracture width (right) under different length ratios of coal rock and mudstone combination.

Figure 21 .
When the injection time is the same, the injection displacement increases, and the maximum fracture width at the injection point F I G U R E 19 Fracture propagation morphology under different construction displacements.

F
I G U R E 20 The injection pressure curve under different construction displacements.