Energy 120 (2017) 20e33

Contents lists available at ScienceDirect

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

Numerical simulation of the heat extraction in EGS with thermal-

hydraulic-mechanical coupling method based on discrete fractures
model
Zhi-xue Sun a, Xu Zhang a, Yi Xu b, Jun Yao a, *, Hao-xuan Wang a, Shuhuan Lv a,
Zhi-lei Sun c, Yong Huang a, Ming-yu Cai a, Xiaoxue Huang d
a
School of Petroleum Engineering, China University of Petroleum (East China), Qingdao 266580, China
b
State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China
c
Qingdao Institute of Marine Geology, Qingdao 266071, China
d
Key Laboratory of Efficient Utilization of Low and Medium Grade Energy, Ministry of Education, Tianjin University, Tianjin 300072, China

a r t i c l e i n f o a b s t r a c t

Article history: The Enhanced Geothermal System (EGS) creates an artificial geothermal reservoir by hydraulic fracturing
Received 9 June 2016 which allows heat transmission through the fractures by the circulating fluids as they extract heat from
Received in revised form Hot Dry Rock (HDR). The technique involves complex thermalehydraulicemechanical (THM) coupling
16 August 2016
process. A numerical approach is presented in this paper to simulate and analyze the heat extraction
Accepted 17 October 2016
process in EGS. The reservoir is regarded as fractured porous media consisting of rock matrix blocks and
Available online 2 January 2017
discrete fracture networks. Based on thermal non-equilibrium theory, the mathematical model of THM
coupling process in fractured rock mass is used. The proposed model is validated by comparing it with
Keywords:
Hot dry rock
several analytical solutions. An EGS case from Cooper Basin, Australia is simulated with 2D stochastically
Enhanced geothermal system generated fracture model to study the characteristics of fluid flow, heat transfer and mechanical response
Fractured rock mass in geothermal reservoir. The main parameters controlling the outlet temperature of EGS are also studied
THM coupling by sensitivity analysis. The results shows the significance of taking into account the THM coupling effects
Numerical simulation when investigating the efficiency and performance of EGS.
© 2016 Elsevier Ltd. All rights reserved.

1. Introduction by traditional techniques. The enhanced geothermal system (EGS)

creates an artificial heat reservoir by hydraulic fracturing process
1.1. Background which has been widely used to develop unconventional oil and gas
[12], this increases the permeability of the reservoir and allows heat
Majority of geothermal energy is stored in the hot (normally at transmission fluids (or working fluids, typically water) to be
temperature up to 150e650
C) dry or with little fluid crystalline circulated in the reservoir to extract heat from HDR economically
basement rocks located within subsurface of 3e10 km depths [1]. [13,14]. EGS represents the key technology for effective HDR heat
Hot Dry rock (HDR) geothermal heat has aroused interests around utilization, which has attracted broad attention in many countries
the world due to its cleanness, reproducibility, large abundance of around the world.
storage and its high temperature. It is considered as one of the most Extraction of heat from an EGS reservoir involves several
promising new energy resources of 21st century [2,3]. A vast coupled physical to geo-mechanical processes in the fractured rock
amount of research [4e11] was conducted all over the world since mass, including heat transfer, water seepage and solid deformation
the successful testing of Fenton Hill geothermal site in 1984. Rocks which are generally termed as “thermalehydraulicemechanical
in the HDR reservoir are dense and nearly impervious as they are (THM) coupling”. The spatial evolution of these multi-physical
deeply buried, making it difficult to extract stored thermal energy processes is complex and affected by a large number of parame-
ters. Therefore, the numerical method has been recognized as a
viable approach to model and analyze the coupling process in EGS
* Corresponding author.
reservoir [15]. Fluid migration, heat exchange and stress evolution
E-mail address: [email protected] (J. Yao). in the fractured rock mass can be effectively modeled using

http://dx.doi.org/10.1016/j.energy.2016.10.046
0360-5442/© 2016 Elsevier Ltd. All rights reserved.
 Z.-x. Sun et al. / Energy 120 (2017) 20e33 21

numerical approach. This plays a crucial role in many key technical model is implemented in the commercial finite element software,
issues for EGS development such as controlling the performance of COMSOL Multiphysics. The proposed model and the numerical
heat extraction, ensuring the stability and security of the wellbore, approach are validated by comparing with some analytical solu-
enhancing the geothermal reservoir recovery and keeping its tions. Finally, the model is used to simulate an EGS case with 2D
service-life for sustainable utilization period. random generated fracture network to study the characteristics of
In recent years, many useful models have been developed for flow, heat transfer and mechanical behaviors in HDR reservoir.
modeling the performance of EGS. Jiang et al. [16,17] presented a
three-dimensional (3D) transient model for EGS subsurface 2. THM coupling method and mathematical model
thermo-hydraulic process by which the geothermal reservoir is
treated as an equivalent porous medium of a single porosity. Xu 2.1. Basic assumptions
et al. [18] proposed a simplified approach to simulate the coupled
hydro-thermal system for EGS, capable of providing a detailed In this work, the following assumptions are made in order to
prediction of fluid flow and heat transfer in geothermal reservoir formulate the mathematical model:
based on an equivalent pipe network model. Shaik et al. [19]
developed a numerical procedure to simulate the heat extraction (i) The rock mass in EGS reservoir is treated as a 2D fractured
from naturally fractured geothermal systems by coupling fluid flow porous media consisting of rock matrix blocks and discrete
with heat transfer between the rock matrix and circulating fluid. It fractures; the rock matrix blocks can be simplified as
provides a dynamic treatment of the characteristic properties continuous porous media with isotropic properties, whose
(aperture, length and orientation) of individual fractures. However, permeability is much lower than that of fractures. The frac-
the effect of fractured rock mass deformation on the coupled tures created by hydraulic stimulation are the main flow
hydro-thermal process in EGS is overlooked in these studies. For a pathways in the reservoir.
better understanding of the performance of EGS, Bahrami et al. [20] (ii) The fractured porous geothermal reservoir is saturated with
studied several self-propped single fracture THM models based on single-phase liquid, i.e. water and the flow of fluid both in the
the poro-elastic theory. Zeng et al. [13] investigated HDR in the rock matrix block and the fractures obeys Darcy's Law.
Desert Peak geothermal field in American through a single vertical (iii) The reservoir pressure is in the range of 68.6 MPae80.4 MPa
and horizontal fractures. Zhao et al. [21] established a 3D THM (the corresponding boiling temperature is 787 K and 811 K,
coupling model of fractured media to simulate the extraction of respectively, referring to the literature [22]). Therefore, water
HDR geothermal energy, by which the fracture area was discretely in the reservoir (maximum temperature 473.15 K) does not
modeled by the Goodman joint element. However, these models vaporize, and is considered to be always in fluid in this work.
look too simplified to be used for modeling heat extraction in rock (iv) The heat exchange between water and rock matrix is ach-
mass containing dense fracture networks. ieved by advection and conduction process.
Although substantial efforts have been made on the numerical (v) The rock matrix in this study is considered to be thermal
modeling of THM coupling process in fractured rock mass, there is elastic based on small strain assumption.
still a challenge to accurately estimate the heat extraction process
in EGS and further theoretical research and numerical simulation
are needed [21]. Fractures and fracture networks are the funda-
2.2. The governing equation
mental components of EGS to determine its technical and economic
viability, whereas relatively fewer investigations have focused on
With the above-mentioned assumptions, the governing equa-
the simulation of THM coupling process in dense fracture networks.
tion of THM coupling model in fractured porous media can be
The deformation and stress variation of rock mass can result in
expressed as follows:
fracture opening and closing which may dramatically influence the
permeability during extracting geothermal energy in EGS. The
coupled impacts of THM process in fractures, especially how me- 2.2.1. Mass conservation equation
chanical process affects the fluid flow, heat exchange between rock The flow of water in deformable and saturated porous media is
and circulating fluid is essential for understanding the heat described by the following mass balance equation [23].
extraction process in EGS, which should be considered in the vp ve
simulation. On the other hand, more numerical models were S þ V$u ¼  (1)
vt vt
developed based on local thermal equilibrium theory other than
local thermal non-equilibrium theory, neglecting the heat ex- where, S is the constrained specific storage of the media, p is
change between rock matrix and fluid and failing to reflect the pressure, t is time, u is the volumetric flow flux, and e is the volu-
actual scenario of heat exchange occurring in the geothermal metric strain. When the flow is assumed to follow the Darcy's Law,
reservoir [17]. Therefore, a realistic THM coupling model that can the volumetric flow flux u can be written as:
adequately describe the fracture-stimulated geothermal reservoir
is critical for the numerical simulation of the enhanced geothermal k

u¼ Vp þ rf gVz (2)
system. h

1.2. Research objectives where, k is intrinsic permeability of porous media; h is fluid dy-
namic viscosity, rf is fluid density, g is the gravitational acceleration,
Considering the local thermal non-equilibrium, a mathematical and z is a unit vector in the direction over which the gravity acts.
model incorporating THM coupling effects is presented in this work The mass balance equation in discrete fractures can be described
for simulating the fractured EGS reservoir. The geothermal reser- as [23].
voir is regarded as a fractured porous media consisting of rock
vp vef
matrix blocks and discrete fractures. Deformation, water seepage df Sf þ Vt $uf ¼ df þ Qf (3)
and heat transfer processes in both of the rock matrix and fractures vt vt
are taken into account, as well as their interactions. The coupling
22 Z.-x. Sun et al. / Energy 120 (2017) 20e33

fluid.
kf

uf ¼ df Vt p þ rf gVt z (4) Assume that the heat exchange between rock matrix and frac-
h ture water follows the Newton's heat transfer law, the heat flow
from the rock to the fracture water per unit area is defined as [18].
where, Sf is the specific storage for the fracture, kf is the fracture's
permeability, df is the thickness of the fracture, ef is the volumetric


qf ¼ h Ts  Tf (11)
strain of the fracture; Vt denotes the gradient operator restricted to
the fracture's tangential plane, and Qf is the flow exchange on the
k vp
fracture surface between rock matrix and the fracture, Qf ¼  hf vn , where, “h” is the convection efficiency. When h is sufficient large,
n represents the normal direction on the fracture surface. the temperature of rock and that of water are equal at the fracture
surface.

2.2.2. Energy balance equation 2.2.3. Equilibrium equation

In this work, the thermal non-equilibrium theory is adopted to The equilibrium equation in the terms of stress tensor of linear
model the heat exchange between the rock matrix and the fluid elastic porous media with fluid saturated can be written as
flowing in the fractures, which reflects the actual scenario of heat
transfer occurring in EGS heat reservoir. It uses two energy equa- sij;j  aB p;i þ Fi ¼ 0 (12)
tions to describe the thermal movement in rock matrix and in the
fluid. The energy balance equation for heat transport in porous rock where, sij,j represents the divergence of the transpose of the Cauchy
matrix is given as [18,24]. stress tensor, “Fi” is the body force per unit volume in the i-coor-
dinate (i ¼ x,y in 2D), and aBp,i denotes the seepage body force
vTs
ð1  εÞrs Cs ¼ ð1  εÞls V2 Ts (5) resulting from the pore pressure. The effective stress law of porous
vt media,s0 ij ¼ sij  aB p is used, where aB is the Biot's constants.
Combined with the linear stress-strain constitutive law, the
vTf deformation equation of the rock matrix [21].
εrf Cf þ rf Cf $uVTf ¼ εlf V2 Tf (6)
vt
mui;jj þ ðl þ mÞuj;ji  aB p;i  K 0 aT Ts;i þ Fi ¼ 0 (13)
where, T is temperature, r is the density, C is the specific heat ca-
pacity, l is the heat conductivity. The subscripts ‘‘s” and ‘‘f” mean where ui is the displacement component; l and m are Lame's con-
solid and fluid respectively. “ε” is the porosity of the rock matrix. stants, defined as l ¼ En/[(1 þ n)(1  2n)], m ¼ E/2(1 þ n), in which E
0
The second term of Eq. (6) denotes the convection effects of and n are elastic modulus and Poisson's ratio respectively; K aTTs,i
water flow on the temperature distribution. Compared with the represents the thermal stress term, in which aT is the coefficient of
flow in fracture, water transport velocity in porous rock matrix is volumetric expansion corresponding to the bulk medium, and
0
much smaller. Therefore, it can be assumed that the temperature of K ¼ E/(1  2n) is the bulk modulus of the porous medium. Note that
porous water is identical to that of rock. The heat advection in the the thermal stress induced by temperature change is relative to the
pores has thus been overlooked. The equivalent energy equation for reference temperature or initial temperature. In this regard, the
modeling heat transfer in porous media considers equality of solid term “temperature increment”,Ts,i is used here to distinguish from
and fluid temperatures as described by the equation below; current temperature state.
The discrete fracture is approximated by a pair of surfaces be-
vTs
ðrCÞeff ¼ leff V2 Ts (7) tween which normal and shear displacements are permissible. The
vt deformation equation of the rock mass fracture can be written as
where, (rC)eff ¼ (1  ε)rsCs þ εrfCf and leff ¼ (1  ε)ls þ εlf are [21].
respectively the effective heat capacity and the effective thermal
un ¼ s0 n =kn ; us ¼ s0 s =ks (14)
conductivity obtained by volume average.
The energy balance equation for heat transfer and heat advec-
tion in the discrete fracture is written as [25,26]. s0 n ¼ sn  aB p; s0 s ¼ ss (15)
0
vTf
 whereu,s,s and k denote the displacement, total stress, effective
df rf Cf þ df rf Cf ,uf Vt Tf ¼ Vt $ df lf Vt Tf (8) stress and stiffness, respectively. The subscripts ‘‘n” and ‘‘s” repre-
vt
sent normal and tangential directions to the fracture plane
At the fissure surface, temperature of rock, Ts and that of water, respectively.
Tf are not equal. Therefore, heat exchange takes place between rock
matrix and water flowing in the fracture. The heat loss of the rock 2.3. Coupling effects
matrix to water,qf can be incorporated into the energy balance
equations as: 2.3.1. HM coupling effect on the evolution of permeability
The mechanical process in EGS can cause severe deformation of
vTs
ðrCÞeff ¼ leff V2 Ts (9) the pore and fracture structure in the rock mass. Therefore, the
vt
permeability of the rock mass, one of the most important param-

 eter controlling the feasibility and the viability of EGS [27] will
vTf
df rf Cf þ df rf Cf $uf Vt Tf ¼ Vt $ df lf Vt Tf þ qf (10) change as it is stress or strain dependent. Particularly, in the rock
vt mass containing dense fracture networks, the more deformable
Notice that the fracture thickness,df is used to calculate the fractures will be the primary flow pathways for the work fluid.
volumetric thermal flux. “qf” is introduced to describe the heat Thus, stress variation can bring a considerable change to the
exchange in the reservoir. The negative sign on this term in Eq. (9) permeability. The understanding of the evolution law of the frac-
means the heat is extracted from the rock and the positive sign on ture's permeability is an important aspect to study the hydraulic-
this term in Eq. (10) means the heat is absorbed by the working mechanical (HM) coupling effect in fractured rock mass. The
 Z.-x. Sun et al. / Energy 120 (2017) 20e33 23

fracture's permeability is strongly related with its thickness or the essential features of thermal extraction through heat conduc-
aperture, which is dependent on effective normal stress. In this tion from the rock matrix to surrounding fracture. A single rect-
work, the permeability takes the form [28,29]: angular, horizontal fracture of constant width 2df separates two
blocks of homogeneous, isotropic, impermeable rock (see Fig. 1). A
kf ¼ k0 expð  as0 n Þ (16) Cartesian coordinate system has been placed such that x¼0 co-
incides with the inlet.
where kf is the fracture's permeability, k0 is a baseline permeability, In the single fracture model, the rock is assumed to extend in y
s0 n is the effective normal stress on the fracture plane, and a is a direction to ±∞ value, the rock matrix and the fracture extend to
normalizing constant, which depends on the fracture's develop- infinity in the x direction. Initially (t¼0), this semi-infinite system is
ment characteristic. at a uniform temperature T0. During the heat extraction phase,
Thus, the fracture's permeability is a function of the stress state. water is injected at a constant temperature Tin at constant flow
Since the rock matrix is assumed to be much more impervious than velocityuf. The water flows parallel through the fracture to the
the fractures and its deformation is relatively smaller, the perme- outlet. Furthermore, the following assumptions are made: The
ability of the rock matrix is assumed to be constant. temperature variation in y corresponding to water is negligible. As
the aperture of the fracture is very small, conduction in the vertical
2.3.2. TH coupling effect on the evolution of fluid properties direction,y in both the fracture and the rock formation is neglected.
The property of the work fluid is another important factor Thermo-physical properties of both water and rock are constant
influencing the coupling effects. Under high pressure and high and independent of temperature. An analytical solution of the
temperature conditions in the deep underground reservoir, the temperature in the fracture can be obtained from the plane-
density of water,rf is no longer a constant value, this can be symmetric solution by Barends [32].
expressed as a function of pressure and temperature [21],
. 0 1
1 rf ¼ 3:086  0:899017ð4014:15  TÞ0:147166 . !
B ls x rf Cf df C x
 0:39ð658:15  TÞ 1:6
ðp  225:5Þ þ d (17) Tf ¼ T0 þ ðTin  T0 ÞerfcB
@ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
C
AU t  u (19)
2 ls =rs Cs ðuf t  x f
where T is the water temperature, p is absolute water pressure,
and d is a function of “p” and “T”; in general, the value of d does not
exceed 6% of 1/rf. The temperature of water also affects its dynamic where Tf is the temperature in the fracture, r is the density, C is the
viscosity,h ¼ yrf in which y is the kinematic viscosity of water. A specific heat capacity, l is the heat conductivity, and the subscripts
commonly used empirical formula for the kinematic viscosity is ‘‘s” and ‘‘f” are the rock matrix and water respectively. In this
[30]. equation, erfc is the complementary error function, and U is a unit
step function.
0:01775 A numerical solution is presented to model this single fracture
y¼ (18) heat extraction problem. In the numerical solution, a finite domain
1 þ 0:033Tf þ 0:000221Tf2
in 100 m  100 m is modeled for the rock matrix. Some parameters
Therefore, the density and viscosity of water are directly related used in this model are listed in Table 1. Fig. 2 compares the nu-
with the temperature, which can dramatically change the coupling merical and analytical results for this model during a 100 day
process of water transport and heat transfer in the reservoir. extraction period. Fig. 2(a) shows the temperature evolution at
three different positions (x ¼ 10 m, 20 m and 50 m) along the
3. THM coupling numerical model and simulation solution fracture. Temperature distribution in the fracture at different times
is plotted in Fig. 2(b). It can be seen that the numerical solution is in
The mathematical model of the THM coupling analysis in frac- good accordance with the analytical one, except a minor deviation.
tured porous EGS reservoir is composed of the aforementioned This deviation may be attributed to the fact that the analytical so-
equations and the corresponding initial and boundary conditions. lution is given based on the assumption that the model is a semi-
The Finite Element Method (FEM) is employed to solve the THM infinite system, while the numerical solution is obtained from a
coupling model. However, the FEM solution is still a huge and model with limited geometry (i.e., 100 m  100 m).
complex work since that there is a large amount of dependent
variables in the coupling model.
In this work, a full coupling solution of the solid deformation,
water seepage and heat transfer is developed in COMSOL Multi-
physics, which is a well-known tool providing full access to solving
partial differential equations [31]. The stiffness matrix of fractures
is added by modifying the weak terms of the so called interior
boundaries which represent the fractures [25,26].

4. Validation of numerical model

4.1. Validation of the model for thermo-hydraulic analysis

In order to examine the accuracy of the proposed model and its

numerical implementation in solving the typical thermo-hydraulic
problem in the fractures, a 2D single fracture model is studied in
this section. Although realistic EGS reservoir consists of complex
fracture networks, a single fracture model can adequately capture Fig. 1. Illustration of a single fracture model.
24 Z.-x. Sun et al. / Energy 120 (2017) 20e33

Table 1 coupling model and the numerical implementation, though no

Model parameters. fractures are involved in this problem.
Parameter Unit Value

rf kg/m3 1000 5. An EGS case study

Cf J/kg/K 4200
rs kg/m3 2700
An EGS case study is carried out based on the published reser-
Cs J/kg/K 1000
ls W/m/K 3 voir information from the Habanero Engineered Geothermal Sys-
uf m/s 0.01 tem in central Australia, which has targeted a conductive-style
m Pa$s 0.001 geothermal resource. High-heat-producing granites of the Big Lake


T0 C 80 Suite are overlain by a thick sedimentary package (the stacked


Tin C 30
Cooper, Eromanga and Lake Eyre basins) that provides sufficient

Fig. 2. Comparison between the analytical and numerical solutions of the single fracture model.

4.2. Validation of the model for THM coupling analysis

In this section, the proposed model is validated by simulating

the 1D thermal consolidation phenomenon to further demonstrate
its capability of fully THM coupling analysis. The thermal consoli-
dation problem involves the mechanical deformation, temperature
variation and pore pressure dissipation of saturated porous media,
which is a typical THM coupling process. Fig. 3 illustrates the
problem of 1D thermal consolidation of a linear elastic soil column
under constant surface compression and constant surface tem-
perature. Bai [33] deduced the analytical solution of the 1D thermal
consolidation problem by finite Fourier transform and its inverse
transform. With the proposed model, a THM coupling analysis is
performed for seeking a numerical solution to this problem.
In the numerical model, the saturated soil column is 7 m in
height, at initial temperature of 10
C and pore pressure of
1  104Pa. The vertical surface compression load on the soil column
is 1  104Pa and the temperature and fluid boundary pressure at
the top surface are 60
C and 0Pa respectively. As shown in Fig. 2, all
the boundaries are assumed to be thermal insulated and imper-
Fig. 3. Illustration of the 1D thermal consolidation problem.
vious with displacements constrained normally, except at the top
surface. Material parameters used in the numerical model are listed
in Table 2. In order to make the numerical solution and the thermal insulation to trap the radiogenic heat, resulting in tem-
analytical one comparable, some parameters for both solid and peratures above 244
C at depths less than 5 km. Development of
fluid properties are set to be equal according to literature [33]. the project began in 2002, and Geodynamics Limited demonstrated
The results of THM coupling analysis are plotted in Figs. 4e6, proof-of-concept in September 2013 through sustained operation
which respectively show the displacement, pore pressure and of a 1 MW pilot plant at surface. The reservoir for the engineered
temperature evolution at different position in the soil column. The system is within the large granite plutons of the Big Lake Suite,
results of the analytical model by Bai [33] are also presented for which are intersected at depths of between 3600 and 3800 m in the
comparison. It can be observed that the numerical results agree Innamincka area. Four wells have been drilled at the Habanero site
well with the analytical model, this verifies the accuracy of the THM
 Z.-x. Sun et al. / Energy 120 (2017) 20e33 25

Table 2
Material parameters used in the thermal consolidation problem.

Parameter Unit Value

Elastic modulus Pa 6.0  107

Poisson ratio 0.4
Biot's coefficient 1.0
Heat conductivity W/m/K 0.836
Linear thermal expansion 1/K 3.0  107
Heat capacity J/kg/K 167.2
Hydraulic conductivity m/s 4.0  106
Porosity 0.2
Density of water kg/m3 1000
Density of soil kg/m3 2000
Soil specific storage 1/Pa 0

Fig. 6. Variation of temperature evolution with time at different positions in the soil
column.

as shown in Fig. 8. In the EGS, hydraulic fracturing is carried out by

injecting water at high pressure into the target region. It allows the
tectonically stressed rock to fail in friction or shear, thus opening,
extending and connecting fractures, enhancing their permeability.
An artificial geothermal reservoir is created and heat can be
extracted from the rock by circulating water. Since the geothermal
reservoir is buried deeply underground, the detailed geometry of
these fractures is difficult to be directly recorded, which presents a
challenge for EGS modeling. Therefore, this study uses stochasti-
cally generated fracture networks to represent the EGS reservoir,
and model the interacting thermal, hydraulic and mechanical
Fig. 4. Displacement evolution with time at different positions in the soil column. process in geothermal heat extraction.

5.1. Computational model

The realistic EGS reservoir spans over thousands of meters and

contains dense fracture networks in 3D space. A 3D simulation is
difficult to carried out at present due to the extremely large
calculation workloads [24]. For the sake of computational efficiency
and accuracy, a 2D stochastically fracture network is modeled to
represent the geothermal reservoir and is used to predict the per-
formance of EGS heat extraction.
Fig. 9 shows the 2D computational model for EGS reservoir (in
300 m  300 m). Two sets of fractures with different dip angles (i.e.,
30
and 110
) are generated. Trace lengths of the fractures follow a
normal distribution with mean value of 30 m and variance of 10 m.
The stochastically generated fracture network model shows a
strong anisotropy, which can represent the fractured rock mass in
the statistical sense and thus reflects the main features of THM
coupling process in EGS [25].
Fig. 5. Pore pressure evolution with time at different positions in the soil column.

5.2. Computational conditions and parameters

e Habanero-1(H-1), Habanero-2 (H-2), Habanero-3 (H-3), and
Habanero-4 (H-4) (see Fig. 7). Multiple hydraulic stimulations have 5.2.1. Initial and boundary conditions
been performed at Habanero to improve the existing fracture The entire process of heat production in EGS is simulated with
permeability in the granite and create a reservoir that enables the aforementioned model. The performed THM coupling analysis
appropriate heat-exchange with the surrounding granite and suf- is spread over a total period of 40a (years) from time water was
ficient fluid-flow between injection and production wells. injected, with each time step of 1d. Initial and boundary conditions
Currently, the total stimulated area is on the order of 4 km2, it has are as follows:
been mapped from micro seismic monitoring during and after the
successive hydraulic stimulations. (1) Seepage field: The given pressure boundary conditions are
Typically, the subsurface components of an EGS contain injec- used to ensure water circulating in the system. The left
tion wells, production wells, geothermal reservoir, and rock boundary is assumed to be the injection well with a given
enclosing the geothermal reservoir (i.e., the base and cap rock) [16] pressure of 80.4 MPa, and the right one is the production well
26 Z.-x. Sun et al. / Energy 120 (2017) 20e33

Fig. 7. Location of the Habanero EGS site: (a) approximate geographic location as indicated by the green box, (b) location of the Habanero wells. (For interpretation of the references
to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 8. Schematic of the enhanced geothermal system.

at 68.6 MPa. For the other external conditions an imperme-

able boundary is selected. The initial water pressure in the Fig. 9. The 2D fracture network model for EGS.
reservoir is assumed as 71.5 MPa before water is injected.
(2) Thermal field: The temperature at the injection well is 20
C.
The top and bottom boundaries are thermal insulated. The fractured rock mass, which may help in understanding their
initial temperature in the reservoir is 200
C for both rock coupling intersections.
and water.
(3) Displacement field: All the external boundaries of the
reservoir are constrained in normal directions. Since engi- 5.2.2. Model parameters
neers are usually more concerned with the disturbance to Computational parameters used in this study are listed in
the stress distribution caused by EGS development and heat Table 3.
production, the initial geo-stress in the reservoir is not taken Eq. (16) is used to describe the permeability evolution of frac-
into account. This study focuses on the effects of pressurized tures affected by HM interaction. For the convenience in analysis,
water injection and heat extraction on the stress field in an assumed value of a as 0.2  106Pa1 is used in this study.
When tensile stress appears on the fracture (i.e., the effective
 Z.-x. Sun et al. / Energy 120 (2017) 20e33 27

Table 3
Computational parameters in EGS modeling.

rf/(kg/m3) h/(Pa$s) g/(m/s2) Cf/(J/kg/K) lf/(W/m/K) h/(W/m2/K) aB

1000 0.001 10 4200 0 3000 1.0

rs/(kg/m3) k/(m2) S/(1/Pa) E/(GPa) v Cs/(J/kg/K) ls/(W/m/K)

2700 1.0  1018 1.0  108 30 0.25 1000 3

aT/(1/K) ε df/(m) kf/(m2) Sf/(1/Pa) kn/(GPa/m) ks/(GPa/m)

5.0  106 0.0001 0.0005 1.0  1011 1.0  109 1200 400

normal stress is positive), the fracture will open and its perme- conduction, the temperature of the rock matrix block also de-
ability increases, otherwise the compressive stress will be creases gradually forming the low-temperature zone. Meanwhile,
accounted and its permeability decreases. the temperature at the production well has no significant change,
indicating that the system can maintain a steady heat production.
5.3. Results and analysis As more heat is extracted from the reservoir, the low-temperature
zone expands. It can be observed that the temperature drop is more
5.3.1. Results of the hydraulic process rapid near some fracture concentrated regions. These connected
A better understanding of the hydraulic process is important for fractures form the main flow pathways for circulating fluid. The
controlling the performance of EGS, since the fluid flow velocity in effect of heat convection is considerably strong in these connected
the reservoir may directly affect the heat exchange process. Fig. 10 fractures due to the high flow velocity. The interactions between
shows the velocity vectors in the reservoir at t ¼ 40a and colors thermal and hydraulic processes are consistent with the results in
indicate water temperature. It can be seen that the injected low- Refs. [18,25].
temperature water flows to the production well along the frac- The numerical results exhibit strong heterogeneity and anisot-
tures under high pressure. ropy of the thermal distribution in EGS. The fractures created by
Fig. 11 shows the pressure distribution in the reservoir at hydraulic stimulation play a vital role in enhancing the conductivity
different times. At the initial running stage, water flows mostly of rock mass and increase the efficiency of heat extraction, which
through the fractures, and the water pressure in the fractures are the key feature of EGS. The proposed model is able to reflect the
rapidly increases to 80.4Mpa. The water pressure of the rock matrix mechanism that the heat in EGS reservoir transfers mainly in
gradually increases since its permeability is extremely lower, it gets fractures, which might help in understanding the thermal-
approximately in accordance with the water pressure of the frac- hydraulic coupling effect.
tures at t ¼ 10a. After about 20a of running the simulation, the
water pressure in the whole system changes only a little and thus 5.3.3. Results of the mechanical process
the hydraulic process nearly remains stable. Water circulates Fig. 13 shows the contours of displacement in the reservoir at
through the fractures with a fixed pressure gradient, keeping the different times. At the initial stage of heat extraction, the injection
EGS geothermal extraction running effectively. of highly pressurized water and the rapid temperature drop near
injection well cause tremendous changes to the displacement dis-
5.3.2. Results of the thermal process tribution of the rock mass. While at the later stage (after 20a),
The temperature distribution in the reservoir at different times nevertheless the water pressure in the system nearly remains sta-
is shown in Fig. 12. At the initial running stage, heat exchange takes ble, the displacement distribution also exhibits considerable
place between the rock matrix and the fluid flowing in the fracture change. This could attribute to that the temperature in the reservoir
due to the injection of low-temperature water then the tempera- varies all the way, which cause thermal stress in the rock mass. As
ture of fluid increase immediately. With the effect of heat the low-temperature zone extends, more rock matrix contract due
to cooling. The evolution of both hydraulic and thermal processes
affects the mechanical process.
On the other hand, the stress variation in the reservoir can in
turn influence the hydraulic process of circulating fluid and thus
the heat extraction. To illustrate this clue, some feature points in
fractures, #1, #2, #3 and #4 are selected as sampling points (see
Fig. 14). Figs. 15 and 16 show the permeability evolution and the
effective normal stress evolution with extraction time at point #1 to
#4 respectively. It can be found that the permeability of these four
points increases considerably as the extraction time increases and
the fractures near the production well exhibit larger permeability.
The evolution of effective normal stress at these points also shows a
similar trend. Under the effect of highly pressurized water injec-
tion, the deformation of the fractures and rock matrix blocks has
enhanced the system's hydraulic conductivity and further acceler-
ated the heat transfer rate. The rapid variation in temperature can
also induce the change of stress distribution in the reservoir. The
three physical fields presented in the EGS reservoir are strongly
related and interact with each other.

5.3.4. The outlet temperature of the EGS system

Fig. 10. Velocity vectors in the reservoir at t ¼ 40a. The outlet temperature in the production well is one of the
28 Z.-x. Sun et al. / Energy 120 (2017) 20e33

Fig. 11. Pressure distribution in the reservoir at different times.

important indexes to assess the production performance of EGS and The summation term represents the fractures while the inte-
its lifespan. Fig. 17 shows the simulated temperature distribution in gration term represents the rock matrix blocks.
the production well at different time. Before t ¼ 20a, temperature in Fig. 18 shows the evolution of the average outlet temperature
the production well is stable at 200
C. At this period, the system with extraction time obtained by different inlet temperature (i.e.,
can maintain a steady heat production. As the low-temperature 10
C, 20
C and 30
C). A stable running condition for about 20a can
zone expands in the reservoir, temperature in the production well be observed for all these three cases. During the period of 20a to
begins to drop after t ¼ 20a. The temperature drop in the produc- 40a, the low-temperature zone expands gradually to the produc-
tion well is not uniform, faster in areas with dense fractures since tion well and the heat production of the system decreases. By
that the heat flux in the fractures is larger. lowering the inlet temperature, the faster the average outlet tem-
Sensitivity analysis with regard to some key model parameters perature drops. At the time of t ¼ 40a, the average outlet temper-
are performed, in order to further evaluate their effects on the EGS ature is approximately 75%e85% of the reservoir's initial
outlet temperature and its running performance. The average temperature, which intensely reduces its efficiency.
outlet temperature is calculated by the following equation [25], The thermal conductivity of rocks is typically in the range of
Z 1e3 W/m/K [34]. Therefore, different values for rock thermal con-
P ductivity are employed to study its effect on the evolution of the
uf df Tf þ uTs dy
Tout ¼ Z (20) average outlet temperature. As shown in Fig. 19, the average outlet
P temperature drops slightly faster when the rock thermal conduc-
uf df þ udy
tivity is lower. This observation shows similar trend to Chen et al.’s
 Z.-x. Sun et al. / Energy 120 (2017) 20e33 29

Fig. 12. Temperature distribution in the reservoir at different times.

study [25]. However, the rock thermal conductivity seems to have a which can lead to the increase in fracture aperture and its
relatively limited effect on the average outlet temperature. conductive properties. The roles that mechanical process plays in
The evolution of the average outlet temperature against the overall heat production are strongly associated with the hy-
different fracture permeability is illustrated in Fig. 20. It can be draulic process in EGS.
found that the average outlet temperature is sensitive to the The evolution of the average outlet temperature, obtained by
permeation properties of the fractures. When the fracture perme- different well distance between injection well and production well
ability is relatively lower at 0.5  1011 m2, the outlet temperature (i.e., 300 m, 360 m and 420 m), is shown in Fig. 22. The production
is observed to be stable during 40a. However, the average tem- temperature differs apparently in the case of different well dis-
perature drops rapidly to 40
C during the 40-year period when the tance. For well distance larger than 420 m, the system can maintain
fracture permeability is larger. Fractures are the main pathways for a stable running condition for about 40 years. However, the average
heat conductivity, therefore are crucial in the thermal evolution outlet temperature drops sharply faster and the stable running
process. period gets shorter when the well distance become smaller.
The effect of different rock elastic modulus is also investigated. Considering that the well distance may influence the residence
As depicted in Fig. 21, the temperature drop varies considerably time as well as heat exchange amount between hot rock and
when the rock elastic modulus is changed. It is accredited that the injected fluid, appropriate well distance is crucial to the perfor-
rock deformation is much more significant with a higher modulus, mance and service-life of EGS.
30 Z.-x. Sun et al. / Energy 120 (2017) 20e33

Fig. 13. Contours of the displacement in the reservoir at different times (unit in m).

6. Conclusions

A 2D numerical model based on FEM is developed for modeling

the thermalehydraulicemechanical coupling process in fractured

Fig. 14. Four sampling points in the fractures. Fig. 15. Permeability evolution of point #1-#4 with extraction time.
 Z.-x. Sun et al. / Energy 120 (2017) 20e33 31

Fig. 16. Effective normal stress evolution of point #1-#4 with extraction time. Fig. 18. Evolution of the average outlet temperature with different inlet temperature.

Fig. 19. Evolution of the average outlet temperature with different rock thermal
conductivity.

Fig. 17. Temperature distribution in the product well at different times.

rock mass, which can be utilized to study the heat extraction by

enhanced geothermal system in HDR. The proposed model and the
numerical approach are validated by comparing with some
analytical solutions. A THM coupling case analysis is presented to
study the main features of heat transfer, water flow and mechanical
behaviors in EGS reservoir.

(1) In the proposed model, rock mass is treated as a fractured

porous media consisting of rock matrix blocks and discrete
fractures, able to effectively model the complex geometry of
EGS reservoir and describe the realistic interactions between
the rock matrix and the circulating fluid.
(2) The presence of discrete fractures created by hydraulic
stimulation is the key feature of EGS since the fractures form
the main flow pathways in the reservoir. The fracture net-
works play a vital role in enhancing the conductivity of rock
mass and increasing heat extraction efficiency necessary to
be explicitly modeled. Fig. 20. Evolution of the average outlet temperature with different fracture
permeability.
32 Z.-x. Sun et al. / Energy 120 (2017) 20e33

improving the procedure of EGS design. Nevertheless, the 2D re-

sults obtained by stochastically generated fracture networks may
be less representative. A 3D numerical model with actual fracture
data should be studied in future work to fully integrate and
describe the actual features in EGS, which might be helpful for the
development of HDR geothermal energy.

Acknowledgements

This study was jointly supported by the National Natural Science

Foundation of China (Grant No. 51404291,
No.41502131&No.51274226). Fundamental Research Funds for the
Central Universities (NO.14CX05024A), Shandong Natural Science
Foundation (NO. 2013ZRE28068& No.ZR2013DL011). We are
grateful to all staff involved in this project, and also wish to thank
the journal editors and the reviewers whose constructive com-
ments improved the quality of this paper greatly.

References
Fig. 21. Evolution of the average outlet temperature with different rock elastic
modulus.
[1] Tester JW, Anderson BJ, Batchelor AS, Blackwell DD, DiPippo R, Drake EM, et al.
The future of geothermal energy, impact of enhanced geothermal systems
(EGS) on the United States in the 21st Century. 2006. Assessment by an MIT-
led interdisciplinary panel (J.W. Tester, Chairman).
[2] Duchane D. Progress in making hot dry rock geothermal energy a viable
renewable energy resource for America in the 21st century. In: Energy con-
version engineering conference; 1996. p. 1628e32.
[3] Wei G, Meng J, Du X, Yang YP. Performance analysis on a hot dry rock
geothermal resource power generation system based on kalina cycle. Energy
Procedia 2015;75:937e45.
[4] Duchane D. Hot dry rock: a realistic energy option. Bull Geotherm Resour
Counc 1990;19(3):83e8.
[5] Duchane D. International programs in hot dry rock technology development.
Bull Geotherm Resour Counc 1991;20(5):135e42.
[6] Brown DW. Recent testing of the HDR reservoir at Fenton Hill New Mexico.
Bull Geotherm Resour Counc 1993;22(9):208e14.
[7] Chamorro CR, Mondejar ME, Ramos R, Segovia JJ, Martin MC, Villamanan MA.
World geothermal power status: energy, environmental and economic study
of high enthalpy technologies. Energy 2012;42(1):10e8.
[8] Chamorro CR, Garcia-Cuesta JL, Mondejar ME, Linares MM. An estimation of
the enhanced geothermal systems potential for the Iberian Peninsula. Renew
Energy 2014;66:1e14.
[9] Zhao YS, Wan ZJ, Kang JR. Introduction to HDR geothermal development.
Beijing: Science Press; 2004.
[10] Wan ZJ, Zhao YS, Kang JR. Forecast and evaluation of hot dry rock geothermal
resource in China. Renew Energy 2005;30:1831e46.
[11] Feng ZJ, Zhao YS, Zhou AC, Zhang N. Development program of hot dry rock
Fig. 22. Evolution of the average outlet temperature with different well distance. geothermal resource in the Yangbajing basin of China. Renew Energy
2012;39:490e5.
[12] Zhu GP, Yao J, Sun H, Zhang M, Xie MJ, Sun ZX, et al. The numerical simulation
of thermal recovery based on hydraulic fracture heating technology in shale
(3) In the EGS case study, the distributions of temperature, water gas reservoir. J Nat Gas Sci Eng 2016;28:305e16.
pressure and deformation are obtained by numerical simu- [13] Zeng YC, Su Z, Wu NY. Numerical simulation of heat production potential from
hot dry rock by water circulating through two horizontal wells at desert peak
lation. The coupling analysis F. Numerical results shows geothermal field. Energy 2013;56:92e107.
strong coupling effects of THM in the heat extraction of EGS. rez MC, Morales MP, D'Amico S, Liarte LA. Enhanced geothermal
[14] Olasolo P, Jua
For a better understanding of the performance of EGS, systems (EGS): a review. Renew Sustain Energy Rev 2016;56:133e44.
[15] McDermott CI, Randriamanjatosoa AR, Tenzer H, Kolditz O. Simulation of heat
further studies on the THM coupling in fractured rock mass extraction from crystalline rocks: the influence of coupled processes on dif-
are required. ferential reservoir cooling. Geothermics 2006;35(3):321e44.
(4) Numerical simulation can be used to evaluate the production [16] Jiang FM, Luo L, Chen JL. A novel three-dimensional transient model for
subsurface heat exchange in enhanced geothermal systems. Int Commun heat
performance of EGS and its lifespan. Sensitivity analysis are
mass Transf 2013;41:57e62.
performed to understand the effects of some model param- [17] Jiang FM, Chen JL, Huang W, Luo L. A three-dimensional transient model for
eters on the average outlet temperature of EGS. The results EGS subsurface thermo-hydraulic process. Energy 2014;72(7):300e10.
show that the average outlet temperature is sensitive to the [18] Xu C, Dowd PA, Zhao FT. A simplified coupled hydro-thermal model for
enhanced geothermal systems. Appl Energy 2015;140:135e45.
conductive properties of fractures, the inlet temperature and [19] Shaik AR, Rahman SS, Tran NH, Tran T. Numerical simulation of fluid-rock
well distance between injection well and production well, coupling heat transfer in naturally fractured geothermal system. Appl
whereas the rock thermal conductivity has a relatively Therm Eng 2011;31(10):1600e6.
[20] Bahrami D, Danko G, Fu P, Guo B, Podgorney R, White M, et al. Poroelastic and
limited effect. self-propped single fracture THM models for EGS studies. In: Proceedings of
fourtieth workshop on geothermal reservoir engineering. Stanford: Stanford
In this study, the numerical simulation and THM coupling University; 2015.
[21] Zhao YS, Feng Z, Yang D, Liang WTHM. (Thermo-hydro-mechanical) coupled
analysis of heat extraction in EGS are limited in two dimensional. mathematical model of fractured media and numerical simulation of a 3D
The proposed 2D model provides the means for more insight into enhanced geothermal system at 573 K and buried depth 6000e7000 M. En-
fluid-heat flow in deep geothermal reservoirs that may assist in ergy 2015;82:193e205.
[22] Thomson GW. The antoine equation for vapor-pressure data. Chem Rev
 Z.-x. Sun et al. / Energy 120 (2017) 20e33 33

1946;38:1e39. [29] Miller SA. Modeling enhanced geothermal systems and the essential nature of
[23] Liang B, Jiang HQ, Li JJ, Gong CC. A systematic study of fracture parameters large-scale changes in permeability at the onset of slip. Geofluids
effect on fracture network permeability based on discrete-fracture model 2015;15(1e2):338e49.
employing Finite Element Analyses. J Nat Gas Sci Eng 2016;28:711e22. [30] Wu C. Hydraulics. second ed. Beijing: Higher Education Press; 1983 (in
[24] Saeid S, Al-Khoury R, Barends F. An efficient computational model for deep Chinese).
low-enthalpy geothermal systems. Comput Geosci 2013;51:400e9. [31] COMSOL Multiphysics User's guide, version 5.0. Massachusetts, USA, Bur-
[25] Chen BG, Song EX, Cheng XH. A numerical method for discrete fracture lington: COMSOL Inc; 2014. www.comsol.com.
network model for flow and heat transfer in two-dimensional fractured rocks. [32] Barends F. Complete solution for transient heat transport in porous media,
Chin J Rock Mech Eng 2014;33(1):43e51 (in Chinese). following Lauwerier//SPE annual technical conference and exhibition. Society
[26] Chen BG, Song EX, Cheng XH. Plane-symmetrical simulation of flow and heat of Petroleum Engineers; 2010.
transport in fractured geological media: a discrete fracture model with [33] Bai B. One-dimensional thermal consolidation characteristics of geotechnical
comsol//multiphysical testing of soils and shales. Berlin Heidelberg: Springer; media under non-isothermal condition. Eng Mech 2005;22(5):186e91 (in
2013. p. 149e54. Chinese).
[27] Dayan DG. Drilling fluid design for geothermal wells. United Nations Uni- [34] McGuinness MJ, Blakeley M, Pruess K, O'sullivan MJ. Geothermal heat pipe
versity; 2014. Geothermal Training Programme. stability: solution selection by upstreaming and boundary conditions. Transp
[28] Rice JR. Fault stress states, pore pressure distributions, and the weakness of Porous Media 1993;11(1):71e100.
the San Andreas fault. Int Geophys 1992;51:475e503.