# A Thermal-Hydraulic-Mechanical Coupling Study of Heat Extraction from the Geothermal Reservoir with a Discrete Fracture Network

Authors: Zhiwei Ye; J. G. Wang
Journal: Geofluids
 (2020)

---

## Abstract

The complex thermal-hydraulic-mechanical (THM) coupling is the key issue to the energy extraction from a geothermal reservoir, where fractures are the main channels for fluid circulation and heat transfer. However, the effects of matrix deformation-induced aperture variation and fracture roughness on heat recovery efficiency are unclear. In this paper, a fully coupling THM model based on a discrete fracture network is proposed to explore these coupling effects. First, the fracture roughness and the fracture aperture variation with effective stress are introduced. Second, the water flow and heat transfer in the matrix and fractures as well as the deformation of the geothermal reservoir are individually formulated for a fractured geothermal reservoir. Third, the model is validated with analytical solution for its thermal-hydraulic (TH) coupling effect and literature data for its hydraulic-mechanical (HM) coupling effect. Finally, the features of heat transfer and fluid flow in the fractured geothermal reservoir are comparatively analyzed through four scenarios. The simulation results indicate that the discrete fracture network severely impacts the pressure distribution and temperature advance. The aperture variation induced by solid deformation can enhance heat transfer efficiency, and the fracture roughness can reduce the heat transfer efficiency.


---

## Body

## 1. Introduction

A general operation method for energy extraction from geothermal reservoirs is to inject cold fluid in and to pump hot fluid out [1]. The fluid circulation and heat exchange are the key processes in energy extraction. The rock matrix in geothermal reservoirs is usually dense and nearly impervious; thus, natural or artificial fractures or joints are the main channels for fluid circulation and heat exchange [2–4]. This energy extraction process involves the interactions among thermal (T), hydro (H), and mechanical (M) processes [5, 6]. The injection and pump of working fluid will redistribute in situ stress, resulting in the variation of fracture aperture and directly impacting the transmissivity of the fractured geothermal reservoir. This may affect the energy extraction efficiency. Therefore, the THM coupling in a fractured geothermal reservoir is the key issue to the heat extraction performance of this reservoir.Different models have been proposed for the THM coupling in an enhanced geothermal system (EGS). First, TH coupling models were developed for the forced convection in geothermal exploitations. They investigated the heterogeneity of porous media [7, 8], the system parameter optimization [9, 10], the heat exchange mechanism [11–13], and the effects of reservoir fracture treatments [14]. These TH coupling models can well evaluate the heat extraction in a rigid reservoir because they ignore the impact of reservoir deformation in the heat extraction process. Second, THM coupling models were developed based on continuum porous media [15–18]. Recent THM coupling models for geothermal reservoirs consider the impact of several isolated discrete fractures on heat extraction performance [19–21]. The fluid circulation in a complex fracture network is still unknown. The fluid circulation and heat transfer in a complex fracture network may heavily affect the extraction efficiency from a fractured geothermal reservoir, but such a study has rarely been observed in public literature.A fractured geothermal reservoir can be modeled by either continuum-based models (CM) or discrete fracture models (DFM) [22, 23]. In the CM, the fractured reservoir is treated as a continuum composed of the fractures and matrix. Typical continuum-based models include an equivalent porous medium model [11, 12, 17] and dual-porosity model [24, 25]. These continuum-based models cannot explicitly express the fractures and thus cannot include the details of fracture characteristics such as location, length, and orientation. In the DFM, the geothermal reservoir is composed of the discrete fracture network (DFN) and rock matrix [6, 26, 27]. These DFM models consider the seepage properties of fractures in a different way from the rock matrix. The complex features of discrete fractures can be individually considered within geothermal reservoirs. However, the complexity of the fracture network largely increases the difficulty in both fracture description and computational cost. A real geothermal reservoir is fractured, and the properties of fractures may play a key role in energy extraction. Thus, the THM coupling analysis is necessary to include the discrete fracture network in the fractured geothermal reservoir.The heat extraction performance of a geothermal system has been evaluated based on DFN models [13, 14, 28]. These models mainly focused on the heat exchange process but overlooked the effect of solid deformation and the fracture roughness. The impacts of aperture evolution and fracture roughness on heat extraction efficiency are still unknown. Hence, a DFN-based THM coupling model is necessary. This model can include the aperture evolution and the fracture roughness in the heat extraction from a geothermal system.This paper develops a THM coupling model to consider the effects of aperture variation and fracture roughness for the heat energy extraction from a geothermal reservoir. Firstly, the discrete fracture network is generated through the Monte Carlo method. The fracture characteristics such as location, length, and orientation are included. Secondly, the aperture variation is described by a function of fracture aperture with effective stress. Further, a roughness factor is introduced into Darcy’s law to describe the effect of fracture roughness on fluid flow in a fracture. Thirdly, a full THM coupling model is established for the heat extraction process from a fractured geothermal reservoir and is validated by two degraded submodels. Finally, the distribution of pressure and temperature in the fractured reservoir is numerically simulated. Their coupling effects are comparatively analyzed.


## 2. Generation of the Fracture Network

### 2.1. Monte Carlo Method for Fracture Generation

A fracture network is constructed by randomly generating fractures through the algorithm of the Monte Carlo method. The trace length, orientation, and location of fractures are the three key parameters to a two-dimensional fracture network. Compared with trace length, fracture aperture is too tiny to be expressed on the graph of the fracture network although it severely impacts the fluid seepage. These three key parameters are determined through the following procedure.


### 2.2. Fracture Trace Lengths

Previous studies have shown that the trace length of fractures in geothermal reservoirs obeys a power-law distribution [29, 30]. This study expresses the trace length of each fracture by
(1)l=lmin−D+Fflmax−D−lmin−D−1/D,where lmin and lmax are the minimum and maximum trace lengths, D is the fractal dimension of length distribution of fractures, and Ff is a random number uniformly distributed in the range of 0 to 1.


### 2.3. Orientations of Fractures

The fractures in the geothermal reservoir can be grouped into different sets according to their orientations [29]. A fracture network may have one, two, or three sets of fractures. Each set of fractures has a mean orientation angle γμ and a standard deviation γσ. The orientation angle of fractures follows a Gaussian distribution. This fracture orientation angle after the Box-Muller transform can be generated by
(2)γ=γμ+γσ−2lnF1cos2πF2,where γ is a fracture orientation angle in a set of fractures and F1 and F2 are random numbers uniformly distributed in the range of 0 to 1.


### 2.4. Location of Fractures

The midpoint of a fracture is used to indicate the location of a fracture. The coordinates of the fracture center is fixed by generating random numbers based on uniform distribution. The fracture centers are sequentially generated. If one fracture is too close to other fractures (<2 m in the domain of300m×300m), the fracture center will be regenerated. The coordinate generation function for the midpoints of fractures is expressed as
(3)x=LFf,y=WFf,where L and W are the length and width of the generation domain, respectively.


### 2.5. Generated Fracture Networks

Two types of fracture networks are generated in Figure1. Table 1 lists the parameters used for the generation of these DFNs. Type (a) consists of a set of horizontal fractures and a set of vertical fractures. The fracture density of (a-1) is much larger than that of (a-2), and the connectivity of (a-1) is better than that of (a-2). Type (b) has two sets of oblique fractures, and each set has a fixed orientation. The fracture density of (b-2) is half of that of (b-1). All the four DFNs are imported into the COMSOL Multiphysics for numerical simulation. The simulation results of the four DFNs are discussed in Section 5 of this paper.Figure 1
Four fracture networks generated by the MATLAB code.Table 1
Geometric parameters used for DFN generation.
DFNParameterFracture lengthFracture numberOrientationCenter pointDistributionPower-law distributionDesignatedGaussian distributionUniform distribution(a-1)Set 140-200 m1200°/02300×300mSet 220-150 m15090°/02300×300m(a-2)Set 140-200 m600°/02300×300mSet 220-150 m7590°/02300×300m(b-1)Set 140-200 m12020°/02300×300mSet 220-150 m150130°/02300×300m(b-2)Set 140-200 m6020°/02300×300mSet 220-150 m75130°/02300×300m


## 2.1. Monte Carlo Method for Fracture Generation

A fracture network is constructed by randomly generating fractures through the algorithm of the Monte Carlo method. The trace length, orientation, and location of fractures are the three key parameters to a two-dimensional fracture network. Compared with trace length, fracture aperture is too tiny to be expressed on the graph of the fracture network although it severely impacts the fluid seepage. These three key parameters are determined through the following procedure.


## 2.2. Fracture Trace Lengths

Previous studies have shown that the trace length of fractures in geothermal reservoirs obeys a power-law distribution [29, 30]. This study expresses the trace length of each fracture by
(1)l=lmin−D+Fflmax−D−lmin−D−1/D,where lmin and lmax are the minimum and maximum trace lengths, D is the fractal dimension of length distribution of fractures, and Ff is a random number uniformly distributed in the range of 0 to 1.


## 2.3. Orientations of Fractures

The fractures in the geothermal reservoir can be grouped into different sets according to their orientations [29]. A fracture network may have one, two, or three sets of fractures. Each set of fractures has a mean orientation angle γμ and a standard deviation γσ. The orientation angle of fractures follows a Gaussian distribution. This fracture orientation angle after the Box-Muller transform can be generated by
(2)γ=γμ+γσ−2lnF1cos2πF2,where γ is a fracture orientation angle in a set of fractures and F1 and F2 are random numbers uniformly distributed in the range of 0 to 1.


## 2.4. Location of Fractures

The midpoint of a fracture is used to indicate the location of a fracture. The coordinates of the fracture center is fixed by generating random numbers based on uniform distribution. The fracture centers are sequentially generated. If one fracture is too close to other fractures (<2 m in the domain of300m×300m), the fracture center will be regenerated. The coordinate generation function for the midpoints of fractures is expressed as
(3)x=LFf,y=WFf,where L and W are the length and width of the generation domain, respectively.


## 2.5. Generated Fracture Networks

Two types of fracture networks are generated in Figure1. Table 1 lists the parameters used for the generation of these DFNs. Type (a) consists of a set of horizontal fractures and a set of vertical fractures. The fracture density of (a-1) is much larger than that of (a-2), and the connectivity of (a-1) is better than that of (a-2). Type (b) has two sets of oblique fractures, and each set has a fixed orientation. The fracture density of (b-2) is half of that of (b-1). All the four DFNs are imported into the COMSOL Multiphysics for numerical simulation. The simulation results of the four DFNs are discussed in Section 5 of this paper.Figure 1
Four fracture networks generated by the MATLAB code.Table 1
Geometric parameters used for DFN generation.
DFNParameterFracture lengthFracture numberOrientationCenter pointDistributionPower-law distributionDesignatedGaussian distributionUniform distribution(a-1)Set 140-200 m1200°/02300×300mSet 220-150 m15090°/02300×300m(a-2)Set 140-200 m600°/02300×300mSet 220-150 m7590°/02300×300m(b-1)Set 140-200 m12020°/02300×300mSet 220-150 m150130°/02300×300m(b-2)Set 140-200 m6020°/02300×300mSet 220-150 m75130°/02300×300m


## 3. Coupling Model for Energy Extraction from a Geothermal Reservoir

### 3.1. Basic Assumptions

This fully coupling model has five separate yet interacting governing equations: a thermo-elastic mechanical equation, two flow equations (for the fractures and rock matrix), two heat transfer equations in the fractures and rock matrix. This model is mainly based on the following assumptions: (1) The rock matrix is continuous, isotropic, and homogeneous. (2) Surrounding rocks are impermeable. The fluid loss is not considered. (3) Both the matrix and fractures are saturated with single-phase liquid. (4) The fluid flow in both the matrix and fractures follows Darcy’s law. (5) The local thermal equilibrium is assumed at the interface between the rock matrix and the working fluid.


### 3.2. Governing Equations for Fluid Flow

#### 3.2.1. Fluid Flow in the Matrix

In the saturated porous matrix, the mass conservation of fluid is(4)∂∂tρfϕ+∇⋅ρfus=Qs,where ρf is the fluid density, ϕ is the porosity of the matrix, t denotes time, us is the flow velocity of the Darcy seepage, and Qs is the source or sink term.The seepage velocityus is expressed by Darcy’s law as
(5)us=−kμ∇p.Furthermore, the storage term at the first term of Equation (4) is
(6)∂∂tρfϕ=ρfS∂p∂t.The porosity dilation-induced fluid variation is treated as the source termQs which is shown as
(7)Qs=−ρfα∂εv∂t,where k is the intrinsic permeability of the rock matrix, μ is the dynamic viscosity of water, S is the storage coefficient of the porous matrix, p is the water pressure, αis Biot’s coefficient of the porous matrix, and εv is the volumetric strain.Thus, the governing equation of water flow in the matrix is(8)ρfS∂p∂t+ρfα∂εv∂t+∇⋅−ρfkμ∇p=0.


#### 3.2.2. Fluid Flow in Fractures

The continuity equation of fluid flow in discrete fractures is(9)b∂ρf∂t+∇T⋅bρfuf=ρfQf,where b is the fracture aperture, uf is the flow velocity in the fracture, ∇T denotes the gradient operator restricted to the tangential plane of a fracture, and Qf is the source or sink term.A fracture consists of two smooth parallel plates; thus, the volumetric flow rateQ is
(10)Q=ufwb=−wb312μ∇Tp,and the permeability of a single fracture is then obtained as
(11)kf=b212,where w is the fracture width and b is the fracture aperture. ∇Tp is the pressure gradient along the fracture.We introducefσ of the roughness factor to consider the effect of fracture roughness on flow:
(12)kf=b212fσ.The following empirical expression is used in our computation:(13)fσ=1+aσrb1.5,where σr is the standard deviation of the fracture surface height and a is a constant (a=17 by Lomize [31] and a=8.8 by Louis [32]). Single fracture flow found that fσ varies from 1.04 to 1.65 [33]. Hence, the fracture flow velocity uf is
(14)uf=−b212μfσ∇Tp.In this study, fractures are treated as the internal boundaries of the matrix seepage model. The fluid pressure in fractures is the Dirichlet boundary condition of fluid flow in the matrix, while the pore pressure of the matrix is the motivation of the source/sink term of fracture flow. The source or sink term in the governing equation for fracture flow is(15)Qf=−kμ∂p∂nu−kμ∂p∂nd,where nu and nd are the normal direction of up and down surfaces of the fracture, respectively.Finally, the governing equation for fracture flow is obtained as(16)b∂ρf∂t+∇T⋅−ρfb312μfσ∇Tp=−ρfkμ∂p∂nu−ρfkμ∂p∂nd.


### 3.3. Governing Equations for Heat Transfer in the Matrix and Fractures

As a kind of porous media, the rock matrix has a higher specific surface area. Fluid can exchange heat with the matrix instantly. Hence, the temperature of pore water is assumed to be identical to that of the rock matrix. Then, the heat transfer equation in the matrix can be written as(17)ρCeff∂T∂t+ρfCfus⋅∇T+∇⋅q=0,where ρCeff=1−ϕρsCs+ϕρfCf is the effective heat capacity, ϕ is the porosity of the rock matrix, ρs is the matrix density, Cs is the specific heat capacity of the matrix, Cf is the specific heat capacity of fluid, T is the matrix temperature, and q is the heat flux in the rock matrix.The heat fluxq can be expressed through Fourier’s law as
(18)q=−λeff∇T,where λeff=1−ϕλs+ϕλf is the effective heat conductivity, λs is the heat conductivity of the rock matrix, and λf is the heat conductivity of fluid.The heat transfer equation in discrete fractures is(19)bρfCf∂T∂t+bρfCfuf⋅∇TT−∇Tbqf=QTf,where qf is the heat flux in the fracture and is
(20)qf=−λf∇TT.QTf is the source or sink term in fracture heat transfer. It is the heat exchange between the matrix and fracture. Based on the local thermal equilibrium assumption, QTf can be expressed as
(21)QTf=−ρfCfkμ∂p∂nuT−λeff∂T∂nu−ρfCfkμ∂p∂ndT−λeff∂T∂nd.


### 3.4. Mechanical Equilibrium

The constitutive equation of poroelasticity for THM coupling is(22)σ=2Gε+2Gν1−2νε:II−αpI−KαTT−T0I.The Navier equation for displacementv can be expressed as
(23)G∇⋅∇v+G1−2ν∇∇⋅v−α∇p−KαT∇T+F=0,where σ is the stress tensor (positive for tension), G is the shear modulus, ε is the strain tensor, ν is Poisson’s ratio, and I is the unit tensor. α is Biot’s coefficient, which depends on the compressibility; p is the pore pressure (negative for suction); K=2G1+ν/31−2ν=E/31−2ν is the bulk modulus of porous media, and E is Young’s modulus. αT is the coefficient of thermal volumetric expansion of porous media. T0 is the initial temperature of the rock mass. F is the body force.


### 3.5. Multiphysical Couplings

#### 3.5.1. Variation of Fracture Aperture with Effective Normal Stress

A geothermal reservoir is composed of a rock matrix and fracture network. The rock matrix is usually of low porosity and low permeability. The slight change of porosity has little effect on water flow and heat transfer; thus, the rock matrix is assumed to have constant porosity. The fracture network is the primary channel for fluid flow. Any change of fracture aperture may severely affect water flow. The fracture aperture is related to effective normal stress as [34]
(24)δ=σnδmσn+kn0δm,where δ=b0−b is the fracture closure, where b0 is the initial aperture with zero effective normal stress; σn is the normal component of the effective stress tensor σ+αpI; δm is the maximum possible fracture closure; and kn0 is an empirical parameter representing the initial normal stiffness. δm and kn0 can be determined with the tested samples under the specific testing conditions. Baghbanan and Jing [35] found that the above two parameters are difficultly evaluated for a set of fractures. They presented some empirical expressions for use: σncMPa=10kn0δm=487b0mm+2.51,δm=0.9b0.


#### 3.5.2. Change of Fluid Property with Temperature

The fluid property includes water densityρf and dynamic viscosity μ and affects water flow and heat transfer. It is found that pressure change does not significantly modify water density, but temperature severely modifies both water density and dynamic viscosity. A piecewise function is used for water density and dynamic viscosity as [36]
(25)ρf=1000×1−T−3.982503570×T+283T+67.26,0°C≤T≤20°C,996.9×1−3.17×10−4×T−25−2.56×10−6×T−252,20°C<T≤175°C,1758.4+103×T+273.15−4.84×10−3+T+273.151.09×10−5−T+273.15×9.85×10−9,175°C<T≤300°C,μ=exp−0.033+1.962×10−4×T×T×1.787×10−3,0°C≤T≤40°C,1+0.015512×T−20−1.572×10−3,40°C<T≤100°C,0.2414×10247.8/T+133.15×10−4,100°C<T≤300°C.


## 3.1. Basic Assumptions

This fully coupling model has five separate yet interacting governing equations: a thermo-elastic mechanical equation, two flow equations (for the fractures and rock matrix), two heat transfer equations in the fractures and rock matrix. This model is mainly based on the following assumptions: (1) The rock matrix is continuous, isotropic, and homogeneous. (2) Surrounding rocks are impermeable. The fluid loss is not considered. (3) Both the matrix and fractures are saturated with single-phase liquid. (4) The fluid flow in both the matrix and fractures follows Darcy’s law. (5) The local thermal equilibrium is assumed at the interface between the rock matrix and the working fluid.


## 3.2. Governing Equations for Fluid Flow

### 3.2.1. Fluid Flow in the Matrix

In the saturated porous matrix, the mass conservation of fluid is(4)∂∂tρfϕ+∇⋅ρfus=Qs,where ρf is the fluid density, ϕ is the porosity of the matrix, t denotes time, us is the flow velocity of the Darcy seepage, and Qs is the source or sink term.The seepage velocityus is expressed by Darcy’s law as
(5)us=−kμ∇p.Furthermore, the storage term at the first term of Equation (4) is
(6)∂∂tρfϕ=ρfS∂p∂t.The porosity dilation-induced fluid variation is treated as the source termQs which is shown as
(7)Qs=−ρfα∂εv∂t,where k is the intrinsic permeability of the rock matrix, μ is the dynamic viscosity of water, S is the storage coefficient of the porous matrix, p is the water pressure, αis Biot’s coefficient of the porous matrix, and εv is the volumetric strain.Thus, the governing equation of water flow in the matrix is(8)ρfS∂p∂t+ρfα∂εv∂t+∇⋅−ρfkμ∇p=0.


### 3.2.2. Fluid Flow in Fractures

The continuity equation of fluid flow in discrete fractures is(9)b∂ρf∂t+∇T⋅bρfuf=ρfQf,where b is the fracture aperture, uf is the flow velocity in the fracture, ∇T denotes the gradient operator restricted to the tangential plane of a fracture, and Qf is the source or sink term.A fracture consists of two smooth parallel plates; thus, the volumetric flow rateQ is
(10)Q=ufwb=−wb312μ∇Tp,and the permeability of a single fracture is then obtained as
(11)kf=b212,where w is the fracture width and b is the fracture aperture. ∇Tp is the pressure gradient along the fracture.We introducefσ of the roughness factor to consider the effect of fracture roughness on flow:
(12)kf=b212fσ.The following empirical expression is used in our computation:(13)fσ=1+aσrb1.5,where σr is the standard deviation of the fracture surface height and a is a constant (a=17 by Lomize [31] and a=8.8 by Louis [32]). Single fracture flow found that fσ varies from 1.04 to 1.65 [33]. Hence, the fracture flow velocity uf is
(14)uf=−b212μfσ∇Tp.In this study, fractures are treated as the internal boundaries of the matrix seepage model. The fluid pressure in fractures is the Dirichlet boundary condition of fluid flow in the matrix, while the pore pressure of the matrix is the motivation of the source/sink term of fracture flow. The source or sink term in the governing equation for fracture flow is(15)Qf=−kμ∂p∂nu−kμ∂p∂nd,where nu and nd are the normal direction of up and down surfaces of the fracture, respectively.Finally, the governing equation for fracture flow is obtained as(16)b∂ρf∂t+∇T⋅−ρfb312μfσ∇Tp=−ρfkμ∂p∂nu−ρfkμ∂p∂nd.


## 3.2.1. Fluid Flow in the Matrix

In the saturated porous matrix, the mass conservation of fluid is(4)∂∂tρfϕ+∇⋅ρfus=Qs,where ρf is the fluid density, ϕ is the porosity of the matrix, t denotes time, us is the flow velocity of the Darcy seepage, and Qs is the source or sink term.The seepage velocityus is expressed by Darcy’s law as
(5)us=−kμ∇p.Furthermore, the storage term at the first term of Equation (4) is
(6)∂∂tρfϕ=ρfS∂p∂t.The porosity dilation-induced fluid variation is treated as the source termQs which is shown as
(7)Qs=−ρfα∂εv∂t,where k is the intrinsic permeability of the rock matrix, μ is the dynamic viscosity of water, S is the storage coefficient of the porous matrix, p is the water pressure, αis Biot’s coefficient of the porous matrix, and εv is the volumetric strain.Thus, the governing equation of water flow in the matrix is(8)ρfS∂p∂t+ρfα∂εv∂t+∇⋅−ρfkμ∇p=0.


## 3.2.2. Fluid Flow in Fractures

The continuity equation of fluid flow in discrete fractures is(9)b∂ρf∂t+∇T⋅bρfuf=ρfQf,where b is the fracture aperture, uf is the flow velocity in the fracture, ∇T denotes the gradient operator restricted to the tangential plane of a fracture, and Qf is the source or sink term.A fracture consists of two smooth parallel plates; thus, the volumetric flow rateQ is
(10)Q=ufwb=−wb312μ∇Tp,and the permeability of a single fracture is then obtained as
(11)kf=b212,where w is the fracture width and b is the fracture aperture. ∇Tp is the pressure gradient along the fracture.We introducefσ of the roughness factor to consider the effect of fracture roughness on flow:
(12)kf=b212fσ.The following empirical expression is used in our computation:(13)fσ=1+aσrb1.5,where σr is the standard deviation of the fracture surface height and a is a constant (a=17 by Lomize [31] and a=8.8 by Louis [32]). Single fracture flow found that fσ varies from 1.04 to 1.65 [33]. Hence, the fracture flow velocity uf is
(14)uf=−b212μfσ∇Tp.In this study, fractures are treated as the internal boundaries of the matrix seepage model. The fluid pressure in fractures is the Dirichlet boundary condition of fluid flow in the matrix, while the pore pressure of the matrix is the motivation of the source/sink term of fracture flow. The source or sink term in the governing equation for fracture flow is(15)Qf=−kμ∂p∂nu−kμ∂p∂nd,where nu and nd are the normal direction of up and down surfaces of the fracture, respectively.Finally, the governing equation for fracture flow is obtained as(16)b∂ρf∂t+∇T⋅−ρfb312μfσ∇Tp=−ρfkμ∂p∂nu−ρfkμ∂p∂nd.


## 3.3. Governing Equations for Heat Transfer in the Matrix and Fractures

As a kind of porous media, the rock matrix has a higher specific surface area. Fluid can exchange heat with the matrix instantly. Hence, the temperature of pore water is assumed to be identical to that of the rock matrix. Then, the heat transfer equation in the matrix can be written as(17)ρCeff∂T∂t+ρfCfus⋅∇T+∇⋅q=0,where ρCeff=1−ϕρsCs+ϕρfCf is the effective heat capacity, ϕ is the porosity of the rock matrix, ρs is the matrix density, Cs is the specific heat capacity of the matrix, Cf is the specific heat capacity of fluid, T is the matrix temperature, and q is the heat flux in the rock matrix.The heat fluxq can be expressed through Fourier’s law as
(18)q=−λeff∇T,where λeff=1−ϕλs+ϕλf is the effective heat conductivity, λs is the heat conductivity of the rock matrix, and λf is the heat conductivity of fluid.The heat transfer equation in discrete fractures is(19)bρfCf∂T∂t+bρfCfuf⋅∇TT−∇Tbqf=QTf,where qf is the heat flux in the fracture and is
(20)qf=−λf∇TT.QTf is the source or sink term in fracture heat transfer. It is the heat exchange between the matrix and fracture. Based on the local thermal equilibrium assumption, QTf can be expressed as
(21)QTf=−ρfCfkμ∂p∂nuT−λeff∂T∂nu−ρfCfkμ∂p∂ndT−λeff∂T∂nd.


## 3.4. Mechanical Equilibrium

The constitutive equation of poroelasticity for THM coupling is(22)σ=2Gε+2Gν1−2νε:II−αpI−KαTT−T0I.The Navier equation for displacementv can be expressed as
(23)G∇⋅∇v+G1−2ν∇∇⋅v−α∇p−KαT∇T+F=0,where σ is the stress tensor (positive for tension), G is the shear modulus, ε is the strain tensor, ν is Poisson’s ratio, and I is the unit tensor. α is Biot’s coefficient, which depends on the compressibility; p is the pore pressure (negative for suction); K=2G1+ν/31−2ν=E/31−2ν is the bulk modulus of porous media, and E is Young’s modulus. αT is the coefficient of thermal volumetric expansion of porous media. T0 is the initial temperature of the rock mass. F is the body force.


## 3.5. Multiphysical Couplings

### 3.5.1. Variation of Fracture Aperture with Effective Normal Stress

A geothermal reservoir is composed of a rock matrix and fracture network. The rock matrix is usually of low porosity and low permeability. The slight change of porosity has little effect on water flow and heat transfer; thus, the rock matrix is assumed to have constant porosity. The fracture network is the primary channel for fluid flow. Any change of fracture aperture may severely affect water flow. The fracture aperture is related to effective normal stress as [34]
(24)δ=σnδmσn+kn0δm,where δ=b0−b is the fracture closure, where b0 is the initial aperture with zero effective normal stress; σn is the normal component of the effective stress tensor σ+αpI; δm is the maximum possible fracture closure; and kn0 is an empirical parameter representing the initial normal stiffness. δm and kn0 can be determined with the tested samples under the specific testing conditions. Baghbanan and Jing [35] found that the above two parameters are difficultly evaluated for a set of fractures. They presented some empirical expressions for use: σncMPa=10kn0δm=487b0mm+2.51,δm=0.9b0.


### 3.5.2. Change of Fluid Property with Temperature

The fluid property includes water densityρf and dynamic viscosity μ and affects water flow and heat transfer. It is found that pressure change does not significantly modify water density, but temperature severely modifies both water density and dynamic viscosity. A piecewise function is used for water density and dynamic viscosity as [36]
(25)ρf=1000×1−T−3.982503570×T+283T+67.26,0°C≤T≤20°C,996.9×1−3.17×10−4×T−25−2.56×10−6×T−252,20°C<T≤175°C,1758.4+103×T+273.15−4.84×10−3+T+273.151.09×10−5−T+273.15×9.85×10−9,175°C<T≤300°C,μ=exp−0.033+1.962×10−4×T×T×1.787×10−3,0°C≤T≤40°C,1+0.015512×T−20−1.572×10−3,40°C<T≤100°C,0.2414×10247.8/T+133.15×10−4,100°C<T≤300°C.


## 3.5.1. Variation of Fracture Aperture with Effective Normal Stress

A geothermal reservoir is composed of a rock matrix and fracture network. The rock matrix is usually of low porosity and low permeability. The slight change of porosity has little effect on water flow and heat transfer; thus, the rock matrix is assumed to have constant porosity. The fracture network is the primary channel for fluid flow. Any change of fracture aperture may severely affect water flow. The fracture aperture is related to effective normal stress as [34]
(24)δ=σnδmσn+kn0δm,where δ=b0−b is the fracture closure, where b0 is the initial aperture with zero effective normal stress; σn is the normal component of the effective stress tensor σ+αpI; δm is the maximum possible fracture closure; and kn0 is an empirical parameter representing the initial normal stiffness. δm and kn0 can be determined with the tested samples under the specific testing conditions. Baghbanan and Jing [35] found that the above two parameters are difficultly evaluated for a set of fractures. They presented some empirical expressions for use: σncMPa=10kn0δm=487b0mm+2.51,δm=0.9b0.


## 3.5.2. Change of Fluid Property with Temperature

The fluid property includes water densityρf and dynamic viscosity μ and affects water flow and heat transfer. It is found that pressure change does not significantly modify water density, but temperature severely modifies both water density and dynamic viscosity. A piecewise function is used for water density and dynamic viscosity as [36]
(25)ρf=1000×1−T−3.982503570×T+283T+67.26,0°C≤T≤20°C,996.9×1−3.17×10−4×T−25−2.56×10−6×T−252,20°C<T≤175°C,1758.4+103×T+273.15−4.84×10−3+T+273.151.09×10−5−T+273.15×9.85×10−9,175°C<T≤300°C,μ=exp−0.033+1.962×10−4×T×T×1.787×10−3,0°C≤T≤40°C,1+0.015512×T−20−1.572×10−3,40°C<T≤100°C,0.2414×10247.8/T+133.15×10−4,100°C<T≤300°C.


## 4. Simulation Implementation and Model Validation

### 4.1. Numerical Simulation Procedure

These governing equations with the initial and boundary conditions formulate a complete THM coupling model. Generally, these partial differential equations are too complex to derive analytical solutions. This study uses a finite element method (FEM) to solve this problem. Particularly, the FEM solutions are numerically obtained through the COMSOL Multiphysics, a powerful solver for partial differential equations. The solution procedure has four steps. Firstly, rock fractures are generated with the Monte Carlo algorithm. Then, the geometric model of fractured rock mass is imported into COMSOL Multiphysics. Thirdly, both governing equations and boundary conditions are assigned to the geometric model, and model parameters and initial conditions are set. Finally, the distributions of temperature and pressure are solved and visualized.


### 4.2. Model Parameters

A realistic geothermal reservoir usually has a three-dimensional (3D) fracture network. A complex 3D fracture network is extremely difficult in implementation and largely increases computation cost. In this study, a 2D fracture network in Figure2 is used to model a geothermal reservoir. This geothermal reservoir is at 3500 m deep and has a computational domain of 300m×300m. The simulation time is 50 years from the beginning of cool water injection. Working fluid of water is injected into the geothermal reservoir. The left side boundary of the reservoir is the injection well, and the right side boundary is the production well. The initial and boundary conditions are described below.Figure 2
FEM mesh model of the geothermal reservoir with initial and boundary conditions.For the reservoir seepage, the hydraulic pressure is kept as 80 MPa at the injection well and as 66 MPa at the right boundary. The top and bottom boundaries are impermeable. The initial pressure of this reservoir is 70 MPa. For the heat transfer in the reservoir, the water temperature is 30°C at the injection well. The top and bottom boundaries are thermally insulated. The initial temperature of the reservoir is 180°C. For the mechanical equilibrium model, the top and right boundaries are applied a constant normal compression of 90 MPa, which is equivalent to the overburden stress at 3500 m deep. The other two boundaries have zero normal displacement and free tangential displacement. Table2 lists the model parameters used in simulations.Table 2
Parameters of the coupling model.
VariableParameterValueUnitρsRock matrix density2700kg·m-3CfSpecific heat capacity of water4200J·kg-1·K-1CsSpecific heat capacity of the rock matrix1000J·kg-1·K-1λfThermal conductivity of water0.7W·m-1·K-1λsThermal conductivity of the rock matrix3W·m-1·K-1ϕPorosity0.051kMatrix permeability1.0×10−18m2EYoung’s modulus30GPaνPoisson’s ratio0.251αBiot’s coefficient11αTCoefficient of volumetric expansion5.0×10−61/KSStorage coefficient of the matrix1.0×10−81/Pab0Initial aperture of fractures0.18mmσrStandard deviation of the fracture surface height0.005mm


### 4.3. Model Validation

This fully coupling THM model is divided into two submodels: a TH coupling model for the heat transfer in a single fracture and an HM coupling model for the variation of fracture aperture. Figure3 is the geometry model for the heat transfer in a single fracture. Cool water is injected into the fracture and heated by the impermeable rock matrix. For a problem similar to Figure 3 but with different boundary conditions, Ye and Wang [37] derived an analytical solution through the Laplace transform and Stehfest method. Following the same method, the temperature distribution of heat transfer in a single fracture can be obtained analytically. The deduction process is presented in the appendix, which is used for model validation.Figure 3
Geometric model of heat transfer in a single fractured rock.The distribution of fracture temperature and the observation line at the 100th day are plotted in Figure4(a), and the variation of temperature at the fixed points is plotted in Figure 4(b). The temperature distribution obtained from the current numerical simulation is in good agreement with the analytical solution. Only small deviation is observed at the observation point, but within the acceptable range.Figure 4
Comparison between analytical solution and present simulation.
(a)
Temperature at the 100th day(b)
Temperature at fixed pointsFigure5 presents a geometry model with two horizontal wells and four fractures. The domain size and boundary conditions are marked. The two wells are located at the left and right boundaries. The four fractures are located inside the domain: fracture (a) has a length of 200 m and is connected with two wellbores, fracture (b) has a length of 100 m and is only connected with the injection wellbore, fracture (c) has a length of 100 m and is only connected with the production wellbore, and fracture (d) has a length of 100 m and is isolated from the two wellbores.Figure 5
Fracture distribution diagram between two horizontal wells.Under the same parameters, both the simulation results in this study and those from Moradi et al. [38] are plotted in Figure 6. Generally, both solutions are in good agreement. Our numerical model can correctly describe the variation of fracture aperture induced by the HM coupling effect.Figure 6
Comparison of fracture apertures at different locations.
(a)(b)(c)(d)


## 4.1. Numerical Simulation Procedure

These governing equations with the initial and boundary conditions formulate a complete THM coupling model. Generally, these partial differential equations are too complex to derive analytical solutions. This study uses a finite element method (FEM) to solve this problem. Particularly, the FEM solutions are numerically obtained through the COMSOL Multiphysics, a powerful solver for partial differential equations. The solution procedure has four steps. Firstly, rock fractures are generated with the Monte Carlo algorithm. Then, the geometric model of fractured rock mass is imported into COMSOL Multiphysics. Thirdly, both governing equations and boundary conditions are assigned to the geometric model, and model parameters and initial conditions are set. Finally, the distributions of temperature and pressure are solved and visualized.


## 4.2. Model Parameters

A realistic geothermal reservoir usually has a three-dimensional (3D) fracture network. A complex 3D fracture network is extremely difficult in implementation and largely increases computation cost. In this study, a 2D fracture network in Figure2 is used to model a geothermal reservoir. This geothermal reservoir is at 3500 m deep and has a computational domain of 300m×300m. The simulation time is 50 years from the beginning of cool water injection. Working fluid of water is injected into the geothermal reservoir. The left side boundary of the reservoir is the injection well, and the right side boundary is the production well. The initial and boundary conditions are described below.Figure 2
FEM mesh model of the geothermal reservoir with initial and boundary conditions.For the reservoir seepage, the hydraulic pressure is kept as 80 MPa at the injection well and as 66 MPa at the right boundary. The top and bottom boundaries are impermeable. The initial pressure of this reservoir is 70 MPa. For the heat transfer in the reservoir, the water temperature is 30°C at the injection well. The top and bottom boundaries are thermally insulated. The initial temperature of the reservoir is 180°C. For the mechanical equilibrium model, the top and right boundaries are applied a constant normal compression of 90 MPa, which is equivalent to the overburden stress at 3500 m deep. The other two boundaries have zero normal displacement and free tangential displacement. Table2 lists the model parameters used in simulations.Table 2
Parameters of the coupling model.
VariableParameterValueUnitρsRock matrix density2700kg·m-3CfSpecific heat capacity of water4200J·kg-1·K-1CsSpecific heat capacity of the rock matrix1000J·kg-1·K-1λfThermal conductivity of water0.7W·m-1·K-1λsThermal conductivity of the rock matrix3W·m-1·K-1ϕPorosity0.051kMatrix permeability1.0×10−18m2EYoung’s modulus30GPaνPoisson’s ratio0.251αBiot’s coefficient11αTCoefficient of volumetric expansion5.0×10−61/KSStorage coefficient of the matrix1.0×10−81/Pab0Initial aperture of fractures0.18mmσrStandard deviation of the fracture surface height0.005mm


## 4.3. Model Validation

This fully coupling THM model is divided into two submodels: a TH coupling model for the heat transfer in a single fracture and an HM coupling model for the variation of fracture aperture. Figure3 is the geometry model for the heat transfer in a single fracture. Cool water is injected into the fracture and heated by the impermeable rock matrix. For a problem similar to Figure 3 but with different boundary conditions, Ye and Wang [37] derived an analytical solution through the Laplace transform and Stehfest method. Following the same method, the temperature distribution of heat transfer in a single fracture can be obtained analytically. The deduction process is presented in the appendix, which is used for model validation.Figure 3
Geometric model of heat transfer in a single fractured rock.The distribution of fracture temperature and the observation line at the 100th day are plotted in Figure4(a), and the variation of temperature at the fixed points is plotted in Figure 4(b). The temperature distribution obtained from the current numerical simulation is in good agreement with the analytical solution. Only small deviation is observed at the observation point, but within the acceptable range.Figure 4
Comparison between analytical solution and present simulation.
(a)
Temperature at the 100th day(b)
Temperature at fixed pointsFigure5 presents a geometry model with two horizontal wells and four fractures. The domain size and boundary conditions are marked. The two wells are located at the left and right boundaries. The four fractures are located inside the domain: fracture (a) has a length of 200 m and is connected with two wellbores, fracture (b) has a length of 100 m and is only connected with the injection wellbore, fracture (c) has a length of 100 m and is only connected with the production wellbore, and fracture (d) has a length of 100 m and is isolated from the two wellbores.Figure 5
Fracture distribution diagram between two horizontal wells.Under the same parameters, both the simulation results in this study and those from Moradi et al. [38] are plotted in Figure 6. Generally, both solutions are in good agreement. Our numerical model can correctly describe the variation of fracture aperture induced by the HM coupling effect.Figure 6
Comparison of fracture apertures at different locations.
(a)(b)(c)(d)


## 5. Simulation Results and Discussion

### 5.1. Pore Pressure Distribution

The distribution of pore pressure has similar characteristics in the four fracture networks. Thus, only the pressure of the fracture network (a-1) is plotted for discussion. As shown in Figure7, the pressure distribution is mainly along the fractures. That is because fractures are the main channels for fluid seepage. In the initial stage of water injection (before 3 years), the fractures beside the injection well have higher pressure than the matrix around the fracture. The fractures beside the production well have lower pressure than the matrix around the fracture. When the pressure of fracture water reaches a stable state, the pressure of matrix water is still changing continuously. After 15 years of water injection, the pore pressure of the reservoir tends to be stable. The area with connected fractures has a lower pressure gradient. That is, the water pressure decreases gently (the dash circle area of Figure 7). The area with nonconnected fractures has a higher pressure gradient, and the water pressure decreases sharply (the solid circle area of Figure 7).Figure 7
Pressure distribution of the reservoir at the 3rd year and 30th year.


### 5.2. Reservoir Temperature Distribution

The temperature distributions of the four fracture networks at the 30th year are presented in Figure8. After cool water is injected into the geothermal reservoir, the water as working fluid is immediately heated by the rock matrix. At the same time, the temperature of the rock matrix decreases quickly. The flow velocity is much higher in fractures than in the matrix. Fracture temperature first decreases and then extends to the matrix. As more and more heat is extracted, the cooling area gradually expands from the injection well to the production well. However, the temperature at the production well does not decrease sharply before 30 years. This implies that the geothermal system can maintain a steady heat production for a quite long period.Figure 8
Temperature distributions of the four fracture networks at the 30th year.If the fractures form a network which connects the injection well with the production well, the cooling area will advance faster along these fractures than in other zones (the dash circle area of Figure8). Higher fracture density increases the probability of fracture intersection, resulting in better connectivity. Thus, the connectivity of (a-1) and (b-1) is better than that of (a-2) and (b-2). Actually, fracture networks (a-2) and (b-2) only have one cluster of fractures that connected both injection and production wells (dashed arrow line in Figure 8). Thus, their thermal front propagation is much slower than that of (a-1) and (b-1). The connected fractures are the main flow channels for fluid circulation. Their fracture flow flux is high, carrying away a majority of heat energy.


### 5.3. Temperature and Flow Velocity at the Outlet

The temperature and flow velocity at the production well are two important parameters to measure EGS performance. Based on the simulation results of the fracture network (a-1), the water temperature in the matrix and fractures along the production well is plotted in Figure9(a). The line markers are for the matrix, and the scatter markers are for fractures. Since the permeability of fractures is much higher than that of the matrix, the flow in fractures is much faster than that in the matrix. The flow velocity in the matrix (lines) and fractures (scatter markers) along the production well is plotted in Figure 9(b).Figure 9
Outlet temperature and flow velocity at the different times.
(a)
Temperature(b)
Flow velocityThe temperature at the production well is basically kept at about 180°C in the first 15 years of water injection. During this period, the EGS can maintain a high-efficiency heat production. After 30 years of cool water injection, the outlet temperature could remain above 140°C. Then, the cooling area gradually expands to the production well. The temperature drop along the production well is not uniform. It is faster in the zones with dense fractures where the faster flow of working fluid can more quickly transfer the thermal heat of the reservoir.Under THM coupling, the aperture of fractures is in the range of 60μm to 90 μm. Thus, the permeability in fractures is several orders of magnitude higher than that in the matrix. As shown in Figure 9(b), the outlet flow velocity in fractures (scatter markers) is about 7 orders of magnitude higher than that in the matrix (line markers). At the initial stage of heat extraction (before 1 year), the flow velocity in the matrix away from fractures is greater than that near fractures. The water in the matrix near fractures firstly flows into the fractures and continues to flow into the production well.


### 5.4. Comparison of Different Coupling Models

The TH coupling is the core for heat extraction. In the EGS operation, solid deformation and some other factors affect the TH coupling, thus modifying the heat extraction efficiency. In this section, four scenarios in the fully coupling THM model are investigated to comparatively analyze the different coupling effects on the outlet temperature and flow velocity. Table3 lists the coupling effect in each scenario.Table 3
Different coupling factors in four scenarios.
Coupling factorFormula termScenario 1 (TH)Scenario 2 THM)Scenario 3 (TH+roughness)Scenario 4 (THM+roughness)HMInitial stress equilibrium√√√√THρfCfus⋅∇TbρfCfuf⋅∇TTρfT, μT√√√√HM for the matrixρfα∂εv/∂t−α∇p✗√✗√HM for the fractureδ=σnδm/σn+kn0δm−α∇p✗√✗√TM−KαT∇T✗√✗√Fracture roughnessfσ=1+8.8σr/b1.5✗✗√√The initial stress equilibrium refers to the initial state of fluid-solid static coupling before cool water injection. Hence, the first step is to calculate the initial stress equilibrium. Then, the four models with different couplings are implemented to investigate their different coupling effects on heat transfer.Based on the simulation results, the temperature distributions of the four scenarios at 30 years are presented in Figure10. The temperature difference of the four scenarios mainly lies in the locations of thermal fronts. The thermal front of the THM model is more ahead than that of the other three models. This is mainly due to the increase in heat transfer efficiency caused by the enlargement of fracture aperture. Since fracture roughness results in the decrease in reservoir permeability, the thermal front of the TH+roughness model is in the most lag behind place.Figure 10
Temperature distributions of the four scenarios at the 30th year.The outlet temperature and flow velocity are plotted in Figure11. Compared to the TH model, the THM model with smooth fracture has higher flow velocity but lower temperature at the production well. This implies that solid deformation can enhance reservoir permeability and finally increases heat transfer efficiency. If the fracture roughness is considered, the outlet velocity decreases but the outlet temperature increases. Hence, the fracture roughness leads to a decrease in reservoir permeability and finally decreases the heat transfer efficiency.Figure 11
Outlet temperature and flow velocity under different scenarios.
(a)
Temperature 25th year(b)
Flow velocity 25th yearThe average outlet temperature is calculated by(26)Tout=∑ufbT+∫usTdy∑ufb+∫usdy.The summation term refers to the fractures while the integration term refers to the rock matrix. The average outlet temperature is a comprehensive indicator for the evaluation of EGS performance. The outlet temperatures at 25 years are plotted in Figure12 for the four scenarios. It is observed that the TH model with fracture roughness has the highest average outlet temperature. The THM model with smooth fracture has the lowest average outlet temperature. The variation of fracture aperture induced by matrix deformation leads to the enhancement of heat transfer efficiency. However, fracture roughness leads to a decrease in heat transfer efficiency. It is observed that the average outlet temperature of the fully coupling model is similar to that of the simple TH model in this case.Figure 12
Average outlet temperatures of the four scenarios.


## 5.1. Pore Pressure Distribution

The distribution of pore pressure has similar characteristics in the four fracture networks. Thus, only the pressure of the fracture network (a-1) is plotted for discussion. As shown in Figure7, the pressure distribution is mainly along the fractures. That is because fractures are the main channels for fluid seepage. In the initial stage of water injection (before 3 years), the fractures beside the injection well have higher pressure than the matrix around the fracture. The fractures beside the production well have lower pressure than the matrix around the fracture. When the pressure of fracture water reaches a stable state, the pressure of matrix water is still changing continuously. After 15 years of water injection, the pore pressure of the reservoir tends to be stable. The area with connected fractures has a lower pressure gradient. That is, the water pressure decreases gently (the dash circle area of Figure 7). The area with nonconnected fractures has a higher pressure gradient, and the water pressure decreases sharply (the solid circle area of Figure 7).Figure 7
Pressure distribution of the reservoir at the 3rd year and 30th year.


## 5.2. Reservoir Temperature Distribution

The temperature distributions of the four fracture networks at the 30th year are presented in Figure8. After cool water is injected into the geothermal reservoir, the water as working fluid is immediately heated by the rock matrix. At the same time, the temperature of the rock matrix decreases quickly. The flow velocity is much higher in fractures than in the matrix. Fracture temperature first decreases and then extends to the matrix. As more and more heat is extracted, the cooling area gradually expands from the injection well to the production well. However, the temperature at the production well does not decrease sharply before 30 years. This implies that the geothermal system can maintain a steady heat production for a quite long period.Figure 8
Temperature distributions of the four fracture networks at the 30th year.If the fractures form a network which connects the injection well with the production well, the cooling area will advance faster along these fractures than in other zones (the dash circle area of Figure8). Higher fracture density increases the probability of fracture intersection, resulting in better connectivity. Thus, the connectivity of (a-1) and (b-1) is better than that of (a-2) and (b-2). Actually, fracture networks (a-2) and (b-2) only have one cluster of fractures that connected both injection and production wells (dashed arrow line in Figure 8). Thus, their thermal front propagation is much slower than that of (a-1) and (b-1). The connected fractures are the main flow channels for fluid circulation. Their fracture flow flux is high, carrying away a majority of heat energy.


## 5.3. Temperature and Flow Velocity at the Outlet

The temperature and flow velocity at the production well are two important parameters to measure EGS performance. Based on the simulation results of the fracture network (a-1), the water temperature in the matrix and fractures along the production well is plotted in Figure9(a). The line markers are for the matrix, and the scatter markers are for fractures. Since the permeability of fractures is much higher than that of the matrix, the flow in fractures is much faster than that in the matrix. The flow velocity in the matrix (lines) and fractures (scatter markers) along the production well is plotted in Figure 9(b).Figure 9
Outlet temperature and flow velocity at the different times.
(a)
Temperature(b)
Flow velocityThe temperature at the production well is basically kept at about 180°C in the first 15 years of water injection. During this period, the EGS can maintain a high-efficiency heat production. After 30 years of cool water injection, the outlet temperature could remain above 140°C. Then, the cooling area gradually expands to the production well. The temperature drop along the production well is not uniform. It is faster in the zones with dense fractures where the faster flow of working fluid can more quickly transfer the thermal heat of the reservoir.Under THM coupling, the aperture of fractures is in the range of 60μm to 90 μm. Thus, the permeability in fractures is several orders of magnitude higher than that in the matrix. As shown in Figure 9(b), the outlet flow velocity in fractures (scatter markers) is about 7 orders of magnitude higher than that in the matrix (line markers). At the initial stage of heat extraction (before 1 year), the flow velocity in the matrix away from fractures is greater than that near fractures. The water in the matrix near fractures firstly flows into the fractures and continues to flow into the production well.


## 5.4. Comparison of Different Coupling Models

The TH coupling is the core for heat extraction. In the EGS operation, solid deformation and some other factors affect the TH coupling, thus modifying the heat extraction efficiency. In this section, four scenarios in the fully coupling THM model are investigated to comparatively analyze the different coupling effects on the outlet temperature and flow velocity. Table3 lists the coupling effect in each scenario.Table 3
Different coupling factors in four scenarios.
Coupling factorFormula termScenario 1 (TH)Scenario 2 THM)Scenario 3 (TH+roughness)Scenario 4 (THM+roughness)HMInitial stress equilibrium√√√√THρfCfus⋅∇TbρfCfuf⋅∇TTρfT, μT√√√√HM for the matrixρfα∂εv/∂t−α∇p✗√✗√HM for the fractureδ=σnδm/σn+kn0δm−α∇p✗√✗√TM−KαT∇T✗√✗√Fracture roughnessfσ=1+8.8σr/b1.5✗✗√√The initial stress equilibrium refers to the initial state of fluid-solid static coupling before cool water injection. Hence, the first step is to calculate the initial stress equilibrium. Then, the four models with different couplings are implemented to investigate their different coupling effects on heat transfer.Based on the simulation results, the temperature distributions of the four scenarios at 30 years are presented in Figure10. The temperature difference of the four scenarios mainly lies in the locations of thermal fronts. The thermal front of the THM model is more ahead than that of the other three models. This is mainly due to the increase in heat transfer efficiency caused by the enlargement of fracture aperture. Since fracture roughness results in the decrease in reservoir permeability, the thermal front of the TH+roughness model is in the most lag behind place.Figure 10
Temperature distributions of the four scenarios at the 30th year.The outlet temperature and flow velocity are plotted in Figure11. Compared to the TH model, the THM model with smooth fracture has higher flow velocity but lower temperature at the production well. This implies that solid deformation can enhance reservoir permeability and finally increases heat transfer efficiency. If the fracture roughness is considered, the outlet velocity decreases but the outlet temperature increases. Hence, the fracture roughness leads to a decrease in reservoir permeability and finally decreases the heat transfer efficiency.Figure 11
Outlet temperature and flow velocity under different scenarios.
(a)
Temperature 25th year(b)
Flow velocity 25th yearThe average outlet temperature is calculated by(26)Tout=∑ufbT+∫usTdy∑ufb+∫usdy.The summation term refers to the fractures while the integration term refers to the rock matrix. The average outlet temperature is a comprehensive indicator for the evaluation of EGS performance. The outlet temperatures at 25 years are plotted in Figure12 for the four scenarios. It is observed that the TH model with fracture roughness has the highest average outlet temperature. The THM model with smooth fracture has the lowest average outlet temperature. The variation of fracture aperture induced by matrix deformation leads to the enhancement of heat transfer efficiency. However, fracture roughness leads to a decrease in heat transfer efficiency. It is observed that the average outlet temperature of the fully coupling model is similar to that of the simple TH model in this case.Figure 12
Average outlet temperatures of the four scenarios.


## 6. Conclusions

This study developed a fully coupling thermal-hydro-mechanical model based on a discrete fracture network. The effects of deformation-induced aperture variation and fracture roughness on heat transfer efficiency in geothermal reservoirs were numerically investigated. Both TH and HM coupling effects in this coupling model were validated by either analytical solution or literature data. The different coupling effects on the distributions of pressure and temperature and heat recovery efficiency were comparatively analyzed through four scenarios. Based on these investigations, the following conclusions can be drawn.(1)
The pressure distribution and temperature advance in the geothermal reservoir are severely impacted by the discrete fracture network. The lower pressure gradient and faster temperature advance usually occur in the zone with denser interconnected fractures(2)
The deformation of the matrix induced by cooler water increases the aperture of the fracture and its permeability. This will enhance the heat transfer efficiency. However, the fracture roughness reduces the fracture flow velocity, thus reducing the heat transfer efficiency(3)
The simple TH model can predict the similar outlet temperature to the fully coupling THM model in our calculation case. This may be due to the fact that the combining coupling effects of deformation and roughness are ignorable in that case. On this sense, a TH model is a simple and practical tool for heat extraction prediction



---

*Source: 8875918-2020-11-23.xml*