Research of fracture initiation and propagation in HDR fracturing under thermal stress from meso-damage perspective

Energy 178 (2019) 508e521

Contents lists available at ScienceDirect

Energy journal homepage: www.elsevier.com/locate/energy

Research of fracture initiation and propagation in HDR fracturing under thermal stress from meso-damage perspective Wei Zhang*, Tian-kui Guo**, Zhan-qing Qu, Zhiyuan Wang School of Petroleum Engineering, China University of Petroleum (East China), Shandong, Qingdao, 266580, China

a r t i c l e i n f o

a b s t r a c t

Article history: Received 12 September 2018 Received in revised form 19 April 2019 Accepted 20 April 2019 Available online 28 April 2019

By HDR (hot dry rock) fracturing the deep buried geothermal energy can be efficiently extracted from the established EGS (enhanced geothermal system). While the fracture initiation and propagation is subjected to the interaction of cryogenic induced thermal stress and liquid pressure. Based on the mesodamage mechanics, elastic thermodynamics and Biot seepage mechanics, a mesoscopic thermo-hydromechanical-damage coupling model (THM-damage) is proposed to analyze the fracturing stimulation in HDR. Firstly, the mathematical model and numerical implementation method is validated by hightemperature granite fracturing experiments. Secondly, the action mechanism of thermal stress in hydraulic fracturing of HDR is discussed. Thirdly, the evolution of multiple physical fields during the initiation and propagation of HDR fracturing is researched. Finally, the effects of various parameters on HDR fracturing process are also studied. The results indicate that when the rock temperature exceeds 200
C the fracture network can be formed by hydraulic fracturing, which extends along the direction perpendicular to the minimum in-situ stress. Increasing rock temperature can reduce fracture initiation pressure and rock failure pressure. The heat transfer coefficient between fracturing fluid and rock and the rock Young’s modulus have influence on the fracture morphology during HDR fracturing under thermal stress. © 2019 Elsevier Ltd. All rights reserved.

Keywords: HDR fracturing Enhanced geothermal system Thermal stress High-temperature fracturing test THM-Damage model

1. Introduction Geothermal energy has the advantages of green, stable, abundant and reproducible compared with fossil energy sources such as oil and gas. Although there is abundant geothermal energy in hot dry rock, the characteristics of large depth, low porosity and low permeability make it difficult to use the high temperature resources [1,2]. Formation of enhanced geothermal systems by artificial fracturing provides the possibility to efficiently extract geothermal energy from hot dry rock [3e5]. During the fracturing process, due to the injection of low temperature fluid (the ground temperature is about 20
C) into the high temperature rock (greater than 150
C), the thermal stress generated by the temperature gradient and the thermal stress generated by the non-uniform expansion of the particles have an important influence on the fracture initiation and propagation [6e9].

* Corresponding author. ** Corresponding author. E-mail addresses: [email protected] (W. Zhang), [email protected] (T.-k. Guo). https://doi.org/10.1016/j.energy.2019.04.131 0360-5442/© 2019 Elsevier Ltd. All rights reserved.

However, the fracture initiation and propagation mechanism in hot dry rock fracturing under the influence of thermal stress is still unclear, which makes it difficult to effectively guide the scheme design of hydraulic fracturing and accurate prediction of heat mining capacity. Some EGS field tests have proved that thermal cracking is the key to the building of EGS. Rutqvist et al. [9], Garcia et al. [10] adopted the micro-seismic monitoring in HDR mining field, the results showed that a large scale of high permeability zones were created by the injected cryogenic fluids. Fracturing research on oil and gas reservoirs started earlier, involving tight reservoirs [11,12], fractured reservoirs with natural fractures [13e15] and glutenite reservoirs [16]. However, since the oil and gas reservoir temperature is usually below 100
C, these studies have not considered the effect of temperature field. Therefore, the study of HDR fracturing affected by thermal stress is also a perfection of the hydraulic fracturing mechanism. In order to study the permeability enhancement mechanism of thermal cracking induced by low-temperature for high temperature tight reservoirs, many scholars only studied the thermal stress cracking caused by rapid cooling of hot rocks without considering

W. Zhang et al. / Energy 178 (2019) 508e521

the effect of injection pressure. Xi et al. launched the experimental research on mechanical properties of water-cooled granite. The results showed that the drastic changes of temperature caused thermal cracking and thermal shock in granite and the rock mechanical properties were degraded [17]. Kumari et al. investigated the mechanical behavior of Australian Strathbogie granite under the increasing temperature (from room temperature to 800
C) followed by two cooling methods (both rapid and slow), besides, the stress concentration zone was confirmed by ARAMIS images [18]. The conclusion similar to Xi was obtained. Kumari et al. heated the granite to 600
C then cooled it by cold water, the rock permeability increased from 1  1019 m2 to 6.5  1015 m2 [19]. From the above experiments, it can be seen that when the injection pressure is not applied, it is required at least a temperature difference of 500e600
C to obtain an effective permeability increase only rely on the thermal cracking caused by temperature gradient [17e19]. Therefore, for high temperature rock of 150e400
C, it’s necessary and efficient to achieve the fractures or fracture networks in HDR by hydraulic fracturing under the combined action of injection pressure and thermal stress [20e22]. Guo proved that hydraulic fracturing at 100
C can form multiple fractures through the prefabricated cement blocks [23]. Zhou et al. [24] and Kumari et al. [25] conducted hydraulic fracturing experiments with cylindrical rock samples under different temperatures and confining pressures. The results showed that the higher the rock temperature can reduce the failure pressure of rock. In addition, the stimulation mechanism of Liquid nitrogen (LN) fracturing is inseparable from low-temperature induced thermal stress. Cha et al. carried out Liquid nitrogen fracturing test on concrete blocks and sandstones, the results showed that the formation of fractures is related to rock properties [26]. Zhang et al. analyzed the fracture morphology formed by liquid nitrogen fracturing in granite with different rock temperature [27]. It is showed that the increase of rock temperature promotes the connectivity of fracture network. However, the triaxial stress was overlooked in the study of Cha and Zhang. Through the high temperature fracturing experiments, we obtain the macroscopic understanding that the increase of rock temperature can reduce the rock failure pressure and form multiple fractures [24e27]. However, laboratory tests cannot describe the process of fracture initiation and propagation during rock failure. In particular, we want to know how the fracture networks are initiated and interconnected during the HDR fracturing, as shown in Fig. 1. As we know, the initiation, propagation and the interactions of multi-cracks in heterogeneous materials is a complex process [29,30]. Scholars have carried out numerical simulations of rock instability and fracture propagation based on the meso-damage mechanics. Thus, many engineering problems are effectively analyzed, such as the hydraulic fracturing in natural fractured shale reservoir [31], glutenite fracturing [32], underground radioactive waste repository [33], deep underground excavation [34] and repeat fracturing [35]. The above studies proved that the thermal stress and thermal cracking caused both by thermal loading and water cooling have key effects on mechanical properties and seepage characteristics of rocks in laboratory-scale. Because of the low-permeability of HDR matrix, the stimulated fractures become the main heat transfer channels for heat carrying medium, therefore, it’s of significance to confirm the fracture evolution rules in HDR fracturing under thermal stress for the accurate prediction of heat mining performance [36e38]. Based on the meso-damage mechanics, elastic thermodynamics and Biot seepage mechanics, a mesoscopic thermo-hydromechanical-damage coupling model (THM-damage) is proposed

509

Fig. 1. Schematic diagram of the fracture morphology after HDR hydraulic fracturing.

and implemented to analyze the fracturing stimulation in hot dry rocks. Firstly, the mathematical model and numerical implementation method using COMSOL combined with MATLAB programming are validated by the high-temperature granite fracturing experiments under the triaxial stress. Secondly, the action mechanism of thermal stress in hydraulic fracturing of hot dry rock is discussed. Thirdly, the evolution of multiple physical fields during the initiation and propagation of HDR fracturing are researched. Finally, the effects of various parameters on the stimulation process of HDR fracturing are also studied. This work reveals the formation of fracture network in HDR fracturing under thermal stress, and provides an effective simulation method for fracture morphology prediction. In addition, more accurate simulation of heat mining performance can be carried out according to the reasonable fracture morphology. 2. Numerical model HDR fracturing stimulation is a multi-physical coupling process of matrix deformation, fracturing liquid flow and heat transfer between the injected liquid and hot matrix. 2.1. Governing equations of thermo-hydro-mechanical coupling model The static equilibrium equation of a linear elastomer when considering the pore fluid pressure and temperature changes [28,29]:

Gui;jj þ

G u  ap;i  K ’ aT T;i þ Fi ¼ 0 1  2n j;ji

(1)

Where, G is the shear modulus, (Pa), G ¼ E=2ð1 þ vÞ, v is the Poisson’s ratio, E is the elastic modulus (Pa); Fi and ui ði ¼ x; y; zÞ are respectively the components of body force and displacement in i 2 2 2 v2 u þ v2 v þ v2 w ;ap is the direction, ui;jj ¼ vvxu2 þ vvyu2 þ vvzu2 , uj;ji ¼ vxvx ;i vxvy vxvz water pressure action item (penetration force under pore pressure), a is Biot coefficient; K 0 aT T;i is the thermal stress action item, aT is the matrix thermal expansion coefficient (
C1), K 0 is the drainage

510

W. Zhang et al. / Energy 178 (2019) 508e521

volume modulus (Pa),K 0 ¼ 2Gð1 þ vÞ=3ð1  2vÞ。 The seepage field governing equation considering the action of mechanical stress and temperature [31,32]:

c1


 vεv vT vp k  c2 þ c3 ¼ V, ðVp þ rl gVzÞ vt vt vt m

(2)

8 K0 > > ¼ 1  c > 1 > > Ks > > > <

Where, εv is volumetric strain; f is medium porosity; al is the volume thermal expansion coefficient of fluid (
C1); as is the volume thermal expansion coefficient of solid (
C1); bl is the volume modulus of porous fluid (Pa); Ks is the effective volume modulus of solid (Pa); k is permeability of continuous medium (m2); m is dynamic viscosity of fluid (Pa$s); rl is fluid density (Kg/m3); g is gravity acceleration (m/s2); z is the coordinates in the vertical direction (m). The temperature field governing equation considering thermal convection and mechanical stress [33]:

vT vε k þ ðT0 þ TÞK 0 aT T þ rl Cl ðT0 þ TÞ Vp ¼ lM V2 T vt vt m

(3)

Where, lM ¼ ls ð1  fÞ þ ll f; ls and ll respectively represents the thermal conductivity coefficient of matrix and fluid (W/(m$K)); T0 is the reference temperature under zero stress state (K); the heat capacity of porous media containing fluid can be described as ðrCÞM ¼ rs Cs ð1  fÞ þ rl Cl f (KJ$m3$K1). Equations (1)e(3) consider the heat and mass transfer in thermodynamics and the contraction of each component under the action of mechanics stress and temperature, then constitute a coupled nonlinear equations for controlling the thermo-elastic response of saturated porous media.

(4) 1 þ sin 4 s 1  sin 4 1

(5)

Where, ft is the rock uniaxial tensile strength, Pa; fc is the rock uniaxial compressive strength (Pa); 4 is the rock internal friction angle (
). Due to the HDR reservoirs are generally granite, the range of value is 45e60
. The elastic-brittle damage judgment model is adopted. When it is satisfied F1  0, the rock is subjected to tensile damage, and the relationship between the damage variable D and the strain can be expressed as [31]:

8 0; ε < εt0 > > > < lε D ¼ 1  t0 ; εt0  ε  εtu > ε > > : 1; εtu  ε

According to the elastic damage theory, the elastic modulus is given by the following relationship [32,33]:

E ¼ ð1  DÞE0

  f ¼ ðf0  fr Þexp  af sv þ fr

(9)

Where, f0 is porosity at zero stress state; af is the stress sensitive coefficient of porosity, 5.0  108 Pa1; fr is the limit value of porosity at high compressive stress state; sv is average effective stress, and can be calculated as follows [28]: a is Biot coefficient. the s1 ,s2 and s3 represses the principal stresses respectively.

sv ¼ ðs1 þ s2 þ s3 Þ=3  ap

(10)

The effect of damage on permeability is as follows:

f f0

3

expðaD DÞ

(11)

where, k0 is permeability under zero stress state (m2); k is the permeability under stress action (m2). The effect of damage on thermal conductivity can be expressed as follows:



ls ðT; DÞ ¼ ls ðTÞexp D=a T

 (12)

2.4. Rock parameter assignment based on statistical distribution To describe the heterogeneity of the material properties, it is considered that the mesoscopic element properties of the rock satisfy the Weibull distribution, which can be defined by the following function [35,39].

f ðuÞ ¼ (6)

(8)

where, E0 and E are the elastic modulus before and after the damage, respectively. Equation (1) embodies the influence of thermal stress caused by temperature difference and fluid pressure on rock deformation process. In addition, the rock porosity is also related to the stress state the rock subjected, as shown in following:

k ¼ k0

When the rock stress state satisfies the maximum tensile stress criterion F1  0 or the Mohr-Coulomb criterion F2  0, the tensile or shear damage will occur, the expression of F1 and F2 can be described as [29]:

F2 ¼ s3  fc þ

(7)

2.3. Effect of damage on solid deformation, fluid flow and heat transfer performance



2.2. The failure criterion

F1 ¼ s1  ft

8 < 0; εc0 < ε D¼ : 1  lεc0 ; ε  ε t0 ε

Where, l is the residual tensile strength coefficient, l ¼ ftr =ft0 , ft0 and ftr expresses the uniaxial tensile strength and residual strength respectively, εt0 is the tensile strain corresponding to the elastic limit, the ultimate tensile strain can be described as εtu ¼ h$εt0 , h is the ultimate strain coefficient.

a K0 c2 ¼ fal þ ð1  fÞas  T > Ks > > > > > f 1  f > > : c3 ¼ þ Ks bl

ðrCÞM

When it is satisfied F2  0, the rock is subjected to shear damage, and the relationship between the damage variable D and the strain can be expressed as:

 
 m  m u m1 u exp  u0 u0 u0

(13)

Where, u represents the value that satisfies the distribution parameters (such as strength, elastic modulus, Poisson’s ratio, etc.); u0 is a parameter related to the average of all element parameters.

W. Zhang et al. / Energy 178 (2019) 508e521

The shape parameter m defines the shape of the Weibull distribution density function. 2.5. Numerical implementation In the thermo-hydro-mechanical coupling, the physical equations are highly nonlinear partial differential equations in space and time domain. A series of solvers provided by COMSOL can achieve efficient numerical solution of the coupled nonlinear equations. The real-time connection between the MATLAB programming language and COMSOL is realized by COMSOL with MATLAB port, making the multi-physical equations and transient solver easily invoked by the designed MATLAB program, thus providing convenience for the establishment of complex coupling equations and the implementation of finite element [40]. In this study, the finite element solver COMSOL is adopted to solve the partial differential equations (PDE) of the THM coupling to analyze the interaction between the increasing liquid pressure and thermal stress during the cryogenic fluid injection process. And COMSOL combined with MATLAB programing is employed to numerically realize the THM-damage process of HDR fracturing stimulation. 3. Model validation The self-designed high-temperature fracturing device is used to verify the established THM-damage coupling and the numerical implementation method of COMSOL combined with MATLAB programming. As shown in Fig. 2(e), the experimental device can load axial pressure in the direction of X, Y and Z to simulate different insitu stress values. The granite samples used in the experiments are 0.1 m  0.1 m  0.12 m. Firstly, the borehole with bottom diameter of 0.008 m and upper diameter of 0.014 m was drilled on the rock sample, as shown in Fig. 2(a), and then the wellbore with perforations was put into the drilled borehole. The annulus between wellbore and rock sample is sealed by high temperature sealant. The specific installation method of wellbore is shown in Fig. 2. In order to ensure that the wellbore installed by high temperature sealant can withstand as high pressure as possible, two metal

511

sealing rings are installed on each wellbore. (Fig. 2(a) (d)). The rock samples used in the experiment are Rizhao granite from Shandong Province. Before the high temperature hydraulic fracturing experiment was launched, rock cores with diameter of 2.5 cm was drilled in granite to carry out the rock physical and mechanical properties tests, including permeability, porosity, Young’s modulus, tensile strength, compressive strength and so on. Then the key parameters obtained from granite will be applied to the validation numerical simulations, as shown in Table 1. Fig. 3(a) shows the boundary conditions in high-temperature fracturing experiments. Concerning the triaxial condition, we adopt thesV ¼ 10 MPa, sH ¼ 7.5 MPa, sh ¼ 5 MPa. Roller boundary is placed on the opposite of the surface subjected to the in-situ stress, and the cryogenic fracturing liquid is continuously injected from the wellbore. The rock temperature of 200
C, 100
C and 50
C were adopted in the fracturing validation experiments. In high temperature experiments, the specimen under in-situ stress first underwent a slow thermal loading process (5
C/min) in order to reduce the damage caused by thermal loading. Then keep them at the designed temperature (200
C, 100
C and 50
C) for 12 h to make the rock evenly heated. Distilled water was used as fracturing medium in the experiments. In addition, after the fracturing experiment the red ink was injected into the fractured rock samples to observe the flow path of the red ink along fractures. Limited by the computer’s ability the two-dimension validation model is used instead of the 3D model, which is intercepted from the three-dimensional model in Fig. 3(b), so the size of the validation numerical model is 0.1 m  0.1 m. Fig. 3(c) shows the boundary conditions in 2D numerical simulations. The 2D model is divided into 100,000 computing units using the free quadrilateral grid. The injection water temperature is set to 20
C and the heat

Table 1 Parameters of the verification model. E/(GPa) 32 st/(MPa) 20.5

km/(m2) 18

0.6  10 sc/(MPa) 190

Fig. 2. Wellbore installation sketch and high-temperature hydraulic fracturing device.

Tm/(
C)

n

m [29]

200/100 aT/(1/K) 6.2  106

0.2 Tinj/(
C) 20

10 mf/(cP) 1

512

W. Zhang et al. / Energy 178 (2019) 508e521

Fig. 3. Schematic diagram of boundary condition settings in high-temperature fracturing experiments; (a) and the 2D Section employing in numerical simulation.

transfer coefficient between water and granite is set to 2000 W/ (m2$K). E elastic modulus, n Passion’s ratio, st tensile strength, sc compressive uniaxial strength, Tm matrix initial temperature, Tinj injection fluid temperature, km matrix permeability, aT thermal expansion coefficient, mf fracturing fluid viscosity, m homogeneity index [29]. After the high-temperature fracturing, the red ink flowing out of the rock can intuitively depict the fracture formed under the combination of thermal stress and injection pressure. Fig. 4(a) and Fig. 5(a) show the fracture morphology of the 100
C/200
C granite after hydraulic fracturing from two perspectives. Figs. 4(b) and Fig5(b) show the outflow of red ink on ①-④ surface. It can be seen that the formed fractures extend perpendicular to the direction of minimum in-situ stress. Besides, when the rock temperature rises, the thermal stress generated by the injection of fracturing fluid is more severe, and the influence of thermal stress on the fracture morphology is greater. As shown in Fig. 5(b), when the granite temperature is 200
C, multiple fractures appear on the surface①

Fig. 4. Fracture morphology obtained by experiments at 100
C; (a) two views of the fractured rock; (b) fracture distribution on surface ①-④

Fig. 5. Fracture morphology obtained by experiments at 200
C; (a) two views of the fractured rock; (b) fracture distribution on surface ①-④

and ②, and the corresponding red ink flow out simultaneously from multiple points on surface① and ②. In the validation numerical simulation, Fig. 6 and Fig. 7 show the damage evolution of the fracturing rock on the two-dimension section (as shown in Fig. 3(b)) at different temperatures. When the rock temperature is 200
C, more fractures initiate along the direction perpendicular to the minimum in-situ stress under the action of low-temperature induced thermal stress. With the increase of injection pressure, fractures continue to extend accompanied by the formation of branching fractures, and then interconnect with each other, finally forming the fracture network of a certain scale. For the rock at 100
C, many fractures are formed along the direction perpendicular to the minimum in-situ stress, but the fractures are not connected. Compared with rock fracturing at 200
C, the number of damaged fractures in the rock fracturing of 100
C decreases obviously. In addition, as the rock temperature increases, the injection pressure required for the initiation and propagation of fractures decreases under the influence of thermal stress. Under the low-temperature induced thermal stress, the

W. Zhang et al. / Energy 178 (2019) 508e521

513

In addition, the validation simulation also shows that the damage process inside the rock can be visually depicted by using the meso-damage mechanics analysis, which provides a powerful tool for studying the rock damage and failure process under the coupling of multi-physical fields.

(a) 14MPa-11

(b) 16MPa-82

4. Mechanism of thermal stress in hydraulic fracturing of hot rock

(c) 18MPa-17

100
C

Fig. 6. Damage evolution of hot rock fracturing at on the 2D Section in validation tests (in Fig 3(b)). (In the illustration, 0 to 1 represents the damage degree, 1 represents the complete damage).

(a) 8MPa-19

(b) 13MPa-22

(c) 15MPa-39

Fig. 7. Damage evolution of hot rock fracturing at 200
C on the 2D Section in validation tests (in Fig 3(b)).

fracture initiation pressure is reduced to 5 MPa when the temperature of granite fracturing is 200
C, while at 100
C, the initiation pressure is 8 MPa. In order to quantitatively analyze the damage evolution in the fracturing numerical simulation, the number of damaged elements in each pressure loading step is recorded in the process of damage calculation. Fig. 8(a) shows the variation of the injection pressure in the fracturing experiment, and Fig. 8(b) shows the variation of the damaged elements with the injection pressure in numerical simulation. It can be seen that the numerical simulation is consistent with the rock failure pressure at different temperatures obtained by laboratory experiments. The corresponding failure pressures at the rock temperatures of 50
C, 100
C and 200
C is respectively 24 MPa, 21 MPa and 16 MPa. As shown in Fig. 8(a), when the temperature is 200
C, the fracture is more likely to initiate and extend under thermal stress, so the slope of rising portion in the pressure curve is lowered. Compared with the fracturing curve, the description of the damaged elements varies with the injection pressure can be more convenient to analyze the initiation and propagation pressure of the rock.

In order to clarify the role of thermal stress in hydraulic fracturing, without considering the impact of damage, the following studies have been carried out by means of the verified coupling model: ① firstly, we analyzed the stress distribution around the wellbore when only the thermal stress exists, that is, the injection pressure in the fracturing process was not considered, but only the influence of thermal stress generated by low-temperature fluid on stress field of hot rock was analyzed; ② after defining the influence of low-temperature induced thermal stress on stress field, to simulate the hot rock fracturing process the increasing injection pressure was added to the wellbore under the condition of thermal stress. Then the variation of the maximum tensile stress in hot rock under the combined action of thermal stress and injection pressure was analyzed. Furthermore, the mechanism of the interaction between thermal stress and injection pressure on rock initiation is also studied. 4.1. Computational conditions and parameters The numerical model is assembled with a 0.3 m  0.3 m 2D square, and a 0.015 m diameter borehole in center. Fig. 9 describes the boundary condition settings. The typical parameters for THM coupling are shown in Table 2, and the homogeneous model is used in this part. Because the temperature field in rock varies with time, the model is solved by a transient solver in COMSOL [40], and time step of numerical calculation is 1s. The injection temperature is set to 20
C and the heat transfer coefficient between water and granite is set to 2000 W/(m2$K). The research is conducted in a 200
C hot rock with 20
C injection liquid, and the increasing liquid pressure is respectively set as 0/4/8/12/16/20 MPa. Meanwhile, the interaction mechanism between liquid pressure and thermal stress as well as its influence on fracture initiation is analyzed. 4.2. Stress distribution only under the cryogenic induced thermal stress The rapid cooling of low-temperature fracturing fluid to hot rock









7L  7L  7L 



'DPDJHG (OHPHQWV

3UHVVXUH 03D



7L  7L  7L 



    

 















7LPH V

(a) Injection pressure curve in fracturing tests























,QMHFWLRQ 3UHVVXUH 03D

(b) damage curve in numerical simulation

Fig. 8. Fracturing curves of granite at different temperatures; (a) Injection pressure curve in fracturing tests; (b) damage curve in numerical simulation.

514

W. Zhang et al. / Energy 178 (2019) 508e521



7HPSHUDWXUH 



    

7LPH 7LPH 7LPH 7LPH 7LPH 7LPH

  

V V V V V V

 







'LVWDQFH P

Fig. 9. Model schematic.

(a) Evolution of temperature distribution Table 2 Typical parameters for THM coupling [33,36,41]. km/(m )

lm/(W/m/K)

E/(GPa)

n

2700 sH/(MPa) 15

1  1018 sh/(MPa) 10

4

37 4 0.01

0.2 Cm/(J/kg/K) 950

2

aT/(1/K)

6  106

E elastic modulus, n Passion’s ratio, sY In-situ stress in Y, sX In-situ stress in X, km matrix permeability, rm matrix density, aT thermal expansion coefficient, lm thermal conductivity, 4 porosity, Cm heat capacity.

results in the formation of temperature gradient around the wellbore, which leads to the generation of thermal stress [17,18]. However, there is a cooling process in the rock matrix surrounding the wellbore by the cooling of the cryogenic fluid. As shown in Fig. 10(a), when the rock temperature is 200
C, the matrix temperature at the wellbore edge drops from an initial temperature of 200
Ce120
C at 1 s and to 73
C at 10 s. With the passage of time, the low temperature gradually transfers from the wellbore to the inside of the hot rock. At 100s, the matrix temperature is still 200
C at 0.05 m away from wellbore edge, and at 1000s, the temperature drops to 165
C, then at 1000s, the temperature at this place falls to 87
C. It can be seen that the process of low temperature transfer from wellbore to rock sample is accompanied by the gradual decrease of temperature gradient. When the rock temperature is 200
C, the evolution of the maximum principal stress along the transversal (as shown in Fig. 9) with time is shown in Fig. 10(b), the tensile stress generated by rapid cooling at t ¼ 1s is 5 MPa. With the increase of cooling time, temperature at the wellbore edge gradually decreases while the thermal stress increases gradually. The tensile stress is 18.6 MPa at t ¼ 10s and 26.1 MPa at t ¼ 100s. However, due to the constant heat transfer between cryogenic fluid and hot rock matrix, the final temperature of rock sample and injection fluid will tend to balance. At this time, there is no temperature gradient, and the thermal stress decreases gradually until it disappears. The distribution of maximum principal stress at different rock temperatures for transient solution t ¼ 1s is shown in Fig. 11. As the rock temperature increases, the tensile stress generated by the cold shrinkage increases. The maximum tensile stress around the wellbore is 1.07 MPa when the rock temperature is 100
C, and 5.25 MPa when the temperature increases to 200
C. Therefore, during the rapid cooling of hot rock, the thermal stress presents as the tensile stress, which increases first and then

0D[LPXP 3ULQFLSOH 6WUHVV 03D



rm/(kg/m ) 3

7LPH 7LPH 7LPH 7LPH 7LPH 7LPH

  



 



V V V V V V

 







 











  









'LVWDQFH P

(b) Evolution of stress distribution Fig. 10. Evolution of temperature and stress distribution along the two-dimensional section in rock sample.

decreases with time. The first increase of tensile stress is due to the continuous cooling around wellbore edge. The reduction phase of tensile stress is owing to the heat exchange between the injection fluid and the rock, then the temperature tends to balance and the temperature gradient gradually disappears. 4.3. Stress distribution under the combination of thermal stress and liquid pressure The hot rock hydraulic fracturing process involves the interaction of injection pressure and low-temperature induced thermal stress. On the basis of defining the law of thermal stress acting on the wellbore surrounding (4.1), influence of the increasing injection pressure on stress field of rock under thermal stress is studied, and the interaction between injection pressure and thermal stress during the fracturing process and its influence on fracture initiation

W. Zhang et al. / Energy 178 (2019) 508e521

515

Fig. 11. Maximum principal stress distribution under different matrix temperatures with only thermal stress.

are analyzed. Fig. 12 shows the distribution of stress field in rock with the increase of injection pressure when the initial temperature of matrix is 150
C and the distilled water is injected at 20
C. When the injection pressure is not loaded, the maximum tensile stress around the wellbore caused by thermal stress is only 1.4 MPa. When the injection pressure increases to 4 MPa, the maximum tensile stress around the wellbore is 5.15 MPa. When the injection pressure increases to 12 MPa, the maximum tensile stress is 15.3 MPa. It can be seen that under the influence of low-temperature induced thermal stress when injecting cryogenic fluid into the

hot rock, the tensile stress of rock is caused by the combined contribution of low-temperature induced thermal stress and injection pressure, and the maximum tensile stress increases with the raise of injection pressure [27]. This is because the thermal stress generated by the cold shrinkage of the rock is presented in tensile stress around the wellbore and the injection pressure also acts in tension stress. The tensile stress caused by the two reasons (injection pressure and thermal stress) is superimposed on each other, so the hot rock is more likely to reach the tensile failure condition. The increase of rock temperature intensifies the lowtemperature induced thermal stress during fracturing. Fig. 13 shows the maximum tensile stress of rock at t ¼ 1s under the

Fig. 12. Effect of the increasing liquid pressure on the distribution of maximum principal stress at a rock temperature of 150
C.

0D[LPXP 3ULQFLSOH 6WUHVV 03D

516

W. Zhang et al. / Energy 178 (2019) 508e521

4.5. Evolution of multiple physical fields

                 

5RRP 7HPSHUDWXUH   





In HDR hydraulic fracturing, the temperature field has an important influence on the initiation and propagation of fractures. Therefore, the evolution of damage field, temperature field, seepage field and stress field in the fracturing process is analyzed for rocks at 100
C and 200
C.







,QMHFWLRQ 3UHVVXUH 03D Fig. 13. Maximum tensile stress under different matrix temperatures with the combination of thermal stress and liquid pressure.

increasing injection pressure. For example, when the injection pressure is 4 MPa, at room temperature the rock is still in the compressive state under the in-situ stress, when at 150
C the rock is subjected to tensile stress of 5.1 MPa along the direction of the maximum in-situ stress, while the rock at 200
C is subjected to the tensile stress of 10.6 MPa. It can be seen that the rock with higher temperature is more likely to reach the tensile failure criterion and begin to damage. The results are consistent with those obtained by Rutqvist [9] and Garcia [10] in field construction experiments. 4.4. Fracturing stimulation process of HDR The coupling of injection pressure and thermal stress during the injection of cryogenic fluid and its effect on fracture initiation are analyzed above. In this part, the established THM-damage model is used to investigate the process of fracture initiation and propagation in hydraulic fracturing of hot rocks under the influence of thermal stress. Because the exact three-dimensional simulation research needs to reach the numerical solving scale of more than ten million freedom degrees [29], limited to the computing power, the model used in this part is a 0.3 m  0.3 m two-dimension square tight rock with a 0.015 m diameter drilling hole in central, as shown in Fig. 9. The maximum and minimum in-situ stresses of the rock are respectively 15 MPa and 10 MPa. The tight rock model is divided by free triangular grid with 26 thousand calculation units with 100,000 computing nodes, and the maximum unit size is 0.003 m. The heterogeneity of rock parameter m ¼ 10, the injection pressure increases by 0.5Mpa per load step. Fig. 9 shows the boundary condition schematic. Table 2 provides the basic model parameters for THM coupling. The parameters related to the simulation of rock damage are illustrated in Table 3.

Table 3 Parameters related to damage model [29,33].

st//(MPa)

sc/(MPa)

hc

ht

l

q/(
)

16

170

150

5

0.1

45

hc ultimate compressive strain coefficient, ht ultimate tensile strain coefficient, lresidual strength factor, q internal friction angle.

4.5.1. Evolution of damage field When the rock temperature is 100
C, only the low-temperature induced thermal stress is insufficient to cause the thermal cracking. Until the injection pressure increases to 10 MPa, the rock damage begins to appear under the combination of thermal stress and injection pressure. The emerging fracture extends perpendicular to the direction of minimum in-situ stress. And the propagation of fractures requires an increasing injection pressure. As shown in Fig. 14, the fractures formed at 20 MPa include a two-wing fracture extending outward from the wellbore and a fracture extending only in the Y direction, but the single fracture along Y direction expands slower than the two-wing fracture. Therefore, under the given simulation conditions, three fractures are formed after rock fracturing at 100
C, but the fractures are not interconnected. When the rock temperature increases to 200
C, the lowtemperature induced thermal stress intensifies, and the rock begins to damage when the injection pressure is 5 MPa. Comparing Fig. 14(a) and (b), it can be seen that more fractures initiate due to the influence of thermal stress. With the increase of injection pressure, more branching fractures continue to form, then the branching fractures cross-propagate in the rock and gradually connect with each other. At the injection pressure of 12 MPa, the fracture network begins to form, and it extends perpendicular to the minimum in-situ stress under the in-situ stress. Before the rock is completely destroyed, the increase in injection pressure will enhance the complexity of the fracture network. For the construction of EGS, the increase of fractures can provide more flow channels for heat-carrying fluids. What’s more, the complex network can make the heat transfer of heat-carrying fluid more efficient. 4.5.2. Evolution of temperature field, seepage field and stress field As the cryogenic fluid injected into the hot rock at beginning of the fracturing, temperature and seepage fields are disturbed, resulting in temperature and pressure gradients. The thermal stress is generated by temperature gradient and the pore pressure is formed by pressure gradient. At the same time, the stress field of rock is disturbed by the combined action of thermal stress and injection pressure. When the failure condition is reached, the rock begins to damage [42e44]. Then fractures or micro-fractures are formed after damage, meanwhile the mechanical properties of the damage elements are reduced and the heat and mass transfer abilities are improved, as illustrated in 2.3. After the fractures emerge, the cryogenic fluid is present in the fractures, so the temperature inside the fractures is lower. Over time, although cryogenic fluid is continuously injected into the well, the cryogenic fluid in the fractures will be heated by the surrounding hot rock. As shown in Fig. 15(a), when the rock temperature is 100
C, the temperature of the fracture tip along -y direction is 54 C in the calculation step of 19MPa-3, and 75
C in the calculation step of 20 MPae34. As the rock temperature increases to 200
C, more fractures are created around the wellbore during fracturing. The cryogenic fluid is heated by surrounding rock in fractures. Due to the increase of rock temperature, the fluid temperature in fractures rises faster. At 10 MPa-13 calculation step, the average temperature of fluid in multiple fractures along -y direction is 79
C, and at 14 MPa-24

W. Zhang et al. / Energy 178 (2019) 508e521

517

(a) rock temperature of 100

(b) rock temperature of 200 Fig. 14. Evolution of damage field for the rock temperature of 100
C and 200
C (The legend represents damage index, 0  D  1 for tensile damage).

calculation step, the average temperature of fluid in multiple fractures along -y direction is 112
C. After the fractures are formed, the injected fluid enters the fractures and produces the pressure gradient in rock [44,45]. As shown in Fig. 15, the injection pressure quickly propagates along the fracture and acts on the matrix around the fracture. Compared with the evolution of temperature field, injection pressure propagates faster and there is no loss in the propagation process, which is due to the heating of fluid by hot rock matrix. Since a large amount of interconnected branch fractures are formed at rock temperature of 200
C, the propagation range of injection pressure is wider. In addition, since the numerical model is based on the physical parameters of granite, the permeability of matrix is extremely low. So after the injection pressure propagates in the fractures there is almost no diffusion of the injection pressure into the matrix. The evolution of the stress field in the fracturing process determines the rock damage and failure. When the stress of the calculation element reaches the failure condition, the element will be damaged. As shown in Fig. 15(a), the maximum tensile stress is increased by the injection pressure. When rock temperature is 100
C, the maximum tensile stress is 21.3 MPa under the injection pressure of 19 MPa, and the maximum tensile stress increases to 26 MPa under the injection pressure of 20 MPa. For the rock temperature of 200
C in Fig. 15(b), due to the more intense low-temperature induced thermal stress, the maximum tensile stress is 20 MPa under the injection pressure of 10 MPa, and the maximum tensile stress increases to 26 MPa under the injection pressure of 14 MPa. When the heterogeneity is m ¼ 10, it can be seen from the analysis in 4.1 that both the low-temperature induced thermal stress and the injection pressure cause damage to the rock by tensile stress. Comparing Fig. 15(a) and (b), it can be seen that the increase of rock temperature enhances the thermal stress induced by low-temperature, and more tensile stress concentration occurs around the wellbore. Simultaneously, more calculation elements will reach the failure criterion and form fractures. In addition, the compressive stress caused by thermal expansion of rock particles

increases with the raise of temperature. 4.6. Failure mechanism in HDR fracturing By analyzing the influence of various parameters on fracture morphology, we try to clarify the failure mechanism of HDR fracturing process [46,47]. The damage scale of rock is quantitatively analyzed by recording the change of damage elements with the injection pressure during the fracturing process. Variation of damage elements with injection pressure during fracturing under different conditions is shown in Fig. 16. 4.6.1. Rock temperature Fracture initiation pressure and rock failure pressure decrease with the increase of rock temperature. As shown in Fig. 16(a), the fracture initiation pressure is 16 MPa when the temperature is 50
C, and 10 MPa when the temperature is 100
C. During the fracture propagation stage, the velocity of damage increases as the temperature raises. The fracture morphology of rock fracturing under different temperature is shown in Fig. 14. 4.6.2. Young’s modulus It can be known from equations (3)e(1) that the bulk modulus of the porous medium in the thermal stress term is proportional to the Young’s modulus of the rock [33]. Therefore, the increase of Young’s modulus causes the greater thermal stress under the same temperature gradient (Fig. 16(b)). As shown in Fig. 17, the increase in Young’s modulus makes the fracture network more complex. It can be seen that in the layer selection of EGS, the selection of hot rock layers with larger Young’s modulus is conducive to constructing a more complex fracture network in EGS. 4.6.3. In-situ stress condition In this part, l represents the differential coefficient of in-situ stress. As shown in Fig. 16(c), with the increase of in-situ stress, the more serious the rock is subjected to compression, and the higher injection pressure is required for the tensile failure of rock.

518

W. Zhang et al. / Energy 178 (2019) 508e521

(a) rock temperature of 100

(b) rock temperature of 200 Fig. 15. Evolution of temperature field, seepage field and stress field during the fracturing for the rock temperature of 100
C and 200
C.

For high temperature fracturing, the greater the difference of insitu stress, the rock will be more likely to crack. However, the large difference of in-situ stress makes the fracture network expand rapidly along the direction perpendicular to the minimum in-situ stress. When the difference of in-situ stress is small, the fracture network will be larger. Taking the minimum in-situ stress of 20 MPa as an example, when the injection pressure is 24 MPa, the damaged element is 4048 when the difference of in-situ stress is l ¼ 1.25, while the damaged element is 2140 when the difference of in-situ stress is l ¼ 1.5, and the number of fractures increases obviously when l ¼ 1.25 as can be seen in Fig. 18. 4.6.4. Heat transfer coefficient With the increase of heat transfer coefficient, the temperature of the rock around the wellbore decreases more in a certain period of time when it encounters the injection fluid of the same temperature, which makes the thermal stress enhance with the increase of

temperature gradient. As shown in Fig. 19, the increase in heat transfer coefficient makes the fracture network formed larger. As shown in Fig. 16(d), when the heat transfer coefficient is 1000 W/(m2$K), the crack initiation pressure is 7 MPa, and when the heat transfer coefficient is 2000 W/(m2$K), the crack initiation pressure is 5 MPa. The change of heat transfer coefficient not only affects the heat exchange between the cryogenic fluid and the hot rock matrix around the wellbore, but also has an important influence on the heat transfer between the fracture wall and the fracturing fluid during the fracturing process. Under the injection pressure of 12 MPa, the number of damaged units is 1876 when the heat transfer coefficient is 1000 W/(m2$K), and the number of damaged units increases to 7199 when the heat transfer coefficient is 2000 W/(m2$K). It can be seen that for a reservoir with a certain temperature, increasing the heat transfer coefficient between the injected fluid and the rock can more effectively utilize the low-temperature

W. Zhang et al. / Energy 178 (2019) 508e521



7L 7L 7L 7L

( *3D ( *3D ( *3D



'DPDJHG (OHPHQWV

'DPDJHG (OHPHQWV





   

   

   



 























,QMHFWLRQ 3UHVVXUH 03D











,QMHFWLRQ 3UHVVXUH 03D

(a) rock temperature

(b) Young's modulus 



   



K K K K

03D 03D 03D 03D

+ :Pg. + :Pg. + :Pg.

 

'DPDJHG (OHPHQWV

'DPDJHG (OHPHQWV

519



   

  





























,QMHFWLRQ 3UHVVXUH 03D













,QMHFWLRQ 3UHVVXUH 03D

(c) in-situ stress condition

(d) heat transfer coefficient

Fig. 16. Variation of damage elements with injection pressure during fracturing.

(a) E=33GPa

(b) E=37GPa

(c) E=41GPa

Fig. 17. Fracture morphology under different Young’s modulus (a) E ¼ 33 GPa, 15 MPa-82; (b) E ¼ 37 GPa, 14 MPa-22; (c) E ¼ 41 GPa, 13MPa-2.

induced thermal stress. 5. Conclusions In order to clarify the fracture morphology after HDR fracturing stimulation and provide the basis for accurate simulation of EGS heat mining, based on the meso-damage mechanics, elastic thermodynamics and Biot seepage mechanics, a mesoscopic THMdamage coupling model is proposed and implemented to analyze the initiation and propagation of HDR fracturing stimulation. Firstly, the mathematical model and numerical implementation method is validated by high-temperature granite fracturing experiments. Secondly, the action mechanism of thermal stress in HDR fracturing is discussed. Thirdly, the evolution of multiple physical fields during the initiation and propagation of HDR

fracturing is researched. Finally, the effects of various parameters on HDR fracturing process are also studied. The following conclusions have been drawn: (1) In the high-temperature granite hydraulic fracturing experiments, the rock failure pressure decreases as the raises of rock temperature. When the temperature raises from 50
C to 200
C, the rock failure pressure decreases by 27.7%. What’s more, the fracture network along the direction perpendicular to minimum in-situ stress can be achieved at 200
C. (2) The numerical simulation of HDR fracturing under thermal stress is in good agreement with the laboratory test both in the fracture morphology and rock failure pressure. The validation simulation also shows that the micro-damage

520

W. Zhang et al. / Energy 178 (2019) 508e521

(6) Under the same minimum in-situ stress condition, the increase of the in-situ stress difference makes the tendency of the fracture extending perpendicular to the minimum in-situ stress more obvious, and ultimately leads to the decrease of the expansion range of the fracture network.

(a)

=1.5

h=10MPa

(b)

=1.25

h=10MPa

In this study, the established mesoscopic THM-damage model provides the possibility for the prediction of HDR fracturing process. Adopting this numerical model to guide the fracturing project can reduce the cost and improve the success rate of EGS construction. Additionally, based on the obtained fracture morphology, it can provide powerful support for reasonable and accurate heat mining simulation of EGS. The current research is limited to smallscale two-dimensional model, and the purpose of which is to clarify the mechanism of fracture initiation and propagation in HDR fracturing under thermal stress. In the follow-up studies, the threedimensional model will be constructed and simulated for the fieldscale EGS project. In addition, for the HDR that already contains natural fractures, the impact of the presence of natural fractures on the fracturing of HDR needs to be investigated in the future. Acknowledgments

(c)

=1.5

h=20MPa

(d)

=1.25

h=20MPa

Fig. 18. Fracture morphology under different in-situ stress conditions (a) l ¼ 1.5 sh ¼ 10 MPa, 14 MPa-22; (b) l ¼ 1.25 sh ¼ 10 MPa, 14MPa-6; (c) l ¼ 1.5 sh ¼ 20 MPa, 28 MPa-34; (d) l ¼ 1.25 sh ¼ 20 MPa, 27 MPa-31.

This research was supported by the National Natural Science Foundation of China (51874338); the Fundamental Research Funds for the Central Universities (17CX06008, 17CX02077); the Shandong Natural Science Fund (ZR2016EEM30); the Applied Basic Research Project for Qingdao (17-1-1-20-jch). References

(a) 1000W/(m2· K)

(b) 2000W/(m2· K)

(c) 3000W/(m2· K)

Fig. 19. Fracture morphology with different heat transfer coefficients (a) Heat transfer coefficient is 1000 W/(m2$K), 16 MPa-34; (b) Heat transfer coefficient is 2000 W/ (m2$K), 14 MPa- 22; (c) heat transfer coefficient is 3000 W/(m2 K), 11 MPa-31.

mechanics analysis of rock failure process can truly depict the internal damage of rock, which provides a powerful tool for the study of rock damage process under the coupling of multiple physical fields. (3) The temperature gradient formed by the contact between cryogenic fracturing fluid and hot rock increases first and then decreases, so the thermal stress increases first and then decreases corresponding. (4) In the HDR hydraulic fracturing, the low-temperature induced thermal stress and the increasing injection pressure both act on the matrix in the form of tensile stress. The initiation and propagation of fractures are the joint action of the thermal stress and injection pressure. Therefore, both the initiation pressure and the failure pressure of hot rock are reduced during hydraulic fracturing. (5) With the increase of rock temperature, heat transfer coefficient and Young’s modulus, the more severe the thermal stress induced by low-temperature is, the lower the fracture initiation and rock failure pressure are, and a larger and more complex fracture network will be formed.

[1] Hofmann H, Weides S, Babadagli T, Zimmermann G, Moeck I, Majorowicz J, Unsworth M. Potential for enhanced geothermal systems in Alberta, Canada. Energy 2014;69:578e91. _ zał _ A. [2] Bujakowski W, Barbacki A, Miecznik M, Paja˛ k L, Skrzypczak R, Sowizd Modeling geothermal and operating parameters of EGS installations in the lower Triassic sedimentary formations of the central Poland area. Renew Energy 2015;80(2):441e53. [3] Zhang YJ, Guo LL, Li ZW, Yu ZW, Xu TF, Lan CY. Electricity generation and heating potential from enhanced geothermal system in Songliao Basin, China: different reservoir stimulation strategies for tight rock and naturally fractured formations. Energy 2015;93:1860e85. [4] Sun ZX, Zhang X, Xu Y, Yao J, Wang HX, Lv SH, Sun ZL, Huang Y, Cai MY, Huang XX. Numerical simulation of the heat extraction in EGS with thermalhydraulic-mechanical coupling method based on discrete fractures model. Energy 2017;120:20e33. [5] Qu ZQ, Zhang W, Guo TK. Influence of different fracture morphology on heat mining performance of enhanced geothermal systems based on COMSOL. Int J Hydrogen Energy 2017;42:18263e78. [6] Zhao Yangsheng, Wan Zhijun, Zhang Yuan, Zhang Ning, Feng Zijun, Dong Fuke, Wu Jinwen, Qu Fang. Experimental study of related laws of rock thermal cracking and permeability. Chin J Rock Mech Eng 2010;29(10): 1970e6 [in Chinese)]. [7] Feng Zijun, Zhao Yangsheng, Zhang Yuan, Wan Zhijun. Real-time permeability evolution of thermally cracked granite at triaxial stresses. Appl Therm Eng 2018;133:194e200. [8] Zhao Zhihong. Thermal influence on mechanical properties of granite: a microcracking perspective. Rock Mech Rock Eng 2016;49:747e62. [9] Rutqvist J, Jeanne P, Dobson PF, Garcia J, Hartline C, Hutchings L, Singh A, Vasco DW, Walters M. The northwest geysers EGS demonstration project, California-Part 2: modeling and interpretation. Geothermics 2015;63(1): 120e38. [10] Garcia J, Hartline C, Walters M, Wright M, Rutqvist J, Dobson PF, Jeanne P. The northwest geysers EGS demonstration project, California-Part 1: characterization and reservoir response to injection. Geothermics 2015;63:97e119. [11] Vishkai Mahta, Gates Ian. On multistage hydraulic fracturing in tight gas reservoirs: montney Formation, Alberta, Canada. J Pet Sci Eng 2019;174: 1127e41. [12] ErfanSaberhosseini Seyed, Ahangari Kaveh, Mohammadrezaei Hossein. Optimization of the horizontal-well multiple hydraulic fracturing operation in a low-permeability carbonate reservoir using fully coupled XFEM model. Int J Rock Mech Min Sci 2019;114:33e45. [13] Jaber Taheri-Shakib, Amir Ghaderi, Mohammad AminSharif Nik. Numerical study of influence of hydraulic fracturing on fluid flow in natural fractures. DOI:10.1016/j.petlm.2018.10.005.

W. Zhang et al. / Energy 178 (2019) 508e521 [14] Li Tianjiao, Li Lianchong, Chun’anTang, et al. A coupled hydraulic-mechanicaldamage geotechnical model for simulation of fracture propagation in geological media during hydraulic fracturing. J Pet Sci Eng 2019;173: 1390e416. [15] Cruz Francisco, Roehl Deane, Vargas Jr Euripedes do Amaral. An XFEM element to model intersections between hydraulic and natural fractures in porous rocks. Int J Rock Mech Min Sci 2018;112:385e97. [16] Ma Xinfang, Zou Yushi, Ning Li, et al. Experimental study on the mechanism of hydraulic fracture growth in a glutenite reservoir. J Struct Geol 2017;97: 37e47. [17] Xi Baoping, Zhao Yangsheng. Experimental research on mechanical properties of water-cooled granite under high temperatures within 600
C. Chin J Rock Mech Eng 2010;29(05):892e8. [18] Kumari WGP, Rangith PG, Perera MSA, Chen BK, Abdulagatov IM. Temperature-dependent mechanical behavior of Australian Strathbogie granite with different cooling treatments. Eng Geol 2017;229:31e44. [19] Kumari WGP, Rangith PG, Perera MSA, Chen BK. Experimental investigation of quenching effect on mechanical, microstructural and flow characteristics of reservoir rocks: thermal stimulation method for geothermal energy extraction. J Pet Sci Eng 2018;162:419e33. [20] Zhao Yangsheng, Feng Zijun, Zhao Yu, et al. Experimental investigation on thermal cracking, permeability under HTHP and application for geothermal mining of HDR. Energy 2017;32:305e14. [21] Li Zheng-Wei, Feng Xia-Ting, Zhang Yan-Jun, et al. Experimental research on the convection heat transfer characteristics of distilled water in manmade smooth and rough rock fractures. Energy 2017;133:206e18. [22] Zhao Yangsheng, Feng Zijun, Feng Zengchao, et al. THM (Thermo-hydro-mechanical) coupled mathematical model of fractured media and numerical simulation of a 3D enhanced geothermal system at 573K and buried depth 6000-7000M. Energy 2015;82:193e205. [23] Guo Liangliang. Test and model research of hydraulic fracturing and reservoir damage evolution in Enhanced Geothermal System[D]. Ji Lin University; 2016 [in Chinese)]. [24] Zhou Changbing, Wan Zhijun, Zhang Yuan, et al. Experimental study on hydraulic fracturing of granite under thermal shock. Geothermics 2018;71: 146e55. [25] Kumari WGP, Ranjith PG, Perera MSA, et al. Hydraulic fracturing under high temperature and pressure conditions with micro CT applications: geothermal energy from hot dry rocks. Fuel 2018;230:138e54. [26] Cha Minsu, Yin Xiaolong, Kneafsey Timothy, et al. Cryogenic fracturing for reservoir stimulation-Laboratory studies. J Pet Sci Eng 2014;124:436e50. [27] Zhang Shikun, Huang Zhongwei, Huang Pengpeng, Wu Xiaoguang, et al. Numerical and experimental analysis of hot dry rock fracturing stimulation with high-pressure abrasive liquid nitrogen jet. J Pet Sci Eng 2018;163:156e65. [28] Li lianchong. Investigation on numerical model of coupled thermo-hydromechanical-damage (THM-D) for rock failure process and associated application [D]. Northeastern University; 2006 [in Chinese)]. [29] Li Lianchong, Yang Tianhong, Tang Chunan, et al. Coupling analysis of hydraulic fracturing process of rock. Chin J Rock Mech Eng 2003;22(7):1060e6 [in Chinese)]. [30] Rayudu Nithin Manohar, Tang Xuhai, Singh Gaurav. Simulating three dimensional hydraulic fracture propagation using displacement correlation

521

method. Tunn Undergr Space Technol 2019;85:84e91. [31] Li Zhichao, Li Lianchong, Huang Bo, et al. Numerical investigation on the propagation behavior of hydraulic fractures in shale reservoir based on the DIP technique. J Pet Sci Eng 2017;154:302e14. [32] Lian-chong LI, Gen LI, Meng Qing-min, Wang Hao, Wang Zhen. Numerical simulation of propagation of hydraulic fractures in glutenite formation. Rock Soil Mech 2013;34(5):1501e7 [in Chinese)]. [33] Wei CH, Zhu WC, Yu QL, Xu T, Jeon S. Numerical simulation of excavation damaged zone under coupled thermalemechanical conditions with varying mechanical parameters. Int J Rock Mech Min Sci 2015;75:169e81. [34] C Zhu W, Wei J, Zhao J, et al. 2D numerical simulation on excavation damaged zone induced by dynamic stress redistribution. Tunn Undergr Space Technol 2014;43:315e26. [35] Xiang Li, Wang Jiehao, Elsworth Derek. Stress redistribution and fracture propagation during restimulation of gas shale reservoirs. J Pet Sci Eng 2017;154:150e60. [36] Pandey SN, Chaudhuri A, Kelkar S. A coupled thermo-hydro-mechanical modeling of fracture aperture alteration and reservoir deformation during heat extraction from a geothermal reservoir. Geothermics 2017;65:17e31. [37] Gelet Rachel, Loret Benjamin, Khalili Nasser. A thermo-hydro-mechanical coupled model in local thermal non-equilibrium for fractured HDR reservoir with double porosity. J Geophys Res 2012;117:1e23. [38] Taron Joshua, Elsworth Derek. Thermalehydrologicemechanicalechemical processes in the evolution of engineered geothermal reservoirs. Int J Rock Mech Min Sci 2009;46:855e64. [39] Sesetty Varahanaresh, Ahmad Ghassemi. Effect of rock anisotropy on wellbore stresses and hydraulic fracture propagation. Int J Rock Mech Min Sci 2018;112:369e84. [40] COMSOL Multiphysics User’s guide, version 5.3. Massachusetts, USA, Burlington: COMSOL Inc; 2017. www.comsol.com. [41] Ghassemi A, Zhang Q. A transient fictitious stress boundary element method for porothermoelastic media. Eng Anal Bound Elem 2004;28(11):1363e73. [42] Guo Tiankui, Li Yanchao, Ding Yong, Qu Zhanqing, Gai Naicheng, Rui Zhenhua. Evaluation of acid fracturing treatments in shale formation. Energy Fuels 2017;31(10):10479e89. [43] Guo Tiankui, Zhang Shicheng, Ge Hongkui, Wang Xiaoqiong, Lei Xin, Xiao Bo. A new method for evaluation of fracture network formation capacity of rock. Fuel 2015;140:778e87. [44] Enayatpour Saeid, Oort Ericvan, Patzek Tad. Thermal cooling to improve hydraulic fracturing efficiency and hydrocarbon production in shales. J Nat Gas Sci Eng 2019;62:184e201. [45] Zhang Shikun, Huang Zhongwei, Li Gensheng, et al. Numerical analysis of transient conjugate heat transfer and thermal stress distribution in geothermal drilling with high-pressure liquid nitrogen jet. Appl Therm Eng 2018;129:1348e57. [46] Guo Tiankui, Zhang Shicheng, Qu Zhanqing, et al. Experimental study of hydraulic fracturing for shale by stimulated reservoir volume. Fuel 2014;128: 373e80. [47] Chen Yun, Ma Guowei, Li Tuo, et al. Simulation of wormhole propagation in fractured carbonate rocks with unified pipe-network method. Comput Geotech 2018;98:58e68.