Numerical Simulation of Hydraulic Fracture Propagation in Coal Seams with Discontinuous Natural Fracture Networks

To investigate the mechanism of hydraulic fracture propagation in coal seams with discontinuous natural fractures, an innovative finite element meshing scheme for modeling hydraulic fracturing was proposed. Hydraulic fracture propagation and interaction with discontinuous natural fracture networks in coal seams were modeled based on the cohesive element method. The hydraulic fracture network characteristics, the growth process of the secondary hydraulic fractures, the pore pressure distribution and the variation of bottomhole pressure were analyzed. The improved cohesive element method, which considers the leak-off and seepage behaviors of fracturing liquid, is capable of modeling hydraulic fracturing in naturally fractured formations. The results indicate that under high stress difference conditions, the hydraulic fracture network is spindle-shaped, and shows a multi-level branch structure. The ratio of secondary fracture total length to main fracture total length was 2.11~3.62, suggesting that the secondary fractures are an important part of the hydraulic fracture network in coal seams. In deep coal seams, the break pressure of discontinuous natural fractures mainly depends on the in-situ stress field and the direction of natural fractures. The mechanism of hydraulic fracture propagation in deep coal seams is significantly different from that in hard and tight rock layers. Record Type: Published Article Submitted To: LAPSE (Living Archive for Process Systems Engineering) Citation (overall record, always the latest version): LAPSE:2018.0417 Citation (this specific file, latest version): LAPSE:2018.0417-1 Citation (this specific file, this version): LAPSE:2018.0417-1v1 DOI of Published Version: https://doi.org/10.3390/pr6080113 License: Creative Commons Attribution 4.0 International (CC BY 4.0) Powered by TCPDF (www.tcpdf.org)


Introduction
Coal seam gas that mainly contains methane is an economical and promising solution for the world's energy crisis [1]. Global coalbed methane (CBM) resources are estimated to range from 84.38 Tm 3 to 262.21 Tm 3 , of which approximately 14.44 Tm 3 are recoverable [2]. China has the world's third largest CBM reserve, totaling approximately 37 Tm 3 [3,4]; however, as they are influenced by coalification, geological structure and depth, most Chinese CBM reservoirs have low permeability [5][6][7], which makes gas extraction difficult, and even leads to coal and gas outbursts [8].
Hydraulic fracturing is an efficient stimulation method that is widely used to enhance CBM production [9][10][11][12][13] by creating a new fracture network to improve the reservoirs' permeability and by disturbing the gas in its absorbed state to promote gas desorption from coal [9,14]. Moreover, hydraulic fracturing has also been applied to destress coal seams to prevent coal and gas outbursts. Other significant applications of hydraulic fracturing include the measurement of geo-stress [15,16], This study proposes a hydraulic fracture propagation simulation method using the ABAQUS Software for discontinuous natural fracture networks and a preliminary analysis of hydraulic fracture morphology in coal seams. 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. A simulation methodology including leak-off and seepage effects was developed. The hydraulic fracture network characteristics, the growth mechanism of secondary hydraulic fractures, the pore pressure distribution in coal seams, and the variation of fracturing liquid pressure were investigated. Overall, this study aims to deepen the understanding of the mechanism of hydraulic fracture propagation in coal seams.

Concept and Methodology of the Hydraulic Fracturing Simulation Using the Cohesive Element Method
The cohesive element method is based on damage mechanics and the cohesive zone model (CZM). By embedding zero-thickness cohesive elements into a mesh, dynamic and unrestricted fracture propagation (especially crack branching, coalescence and swerve) can be modeled. Figure 1 illustrates the process of generating zero-thickness cohesive elements into a solid element mesh. This algorithm was coded with the Python programing language and implemented as an ABAQUS plugin to generate cohesive elements. By executing the plugin, cohesive elements are generated at every solid element interface. Therefore, the final model consists of solid elements and cohesive elements. The solid element and its neighboring cohesive element are connected by their sharing two common nodes, and the two adjacent cohesive elements connect to each other by sharing one node, as shown in Figure 1. Under internal/external traction loads, cohesive elements experience loading, damage, stiffness degradation and cracking. Crack initiation and propagation are represented by the damage and failure of cohesive elements. The coal/rock micro-elastic-plastic deformation and seepage behaviors are dominated by solid elements.
Processes 2018, 6, x FOR PEER REVIEW 4 of 25 hydraulic fracturing based on the cohesive element method, featuring global cohesive element embedding and global pore pressure node sharing techniques. A simulation methodology including leak-off and seepage effects was developed. The hydraulic fracture network characteristics, the growth mechanism of secondary hydraulic fractures, the pore pressure distribution in coal seams, and the variation of fracturing liquid pressure were investigated. Overall, this study aims to deepen the understanding of the mechanism of hydraulic fracture propagation in coal seams.

Concept and Methodology of the Hydraulic Fracturing Simulation Using the Cohesive Element Method
The cohesive element method is based on damage mechanics and the cohesive zone model (CZM). By embedding zero-thickness cohesive elements into a mesh, dynamic and unrestricted fracture propagation (especially crack branching, coalescence and swerve) can be modeled. Figure 1 illustrates the process of generating zero-thickness cohesive elements into a solid element mesh. This algorithm was coded with the Python programing language and implemented as an ABAQUS plugin to generate cohesive elements. By executing the plugin, cohesive elements are generated at every solid element interface. Therefore, the final model consists of solid elements and cohesive elements. The solid element and its neighboring cohesive element are connected by their sharing two common nodes, and the two adjacent cohesive elements connect to each other by sharing one node, as shown in Figure 1. Under internal/external traction loads, cohesive elements experience loading, damage, stiffness degradation and cracking. Crack initiation and propagation are represented by the damage and failure of cohesive elements. The coal/rock micro-elastic-plastic deformation and seepage behaviors are dominated by solid elements. Note that 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.
In this study, the cohesive pore pressure element (CPPE), a special type of cohesive element, was used to simulate hydraulic fracture propagation and leak-off behavior from hydraulic fractures to the porous coal matrix. As illustrated in Figure 2a, a 2D zero-thickness CPPE has six nodes and two integration points. The connection of nodes 1 and 2 forms the bottom face of a cohesive element, and nodes 3 and 4 form the top face. The relative displacements between the top and bottom faces (both in 1-direction and 2-direction) represent crack opening and slipping. Therefore, Nodes 1 to 4 are used to calculate the relationship between tractions and separations. Note that 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.
In this study, the cohesive pore pressure element (CPPE), a special type of cohesive element, was used to simulate hydraulic fracture propagation and leak-off behavior from hydraulic fractures to the porous coal matrix. As illustrated in Figure 2a, a 2D zero-thickness CPPE has six nodes and two integration points. The connection of nodes 1 and 2 forms the bottom face of a cohesive element, and nodes 3 and 4 form the top face. The relative displacements between the top and bottom faces (both in 1-direction and 2-direction) represent crack opening and slipping. Therefore, Nodes 1 to 4 are used to calculate the relationship between tractions and separations.
Nodes 5 and 6 are the pore pressure nodes, to calculate the fracturing liquid flow from one CPPE to the next CPPE (tangential flow) and leak-off from cohesive elements to solid elements (normal flow). The CPPEs are zero-thick, so the pore pressure nodes of the discrete adjacent CPPEs occupy the same coordinates after generating the discrete CPPEs in a mesh, as shown in Figure 2a (pore pressure nodes p1, p2, p3, p4 have the same coordinates P (x, y)). The fracturing liquid flow and transmission need to go through the pore pressure nodes. All the adjacent CPPEs need to share one common pore pressure node to connect to each other and to allow the fracturing liquid to flow in the fracture from one CPPE to the next. Therefore, after creating the discrete CPPEs, the pore pressure nodes at the same coordinates were merged, as shown in Figure 2b. Nodes 5 and 6 are the pore pressure nodes, to calculate the fracturing liquid flow from one CPPE to the next CPPE (tangential flow) and leak-off from cohesive elements to solid elements (normal flow). The CPPEs are zero-thick, so the pore pressure nodes of the discrete adjacent CPPEs occupy the same coordinates after generating the discrete CPPEs in a mesh, as shown in Figure 2a (pore pressure nodes p1, p2, p3, p4 have the same coordinates P (x, y)). The fracturing liquid flow and transmission need to go through the pore pressure nodes. All the adjacent CPPEs need to share one common pore pressure node to connect to each other and to allow the fracturing liquid to flow in the fracture from one CPPE to the next. Therefore, after creating the discrete CPPEs, the pore pressure nodes at the same coordinates were merged, as shown in Figure 2b.
(a) Schematic of the CPPE and the mesh structure by embedding discrete CPPEs; (b) the mesh structure after merging the pore pressure nodes at the same coordinates and the concept of fracturing liquid flow in fractures. Note that the CPPEs are zero-thick, and the pore pressure nodes p1, p2, p3, p4 have the same coordinates P (x, y). The fracturing liquid flow has two components: the tangential flow from one cohesive element to the next cohesive element, and the normal flow from cohesive elements to solid elements.

Discontinuous Natural Fracture Networks
As discussed in the introduction, natural fractures in coal seams usually have obvious directionality, thus coal seams are divided into regular hexahedral masses by the discontinuous natural fracture network. In fact, the fracture direction varies widely. This study assumed that the angle between the two sets of natural fractures was 90°, and its angle bisector was parallel to the maximum principal stress, as shown in Figure 3. The discontinuous fracture network was created as follows: (1) A 2D plane model was created, then this plane was partitioned by two sets of lines that are orthogonal to each other. (2) This plane was meshed to generate solid elements.
(3) CPPEs were embedded in this mesh model and numbered. The serial numbers of the CPPEs that were on the two sets of partition lines into set A were picked. (4) A certain proportion of elements from set A was randomly selected. By assigning very low mechanical properties, these selected elements were used to represent the discontinuous fractures. The above process for the creation of a discontinuous fracture network was executed using a Python script program in ABAQUS.   Figure 2. (a) Schematic of the CPPE and the mesh structure by embedding discrete CPPEs; (b) the mesh structure after merging the pore pressure nodes at the same coordinates and the concept of fracturing liquid flow in fractures. Note that the CPPEs are zero-thick, and the pore pressure nodes p1, p2, p3, p4 have the same coordinates P (x, y). The fracturing liquid flow has two components: the tangential flow from one cohesive element to the next cohesive element, and the normal flow from cohesive elements to solid elements.

Discontinuous Natural Fracture Networks
As discussed in the introduction, natural fractures in coal seams usually have obvious directionality, thus coal seams are divided into regular hexahedral masses by the discontinuous natural fracture network. In fact, the fracture direction varies widely. This study assumed that the angle between the two sets of natural fractures was 90 • , and its angle bisector was parallel to the maximum principal stress, as shown in Figure 3. The discontinuous fracture network was created as follows: (1) A 2D plane model was created, then this plane was partitioned by two sets of lines that are orthogonal to each other. (2) This plane was meshed to generate solid elements.
(3) CPPEs were embedded in this mesh model and numbered. The serial numbers of the CPPEs that were on the two sets of partition lines into set A were picked. (4) A certain proportion of elements from set A was randomly selected. By assigning very low mechanical properties, these selected elements were used to represent the discontinuous fractures. The above process for the creation of a discontinuous fracture network was executed using a Python script program in ABAQUS.
Processes 2018, 6, x FOR PEER REVIEW 6 of 25 Figure 3. The process of creating the discontinuous fracture network.

Seepage-Stress Coupling Equation for the Coal Matrix
Coal matrix is a typical porous medium that contains solid coal materials, gas-filled pores, and fracturing leak-off liquid that permeates the coal matrix during fracturing. In the seepage-stress coupling analysis, the effective stress principle was used, given by [86]: where σ is the effective stress matrix at a point, Pa; σ is the total stress matrix, Pa;  is the Biot coefficient; w u is the pore pressure, Pa; I is a unit matrix.
The porosity of the coal matrix, n, is defined as: where v dV is the void volume, m 3 ; and dV is the total volume in the current configuration, m 3 .
The void ratio e is defined as: Saturation, s, is the ratio of the total free seepage volume of liquid and gas, w dV , to the void volume, and is defined as:

Discretized Equilibrium Equation
Based on the virtual work principle, the stress equilibrium equation [87] in the current configuration can be expressed as: (4) where  v is the virtual velocity field, m/s;  ε is the virtual rate of deformation tensor, s −1 ; t is the surface traction tensor, Pa; f is the body force except the seepage liquid weight, N/m 3 ;  w is the density of the seepage liquid, kg/m 3 ; g is the gravity acceleration, m/s 2 ; V is the integration domains of volume, m 3 ; S is the integration domains of area, m 2 . Note that the last term in Equation (5) is the weight of the seepage liquid.

Partition
Mesh Discontinuous fracture network

Seepage-Stress Coupling Equation for the Coal Matrix
Coal matrix is a typical porous medium that contains solid coal materials, gas-filled pores, and fracturing leak-off liquid that permeates the coal matrix during fracturing. In the seepage-stress coupling analysis, the effective stress principle was used, given by [86]: where σ is the effective stress matrix at a point, Pa; σ is the total stress matrix, Pa; χ is the Biot coefficient; u w is the pore pressure, Pa; I is a unit matrix. The porosity of the coal matrix, n, is defined as: where dV v is the void volume, m 3 ; and dV is the total volume in the current configuration, m 3 . The void ratio e is defined as: Saturation, s, is the ratio of the total free seepage volume of liquid and gas, dV w , to the void volume, and is defined as:

Discretized Equilibrium Equation
Based on the virtual work principle, the stress equilibrium equation [87] in the current configuration can be expressed as: where δv is the virtual velocity field, m/s; δε is the virtual rate of deformation tensor, s −1 ; t is the surface traction tensor, Pa; f is the body force except the seepage liquid weight, N/m 3 ; ρ w is the density of the seepage liquid, kg/m 3 ; g is the gravity acceleration, m/s 2 ; V is the integration domains of  (5) is the weight of the seepage liquid.
In finite element analysis, the equilibrium equation is discretized into a set of equations with interpolation functions. The virtual velocity field and the virtual rate of deformation are interpolated as: 6) where N N is the matrix of interpolation functions defined with respect to the material coordinates; and x is the unit coordinate vector, m. By substituting Equation (6) into Equation (5), the virtual work equation is discretized as Note that the term V β N : σdV is the internal force array, I N , and S N N · tdS + V N N · f dV + V snρ w N N · gdV is the external force array, P N . Therefore, the equilibrium Equation (7) can be expressed as:

Continuity Equation of Seepage
The porous coal matrix can be modeled approximately by attaching the solid element mesh to the solid materials. Thus, the liquid flows through the mesh. The time rate of the change of liquid mass in a control volume is: where R is the time rate of change of liquid mass in a control volume, kg/s; t is time, s; J is the rate of volume change of the porous coal matrix, defined as the ratio of volume in the current configuration, dV, to the volume in the reference configuration, dV 0 : The liquid mass which crosses a surface, S, and flows into the control volume per unit time is: (11) where v w is the liquid seepage velocity, m/s; n is the outward normal vector to S. Obviously, the liquid mass change per unit time, R, is equal to the liquid mass flowing into the control volume per unit time, F, that is [85,88]: According to the divergence theorem, Equation (12) can be transformed to an equivalent weak form: Equation (13) is the continuity equation of the seepage. In this study, the seepage in the porous coal matrix obeys Darcy's law, which is given by [89]: Processes 2018, 6, 113 8 of 25

Flow Equation of the Fracturing Liquid Flow in CPPE
The constitutive response of the fracturing liquid flow in the CPPE combines tangential flow within the fracture and normal flow across the fracture, as illustrated in Figure 4. The fracturing liquid is assumed to be incompressible Newtonian. The tangential flow within the fracture is governed by [74,90]: where q is the fracturing liquid flow rate per unit area, m/s; w is the opening width of the hydraulic fracture, m; µ is the viscosity of the fracturing liquid, Pa·s; ∇p w is the pressure gradient along the CPPE, Pa/m. The term w 3 12µ is equivalent to the permeability or the resistance to fluid flow in the fracture.

Flow Equation of the Fracturing Liquid Flow in CPPE
The constitutive response of the fracturing liquid flow in the CPPE combines tangential flow within the fracture and normal flow across the fracture, as illustrated in Figure 4. The fracturing liquid is assumed to be incompressible Newtonian. The tangential flow within the fracture is governed by [74,90]: where q is the fracturing liquid flow rate per unit area, m/s; w is the opening width of the hydraulic fracture, m; μ is the viscosity of the fracturing liquid, Pa·s;  w p is the pressure gradient along the CPPE, w is equivalent to the permeability or the resistance to fluid flow in the fracture. Normal flow, representing the leak-off behavior from fractures into porous coal matrix, is defined as [74,85]: where t q and b q are the leak-off flow rates in the normal direction from the mid-face into the top face and bottom face of CPPE respectively, m/s; t c and b c are the corresponding leak-off coefficients, m/Pa·s; m p is the mid-face pressure, Pa; t p and b p are the pore pressures on the top and bottom faces respectively, Pa. The leak-off coefficient can be interpreted as the permeability of the hydraulic fracture surfaces (i.e., the top and bottom faces of CPPE). According to Equation (16), the normal flow rate depends on the leak-off coefficient, and the pressure difference between the fracture (CPPE) and the porous coal matrix (solid element).

Damage Initiation and Evolution Law of Hydraulic Fractures
Crack initiation and propagation are governed by the traction-separation law. Figure 5 illustrates the fracture propagation driven by the fracturing liquid based on the cohesive element method. The traction-separation law defines the relationship between traction force, T, and the relative separation displacement between the top and bottom faces, δ. In the red cohesive zone in  Normal flow, representing the leak-off behavior from fractures into porous coal matrix, is defined as [74,85]: where q t and q b are the leak-off flow rates in the normal direction from the mid-face into the top face and bottom face of CPPE respectively, m/s; c t and c b are the corresponding leak-off coefficients, m/Pa·s; p m is the mid-face pressure, Pa; p t and p b are the pore pressures on the top and bottom faces respectively, Pa. The leak-off coefficient can be interpreted as the permeability of the hydraulic fracture surfaces (i.e., the top and bottom faces of CPPE). According to Equation (16), the normal flow rate depends on the leak-off coefficient, and the pressure difference between the fracture (CPPE) and the porous coal matrix (solid element).

Damage Initiation and Evolution Law of Hydraulic Fractures
Crack initiation and propagation are governed by the traction-separation law. Figure 5 illustrates the fracture propagation driven by the fracturing liquid based on the cohesive element method. The traction-separation law defines the relationship between traction force, T, and the relative separation displacement between the top and bottom faces, δ. In the red cohesive zone in Figure 5, the cohesive elements are in an elastic state or in an incomplete damage state. state, their tensile stiffness cannot recover. The concept of the crack tip was referenced from the literature [91]. The point of δ = δ f m is the material crack tip, while the point of δ = δ o m is the cohesive crack tip. The cohesive elements between the two crack tips are in the damaged state. The traction force, T, is the resultant force of hydraulic pressure in the cohesive element, pore pressure in the adjacent solid element, and geo-stress. Therefore, the hydraulic fracture is affected by fracturing pressure, the discontinuous natural fracture network, the seepage process in the coal matrix and geo-stress.
To model crack initiation and propagation in different materials using cohesive elements, various traction-separation models were developed. In this work, the irreversible bilinear model was used to simulate hydraulic fracture opening and closing behaviors. This model was also used by Chen [74] and Zhao et al. [85] to model hydraulic fracturing in rock. The model assumes that the mechanical response of cohesive elements satisfies a reversible linear elastic behavior prior to the initial damage point. The bearing capacity of cohesive elements decreases linearly in the damage evolution stage, and their stiffness cannot recover during unloading. The irreversible bilinear model was introduced as follows.  o m is the cohesive crack tip. The cohesive elements between the two crack tips are in the damaged state. The traction force, T, is the resultant force of hydraulic pressure in the cohesive element, pore pressure in the adjacent solid element, and geo-stress. Therefore, the hydraulic fracture is affected by fracturing pressure, the discontinuous natural fracture network, the seepage process in the coal matrix and geo-stress.
To model crack initiation and propagation in different materials using cohesive elements, various traction-separation models were developed. In this work, the irreversible bilinear model was used to simulate hydraulic fracture opening and closing behaviors. This model was also used by Chen [74] and Zhao et al. [85] to model hydraulic fracturing in rock. The model assumes that the mechanical response of cohesive elements satisfies a reversible linear elastic behavior prior to the initial damage point. The bearing capacity of cohesive elements decreases linearly in the damage evolution stage, and their stiffness cannot recover during unloading. The irreversible bilinear model was introduced as follows.

Damage Initiation
The tension component of traction force is mainly induced by hydraulic pressure, while the shear component of traction force may be caused by geo-stress differences and the local shear stress concentration around the natural fracture. Therefore, it is significant and necessary to consider the combined effects of the tension force and the shear force when modeling hydraulic fracture initiation and propagation in coal seams. Thus, a quadratic stress criterion that incorporates the tension and shear mixed modes was used to predict damage initiation occurrence. The 2D crack damage initiates when the stress state of the cohesive elements satisfies [33,85]: where  t and  s are tensile and shear strength, respectively, Pa; n t and s t represents normal traction (along the local 2-direction in Figure 2a) and shear traction (along the local 1-direction in Figure 2a), respectively, Pa.

Damage Evolution
A scalar damage variable, D, is defined to represent the overall damage degree in one cohesive element. D evolves from 0 to 1 monotonically upon further loading after damage initiation. Stress components are affected by the damage according to [92]:

Damage Initiation
The tension component of traction force is mainly induced by hydraulic pressure, while the shear component of traction force may be caused by geo-stress differences and the local shear stress concentration around the natural fracture. Therefore, it is significant and necessary to consider the combined effects of the tension force and the shear force when modeling hydraulic fracture initiation and propagation in coal seams. Thus, a quadratic stress criterion that incorporates the tension and shear mixed modes was used to predict damage initiation occurrence. The 2D crack damage initiates when the stress state of the cohesive elements satisfies [33,85]: where σ t and σ s are tensile and shear strength, respectively, Pa; t n and t s represents normal traction (along the local 2-direction in Figure 2a) and shear traction (along the local 1-direction in Figure 2a), respectively, Pa.

Damage Evolution
A scalar damage variable, D, is defined to represent the overall damage degree in one cohesive element. D evolves from 0 to 1 monotonically upon further loading after damage initiation. Stress components are affected by the damage according to [92]: where t n and t s are stress components predicted by the linear elastic law at the current separation displacement, Pa, as shown in Figure 6a. Note that in Equation (18), the tensile stiffness of cohesive elements is degraded at the damage stage, whereas the compressive stiffness is unchanged.
For the linear damage softening evolution, D is expressed as [33,85,93]: where δ max m is the maximum value of the relative separation displacement recorded during the loading history, m. As illustrated in Figure 6b, the loading history is the loading along the red path (path 1) or reloading via the green path (path 2 where n t and s t are stress components predicted by the linear elastic law at the current separation displacement, Pa, as shown in Figure 6a. Note that in Equation (18), the tensile stiffness of cohesive elements is degraded at the damage stage, whereas the compressive stiffness is unchanged. For the linear damage softening evolution, D is expressed as [33,85,93]: where  max m is the maximum value of the relative separation displacement recorded during the loading history, m. As illustrated in Figure 6b, the loading history is the loading along the red path The power law criterion of fracture energy is widely used to predict the complete damage state, given by [94]: where I G and II G are energy components of tension and shear respectively, Pa·m; IC G and IIC G are I-mode and II-mode fracture energies respectively, Pa·m. When I G and II G satisfy Equation (20), the damage variable D reaches 1.
Note that the stress state and the corresponding displacement at the initial damage point is determined by Equation (17) (19), and then the stress components are updated by substituting D into Equation (18). Hydraulic fracture propagation is simulated using this damage law. The power law criterion of fracture energy is widely used to predict the complete damage state, given by [94]: where G I and G I I are energy components of tension and shear respectively, Pa·m; G IC and G I IC are I-mode and II-mode fracture energies respectively, Pa·m. When G I and G I I satisfy Equation (20), the damage variable D reaches 1.
Note that the stress state and the corresponding displacement at the initial damage point is determined by Equation (17)  using Equation (19), and then the stress components are updated by substituting D into Equation (18). Hydraulic fracture propagation is simulated using this damage law.

The Mesh Model of Coal Seams with Discontinuous Natural Fracture Networks
By setting the out-of-plane dimension as 3 m, the 2D numerical coal seam model was equivalent to a 3D model with a thickness of 3 m. As shown in Figure 7, the model is a square with the dimensions 400 m 2 . To avoid the boundary effect, the discontinuous natural fracture network zone is at the center of the model, with a size of 10 m × 10 m. The two fracture sets are orthotropic with each other, and their bisector is parallel to the direction of the principle stress. The fracture spacing was set to 0.

Mechanical Properties of the Coal Matrix, Discontinuous Natural Fracture and Fracturing Liquid
Many studies show that hydraulic fracturing has significant effects on the permeability, leak-off coefficient, and porosity of coal seams [85,95−97]. According to Xu et al. [38], the dynamic permeability k, the dynamic porosity n and the dynamic leak-off coefficient c can be calculated as follows:

Mechanical Properties of the Coal Matrix, Discontinuous Natural Fracture and Fracturing Liquid
Many studies show that hydraulic fracturing has significant effects on the permeability, leak-off coefficient, and porosity of coal seams [85,[95][96][97]. According to Xu et al. [38], the dynamic permeability k, the dynamic porosity n and the dynamic leak-off coefficient c can be calculated as follows: Here, k 0 is the initial permeability, m/s; ∆p is the pressure difference between the fracture (CPPE) and the porous coal matrix (solid element), Pa; n 0 is the initial porosity; c 0 is the initial leak-off coefficient, m/Pa·s; C f is the pore compressive coefficient; and C f = (n − n 0 )/(n · ∆p), Pa −1 .
The user subroutine was used to model the dynamic variation of the coal seepage properties. In one numerical increment step, if ∆p is determined, n can be calculated using Equation (22), and then C f can be obtained. k and c are then calculated using Equations (21) and (23). Until this point, k, n and c were updated dynamically. In the next incremental step, these parameters were set as the updated initial values and substituted into the governing equations for the seepage to model the dynamic fracturing process.
The mechanical properties used in this numerical model are listed in Table 1. Note that the mechanical strengths of the discontinuous natural fractures are significantly weaker than those of the coal matrix.

Initial Conditions
In this study, the injection rate of the fracturing liquid via the injection node was 0.003 m 3 /s; the initial gas pore pressure in the coal seams was 5 MPa. To analyze the hydraulic fracture propagation and its interaction with natural fractures under different ground stress conditions, four different stress conditions were applied, as shown in Table 2.  Figure 8 illustrates the pre-processing and the solving procedures of the ABAQUS/Standard solver. The model first calibrates geostatic stress to a balance state and then computes hydraulic fracturing. For the two steps above, the procedure types were 'geostatic' and 'soils', respectively. The module 'soils' in ABAQUS allows users to model hydraulic fracturing and fluid flow through porous media. In the 'soils' module, the transient analysis method and the direct backward difference operator were used to optimize the accuracy of the continuity equation. The solution technique used was the full Newton method for nonlinear equilibrium equations. In pore liquid flow analysis, pore pressure, w u , was selected as a field index to evaluate the accuracy of the Jacobian matrix solution. According to the ABAQUS documentation [98], most nonlinear calculations are sufficiently accurate if the error of the largest field index residual is less than 0.005. In this study, the convergence criterion was also used: To guarantee the accuracy of the fracture propagation when using CPPE, the cohesive element size should be smaller than the cohesive zone length c d , which is determined by the material mechanical properties, and is calculated as [74]:

Solving Procedures, Convergence Criterion and Error Control
where E is the Young's modulus, Pa; v is the Poisson's ratio. By substituting the parameter values from Table 1 into Equation (25), c d = 0.182 m. The cohesive element size used in this study was 0.1 m, which is significantly less than dc. This suggests good error control for modeling hydraulic fracture propagation.

Hydraulic Fracture Network Characteristics
The CPPEs with fracturing liquid flowing were selected and marked as hydraulic fractures during the fracturing process. Figure 9 shows the hydraulic fracture network characteristics under different stress conditions. According to this figure, hydraulic fractures in coal seams are not just a single crack, but a complex fracture network. Micro-seismic monitoring in the field also confirms that hydraulic fractures in coal seams usually form a network of fracture branches using the horizontal well technique [99]. The solution technique used was the full Newton method for nonlinear equilibrium equations. In pore liquid flow analysis, pore pressure, u w , was selected as a field index to evaluate the accuracy of the Jacobian matrix solution. According to the ABAQUS documentation [98], most nonlinear calculations are sufficiently accurate if the error of the largest field index residual is less than 0.005. In this study, the convergence criterion was also used: where r u w max is the largest residual in the pore pressure balance equation, m 3 /s; q u w is the overall time-averaged value of a typical flux for u w until this point during this step, including the current increment, m 3 /s.
To guarantee the accuracy of the fracture propagation when using CPPE, the cohesive element size should be smaller than the cohesive zone length d c , which is determined by the material mechanical properties, and is calculated as [74]: where E is the Young's modulus, Pa; v is the Poisson's ratio. By substituting the parameter values from Table 1 into Equation (25), d c = 0.182 m. The cohesive element size used in this study was 0.1 m, which is significantly less than d c . This suggests good error control for modeling hydraulic fracture propagation.

Hydraulic Fracture Network Characteristics
The CPPEs with fracturing liquid flowing were selected and marked as hydraulic fractures during the fracturing process. Figure 9 shows the hydraulic fracture network characteristics under different stress conditions. According to this figure, hydraulic fractures in coal seams are not just a single crack, but a complex fracture network. Micro-seismic monitoring in the field also confirms that hydraulic fractures in coal seams usually form a network of fracture branches using the horizontal well technique [99].  Figure 9 shows that the hydraulic fracture network has a main fracture trunk with many fracture branches. However, the main fractures are not strictly defined. Generally, this usually refers to the more prominent fractures with a relatively large opening. In hydraulic fracturing, the main fractures  Figure 9 shows that the hydraulic fracture network has a main fracture trunk with many fracture branches. However, the main fractures are not strictly defined. Generally, this usually refers to the more prominent fractures with a relatively large opening. In hydraulic fracturing, the main fractures are the main liquid flow fracture channels. The fracture opening degree and the stiffness of the CPPEs are indicators of main hydraulic fractures. The main hydraulic fractures were identified by magnifying the hydraulic fracture opening 20 times, as shown in Figure 10. During the initial hydraulic fracturing stage, the total fracture volume is small, thus a reticulated secondary fracture structure is created around the injection node (i.e., the injection borehole) when subject to the impact of high-pressure fracturing liquid, as shown in Figure 9. Then the total volume During the initial hydraulic fracturing stage, the total fracture volume is small, thus a reticulated secondary fracture structure is created around the injection node (i.e., the injection borehole) when subject to the impact of high-pressure fracturing liquid, as shown in Figure 9. Then the total volume of fractures increases gradually, which in turn increases the fracturing liquid leak-off and weakens fluid impact. During the later stage, the secondary fractures grow from the main fracture trunk like branches.
In Figure 9, the maximum width of the hydraulic fracture network gradually narrows as the stress difference increases, indicating that the general direction of the fracture network follows σ 1 under high stress difference conditions. However, when σ 1 − σ 3 = 2.5 MPa, the main hydraulic fractures propagate straight along the natural fractures. These results suggest that the general direction of hydraulic fracture networks in coal seams is jointly influenced by geo-stress differences and pre-existing natural fracture networks. When the ground stress difference is insignificant, the natural fracture network dominates the direction of the main hydraulic fractures. On the contrary, the maximum principal stress shows strong control over the general direction of the hydraulic fracture network, but the hydraulic fracture network still has many secondary fractures.
When the stress difference is 7.5 MPa or 10 MPa, the hydraulic fracture network is spindle-shaped. As the fractures propagate, the general direction of the hydraulic fracture network gradually converges to the direction of σ 1 , and the number of secondary fractures gradually decreases. To better understand how the stress field and natural fractures affect the hydraulic fracture morphology, tests that consider the direction and discontinuous degree of the natural fractures must be carried out in the future. The improved cohesive element method used in this study is leading to continuing research on these topics.
The total lengths of hydraulic fractures under different stress conditions were counted, as illustrated in Figure 11. According to the figure, the total lengths of the hydraulic fractures were between 82.4 m (when σ 1 − σ 3 = 2.5 MPa) and 100.7 m (when σ 1 − σ 3 = 5 MPa). There was no obvious correlation between the total length and stress difference. In each stress condition, the total length of the secondary fractures was significantly bigger than that of the main fracture. The ratio of secondary fracture total length to that of the main fracture was between 2.11 and 3.62, which shows no obvious correlation with stress difference. The results indicate that: (1) hydraulic fracturing can create a large number of secondary fractures in coal seams; (2) stress difference has no significant effect on the total length of the hydraulic fracture and the ratio of the secondary fracture length to the main fracture length. of fractures increases gradually, which in turn increases the fracturing liquid leak-off and weakens fluid impact. During the later stage, the secondary fractures grow from the main fracture trunk like branches.
In Figure 9, the maximum width of the hydraulic fracture network gradually narrows as the stress difference increases, indicating that the general direction of the fracture network follows  1 under high stress difference conditions. However, when    1 3 = 2.5 MPa, the main hydraulic fractures propagate straight along the natural fractures. These results suggest that the general direction of hydraulic fracture networks in coal seams is jointly influenced by geo-stress differences and pre-existing natural fracture networks. When the ground stress difference is insignificant, the natural fracture network dominates the direction of the main hydraulic fractures. On the contrary, the maximum principal stress shows strong control over the general direction of the hydraulic fracture network, but the hydraulic fracture network still has many secondary fractures.
When the stress difference is 7.5 MPa or 10 MPa, the hydraulic fracture network is spindleshaped. As the fractures propagate, the general direction of the hydraulic fracture network gradually converges to the direction of  1 , and the number of secondary fractures gradually decreases. To better understand how the stress field and natural fractures affect the hydraulic fracture morphology, tests that consider the direction and discontinuous degree of the natural fractures must be carried out in the future. The improved cohesive element method used in this study is leading to continuing research on these topics.
The total lengths of hydraulic fractures under different stress conditions were counted, as illustrated in Figure 11. According to the figure, the total lengths of the hydraulic fractures were between 82.4 m (when    1 3 = 2.5 MPa) and 100.7 m (when    1 3 = 5 MPa). There was no obvious correlation between the total length and stress difference. In each stress condition, the total length of the secondary fractures was significantly bigger than that of the main fracture. The ratio of secondary fracture total length to that of the main fracture was between 2.11 and 3.62, which shows no obvious correlation with stress difference. The results indicate that: (1) hydraulic fracturing can create a large number of secondary fractures in coal seams; (2) stress difference has no significant effect on the total length of the hydraulic fracture and the ratio of the secondary fracture length to the main fracture length.

Growth Process of the Secondary Hydraulic Fractures
According to Figures 9 and 10, the secondary fractures contain new fractures and natural fractures. There are two causes of secondary cracks: change to local stress fields caused by hydraulic fracture propagation and seepage in the porous coal matrix, and pre-existing low-intensity natural fractures. Figure 12 shows the growth process of secondary fractures during hydraulic fracturing when σ 1 − σ 3 = 7.5 MPa. As the main fracture propagates, the secondary fractures gradually branch out of the main fracture trunk. At t = 8.380 min, fracture branch 1 appears. At t = 13.244 min, fracture branches 2 to 6 appear, and branch 1 grows longer due to the fracturing liquid. At t = 21.696 min, new fractures 7 to 13 appear as the main fracture propagates. Note that fractures 7, 9 and 10 grow on the secondary fractures and are considered micro-secondary fractures. The hydraulic fracture network shows a multi-level branch structure.
The main fracture in Figure 12 propagates from natural fractures to the coal matrix, and then swerve back along the natural fracture. Affected by the stress field and discontinuous natural fracture network, some secondary fractures, such as fracture branches 5, 6, 8, 12, propagate along the discontinuous natural fractures, whereas the other secondary fractures start in the coal matrix. Although the tensile strength of the coal matrix is higher than that of the natural fracture, the difficulty of secondary fracture initiation in the coal matrix is almost equal to that in natural fractures due to the high geo-stress conditions. This result suggests that differences in geo-stress play an important role in inducing secondary fractures.

Growth Process of the Secondary Hydraulic Fractures
According to Figures 9 and 10, the secondary fractures contain new fractures and natural fractures. There are two causes of secondary cracks: change to local stress fields caused by hydraulic fracture propagation and seepage in the porous coal matrix, and pre-existing low-intensity natural fractures. Figure 12 shows the growth process of secondary fractures during hydraulic fracturing when    1 3 = 7.5 MPa. As the main fracture propagates, the secondary fractures gradually branch out of the main fracture trunk. At t = 8.380 min, fracture branch 1 appears. At t = 13.244 min, fracture branches 2 to 6 appear, and branch 1 grows longer due to the fracturing liquid. At t = 21.696 min, new fractures 7 to 13 appear as the main fracture propagates. Note that fractures 7, 9 and 10 grow on the secondary fractures and are considered micro-secondary fractures. The hydraulic fracture network shows a multi-level branch structure.
The main fracture in Figure 12 propagates from natural fractures to the coal matrix, and then swerve back along the natural fracture. Affected by the stress field and discontinuous natural fracture network, some secondary fractures, such as fracture branches 5, 6, 8, 12, propagate along the discontinuous natural fractures, whereas the other secondary fractures start in the coal matrix. Although the tensile strength of the coal matrix is higher than that of the natural fracture, the difficulty of secondary fracture initiation in the coal matrix is almost equal to that in natural fractures due to the high geo-stress conditions. This result suggests that differences in geo-stress play an important role in inducing secondary fractures. Figure 12. The numerical simulation results for the secondary hydraulic fracture growth process and pore pressure distribution. The serial numbers represent the order of secondary fracture initiation. Figure 12 also shows pore pressure distribution in the coal matrix and in hydraulic fractures. Note that there are two low pore pressure zones (the blue zones) in front of the two main fracture tips and two high pore pressure zones (the yellow zones) at both sides of the main fractures. Figure  13 shows pore pressure variation in the coal matrix on both sides of the hydraulic fractures. In Figure   Figure 12. The numerical simulation results for the secondary hydraulic fracture growth process and pore pressure distribution. The serial numbers represent the order of secondary fracture initiation. Figure 12 also shows pore pressure distribution in the coal matrix and in hydraulic fractures. Note that there are two low pore pressure zones (the blue zones) in front of the two main fracture tips and two high pore pressure zones (the yellow zones) at both sides of the main fractures. Figure 13 shows pore pressure variation in the coal matrix on both sides of the hydraulic fractures. In Figure 13, the distance between the pore pressure peak and main fractures is approximately 1.7 m. The spindle-shaped zone, where the hydraulic fracture network is located, can be regarded as a coal damage zone that changes the transfer path of σ 1 from the far field to near field, as shown in Figure 14. The compressive stress transfer path bypasses the damage zone, resulting in a compressive stress concentration on both sides of the damage zone, which increases pore pressure. In front of the main hydraulic fractures, the coal matrix is subjected to tensile stress by hydraulic fracturing, therefore its volume tends to expand, lowering the pore pressure.

Pore Pressure Distribution Characteristics
Processes 2018, 6, x FOR PEER REVIEW 18 of 25 13, the distance between the pore pressure peak and main fractures is approximately 1.7 m. The spindle-shaped zone, where the hydraulic fracture network is located, can be regarded as a coal damage zone that changes the transfer path of σ1 from the far field to near field, as shown in Figure  14. The compressive stress transfer path bypasses the damage zone, resulting in a compressive stress concentration on both sides of the damage zone, which increases pore pressure. In front of the main hydraulic fractures, the coal matrix is subjected to tensile stress by hydraulic fracturing, therefore its volume tends to expand, lowering the pore pressure.

Effect of the Injection Rate on Bottomhole Pressure
Due to the large scale of the model (400 m 2 ), it is difficult to consider the injection borehole profile, which is usually less than 0.04 m 2 in the field; thus an injection node was used instead. The bottomhole pressure variation during hydraulic fracturing was captured by the injection node. Figure  15 illustrates the bottomhole pressure variations at different injection rates under the stress condition of  1 = 17.5 MPa and  3 = 10 MPa. According to Figure 15, within the first two minutes of hydraulic fracturing, the bottomhole pressure increased rapidly due to the injection fluid impact    13, the distance between the pore pressure peak and main fractures is approximately 1.7 m. The spindle-shaped zone, where the hydraulic fracture network is located, can be regarded as a coal damage zone that changes the transfer path of σ1 from the far field to near field, as shown in Figure  14. The compressive stress transfer path bypasses the damage zone, resulting in a compressive stress concentration on both sides of the damage zone, which increases pore pressure. In front of the main hydraulic fractures, the coal matrix is subjected to tensile stress by hydraulic fracturing, therefore its volume tends to expand, lowering the pore pressure.

Effect of the Injection Rate on Bottomhole Pressure
Due to the large scale of the model (400 m 2 ), it is difficult to consider the injection borehole profile, which is usually less than 0.04 m 2 in the field; thus an injection node was used instead. The bottomhole pressure variation during hydraulic fracturing was captured by the injection node. Figure  15 illustrates the bottomhole pressure variations at different injection rates under the stress condition of  1 = 17.5 MPa and  3 = 10 MPa. According to Figure 15, within the first two minutes of hydraulic fracturing, the bottomhole pressure increased rapidly due to the injection fluid impact

Effect of the Injection Rate on Bottomhole Pressure
Due to the large scale of the model (400 m 2 ), it is difficult to consider the injection borehole profile, which is usually less than 0.04 m 2 in the field; thus an injection node was used instead.
The bottomhole pressure variation during hydraulic fracturing was captured by the injection node. Figure 15 illustrates the bottomhole pressure variations at different injection rates under the stress condition of σ 1 = 17.5 MPa and σ 3 = 10 MPa. According to Figure 15, within the first two minutes of hydraulic fracturing, the bottomhole pressure increased rapidly due to the injection fluid impact effect, and then decreased rapidly due to hydraulic fracture initiation. The peak pressure is called the break-down pressure. After three minutes, the bottomhole pressure stabilized. The stable bottomhole pressure and the break-down pressure increased as the injection rate increased, indicating that hydraulic fracture propagation is a dynamic process. Under high injection rates, the fracture volume rate increase is slower than the injection liquid volume rate increase. This explains why the bottomhole pressure increased as the injection rate increased. effect, and then decreased rapidly due to hydraulic fracture initiation. The peak pressure is called the break-down pressure. After three minutes, the bottomhole pressure stabilized. The stable bottomhole pressure and the break-down pressure increased as the injection rate increased, indicating that hydraulic fracture propagation is a dynamic process. Under high injection rates, the fracture volume rate increase is slower than the injection liquid volume rate increase. This explains why the bottomhole pressure increased as the injection rate increased.

Break Pressure of the Discontinuous Natural Fracture
As shown in Figure 16, when the hydraulic fracture propagates along discontinuous natural fractures, the break pressure of the fracture initiation p f should overcome the sum of the normal compressive stress on discontinuous natural fractures p n and tensile strength T f . That is: Based on the stress balance principle, p n can be calculated as: where θ is the angle between the normal direction of discontinuous natural fractures and the σ 1 direction. σ 1 and σ 3 are the effective stresses corresponding to σ 1 and σ 3 , respectively. The effective stresses can be calculated using Equation (1). By substituting Equation (27) into Equation (26), p f can be calculated as: Equation (28) clearly states that the break pressure, p f , depends on the ground stress conditions, tensile strength and direction of discontinuous natural fractures.
When σ 1 = 10 MPa, σ 3 = 17.5 MPa, u w (pore pressure) = 5 MPa, and θ = 45 • , p f was approximately 8.755 MPa for the natural fractures, and 9.52 MPa for other parts of the discontinuous natural fractures that have the same tensile strength as the coal matrix. Compared to the ground stress, the coal seam tensile strength was very low. According to Equation (28), in deep coal seams, σ 1 and σ 3 are much greater than T f , so the break pressure of discontinuous natural fractures is mainly determined by the ground stress conditions and the direction of the discontinuous natural fractures.
Hard and tight rock layers, such as geothermal energy reservoirs, have an intact rock tensile strength significantly higher than that of the joints, so hydraulic fractures tend to propagate along natural fractures. When the fracturing liquid is full of natural fractures, the fracture shear strength is drastically reduced and eventually fails. This is why the PSS model is widely used in EGS to predict fracture propagation, but is not suitable for the prediction of hydraulic fractures in deep coal seams.
Note that the above analytical process is based on the assumption of static fracturing, which is different from true dynamic fracturing. Comparing the above analytical results with the numerical simulation results in Figure 15, it can be seen that the simulated bottomhole pressure was higher than the static analytical break pressure result. Hard and tight rock layers, such as geothermal energy reservoirs, have an intact rock tensile strength significantly higher than that of the joints, so hydraulic fractures tend to propagate along natural fractures. When the fracturing liquid is full of natural fractures, the fracture shear strength is drastically reduced and eventually fails. This is why the PSS model is widely used in EGS to predict fracture propagation, but is not suitable for the prediction of hydraulic fractures in deep coal seams.
Note that the above analytical process is based on the assumption of static fracturing, which is different from true dynamic fracturing. Comparing the above analytical results with the numerical simulation results in Figure 15, it can be seen that the simulated bottomhole pressure was higher than the static analytical break pressure result.

Conclusions
In this study, an innovative finite element meshing scheme using cohesive elements (including global cohesive element embedding and global pore pressure node sharing techniques) was used to model hydraulic fracturing. Fracture propagation and interaction with discontinuous natural fracture networks in coal seams were simulated using the cohesive element method. The hydraulic fracture network characteristics under different stress conditions, the growth process of secondary hydraulic fractures, the pore pressure distribution in coal seams, and the variation of fracturing liquid pressure were investigated in detail. The numerical simulation results show that the improved cohesive element method in this study is capable of modeling hydraulic fracturing in naturally fractured formations.
In coal seams, hydraulic fracturing generally creates a complex hydraulic fracture network that is composed of a main fracture trunk and many secondary fracture branches. The morphology of the fracture network is determined by the geo-stress field and the discontinuous natural fracture network. Under high stress difference conditions, the hydraulic fracture network is spindle-shaped, and shows a multi-level branch structure. The higher the geo-stress difference, the narrower the width of the spindle-shaped hydraulic fracture network.
The ratio of secondary fracture total length to main fracture total length was 2.11~3.62, indicating that the secondary fractures are an important part of the hydraulic fracture network. The geo-stress difference had no significant effect on the ratio of secondary fracture length to main fracture length.
The stable bottomhole pressure and the break-down pressure increased with the injection rate, indicating that fracture propagation is a dynamic process in coal seams. In deep coal seams, the break pressure of discontinuous natural fractures mainly depends on the in-situ stress field and the direction of natural fractures. The mechanism of hydraulic fracture propagation in deep coal seams is significantly different from that in hard and tight rock layers. The PSS model, which is more frequently used in EGS to predict hydraulic fracture propagation, is not suitable for deep coal seams.  Figure 16. Schematic of the analytical model for calculating the break pressure of discontinuous natural fractures. p f is the break pressure; p n and p s are the normal compressive stress and the tangential shear stress of discontinuous natural fractures.

Conclusions
In this study, an innovative finite element meshing scheme using cohesive elements (including global cohesive element embedding and global pore pressure node sharing techniques) was used to model hydraulic fracturing. Fracture propagation and interaction with discontinuous natural fracture networks in coal seams were simulated using the cohesive element method. The hydraulic fracture network characteristics under different stress conditions, the growth process of secondary hydraulic fractures, the pore pressure distribution in coal seams, and the variation of fracturing liquid pressure were investigated in detail. The numerical simulation results show that the improved cohesive element method in this study is capable of modeling hydraulic fracturing in naturally fractured formations.
In coal seams, hydraulic fracturing generally creates a complex hydraulic fracture network that is composed of a main fracture trunk and many secondary fracture branches. The morphology of the fracture network is determined by the geo-stress field and the discontinuous natural fracture network. Under high stress difference conditions, the hydraulic fracture network is spindle-shaped, and shows a multi-level branch structure. The higher the geo-stress difference, the narrower the width of the spindle-shaped hydraulic fracture network.
The ratio of secondary fracture total length to main fracture total length was 2.11~3.62, indicating that the secondary fractures are an important part of the hydraulic fracture network. The geo-stress difference had no significant effect on the ratio of secondary fracture length to main fracture length.
The stable bottomhole pressure and the break-down pressure increased with the injection rate, indicating that fracture propagation is a dynamic process in coal seams. In deep coal seams, the break pressure of discontinuous natural fractures mainly depends on the in-situ stress field and the