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

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 THMmodel 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.


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][3][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 geo-thermal 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][12][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][16][17][18]. Recent THM coupling models for geothermal reservoirs consider the impact of several isolated discrete fractures on heat extraction performance [19][20][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. Com-pared 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.

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 where l min and l max are the minimum and maximum trace lengths, D is the fractal dimension of length distribution of fractures, and F f is a random number uniformly distributed in the range of 0 to 1.

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 where γ is a fracture orientation angle in a set of fractures and F 1 and F 2 are random numbers uniformly distributed in the range of 0 to 1.

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 of 300 m × 300 m), the fracture center will be regenerated. The coordinate generation function for the midpoints of fractures is expressed as 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 Figure 1. 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.

Coupling Model for Energy Extraction from a Geothermal Reservoir
where ρ f is the fluid density, ϕ is the porosity of the matrix, t denotes time, u s is the flow velocity of the Darcy seepage, and Q s is the source or sink term. The seepage velocity u s is expressed by Darcy's law as Furthermore, the storage term at the first term of Equa-tion (4) is The porosity dilation-induced fluid variation is treated as the source term Q s which is shown as 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 where b is the fracture aperture, u f is the flow velocity in the fracture, ∇ T denotes the gradient operator restricted to the tangential plane of a fracture, and Q f is the source or sink term. A fracture consists of two smooth parallel plates; thus, the volumetric flow rate Q is and the permeability of a single fracture is then obtained as where w is the fracture width and b is the fracture aperture. ∇ T p is the pressure gradient along the fracture. We introduce f σ of the roughness factor to consider the effect of fracture roughness on flow: The following empirical expression is used in our computation: 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 3 Geofluids from 1.04 to 1.65 [33]. Hence, the fracture flow velocity u f is 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 where n u and n d are the normal direction of up and down surfaces of the fracture, respectively. Finally, the governing equation for fracture flow is obtained as 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 where ðρCÞ eff = ð1 − ϕÞρ s C s + ϕρ f C f is the effective heat capacity, ϕ is the porosity of the rock matrix, ρ s is the matrix density, C s is the specific heat capacity of the matrix, C f is the specific heat capacity of fluid, T is the matrix temperature, and q is the heat flux in the rock matrix.
The heat flux q can be expressed through Fourier's law as 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 where q f is the heat flux in the fracture and is Q Tf 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, Q Tf can be expressed as 3.4. Mechanical Equilibrium. The constitutive equation of poroelasticity for THM coupling is The Navier equation for displacement v can be expressed as 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 = 2Gð1 + νÞ/3ð1 − 2νÞ = E/3ð1 − 2νÞ is the bulk where δ = b 0 − b is the fracture closure, where b 0 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 k n0 is an empirical parameter representing the initial normal stiffness. δ m and k n0 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: σ nc ðMPaÞ = 10k n0 δ m = 487 b 0 ðmmÞ + 2:51, δ m = 0:9 b 0 .

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] 4. Simulation Implementation and Model Validation

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 Figure 2 is used to model a geothermal reservoir. This geothermal reservoir is at 3500 m deep and has a computational domain of 300 m × 300 m. 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.
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. Table 2 lists the model parameters used in simulations.

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. Figure 3 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   6 Geofluids 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.
The distribution of fracture temperature and the observation line at the 100th day are plotted in Figure 4(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.

Geofluids
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.

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 Figure 7, 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).

Reservoir Temperature Distribution.
The temperature distributions of the four fracture networks at the 30th year are presented in Figure 8. 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. Frac-ture 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.
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 Figure 8). 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.

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 Figure 9(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).
The 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.  Table 3 lists the coupling effect in each scenario. 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 Figure 10. 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.
The outlet temperature and flow velocity are plotted in Figure 11. 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.
The average outlet temperature is calculated by 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   11 Geofluids are plotted in Figure 12 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. TH HM for the fracture

Conclusions
This study developed a fully coupling thermal-hydromechanical 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  Figure 3. All the properties of rock and water are kept as constants, and the heat conduction of water is omitted. For the convenience of symmetrical derivation, the fracture aperture b is equivalent to 2b s . Based on the cubic law of fracture flow, the flow velocity in the fracture is Since all the parameters in Equation (A.1) are kept constant, the flow velocity is constant in this context. The heat transfer in the fracture can be described as