Numerical study on heat extraction performance of multistage fracturing Enhanced Geothermal System

Journal Pre-proof Numerical study on heat extraction performance of multistage fracturing Enhanced Geothermal System Songcai Han, Yuanfang Cheng, Qi Gao, Chuanliang Yan, Jincheng Zhang PII:

S0960-1481(19)31611-8

DOI:

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

Reference:

RENE 12483

To appear in:

Renewable Energy

Received Date: 30 June 2019 Revised Date:

23 September 2019

Accepted Date: 19 October 2019

Please cite this article as: Han S, Cheng Y, Gao Q, Yan C, Zhang J, Numerical study on heat extraction performance of multistage fracturing Enhanced Geothermal System, Renewable Energy (2019), doi: https://doi.org/10.1016/j.renene.2019.10.114. 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. © 2019 Published by Elsevier Ltd.

Numerical study on heat extraction performance of multistage fracturing Enhanced Geothermal System Songcai Hana, Yuanfang Chenga*, Qi Gaoa,b*, Chuanliang Yana, Jincheng Zhangc a

School of Petroleum Engineering, China University of Petroleum (East China), Qingdao, 266580, China School of Engineering, The University of Western Australia, 35 Stirling Highway, Perth, WA, 6009, Australia c Research Institute of Petroleum Engineering, Sinopec, Beijing, 100101,China *Corresponding authors, E-mail address： [email protected] (Y. Cheng); [email protected] (Q. Gao) b

Overburden layer

Overburden layer

Injection well with perforated casing completion

Production wells with openhole Stimulated completion HDR reservoir Stimulated reservoir volume (SRVs)

Injection well with perforated casing completion Induced microseismic events

HDR reservoir Hydraulic fractures

Underburden layer Multistage fracturing in injection well

Underburden layer Heat extraction from EGS

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

Abstract: A conceptual Enhanced Geothermal System (EGS) with triplet horizontal well layout is proposed, in which the injection well is flanked by two production wells. The injection well is completed with perforated casing, followed by multistage fracturing with casing packers. Two production wells are drilled through the stimulation zone induced by multistage fracturing and completed openhole to maximize contact with the reservoir. Based on a fully coupled hydrothermal (HT) model accounting for local thermal non-equilibrium (LTNE), sensitivity analysis and optimization for EGS design parameters are conducted. Simulation results show that increasing circulation rate only facilitates the heat extraction rate before thermal breakthrough, but cannot improve the heat extraction performance at the end of EGS operation. Compared with the fracture aperture, the increase in the stimulation zone permeability is more conducive to reducing the injection pressure and increasing the cumulative output thermal power. Insufficient or excessive fracturing stages are detrimental to economic heat extraction. It is found that more uniform flow into each perforation can achieve better thermal sweep efficiency and prevent local thermal breakthrough. The stimulation zone has an optimal reservoir permeability to match the fracture aperture. The optimum number of fracturing segments for a specified wellbore length is obtained. The preferable zone for the production well is between the fracture front and the outermost periphery of the stimulation zone. Keywords: EGS; horizontal wells; multistage fracturing; HT coupling; heat extraction performance; optimization

22

1. Introduction

23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41

Geothermal resources due to their renewability, cleanness and universality have attracted broad attention around the world [1-2]. Generally, the hot dry rock (HDR) reservoir possesses extra-low porosity and ultra-low permeability, which results in the difficulty of obtaining economic heat efficiency before stimulation treatment [2-3]. Consequently, permeability enhancement techniques play a critical role in these nearly impervious reservoirs [3]. Enhanced geothermal system (EGS) is widely regarded as an efficient method to extract heat from the HDR [1-3]. The key is to establish sufficient hydraulic channels for fluid circulation by a series of stimulation techniques, such as hydraulic fracturing, thermal shock, and chemical treatment [4-7]. Currently, most practical or conceptual EGS projects are designed as a single-stage fracturing in nearly vertical wells with openhole completion [5-6, 8-12]. On the one hand, this is mainly because horizontal drilling in hard and high-temperature rocks is not only technically challenging, but also extremely costly. On the other hand, EGS design focuses more on shear stimulation. It is generally believed that the shear slip of natural fractures can improve the hydraulic transmissivity of the overall reservoir due to the self-propping effect. However, mining experience has revealed limitations to the conventional EGS design [13-14]. An important uncertainty is whether the slipping fracture would be capable of self-propping under in-situ stresses. It has been found that only some natural fractures 1

42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59

in the stimulation zone have the strong conductivity, resulting in insufficient connectivity and injectivity of the reservoir. In addition, the flow generally tends to localize into highly conductive fractures around the wellbore, which will bring about low heat sweep efficiency. There were of course some exceptions, such as Paralana project and Salavatli project, which adopted the casing perforation completion technology [15-16]; Groβ Schönebeck project and Klaipėda project, which used the horizontal (or deviated) and multilateral drilling technology, respectively [17-18]; Groβ Schönebeck project, Newberry project and Qiabuqia project, which performed multiple stage stimulation using packers [17, 19-20]. The locations of these EGS sites are shown in Fig.1, and their characteristics are generalized in Table 1. More detail discussions about the feasibility of horizontal drilling and multistage fracturing in geothermal reservoirs can be found in Refs. [21-22]. In fact, EGS designs combining horizontal drilling with multistage fracturing may have many potential advantages: (1) creating more adequate stimulated reservoir volume (SRV), (2) improving the connectivity of the wellbore with HDR, (3) increasing the circulation rate through the system, (4) inhibiting thermal short circuiting and delaying thermal breakthrough, (5) achieving a better heat sweep efficiency, (6) inducing more thermal fractures perpendicular to transverse fractures.

60 61 62

Fig.1. Illustration of EGS sites distribution [8-9, 15-20] Table 1 Country Australia Turkey Lithuania

Germany

Description of EGS sites with different completion and stimulation methods [8-9, 15-20]

Project (Development period)

Paralana (2005-present) Salavatli (2003-present) Klaipėda (1996-present) Groβ Schönebeck (2007-present)

Well name

Depth (m)

Temperature (℃)

Well type

Completion method

Stimulation method

Paralana-2

4003

180~200

deviated well

perforated casing

multistage fracturing

AS-9

2500

160~173

perforated casing

acid treatment

-

>1200

50~80

vertical well multilateral well

open hole

radial jet

zonal isolation materials

multistage fracturing + an acid treatment

GrSk-4

4400

>150

2

deviated well

63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90

USA

Newberry (2009-present)

NWG 55-29

3067

331

China

Qiabuqia (2007-present)

GR-1

3600

>180

nearly horizontal well horizontal well

zonal isolation materials perforated casing

multistage fracturing multistage fracturing

Evaluating the heat extraction performance of an EGS involves precisely solving a hydraulic-thermal (HT) coupling process. Li et al [21] used an analytical model to perform a sensitivity analysis and optimization for an EGS with two horizontal wells and multiple fracturing stages. Still, numerical simulation is the most efficient way to analyze the heat extraction from an EGS. Based on regional geological models and field data, Hofmann et al [23-24] conducted some numerical studies to identify the production potential of an EGS with difference well and fracture arrangements. Up to now, different conceptual EGS models have been designed, modeled, and evaluated in many works, such as doublet vertical well layout [10-11], doublet horizontal well layout in the vertical direction [25-26], triplet vertical well layout [12, 27], triplet horizontal well layout in the horizontal direction [20, 28], quintuplet vertical well layout [12, 29], and multilateral well layout in the vertical direction [30-31]. In addition, the adaptability of different circulating fluids (such as water and supercritical CO2) in EGS has been investigated in some studies [28, 32]. These efforts can provide insights on how to optimize the long-term heat extraction performance of an EGS. However, in many models, it was assumed that the local thermal equilibrium (LTE) between the rock matrix and the pore fluid could be achieved instantaneously. Many studies have demonstrated that the local thermal non-equilibrium (LTNE) theory is more appropriate if there is a rapid heat transfer or large temperature difference between the rock matrix and the flowing fluid [33-35]. In addition, the effect of temperature change on pore pressure and fluid properties should also be incorporated into the fully coupled HT model. The purpose of this paper is to investigate the heat extraction performance of a novel EGS with triplet horizontal well layout (one injection well with multistage fracturing and two production wells with openhole completion), as shown in Fig.2. With a sensitivity study of the heat extraction performance to main design parameters, these results can provide an optimal strategy for how to achieve better heat extraction performance for an EGS combining horizontal wells with multistage fractures.

3

91 92

Fig. 2.

Schematic of heat extraction from an EGS with horizontal wells and multistage fracturing

93

2. Mathematical model

94

2.1. HT coupling equations

95 96 97 98 99 100 101 102 103 104 105 106 107 108 109

Simulating heat extraction from a stimulated geothermal reservoir has to account for the dominant physical phenomena including fluid flow and heat transfer in fractures and porous media. The process involves a complex HT coupling problem that usually gives rise to formidable challenges in developing a proper mathematical model. Some assumptions have to be adopted to simplify the complex process. (1) Multistage fracturing in the injection well can create equidistant and equal-sized permeability enhancement zones (namely SRV) that are enclosed by the unstimulated HDR. (2) The SRV and HDR are simplified as the equivalent continuous porous media with homogeneous and isotropic hydrothermal properties, but the permeability of the SRV is much greater than that of the HDR. (3) The reservoir is fully saturated with water, and the pore water does not vaporize throughout the simulation. (4) Fluid flow in porous media and fractures obeys Darcy’s law and Cubic law, respectively. (5) The heat transfer is dominated by conduction and convection processes, and local thermal non-equilibrium is considered between the solid phase and the fluid phase. Mass conservation equations in porous media and fractures [36-40]:

110

ρl C1m

111

ρl C1 f eh

112

ρk  ∂p ∂T − ρl C2 m = ∇ ⋅  l m (∇p + ρl g ∇z )  ∂t ∂t  µ 

(1)

 eh ρl k f  ∂p ∂T − ρl C2 f eh = ∇T ⋅  (∇ T p + ρl g ∇ T z )  + n ⋅ Qm ∂t ∂t  µ 

(2)

Energy balance equations in porous media and fractures [36-40]:

113

(1 − φm ) ρm Cm

114

φm ρl C p,l

∂Tm = ∇ ⋅ ( (1 − φm )λm∇Tm ) + qml (Tl − Tm ) ∂t

∂Tl + φm ρl C p,l um ⋅ ∇Tl = ∇ ⋅ (φm λl ∇Tl ) + qml (Tm − Tl ) ∂t 4

(3) (4)

115

(1 − φ f )eh ρ f C f

∂Tm = ∇T ⋅ ( (1 − φ f )eh λ f ∇TTm ) + eh q fl (Tl − Tm ) + n ⋅ qm ∂t

(5)

117

∂Tl + φ f eh ρl C p ,l u f ⋅ ∇TTl = ∇T ⋅ (φ f eh λl ∇TTl ) + eh q fl (Tm − Tl ) + n ⋅ ql (6) ∂t where p , T l and Tm are the pore pressure, fluid temperature, and rock matrix

118 119 120

temperature, respectively; C1m = φm Sl + (1 − φm ) Sm represents the effective storage coefficient of porous media; C2 m = α m β m + φm ( βl − β m ) represents the effective thermal expansion coefficient of porous media; C1 f = φ f Sl + (1 − φ f ) S f represents the effective

121

storage coefficient the fracture; C2 f = α f β f + φm ( β l − β f ) represents the effective

122

thermal expansion coefficient of the fracture; Sl , Sm and S f represent the storage

123 124 125

coefficients of the fluid, rock matrix and fracture, respectively; β l , β m and β f represent the thermal expansion coefficients of the fluid, rock matrix and fracture, respectively; ρl , ρ m and ρ f are the density of the fluid, rock matrix, and fracture,

126

respectively; C p ,l , Cm and C f are the specific heat capacity of the fluid, rock matrix

127 128

and fracture, respectively; λl , λm and λ f are the heat conductivity of the fluid, rock matrix and fracture, respectively; α m and α f represent the Biot’s coefficients of

129

porous media and fractures, respectively; φm and φ f are the porosity of porous

130 131 132 133

media and fractures, respectively; k m and k f are the permeability of porous media and fractures, respectively; µ is the fluid dynamic viscosity; g is the gravitational acceleration; z represents elevation in the z direction; qml represents the matrix-fluid interface heat transfer coefficient; q fl represents the fracture-fluid interface heat

134 135 136 137 138 139 140 141

transfer coefficient. k f = eh2 12 f , where eh is the hydraulic aperture between fracture surfaces; f is the roughness factor, characterizing the tortuosity and roughness of the fracture. n ⋅ Q m = n ⋅ ( − ρ l k m µ ∇ p ) is the mass flux exchange between porous media and fractures. n ⋅ q m = n ⋅ ( − (1 − φ m ) λ m ∇ Tm ) represents the heat flux exchange between the rock matrix and the fracture. n ⋅ ql = n ⋅ ( −φm λl ∇ Tl ) denotes the heat flux exchange between the pore fluid and the fracture fluid. um = − km µ (∇p + ρl g∇z ) is the Darcy’s rate in porous media. u f = − k f µ (∇T p + ρl g∇T z ) is the Darcy’s rate in fractures. ∇ T denotes the gradient operator restricted to the fracture’s tangential plane.

142

2.2. Temperature-dependent thermodynamic properties

143 144 145 146 147

Temperature-dependent fluid thermodynamic properties are also incorporated into the aforementioned mass and energy conservation equations to account for the bidirectional HT coupling. The changes in fluid viscosity, specific heat capacity, density, heat conductivity and thermal expansion coefficient with the temperature are described by some empirical formulas as follows [28]

116

φ f eh ρl C p,l

5

148

 µl = 1.3799 − 0.0212Tl + 1.3604 × 10−4 Tl 2 − 4.6454 × 10−7 Tl 3 + 8.9043  × 10−10 Tl 4 − 9.0791× 10−13 Tl 5 + 3.8457 × 10−16 Tl 6 , Tl ∈[273K, 413K]   −5 −8 2 −11 3  µl = 0.004 − 2.1075 × 10 Tl + 3.8577 × 10 Tl − 2.3973 × 10 Tl , Tl ∈ [413K,553K]  2 −4 3 −7 4 C p,l = 12010 − 80.4Tl + 0.3Tl − 5.4 × 10 Tl + 3.6 × 10 Tl , Tl ∈[273K,553K] (7)  2 −7 3  ρl = 838.4661 + 1.4005Tl − 0.003Tl − 3.7182 × 10 Tl , Tl ∈ [273K,553K] λ = −0.8691 + 0.0089T − 1.5837 × 10−5 T 2 + 7.9754 × 10−9 T 3 , T ∈[273,553] l l l l  l  βl = − (1 ρl ) ⋅ ( d ρl dTl ) , Tl ∈ [273K,553K]

149

3. Numerical method and Validation

150

3.1. Numerical implementation

151 152 153 154 155 156 157 158 159 160 161 162

The governing equations describing the heat extraction process, Eqs.(1)-(7), combined with the corresponding initial and boundary conditions can be solved by the finite element method. COMSOL Multiphysics, which is a well-known commercial tool, is adopted to solve the partial differential equations. To accelerate convergence, the segregated approach, which solves each physics sequentially until convergence, is adopted to solve the HT coupling model. Although more iterations may be required for the same problem, each loop through the segregated solution approach can be much faster than the Newton-Raphson step required for the fully coupled approach. It is also noticeably profitable to model all fractures as interior boundaries, because it eliminates the need to create a geometry with a large aspect ratio. Otherwise, the very long and narrow fracture domain would bring about a dense mesh system consisting of great number of tiny elements.

163

3.2. Validation of heat extraction

164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179

It should be noted that there are currently no on-site or experimental data available to validate the numerical model proposed in this paper. The applicability of the hydrothermal simulation method has been confirmed in our fully coupled thermal-hydraulic-mechanical model [37]. To further testify the validity of the heat extraction model considering LTNE theory, a simple two-dimensional geometry model embedded a smooth fracture with constant aperture is used to model the heat exchange between the fracture wall and flowing fluid (as shown in Fig.3). The solid temperatures on the up and down outer boundaries are maintained constant (T0), however, the left and right boundaries are thermally insulated. The initial temperature of the fluid in the fracture is identical with the initial temperature of the solid (Tf0= Ts0= T0). It is assumed that the fluid injected with a constant rate flows parallel through the impermeable fracture walls. Temperature-dependent fluid properties are ignored. A constant heat transfer coefficient is assigned along the fracture plane. The low temperature injection fluid is gradually heated to a high temperature fluid by heat convection. Another reasonable assumption for the thin fracture is that a uniform fluid temperature crosses the fracture section. The parameters used in this simulation are 6

180

181 182 183

listed in Table.2.

Fig. 3.

A conceptual diagram of heated flow-through simulation in a rock fracture

Table 2

Input parameters for heat transfer validation in a fracture [35]

Parameters

Symbol

Unit

Value

Initial temperature

T0

K

363

Injection fluid temperature

Tin

K

315

Injection velocity

u

mm/s

10.63

Fluid density

ρl

3

kg/m

1000

Fluid viscosity

µl

Pa·s

0.001

Cp,l

J/(kg·K)

4200

Fluid heat conductivity

λl

J/(m·K·s)

0.6

Rock heat conductivity

λm

J/(m·K·s)

3.5

Rock heat capacity

Cp,m

J/(kg·K)

790

Hydraulic aperture

eh

µm

19.17

Fluid heat capacity

184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200

There exists an analytical solution for heat convection in a single fracture. The fluid temperature in the fracture can be obtained as follows [41] T f ( x ) = T0 + (Tin − T0 ) exp( − x

4 h λm ) u ρ l cl eh (2 λ m + hH )

(8)

where x is the distance from the inlet end to the outlet end; h is the fracture heat transfer coefficient; H is the height of the model. Meanwhile, the change in the solid temperature due to heat conduction satisfies [41] Tm ( x , y ) = T0 +

2h H (T0 − T f ( x ))( y − ) 2 λm + hH 2

(9)

where x and y denote the coordinates in the two-dimensional plane. It is noteworthy that Eq.9 is derived based on one-dimensional heat conduction, perpendicular to the fracture plane. Therefore, there may exist a slight difference between the analytical solution and the numerical solution. Fig.4 demonstrates a comparison between the analytical solutions and the numerical solutions for the fluid temperature in the fracture and the matrix temperature along line AB. It can be clearly seen that the remarkably good agreements for different heat transfer coefficients demonstrate the accuracy of the numerical algorithm. As the heat transfer coefficient increases, the fluid temperature rises faster 7

201 202 203

and the matrix temperature drops faster. This result indicates that the heat convection process should be taken into account when investigating the heat extraction of the fractured reservoir. 363 360

Matrix temperature (K)

Fluid temperature (K)

362 350

340

330

Analytical solution Numerical solution

h=7 W/(m2·K) h=12 W/(m2·K) h=25 W/(m2·K) h=50 W/(m2·K)

320

310 0

20

40

60

80

361 360 359 Analytical solution Numerical solution h=7 W/(m2·K) h=12 W/(m2·K) h=25 W/(m2·K) h=50 W/(m2·K)

358 357 356 0

100

Distance from inlet piont (mm)

5 10 15 20 Distance perpendicular to the fracture (mm)

25

204 205 206

Fig. 4.

207

4. A conceptual EGS model

208

4.1. Model description

209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233

The idealized EGS consists of a horizontal injection well and two horizontal production wells (see Fig.2). The injection well is completed with perforated casing. Sequentially, multistage fracturing in the injection well is performed by casing packers. Transverse fractures are evenly generated from the toe to heel of the wellbore. Meanwhile, a permeability enhancement zone is also induced around each main fracture. Two production wells are located on both sides of the injection well, and their positions could be determined based on microseismic observations. The production wells are completed openhole to maximize contact with the reservoir. During fluid circulation, the cold fluid is only injected through perforations. Conversely, the heated fluid can be withdrawn through the entire production well. The two-dimensional schematic of the well and fracture arrangement is demonstrated in Fig.5(a). In this study, the computational domain with a size of 1600 m×1200 m×600 m cube is located at a depth of 3900-4500 m, which contains unstimulated HDR, SRVs, three horizontal wells and transverse fractures. The heat extraction zone of interest is a 910 m×800 m×200 m hexahedral volume and located at the center of the computational domain. These horizontal wells are 820 m in length, and 300 m away from the top of the computational domain. To facilitate fracturing implementation, the first and last fractures are offset ∆xe/2(=20 m) from the toe and heel of the wellbore, respectively. It is assumed that the aperture of transverse fractures complies with the PKN model [42-43] (Fig.6), and their length and height are 700 m and 200 m, respectively. These artificial transverse fractures are the primary conduits for fluid circulation and heat extraction. The scope of the induced SRV is closely related to fracturing parameters and reservoir properties, which can also be assessed by microseismic events. It is assumed that each SRV has the same size of 130 m×800 m×200 m.

Comparison of numerical solutions with analytical solutions: (a) fluid temperature in fracture; (b) matrix

temperature along line AB.

8

234 235 236 237 238 239 240 241 242 243 244 245

In the basic case, the spacing between injection and production wells is equal to the half-length of the fracture, which indicates that the production well is fully connected to the injection well. Insufficient fracturing stages (e.g., Nf=5 in the basic case) can pose inadequate stimulation. It is also noteworthy that excessive fracturing stages will result in overlapping stimulation zones. Considering the symmetry of the model, only half of the computational domain is simulated, as shown Fig.5(b). Numerical experiments have been conducted to investigate the effect of the number of the mesh on the average production temperature. Accuracy indicates that the change in production temperature is less than 1% from a mesh of about 380,000 elements to one with about 700,000. However, an increase in the number of elements will cause a significant increase in computation time. Therefore, the mesh system with about 380,000 tetrahedral elements is adopted in this study, as shown Fig.5(c). lh = ∆xe + ( N f − 1)∆ye

246

247 248 249

250 251

Fig. 5.

Computional model and mesh system. (a) 2-D schematic of the arrangement of wells and fractures;

(b) half of the geometry model of a conceptual EGS; (c) the mesh scheme of the simulation zone.

Fig. 6.

Distribution of the fracture aperture 9

252

4.2. Initial and boundary conditions

253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271

A 40-year production period is designated for simulating heat extraction in the conceptual EGS. (1) The initial conditions for the pore pressure and temperature fields are prescribed as follows: p0 with a hydraulic gradient of 0.85MPa/100m and a ground pressure of 0MPa; T0 with a thermal gradient of 4K/100m and a ground temperature of 300K. (2) For infinite peripheral boundaries, the pore pressure and temperature obey the Dirichlet boundary condition with pout=p0 and Tout=T0, respectively. (3) The top and bottom boundaries are considered impervious and thermally insulated. (4) For internal boundary conditions, the total injection rate is assigned as Qtot ( rw , t ) = 160 kg s . The rate flowing into each transverse fracture can be adjusted by the perforation pressure drop to achieve a good heat sweep efficiency (a strategy called “limited entry completion” [44]). For the basic case, the total injection rate is evenly allocated to each perforation stage, namely Qi ( rw , t ) = Qtot (2 N f )=16kg s . Meanwhile, a point heat source can be defined as follows based on the injection rate and injection temperature, q ( rw , t ) = C p , l Qi ( rw , t ) (Tin − T (t ))W m , where the injection temperature Tin=303 K, and T(t) represents the current temperature at the injection point. The production wells extract hot fluid from the geothermal reservoir at a constant pressure of 24 MPa. It should be noted that all artificial fractures within the domain are regarded as internal boundaries, implicitly considering the mass and energy exchange between porous media and fractures.

272

4.3. Basic input parameters

273 274 275 276 277 278

Most of physical parameters required for the fully coupled HF model are collected from Refs. [10, 28, 31, 37-39] and listed in Table.3. All parameters used in this simulation are reasonable and consistent with the physical properties of typical geothermal reservoirs. Temperature-dependent fluid thermodynamic properties are described by some empirical formulas (see Eq. 7). Table 3

Input parameters for HF coupling analysis of an EGS

Parameters

Symbol

Unit

Value 0.8 [39]

Biot’s coefficient of unstimulated zone

αm

-

Matrix porosity of unstimulated zone

φm

-

0.06 [31]

Matrix permeability of unstimulated zone

km0

m

1×10-17 [31]

Matrix density of unstimulated zone

ρm

kg/m3

2700 [31]

Cp,m

J/(kg·K)

1000 [31]

Matrix thermal diffusivity of unstimulated zone

λm

W/(m·K)

2.8 [28, 31]

Matrix compressibility of unstimulated zone

Sm

1/MPa

8.4×10-3 [37-39]

Thermal expansion coefficient of unstimulated zone

βm

1/K

3.2×10-6 [37-39]

Biot’s coefficient of stimulated zone

αm

-

0.8 [39]

Matrix porosity of stimulated zone

φm

-

Matrix heat capacity of unstimulated zone

2

0.15 [31]

Matrix permeability of stimulated zone

km0

m

1×10-15 [31]

Matrix density of stimulated zone

ρm

kg/m3

2600 [31]

10

2

Cp,m

J/(kg·K)

1000 [28]

Matrix thermal diffusivity of stimulated zone

λm

W/(m·K)

2.8 [31]

Matrix compressibility of stimulated zone

Sm

1/MPa

8.4×10-3 [37-39]

Thermal expansion coefficient of stimulated zone

βm

1/K

3.2×10-6 [37-39]

Porous media heat transfer coefficient

qml

W/(m3·K)

1.0 [10]

Biot’s coefficient of fracture

α

-

0.8 [39]

Matrix heat capacity of stimulated zone

f

3

Fracture density

ρ

kg/m

2000 [31]

Fracture porosity

φf

-

0.95 [31]

Fracture heat capacity

Cp,f

J/(kg·K)

850 [28, 31]

Fracture thermal diffusivity

λf

W/(m·K)

1.8 [28, 31]

Fracture bulk compressibility

Sf

1/MPa

6.4×10-3 [37]

Fracture thermal expansion coefficient

βf

1/K

5.0×10-6 [37]

Fracture surface heat transfer coefficient

qfl

W/(m3·K)

1.0 [10]

Fracture aperture at injection point

e0

cm

0.7 [28]

Fracture friction factor

f

-

15

279

5. Simulation results and analysis

280

5.1. Definition of evaluation parameters

281 282 283 284 285 286 287 288

To optimize the heat extraction performance, some evaluation parameters are defined to characterize the long-term heat extraction process. It is noteworthy that only half of the numerical model is simulated due to its symmetry. Therefore, cumulative evaluation parameters should be multiplied by 2 to account for this entire model of interest. In the LTNE theory, there is a temperature difference between the solid phase and the fluid phase. The equivalent temperature Teq is defined to account for both the temperature of the fluid and the solid, as follows

289 290 291 292 293 294 295 296 297 298

Teqi =

(1 − φi )ρi Cp,iTi + φi ρl Cp,lTl (1 − φi )ρi Cp,i + φi ρl Cp,l

(10)

where i=m and f represent porous media and fractures, respectively. The average production temperature is the line average of fluid temperature along the production wellbore Tav =

∫

lh

Tl (t ) dl lh

(11)

where lh represents the length of the production wellbore. The local heat extraction ratio ηl is defined to reveal the local heat transfer efficiency of the EGS [11, 37] ηl =

T0 − Tm ( t ) T0 − Tinj

(12)

where T0 is the initial temperature of the reservoir; Tm (t ) is the reservoir temperature 11

299

at time t; Tinj is the temperature of the injection fluid.

300 301

The total heat extraction ratio ηt is defined as the cumulative heat extracted by the heat transfer divided by the total heat stored in the EGS [31, 37] 2 ∫ ρ m C p , m (T0 − Tm (t )) dV

302

ηt =

Vs

2 ∫ ρ m C p , m (T0 − Tin )dV

(13)

Vs

303 304 305 306 307 308 309 310 311 312

where Vs is equal to half of the heat extraction zone of interest. The output thermal power P is defined as the heat extraction rate from the reservoir P = 2∫ Qout C p,l (Tl − Tin )dl l

(14)

where Qout represents the output rate of each production well. According to Ref. [26], it is assumed that no further thermal power can be produced when the service time of the EGS is arrived. The service time ts will be defined in the following. The cumulative output thermal power W is the total extraction heat from the reservoir during the service time of the EGS t (15) W = 2 ∫ ∫ Qout C p ,l (Tl − Tin ) dldt s

0

lh

313 314 315 316 317 318 319 320 321

Thermal breakthrough time (tb) is defined as a critical time, after which the temperature of the production fluid begins to decline obviously. The service time or lifetime of the EGS (ts) is defined as the operation time, at which the production fluid temperature reaches the abandonment temperature. Currently, different criteria for abandonment temperature have been adopted in some literature. These criteria generally range from 5% to 15% decline of the production temperature compared with its maximum value [10-11]. According to Ref.[11], a dimensionless temperature is defined to characterize the production temperature drop θ (t ) =

Tav ,max − Tav (ts ) Tav ,max − Tin

(16)

322 323

where Tav ,max represents the maximum temperature of the production fluid. θ=10% is adopted in this study to determine the service time of the EGS.

324

5.2. Evolution of pressure and temperature

325 326 327 328 329 330 331 332 333 334

Fig.7(a) illustrates the temporal and spatial evolution of pore pressure and flow rate. It can be seen that the pore pressure gradually increase with time. This is because the hot fluid occupying the pore space is gradually replaced by the injected cold fluid. The increase in fluid viscosity will result in a continuous increase in the seepage resistance. As the injected fluid diffuses around transverse fractures, the low permeable reservoir is insufficient to accommodate the continuous injection fluid, which will induce pressurized zones around transverse fractures. The pore pressure in inner zones is greater than that on both sides of the simulated zone. This is mainly because the seepage interference between inner fractures results in a higher flow resistance. Therefore, the fluid flow velocity in outer cracks is slightly larger than that 12

335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353

in internal cracks. It is taken for granted that the overall reservoir will eventually reach a steady state with the circulation of the injected fluid. Fig.7(b) demonstrates the evolution of the equivalent temperature. It is evident that transverse fractures are the main pathways for heat transport. With the continuous injection of the cold fluid, the low-temperature zones generated by heat conduction and heat convection gradually expand toward the length and height directions of the injection points. It can be found that these cooled zones enclosing transverse fractures are somewhat non-uniformly even at the same injection rates. The heat transfer in outer transverse fractures is slightly faster that in inner transverse fractures. These results further reveal that the seepage interference among inner fractures has a certain effect on the heat extraction. Fig.7(c) illustrate the evolution of local heat extraction ratio of the reservoir. It can be seen that the heat extraction profile advances in an approximately spindle shape. This is due to the opening of the transverse fracture conforms to the morphology of the PKN model. The local heat extraction ratio in each SRV is slightly different. The closer to both ends of the wellbore, the higher the heat extraction ratio of the SRV enclosing the fracture. Overall, the total heat extraction ratio of the reservoir is relatively low for this case (See Fig.8(a)). Therefore, the number of fracturing stages needs to be optimized. (a) Pore pressure and fow rate

(b) Reservoir equivalent temperature

(c) Reservoir local heat extraction ratio

t=1a

t=10a

t=20a

t=40a

m/s K MPa

354 13

355 356 357 358

Fig. 7. (a) Temporal and spatial evolution of pore pressure and flow velocity. The arrows indicate the direction of fluid flow, and the color of the arrow represent the relative magnitude of the velocity. (b) Temporal and spatial evolution of the reservoir equivalent temperature. (c) Temporal and spatial evolution of the local heat extraction ratio.

359

5.3 Optimization of heat extract efficiency

360

5.3.1 Effect of injection rate on heat extraction

361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389

Fig.8(a) illustrates the change in production temperature and heat extraction ratio with production time at different injection rates. It should be noted that the injection rate here is equal to half of the total injection rate due to symmetry. As the injection rate increases, the production temperature declines faster, which results in an earlier thermal breakthrough and a shorter service time. The service times of the EGS are 32a, 20a, 13.8a and 10.3a at different injection rates, e.g., 60kg/(m·s), 80kg/(m·s), 100kg/(m·s) and 120kg/(m·s), respectively. Although the heat extraction ratio increases faster at higher injection rates, the final heat extraction ratio is approximately equal (about 0.12) when arriving at the abandonment temperature. The result reveals that the heat extraction ratio of the overall reservoir cannot be improved by changing injection rate. Fig.8(b) illustrates the change in output thermal power and cumulative output thermal power. It can be seen that the heat extraction process shows three distinct stages during the EGS operation time. The duration of the first stage is about 1 year, during which the output thermal power drops sharply, and the corresponding cumulative output thermal power increases rapidly. This is due to the rapid decline in early productivity at constant production pressure. There is no obvious difference in the first stage at different injection rates. The second stage is the stable production stage or the ascending stage, depending on to the injection rate. The duration of this stage is closely related to the thermal breakthrough time of the EGS (e.g., 14a, 11a, 7a, 5a and 4a for five injection rates, respectively). After the second stage, the output thermal power drops at a relatively fast rate. The greater the injection rate, the faster the declining rate. The changes in production temperature and output thermal power versus EGS operation time are qualitatively similar with some Refs.[10-12, 31]. When arriving at the service time of the EGS, the cumulative output thermal powers are 3.27×1016J, 3.17×1016J, 2.55×1016J, 2.16×1016J and 1.91×1016J, respectively. Therefore, increasing the injection rate only facilitates the heat extraction rate before thermal breakthrough. It can be inferred that the optimal injection rate is 60kg/(m·s) in this case.

14

0.6

Q=50 kg/(m·s) Q=60 kg/(m·s) Q=80 kg/(m·s) Q=100kg/(m·s) Q=120kg/(m·s)

420 410 400 390 380

0.4

0.2

(a)

370

0.0

360 0

5

10

15

20

25

30

35

40

Heat extraction ratio

Temperature (K)

430

Output thermal power (MW)

0.8

440

60 50

3.5x1016 3.0x1016 2.5x1016

40

2.0x1016 30

1.5x1016

20

1.0x1016

(b)

10

5.0x1015 0.0

0 0

Time (a)

Fig. 8.

4.0x1016

Q=50kg/(m·s) Q=60 kg/(m·s) Q=80 kg/(m·s) Q=100kg/(m·s) Q=120kg/(m·s)

70

450

390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406

4.5x1016

80

460

5

10

15

20

25

30

35

Cumulative thermal extraction power (J)

1.0

470

40

Time (a)

Effect of the injection rate on heat extraction. (a) Production temperature and total heat

extraction ratio, (b) output thermal power and cumulative output thermal power.

Fig.9 illustrates the effect of the injection rate distribution on the heat extraction process. It can be seen that although the injection rate is 80kg/s for both cases, the local heat extraction ratio is completely different. The more uniformly distributed injection rate contributes to evenly cooling the heat rock mass (Fig.9(b)). Conversely, the non-uniform injection mode will result in earlier local thermal breakthrough and lower local heat extraction ratio (Fig.9(a)). Generally, due to the flow resistance along the horizontal wellbore, more fluid is allocated to the perforation section near the heel of the horizontal wellbore, where the thermal breakthrough typically occurs earlier. In order to achieve better thermal sweep efficiency and prevent thermal short circuiting during long-term circulation, it would desirable for flow to be as uniform as possible across each stage. Referring to a strategy called “limited entry completion”, it is possible to intentionally create perforations with relatively small diameters to increase the perforation pressure drop and promote a more uniform flow into the perforations [44].

407 408 409

Fig. 9.

Effect of the injection rate on the local heat extraction. (a) Linearly distributed injection

rate, (b) uniformly distributed injection rate

15

410

5.3.2 Effect of hydraulic transmissivity on heat extraction

411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431

Fig.10(a) illustrates the effect of the maximum fracture aperture at the injection point on the production temperature and total heat extraction ratio. It can be found that as the fracture aperture increases, the decline in production temperature becomes slower, and the increase in the heat extraction ratio is slower. This is because the smaller fracture aperture will cause a higher injection pressure at the same injection rate, resulting in a higher flow rate in the narrow fracture channel. It can be seen from Fig.10(b) that the increase in the fracture aperture contributes to reducing the injection pressure and increasing the cumulative output thermal power. When the fracture aperture increases by 900% (from 0.1cm to 1cm), the cumulative output thermal power increases by 203.5%, and the injection pressure decreases by 46.2%. It can be concluded from the research that increasing or maintaining the fracture aperture by proppants can improve the heat extraction performance. Mining experience has shown that proppants used in EGS have consistently led to improved productivity and injectivity [14, 17, 45]. In fact, there has been controversy about whether to use proppants in geothermal reservoirs. A skepticism is whether newly created fractures would be capable of self-propping. Some experiments and practices have demonstrated that relying solely on the self-propping effect of fractures cannot maintain high conductivity in long-term production [46-48]. Another potential concern is that using proppants may bring about the thermal short-circuit pathway. If a thermal short circuit develops, it could be easily mitigated by the zonal isolation technology [22]. 80

Temperature (K)

440 430

0.6

420 Fracture aperture at injection point 410

0.1 cm 0.4 cm 0.7 cm 1 cm

400 390 380 370

0.4

5

10

15

2.0x1016

60

1.5x1016 50 1.0x1016

20

25

30

35

(b)

0.0

40

5.0x1015 0.2 0.4 0.6 0.8 1.0 Fracture aperture at injection point (cm)

Time (a)

Fig. 10. Effect of fracture aperture on heat extraction. (a) Production temperature and total heat 1.0

Temperature (K)

450 440

Permeability in stimulation zones -16

0.6

2

2×10 m 6×10-16 m2 1×10-15 m2 2×10-15 m2 3×10-15 m2

430 420 410 400 390

0.4

Heat extraction ratio

0.8

Injection pressure (MPa)

460

0.2

(a) 10

15

20

25

30

35

90

4.5x1016 4.0x1016

80

40

Time (a)

3.5x1016

70

3.0x1016

60

2.5x1016 2.0x1016

50

30 0.0

0.0 5

5.0x1016

40

380 0

100

(b) 6.0x10-16 1.2x10-15 1.8x10-15 2.4x10-15 3.0x10-15 Permeability of stimulated zone (m2)

1.5x1016 1.0x1016

Cumulative output thermal power (J)

extraction ratio, (b) injection pressure and cumulative output thermal power. 470

435 436 437

2.5x1016

70

40 0.0

350 0

3.0x1016

0.2

(a)

360

432 433 434

Heat extraction ratio

0.8

450

Injection pressure (MPa)

460

Cumulative output thermal power (J)

3.5x1016

1.0

470

5.0x1015

Fig. 11. Effect of the permeability of the stimulation zone on heat extraction. (a) Production temperature and total heat extraction ratio, (b) injection pressure and cumulative output thermal 16

438 439 440 441 442 443 444 445 446 447 448 449 450 451

Fig.11(a) illustrates the effect of the permeability of the stimulation zone on the production temperature and total heat extraction ratio. As the permeability increases, the decline in production temperature becomes slower, which contributes to prolong the service time of the EGS. Although the total heat extraction ratio of the low permeable reservoir is slightly higher, this will result in a very high injection pressure to achieve the requirement of a constant injection rate (Fig.11(b)). It can be inferred from Fig.11(b) that there may be an optimal reservoir permeability to match the specified fracture aperture (e.g., in this case, km=2×1015m2 for 0.7m of fracture aperture). When the stimulation zone permeability increases by 900% (from 2×1016m2 to 2×1015m2), the cumulative output thermal power increases by 519.0%, and the injection pressure decreases by 58.6%. Therefore, in multistage fracturing designs, more attention should be paid to shear stimulation to enhance the permeability of the HDR reservoir enclosing transverse fractures.

452

5.3.3 Effect of fracture number on heat extraction

453 454 455 456 457 458 459 460 461 462 463

Fig.12(a) illustrates the effect of the number of fractures (Nf) on the production temperature and total heat extraction ratio. In this paper, the length of the horizontal wellbore and the size of each SRV remain constant. When the number of the fracturing stage is 7, the geothermal reservoir of interest is fully stimulated. It can be found that when the fracture number is between 7 and 15, there is only a slight difference in production temperature. However, the total heat extraction ratio increases with the fracture number. Fig.12(b) demonstrates the effect of the number of fractures on the service time and cumulative output thermal power of the EGS. It can be seen that there are three representative regions with the increase of the fracture number. The first region is regarded as inadequate stimulation (I: Nf <7), where the service time of EGS drops sharply, but the cumulative output thermal power rises

464

rapidly. The second region is called as full stimulation (II: 7≤Nf≤15), in which the

465 466 467 468 469 470 471

service time has a slight increase, and the cumulative output thermal power continues to increase at a relatively slow rate. The third region is considered to be over-stimulated (III: Nf >15), in which the service time begins to decrease again, but the cumulative output thermal power still increases slightly. In addition, excessive fracturing stages not only intensify the inter-fracture interference, but also substantially increase the capital investment. It can be concluded that in this case, the optimum number of fracturing stages is 15.

power.

17

473 474

6.4x1016

26

I

0.6 0.4

Service time (a)

0.8

0.2

20 25 Time (a)

30

35

II

III

4.8x1016

22

4.0x1016

20

3.2x1016 2.4x1016

18 16

1.6x1016

(b)

8.0x1015

14

0.0 40

5.6x1016

24

2

4

6

cumulative output thermal power (J)

1.0

Heat extraction ratio

Temperature (K)

472

470 460 450 440 Number of fractures 430 3 7 9 11 420 13 15 410 21 400 390 380 (a) 370 360 0 5 10 15

8 10 12 14 16 18 20 22 Number of fracturing stages

Fig. 12. Effect of the number of fracturing stages on heat extraction. (a) production temperature and total heat extraction ratio, (b) service time and cumulative output thermal power after lifetime. 480

T0=468 K

470 460 Temperature (K)

450 440 430 420 410 400

3 fractures 5 fractures 7 fractures 15 fractures

390 380 370

475 476 477 478 479 480 481 482 483 484 485 486 487 488 489

0

80

160

240 320 400 480 560 640 Distance from the production well (m)

720

800

Fig. 13. Effect of the number of fracturing stages on production temperature distribution.

Fig.13 shows the fluid temperature distribution along the production well at different number of fractures. When the fracturing stage is insufficient, the heat extraction around main fractures has a significant fingering phenomenon. When there are only three transverse fractures, local thermal breakthrough would occur earlier, resulting in a lower total heat extraction ratio. This is because the SRV induced by each main fracture is limited (130 m×800 m in plane). The permeability of unstimulated zones between main fractures is extremely low, where the cooling of the heat rock mass is dominated by the heat conduction regime. Therefore, the heat extraction efficiency is low in these unstimulated zones. As the fracture number increases, the overall reservoir of interest can get fully stimulated (more than 7 fractures). In permeability enhancement zones, the thermal transport by convection and conduction is the dominant process. It can be found that for 15 transverse fractures, the overall reservoir can be cooled down more uniformly.

490

5.3.4 Effect of well spacing on heat extraction

491 492 493 494 495 496 497

Fig.14(a) illustrates the effect of the well spacing on the production temperature and total heat extraction ratio. In this paper, it is assumed that the half-length of the transverse fracture is 350m, and the outermost periphery of the stimulation zone is 400m away from the injection wellbore. When the production well is placed outside the zone between the fracture tip and the outermost periphery of the stimulation zone (i.e., the well spacing is less than 350m or greater than 400m), the service time of the EGS can be prolonged. When the well spacing is less than 350m, the total heat 18

extraction ratio of the EGS is extremely low, only 0.11 after operation time of 40a. When the well spacing is greater than 400m, the total heat extraction ratio increases obviously, up to 0.61 after operation time of 40a. Fig.14(b) shows the changes in injection pressure and cumulative output thermal power with the distance (l) between the production well and the fracture front. When the production well intersects transverse fractures (l<0m), the injection pressure is relatively low. As the distance l increases, the injection pressure first increases moderately (e.g., ∆p=10.6MPa from -10m to 0m) and then increases significantly (e.g., ∆p=62.4MPa from 50m to 60m). More importantly, it can be seen from the cumulative output thermal power that the

507

heat extraction from the EGS is economically acceptable only when 0m≤l≤50m. After

508 509

stimulation treatment, the optimal zone of the production well can be determined by the fracture length and the microseismic cloud map. 4.0x1016

1.0 470

Temperature (K)

Well spacing

450

∆d=330m ∆d=350m ∆d=360m ∆d=410m

440 430 420

0.6

∆d=340m ∆d=360m ∆d=400m ∆d=420m

0.4 0.2

(a)

Heat extraction ratio

0.8

460

Injection pressure (MPa)

120 3.2x1016 100

410 0

510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525

5

10

15

20 25 Time (a)

30

35

Optimal zone for production well

80

1.6x1016

60 40

2.4x1016

(b)

8.0x1015

0.0 20 -30 -20 -10 0 10 20 30 40 50 60 70 80 Distance from from the fracture front (m)

0.0 40

Cumulative output thermal power (J)

498 499 500 501 502 503 504 505 506

Fig. 14. Effect of well spacing on heat extraction. (a) production temperature and total heat extraction ratio, (b) injection pressure and cumulative output thermal power after 40a.

Fig.15 shows the fluid temperature distribution along the production well at different well spacing. It can be clearly seen that when the production wellbore intersects transverse fractures (∆d=340m), the temperature drop is localized around these main fractures. This reveals the reason why the total heat extraction ratio is low for this case in Fig.14(a). When the position of the production well exceeds the fracture front, the cooling zone becomes larger, approximately equal to the width of each SRV (130m in this paper). Compared with the well spacing of 360m, the production temperature drop for the well spacing of 410m is slower before thermal breakthrough, after which the production temperature drops faster. This is because when the production well is placed outside the stimulation zone, higher injection pressures will be generated due to the low permeability of the unstimulated zone. At high injection pressures, the production temperature will drop rapidly after thermal breakthrough.

19

480

Temperature (K)

460 440 420 400 380 ∆d=340m (10a) ∆d=360m (10a) ∆d=410m (10a)

360 340

0

80

160

240

∆d=340m (40a) ∆d=360m (40a) ∆d=410m (40a)

320

400

480

560

640

720

800

526 527

Fig. 15. Effect of well spacing on production temperature distribution.

528

6. Conclusions

529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558

A conceptual EGS model combining horizontal wells with multistage fractures is proposed. Based on the fully coupled HT model developed in this paper, a series of numerical studies are conducted to investigate the effect of potential design parameters on the heat extraction process. The heat extract performance of the EGS is optimized by considering key evaluation parameters such as cumulative output thermal power, service time and injection pressure. Some remarkable conclusions are summarized as follows. (1) The cooling of the heat rock mass advances in an approximately spindle shape. The local heat extraction profile around each fracture is somewhat non-uniform even at the same injection rate. Due to the flow interference among inner fractures, the closer to both ends of the wellbore, the higher the local heat extraction ratio of the SRV. The inter-fracture interference may be mitigated by increasing the fracture spacing. (2) Increasing the injection rate will result in an earlier thermal breakthrough, and is ineffective for improving the heat extraction ratio and cumulative output thermal power at the end of the EGS lifetime. Compared to non-uniformly distributed injection rates, the equal injection rate for each segment contributes to achieving a good thermal sweep efficiency and preventing the local thermal short circuiting. (3) Compared to the effect of the fracture aperture, the increase in the permeability of the stimulation zone is more conducive to reducing the injection pressure and increasing the cumulative output thermal power. In multistage fracturing designs, more attention should be paid to shear stimulation to enhance the permeability of the surrounding reservoir. There may be an optimal reservoir permeability to match a certain fracture aperture (e.g., km=2×1015m2 for 0.7m of fracture aperture in this study). (4) Increasing the number of fracturing stages contributes to improving the cumulative output thermal power, but overall results in a decrease in the service time of EGS. Insufficient or excessive stimulation is detrimental to the economic development of EGS. For example, for a wellbore length of 820m, the optimum

Distance from the production well (m)

20

559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600

number of fracturing stages is 15, and the corresponding fracture spacing is 55.7m. (5) When the production well is located behind the fracture tip, heat extraction is localized around transverse fractures, resulting in an extremely low total heat extraction ratio and cumulative output thermal power. When the production well exceeds the outermost periphery of the stimulation zone, although the total heat extraction ratio is improved, the injection pressure is extremely high and the cumulative output thermal power is also undesirable. The preferable zone for the production well should be between the fracture tip and the periphery of the stimulation zone, which could be determined through the microseismic cloud map created by fracturing. Currently, it is assumed that the induced SRV is an equal-sized permeability enhancement zone. On-site hydraulic fracturing typically generates irregular SRVs with different sizes, which should be delineated by the scope of microseismic events. Moreover, the stimulation zone is simplified as an equivalent continuous porous medium that is inconsistent with the actual SRV containing a large number of activated natural fractures. We have developed a fully coupled thermal-hydraulic-mechanical (THM) model to investigate the heat extraction of naturally fractured geothermal reservoirs [37]. In future study, the HTM model combined with random fractures will be incorporated into the conceptual EGS model proposed in this paper. Acknowledgements This work was supported by the National Natural Science Foundation of China (Grant No.U1762216, Grant No.51574270, and Grant No.51504280), the Program for Changjiang Scholars and Innovative Research Team in University (IRT_14R58) and the Key Research and Development Program of Shandong Province (2019GGX103025). References [1] Panel M.L. The future of geothermal energy: Impact of enhanced geothermal systems (EGS) on the United States in the 21st Century. Geothermics 2006, 17:881-882. [2] Breede K., Dzebisashvili K., Liu X.L., Falcone G. A systematic review of enhanced (or engineered) geothermal systems: past, present and future. Geotherm. Energy 2013, 1(1):1-27. [3] Kumari W G P, Ranjith P G, Perera M S A, Li X, Li L H, Chen B K, Avanthi Isak B L, De Silva V R S. Hydraulic fracturing under high temperature and pressure conditions with micro CT applications: geothermal energy from hot dry rocks. Fuel 2018, 230:138-154. [4] Vidal J, Genter A, Schmittbuhl J. Pre- and post-stimulation characterization of geothermal well GRT-1, Rittershoffen, France: insights from acoustic image logs of hard fractured rock. Geophys. J. Int. 2016, 206(2): 845-860. [5] Lengliné O, BoubacarM, Schmittbuhl J. Seismicity related to the hydraulic stimulation of GRT1, Rittershoffen, 21

601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642

France. Geophys. J. Int. 2017, 208(3):1704-1715. [6] Bradford J, McLennan J, Tiwari S, Moore J, Podgorney R, Plummer M, Majer E. Application of hydraulic and thermal stimulation techniques at Raft River, Idaho: a doe enhanced geothermal system demonstration project. In: 50th U.S. Rock Mechanics/Geomechanics Symposium. American Rock Mechanics Association 2016. [7] Han S, Cheng Y, Gao Q, Yan C, Han Z. Experimental study of the effect of liquid nitrogen pretreatment on shale fracability. J Nat Gas Sci Eng 2018; 60:11-23. [8] Lu, S M. A global review of enhanced geothermal system (EGS). Renew Sust Energ Rev 2018, 81: 2902-2921. [9] Cladouhos T T, Petty S, Swyer M W, Uddenberg M E, Grasso K, Nordin Y. Results from Newberry Volcano EGS Demonstration, 2010–2014. Geothermics 2016 63: 44-61. [10] Cao W.J, Huang W.B, JiangF.M. A novel thermal–hydraulic–mechanical model for the enhanced geothermal system heat extraction. Int. J. Heat Mass Transfer 2016, 100:661-671. [11] Huang W.B, Cao W.J, Jiang F.M. Heat extraction performance of EGS with heterogeneous reservoir: A numerical evaluation. Int. J. Heat Mass Transfer 2017, 108:645-657. [12] Chen J, Jiang F. Designing multi-well layout for enhanced geothermal system to better exploit hot dry rock geothermal energy. Renew. Energ 2015, 74:37-48. [13] Dezayes C, Genter A, Valley B. Structure of the low permeable naturally fractured geothermal reservoir at Soultz. C.R. Geosci. 2010, 342(7-8):517-530 [14] Brown D.W., Duchane D.V., Heiken G., Hriscu, V.T. Mining the earth’s heat: hot dry rock geothermal energy. Springer Berlin Heidelberg. 2012. [15] Bendall B, Hogarth R, Holl H, McMahon A., Larking A., Reid P. Australian experiences in EGS permeability enhancement - a review of 3 case studies. In: Thirty-Ninth Workshop on Geothermal Reservoir Engineering, Stanford University 2014. [16] Serpen U, Aksoy N. Successful perforation operation experience in a geothermal well of Salavatli geothermal field. In: Fourtieth Workshop on Geothermal Reservoir Engineering, Stanford University 2015. [17] Zimmermann G, Reinicke A, Blöcher G, Moeck I, Kwiatek G, Brandt W, Regenspurg S, Schulte T, Saadat A, Huenges E. Multiple fracture stimulation treatments to develop an enhanced geothermal system (EGS) conceptual design and experimental results. In: the Word Geothermal Congress, Bali, Indonesia 2010. [18] Nair R, Peters E, Šliaupa S, Valickas R, Petrauskas S. A case study of radial jetting technology for enhancing geothermal energy systems at Klaipėda geothermal demonstration plant. In: 42nd workshop on Geothermal Reservoir Engineering, Stanford University 2017. [19] Petty S, Nordin Y, Glassley W, Cladouhos T T, Swyer M. Improving geothermal project economics with multi-zone stimulation: results from the Newberry Volcano EGS demonstration. In: Thirty-Eighth Workshop on Geothermal Reservoir Engineering, Stanford University 2013. [20] Lei Z H, Zhang Y J, Zhang S Q, Fu L, Hu Z J, Yu Z W, Li L Z, Zhou J. Electricity generation from a three-horizontal-well enhanced geothermal system in the Qiabuqia geothermal field, China: Slickwater fracturing treatments for different reservoir scenarios. Renew. Energ 2020, 145:65-83. [21] Li T, Shiozawa S, Mcclure M W. Thermal breakthrough calculations to optimize design of a multiple-stage Enhanced Geothermal System. Geothermics, 2016, 64:455-465. [22] Shiozawa S and McClure M. EGS designs with horizontal wells, multiple stages, and proppant. In: Thirty-Ninth Workshop on Geothermal Reservoir Engineering, Stanford University 2014. [23] Hofmann H, Weides S, Babadagli T, Zimmermann Günter, Moeck I, Majorowicz J, Unsworth M. Potential for enhanced geothermal systems in Alberta, Canada. Energy 2014, 69:578-591. 22

643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684

[24] Hofmann H, Babadagli T, Zimmermann Günter. Hot water generation for oil sands processing from enhanced geothermal systems: Process simulation for different hydraulic fracturing scenarios. Appl. Energ 2014, 113:524-547. [25] Zeng Y C, Su Z, Wu N Y. Numerical simulation of heat production potential from hot dry rock by water circulating through two horizontal wells at Desert Peak geothermal field. Energy 2013, 56:92-107. [26] Xia Y, Plummer M, Mattson E, Podgorney R, Ghassemi A. Design, modeling, and evaluation of a doublet heat extraction model in enhanced geothermal systems. Renew. Energ 2017, 105:232-247. [27] Jiang F, Chen J, Huang W, Luo L. A three-dimensional transient model for EGS subsurface thermo-hydraulic process. Energy 2014, 72:300-310. [28] Qu Z Q, Zhang W, Guo T K. Influence of different fracture morphology on heat mining performance of enhanced geothermal systems based on COMSOL. Int. J. Hydrogen Energ 2017, 42:18263-18278. [29] Zeng Y, Tang L, Wu N, Cao Y. Numerical simulation of electricity generation potential from fractured granite reservoir using the MINC method at the Yangbajing geothermal field. Geothermics 2018, 75:122-136. [30] Jie Zhang , Jingxuan Xie, Xueling Liu. Numerical evaluation of heat extraction for EGS with tree-shaped wells. Int. J. Heat Mass Transfer 2019. 134:296-310. [31] Song X, Shi Y, Li G, Yang R, Wang G, Zheng R, Li J, Lyu Z. Numerical simulation of heat extraction performance in enhanced geothermal system with multilateral wells. Appl. Energ 2018, 218:325-337. [32] Zhang F Z, Jiang P X, Xu R N. System thermodynamic performance comparison of CO2-EGS and water-EGS systems. Appl. Therm. Eng 2013, 61(2):236-244. [33] Bai B, He Y Y, Li X C, Li J, Huang X X, Zhu J L. Experimental and analytical study of the overall heat transfer coefficient of water flowing through a single fracture in a granite core, Appl. Therm. Eng. 2017, 116:79-90. [34] Baytas A.C. Thermal non-equilibrium natural convection in a square enclosure filled with a heat-generating solid phase, non-Darcy porous medium. Int. J. Energy Res. 2003, 27:975-988. [35] Chen Y, Ma G W, Wang H D. Heat extraction mechanism in a geothermal reservoir with rough-walled fracture networks, Int. J. Heat Mass Transfer 2018, 126:1083-1093. [36] COMSOL. Reference Manual Version 5.3a. COMSOL Corporation, Stockholm, Sweden, 2017. [37] Han S, Cheng Y, Gao Q, Yan C, Han Z, Zhang J. Investigation on heat extraction characteristics in randomly fractured geothermal reservoirs considering thermos-poroelastic effects. Energy Sci Eng 2019, 00:1-22. [38] Sun Z X, Zhang X, Yao J, Wang H X, Lv S H, Sun Z L, Huang Y, Cai M Y, Huang X Y. Numerical simulation of the heat extraction in EGS with thermal-hydraulic-mechanical coupling method based on discrete fractures mode. Energy 2017, 120:20-33. [39] Ghassemi A, Zhang Q. A transient fictitious stress boundary element method for porothermoelastic media. Eng Anal Bound Elem. 2004, 28(11):1363‐1373 [40] Huang N, Liu R, Jiang Y, Cheng Y, Li B. Shear-flow coupling characteristics of a three-dimensional discrete fracture network-fault model considering stress-induced aperture variations. J Hydrol 2019, 571: 416-424. [41] Zhao Z H. On the heat transfer coefficient between rock fracture walls and flowing fluid. Comput. Geotech. 2014, 59:105-111. [42] Perkins T K, and Kern L R. Widths of Hydraulic Fractures. J. Petrol. Technol. 1961: 937-949. [43] Nordgren R P. Propagation of a Vertical Hydraulic Fracture. SPE J 1972, 253: 306-314. [44] Lecampion B, Desroches J. Robustness to formation geological heterogeneities of the limited entry technique for multi-stage fracturing of horizontal wells. Rock Mech. Rock Eng. 2015, 48(6):2637-2644. 23

685 686 687 688 689 690 691 692 693

[45] Bennett T S, and Barker K., Operational aspects of placing proppant in a naturally fracture geothermal reservoir, GRC Transactions 1989, 13: 359-365. [46] Evans K F. Permeability creation and damage due to massive fluid injections into granite at 3.5 km at Soultz: 2. Critical stress and fracture strength. J. Geophys. Res. 2005, 110 (B4):B04204. [47] Dezayes C, Genter A, Valley B. Structure of the low permeable naturally fractured geothermal reservoir at Soultz. C R Geosci 2010, 342(7-8):517–530. [48] Stoddard T, McLennan J D, Moore J. Residual conductivity of a bauxite-propped geothermal system influence of self-propping, time, and closure stress. In: 46th US Rock Mechanics/Geomechanics Symposium. American Rock Mechanics Association 2012.

24

Highlights: A novel EGS model for horizontal wells with multistage fractures is designed. An EGS hydrothermal simulation method considering local thermal non-equilibrium is developed. Sensitivity analysis of potential factors to heat extraction performance is performed. Optimization results for key design parameters are presented.