Geofluids
Volume 2019 (2019), Article ID 2624716, 13 pages
https://doi.org/10.1155/2019/2624716

Research Article

Numerical Simulations of Fracture Propagation in Jointed Shale Reservoirs under CO₂ Fracturing

Qi Zhang,1,2 Dan Ma,1,2 Jiangfeng Liu,² Jiehao Wang,³ Xibing Li,¹ and Zilong Zhou¹

¹School of Resources and Safety Engineering, Central South University, Changsha, 410083 Hunan, China
²State Key Laboratory for Geomechanics and Deep Underground Engineering, China University of Mining and Technology, Xuzhou, 221116 Jiangsu, China
³Department of Energy and Mineral Engineering, EMS Energy Institute and G3 Center, Pennsylvania State University, University Park, PA 16802, USA

Correspondence should be addressed to Dan Ma; dan.ma@csu.edu.cn

Received 8 November 2018; Revised 5 February 2019; Accepted 5 March 2019; Published 7 April 2019

Guest Editor: Qianbing Zhang

Copyright © 2019 Qi Zhang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Water-based hydraulic fracturing for the exploitation of shale gas reservoirs may be limited by two main factors: (1) water pollution and chemical pollution after the injection process and (2) permeability decrease due to clay mineral swelling upon contact with the injection water. Besides, shale rock nearly always contains fractures and fissures due to geological processes such as deposition and folding. Based on the above, a damage-based coupled model of rock deformation and gas flow is used to simulate the fracturing process in jointed shale wells with CO₂ fracturing. We validate our model by comparing the simulation results with theoretical solutions. The research results show that the continuous main fractures are formed along the direction of the maximum principal stress, whilst hydraulic fractures tend to propagate along the preexisting joints due to the lower strength of the joints. The main failure type is tensile damage destruction among these specimens. The preexisting joints can aggravate the damage of the numerical specimens; the seepage areas of the layered jointed sample, vertical jointed sample, and orthogonal jointed sample are increased by 32.5%, 29.16%, and 35.05%, respectively, at time  s compared with the intact sample. The preexisting horizontal joints or vertical joints promote the propagation of hydraulic fractures in the horizontal direction or vertical direction but restrain the expansion of hydraulic fractures in the vertical or horizontal direction.

1. Introduction

In recent years, the application of hydraulic fracturing has significantly increased the production of oil and natural gas. It is estimated that hydraulic fracturing has improved recoverable reserves of oil by at least 30% and of natural gas by 90% [1]. In addition, more than 60% of oil and gas wells need to be fractured first, especially for unconventional gas resource deposited in deep underground shale layers with extremely small permeability (usually less than 1 mD) [2, 3]. Hydraulic fracturing may increase permeability of shale by three to five orders of magnitude [4]. Thus, the production of the fractured gas wells increases dramatically. In all, hydraulic fracturing introduced in oil and gas fields has revolutionized gas production around the world [5, 6].

However, the water-based fracturing technology has some limitations and environmental concerns. First, hydraulic fracturing consumes large amounts of water. According to a report [7], the water use for hydraulic fracturing accounts for 9% of total freshwater consumption in Texas. The large consumption of water will restrict the oil and gas reservoir developments in a water-deficient area [8, 9]. Second, due to the large amount of water and chemical reagents used in the hydraulic fracturing, it may cause potential water pollution and chemical pollution if the treatment of flowback fluids with chemical reagents is insufficient [10–12]. Third, for the reservoirs containing clay minerals, the permeability may decrease after water injection, thereby decreasing the production of gas reservoirs [13–15]. The main reason is that, when the hydration minerals meet injection water, clay minerals can swell and result in blockage of seepage channels [16]. All of these disadvantages of hydraulic fracturing promote the study and development of waterless fracturing [17].

Several waterless fracturing technologies have been introduced in oil and gas industries over the past few decades, including oil-based fracturing, N₂ fracturing, and CO₂ fracturing [18]. Oil-based fracturing was first used in Colorado, Texas, and Kansas in late 1940s [19]. Compared with hydraulic fracturing, it could be conducted in frozen areas. However, oil-based fracturing is expensive and may impair the effective permeability of wells [20]. N₂ fracturing and CO₂ fracturing are the two most popular fracturing methods because they are more economical and efficient compared with hydraulic fracturing [21]. According to engineering production data [22], the production of the reservoirs stimulated by CO₂ is 1.9 times as much as the production of those stimulated by N₂. The laboratory experiment results indicate that the gas with lower viscosity, higher diffusivity, and lower surface tension can penetrate into smaller pore space to create more complex fracture networks compared with hydraulic fracturing [23, 24]. In addition, the fracture surface created by gas fracturing has a larger roughness and complexity, resulting in a greater increase in permeability [4].

In addition, shale reservoirs always contain kinds of joints caused by the geological deposition and folding [25–27]. The existing joints have a significant influence on the initiation and propagation of the induced hydraulic fractures [28–30]. The fracture networks of jointed reservoirs may be very complex due to the reopening of the existing joints, the expansion of hydraulic fractures, and the intersection between joints and hydraulic fractures [31]. Nitrogen fracturing experiments were conducted on shale samples vertical and parallel to the bedding plane; the results indicated that a relative complex fracture surface is formed in the shale sample vertical to the bedding plane [32]. He et al. [33] performed hydraulic fracturing on shale with bedding planes; the results showed that the bedding planes in shale formation have a significant influence on the propagation of hydraulic fractures. However, the mechanism of the fracture initiation and propagation in kinds of jointed reservoirs is not well investigated. It is important to learn the distribution of fracture networks for the successful design of stimulation in jointed reservoirs.

To this end, the numerical tools COMSOL and MATLAB are used to simulate the hydraulic fracture propagation driven by injection fluids in several jointed reservoirs. The distribution of fracture networks and the development of horizontal and vertical fracture radii are studied in this work.

2. Governing Equations

In the numerical simulation, CO₂ is injected into the borehole. Then, the rock mass begins to fracture with the increasing injection pressure. The process of CO₂ fracturing involves solid deformation and fluid seepage. In this part, a series of governing equations are set up for solid mechanic field and flow field. Besides, damage equations are introduced to describe the destruction of the calculation elements.

2.1. Rock Deformation and Damage Evolution Equations

In this work, shale rock is assumed as an elastic continuum material, whose constitutive relation satisfies with the physical equation of elasticity. It should be noted that the influence of pore pressure on stress distribution is also considered in the equation. Thus, the modified physical equation can be induced as

where is the total stress tensor, is the total strain tensor, is the shear modulus of rock, is the elastic modulus of the rock, is Poisson’s ratio of the rock, is the volumetric strain, is the Kronecker delta, is the Biot coefficient, and is the pore pressure.

The relationship between strain and displacement is expressed by a geometric equation as follows:

where and are the components of displacement in and directions, respectively.

Substituting the modified physical equation (1) and the geometric equation (2) into the equilibrium equation, then the modified Navier-type equation is induced as

where is the component of the net body force.

Since the initiation and propagation of hydraulic fractures are studied in this work, a damage model is introduced to characterize the damage condition during the injection process. The damage model is used to determine whether shale damage occurs after every calculation step. For the calculation element, when the stress state meets the maximum tensile stress criterion or the Mohr-Coulomb criterion, the tensile crack or shear crack occurs. It should be noted that the tensile crack is first generated, because the compressive strength is ten times greater than the tensile strength. Equations (4) and (5) are the maximum tensile stress criterion and the Mohr-Coulomb criterion, respectively:

where and are the first and third principle stresses; and are the tensile strength and compressive strength of rock, respectively; and is the internal friction angle.

When elements start to be damaged, the elastic modulus reduces correspondingly according to damage theory [34]. The evolution of elastic modulus is defined as

where is the initial elastic modulus of rock and is the damage variable and is calculated as [35–37]

where and are the first and third principal strains and and are the tensile strain and compressive strain, respectively.

2.2. Gas Flow Equation

Gas flow equations are defined to describe the injection gas flow in this part. The gas continuity equation during gas transportation is defined as

where is the gas mass per volume of rock, is the density of the injection gas, is the seepage velocity of the gas, is the source origin, and is the time variable.

On the basis of Darcy’s law, the seepage velocity of gas is shown as

where is the permeability of the rock and is the dynamic viscosity coefficient.

Assuming that shale rock is saturated by CO₂ after the injection, the gas content per volume of rock can be defined as , and is the porosity of the rock. Injected CO₂ gas enters the supercritical state when the pressure exceeds 7.56 MPa at the temperature 76.8 degrees Celsius [38]. When the injection CO₂ is transformed from the gaseous state to the supercritical state, density and viscosity change dramatically under the different pressure. The evolution of density and viscosity of CO₂ varying with pressure is shown in Figure 1. The relationship between density, viscosity, and pressure can be described by interpolating function in the model calculation. Thus, the first item in equation 7 is induced to equation (10) [39] as

Figure 1: The evolution of density and viscosity of CO₂ with pressure (°C).

The density of CO₂ changes rapidly with pressure whilst the porosity of shale changes slightly during hydraulic fracturing. Thus, is ignored since is much smaller than , where is the compressibility coefficient of CO₂ and can be calculated from Figure 1.

Putting equations (9) and (10) into equation (8), the gas continuity equation can be induced into equation (11) as

Considering the influence of solid mechanics on the evolution of porosity, the dynamic evolution of porosity can be described by [40]

where is the initial porosity of shale rock, is the initial volumetric strain, is the initial pore pressure, and is the bulk modulus of rock grains.

2.3. Characterization of Rock Heterogeneity

Since rock is heterogeneous material with natural fissures and induced fissures [41], Weibull distribution function is introduced to represent the heterogeneity of shale rock in this part [42]. In this work, parameters, such as elastic modulus and strength, are assumed to conform to the Weibull distribution and produced in MATLAB. Weibull statistical distribution and the probability density function are defined as

where is the mechanical parameter of rock, is the average value of , and is the coefficient of heterogeneity.

3. Numerical Simulation Setup

3.1. Model Geometry

Before the simulation, the models for gas fracturing, boundary settings, and parameters are introduced in this part. The 2D calculation model here is a 2D plane cross-section through a horizontal well, which is a 2D plane square with a length of 0.2 m and a borehole with a diameter of 0.04 m drilled in the center. The distribution of preexisting horizontal joints, vertical joints, and orthogonal joints is showed in Figure 2. The thickness of the joints is 1 mm among these models. For mechanical boundary settings, the load applied on the left boundary is equal to that on the upper boundary, the value is 5 MPa, whilst the right boundary and lower boundary are roller boundaries; these boundary settings are the same for the intact sample, layered jointed sample, vertical jointed sample, and orthogonal jointed sample. The gas is injected from the inner boundary with the injection rate of 0.0106 m³/s. As for seepage boundaries, there is no flow boundary on the outer boundaries. The parameters for the rock matrix and joints are listed in Table 1; it should be noted that the elastic modulus, tensile strength, and compressive strength of joints are only half of those of the rock matrix.

(a)

(b)

(c)

(d)

Figure 2: Numerical samples with different kinds of joints: (a) intact sample, (b) sample with layered joints, (c) sample with vertical joints, and (d) sample with orthogonal joints.

Table 1: Mechanical parameters for the simulations.

Mechanical parameters
Rock matrix
Joint
Unit

Elastic modulus
36
18
GPa

Poisson’s ratio
0.225
0.25
—

Rock density
2600
2600
kg/m³

Compressive strength
62
31
MPa

Tensile strength
8.2
4.1
MPa

Initial rock permeability
10⁻¹⁸

m²

Initial rock porosity
0.01
0.015
—

Biot coefficient
0.2
0.3
—

Internal fraction angle
0.117
0.117
rad

3.2. Numerical Implementation of the Model

In the front parts, the coupling equations and model geometry are established. Due to the complexity of the coupling relationship between solid deformation and fluid flow, these equations are difficult to be directly calculated. Thus, the finite element method is adopted to solve coupling equations via COMSOL Multiphysics and MATLAB. The distribution of stress and pore pressure is obtained to further discuss the condition of the numerical sample. The flow chart of numerical calculation procedures is showed in Figure 3.

Figure 3: Calculation procedures of the numerical model.

3.3. Model Validation

It is important to validate the effectiveness of the model before it is used to simulate the fracture propagation in the jointed reservoirs. The classical theoretical solution for forecasting the breakdown pressure [43] is shown as

where is the breakdown pressure, is the tensile strength of the rock, and are the maximum and minimum principal stresses, respectively, and is the initial pore pressure.

The numerical specimen used in this part is the same with the intact sample shown in Figure 2(a). In this part, the initial pore pressure is 1 MPa, the average tensile strength and compressive strength of the specimen are 5.9 MPa and 59.4 MPa, respectively, and the average elastic modulus of the numerical specimen is 35.9 GPa. The horizontal stress is fixed at 25 MPa, and the vertical stress is varying from 10 MPa to 25 MPa. A tectonic stress coefficient is defined in this part, and the breakdown pressure versus the tectonic stress coefficient is calculated in this part.

A comparison of the breakdown pressure with different tectonic stress coefficients by the theoretical solution and numerical solution is shown in Figure 4. The results show that the numerical simulation results agree well with results obtained from the theoretical equation, indicating that the proposed model is suitable to simulate the hydraulic fracture propagation under CO₂ fracturing.

Figure 4: Comparison of the theoretical solution and numerical solution of breakdown pressure under different tectonic stress coefficients.

4. Results and Discussion

4.1. Fracture Propagation in the Intact Sample

The hydraulic fracture initiation and propagation in an intact sample without joints are shown in Figure 5. To distinguish tensile cracks and shear cracks, the damage in tension is defined as negative whilst that in shear is defined as positive. It can be seen that the fractures first appear around the borehole and propagate gradually to the surrounding rock with the increase of injection gas pressure [44]. Eventually, several main fractures are formed uniformly in the rock sample. As shown in Figure 5, hydraulic fractures are principally formed in the tensile mode in the intact sample during gas fracturing, and the number of damaged elements in the tensile mode accounts for 93.6% of the total number of damaged elements at time s.

Figure 5: Distribution and evolution of the fracture in the intact sample.

Figure 6 shows the development of the seepage area versus injection time. It should be noted that the cracks (damage area) are generated from the destruction of the elements. The seepage area keeps as 0 m² due to the low gas pressure from  s to  s. Then, the seepage area increases slowly from 3.27⁻⁵ m² to 3.44⁻³ m² from  s to  s with fracture initiation and gradual propagation. Besides, it increases dramatically from 3.44⁻³ m² to 6.93⁻³ m² from  s to  s, indicating that unstable fracture propagation occurs. The whole process of the evolution of the seepage area can be fitted by an exponential function as shown in Figure 6.

Figure 6: Development of the seepage area versus injection time in the intact sample.

The development of horizontal and vertical fracture radii versus injection time in the intact sample is presented in Figure 7. It can be seen that hydraulic fractures initiate at time  s. As the injection gas pressure in the borehole further increases, the fracture radii in horizontal and vertical directions increase at a similar rate, reaching 6.3⁻² m and 6.49⁻² m at time  s, respectively. The slight difference of two curves shown in Figure 7 may result from the heterogeneity of mechanical parameters assigned to the sample.

Figure 7: Development of horizontal and vertical fracture radii versus injection time in the intact sample.

4.2. Fracture Propagation in the Layered Jointed Sample

In the shale gas wells, the surrounding rock contains kinds of joints and fissures caused by the geological deposition and folding. The existence of these joints and fissures has a significant influence on the propagation of hydraulic fractures. Based on the above, hydraulic fracture propagation in kinds of joint samples will be discussed in the following parts. The fracture propagation in layered jointed samples is researched in this part.

The distribution and evolution of the fracture in the layered jointed sample are presented in Figure 8. At time  s, the sample begins to fracture. Hydraulic fractures first emerge around the drilling hole at time  s due to stress concentration in the inner boundary. With the increase of the injection time, hydraulic fractures propagate to the surrounding rock. Finally, fractures propagate along the diagonal direction and horizontal direction forming several main fractures, as shown in Figure 8. The distribution of the fracture is different from that in the intact sample. This result from the fracture propagation favors the horizontal preexisting joints along which the least energy is dissipated due to their low strength and elastic modulus.

Figure 8: Distribution and evolution of the fracture in the layered jointed sample.

The seepage area begins to increase from time  s; then, it increases slowly from 7.92⁻⁶ m² to 3.22⁻³ m² from time  s to  s due to the low injection gas pressure at the initial stage. With the increase of injection gas pressure, the seepage area increases rapidly from 3.22⁻³ m² to 9.18⁻³ m² from time  s to  s. The seepage area increases by 32.5% in this case compared with that in the intact sample at time  s. An exponential function is used to describe the development of the seepage area (Figure 9).

Figure 9: Development of the seepage area versus injection time in the layered jointed sample.

Figure 10 is the development of horizontal and vertical fracture radii versus injection time in the layered jointed sample. Both the horizontal radius and vertical radius increase with the injection time. It can be seen that the horizontal radius is larger than the vertical radius during the whole fracturing process, due to the low strength of the joint in the horizontal direction. The horizontal radius and vertical radius are 8⁻² m and 5.55⁻² m at time  s, respectively. Compared with the intact sample (horizontal radius and vertical radius are 6.3⁻² m and 6.49⁻² m, respectively, at time  s), the preexisting joints promote the expansion of cracks in the horizontal direction but restrain the crack expansion in the vertical direction.

Figure 10: Development of horizontal and vertical fracture radii versus injection time in the layered jointed sample.

4.3. Fracture Propagation in the Vertical Jointed Sample

The distribution and evolution of hydraulic fractures in the vertical jointed sample are presented in Figure 11. It can be seen that fractures are first generated around the borehole due to stress concentration. Then, fractures propagate along the diagonal direction and vertical direction with the increase of injection pressure, probably resulting from the equal load applied on the horizontal and vertical directions, and the preexisting joints in the vertical direction. Eventually, several main fractures are formed along the diagonal and vertical directions. Tensile damage is the main model of destruction, and the number of tensile cracks accounts for 93.9% of the total cracks.

Figure 11: Distribution and evolution of the fracture in the vertical jointed sample.

It can be seen from Figure 12 that the seepage area increases slowly from 3.96⁻⁶ m² to 2.18⁻³ m² with the injection time varying from  s to  s, owing to only a small amount of crack propagation at the low injection pressure stage. Then, the seepage area increases sharply from 2.18⁻³ m² to 8.95⁻³ m² with the injection time varying from  s to  s, due to the propagation of several main fractures at the same time in this period. Besides, the development of the seepage area is fitted by an exponential function, as shown in Figure 12.

Figure 12: Development of the seepage area versus injection time in the vertical jointed sample.

The development of horizontal and vertical fracture radii versus injection time in the vertical jointed sample is presented in Figure 13. It can be seen from the curve that fractures begin to propagate at time  s. With the increase of injection gas pressure, both the horizontal and vertical radii increase gradually. It should be noted that the vertical radius is greater than the horizontal radius in the whole injection process. The horizontal and vertical radii are 5.1⁻² m and 8⁻² m, respectively, at time  s. Compared with the intact sample (horizontal radius and vertical radius are 6.3⁻² m and 6.49⁻² m, respectively, at time  s), the preexisting vertical joints promote the expansion of cracks in the vertical direction but restrain the crack expansion in the horizontal direction.

Figure 13: Development of horizontal and vertical fracture radii versus injection time in the vertical jointed sample.

4.4. Fracture Propagation in the Orthogonal Jointed Sample

Orthogonal joints are also common in naturally fractured reservoirs. The distribution and evolution of hydraulic fractures after gas fracturing in the orthogonal jointed sample are shown in Figure 14. Fractures appear around the borehole at the initial stage (at time s), and it propagates to surrounding rock along the horizontal and vertical directions at time s. The preexisting orthogonal joints play a significant role in the propagation direction of the fracture: many generated fractures are extended in the direction of the preexisting natural fractures, as shown in Figure 14. During this period, fractures begin to expand gradually to form continuous cracks; the seepage area increases gradually from 1.98⁻⁵ m² to 3.13⁻³ m². The main fractures are connected to form the crush area near the injection hole at time s. As shown in Figure 15, the seepage area increases dramatically from 3.13⁻³ m² to 9.36⁻³ m², and a complex fracture network is developed eventually. The seepage area increases 35.05% compared with that of the intact sample at time s. Tensile damage is the primary damage mode in this situation with 93.02% of the damaged elements being in the tensile mode.

Figure 14: Distribution and evolution of the fracture in the orthogonal jointed sample.

Figure 15: Development of the seepage area versus injection time in the orthogonal jointed sample.

Figure 16 is the development of horizontal and vertical fracture radii with injection time in the orthogonal jointed sample. The fractures begin to propagate into the further surrounding rock at time  s. The horizontal radius and vertical radius are almost identical before time  s; then, the vertical radius becomes greater than the horizontal radius in the following process, which may result from the heterogeneous mechanical parameters used in this model. The horizontal radius and vertical radius are 7.3⁻² m and 8⁻² m at time  s, respectively. Compared with the intact sample (the horizontal radius and vertical radius are 6.3⁻² m and 6.49⁻² m, respectively, at time  s), the preexisting joints promote the expansion of fractures in the vertical direction.

Figure 16: Development of horizontal and vertical fracture radii versus injection time in the orthogonal jointed sample.

Based on the numerical simulations conducted on the intact sample, layered jointed sample, vertical jointed sample, and orthogonal jointed sample, the numerical results can be concluded in Table 2. The numerical results indicate that the preexisting horizontal joints or vertical joints promote the propagation of hydraulic fractures in the horizontal direction or vertical direction but restrain the expansion of hydraulic fractures in the vertical or horizontal direction. And the preexisting orthogonal joints promote the propagation of hydraulic fractures both in the horizontal direction and vertical direction. The similar results can also be observed in the research of Wang et al. [28] as shown in Figure 17. In addition, the preexisting joints can aggravate the damage of the numerical specimens; the seepage areas of the layered jointed sample, vertical jointed sample, and orthogonal jointed sample are increased by 32.5%, 29.16%, and 35.05%, respectively, at time  s compared with the intact sample.

Table 2: Numerical results of the intact sample, layered jointed sample, vertical jointed sample, and orthogonal jointed sample under CO₂ fracturing.

Sample
Damage area (10⁻² m²)
Horizontal radius (m)
Vertical radius (m)

Intact
0.693
0.063
0.0649

Layered jointed
0.918
0.08
0.0555

Vertical jointed
0.895
0.051
0.08

Orthogonal jointed
0.936
0.073
0.08

(a)

(b)

(c)

Figure 17: The distribution of hydraulic fractures of jointed samples: (a) layered jointed sample, (b) vertical jointed sample, and (c) orthogonal jointed sample [28].

5. Conclusions

In this work, a series of numerical simulations are performed on jointed samples under the CO₂-based hydraulic fracturing. A damage model is introduced to describe the initiation and propagation of microcrack in the fracturing process. The mechanical mechanism of the propagation of hydraulic fractures in several kinds of joint reservoirs is researched in this work. The numerical results of the intact sample, layered jointed sample, vertical jointed sample, and orthogonal jointed sample under CO₂ fracturing are shown in Table 2. Based on the results, the following conclusions are obtained.

It is shown that several continuous main fractures are formed uniformly in the intact sample. For the jointed samples, hydraulic fractures mainly propagate along the preexisting joints due to the lower strength and elastic modulus of joints compared with the rock matrix. Hydraulic fractures intersect with preexisting joints to form complex fracture networks. Besides, tensile damage is the main failure model in these numerical samples.

The variations of the seepage area versus injection time among these conditions are similar. The seepage area increases smoothly at the beginning since hydraulic fractures propagate slowly under low injection pressure, and then, it increases rapidly with the propagation of a large amount of hydraulic fractures to the surrounding rock. The whole process can be described by an exponential function. The seepage areas of the layered jointed sample, vertical jointed sample, and orthogonal jointed sample are increased by 32.5%, 29.16%, and 35.05%, respectively, at time  s compared with the intact sample. For the complexity of the fracture networks, it can be described by the development of horizontal and vertical fracture radii. The preexisting horizontal joints or vertical joints promote the propagation of the fracture in the horizontal direction or vertical direction but restrain the fracture expansion in the vertical or horizontal direction. In addition, the hydraulic fracture propagation is promoted along the preexisting orthogonal joints to form the complex fracture networks.

Data Availability

The numerical simulation data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (51804339), the China Postdoctoral Science Foundation (2018M640760), and the Innovation-Driven Project of Central South University. The corresponding author would like to thank the financial support by the State Key Laboratory for GeoMechanics and Deep Underground Engineering, China University of Mining and Technology (SKLGDUEK1805). We would like to thank ShiningStar Translation (email: shiningstartrans@foxmail.com) for providing linguistic assistance during the preparation of this manuscript.

References

C. T. Montgomery and M. B. Smith, “Hydraulic fracturing: history of an enduring technology,” Journal of Petroleum Technology, vol. 62, no. 12, pp. 26–40, 2010.

A. Reinicke, E. Rybacki, S. Stanchits, E. Huenges, and G. Dresen, “Hydraulic fracturing stimulation techniques and formation damage mechanisms—implications from laboratory testing of tight sandstone–proppant systems,” Chemie der Erde-Geochemistry, vol. 70, pp. 107–117, 2010.

V. Grechka, P. Mazumdar, and S. A. Shapiro, “Predicting permeability and gas production of hydraulically fractured tight sands from microseismic data,” Geophysics, vol. 75, no. 1, pp. B1–B10, 2010.

Y. Jia, Y. Lu, D. Elsworth, Y. Fang, and J. Tang, “Surface characteristics and permeability enhancement of shale fractures due to water and supercritical carbon dioxide fracturing,” Journal of Petroleum Science and Engineering, vol. 165, pp. 284–297, 2018.

D. Rahm, “Regulating hydraulic fracturing in shale gas plays: the case of Texas,” Energy Policy, vol. 39, no. 5, pp. 2974–2981, 2011.

D. Ma, M. Rezania, H. S. Yu, and H. B. Bai, “Variations of hydraulic properties of granular sandstones during water inrush: effect of small particle migration,” Engineering Geology, vol. 217, pp. 61–70, 2017.

J. P. Nicot and B. R. Scanlon, “Water use for shale-gas production in Texas, US,” Environmental Science & Technology, vol. 46, no. 6, pp. 3580–3586, 2012.

A. Josifovic, J. J. Roberts, J. Corney, B. Davies, and Z. K. Shipton, “Reducing the environmental impact of hydraulic fracturing through design optimisation of positive displacement pumps,” Energy, vol. 115, pp. 1216–1233, 2016.

I. Mrdjen and J. Lee, “High volume hydraulic fracturing operations: potential impacts on surface water and human health,” International Journal of Environmental Health Research, vol. 26, no. 4, pp. 361–380, 2016.

A. Vengosh, R. B. Jackson, N. Warner, T. H. Darrah, and A. Kondash, “A critical review of the risks to water resources from unconventional shale gas development and hydraulic fracturing in the United States,” Environmental Science & Technology, vol. 48, no. 15, pp. 8334–8348, 2014.

D. Ma, H. Duan, J. Liu, X. Li, and Z. Zhou, “The role of gangue on the mitigation of mining-induced hazards and environmental pollution: an experimental investigation,” Science of the Total Environment, vol. 664, pp. 436–448, 2019.

D. Ma, X. Cai, Q. Li, and H. Duan, “In-situ and numerical investigation of groundwater inrush hazard from grouted karst collapse pillar in longwall mining,” Water, vol. 10, no. 9, p. 1187, 2018.

A. Lisjak, B. S. A. Tatone, O. K. Mahabadi et al., “Hybrid finite-discrete element simulation of the EDZ formation and mechanical sealing process around a microtunnel in Opalinus Clay,” Rock Mechanics and Rock Engineering, vol. 49, no. 5, pp. 1849–1873, 2016.

Z. Zhou, X. Cai, D. Ma et al., “Water saturation effects on dynamic fracture behavior of sandstone,” International Journal of Rock Mechanics and Mining Sciences, vol. 114, pp. 46–61, 2019.

Z. Zhou, X. Cai, D. Ma, L. Chen, S. Wang, and L. Tan, “Dynamic tensile properties of sandstone subjected to wetting and drying cycles,” Construction and Building Materials, vol. 182, pp. 215–232, 2018.

J. D. Arthur, B. K. Bohm, B. J. Coughlin, M. A. Layne, and D. Cornue, “Evaluating the environmental implications of hydraulic fracturing in shale gas reservoirs,” in SPE Americas E&P environmental and safety conference, San Antonio, TX, USA, March 2009.

Q. Zhang, D. Ma, Y. Wu, and F. Meng, “Coupled thermal-gas-mechanical (TGM) model of tight sandstone gas wells,” Journal of Geophysics and Engineering, vol. 15, no. 4, pp. 1743–1752, 2018.

L. Wang, B. Yao, M. Cha et al., “Waterless fracturing technologies for unconventional reservoirs-opportunities for liquid nitrogen,” Journal of Natural Gas Science and Engineering, vol. 35, pp. 160–174, 2016.

J. B. Clark, “A hydraulic process for increasing the productivity of wells,” Journal of Petroleum Technology, vol. 1, no. 01, pp. 1–8, 1949.

C. F. Smith, “Gas well fracturing using gelled non-aqueous fluids,” in Fall Meeting of the Society of Petroleum Engineers of AIME, Las Vegas, NV, USA, September 1973.

S. Kothare, Economics and Applicability of Nitrogen for Fracking, Air Products and Chemicals, Inc, 2012.

A. B. Yost, “Analysis of production response to CO₂/sand fracturing: a case study,” in SPE Eastern Regional Meeting, Charleston, WV, USA, November 1994.

J. Zhou, G. Liu, Y. Jiang et al., “Supercritical carbon dioxide fracturing in shale and the coupled effects on the permeability of fractured shale: an experimental study,” Journal of Natural Gas Science and Engineering, vol. 36, pp. 369–377, 2016.

Y. Chen, Y. Nagaya, and T. Ishida, “Observations of fractures induced by hydraulic fracturing in anisotropic granite,” Rock Mechanics and Rock Engineering, vol. 48, no. 4, pp. 1455–1461, 2015.

X. Zhang, B. Wu, R. G. Jeffrey, L. D. Connell, and G. Zhang, “A pseudo-3D model for hydraulic fracture growth in a layered rock,” International Journal of Solids and Structures, vol. 115-116, pp. 208–223, 2017.

Q. Zhang, D. Ma, J. Liu, K. Zhang, and Z. Fan, “Numerical studies on the failure process of heterogeneous rock material with preexisting fracture under uniaxial compression,” Advances in Civil Engineering, vol. 2018, Article ID 9203549, 12 pages, 2018.

D. Ma, H. Duan, X. Cai, Z. Li, Q. Li, and Q. Zhang, “A global optimization-based method for the prediction of water inrush hazard from mining floor,” Water, vol. 10, no. 11, p. 1618, 2018.

T. Wang, W. Hu, D. Elsworth et al., “The effect of natural fractures on hydraulic fracturing propagation in coal seams,” Journal of Petroleum Science and Engineering, vol. 150, pp. 180–190, 2017.

J. Wang, D. Elsworth, and M. K. Denison, “Hydraulic fracturing with leakoff in a pressure-sensitive dual porosity medium,” International Journal of Rock Mechanics and Mining Sciences, vol. 107, pp. 55–68, 2018.

X. Zhang, B. Wu, L. D. Connell, Y. Han, and R. G. Jeffrey, “A model for hydraulic fracture growth across multiple elastic layers,” Journal of Petroleum Science and Engineering, vol. 167, pp. 918–928, 2018.

A. Riahi and B. Damjanac, “Numerical study of interaction between hydraulic fracture and discrete fracture network,” in Effective and Sustainable Hydraulic Fracturing, InTech, 2013.

P. Hou, F. Gao, Y. Ju et al., “Experimental investigation on the failure and acoustic emission characteristics of shale, sandstone and coal under gas fracturing,” Journal of Natural Gas Science and Engineering, vol. 35, pp. 211–223, 2016.

J. He, C. Lin, X. Li, and X. Wan, “Experimental investigation of crack extension patterns in hydraulic fracturing with shale, sandstone and granite cores,” Energies, vol. 9, no. 12, p. 1018, 2016.

W. C. Zhu and C. A. Tang, “Micromechanical model for simulating the fracture process of rock,” Rock Mechanics and Rock Engineering, vol. 37, no. 1, pp. 25–56, 2004.

C. H. Wei, W. C. Zhu, Q. L. Yu, T. Xu, and S. Jeon, “Numerical simulation of excavation damaged zone under coupled thermal–mechanical conditions with varying mechanical parameters,” International Journal of Rock Mechanics and Mining Sciences, vol. 75, pp. 169–181, 2015.

W. C. Zhu, D. Gai, C. H. Wei, and S. G. Li, “High-pressure air blasting experiments on concrete and implications for enhanced coal gas drainage,” Journal of Natural Gas Science and Engineering, vol. 36, pp. 1253–1263, 2016.

D. Ma, X. Miao, H. Bai et al., “Effect of mining on shear sidewall groundwater inrush hazard caused by seepage instability of the penetrated karst collapse pillar,” Natural Hazards, vol. 82, no. 1, pp. 73–93, 2016.

WebBook, NIST, “Thermophysical properties of fluid systems,” 2016, http://webbook.nist.gov/chemistry/fluid/.

J. Wang, D. Elsworth, Y. Wu, J. Liu, W. Zhu, and Y. Liu, “The influence of fracturing fluids on fracturing processes: a comparison between water, oil and SC-CO2,” Rock Mechanics and Rock Engineering, vol. 51, no. 1, pp. 299–313, 2018.

T. Xia, F. Zhou, J. Liu, J. Kang, and F. Gao, “A fully coupled hydro-thermo-mechanical model for the spontaneous combustion of underground coal seams,” Fuel, vol. 125, pp. 106–115, 2014.

Y. H. Huang, S. Q. Yang, W. L. Tian, J. Zhao, D. Ma, and C. S. Zhang, “Physical and mechanical behavior of granite containing pre-existing holes after high temperature treatment,” Archives of Civil And Mechanical Engineering, vol. 17, no. 4, pp. 912–925, 2017.

R. P. Covert and G. C. Philip, “An EOQ model for items with Weibull distribution deterioration,” AIIE Transactions, vol. 5, no. 4, pp. 323–326, 1973.

M. Hubbert and D. Willis, “Mechanics of hydraulic fracturing,” Transactions of Society of Petroleum Engineers of AIME, vol. 210, pp. 153–168, 1957.

D. Ma, Q. Li, M. R. Hall, and Y. Wu, “Experimental investigation of stress rate and grain size on gas seepage characteristics of granular coal,” Energies, vol. 10, no. 4, p. 527, 2017.