Stratified rock hydraulic fracturing for enhanced geothermal system and fracture geometry evaluation via effective length

Journal Pre-proof Stratified rock hydraulic fracturing for enhanced geothermal system and fracture geometry evaluation via effective length Likui Yu, Xiaotian Wu, Yadan Wang, Weiwu Ma, Gang Liu PII:

S0960-1481(20)30119-1

DOI:

https://doi.org/10.1016/j.renene.2020.01.097

Reference:

RENE 12960

To appear in:

Renewable Energy

Received Date: 29 October 2019 Revised Date:

6 January 2020

Accepted Date: 22 January 2020

Please cite this article as: Yu L, Wu X, Wang Y, Ma W, Liu G, Stratified rock hydraulic fracturing for enhanced geothermal system and fracture geometry evaluation via effective length, Renewable Energy (2020), doi: https://doi.org/10.1016/j.renene.2020.01.097. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Ltd.

Credit Author Statement Lkui Yu: Conceptualization, Validation, Writing-Original draft. Xiaotian Wu: Conceptualization, Methodology, Writing-Original draft, Formal analysis. Yadan Wang: Software, Investigation . Weiwu Ma: Data curation, Investigation. Gang Liu: Supervision, Writing-review, Funding acquisition.

Manuscript Submission Renewable energy

Stratified rock hydraulic fracturing for enhanced geothermal system and fracture geometry evaluation via effective length

Likui Yu, Xiaotian Wu, Yadan Wang, Weiwu Ma, Gang Liu* (School of Energy Science and Engineering, Central South University, Changsha, 410083, China)

*Corresponding author: Gang Liu Tel: +86 180 7516 9583 Fax: 0086 0731 88879863 Email: [email protected] Address: Energy building, School of Energy Science and Engineering, Central South University, Changsha, Hunan, China, 410083.

Submitted: January 25, 2020

Abstract Rock stratification and anisotropy extensively exist in geothermal reservoir, and have great effects on the results of hot dry rock hydraulic fracturing. Therefore, the understanding of fracturing mechanism is critical to successfully building an engineered reservoir with high permeability for enhanced geothermal 5

system. Some previous researchers have attempted to reveal hydraulic fracturing process of stratified rock, but the rock hydraulic fracturing propagation and fracture quality evaluation are still not clear enough. We established a numerical model to study the propagation of non-planar fractures in stratified rock by using extended finite element method. For fracture quality evaluation, we proposed fracture effective length as an index, which is the length of the segment where fracture width greater than 1 millimetre, considering the

10

volume of sand within injected proppant. It is found that fractures generated from stratified rock is more conducive to heat extraction process. A large in-situ stress difference is more conducive to an increase in the fracture effective length. The methodology and results would guide the perforation location selection for stractified rock hydraulic fracturing process to achieve a reasonable fracture morphology. Keywords: Enhanced geothermal system; Hydraulic fracturing; Hot dry rock; Stratified rock; Fracture

15

morphology; In-situ stress difference

Preprint submitted to

January 25, 2020

1. Introduction Hot dry rock (HDR) has a high density and a low porosity and permeability, and it must be artificially modified to perform heat extraction [1, 2]. Enhanced geothermal system (EGS) is an artificial system 20

in which engineering technology is used to artificially construct fracture networks or improve original fracture networks in HDR to form geothermal reservoirs [3, 4]. Hydraulic fracturing is the most important technology for building artificial geothermal reservoirs. Mastering the fracturing and spreading laws of hydraulic fractures plays a crucial role in the utilization of deep geothermal energy [5]. Hydraulic fracturing is a key technique for heat extraction from HDR, and its purpose is to create

25

fracture networks as flow paths in geothermal reservoirs between injection and production wells [6, 7]. Hydraulic fracturing involves complex coupling of multiple physics and different scale problems [8, 9]. Even under the premise of many assumptions, it is difficult to theoretically clarify the mechanism of hydraulic fracturing [10]. Due to the difficulties in monitoring internal fracture network of rock, numerical simulation is still an important tool for engineers to estimate the fracture network geometry [11, 12].

30

In recent years, most widely used hydraulic fracturing modeling methods focused on the displacement discontinuity method (DDM) [13, 14], the discrete element method (DEM) [15–17], the finite element method (FEM) [18, 19] and the extended finite element method(XFEM) [20, 21]. DDM uses the discontinuous displacement as a function to be sought, which can describe the fracture elements naturally. Hydraulic fracturing would break the equilibrium state of in-situ stress, resulting in stress

35

redistribution [22]. The limitations of DDM itself make the analysis objects limited to linear elastic materials and cannot consider the coupling of heat conduction and deformation. The model based on DEM can reflect the discontinuity and joint rationality of rock [16, 23]. However, the theory of DEM is not yet systematic and lacks the support of mathematical theory. The calculation speed is closely related to the partition of elements. When the number of DEM computing elements is large, the calculation cost

40

would get expensive [24]. Based on the continuum theory, FEM divides the object into finite elements and links nodes with mathematical equations [25]. When the rock contains joint fractures, the FEM calculation results have some errors. Besides, using this method to simulate hydraulic fracturing is very time-consuming. The combination of FEM and cohesive zone method causes the fractures to spread only

2

along the edge of elements, thereby reducing the computational cost [19, 26]. 45

But the fracture

propagation path is also restricted in this case. In order to overcome the shortcomings of FEM, XFEM is proposed, which can be seen as an improvement of FEM. When XFEM is used to simulate the fracture propagation, the meshes do not need to be reconstructed and fractures are allowed to propagate along any path [27, 28]. Therefore, compared with the traditional FEM, using XFEM can significantly reduce the calculation cost.

50

Because of those superiorities of XFEM, many scholars have used this method to study the morphology and propagation process of hydraulic fractures. For example, the effect of natural fractures on hydraulic fracture propagation varies due to formation conditions differences [29]. When horizontal wells are multiperforated for simultaneous fracturing or sequential fracturing, stress shadows between fractures can affect the formation of complex fracture networks [30, 31]. The percolation and suction of fracturing fluid in

55

medium are investigated. Sensitivity analyses are carried out for the cases with different injection rates [32]. However, rock anisotropy, which has an important influence on the morphology of hydraulic fractures, is ignored in these studies. In fact, rock properties are generally anisotropic [33, 34]. Therefore, XFEM has also been used to study hydraulic fracturing in the case of rock anisotropy. For example, the influence factors of hydraulic fracture morphology are analyzed under the rock anisotropy conditions (such as: in-

60

situ stress difference, ratio of Young’s modulus, fracture toughness, rock material angle, etc.) [35, 36]. Rock material angle is the angle between initial fracture direction and anisotropic material axis. The propagation direction of hydraulic fractures is affected by material angle and in-situ stress [37]. Although predecessors have studied the effect of rock properties on hydraulic fracturing results, the situation of rock stratification was ignored in previous studies. Actually, rock anisotropy is only one aspect of rock

65

complexity. HDR may be stratified due to long-term geological loads (Fig. 1(a)). As a result, rock properties change greatly in different strata [38–40]. The fracture initiation pressure is mainly determined by the material stiffness and failure coefficient. When hydraulic fracturing is performed on stratified rock, the micro-fractures first appear in the part with small material failure coefficient, which increases the rock permeability. Moreover, due to the effect of high pressure on stratified rock, larger deformations occur in

70

the part with smaller material stiffness. Therefore, the rock stratification has a great influence on fracture

3

morphology and initiation pressure. Due to the lack of previous researches in this aspect, it is necessary to study the initiation and propagation mechanism of hydraulic fractures in stratified rock. In this study, the entire HDR is divided into rock mass I and II (Fig. 1(b)). Although they are isotropic in themselves, their properties are different from each other. The material stiffness and material failure parameters of 75

rock mass II are 1/10 of that of rock mass I. The fluid filtration rate and permeability coefficient of rock mass II are consistent with that of rock mass I.

(a)

(b)

Fig. 1. Diagram of stratified rock. (a) Texture structure of rock stratification; (b) 2D geometric model of rock stratification.

In recent years, many studies attempted to reveal the hydraulic fracturing mechanism of porous rock. For example, Amir et al. studied the flow of production wells with different apertures and fracture lengths [17]. The results show that the effect of fracture length on production is much greater than the aperture. 80

Paul et al. established a fluid-solid coupling model to study the interaction between adjacent fractures [27]. However, the analysis of hydraulic fracturing results in previous studies was limited to understand the impacts of fracture length and width on practical reservoir performance, which is quite different from practical projects. In practical projects, a proppant composed of water, sand and chemicals is injected into hydraulic fractures to prevent those from closing [41]. As the sand clogs the tail of fractures, only a

85

portion of fractures can be used for heat extraction. In this study, we defined the fracture length of this portion as the fracture effective length. The effective fracture length, rather than total length, is used to evaluate the fracture quality.

4

In summary, a fluid-solid coupling model based on XFEM is used to study the hydraulic fracturing process of stratified rock. By comparing the hydraulic fracturing results of isotropic and stratified rock, 90

the initiation and propagation mechanisms of hydraulic fractures in stratified rock are analyzed. Then, the model is used to study impact factors of hydraulic fracturing, including the in-situ stress difference and the distance between rock mass II and perforation. Among them, in-situ stress difference is one of the key parameters of hydraulic fracturing, and the distance between rock mass II and perforation is one of the main factors affecting hydraulic fracturing of stratified rock. Finally, by considering the injection of

95

proppant, the fracture effective part is used as an indicator to evaluate the fracture quality.

2. Mathematical framework and computer implementation Hydraulic fracturing is a dynamic coupling process between fracturing fluid flow and rock deformation. The mathematical framework of hydraulic fracturing consists of fluid-solid coupling equations, the fracture initiation and propagation criteria, and the extended finite element method.

100

2.1. Fluid-solid coupling equations Due to the existence of fluid pressure in pores, the rock stress is shared by the rock skeleton and void fluid, and only the effective stress transmitted through the rock skeleton can cause this rock to deform. This is the Terzaghi effective stress principle. At some points, the virtual work on rock mass is equal to the virtual work generated by the force acting on rock mass. The following equation is the stress equilibrium

105

equation of rock mass, also known as virtual work equation. Expressed by eq. 1:

Z V

δT dδdV −

Z

δuT df dV −

V

Z

δuT dtdS = 0

(1)

S

where t is the surface force, f is the physical force; δu is the virtual strain, and δ is the virtual displacement. For rock with constant volume, the following conclusion can be obtained according to the law of conservation of mass. Within a certain period of time, the fluid mass flowing into rock is equal to the sum

5

of increment of fluid mass in rock and the fluid mass flowing out of rock [42]. Expressed by eq. 2:

∂ρw φ + 5 · (ρw φvw ) = 0 ∂t 110

(2)

where ρw is the fluid density in porous medium; φ is the porosity; vw is the percolation velocity. The fluid flow in hydraulic fractures includes tangential flow and normal filtration. The fracturing fluid is assumed to be incompressible Newtonian fluid. The tangential velocity depends on the lubrication equation [43], which can be written as follows:

q=−

w3 5 pf 12µ

(3)

where pf is the fluid pressure in fracture, i.e. the traction force of fracture propagation; w is the fracture 115

width; µ is the fluid viscosity. While the fluid in fractures flows tangentially, part of the fluid is lost to stratigraphic medium through the upper and lower surfaces of fractures (Fig. 2). The normal flow formulas of fluid perpendicular to the upper and lower surfaces of fractures are as follows:     qt = ct (pf − pt )

(4)

  qb = cb (pf − pb ) where pf and pb are the pore pressures of rock mass near fracture surfaces; cb and ct are the fluid loss 120

coefficients.

6

Fig. 2. Schematic diagram of fluid flow inside the fracture.

2.2. Criteria for fracture initiation and propagation Damage initiation is the starting point for the reduction of material mechanical properties. When the stress or strain relationship satisfies specified initial damage criterion, the material stiffness begins to degrade and load carrying capacity decreases. Thus, the fracture initiation occurs. In this paper, the 125

criterion used in numerical simulation is the maximum principal stress criterion. The criterion is that once the maximum principal stress reaches critical stress, the material would be damaged. Its expression is shown as follows [44, 45]:

f=

hσmax i 0 σmax

(5)

0 is the maximum principal stress; h i is the Macaulay bracket, indicating that compressive where σmax

stress does not damage material. When value of f is 1, the rock mass begins to appear damaged. 130

For the expansion process of composite fractures, we used the critical energy release rate criterion (B-K criterion) for fracture propagation proposed by Benzeggagh and Kenane [45, 46].

C C GC n + (Gs − Gn )



GS GT

η

= GC

(6)

where GS = Gs + Gt , GT = Gn + Gs , Gcn is the critical strain energy release rate of normal fracture; GC C is the critical fracture energy release rate of composite fracture; GC s , Gt are the critical energy release

7

rates of first tangential and second tangential fractures, respectively; η is a constant related to material 135

properties.

2.3. Extended finite element method Belytschko and Black [47, 48] first proposed the extended finite element method (XFEM) in 1999 to deal with the discontinuity mesh problem in numerical simulation. XFEM is more efficient than other numerical methods in solving discontinuous mechanics problems. Based on the theory of FEM, XFEM 140

effectively inherits all advantages of FEM based on the partition of unity method, and does not need to re-mesh geometric interface generated by the existence of fractures in grid elements. The displacement and pore pressure of nodes in elements can be calculated by the modified interpolation function.                  h       a1j  ui  u (x)                    X    X X T h + Nm (x) [L] + N (x)H(x) = N (x) j i a v v (x) 2j i             m∈ΩΛ j∈ΩΓ i∈Ω                        pi    ph (x)   a3j 

     tip   um          tip vm               ptip m

(7)

where H(x) is the Heaviside function of fracture surfaces. It can be written as:

H(x) =

   +1 x ≥ 0

(8)

  −1 x < 0 145

where N i is the shape function; Ω is the area where the entire element is located; ΩΓ is the mesh area through which fracture penetrates; ΩΛ is the element where the fracture tip is located; (ui , vi ) is the displacement of joint; (a1j , a2j ) is the enhanced degree of freedom associated with Heaviside function tip ; (utip m , vm ) is the progressive displacement field of fracture tip; [L] is the transition between the local

coordinates of fracture tip and global coordinates.

150

2.4. Computer implementation The main processes of simulating hydraulic fracturing using the extended finite element method is shown in Fig. 3. Material properties and initial fractures locations are defined to establish HDR geometric

8

models. The relative displacement of fracture surfaces in this model is solved. If this process does not converge, the models need to be optimized. Then, the energy release rate is solved. According to Eqs. 5 155

and 6, when f is 1, the fractures begin to propagate along the direction of maximum in-situ stress. Finally, if the calculation time does not reach the set time, new fracture tip positions is calculated to perform the above process again. When the calculation time reaches the set time, the calculation ends and results are output.

Fig. 3. Algorithm and flowchart of extended finite element method for hydraulic fracturing.

9

3. Model introduction and verification

160

3.1. Model introduction In this study, a two-dimensional fluid-solid coupling numerical model for the propagation of non-planar fractures in porous media is established by using XFEM, taking into account of rock stratification. The length and width of HDR model are 40 m, and the thickness of rock mass II is 2 m. The Young’s modulus of rock mass I is 30 GPa and the tensile strength is 6 MPa. The value of material stiffness and material

165

failure parameters of rock mass II are 1/10 of that of rock mass I. The fluid loss rate and permeability coefficient are the same as those of rock mass I. In order to make the simulation results accurate and convergent, the coarse element is used in the far field and fine element is used in the near field. For studying the impacts of fracture morphology on heat extraction, we chose the research time point is at 120 seconds, which is equal to the total injection time. This is because heat extraction only occurs

170

when hydraulic fracturing finished. For studying the influence of factors (rock stratification, in-situ stress difference, the distance between rock mass II and perforation) on fracture initiation and propagation, we set the research time resolution as 10 seconds. When fracturing, the processes of initiation, expansion and mergence of micro-cracks constantly accumulated with time. The time resolution of 10-second is accurate enough to demonstrate these phenomenon.

175

3.2. Model validation The accuracy of model in this study is verified by comparing numerical solution of hydraulic fracturing with experimental and analytical solutions. First, the numerical solution obtained from this model is compared with the experimental solution in existing literature [49, 50]. Lv carried out experiments on the hydraulic fracturing process of granite. The properties of granite are basically the same as those

180

of HDR. In order to be consistent with the experiment, the injection time is increased in the numerical simulation. The fracture internal pressure in numerical and experimental solutions is analyzed. Then, the numerical solution of this model is compared with the analytical solution of Perkins-Kern-Nordgren (PKN) model [51, 52]. The PKN model is a more mature two-dimensional model for hydraulic fracturing, and it is suitable for the situations of large aspect ratio. The fracture length and width can be accurately

10

185

determined by the PKN model. The maximum fracture width between numerical and analytical solutions is compared. Due to the low permeability of HDR, additional pumping pressure is required to cause the fracture initiation. The peak pressure in Fig. 4(a) is initiation pressure of fractures. After fractures initiation, the fracture extension does not require additional pumping pressure. The extension pressure of fractures is less

190

than the one of initiation. The fracture internal pressure changes of numerical and experimental solutions are approximately the same. Around the peak pressure, the numerical solution is somewhat different from the experimental one. This phenomenon is mainly because the model is two-dimensional, and the rock in the experiment is three-dimensional. There are pressure build-up and relief phenomena during the actual fracturing process. The fracture internal pressure is slightly increased during the pressure build-

195

up process, and it is slightly reduced during the pressure relief process. Therefore, the curve of fracture internal pressure with time is not smooth (Fig. 4(a)). In early stage of hydraulic fracturing, the maximum fracture width of numerical solution is basically consistent with that of analytical solution. The maximum fracture width of numerical solution at the end of hydraulic fracturing is slightly larger (Fig. 4(b)). This is mainly because PKN model assumes that the rock is infinite, and the model in numerical simulation has

200

boundaries. The fracture is basically unaffected by boundaries in early stage. As the fracture expands, the influence of boundary effect on it cannot be ignored. Due to the existence of pressure build-up and relief processes, the curve of fracture maximum width with time in actual process is fluctuating, and the overall trend is increasing. In summary, the numerical solution is generally consistent with experimental and analytical solutions, in spite of small errors which are fully acceptable. Therefore, the simulation

205

results are reliable and can be used as a reference for fracturing design.

11

(a)

(b)

Fig. 4. Comparison of hydraulic fracturing results between numerical solution and experimental and analytical solutions. (a) Relationship between fracture internal pressure and time; (b) relationship between maximum fracture width and time.

4. Results and discussions 4.1. Effects of rock stratification on hydraulic fracturing results The results of hydraulic fracturing in two conditions of rock isotropy and stratification are compared. The maximum and minimum in-situ stress in fracturing process are 40 MPa and 38 MPa, respectively. 210

The flow rate of fracturing fluid is 0.001 m3 /s. In the fracturing simulation of stratified rock, the distance between rock mass II and fracture perforation is 2 m. The results obtained are shown as follows. During the hydraulic fracturing process, as the fracturing fluid is injected, the circumferential stress on bottom rock is changed from a compressed state to a tensile state. When the tensile stress is large enough to overcome rock tensile strength, the rock original fractures would be opened and extended or

215

new fractures would be formed. Figs. 5(a) and 5(b) show the results obtained by hydraulic fracturing of isotropic rock. It has been assumed that the rock is a homogeneous, isotropic, linear elastic porous material with no stratification. The hydraulic fracture extends in the direction of maximum in-situ stress without deflection.

12

(a)

(b)

(c)

(d)

Fig. 5. Comparison of hydraulic fracturing results between rock stratification and isotropy. (a) Pore pressure of isotropic rock; (b) Deformation distribution of isotropic rock; (c) pore pressure of stratified rock; (d) deformation distribution of stratified rock.

The generation of fractures in rock is not a one-step process, but a process of gradual deterioration. 220

Most rock masses have produced micro-cracks and micro-holes before macro-cracks appear. Due to the influence of various external factors, micro-cracks and micro-holes are constantly fracturing, expanding, merging and penetrating, which deteriorates the performance of rock mass and causes damage. There is a process from quantitative change to qualitative change in the damage of rock masses. The material stiffness and material failure parameters of rock mass II are both 1/10 of that of rock mass I, which means

225

that the initiation pressure of rock mass II is far less than that of rock mass I. As fracturing fluid leaks into formation, the rock mechanical effects around formation increase. When the pressure on rock mass II is

13

greater than the fracture initiation pressure and rock mass I is not, micro-cracks are first formed and large deformation occurs in rock mass II (Fig. 5(d)). Due to large deformation of rock mass II, the fracture would bend toward rock mass II while extending along the direction of maximum in-situ stress (Fig. 6). 230

Moreover, the initiation pressure of stratified rock is significantly lower than that of isotropic rock.

Fig. 6. Comparison of fracture morphology between rock stratification and isotropy.

As the fracturing fluid leaks into formation, the variation of pore pressure would change the effective stress and fluid leakage speed. The direction of fracture propagation is also changed. The variation of fracture width with fracture length is shown in Fig. 7. In the case of rock isotropy, the fracture continues to extend along direction of maximum in-situ stress. For hydraulic fracturing of stratified rock, the fracture 235

would bend toward rock mass II while extending along direction of maximum in-situ stress. The simulation results show that the fracture of stratified rock is shorter and wider than that of isotropic rock. Compared with isotropic rock, the fracture length in stratified rock decreased by 57.7%, and the maximum fracture width increased by 110.8%. This is due to the fact that in process of injecting fluids into stratified rock, the smaller value of material stiffness and material failure of rock mass II would lead to micro-cracks and

240

large deformation. A large deformation in Y direction would result in fractures with greater width. In addition, the presence of micro-cracks causes more fracturing fluid to leak into formation rock, resulting

14

in less fluid to squeeze the rock. This is also the reason why the fracture length in stratified rock is much smaller than that in isotropic rock. In the process of two kinds of hydraulic fracturing, fracturing fluid is injected into HDR at a certain rate within a certain period of time, and the total volume of fracturing 245

fluid is the same. When fracture width is large, the length is relatively small. So two results obtained by simulation are also accord with actual situation of projects. The mechanism of initiation and propagation of hydraulic fractures in stratified rock is studied by comparing isotropic rock. The research results can provide a new perspective for the prediction of fracture morphology during construction process.

Fig. 7. Relationships between fracture width and length.

The hydraulic fracturing process involves the injection of fracturing fluid and proppant. During 250

proppant injection process, the fracture width must be larger than particle size in proppant. Otherwise, it would cause sand plugging, and it is difficult to extract heat. In previous studies, scholars used the fracture total length to evaluate the quality of fractures, which differed greatly from actual engineering. The part of fracture length that can be normally injected into the proppant should be used as an evaluation standard instead of the total fracture length. In this study, the minimum limit of fracture width that can

255

be normally injected into proppant is set to 1 mm. The working fluid can only perform convective heat transfer in the part where the fracture width is greater than 1 mm. The remaining part has no effect on heat transfer due to sand plugging. The part with a fracture width greater than 1 mm was defined as

15

fracture effective part, and its length is fracture effective length. The fracture effective length in isotropic rock and stratified rock are 3.7 m and 6.6 m, respectively. The purpose of hydraulic fracturing is to extract 260

deep geothermal energy, and only the fracture effective part can normally exchange heat. Therefore, the fractures in stratified rock can meet the requirements of actual engineering better, although the fracture total length in isotropic rock is larger. The fracture effective length proposed in this study makes a more accurate evaluation of the heat extraction of different fractures, and provides a scientific basis for the optimization of fracture morphology in actual engineering.

265

4.2. Effects of the distance between rock mass II and perforation (d) on hydraulic fracturing results The results of hydraulic fracturing are compared by changing the distance between rock mass II and perforation. During the numerical simulation, the perforation position is kept constant. The value of d is changed by changing the position of rock mass II, where d is the distance between rock mass II and perforation. Let d be 2 m, 4 m, 6 m, 8 m, 10 m respectively. The maximum and minimum in-situ stresses

270

in fracturing process are 40 MPa and 38 MPa, respectively. The fracturing fluid flow rate is 0.001 m3 /s. The results obtained are reproduced below. By changing the value of d, the comparison of rock pore pressure distribution after hydraulic fracturing is shown in Fig. 8. It can be seen from the figure that the pore pressure in rock mass II is generally greater than that in rock mass I. The maximum value of pore pressure is also present in rock mass II. This situation

275

is due primarily to the fact that the stiffness and failure coefficient of rock mass II is much smaller than that of rock mass I. During the injection of fracturing fluid, the rock mass II would produce micro-cracks and deform early, and the liquid would enter the rock by seepage. Therefore, the pore pressure in rock mass II is increased, which is significantly larger than that in rock mass I.

16

(a)

(b)

(c)

(d)

Fig. 8. Effects of d values on fracture results. (a) d = 4 m; (b) d = 6 m; (c) d = 8 m; (d) d = 10 m.

It can be seen from Fig. 10(b) that as value of d increases, the fracture length increases, the maximum 280

fracture width decreases, and the angle at which fracture deflects toward rock mass II is significantly reduced. The fracture morphology tends to be “narrow, long and straight” as the value of d increases (Fig. 9). When the value of d increases from 2 m to 10 m, the fracture length increases by 101.2% and the fracture width decreases by 46.1%. Comparing different values of d in Fig. 8, it can be seen that the rock pore pressure generally decreases with the value of d. This is mainly because the stress effect on rock

285

mass II decreases with the value of d, which leads to the reduction of rock deformation. The working fluid leaking into rock is also reduced, resulting in a decrease in rock pore pressure. The decrease of rock mass II deformation leads to the reduction of total rock deformation, so that the fracture width and deflection 17

angle are reduced. Since the fracturing fluid leaking into rock is reduced, more fracturing fluid can be used to crack the rock, increasing the fracture length.

Fig. 9. Relationship between fracture morphology and d values.

290

In Fig. 10(a), the fracture effective part is above horizontal dashed line. The fracture effective part can normally inject proppant and extract heat. We can see from Fig. 10(a) that the fracture effective length increases first and then decreases with the value of d. When the value of d is 6 m, the fracture effective part is the longest. In the other word, when the value of d is 6 m, the fracture morphology is more in line with the need for heat extraction in actual engineering. The research results can provide an

295

engineering guidance for the selection of perforation positions in stratified rock hydraulic fracturing.

18

(a)

(b)

Fig. 10. Effects of d values on hydraulic fracturing results. (a) Relationship between fracture width and fracture length; (b) relationship between fracture length and d values, maximum fracture width and d values.

4.3. Effects of in-situ stress difference (∆σ) on hydraulic fracturing results The maximum in-situ stress is maintained at 40 MPa and minimum in-situ stress is changed so that the value of ∆σ is 0 MPa, 2 MPa, 4 MPa, 6 MPa, 8 MPa, and 10 MPa, respectively. The results of hydraulic fracturing are compared and analyzed. In the process of fracturing, the distance between rock 300

mass II and perforation is 2 m. Fracturing fluid flow rate is 0.001 m3 /s. The results are as follows. The rock initiation pressure decreases with the value of ∆σ (Tab. 1). This means that in the case of large ∆σ, rock can rupture under lower pressure. This situation occurs mainly because the fluid pressure must overcome in-situ stress to make it possible to expand the fractures. Therefore, in-situ stress is a hindrance for fracture extension. When the value of ∆σ increases, that is, when the minimum in-situ

305

stress decreases, the pressure required for fracture extension decreases, and the rock is more susceptible to cracking. Table 1 Relationship between initiation pressure and ∆σ values.

∆σ (MPa)

0

2

4

6

8

10

Initiation pressure (MPa)

25.81

24.41

22.99

22.49

21.98

20.58

The fracture morphology can be observed in Fig. 11. As the value of ∆σ increases, the angle at 19

which the fracture deflects toward rock mass II is significantly reduced, and the fracture tendency to expand in the direction of maximum in-situ stress is more pronounced. The variation of fracture width 310

with fracture length under different values of ∆σ is shown in Fig. 12(a). As shown in Fig. 12(b), the fracture length increases and the maximum fracture width decreases with the value of ∆σ. When the value of ∆σ increases from 0 MPa to 10 MPa, the fracture length increases by 53.3% and the fracture width decreases by 20.4%. This phenomenon is mainly due to the variation of ∆σ, which would affect the difficulty of fracture propagation, and then affect the fracture length and width. When the value

315

of ∆σ is increased, the rock initiation pressure is reduced. The fracture is more likely to extend in the direction of the maximum in-situ stress, and fracture length increases. In the case of large ∆σ, the fracture internal pressure is smaller. The fracturing fluid leaking into surrounding rock during the injection process is reduced so that the deformation of surrounding rock (especially rock mass II) is small. The reduction of rock deformation in Y direction means that the fracture width is reduced and the angle at which the

320

fracture deflects toward rock mass II is also reduced. Due to the same injection time and injection rate, the fluid loss coefficient is a certain value, and the amount of fracturing fluid in fractures is basically the same. The volume of two-dimensional model fracture is determined by the fracture length and width. Increasing in-situ stress difference, the fracture length increases and the width decreases, which is consistent with actual engineering.

20

Fig. 11. Relationship between fracture morphology and ∆σ values.

(a)

(b)

Fig. 12. Effects of ∆σ values on hydraulic fracturing results. (a) Relationship between fracture width and fracture length; (b) relationship between fracture length and ∆σ values, maximum fracture width and ∆σ values.

325

When the value of ∆σ is 10 MPa, the fracture effective part is the longest (Fig. 12(a)). After ∆σ is larger than 6 MPa, the fracture effective length does not change significantly with the increase of ∆σ. It can be seen from Tab. 1 that the rock initiation pressure increases with the value of ∆σ. In actual engineering, we hope that the rock initiation pressure is as small as possible to meet the economic requirements. In summary, when the value of ∆σ is 6 MPa, the fracture morphology is more accord with the requirements

21

330

of practical engineering. Although all factors have an effect on hydraulic fracturing, the degree is not the same. In the previous study, we proposed a sensitivity indicator to evaluate the impact factors [19]. However, this indicator has a limitation, which is only applicable to the case where the fracture length and width have the same trend. In this study, we optimized the solution formula of the sensitivity indicator to make it more applicable.

335

This indicator is obtained by the following formula.

QA =

  φf    φl ·φw        φf ·φw

(φl ≥ 1 ∪ φw ≥ 1) (φl ≥ 1 ∪ φw < 1)

φl

(9)

 φf ·φl   (φl < 1 ∪ φw ≥ 1)  φw       φf · φl · φw (φl < 1 ∪ φw ≥ 1) where QA is the sensitivity indicator of each factor; φf is the variation multiple of each factor; φl is the variation multiple of fracture length; φw is the variation multiple of fracture width. Analyzing the above hydraulic fracturing results, the sensitivity of each factor is shown in (Tab. 2). The smaller QA value is, the more obvious hydraulic fracturing is affected by this factor, and vice versa. 340

It can be seen from Tab. 2 that Qd < Q4σ . This situation means that d has a greater effect on fracture morphology than 4σ. According to the values of φl and φw , when the value of d is changed, the variation multiples of fracture length and width are larger. We can see from the sensitivity indicator that the fracture morphology is most sensitive to the variation of d values. The improved formula in this study is suitable for multi-factor sensitivity analysis, and the influence degree of various factors on hydraulic fracturing can

345

be studied by this formula. The factors with a large sensitivity indicator are optimized to obtain the ideal fractures that meet the heat transfer requirements. Table 2 Sensitivity analysis of impact factors

Factor

φf

φl

φw

QA

d

5.00

2.01

0.54

1.34

4σ

5.00

1.53

0.80

2.61

22

5. Conclusions Based on XFEM, a fluid-solid coupling model is established to study the hydraulic fracturing process with stratified rock. The model is used to study the factors affecting hydraulic fracturing results. By 350

considering the injection of proppant, the fracture effective part is used as an indicator to evaluate the fracture quality. We improved the formula proposed in our previous study to analyze the sensitivity of the factors affecting hydraulic fracturing. According to simulation results, the following conclusions can be obtained. (1) Rock stratification has a significant effect on the pore pressure distribution and fracture

355

morphology. Compared with isotropic rock, the fracture length in stratified rock decreased by 57.7% and the maximum fracture width increased by 110.8%. Hydraulic fractures in stratified rock deviate from the desired orthogonal path, which is not conducive to the formation of ideal fractures. However, it can be seen from the fracture evaluation indicator that the fracture effective length in stratified rock is 6.6 m, which is longer than that of the isotropic rock. Thus, hydraulic fractures in stratified rock have a larger

360

convective heat transfer area, which is beneficial for the extraction of geothermal energy. (2) As the value of d increases, the fracture morphology tends to be “narrow, long and straight”. When the value of d increases from 2 m to 10 m, the fracture length increases by 101.2% and the fracture width decreases by 46.1%. A large value of d reduces the rock pore pressure. The fracture effective length increases first and then decreases with the value of d. In actual engineering, engineers can determine the

365

perforation position by the optimal d value (6 m). (3) A large value of ∆σ would reduce the rock initiation pressure and increase the tendency of fractures to propagate along the ideal orthogonal path. The fracture effective length increases with the value of ∆σ. Considering the requirements of economy, when the value of ∆σ is 6 MPa, the fracture morphology is more in line with the needs of practical project. When the value of d is changed, the variation multiples

370

of fracture length and width are larger. We can see from the sensitivity indicator (QA ) that the fracture morphology is most sensitive to the variation of d values. We studied the initiation and propagation mechanism of hydraulic fractures in stratified rock, which provides a new perspective for other scholars to explore hydraulic fracturing of rocks with different

23

properties. Engineers can determine what kind of fractures should be obtained in practical project from 375

the fracture evaluation indicator (effective length) proposed in this study. Besides, the research results can guide the selection of perforation locations and prediction of fracture morphology during construction process.

Acknowledgement The work was supported by the National Natural Science Foundation of China (51606225).

380

References [1] Z.Feng, Y.Zhao, A.Zhou, N.Zhang. Development program of hot dry rock geothermal resource in the Yangbajing Basin of China. Renew. Energ., 39:490–495, 2012. [2] J.Haraden. The status of hot dry rock as an energy source. Energy, 17:777–786, 1992. [3] P.Asai, P.Panja, J.McLennan, J.Moore. Efficient workflow for simulation of multifractured enhanced

385

geothermal systems (EGS). Renew. Energ., 131:763–777, 2019. [4] Y.Shi, X.Song, G.Wang, J.McLennan, B.Forbes, X.Li, J.Li. Study on wellbore fluid flow and heat transfer of a multilateral-well CO2 enhanced geothermal system. Appl. Energ., 249:14–27, 2019. [5] M.Li, L.Zhang, G.Liu. Estimation of thermal properties of soil and backfilling material from thermal response tests (TRTs) for exploiting shallow geothermal energy: Sensitivity, identifiability, and

390

uncertainty. Renew. Energ., 132:1263–1270, 2019. [6] W.Zhang, Z.Qu, T.Guo, Z.Wang. Study of the enhanced geothermal system (EGS) heat mining from variably fractured hot dry rock under thermal stress. Renew. Energ., 143:855–871, 2019. [7] S.Li, X.Feng, D.Zhang, H.Tang.

Coupled thermo-hydro-mechanical analysis of stimulation and

production for fractured geothermal reservoirs. Appl. Energ., 247:40–59, 2019. 395

[8] Y.Zhang, Y.Ma, Z.Hu, H.Lei, L.Bai, Z.Lei, Q.Zhang.

An experimental investigation into the

characteristics of hydraulic fracturing and fracture permeability after hydraulic fracturing in granite. Renew. Energ., 140:615–624, 2019. 24

[9] C.Xu, P.Dowd, F.Zhao. A simplified coupled hydro-thermal model for enhanced geothermal systems. Appl. Energ., 140:135–145, 2015.

400

[10] T.Ingrid, M.Gutierrez. Coupled hydro-thermo-mechanical modeling of hydraulic fracturing in quasibrittle rocks using BPM-DEM. J. Rock Mech. Geotech. Eng., 9:92–104, 2017. [11] K.Olga, X.Weng. Numerical modeling of 3D hydraulic fractures interaction in complex naturally fractured formations. Rock Mech. Rock Eng., 51:1–19, 2018. [12] Y.Shi, X.Song, G.Wang, J.Li, L.Geng, X.Li. Numerical study on heat extraction performance of a

405

multilateral-well enhanced geothermal system considering complex hydraulic and natural fractures. Renew. Energ., 141:950–963, 2019. [13] G.E.Blandford, A.R.Ingraffea, J.A.Liggett. Two-dimensional stress intensity factor computations using the boundary element method. Int. J. Numer. Meth. Eng., 17:387–404, 1981. [14] D.Kumar, A.Ghassemi.

410

Three-dimensional poroelastic modeling of multiple hydraulic fracture

propagation from horizontal wells. Int. J. Rock Mech. Min., 105:192–209, 2018. [15] S.Salimzadeh, M.Grandahl, M.Medetbekova, H.Nick. A novel radial jet drilling stimulation technique for enhancing heat recovery from fractured geothermal reservoirs. Renew. Energ., 139:395–409, 2019. [16] F.Zhang, E.Dontsov. Modeling hydraulic fracture propagation and proppant transport in a two-layer formation with stress drop. Eng. Fract. Mech., 199:705–720, 2018.

415

[17] G.Amir, T.Jaber, S.Amin. The distinct element method (DEM) and the extended finite element method (XFEM) application for analysis of interaction between hydraulic and natural fractures. J. Petrol. Sci. Eng., 171:422–430, 2018. [18] Y.Feng, L.Chen, A.Suzuki, T.Kogawa, J.Okajima, A.Komiya, S.Maruyama.

Enhancement of

gas production from methane hydrate reservoirs by the combination of hydraulic fracturing and 420

depressurization method. Energ. Convers. Manage., 184:194–204, 2019. [19] W.Ma, Y.Wang, X.Wu, G.Liu. Hot dry rock (HDR) hydraulic fracturing propagation and impact factors assessment via sensitivity indicator. Renew. Energ., 146:2716–2723, 2019. 25

[20] Y.Feng, K.E.Gray. XFEM-based cohesive zone approach for modeling near-wellbore hydraulic fracture complexity. Acta Geotech., 14:1–26, 2018.

425

[21] B.He, X.Zhuang. Coupled discrete crack and porous media model for hydraulic fractures using the XFEM. KSCE J. Civ. Eng., 23:1017–1027, 2019. [22] R.Safari, A.Ghassemi. 3D thermo-poroelastic analysis of fracture network deformation and induced micro-seismicity in enhanced geothermal systems. Geothermics, 58:1–14, 2015. [23] Z.Sun, X.Zhang, Y.Xu, J.Yao, H.Wang, S.Lv, Z.Sun, Y.Huang, M.Cai, X.Huang.

430

Numerical

simulation of the heat extraction in EGS with thermal-hydraulic-mechanical coupling method based on discrete fractures model. Energy, 120:20–33, 2017. [24] F.Boukouvala, Y.Gao, F.Muzzio, M.G.Ierapetritou. Reduced-order discrete element method modeling. Chem. Eng. Sci., 95:12–26, 2013. [25] M.Aliyu, H.Chen. Optimum control parameters and long-term productivity of geothermal reservoirs

435

using coupled thermo-hydraulic process modelling. Renew. Energ., 112:151–165, 2017. [26] S.Cuvilliez, F.Feyel, E.Lorentz, M.Sylvie. A finite element approach coupling a continuous gradient damage model and a cohesive zone model within the framework of quasi-brittle failure. Comput. Method. Appl. M., 237:244–259, 2012. [27] B.Paul, M.Faivre, P.Massin, R.Giot, D.Colombo , F.Golfier, A.Martin. 3D coupled HM-XFEM

440

modeling with cohesive zone model and applications to non planar hydraulic fracture propagation and multiple hydraulic fractures interference. Comput. Method. Appl. M., 342:321–353, 2018. [28] M.Vahab, N.Khalili. X-FEM modeling of multizone hydraulic fracturing treatments within saturated porous media. Rock Mech. Rock Eng., 51:1–21, 2018. [29] F.Cruz, D.Roehl, E.Vargas. An XFEM element to model intersections between hydraulic and natural

445

fractures in porous rocks. Int. J. Rock Mech. Min., 112:385–397, 2018. [30] D.Kumar, A.Ghassemi. A three-dimensional analysis of simultaneous and sequential fracturing of horizontal wells. J. Petrol. Sci. Eng., 146:1006–1025, 2016. 26

[31] H.Vik, S.Salimzadeh, H.Nick. Heat recovery from multiple-fracture enhanced geothermal systems: The effect of thermoelastic fracture interactions. Renew. Energ., 121:606–622, 2018.

450

[32] S.Salimzadeh, N.Khalili. A three-phase XFEM model for hydraulic fracturing with cohesive crack propagation. Comput. Geotech., 69:82–92, 2015. [33] G.Liu, H.Pu, Z.Zhao, Y.Liu.

Coupled thermo-hydro-mechanical modeling on well pairs in

heterogeneous porous geothermal reservoirs. Energy, 171:631–653, 2019. [34] N.Barton, E.Quadros. Anisotropy is everywhere, to see, to measure, and to model. Rock Mech. Rock 455

Eng., 48:1323–1339, 2015. [35] J.Santos, J.Carcione. Finite-element harmonic experiments to model fractured induced anisotropy in poroelastic media. Comput. Method. Appl. M., 283:1189–1213, 2015. [36] M.Sheng, G.Li, D.Sutula, S.Tian, S.Bordas. XFEM modeling of multistage hydraulic fracturing in anisotropic shale formations. J. Petrol. Sci. Eng., 162:801–812, 2018.

460

[37] X.Wang, F.Shi, H.Liu, H.Wu. Numerical simulation of hydraulic fracturing in orthotropic formation based on the extended finite element method. J. Nat. Gas Sci. Eng., 33:56–69, 2016. [38] H.Eekelen. Hydraulic fracture geometry: Fracture containment in layered formations. Soc. Petrol. Eng. J., 19:127–127, 1982. [39] J.Guo, B.Luo, C.Lu, J.Lai, J.Ren. Numerical investigation of hydraulic fracture propagation in a

465

layered reservoir using the cohesive zone method. Eng. Fract. Mech., 186:195–207, 2017. [40] Y.Muceku, O.Korini, A.Kuriqi. Geotechnical analysis of Hill’s slopes areas in Heritage town of Berati, Albania. Period. Polytech-Civ., 60:61–73, 2016. [41] F.Wang, B.Li, Q.Chen, S.Zhang. Simulation of proppant distribution in hydraulically fractured shale network during shut-in periods. J. Petrol. Sci. and Eng., 178:467–474, 2019.

470

[42] L.Malvern. Introduction to the mechanics of a continuous medium. Lancet, 2:443, 1969. [43] G.K.Batchelor. An introduction to fluid dynamics. Phys. Today, 12:36–38, 1959. 27

[44] C.Liu, X.Wang, D.Deng, Y.Zhang, H.Liu. Optimal spacing of sequential and simultaneous fracturing in horizontal well. J. Nat. Gas Sci. Eng., 29:329–336, 2016. [45] X.Wang, C.Liu, H.Wang, H.Liu, H.Wu. Comparison of consecutive and alternate hydraulic fracturing 475

in horizontal wells using XFEM-based cohesive zone method. J. Nat. Gas Sci. Eng., 143:14–25, 2016. [46] M.Benzeggagh, M.Kenane.

Measurement of mixed-mode delamination fracture toughness of

unidirectional glass/epoxy composites with mixed-mode bending apparatus. Compos. Sci. Technol., 56:439–449, 1996. [47] T.Y.Belytschko, T.Black. Elastic crack growth in finite elements with minimal remeshing. Int. J. 480

Numer. Meth. Eng., 45:601–620, 1999. [48] C.Daux, N.Moes, J.Dolbow, N.Sukumar, T.Belytschko. Arbitrary branched and intersecting cracks with the extended finite element method. Int. J. Numer. Meth. Eng., 48:1741–1760, 2000. [49] T.Lv. Numerical simulation and experimental study of hydraulic fracturing in granite based on particle flow code. Jilin Univ., 2018.

485

[50] L.Guo, Y.Zhang, Y.Zhang, Z.Yu, J.Zhang. Experimental investigation of granite properties under different temperatures and pressures and numerical analysis of damage effect in enhanced geothermal system. Renew. Energ., 126:107–125, 2018. [51] B.Lecampion, A.Bunger, X.Zhang. Numerical methods for hydraulic fracture propagation: A review of recent trends. J. Nat. Gas Sci. Eng., 49:66–83, 2018.

490

[52] P.Wasantha, H.Konietzky, C.Xu. Effect of in-situ stress contrast on fracture containment during single- and multi-stage hydraulic fracturing. Eng. Fract. Mech., 205:175–189, 2019.

28

Highlights 1. Effective length is the length of the segment where fracture width greater than 1 mm 2. Rock stratification increases the fracture effective length 3. Relative position of perforation significantly affects fracture morphology 4. A large value of in-situ stress difference would reduce the rock initiation pressure 5. Perforation position is the most sensitive impact factor on fracture morphology

October 29, 2019

Conflict of interest statement We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work, there is no professional or other personal interest of any nature or kind in any product, service and/or company that could be construed as influencing the position presented in, or the review of, the manuscript entitled.

1