Numerical Simulation of Fracturing in Coals Using Water and Supercritical Carbon Dioxide with Potential-Based Cohesive Zone Models

The Park-Paulino-Roesler (PPR) cohesive zone model (CZM) for coal was established for analyzing mixed-mode I/II fractures using semicircular specimens under punch-through shear (PTS) and three-point bending (SCB) tests. In these methods, the main parameters of the fracture were obtained through SCB tests and PTS tests. And according to the experimental results, the coal specimens show obvious characteristics of ductile fracture under mode I and II loading. Moreover, hydraulic and supercritical carbon dioxide (ScCO2) fracture tests were conducted, and accordingly, it was found that the crack initiation pressure of coal specimens for hydraulic fracturing is 17.76MPa, about 1.59 times that driven by ScCO2. And the crack initiation time of coal with ScCO2 fracturing is 123.73 s, which is 1.58 times that for hydraulic fracturing. A macrocrack eventually formed in the coal specimen due to the hydraulic drive, which penetrated through the entire specimen. Yet, there was no crack penetrating the whole fracture specimen and several widely distributed secondary cracks in the fractured coal specimens by ScCO2. Furthermore, zero-thickness pore pressure cohesive elements were utilized to investigate multicrack propagation in coals undergoing hydraulic and ScCO2 fracturing. The constitutive relationships of the established PPR CZM were introduced into the cohesive elements. The obtained results are consistent with the hydraulic and ScCO2 fracturing experiment results for the coal specimens. This indicates that the established PPR CZMs can accurately represent the crack propagation behavior in coals for hydraulic and ScCO2 fracturing.


Introduction
As an essential type of clean energy, the exploitation of coalbed methane (CBM) is of significant importance to increase the supply of clean energy, thereby decreasing concerns about greenhouse gases and realizing safe coal mining [1,2]. Studies show that the low permeability of coal seams is one of the main challenges for the efficient exploitation of CBM. Generally, coal permeability in Chinese mines is less than 1 mD [3], which is much lower than that in the United States, Australia, and other countries. In this regard, hydraulic fracturing is a widely adopted technology to improve the CBM permeability and, therefore, production by injecting a large volume of water-based fluid to create and extend fracture networks [4,5]. Hence, in the process of hydraulic fracturing, the crack propagation behavior in coals will have a direct influence on the effect of CBM exploitation. However, there are some drawbacks to the hydraulic fracturing technology; for example, it will cause a lot of waste and pollution of water resources, and the fracturing fluid will cause "water sensitive" and "water lock" influence on coalbed methane reservoir [6]. In order to resolve these shortcomings, numerous nonaqueous fracturing technologies have been proposed [7,8], among which the supercritical carbon dioxide (ScCO 2 ) fracturing technology has attracted much attention [9]. ScCO 2 refers to a special state of CO 2 when its temperature and pressure exceed 31.1°C and 7.38 MPa, respectively, which has unique physical and chemical characteristics, including low viscosity, high diffusion coefficient, and high density [10]. Some researchers have conducted experiments on rocks with ScCO 2 fracturing, and the results show that ScCO 2 fracturing can produce more widely distributed and complex fracture networks in rocks than hydraulic fracturing, which can significantly increase reservoir permeability [11][12][13]. In addition, ScCO 2 also has a good displacement effect on methane adsorbed in the coal seam, which will be beneficial to improve the yield of coalbed methane [14]. Therefore, ScCO 2 fracturing technology can promote the efficient exploitation of coalbed methane, and the most important thing is the research on the crack propagation law of coal driven by ScCO 2 .
Since the groundbreaking work of Irwin [15] and Griffith [16], linear elastic fracture mechanics (LEFM) was established, becoming a highly effective theory framework for analyzing crack propagation in brittle solid materials. Bieniawski [17,18] systematically introduced the LEFM into the research of crack propagation behavior in rocks, and since then, rock fracture mechanics has been widely used in rock materials. Generally, the fracture toughness (K c ) is applied as an indicator, reflecting the crack propagation in natural materials [19][20][21]. Nevertheless, LEFM is mostly limited to investigating crack propagation in brittle rocks. Yet, some soft rocks, such as coal, exhibit generally ductile failure behaviors, represented by an obvious strainsoftening stage after the peak stress when the crack initiates [22,23]. It is because of the fracture process zone (FPZ) [24] of these soft rocks, i.e., the particular region in front of the crack tip, where a series of nonlinear softening behaviors including microcrack initiation, and plastic strain and mineral crystal friction occurrence, are sizable. It is also nonnegligible relative to the size of the rock specimen and the size of the crack. In comparison to brittle rocks, abundant primary pores and microfissures exist in the coal body [25], causing the ductile fracture characteristics of coals to be more prominent. Thus, the theory of LEFM does not apply to the study of the fracture behavior of coals.
The CZM inspired by the studies of Barenblatt [26], Dugdale [27], and Hillerborg et al. [28] has been used with success to represent the crack propagation behavior in nonlinear FPZ of ductile materials. In this theory, the FPZ is simplified hypothetically to a discrete line or plane corresponding to either a two-dimensional or three-dimensional case, respectively, in which the hypothetical cohesive stress causes the virtual crack to close (see Figure 1). The constitutive relation of CZM is represented by the relationship between the cohesive stress and relative displacement across this line or plane. The above constitutive relation is usually nonlinear and depends on the form of stress and the evolution characteristics of damage variables. When the material in this region is completely damaged, the cohesive stress will be lost, which means that a new macrocrack surface is generated. The energy consumed in this damaging process is the fracture energy of the material. Accordingly, it is concluded that the cohesive crack model can effectively characterize the ductile fracture behavior of coal.
Based on CZM, the numerical crack propagation model for hydraulic fracturing in a rock is established by the finite element method [29], extended finite element method [30], etc. However, the constitutive relationship of the softening curve has a huge impact on fracture behaviors [31,32], and the linear or bilinear constitutive relationships of CZMs have been adopted in previous numerical models. Hence, it is necessary to establish the CZM of coals to provide an accurate numerical model to predict crack propagation. Mixed-mode I/II crack propagation is prone to occurrence in coals under engineering conditions, especially in supercritical carbon dioxide fracturing [33]. Reviewing the literature indicates that since the LEFM method does not reflect ductile fractures, it is the most widely used scheme to investigate the crack propagation in coal [34,35]. On the other hand, cohesive interactions between fractured surfaces are the main failure mechanisms in the mixed-mode I/II CZM. It should be indicated that these interactions can be expressed through stress-strain equations in fractured surfaces. Nonpotential-based models [36][37][38] were established to characterize the ductile fracture behavior of materials. Considering symmetric systems in cohesive interactions, these models can be simply developed.
Nevertheless, the main drawback of the nonpotentialbased model is that one model cannot explain all possible separations in ductile materials. Furthermore, asymmetric tangential stiffness of material increases the computational expenses. An effective solution for this problem is to apply potential-based models to utilize the initial derivative of the fracture potential energy function [39]. This scheme is based on the cohesive stress over the fractured surfaces, while the second derivative reflects the constitutive association. Based on potential-based models, Park et al. [40] proposed the Park-Paulino-Roesler (PPR) model to simulate the cohesive fracture [41,42]. In this model, fracture energy (including modes I and II) and different initial slopes and cohesive strengths are considered. Meanwhile, corrective variables are defined to cover a wide range of failures in different ductile materials. This model resolves the disadvantages of traditional potential-based models.
In this work, we performed semicircular specimens under PTS and SCB tests to calculate the fracture parameters    Figure 2, coal samples were prepared into semidisk-like specimens, and an artificial crack was prefabricated along the symmetric center starting from the center of the bottom edge of the specimen, and vertical loads were applied on the top of the arc to cause mode I fracture of the specimen. The diameter (2R) and thickness (B) of the coal rock SCB specimen were set as 70 mm and 25 mm, respectively. The ratio (α/R) of the preset crack length to the specimen radius was set as 0.35, and the ratio (S/R) of the base supporting roller span to the specimen diameter was 0.5. The loading mode is displacement control, and the loading speed is 0.02 mm/min. In addition, the crack tip opening displacement (CTOD) of the specimen was measured by the fiber grating (FBG) technique (with an accuracy of 0.5 microstrains) throughout the experiment. In this study, three groups of effective SCB tests were conducted on the coal specimens. In this regard, Figure 3 illustrates the load-CTOD curve of coal SCB specimens during the whole experiment, and according to the various characteristics of the experimental curve, the experimental process is generally divided into four steps, including the compaction, elastic deformation, peak load stage, and postpeak damage stages. When crack initiation occurs in coal specimens, the accumulated energy is not released instantaneously, and there is a nonlinear damage process in the postpeak loading stage. The CTOD of the three coal SCB specimens increased by 0.1192 mm, 0.1153 mm, and 0.0895 mm, respectively, with an average value of 0.108 mm, from the beginning of the specimen subjected to the force to the formation of a new crack surface, that is, from the intact specimen to the fracture of the specimen.

PTS Test.
This method was first proposed by Backers et al. [43] to investigate the fracture of materials in mode II loading conditions. The PTS test was used here to simulate mode II fracture experiments in coal, as it is easy to process, and the experimental results are reliable. As shown in Figure 4, the specimen was a circular cylinder with a diameter D and had circular notches with a diameter ID drilled into the upper and lower end faces along the central axis of the cylinder. Circular cylinders with diameter D and height L were prepared. Moreover, two circular notches were prepared with depths a and b near the upper end and lower end of the specimen, respectively. The width of notches was set to t, and the effective shear length was IP = L-a-b. The experiment was performed using a loading cylinder to apply a vertical shear force (P). It is found that as the applied shear stress increases, the crack propagates along the notches parallel to the axis of the cylinder, as well as mode II fracture characteristics which were then able to be acquired from the experimental results. The diameter D and height L of the specimens was set to 50 mm, and the coals were cut into the PTS specimens using diamond wire cutting under the CNC machine tools. This can limit the micromechanical damage to the coal specimens and improves machining accuracy. In addition, notches with a diameter of ID = 35 mm were prefabricated along the same central axis utilizing a 0.5-millimeter-thick diamond bit with the CNC machine tools. Parameters a and b of these notches were set to 10 mm and 30 mm, respectively. Moreover, the Figure 2: Coal SCB specimen [20]. 3 Geofluids length of IP was set to 10 mm. Over the experiment, the specimen of coal PTS was located in advance on the bottom support, which possesses cylindrical grooves with a diameter of 35 mm and a depth of 15 mm. A shear load was employed to the specimen of coal using a loading cylinder with a diameter of 35 mm. Finally, the test was in the mode of displacement control with a constant rate of 0.02 mm/min to ensure stable crack propagation. Three experiments were executed each for the coal. The stressstrain curves were recorded. Figure 5 shows the experimental curve of coal shear load and tangential displacement, which represents the typical coal type II fracture characteristics. In the initial stage, the shear load has a linear correlation with the shear displacement. When the critical value was obtained for the shear displacement, mode II cracks begin to occur in the FPZ of the coal sample and local damage occurs. Within the postpeak stage, the shear load progressively reduces with the enhancement of shear displacement, and the coal sample presents the characteristics of ductile fracture. The average maximum tangential displacement of the PTS coal specimen is 0.055 mm, and the nonlinear damage softening stage appears in the postpeak stage of the shear process of the PTS coal specimen. In addition, by calculating the ratio of 4 Geofluids the peak shear load to the effective shear area, the shear strength (τ t ) of coal can be obtained directly. The calculation formula is as follows: where P scr is the peak shear load. According to the above formula, the shear strength of the coal sample is 2.36 MPa.

PPR CZM for Coals
In the PPR model, the normal and tangential cohesive interactions (T n , T t ) are functions of the normal or tangential separation (Δ n , Δ t ), respectively. It should be indicated that T n approaches zero when Δ t reaches the tangential conjugate final crack opening displacement ( δ t ) or Δ n reaches the maximum normal crack opening width (δ n ). This represents complete normal failure. Similarly, when Δ t reaches its maximum displacement of tangential crack opening (δ t ) or Δ n attains the normal conjugate final crack opening width ( δ n ), full tangential failure takes place. The expressions are as follows: When Δ n ðΔ t Þ reaches the critical opening width δ nc ðδ tc Þ, the value of T n ðT t Þ is the maximum normal cohesive strength (σ max ). This is shown as follows: The mode I and mode II fracture energy (Φ n , Φ t ) can be calculated by the area underneath the cohesive interactions, as follows: In this study, the mode I and II fracture energies of coal samples were calculated, respectively, by the unit integral area under the load-relative opening displacement curves in Figures 3 and 5 of the above two kinds of tests. The specific shape of the softening response, i.e., the constitutive relationship of the softening process, remarkably affects the crack propagation. Therefore, nondimensional shape parameter indices (α, β) are introduced into the PPR model. When the shape parameter indices are equal to 2, the gradient of the softening process represents a nearly linear relationship. If the shape criteria are lower than 2, the cohesive interactions have concave softening trends. Conversely, if the shape indi-ces are higher than 2, the gradient of the softening procedure demonstrates a convex shape.
Considering the foregoing macroscopic fracture criteria and boundary conditions, the potential energy function can then be mathematically expressed in the form below: The cohesive interactions T n and T t are obtained by taking the first derivative of the PPR model along with the normal vector and tangential vector, respectively, as follows: where <· > is the Macaulay bracket function, whose calculation is as follows: where m and n are the nondimensional exponents, which are determined by the shape parameter indices (α, β) and the boundary conditions of the critical separations (Equations (4) and (5)). m and n are determined by where λ n and λ t are the initial slope indicators, i.e., the ratio of the critical crack opening displacement to the maximum crack separation displacement, as determined by where Γ n and Γ t are considered energy constants, which are functions of mode I and II fracture energy (Φ n , Φ t ). When Φ n is different from Φ t , the formulas of the energy constants are as follows: When the values between Φ n and Φ t are equal, the simplification of energy constants to the following expression is possible: The unknown parameters of the PPR model needed to establish the PPR CZMs for the coals can be calculated based upon data acquired from the SCB tests and PTS tests. The maximum normal and tangential crack opening width can be determined by considering the boundary conditions of the cohesive strength (4) and (5) and fracture energy (6). The equations are as follows: where the parameters of δ n , δ t , Φ n , Φ t , σ max , and τ max have already been determined. The initial slope indicators λ n and λ t can be calculated by Equation (12). The calculation results of the different coals are listed in Table 1. Finally, substituting Equation (11) into Equations (15) and (16), the nondimensional shape parameter indices (α, β) could be ascertained by solving the above equations. The values of α ðβÞ for the three different coals are shown in Table 1 as well.
Based on the experimental results, Φ n is different from Φ t for the coals; hence, the energy constants Γ n and Γ t were calculated using Equation (12). The results are listed in Table 1. In addition, in order to determine the cohesive interaction region, the last displacements of conjugate crack opening δ t and δ n can be calculated by Equations (14) and (15). The cohesive interaction region of mode I is (0.108, 0.009), and the cohesive interaction region of mode II is (0.108, 0.055). When the normal or tangential separation displacement (Δ n , Δ t ) exceeded the region, the cohesive stresses of mode I and mode II were set to zero.

Crack Propagation Experiment of Coals for
Hydraulic and ScCO 2 Fracturing 4.1. Experimental Preparation and Process. In this section, trimmed samples with a diameter of 50 mm and a length of 100 mm were used in the experiment. Figure 6 illustrates the configuration of samples, indicating that there is a central borehole in the upper surface of the specimen. Then, a 3 mm steel pipe was inserted into the borehole to inject the fluid to simulate the fracturing well. The schematic diagram of the hydraulic and ScCO 2 fracturing experimental device is shown in Figure 7. In order to perform the fracturing tests,   6 Geofluids the prepared specimens were placed in a pressurized kettle. In order to prevent damage to the specimens, the pressure increment rate was set to 1 MPa min −1 . Meanwhile, σ a and σ c were set to 10 MPa and 8 MPa, respectively. In addition, the temperature of the triaxial pressure kettle was set at 40°C and maintained for 3 hours prior to the fracturing test to ensure that the coal sample was fully heated, and the injection flow of the fracturing fluid was set to 20 mL min −1 . When ScCO 2 fracturing is performed, the steel injection pipe was heated to 40°C in advance so that the temperature of carbon dioxide is above its critical temperature (31.1°C). Sample breakdown was identified to have occurred once the fracturing fluid pressure reduces suddenly and simultaneously; the confining pressure increases abruptly.

Hydraulic and ScCO 2 Fracturing Experimental Results.
Two samples were run for each experimental condition. Figure 8 shows the four fluid pressure-time curves of coal specimens by the two kinds of fluid fracturing. The fracturing process of all coal specimens can be divided into three stages. The first stage consists of fluid pressure rising. Fracturing fluid is continuously injected into the fracture specimen by a high-pressure pump to resist the strength of the specimen under confining pressure. At the beginning of this stage, the growth rate of the fluid pressure in each coal specimen is very low, especially for ScCO 2 fracturing, because it takes time to fill the anhydrite section of the coal specimen after the fracturing fluid injection. In addition, fracturing fluid injection into the coal body immerses and infiltrates the specimen. In this early stage, the slow increase in fluid pressure becomes significant for ScCO 2 fracturing; this is attributed to the relatively developed fracture structure in the center hole of the coal body, and the infiltration effect of ScCO 2 is noticeable than water. The second stage consists of crack initiations. When the injected fluid pressure reaches a specific critical value, the critical condition for the crack propagation of the coal specimen is reached, and the coal specimen ruptures. The critical pressure is called the crack initiation pressure of fracturing, and the critical fracturing time it takes to reach the crack initiation pressure is called the crack initiation time. The two kinds of fluid fracturing results of the coal specimens are shown in Table 2. The average crack initiation pressure of coal specimens for hydraulic fracturing is 17.76 MPa, about 1.59 times that driven by ScCO 2 fracturing. The third stage is the pressure drop stage; when fracturing occurs in the coal sample, the water pressure accordingly decreases. For hydraulic fracturing in coals, the fluid pressure obviously decreases after fracturing. And there is a significant fluctuation in the fluid pressure after ScCO 2 fracturing. This is because the crack in the coal sample does not completely penetrate the sample. Due to the surrounding rock and axial pressure, the crack closes once again in the coals, and the continuously injected ScCO 2 will drive crack propagation in the specimen repeatedly until the specimen is completely broken. Figure 9 shows the final crack propagation morphology in the cylindrical coal specimens. The coal specimen eventually formed a macrocrack under the hydraulic drive, which penetrated through the whole cylindrical specimen. Yet, there was no crack penetrating the whole fracture specimen and several widely distributed secondary cracks in the fractured coal specimens by ScCO 2 . This is because water has greater viscosity and density, which is easy to produce tensile failure in coal specimens, and eventually forms a single penetrating crack. On the other hand, because of the large diffusion coefficient and strong permeability of ScCO 2 , the influence range in the coal is large, so it is easy to form a

Numerical Simulation of Fracturing in Coals
Based on the PPR Model 5.1. Governing Equations. Hydraulic and ScCO 2 fracturing of coal is a complex, multifield coupling process. Compared with the multiphysical coupling problems in rock mechanic engineering, this crack propagation process in a solid is coupled with fracturing, causing it to be more challenging to model and calculate. The physical model of the fracturing process in coals is shown in Figure 10. Ω represents the entire range of the hydraulic fracturing models, and the fracture width w is located in the center of the model. The two sides of the fracture consist of coal materials. The coal is a porous medium, containing the solid skeleton and pores.
In the fracture, Q represents the flow rate of the injected fracturing fluid. Some of the fluid will be filtered along the upper and lower crack surfaces and permeate into the coal through the pores. The pressure generated by the injected fracturing fluid in the fracture is defined as P f . When the fluid pressure reaches a critical value, the crack in the coal body will expand. The coal fracturing numerical simulation was performed through the pore pressure cohesive element. Figure 11 shows the pore pressure cohesive element with the fluid pressure node. When the cohesive element is affected by the external force, the upper node (1, 2) and the lower node (3,4) in the element are relatively displaced, which damages the cohesive element. Once the critical condition is reached, the cohesive element is destroyed, fracturing the material. During this process, the normal cohesion of the element also changes with the change in normal opening displacement. Also, the tangential cohesion of the element changes with the change of the tangential displacement. The relationship between the two cohesive forces and the displacement of the element nodes is the constitutive relationship of the cohesive crack. As shown in Figure 11, the pore pressure cohesive element consists of adding a group of pore pressure injection nodes (5,6) in the center of the original cohesive element. The injected fluid pressure is already included in the cell calculation model through this pressure node, and the injected fluid pressure causes the relative displacement on the lower and upper surfaces of the cohesive element, resulting in continuous damage to the cohesive element until it is destroyed. This represents crack growth in the fracturing model. In this multifield coupled numerical model, there are four governing equations, including the deformation equation of porous media in coal, the pore seepage equation of porous media in coal, the fracture flow equation in coal, and the constitutive equation of a cohesive crack in coal.
(1) In this model, if the pores in the coal body are filled with a single liquid (water or ScCO 2 ), the deformation of the coal body includes deformation of the    Figure 11: The pore pressure cohesive element. 8 Geofluids solid skeleton and deformation of the liquid in the pores. According to the momentum conservation equation, the following formula can be obtained: where σ is the total stress tensor of the porous media, b is the physical force of the porous medium, and ρ is the density of the coal. The coal density has the following formula: where n is the porosity of the coal, ρ s is the solid skeleton density of coal, and ρ f is the fluid density in the coal pore. According to the Biot porous elasticity theory [44] and the Terzaghi theory [45], the effective stress in coal can be calculated as follows: where δ ij is the Kronecker-delta symbol, p w is the pore pressure in the coal, and α is the Biot coefficient. The Biot coefficient is defined as follows: where K b is the total volume modulus of the porous coal media and K S is the modulus of the solid skeleton in the coal. For the incompressible solid material, K s = ∞, α = 1.
If the Biot coefficient (α) is equal to 0, the porous media material will degenerate into a dense linear elastic solid material.
Based on the hypothesis of small deformation, the formula is as follows: where ε i,j is the strain tensor, u i is the displacement vector, and u i,j and u j,i are the partial derivatives of displacement. The constitutive relation between the stress and strain of coal can be expressed as follows: where D ijkl is the elastic tensor of coal.
(2) Pore seepage in the coal body should satisfy the following mass conservation equation: where 1/Q is the compressibility coefficient of the fluid, p w is the pore pressure of coal, and _ w w is the velocity vector of Darcy flow. The compressibility coefficient 1/Q can be calculated as follows: where n is the porosity of coal and K w is the modulus of the fluid. In Equation (16), the flow rate and fluid pressure gradient in porous media satisfy Darcy's law, and the expression is as follows: where ρ w is the density of the fracturing fluid and k w is the permeability coefficient of the fracturing fluid, which can be calculated as where μ w is the viscosity of the fracturing fluid and k is the permeation matrix. The effect of the inertia term of the fracturing fluid and the roughness of the crack surface are not considered. The fracturing fluid in the fracture itself can be divided into the tangential flow and normal flow. According to the mass conservation theorem, the fluid in the fracture should satisfy the following equation: where 1/W f is the compressibility of the fluid in the fracture, p f is the fluid pressure of the fracturing fluid in the fracture, w is the crack opening, s is the coordinate of the tangential direction along the fracture surface, q t is the flow rate of the fracturing fluid filtered from the upper surface of the fracture into the porous media, q b is the flow rate of the fracturing fluid infiltrating from the lower surface of the fracture into the porous media, QðtÞ is the flow rate of the source term of the fluid, and δðx, yÞ is the Dirac-delta function. If the fracturing fluid is incompressible, the first fluid compression term in the above equation can be ignored. The tangential flow and pressure gradient of the fracture fluid in the fracture satisfy the cubic seepage model [46,47]. They can be related by the following expression: where μ f is the viscosity of the fracture fluid. Some of the fracturing fluid in the fracture infiltrates into the coal body through the fractures. In this numerical model, the fluid flow rate is related to the gradient between the pore pressure in 9 Geofluids the coal and the pressure of the fracturing fluid in the fracture, calculated as follows: where p t and p b are the pore pressures on the upper and lower surfaces of the crack, respectively, and k t and k b are the fluid filtration coefficients on the upper and lower surfaces of the crack, respectively. The equation of the tangential fluid flow in the fracture and the equation of normal fluid flow are taken into the mass conservation equation and expressed as follows: (3)The constitutive relationship of the cohesive element in the numerical simulation of hydraulic fracturing is derived from the mixed-mode I/II PPR potential energy function. In Section 3, the fracture parameter values in the PPR CZMs of the coal were determined through the SCB tests and the PTS tests. The constitutive equations of mixed-mode I/II for a cohesive crack of the different coals were established. The normal and tangential cohesions were obtained by taking the first derivative of the normal displacement and tangential displacement, respectively, by the PPR potential function, and the stiffness matrix of the cohesive element is as follows: The stiffness components of the PPR cohesive constitutive equation are as follows:

Numerical Models of Hydraulic Fracturing for the Coals.
In order to compare the experimental fracturing results of the coals, the coal material parameters and boundary conditions of the coals in the fracturing numerical models were set to be the same as those in the experiments. Figure 12 shows the boundary conditions on the coal fracturing numerical geometric model with a size of 50 mm × 50 mm, as well as the meshing conditions for a section of the model. In the fracturing numerical model, the coal materials were characterized by triangular solid elements, and the pore pressure cohesive elements were inserted between the triangular solid elements to simulate multiple crack propagation driven by hydraulic or ScCO 2 fracturing. To avoid the influence of the overall stiffness of the model after a large number of cohesive elements were inserted between the solid elements (see Figure 12), the corresponding upper and lower nodes in the cohesive element and their intermediate fluid pressure nodes were defined at the same position in the local coordinate system. This numerical simulation method is called the zero-thickness element method [48]. And the fluid injection point was set at the center of the numerical model, and the fluid injection rate was set to 20 mL/min. The involved numerical simulation parameters are given in Table 3, and the PPR model parameters of the three types of coals are listed in Table 1 for the numerical simulation for the hydraulic and ScCO 2 fracturing of the coal specimens. In addition, the contrast numerical simulation of the fracturing in coals was also carried out in which the constitutive relationship of the pore pressure cohesive elements was represented by the common linear elastic fracture mechanics (LEFM). Figure 13 shows the numerical simulation results of hydraulic and ScCO 2 fracturing of the coals. The fracture initiation pressure of the coal with water or ScCO 2 is 18. 45       11 Geofluids in the coals are shown in Figure 14 based on the LFEM, and there are obvious deviations between the simulation results and the test results. Figure 15 shows the simulation results of crack growth for hydraulic and ScCO 2 fracturing in the coal samples. The crack induced by hydraulic fracturing expands along the direction of the maximum principal stress, and at the same time, secondary crack propagation occurs near the main crack. And the multiple crack propagation appears in the coal model with ScCO 2 fracturing. This is also consistent with the experimental results of crack propagation in coal specimens caused by hydraulic and ScCO 2 fracturing. Therefore, this demonstrates that the established PPR CZMs can accurately describe crack propagation behavior in varying coal types for hydraulic and ScCO 2 fracturing.

Conclusion
In this research, the mixed-mode I/II PPR cohesive zone model (CZM) of coals was determined using PTS and SCB tests. The constitutive relationships of the established PPR CZMs were introduced into the pore pressure cohesive elements to simulate crack growth in coals caused by hydraulic and ScCO 2 fracturing. In addition, hydraulic and ScCO 2 fracturing experiments on the coal specimens were performed, and the numerical simulation outcomes were compared with the corresponding experimental results. The following conclusions can be drawn: (1) Several key fracture parameters, including the maximum normal open displacement (δ n ), the maximum tangential open displacement (δ t ), mode I fracture energy (Φ n ), and mode II fracture energy (Φ t ), were obtained through SCB tests and PTS tests. According to the experimental results, there are visible nonlinear damage processes in the stage of postpeak loading, and the coal specimens show obvious characteristics of ductile fracture under mode I and II loading. In addition, the mode II fracture energy of coal type II is 51.62 J/m 2 , which is considerably greater compared with fracture energy of mode I for coal (22.16 J/m 2 ); this shows that the mode II crack propagation in coals will use remarkable energy in coals (2) In the hydraulic and ScCO 2 fracturing experiments of coals, the crack initiation pressure of coal specimens for hydraulic fracturing is 17.76 MPa, about 1.59 times that driven by ScCO 2 fracturing. And the crack initiation time of coal with ScCO 2 fracturing is 123.73 s, which is 1.58 times that for hydraulic fracturing. A macrocrack eventually formed in the coal specimen due to the hydraulic drive, which penetrated through the entire specimen, whereas there was no crack penetrating the whole fracture specimen and several widely distributed secondary cracks in the fractured coal specimens by ScCO 2 . This is because water has greater viscosity and density, which is easy to produce tensile failure in coal specimens, and eventually forms a single penetrating crack. On the other hand, because of the large diffusion coefficient and strong permeability of ScCO 2 , the influence range in the coal is large, so it is easy to form a wide range of tensile and shear mixed-mode cracks in the coal specimens (3) The PPR CZMs of the coal were established using PTS and SCB tests for analyzing the mixed-mode I/II crack propagation. Zero-thickness pore pressure cohesive elements were used to simulate multicrack propagation in coals caused by hydraulic and ScCO 2 fracturing. The constitutive relationships of the established PPR CZM were introduced into the cohesive elements. Overall, the numerical simulation results are consistent with the hydraulic and ScCO 2 fracturing experimental results for the coal specimens. This indicates that the established PPR CZMs can accurately represent crack propagation behavior in coals caused by hydraulic and ScCO 2 fracturing

Data Availability
The underlying data and figures can be found in the manuscript.

Conflicts of Interest
The authors declare that they have no conflicts of interest. 12 Geofluids