#### Abstract

High-energy gas fracturing (HEGF) and gas fracturing (GF) are considered to be efficient to enhance the permeability of unconventional gas reservoir. The existing models for HEGF mainly focus on the dynamic loading of stress wave or static loading of gas pressurization, rather than on the combined actions of them. Studies on the combination of HEGF and GF (HEGF+GF) are also few. In this paper, a damage-based stress wave propagation-static mechanical equilibrium-gas flow coupling model is established. Numerical model and determination of mesomechanical parameters in finite element analysis are described in detail. Numerical simulations on crack evolution under HEGF, GF, and HEGF+GF are carried out, and the impact of in situ stress conditions on crack evolution is discussed further. A total of 11 cracks with length of 2.3-4 m in HEGF, 4 main cracks with length of 6.5–8 m in GF, and 11 radial cracks with length of 2–11.5 m in HEGF+GF are produced. Many radial cracks around the borehole are formed in HEGF and extended further in GF. The crustal stress difference is disadvantageous for crack complexity. This study can provide a reference for the application of HEGF+GF in unconventional gas reservoirs.

#### 1. Introduction

As an efficient clean resource, unconventional gas, such as shale gas and coal bed methane, is of great significance to the sustainable development of economy and energy [1–5]. However, the permeability of unconventional gas reservoirs is usually extremely low, greatly prohibiting the efficient extraction of gas. Therefore, the economic viability of developing unconventional gas depends on the effective stimulation of the reservoirs [6–10]. Currently, hydraulic fracturing (HF) has been used widely to enhance production and has been proven successful in most oil, conventional gas and unconventional gas reservoirs. However, hydraulic fracturing has also caused many problems [6], such as the use of large volumes of water [11], polluting of groundwater [12], and clay swelling and water locking effect in formations with high clay content [13]. Alternatives to water for well fracturing, such as explosive fracturing (EF), high-energy gas fracturing (HEGF, high-pressure gas produced from burning a propellant), and gas fracturing (GF, high-pressure gas produced from compressing air) are considered to be potentially efficient to enhance shale gas and coal gas extraction.

The principal differences between EF, HEGF, and GF are peak pressure and duration time [14]. The peak pressure of GF is usually lower than 35 MPa with duration longer than 1 hour, which is similar to HF. The peak pressure of EF can reach 1000-7000 MPa with duration within 1 ms. And the peak pressure of HEGF ranges 70-350 MPa with duration of 5-30 ms. As a result, GF produce fractures usually in the form of two fractures perpendicular to the minimum principal stress orientation, and EF can produce several radial fractures but cause considerable crushing zone around the wellbore, while HEGF can produce multiple radial cracks without causing substantial damage to the rock around the wellbore [15]. Therefore, HEGF could potentially avoid the limitations in both GF and EF. Up to now, HEGF has been used in many unconventional gas reservoirs [16], such as tight sandstone reservoir [17], shale reservoir [18], coal bed methane reservoir [19], and carbonate reservoir [20].

It is generally agreed that two types of loading are applied on the borehole during HEGF, which is similar to explosive fracturing. First, a stress wave loading that travels outward from the borehole is followed by a longer duration gas pressurization loading [21]. The role of the stress wave is to create initial cracks, while the gas pressure leads to crack propagation. In fact, the stress wave can only initiate limited cracking and crushing of the rock near the borehole which would not exceed more than several borehole diameters [15].

A range of theoretical and computational models have also been developed to simulate the dynamic propagation of fractures in HEGF [15, 16]. The most widely used model was the linear elastic fracture mechanics-based models proposed by Nilson et al. [22] to analyze the gas-driven fracture propagation phase of wedge- and penny-shaped cracks. The Nilson model was then extended by Daehnke et al. [23] to analyze gas flow in propagating conical cracks. The Nilson model was also implemented by Cho et al. [24] to investigate the dynamic fracture process analysis of rock subjected to stress wave and gas pressurization based on dynamic FEM code and one-dimensional gas flow. Ma and An [25] implemented the Johnson–Holmquist (J–H) material model into the commercial software LS-DYNA to simulate the blasting-induced rock fractures. Mohammadi and Pooladi [26] developed a 2D coupled two-mesh interaction model for blast gas flow through fractured and fragmented solid media, in which the solid domain is represented by a combined finite/discrete element technique to model progressive cracking and fragmentation and any potential normal and frictional contacts during the cracking and post-cracking phases. Using one-dimensional transient flow model for gas flow in the cracks and displacement extrapolation method for the stress intensity factor calculation, Goodarzi et al. [15] proposed an extended finite element method-based model which can quantitatively simulate fractures propagation around the borehole as a quasistatic phenomenon.

In order to overcome the limitations of each fracturing method, combination of different fracturing techniques have also been studied. Hou et al. [6] proposed the pulse gas fracturing experiment of low-permeability coal on the laboratory scale, which shows that pulse gas fracturing can improve the fracture propagation in coal seams because low pulse loading can produce more secondary cracks. Experiment of hydraulic fracturing after water pressure control blasting on cement mortar samples was carried out by Huang et al. [27], which indicated that this method can produce not only multiple macroscopic main cracks along borehole axial and radial directions but also local hydraulic cracks in multiple directions and of different types.

Despite the disadvantages of short crack lengths, HEGF has its advantages including low costs, the possibility of having multiple fractures without causing an extensive damage, and the creation of radial fractures with fracture orientations independent of formation stress anisotropy and heterogeneity [15, 28]. The multipulse HEGF could also bring multiplier fracture length by concatenating series propellants with different burning rate [16]. Therefore, HEGF can be used as a valuable prefracturing tool to reduce surface pressure requirements during HF and GF, due to its effectiveness to reduce total friction near wellbore [14].

As mentioned above, the existing models of HEGF focus on the dynamic loading or gas pressurization loading rather than the combined actions of them. Researches on the combined actions of HEGF and GF are few. In view of this, the major objective of this paper is developing a damage-based stress wave propagation-static mechanical equilibrium-gas flow coupling model to study the crack evolution in HEGF. The outline of this paper is as follows. In Section 2, the general governing equations for the coupled gas flow-mechanical damage model for HEGF is established, in which dynamic stage for stress wave and quasistatic stage for gas pressurization are included. In Section 3, the numerical model and determination of mesomechanical parameters in finite elements are described in detail. In Section 4, numerical simulations on crack evolution under HEGF, GF, and HEGF+GF are carried out, and the impact of in situ stress conditions on crack evolution is discussed further.

#### 2. Governing Equations

In this section, the governing equations for the coupled dynamic mechanical-quasistatic mechanical-gas flow system are derived, based on momentum conservation and mass conservation, respectively. The rock damage evolution under the coupled conditions and its impact on mechanical parameters are also considered. These derivations are based on the following assumptions: (a) The stimulated reservoir is elastic and isotropic. (b) The dynamic loading in HEGF is simplified as elastic stress wave. (c) The mechanical stress in rock is calculated based on finite element method. (d) Damage-based elastic model is used for damage evolution of each element. (e) The gas flow in the intact reservoir and cracks follows Darcy’s law. (f) The gas released in HEGF is regarded as ideal and the temperature variation is ignored.

##### 2.1. Mechanical Equilibrium Equation

During explosive blasting, high-pressure air blasting, and high-energy gas fracturing, two types of loadings are applied on the borehole wall, first dynamic stress wave with short duration, followed by quasistatic gas pressurization with longer duration. The modified Navier’s equation (Equation (1)) is used here to express the mechanical equilibrium of rock subjected to dynamic loading [19]:
where is elastic modulus (Pa), is Poisson’s ratio, () is dynamic displacement (m), is rock density (kg m^{−3}), and is time (s).

Based on poroelastic theory, the modified Navier’s equation (Equation (2)) is used here to express the mechanical equilibrium of rock subjected to quasistatic gas pressurization:
where () is static displacement (m), is gas pressure, is Biot’s coefficient, and is the components of the net body force in the direction (N m^{−3}).

Equation (2) can also be used to express the mechanical equilibrium of underground rock subjected to in situ stress loading with gas pressure of zero. Therefore, for underground rock subjected to both dynamic loading induced by blasting stress wave and quasistatic loading induced by in situ stress, Equations (1) and (2) can be used together to obtain the total displacement, strain, and stress based on superposition principle.

##### 2.2. Gas Flow Equation

During explosive blasting, high-pressure air blasting, and high-energy gas fracturing, a huge volume of gas with high pressure is released instantaneously. The gas would then propagate into the porous rock and blasting-induced cracks, which is extremely complicated. Based on continuum mechanics, the gas propagation is simplified as gas flow process in porous medium with modified porosity and permeability of cracks.

If the gas flow follows Darcy’s law,

The mass conservation equation for the gas flow process in porous rock can be expressed as [29]
where is the porosity of rock, is the dynamic viscosity of gas (N s m^{−2}), is the permeability of rock (m^{2}), is the gas density (kg m^{−3}) at standard condition, is the gas pressure (Pa) at standard condition, and is a source term (kg m^{-3} s^{-1}).

##### 2.3. Damage Evolution Equation

As illustrated in Figure 1, the elastic damage-based constitutive law is used to evaluate the mesomechanical damage of elements with the maximum tensile stress criterion for tensile damage and the Mohr–Coulomb criterion for shear damage. The criteria can be expressed as [30–35] where , , , and are the first principal stresses (positive for tensile), the third principal stresses (negative for compressive), the uniaxial tensile strength, and the uniaxial compressive strength, respectively. is the internal frictional angle (°).

According to Figure 1, the damage variable can be calculated as [30–35]

More details can be found in previous publications [30–35].

##### 2.4. Effect of Damage on Mechanical Parameters

The impact of elements damage on mesomechanical parameters is also considered. The elastic modulus and strength of a mesoelement degrade monotonically as damage evolves and can be expressed as [30–35] where subscript “0” represent parameter of the undamaged element.

The permeability and porosity are also correlated to the damage variable according to the following exponential function:
where and are the permeability (m^{2}) and porosity of the undamaged element, respectively.

It is noted that the permeability and porosity are actually not constant when is not equal to zero. The permeability will increase a couple of orders of magnitude when damage occurs, which is related to the aperture of cracks. In fact, governing equations based on the theory of gas dynamics may be more appropriate for gas moving at high velocity in cracks, in which the break of elements should be considered and more complexity will be introduced in numerical simulations. The gas fracturing model proposed here is based on rock damage theory and mainly focus on the propagation of cracks, which is induced by the existence of high gas pressure in cracks. Hence, the permeability and porosity with a relatively high and constant value when damage occurs are used here to simulate gas moving at high velocity in cracks and obtain gas pressure distributions in cracks based on finite element method, damage theory, and gas flow equations.

The rock damage evolution is assumed to be isotropic here. So, mechanical parameters, including , , , , and , are all scalars. It is noted that the isotropic damage model is used for damage evolution of each element. When the mesoscale element damages, the permeability of the damaged element will increase dramatically. When a mass of mesoscale elements damage, a macro- and anisotropic fracture for gas flow will be formed. Hence, the isotropic damage model for mesoscale elements is used which can result in macro- and anisotropic permeability in a specified direction.

#### 3. Model Setup

##### 3.1. Numerical Model

A representative model geometry is shown in Figure 2 containing a wellbore with the diameter of 0.2 m set centrally within a shale domain. For high-energy gas fracturing, a dynamic stress wave is applied first on internal borehole surface, followed by quasistatic gas pressure . For gas fracturing, however, only quasistatic gas pressure is applied on the borehole surface. Crustal stresses of and are applied on the outer boundary of the model.

##### 3.2. Determination of Mesomechanical Parameters

The macromechanical parameters of shale reported by Hou et al. [6] with an elastic modulus of 17 GPa, uniaxial compressive strength (UCS) of 107 MPa, tensile strength of 10 MPa, and Poisson’s ratio of 0.32 are chosen in this study. In addition, a Weibull distribution as defined in the following probability density function is used here to characterize the heterogeneity of shale [30–35]: where is the homogeneity index reflecting the homogeneity degree of shale, is the scale parameter related to the average value of the mechanical parameter (for example, tensile strength), and is the mesomechanical parameter of each mesoscopic element in the numerical model. More detailed description about this function can be referred in [30–35].

For a numerical sample (50 mm wide and 100 mm high) with mesomechanical parameters listed in Table 1, its elastic modulus evolution and complete stress–strain curve under uniaxial compressive condition is shown in Figures 3 and 4, respectively. The heterogeneous elastic modulus distribution is specified initially according to Weibull distribution (Figure 3, step 0). As the axial displacement increases, the elastic modulus decreases gradually (Figure 3, steps 50–70), leading to the formation of macrofracture and unstable failure of the sample (Figure 3, steps 75–80). Finally, the macroelastic modulus of 16.9 GPa and UCS 106.5 MPa were obtained (Figure 4) which are very close to the macromechanical parameters. Therefore, the mesomechanical parameters listed in Table 1 are used in later simulations.

##### 3.3. Determination of Loading Curves

Figure 5 shows the loading curves for high-energy gas fracturing and gas fracturing.

**(a)**

**(b)**

**(c)**

**(d)**

The loading curve for dynamic loading stage in high-energy gas fracturing (HEGF) as shown in Figure 5(a), with peak pressure of 300 MPa and duration of 10 ms, is a simplification of the actual situation based on reference [14]. The dynamic load applied on the borehole wall in HEGF will propagate gradually from the borehole wall into the surrounding rock with time, resulting in crack initiation and propagation around the borehole. In order to reproduce crack evolution in numerical simulations, the propagation of stress wave is divided into many steps. The stress distribution at each step is calculated based on governing equations (Equations (1) and (2)), and the damage variable at each step is calculated based on stress distribution and damage criteria in Equations (5) and (6). So the crack evolution is reproduced step by step. Apparently, the more steps the loading curve is divided into, the more precise the crack evolution is. In current simulations, a total time of 20 ms is calculated to ensure the entire stress wave can propagate from the borehole to the outer boundary of model. Considering the computer capacity, 20 ms is divided into 100 steps to simulate crack evolution.

The loading curve for quasistatic loading stage in HEGF is shown in Figure 5(b). For the similar reason mentioned above, a total of 100 steps (steps 101-200) are calculated to reproduce crack evolution at this stage. The actual gas pressure in the borehole, as shown in Figure 5(b), is calculated based on the initial gas pressure (the peak pressure of stress wave, 300 MPa, is used here), crack distribution (for example, Figure 6, step 100) and governing equation (Equation (4)). Then, the stress state at each step is calculated based on governing equation (Equation (2)), and the damage variable at each step is calculated based on stress distribution and damage criteria in Equations (5) and (6). So the crack evolution is reproduced step by step.

The loading curve for gas fracturing (GF) is shown in Figure 5(c). For the similar reason mentioned above, a total of 100 steps, numbered in order as 201-300, are calculated to reproduce crack evolution in GF. For simplicity in simulations, gas pressure with a constant increasing rate of 5 MPa/step is applied in the borehole, which is equivalent to a constant gas injection rate at atmospheric pressure. Apparently, the injected gas will enter into the cracks around the borehole rapidly and the gas pressure in the borehole will decrease accordingly. The actual gas pressure value in the borehole, as shown in Figure 5(c), is calculated based on gas injection rate and crack distribution (for example, Figure 7, step 250).

The loading curve for HEGF+GF (combined actions of high-energy gas fracturing and gas fracturing) is shown in Figure 5(d). The loading curves at steps 1-200 are the same as the curves in Figures 5(a) and 5(b). The loading curve at steps 201-300 is different from the curve in Figure 5(c) because of different crack distributions.

##### 3.4. Boundary Conditions

Based on the model setup shown in Figure 2, the crustal stress conditions of and applied on the outer boundary of the model in mechanical field can be expressed as

The dynamic load shown in Figure 5(a), which is applied on the borehole wall in the mechanical field, can be expressed as

The quasistatic load shown in Figures 5(b) and 5(c), which is applied on the borehole wall in the mechanical field, can be expressed as

The gas pressure shown in Figures 5(b) and 5(c), which is applied on the borehole wall in the gas flow field, can be expressed as

Equations (1–16) represent a coupled nonlinear model for HEGF. More details about the implementation of these partial differential equations with COMSOL Multiphysics [36] have been presented in previous publications [37].

#### 4. Simulation Results and Discussion

##### 4.1. Crack Evolution in HEGF

In this section, loading curves at steps 1-200 shown in Figure 5 are applied to simulate the damage evolution during HEGF.

Figure 6 shows crack evolution induced by stress wave at the dynamic stage (steps 1-100) of HEGF. The maximum principal stress is shown in the first row, where positive value is for tension and negative value for compression. The stress wave applied on internal borehole surface can result in circumferential tensile stress around the borehole (step 25). The circumferential tensile stress then propagates gradually deep into the surrounding rock (steps 50 and 100) with its value decreasing accordingly. The damage zone distribution is shown in the second row, where damage variable for no damage, for tensile damage, and for shear damage. With stress wave propagation from the borehole wall into the surrounding rock, a tensile damage zone is created first around the borehole (step 20) and expands gradually to 3-5 times of borehole radius (steps 50 and 100). The elastic modulus distribution is shown in the third row, where the initial elastic modulus distribution (at step 10. Elastic modulus distribution at step 10 is the same as the one at step 1 because no damage occurs before step 10) is specified according to Weibull distribution with mesomechanical parameters listed in Table 1. The elastic modulus degrades monotonically as damage evolves (steps 25-100).

Figure 8 shows crack evolution at the quasistatic gas pressurization stage of HEGF (at steps 101-200), where the distributions of permeability, gas pressure, maximum principal stress, damage zone, and elastic modulus are presented.

At dynamic loading stage of HEGF, damage zone is created around the borehole (Figure 6, step 100) in which permeability increase dramatically (Figure 8, the first raw). Then the high-pressure gas (Figure 8, the second raw), produced from burning a propellant during HEGF, can squeeze into the damage zone and cause circumferential tensile stress around the damage zone or tensile stress concentration at the crack tips (Figure 8, the third raw, the maximum principal stress distributions), leading to the producing of new cracks or extending of existing cracks (Figure 8, the forth raw and fifth raw). The gas pressure decrease gradually (Figure 5, steps 101-200) with its squeezing and expanding into the cracks. Finally, a total of 11 radial cracks with length of 2–4 m are produced (Figure 9(a)).

**(a)**

**(b)**

**(c)**

Figure 10 shows the variation of gas pressure in the borehole and damage zone area at the quasistatic gas pressurization stage of HEGF. Gas pressure with initial value of 300 MPa, being identical to the peak pressure of the stress wave, is applied on the borehole. Then, the gas expands instantly into the damage zone around the borehole and the gas pressure declines dramatically. Hence, the gas pressure at step 101 is 122 MPa rather than 300 MPa. The damage zone area increases rapidly and the gas pressure in the borehole decreases accordingly from steps 101 to 110. After step 134, the damage zone area does not increase anymore and the gas pressure keeps at 48 MPa.

##### 4.2. Crack Evolution in GF

In this section, loading curves at steps 201-300 shown in Figure 5 are applied to simulate the damage evolution during GF. Here, the initial step is set to 201 to be consistent with the loading curves in Figure 5**.**

The distributions of maximum principal stress, permeability, gas pressure, damage zone, and elastic modulus are presented in Figure 7.

The static gas pressure applied on the internal borehole surface can cause circumferential tensile stress and results in tensile damage around the borehole (Figure 7, step 225). Then, the gas squeezes into the damage area, causing tensile stress concentration at the crack tips and creation of tensile cracks (Figure 7, steps 250 and 275). Finally, a total of 4 main cracks with length of 6.5–8 m are produced (Figure 7, step 300). As shown in Figure 9(b), crack No.IV includes another two branch cracks, No.IV-1 and No.IV-2.

Figure 11 shows the variation of gas pressure in the borehole and damage zone area under gas fracturing. At steps 201-215, the gas pressure increases linearly as a result of the constant gas injection rate in the borehole. At step 216, damage zone and radial cracks occurs around the borehole, providing space for gas expanding and leading to the rising rate of gas pressure decreases. With cracks extending, the gas pressure decreases rapidly from step 221 to 232. At this stage, gas injection rate is much less than crack propagation speed. From step 233 to 300, gas pressure is in the decline and fluctuate stage, and at this stage, the damage zone area increases linearly on the whole. Finally, the gas pressure decrease to 56 MPa at step 300.

##### 4.3. Crack Evolution under Combined Actions of High-Energy Gas Fracturing and Gas Fracturing (HEGF+GF)

In this section, damage zone evolution under combined actions of high-energy gas fracturing and gas fracturing (HEGF+GF) is simulated. Gas is injected into the borehole after HEGF to fracture the shale reservoir further. Gas injection rate is same as the one of gas fracturing in Section 4.2. The distributions of permeability, damage zone, and elastic modulus at step 200 in Figure 8 is specified as the initial status of gas fracturing in HEGF+GF.

Figure 12 shows the evolution of maximum principal stress, permeability, gas pressure, damage zone, and elastic modulus, respectively. Similar to the damage evolution at the quasistatic stage of HEGF in Figure 8 and GF in Figure 7, the gas injected into the borehole can squeeze into the connected cracks with higher permeability and cause tensile stress concentration at the crack tips (Figure 12, step 225), leading to extension of existing cracks (steps 250 and 275). Finally, a total of 11 radial cracks with length of 2-11.5 m are produced (Figure 12, step 300; Figure 9(b)).

Figure 13 shows the variation of gas pressure in the borehole and damage zone area at the gas fracturing stage of HEGF+GF. The gas produced in HEGF is assumed to be still in the fractured shale reservoir and the gas distribution at step 200 in Figure 8 is regarded as the initial status of gas pressure. Hence, the initial gas pressure in Figure 13 is 48 MPa. Then, gas is injected into the borehole with a constant injection rate and gas pressure increases slowly. At steps 201-230, the damage zone increase rate rises, while the gas pressure decrease rate declines gradually. After step 230, the damage zone area increases linearly, while the gas pressure keeps at 57 MPa.

Figure 9 shows the final cracking pattern in HEGF, GF, and HEGF+GF, respectively.

In HEGF, a total of 11 cracks, with length of 2.3-4 m for the main cracks Nos.1-5 and 1.5-1.9 m for the other minor cracks Nos.6-11, are produced.

In GF, a total of 4 main cracks Nos.I-IV with length of 6.5–8 m, including two branch cracks Nos.IV-1 and IV-2 from crack No.IV, are produced.

In HEGF+GF, the main cracks Nos.1-5 formed in HEGF extend further with their length increasing to 5.2-11.5 m, while the minor cracks Nos.7, 8, 10, and 11 formed in HEGF do not extend anymore. The other two minor cracks Nos.6 and 9 formed in HEGF also extend further with their length increasing to 3.9 and 3.2 m, respectively.

The comparison of final cracking pattern in Figure 9 reveals that HEGF produce more cracks while the crack length is relatively short. GF can produce longer cracks while the amounts of crack are fewer. However, more cracks with longer length are produced in HEGF+GF, which enhances fracturing effect significantly. Apparently, the high-pressure gas applied on the borehole at HEGF stage is advantageous for the formation of more minor cracks around the hole, and the injection of huge amounts of gas at GF stage is advantageous for further extension of minor cracks.

##### 4.4. Impact of In Situ Stress Conditions on Crack Evolution

Crustal stresses of and are applied on the outer boundary of the model in Sections 4.1–4.3. In this section, crustal stress conditions of , and , are applied, respectively, to study the effect of stress conditions on crack evolution under combined actions of HEGF+GF.

As shown in Figure 14, although a lot of minor cracks are formed around the borehole at the dynamic stage (step 100) and quasistatic gas pressurization stage (step 200) of HEGF, only cracks in the direction extend further at the GF stage (step 300). Finally, only 3 cracks with a length of 9–10 m in the direction, which is the same as the maximum crustal stress direction of here, are produced.

Similarly, under crustal stress conditions of , (shown in Figure 15), although a lot of minor cracks are formed around the borehole at dynamic stage (step 100) and quasistatic gas pressurization stage (step 200) of HEGF, only cracks in the direction extend further at the GF stage (step 300). Finally, only 2 cracks with length of 10 m in the direction, which is the same as the maximum crustal stress direction of here, are produced.

Figures 14 and 15 indicate that the main cracks extend only in the maximum crustal stress direction, while the development of minor cracks in the other directions are repressed.

Figure 16 shows crack evolution in HEGF+GF under crustal stress conditions of and , in which stress waves of 450 MPa and 600 MPa are applied, respectively. Compared with crack evolution in Figure 15 in which stress wave of 300 MPa is applied, although the peak pressure of stress wave in Figure 16 increases and more minor cracks are formed around the borehole at dynamic stage (step 100) and quasistatic gas pressurization stage (step 200) of HEGF, still only cracks in the direction extend further at the GF stage (step 300). Apparently, the increment of loading force in the borehole has little effect on crack development in directions other than the direction of maximum crustal stress.

Similar results can be found in rock blasting and hydraulic fracturing [21, 25, 27]. Multiple fracturing technology, which is used to obtain maximum stimulated reservoir volume (SRV) in hydraulic fracturing [38, 39], may be also an effective way to produce more cracks in HEGF, GF, and HEGF+GF under heterogeneous crustal stresses conditions.

#### 5. Conclusions

In this work, a damage model for high-energy gas fracturing, including dynamic stage and quasistatic gas pressurization stage, is established. Numerical simulations on crack evolution under high-energy gas fracturing (HEGF), gas fracturing (GF), and combined actions of high-energy gas fracturing and gas fracturing (HEGF+GF) are carried out, respectively. The simulation results indicate that (1) HEGF produce more cracks while the crack length is relatively short. GF can produce longer cracks while the amounts of crack are fewer. However, more cracks with longer length are produced in HEGF+GF, which enhances fracturing effect significantly. (2) The high-pressure gas applied on the borehole at HEGF stage is advantageous for the formation of more minor cracks around the hole, and the injection of huge amounts of gas at the GF stage is advantageous for further extension of minor cracks. (3) The main cracks extend only in the maximum crustal stress direction, while the development of minor cracks in the other directions are repressed. (4) The increment of loading force in the borehole has little effect on crack development in directions other than the direction of maximum crustal stress.

#### Data Availability

The data used in the manuscript is from the simulation results by MATLAB and COMSOL Multiphysics. And we wrote MATLAB codes independently to conduct these simulations. Therefore, you can contact Yu Bai (baiyu1025@163.com) if you want to learn more about our work.

#### Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

#### Acknowledgments

This research was funded by the National Key R&D Program of China (grant number 2017YFC1503103), National Natural Science Foundation of China (grant numbers 51578347 and 51878420), Key Research and Development Program of Liaoning Province (grant number 2017231010), and Shenyang Science and Technology Program (grant number 18-013-0-31).