Investigation on Heat Extraction Performance of Fractured Geothermal Reservoir Using Coupled Thermal-Hydraulic-Mechanical Model Based on Equivalent Continuum Method

Abstract

Natural fractures and artificial fractures in a tight rock matrix of an enhanced geothermal system make flow and heat transfer become seriously anisotropic. In this paper, a field-scale fractured heterogeneous geothermal reservoir model is built to study the heat transfer process. Based on an equivalent continuum method and local thermal non-equilibrium model, an equivalent permeability tensor is mapped from discrete fractures and a coupled thermal-hydraulic-mechanical mathematical model is adopted in which logarithmic stress sensitivity model is used to couple effective stress and permeability. From numerical simulation results at different injection rates, the contour results show the heterogeneity of flow, heat transfer and stress sensitivity are dominated by fractures distribution. Temperature contours reveal that the heat convection between water and rock in a fracture is more intense than the heat conduction between rock under different temperatures. The predicted power generation of a geothermal plant reveals the adverse effect on heat conversion efficiency, which is caused by the temperature drop at high injection rates.

1. Introduction

Geothermal is one of the promising renewable energy resources around the world. Hot dry rock (HDR) has great possibility to provide large amounts of energy due to the extensive distribution of the heat resource. The assessment report published by the MIT-led interdisciplinary panel forecast an enhanced geothermal system (EGS) could provide 100 GW electricity or more of cost-competitive generating capacity in the next 50 years in USA [1]. The technology challenge of extracting heat from tight rock buried at depths of more than 2 km is to create connected pathways by hydraulic fracturing through which fluid can be induced to flow. Although the Landau plant in Germany, the Soultz plant in France and the recently commissioned Habanero pilot plant prove the commercial feasibility of geothermal power plants [2,3], there are still some complex and unresolved issues in the commercial production of energy from EGS. One of the most critical is difficulty of predicting real-time underground flow and heat transfer condition when considering intricate fracture networks which leads to a lack of clear understanding of system behavior, output and performance.

In practice, the artificial fractures generated by hydraulic fracturing in HDR can induce large-scale flows through a tight rock matrix. Therefore, it is necessary to consider fracture networks in the process of numerical simulation. Despite recent development of fracture network modeling, the complexity of fracture geometry and the fracture connectivity poses further challenges for simulating fluid flows through the HDR. Even if for a conceptual model, when taking explicit discrete fracture networks into consideration, the complexity of simulation increases drastically.

In recent years, remarkable advances have been made in the development of thermal-hydraulic-mechanical (THM) coupling models. Xu et al. [4] proposed a simplified approach to simulate the coupled hydro-thermal system for EGS, capable of providing a detailed prediction of fluid flow and heat transfer in a geothermal reservoir based on an equivalent pipe network model. Jiang et al. [5] presented a three-dimensional transient model for EGS subsurface thermo-hydraulic process by which the geothermal reservoir is treated as an equivalent porous medium of a single porosity. Shaik et al. [6] developed a numerical program to simulate the heat extraction from naturally fractured geothermal systems by coupling fluid flow with heat transfer between the rock matrix and circulating fluid. This provides a dynamic treatment of the characteristic properties (aperture, length and orientation) of individual fractures. Bahrami et al. [7] studied several self-propped single fracture THM models based on the poro-elastic theory. Zhao et al. [8] established a 3D THM coupling model of fractured media to simulate the extraction of HDR geothermal energy, by which the fracture area was discretely modeled by the Goodman joint element. Guo et al. [9] investigated the effects of aperture heterogeneity on flow pattern evolution in a single fracture in a low-permeability crystalline formation. Sun et al. [10] simulated with 2D stochastically generated fracture models to study the characteristics of fluid flow, heat transfer and mechanical response in geothermal reservoirs. Most of the above research was fundamental and based on simplified conceptual models, and the lack of consideration for the complexity of geological structures and fracture systems, results in limitations upon application to field-scale issues.

To utilize a THM coupling model in field-scale cases, a numerical method which fully integrates the anisotropy of fractures and heat transfer between water and rock should be developed to investigate the actual performance of EGS. In this paper, a field-scale fractured geothermal reservoir model is built according to an actual geologic model. To increase calculation efficiency, we adopt the equivalent continuum method to simulate flow, heat transfer and mechanical response, in a which logarithmic stress sensitivity model is employed to couple effective stress and permeability, and the equivalent permeability tensor is mapped from discrete fractures. Considering the local thermal non-equilibrium theory, two energy equations, which can describe heat transfer between rock fabric and pore-fluid, are defined respectively. We focus on numerical simulation results at different injection rates and analyze the temperature evolution. We also predict power generation of a geothermal plant according to the simulation result. This model can better reflect the anisotropy of a fractured reservoir and the simulation results can provide a better representation of the real operational status of a geothermal reservoir compared with simplified conceptual model.

2. Methodology

2.1. Simulation Methods of Flow in Fractured Rock

Simulation methods can be broadly classified into two kinds of approaches by the method of dealing with fractures, namely, discrete fracture network (DFN) and equivalent continuum model (ECM). In the DFN approach, the discrete medium model makes an explicit description of fractures in the rock matrix. Two different methods, which include different governing equations, initial and boundary condition, are employed to simulate flow, heat transfer and mechanical response in the fracture and rock matrix respectively. Sun et al. [10] employed the DFN approach to simulate a stochastically generated 2-D fracture network model with a full THM coupling method. The DFN method is suitable for fundamental research because the parameters of fractures, such as direction and width, are defined accurately during the process of numerical simulation. However, some difficulties on complex fracture network modeling and high computational overhead limit its application to field-scale cases.

To reduce computational overheads and increase availability for field-scale cases, we adopt the equivalent continuum method in this work. Being different from the DFN approach, the equivalent medium model regards the overall fractured reservoir as a continuous porous system. We can represent its heterogeneity by the corresponding equivalent parameters. This model has a high calculation efficiency and a simple requirement for parameters. It has gained a far-reaching development in rock hydraulics. The permeability tensor represents the anisotropic flow resistance caused by the rock fabric, the resistance is added to the momentum equation as a pressure drop (Darcy term in Equation (1)) when fluid flows in porous media. The logarithmic stress sensitivity model describes the coupled relationship between effective stress and permeability [11]. The finite volume method is adopted to discretize model because of its adaptation to the geothermal reservoir boundaries and advantages of high computing efficiency for field-scale models.

2.2. Basic Assumptions

In this work, the following assumptions are made to formulate the mathematical model:

• The fluid flow in heat reservoir is in laminar regime and the Reynolds number is sufficiently low to meet the conditions of the Darcy law. Ignoring inertial loss, the pressure drop caused by porous resistance is proportional to velocity. This greatly simplifies the complexity of pores and throats in rock.
• The heat transmission fluid (i.e., water) always remains liquid. The temperature is always lower than 500 K, yet below the supercritical point of water (647 K, 22.064 MPa) [12], therefore it is reasonable to assume the water is always in liquid state.
• The porous heat reservoir is always fluid saturated, no immiscible gas is trapped in the rock and only single-phase flow is considered.
• There do not exist any fluid-structure interactions, including chemical dissolution, deposition, physical pressing, etc.
• The heat transmission fluid (i.e., water) is incompressible. The rock in this study is considered as thermo-poroelastic systems [13] based on a small strain assumption and obeys the generalized Hooke’s law.

2.3. The Governing Equations

With the aforementioned assumptions, the governing equation in mathematical model is focused on describing the subsurface thermal-hydraulic-mechanical process in EGS, three kinds of governing equations in THM coupling mode can be expressed as follows:

2.3.1. Basic Flow Governing Equations

In the equivalent continuum method, the superficial velocity porous formulation generally gives good representations of the pressure loss through a porous region, and superficial velocity is calculated based on volumetric flow rate. Darcy’s law is regarded as a source term added to Navier–Stokes equations [14].

a. Continuity equation

$$\frac{\partial \left(\epsilon \rho \right)}{\partial t}+\nabla ·\left(\rho \mathit{u}\right)=0$$

(1)

b. Momentum equation

$$\frac{\partial \left(\rho \mathit{u}\right)}{\partial t}+\rho \mathit{u}\nabla \mathit{u}-\mu {\nabla }^2\mathit{u}=-\nabla P+\frac{\mu }{\mathit{k}}\mathit{u}+\rho g$$

(2)

where $\epsilon$ is porosity, $\mathit{u}$ is superficial velocity and equal to $\epsilon {\mathit{u}}_{\mathit{p}\mathit{h}\mathit{y}\mathit{s}\mathit{i}\mathit{c}\mathit{a}\mathit{l}}$, $\mathit{k}$ is permeability tensor, Equation (2) is the standard format Navier–Stokes equation, in which it goes in order from left to right: unsteady term, convection term, diffusion term, source terms respectively.

2.3.2. Energy Equations

Based on local thermal non-equilibrium theory, two energy equations, which, respectively, describe heat transport in rock fabric and pore-fluid, are defined as follows [14]:

a. Fluid energy equation

$$\frac{\partial \left[\epsilon \rho C_f{T}_f\right]}{\partial t}+\mathit{u}·\nabla \left[\rho C_f{T}_f\right]=\nabla ·\left(\epsilon k_f\nabla T_f\right)+\mathrm{ha}\left(T_r-T_f\right)$$

(3)

b. Rock energy equation

$$\frac{\partial \left[\left(1-\epsilon \right)\rho C_r{T}_r\right]}{\partial t}=\nabla ·\left(\left(1-\epsilon \right)k_s\nabla T_r\right)-\mathrm{ha}\left(T_r-T_f\right)$$

(4)

where, $h$ is convective heat transfer coefficient and equal to $\frac{k_{f}Nu}{2d}$, $\mathrm{a}$ is specific surface area and equal to $\frac{2\epsilon }{d}$, $k_f$ and $k_s$ are thermal conductivity of fluid and rock respectively.

2.3.3. Geo-Mechanical Equations

Under the assumption of linear elasticity with small strains for a thermo-poroelastic system, the stress equilibrium equation can be expressed as follows [13]:

$$\frac{3\left(1-v\right)}{\left(1+v\right)}{\nabla }^2\sigma =\frac{2\alpha \left(1-2v\right)}{\left(1+v\right)}{\nabla }^2\mathrm{P}+\frac{6\beta B\left(1-2v\right)}{\left(1+v\right)}{\nabla }^2{T}_s-\nabla \mathit{F}$$

(5)

$$B=\frac{3\left(1-2v\right)}{E}$$

(6)

where $\sigma$ is mean total stress of rock, $v$ is rock Poisson ratio, $\alpha$ is Biot coefficient, $\beta$ is rock linear thermal expansion coefficient, $\mathit{F}$ is the external body force, $B$ is bulk modulus and Equation (6) is the relationship between bulk modulus and Young’s modulus. Based on the classical consolidation theory developed by BIOT, Equation (5) indicates that rock stress changes with pore pressure and rock temperature.

2.4. Coupling Method

Effective stress is represented the force that keeps a collection of rock fabric. The relationship between effective stress, the mean total stress and pore pressure is shown as Equation (7).

(1) Effective stress equation

$${\sigma }_{eff}=\sigma -\alpha p$$

(7)

(2) Logarithmic stress sensitivity model

$$\frac{K}{K_0}=1-S\ln \left(\frac{{\sigma }_{eff}}{{\sigma }_0}\right)$$

(8)

where, ${\sigma }_{eff}$ is effective stress, $K_0$ and ${\sigma }_0$ is primary permeability and effective stress, S is stress sensitivity coefficient.

The fractured granite is categorized as strong stress sensitivity rock [13]. The logarithmic stress sensitivity model (Equation (8)) is suitable to calculate the permeability with the change of effective stress.

2.5. Solving Method

The aforementioned governing equations are solved in the commercial CFD solver, FluentV6. Continuity equation (Equation (1)) and Momentum equation (Equation (2)) are solved by AAMG (aggregative Algebraic Multigrid) solver. Two energy equations (Equations (3) and (4)) and geo-mechanical equation (Equation (5)) are specified as three UDS (User Defined Scalars) equations respectively and solved by SAMG (selective algebraic multigrid) solver. Dual mesh approach is applied for local thermal non-equilibrium model. Two separate but overlapped meshes for fluid and solid zone are created, which represent fluid zone and solid zone respectively. Basic flow equations and the fluid energy equation are discretized on the fluid zone. The rock energy equation and geo-mechanical equation are discretized on the solid zone. The logarithmic stress sensitivity model is coupled with flow, energy and geo-mechanical equations indirectly through a UDF (user defined function). The mature SIMPLE (semi-implicit method for pressure linked equation) algorithm is employed to address the pressure-velocity coupling. The first order upwind differencing scheme is generally adopted to discretize the spatial-derivative terms and a fully implicit scheme is adopted for the discretization of the transient terms.

3. Model Analysis

3.1. Background of Field-Scale Model

Bongor basin has the area of $1.6\times {10}^4 {\mathrm{km}}^2$ located in the southwest of Chad, a central African country [15]. The basin extends from east to west for about 280 km in length and 40–80 km wide [16]. Abundant HDR is discovered in the granite basement in early 2013. Simulation case in this paper was taken from Baobab buried-hill (Figure 1) which mainly consists of Precambrian Monzonite, Amphibolite and Gneisses.

3.2. Boundary Conditions and Geometrical Parameters

As Figure 2 shows, the model geometry is built based on the reservoir geology, it is 2730m long, 1660 m wide, 290 thickness and with a total volume of $6.3\times {10}^8{\mathrm{m}}^3$ (

Table 1. Model geometric dimension and geophysics parameters.

Parameter	Value
Injection well depth, ${\mathrm{h}}_{\mathrm{in}}$($\mathrm{m}$)	3460
Production well depth, ${\mathrm{h}}_{\mathrm{pro}}$($\mathrm{m}$)	3590
Geothermal reservoir length, ${\mathrm{L}}_{\mathrm{l}}\left(\mathrm{m}\right)$	2734
Geothermal reservoir width, ${\mathrm{L}}_{\mathrm{w}}\left(\mathrm{m}\right)$	1660
Geothermal reservoir thickness, ${\mathrm{L}}_{\mathrm{t}}\left(\mathrm{m}\right)$	290
Interwell separation, ${\mathrm{d}}_{\mathrm{ws}}\left(\mathrm{m}\right)$	980
Borehole size, ${\mathrm{d}}_b\left(\mathrm{m}\right)$	0.3
Average temperature, ${\mathrm{T}}_{\mathrm{ave}}\left(\mathrm{K}\right)$	$454$
Initial average permeability, $k_{\mathrm{ave}}\left({\mathrm{m}}^2\right)$	$5.32\times {10}^{-17}$
Initial average porosity, ${\epsilon }_{\mathrm{ave}}(\%)$	$0.23$
Reservoir total heat, ${\mathrm{h}}_t\left(\mathrm{J}\right)$	$1.52\times {10}^{18}$

).

The well layout depends on geological conditions and physical parameters of the reservoir. Doublet well layout is employed because of the extremely low permeability of the rock matrix and the heterogeneous distribution of fractures. A production well is deployed 980 m away from the injection well according to the fractures’ distribution. The injection and production wells both penetrate through the reservoir and boreholes are connected to the reservoir in the thickness range of 290 m. Reservoir boundary is regarded as thermal compensation wall and the change of surface temperature is calculated by analytic approach specified as a UDF.

3.3. Physics Parameters and Properties

3.3.1. Initial Distribution of Permeability

Permeability and porosity are crucial factors influencing flow and heat transfer. These are built from well logging data and shown in Figure 3. In order to keep the isotropic initial permeability of matrix and equivalent permeability tensor of fractures in a consistent form, the former is transformed to a spherical tensor. The reservoir average initial permeability and porosity are $5.32\times {10}^{-17} {\mathrm{m}}^2$ (0.0532 md) and 0.23% respectively. Because the reservoir initial rock is so tight that the heat transmission fluid is barely able to flow, hydraulic stimulation is required to induce large-scale flow.

3.3.2. Physics Parameters of Rock and Fluid

The liquid water is regarded as incompressible fluid in this case, and detailed parameters are listed in

Table 2. Parameters of rock and fluid.

Parameter	Value
Rock thermal conductivity, ${\mathsf{\lambda }}_r$($\mathrm{W}/\mathrm{m}/\mathrm{k}$)	3.42
Rock density, ${\mathsf{\rho }}_r$($\mathrm{kg}/{\mathrm{m}}^3$)	2700
Rock specific heat, $c_r\left(\mathrm{J}/\mathrm{kg}/\mathrm{K}\right)$	1125
Elastic modulus, $\mathrm{E}\left(\mathrm{GPa}\right)$	41
Poisson ratio, $\mathsf{\nu }$	0.23
Biot coefficient, $\mathsf{\alpha }$	1
Rock linear thermal expansion, $\mathsf{\beta }\left(1/\mathrm{K}\right)$	$5\times {10}^{-6}$
Fluid thermal conductivity, ${\mathsf{\lambda }}_w$($\mathrm{W}/\mathrm{m}/\mathrm{k}$)	0.6
Fluid density, ${\mathsf{\rho }}_w$($\mathrm{kg}/{\mathrm{m}}^3$)	1000
Fluid specific heat, $c_w\left(\mathrm{J}/\mathrm{kg}/\mathrm{K}\right)$	4200

. Being different from water, rock fabric is assumed as the thermo-poroelastic media. Properties of rock fabric change with rock temperature and pore pressure. This assumes that rock fabric can move as an elastic material and obeys the generalized version of Hooke’s law. Under the assumption of linear elasticity with small strains for a thermo-poroelastic system, the equations can be expressed as Equation (5) and the parameters associated with the Geo-mechanical equations (Equation (5)) are listed in the Table 2.

3.3.3. Initial Temperature Distribution

The geothermal gradient is the crucial parameter to evaluate reservoir commercial development potential. Higher geothermal gradient not only means higher temperature and more extractable heat, but also higher electricity conversion efficiency of geothermal power plants [18]. Figure 4 shows borehole temperature distribution taken from well logging data. Although the average geothermal gradient is $41.6°\mathrm{C}/\mathrm{km}$, the basement geothermal gradient reach $79.1°\mathrm{C}/\mathrm{km}$ at a depth of 3100 m. HDR is between depths of 3170.1–3460.6 m and with a range of borehole temperature from $165.1°\mathrm{C}(438.1\mathrm{K}$) to $192.5°\mathrm{C}\left(465.5\mathrm{K}\right)$.

The geothermal gradient is a considerable factor in the process of reservoir modeling because the thickness reaches up to 290 m. The vertical temperature distribution (Figure 5) in the model is taken from the previous well logging data shown in Figure 3. Total heat in the reservoir is determined by the temperature, the specific heat of rock and total volume. Average temperature and total heat are $181°\mathrm{C}$ ($454\mathrm{K}$) and $1.52\times {10}^{18}\mathrm{J}$, respectively (Table 2).

3.4. Equivalent Continuum Method for Discrete Fractures

The hydraulic stimulation is the essential measure before extracting heat from HDR. Fractures in the rock create connected pathways through which fluid can be induced to flow and extract heat. In general, most hydraulic fracture apertures are in the range of 0.1 mm to 1.0 mm [9]. At present, there are still some difficulties in fracture modeling and simulating the flow of complex fractures. Besides following the minimum principal stress criterion, there is a certain randomness in fractures extension. In addition, there is very large computational overhead to simulate a 3D field-scale case with explicit fractures.

3.4.1. Equivalent Continuum Method

In view of above difficulties, we try to characterize fractures by the equivalent continuum method. Based on macro scale continuum theory, the fractures in a rock matrix can be equivalent to a permeability tensor to represent the anisotropic flow in fractured porous media (Figure 6). The permeability tensor is calculated based on the fracture attributes (e.g., the fracture size, aperture and orientation in Figure 6a) [19,20] as:

$${\mathit{P}}_{\mathit{i}\mathbf{,}\mathit{j}}={\rho }_{fra}\frac{\pi }{4}{{\displaystyle \int }}_0^{\mathit{f}\mathbf{\left(}\mathit{t}\mathbf{\right)}}{{\displaystyle \int }}_0^{\mathit{f}\mathbf{\left(}\mathit{r}\mathbf{\right)}}{{\displaystyle \int }}_{\mathit{\Omega }}{r}^2{b}^3{\mathit{m}}_{\mathit{i}\mathbf{,}\mathit{j}}E\left(n,b,t\right)\mathrm{d}\mathit{\Omega }\mathrm{d}r\mathrm{d}b$$

(9)

where, ${\mathit{P}}_{\mathit{i}\mathbf{,}\mathit{j}}$ is the fracture tensor that considers the position, shape, aperture and orientation of fractures in a representative elementary volume and averages these features in each arbitrary direction; ${\rho }_{fra}$ is the density of fracture planes in the control volume; $E\left(n,b,t\right)$ is a probability density function that describes the number of fractures with size; ${\mathit{m}}_{\mathit{i}\mathbf{,}\mathit{j}}$ is the unit vector to the fracture plane oriented in a small solid angle $\mathrm{d}\mathit{\Omega }$. This concept has been extended [21] to represent a permeability tensor based on the assumption that the fluid is channeled in parallel fracture planes with volumetric flow rate proportional to $b^3$. Thus, the permeability tensor ${\mathit{k}}_{\mathit{i}\mathbf{,}\mathit{j}}$ is represented as:

$${\mathit{k}}_{\mathit{i}\mathbf{,}\mathit{j}}=\lambda \left(P_{kk}{\delta }_{i,j}-{\mathit{P}}_{\mathit{i}\mathbf{,}\mathit{j}}\right)$$

(10)

where, λ is a dimensionless constant associated with fractures interconnectivity and is restricted between 0 and 1/12; ${\delta }_{i,j}$ is the Kronecker delta; $\mathit{i}$ and $\mathit{j}$ represent Cartesian coordinate directions x, y, z. For a 2D model, $\mathit{i}$ and $\mathit{j}$ are defined as x and y (Figure 6b).

3.4.2. Mapping Permeability Tensor from Discrete Fractures

Stochastic fracture networks are modeled based on the aforementioned reservoir geology condition to represent the fractured geothermal reservoir (Figure 7) and are used to predict the performance of EGS heat extraction. Two sampling lines are established to extract internal variables, the sample line 1 is along the two wells and the sample lines 2 is vertical to the sample line 1. The aperture of 14 fractures is in the range of 0.5 mm to 1 mm.

The permeability tensor is calculated as a second-order symmetric tensor by the equivalent continuum method (Equations (9) and (10)). In order to show the size of permeability tensor of control volume in a 3D contours plot, $\overline{k}$ is an invariant scale of permeability tensor and is defined as:

$$\overline{k}=\frac{\mathit{t}\mathit{r}\mathbf{\left(}{\mathit{k}}_{\mathit{i}\mathbf{,}\mathit{j}}\mathbf{\right)}}{\mathbf{3}}$$

(11)

where, $\mathit{t}\mathit{r}\mathbf{\left(}{\mathit{k}}_{\mathit{i}\mathbf{,}\mathit{j}}\mathbf{\right)}$ is the trace of tensor and equal to the sum of three diagonal elements.

Figure 8 shows the $\overline{k}$ of the control volume permeability tensor after calculation. Fourteen fractures are discretized in the reservoir. $\overline{k}$ of 14 fractures is between $2.896\times {10}^{-13}{\mathrm{m}}^2 \left(289\mathrm{md}\right)$ and $1.061\times {10}^{-12}{\mathrm{m}}^2\left(1061\mathrm{md}\right)$.

5. Results and Discussion

5.1. Analysis of Hydraulic Process

A clear understanding of flow status is important to monitor EGS performance. The water flows in fractures, carrying heat from rock to the production borehole. Therefore, the flow velocity directly determines the intensity of heat transfer. Figure 9 shows the velocity (averaged superficial velocity in three directions) distribution and streamlines at the injection rate of 81.9 kg/s when the flow evolves semi-steady state. There is an obvious difference of the velocity distribution between fractures and rock matrix. The velocity decreases rapidly around the injection well. The velocity contours shows that the water barely flows in the rock matrix in which the velocity falls to $1.12\times {10}^{-11} \mathrm{m}/\mathrm{s}$. Streamlines indicate that fractures are preferential pathways of water flow. Fractures are proved to be the dominant regions of subsurface flow in Figure 9; the distribution of fractures has a principal impact on the EGS operation.

5.2. Analysis of Heat Transfer in Different Conditions

In order to study on the heat transfer under different conditions, three different injection rates ($q_{in}$), 81.9 kg/s, 136.6 kg/s, 191.2 kg/s, are chosen as borehole boundary conditions, respectively.

Figure 10, Figure 11 and Figure 12 show the evolution of temperature in seven profiles at three different flow rates. At the initial running stage, the cold-water flows along fractures while there is a rapid drop in rock temperature around borehole. The water flows into the reservoir, transferring heat from rock until the water temperature rises to the same level as the rock. The production water, therefore, maintains a constant high temperature at initial running stage. With heat being extracted from rock, a low temperature zone develops and expands towards the production well. As fractures are dominant flow regions, heat convection between the water and the rock mainly occurs in these regions. In the tight original rock, conduction is the only heat transfer mode. The temperature in original rock region only drops slowly on account of the low thermal conductivity of the rock.

Two sample lines are set in the model (Figure 7) to extract the internal temperature data. Figure 13 and Figure 14 show the temperature distribution of the 136.6 kg/s case at different times, the intersection postion of fractures and sample line is marked in figures. The temperature drops rapidly in fracture locations and relatively slowly in the rock matrix.

The strong heterogeneous temperature distribution is depicted in the contours at different injection rates because of the heterogeneity of fracture distributions and the anisotropic flow conductivity of fractures. Convection and conduction are two basic modes of the heat transfer process in the EGS operation. The heat convection between the water and the rock is the dominant mode, which leads to a rapid temperature drop in the fractures. According to the convective heat transfer theory, the intensity of convection mainly depends on the flow velocity in porous media [22]. Conduction is the subordinate mode and it conducts heat from the high-temperature original rock to the low-temperature fractures. For the case with lower injection rate, the temperature in fractures drops more slowly because of a lower proportion between convection intensity and conduction intensity.

5.3. Analysis of Stress Sensitivity

The logarithmic stress sensitivity model (Equation (8)) calculates the variation of permeability with the change of effective stress. Figure 15 shows the evolution of permeability at two separate times at the injection rate of 136.6 kg/s. The rising of pore pressure and the dropping of rock temperature both have an impact on the permeability. The rising of pore pressure caused by water injection reduces the rock effective stress, which has a positive effect on the increase of permeability. In addition, the thermal contraction of rock reduces the mean total stress due to the temperature dropping, which further decrease the effective stress. The permeability around the injection well and in fractures increases due to the combined effect of the above two factors.

5.4. Analysis of Enhanced Geothermal System (EGS) Operation Prediction

For a quantitative description of the EGS performance, three parameters relevant to the EGS performance, the average production temperature ($T_{pro}$), the average reservoir temperature ($T_{ave}$) and the heat extraction ratio (${{\rm Y}}_{ra}$), are defined through Equations (9), (10), and (11) respectively. ${{\rm Y}}_{ra}$ is the proportion of total extracted heat to the total heat in the reservoir.

$$T_{pro}\left(t\right)=\frac{\int T_{well}\left(t\right)dl}{l}$$

(12)

$$T_{ave}\left(t\right)=\frac{\int T_{cell}\left(t\right)dv}{V}$$

(13)

$${{\rm Y}}_{ra}\left(t\right)=\frac{\int c_r\Delta T_{cell}\left(t\right)dv}{\int c_r\Delta T_{A}dv}$$

(14)

where $l$ is the wellbore length in the depth range of the reservoir, $T_{well}$ and $T_{cell}$ are the instant temperature in the wall of the borehole and the temperature in the reservoir cell. $\Delta T_{cell}$ is the temperature difference between an instant time (t) and the initial time, $\Delta T_A$ is the difference between the initial temperature and the injection temperature.

Figure 16, Figure 17 and Figure 18 show the changes of $T_{pro}$, $T_{ave}$ and ${{\rm Y}}_{ra}$ with time at different injection rates; three parameters change more rapidly at the higher injection rate. Comparing the changes of the average production temperature ($T_{pro}$) at three different injection rates, the drop of lower injection rate case is relatively slow, and is caused by the lower proportion between heat convection and heat conduction.

As shown in Figure 16, a reference production temperature (130 °C) is chosen to study on the EGS performance. In addition, the $T_{ave}$ and the ${{\rm Y}}_{ra}$ at the reference temperature are marked in Figure 17 and Figure 18. The total extracted heat can be calculated by multiplying the heat extraction ratio (${{\rm Y}}_{ra}$) and the reservoir total heat ($h_t)$ in Table 1. When the production temperature drops to 130 °C, the case with lower injection rate (81.9 kg/s) extracts 20% more heat than the case with higher injection rate (191.2 kg/s), it is because more heat transfer occurs original high-temperature rock over a long duration of operation.

5.5. Analysis of Power Generation Prediction

Power generation is the major way in the multistage utilization process of geothermal. Generally, the conversion efficiency of geothermal power developments is lower than that of all conventional thermal power plants due to the production temperatures. The conversion efficiency for each geothermal power plant is varies with the power plant design (single or double flash, triple flash, dry steam, binary, or hybrid system) and is mainly influenced by non-condensable gas content, parasitic load, heat loss, turbine efficiency, and generator efficiency [23]. The highest reported conversion efficiency is approximately 21%, with an average of about 12% based on the data of 94 worldwide geothermal power plants [24]. Zarrouk and Moon et al. [25] produce a generic model for the conversion efficiency as a function of enthalpy (Equation (12)) though fitting all the available data from 94 geothermal power plants with one curve. Equations (13) and (14) calculate the instant output power and the total electricity generation of the geothermal plant, respectively.

$${\eta }_{act}=7.8795\ln (h_{in})-45.651$$

(15)

$$P_e={\eta }_{act}q{c}_w\left(T_{in}-T_n\right)$$

(16)

$$W_c\left(t\right)=\int P_{e}dt$$

(17)

where $h_{in}$ is the enthalpy of water entering power plant and equal to the enthalpy of production water, the same is true with $T_{in}$, $T_n$ is the water temperature flowing from the plant, $q$ is the water flux, $c_w$ is the specific heat of water listed in Table 2.

As shown in Figure 18, a higher injection rate means a higher power output at initial operation stage. Nevertheless, temperature dropping and conversion efficiency declining of the higher injection rate case reduce the power output below other cases. Comparing these three cases, power output of the lower injection rate stabilizes in a longer period and declines slowly. The influence of different conversion efficiency can be observed in Figure 19; total electricity generation of lower injection rate case (136.6 kg/s) exceeds the case with a higher injection rate (191.2 kg/s) in the late stage of EGS operation, even if the heat extraction ratio of that case is superior in Figure 20.

6. Conclusions

In this paper, a geothermal reservoir model with 14 discrete fractures has been built based on a field-scale geological model. The equivalent continuum method is adapted to represent the fractures and map the permeability tensor. The coupled thermal-hydraulic-mechanical process is simulated. We analyzed numerical simulation results at different injection rates.

• The equivalent continuum method is adopted in order to reduce computational overhead. The permeability tensor calculated from discrete fractures is employed to characterize the fractures. The thermal-hydraulic-mechanical process is coupled by a logarithmic stress sensitivity model.
• Strong heterogeneity is shown in the THM process because of the heterogenous distribution of fractures. In the heat transfer process, fractures are the dominant regions in which water transfers heat from the rock matrix by convection, in addition, the fractures have an impact on the pore pressure distribution, which increase the heterogeneity of geo-mechanical process.
• The EGS operation prediction data has a quantitative description of the subsurface evolving with time in different conditions. The case with high injection rates can rapidly extract heat, and has a relatively lower efficiency and unstable production temperature, however, occur in the late stage of operation. Heat conduction is not able to sufficiently transfer heat from the tight original rock to the fractures. This illustrates the necessity of well pattern optimization and intensive stimulation in the late stage of EGS operation.
• Most of the parameters used in the simulation are taken from actual field data, which includes geometrical data, geothermal gradient data, rock mechanical properties. Therefore, this model can better reflect actual geological physical properties and the simulation results are closer to actual reservoir operation.

Although there may be differences between simulation results and the actual operation of an EGS due to the assumptions and uncertainty of the geological representation, the approach adopted here offers advantages over the simplified conceptual model and DFN method. This method better reflects the complexity of the geology and provides a fast approach to simulating fractured geothermal reservoir performance. These simulation results should offer a closer representation of real reservoir operation and should prove helpful in the development of EGS systems.