Investigation on heat extraction characteristics in randomly fractured geothermal reservoirs considering thermo‐poroelastic effects

A fully coupled thermal‐hydraulic‐mechanical (THM) model was developed to investigate the underlying response mechanisms during heat extraction in fractured geothermal reservoirs. The random fracture network in the stimulated zone was reproduced based on the fractal theory. The coupled model accounts for the dominant physical phenomena including (a) fluid flow, heat transport, and solid deformation in porous media and fractures; (b) local thermal nonequilibrium (LTNE) between rock matrix and flowing fluid; and (c) temperature‐dependent fluid thermodynamic properties and stress‐dependent pore and fracture permeability. The proposed model was validated by several analytic solutions. Sequentially, the evolution of pore pressure, equivalent temperature, effective stress, and reservoir permeability was analyzed. The sensitivity of the heat extraction performance to fracture network morphology was discussed. Results show that interconnected large‐scale fractures dominate mass and heat transport. Widely distributed small‐scale fractures contribute to the cooling of the heat rock mass along many flow paths in parallel. The change in effective stress associated with fully coupled thermo‐poroelastic effects may induce fracture shear dilation and pore expansion, resulting in an enhancement of the overall reservoir permeability. There is a significant temperature difference between the solid phase and the fluid phase in fractures, but not in porous media. It is more reasonable to use the LTNE theory to analyze the heat extraction of fractured geothermal reservoirs. The fracture network morphology has a profound effect on heat extraction efficiency and injectivity. Undeveloped fracture networks will result in a higher injection pressure, earlier thermal breakthrough, and shorter lifetime, but higher heat extraction ratio. Properly increasing the fracture density can delay the thermal breakthrough time, prolong the service life, and improve the injectivity. Fully connected fracture networks may result in thermal short‐circuiting and earlier thermal breakthrough. Generating a complex and scattered fracture network but without preferential channels is conducive to extracting more heat from geothermal reservoirs.


| INTRODUCTION
Due to the renewability, cleanness, and universality, geothermal resources have attracted broad attention around the world. The hot dry rock (HDR) geothermal resource accounts for nearly 90% of the total geothermal resources, so it has great development potential and prospects. 1 Generally, the HDR reservoir possesses extralow porosity and ultralow permeability, making it difficult to achieve economical thermal efficiency without stimulation treatment. 2 Consequently, permeability enhancement techniques play a critical role in the development of these nearly impervious reservoirs. 3 Enhanced geothermal systems (EGS) are widely regarded as an efficient method to extract heat from the HDR. 2,4 The key is to create a complex fracture network for fluid circulation between the injection well and the production well by a series of stimulation techniques, such as hydraulic fracturing, waterless fracturing, and chemical treatment. Recently, cryogenic fracturing using liquid nitrogen or liquid carbon dioxide has raised many researchers' interest due to the thermal shock effect. [5][6][7] Some mining practices have confirmed that the combination of thermal-chemical-hydraulic (TCH) fracturing technology can induce a more pronounced stimulation effect in geothermal reservoirs than using any single stimulation method. [8][9][10] Generally, EGS design is more focused on shear stimulation. [11][12][13] It is believed that fluid injection can induce fracture shear slip, which contributes to the permeability enhancement of the whole reservoir. Micro-seismic monitoring events have confirmed that complex fracture networks can be created with appropriate stimulation in HDR reservoirs. 9,10,12,14 The artificial EGS with embedded fracture networks can obviously improve the contact area between the rock matrix and the circulating fluid, thereby enhancing the heat extraction efficiency. Revealing the mass and energy transport processes in an EGS with different fracture distributions is essential to evaluating the production efficiency and lifetime of a geothermal reservoir. In addition, some pre-existing discontinuities, such as faults, beddings, and natural fractures, may be reactivated during fluid circulation, which will potentially induce formation instability, surface subsidence, or even earthquakes. Recently, injection-induced earthquakes have become a focus of discussion as the broad application of hydraulic fracturing to tight formations. [14][15][16] Therefore, investigation on the thermo-poroelastic responses is particularly important for heat extraction in HDR reservoirs.
In fact, the aforementioned problems involve solving a complex thermal-hydraulic-mechanical (THM) coupling process. Some small-scale experiments have been carried out to reveal the multiphysical-interaction mechanisms during heat extraction. [17][18][19] At present, numerical simulation is still the most efficient method to analyze the THM coupling problem. Based on the thermal-hydraulic (TH) coupling theory, the effects of different fracture morphology, 20,21 circulating fluids, 21,22 and conceptual models or well layouts 23-25 on heat extraction performance were investigated, respectively. In their models, it was assumed that the local thermal equilibrium (LTE) between the rock matrix and the pore fluid can be achieved instantaneously. Numerous studies have demonstrated that the local thermal nonequilibrium (LTNE) theory is more appropriate if there is a rapid heat transport or a large temperature difference between the rock matrix and the flowing fluid. [26][27][28] Considering the LTNE effect, Chen et al 29 developed a unified pipe-network method (UPM) to reveal the heat extraction mechanisms in geothermal reservoirs with rough-walled fracture networks. However, the mechanical effect was not incorporated into the above model. Fully coupled THM models have been put forward to investigate the heat extraction process of geothermal reservoirs. [30][31][32] An alternative THM coupling method has been developed by introducing the mean total stress to reduce the computational cost. 33 For fractured geothermal reservoirs, the mass and heat transport in fractures is faster than in the surrounding porous media. In this case, the effect of pre-existing fractures on the heat extraction process should be taken into consideration. Generally, there exist three models used to simulate fluid flow in fractured reservoirs: equivalent continuum medium model, dual medium model, and discrete fracture network model. [34][35][36] The discrete fracture network model is more suitable for describing the heterogeneity and anisotropy of fractured reservoirs. It can be used to approximately map the microseismic events associated with fracturing stimulation. More importantly, discrete fracture network simulations yield more pessimistic (and more realistic) predictions and can reproduce the flow-channeling phenomenon. 13 Based on the dual porosity model 37 or the discrete fracture network model, [38][39][40][41][42][43] some studies have been conducted to analyze the heat extraction in 2D or 3D EGS. These contributions can provide some Generating a complex and scattered fracture network but without preferential channels is conducive to extracting more heat from geothermal reservoirs.

K E Y W O R D S
fracture network morphology, fractured geothermal reservoirs, heat extraction performance, thermoporoelastic responses, THM coupling | 1707 HAN et Al.
reliable foundation for geothermal development. However, there are also some limitations in these coupling models such as ignoring the effect of temperature change on fluid flow, the stress-dependent fracture aperture, or the pressure-dependent pore permeability. In addition, there is a lack of detailed and comprehensive discussions on the thermo-poroelastic responses during heat extraction of geothermal reservoirs with different fracture network morphology.
In this paper, we presented a fully coupled THM model to investigate the long-term THM responses of heat extraction in an ideally fractured EGS reservoir. The effect of fracture network morphology on heat extraction efficiency was also discussed in detail. This study can provide knowledge of the underlying mechanisms of heat extraction and contribute to better heat extraction performance of the realistic EGS.

NUMERICAL METHOD
It is assumed that an idealized EGS consisting of an injection well and a production well contains a large number of induced fractures after stimulation (as shown in Figure  1). These fractures are the primary conduits for fluid circulation and heat extraction. Actually, the stimulation treatment will create a complex three-dimensional fracture network, which poses a formidable challenge in modeling the fracture distribution and solving the THM coupling equation. For the sake of the computational efficiency, a two-dimensional planar model with embedded random fractures is adopted to investigate the thermo-poroelastic responses during heat extraction (as shown in Section 4). Although the simplification may lead to inconsistency with the realistic EGS, the research results can provide some reliable foundation for geothermal development. In addition, some assumptions have to be adopted to formulate the complex process. (a) The pore space is fully saturated with water, and the pore water does not vaporize throughout the simulation; (b) Fluid flow in porous media and fractures obeys Darcy's law and Cubic law, respectively; (c) Heat transfer is mainly dominated by conduction and convection processes, and local thermal nonequilibrium between the solid phase and the fluid phase is considered; (d) The reservoir matrix is assumed to be a homogeneous isotropic medium, and the rock deformation satisfies the linear elastic planar strain theory; (e) The length of the fracture follows the fractal distribution. All pre-existing fractures within the reservoir are regarded as internal boundaries. All conservation equations are established in the porous medium and fracture domains, respectively. Changes in pore and fracture permeability associated with the thermo-poroelastic effect are elaborated by different empirical formulas.

| Modeling of random fracture distribution
Simulating heat extraction from EGS needs to develop a model to describe the distribution of pre-existing fractures in the stimulated reservoir. Regrettably, it is impossible to reproduce the actual fracture system of the buried reservoir. Geophysical investigations have confirmed that stochastics and the fractal theory are two applicable methods for approximately characterizing the fracture system. 44 Given the complexity and statistical self-similarity of the fracture distribution, the fractal geometry theory can be considered a feasible method to describe the fracture networks. 45

-D c u t p l a n e
Studies have shown that the number of fractures with a length equal or greater than l can be formulated as follows 47,48 where C is a proportional constant; D f is the fracture fractal dimension that can be measured from core observations. It is noteworthy that some studies have found that the fracture fractal dimension D f is related to the b value of Gutenberg-Richter (G-R) relationship. 49,50 The G-R model in seismology is a well-known formula used to describe the frequency-magnitude distribution of seismic events. Therefore, the fracture fractal dimension can also be obtained from microseismic event statistics.
Differentiating Equation (1) with respect to l can lead to Assume that the minimum and maximum lengths of the fractures in the zone of interest are l min and l max , respectively. The corresponding fracture number, N max min , between the specified supremum and infimum can be estimated as Based on Equation (2) and (3), a probability density function about the number of fractures can be defined as follows A distribution function, F(l), represents the ratio of the total number of fractures between l min and l to the total number of fractures counting from l min upward to l max , is defined. The function, F(l), can be written as follows Provided that l ∈ (l min , l max )√2, we can infer F(l) ∈ (0,1). It is assumed that the function F(l) conforms to a certain distribution pattern, such as uniform distribution. The maximum and minimum length of fractures, l max and l min are prescribed, respectively, then the length of the random fracture can be determined as follows In addition, the azimuth and location of all fractures also need to be assigned according to a certain distribution pattern. The next step is to define the density of the fracture network in a planar domain of interest. Finally, the visualization of the fracture distribution can be realized by using the Monte-Carlo method.

| Governing equation for mass conservation
The effect of changes in pore volume and fluid temperature on mass transport is incorporated into the mass balance equation. The final governing equation describing fluid flow in porous media can be written as follows 39,51,52 where p, T, and ε V are the fluid pressure, fluid temperature, and pore volumetric strain in porous media, respectively; C 1m = m S l + (1 − m )S m represents the storage coefficient of porous media; C 2m = m m l + (1 − m ) m represents the thermal expansion coefficient of porous media; C 3m = m represents Biot's coefficient of porous media; m is the reservoir porosity; S l and S m represent the storage coefficients of the fluid and rock matrix, respectively; l and m represent the thermal expansion coefficients of the fluid and rock matrix, respectively; l and μ are the density and dynamic viscosity of the fluid, respectively; and k m is the pressure-dependent pore permeability.
Pre-existing fractures in porous media can be regarded as internal boundaries. Comparing the fracture width (in millimeter) with the fracture length (in meter), the flow along the fracture width direction can be ignored. The tangential derivatives are used to define the flow along the internal fracture 39,51,52 where p, T, and ε V are the fluid pressure, fluid temperature, and pore volumetric strain in the fracture, respectively; C 1f = f S l + (1 − f )S f represents the storage coefficient of the fracture; C 2f = f ( f l + (1 − f ) f ) represents the thermal expansion coefficient of the fracture; C 3f = f represents Biot's coefficient of the fracture, where f is assumed to be equal to m ; f is the fracture porosity; S f represents the storage coefficient of the fracture; f represents the thermal expansion coefficient of the fracture; e h is the hydraulic aperture between the two fracture surfaces; k f is the stressdependent fracture permeability; n ⋅ Q m = n ⋅ ( − k m ∕ ∇p) is the mass flux exchange between porous media and the fracture; and ∇ T denotes the gradient operator restricted to the fracture's tangential plane.

| Governing equation for energy balance
The LTNE theory is adopted to simulate the heat exchange between the rock matrix and the flowing fluid. For rock matrix, the energy transfer process is mainly dominated by the heat conduction and the heat exchange with pore fluid. The governing equation for rock matrix satisfies 29,53 where T m and T l are the matrix and fluid temperatures in porous media, respectively; m is the density of the rock matrix; C p,m is the specific heat capacity of the rock matrix; m is the heat conductivity of the rock matrix; and q ml represents the rock matrix-pore fluid interface heat transfer coefficient. From the point of view of the domain, the ingoing heat flux is received from the thin fracture. Conversely, the outgoing heat flux n ⋅ q m leaves the domain and is received by the adjacent fracture in the source term 29,53 where T m and T l are the matrix and fluid temperatures in the fracture, respectively; f is the fracture density; C p,f is the specific heat capacity of the fracture; f is the heat conductivity of the fracture; q fl represents the fracture-fluid interface heat transfer coefficient, related to the fracture aperture 29 ; and n ⋅ q m = n ⋅ ( − (1 − m ) m ∇T m ) represents the heat flux exchange of the solid between the rock matrix and the fracture.
For pore fluid, the contribution of heat convection should be incorporated into the aforementioned energy balance equation 33,53 where C p,l is the heat capacity of the fluid at constant pressure; l is the heat conductivity of the fluid; and u m = −k m ∕ ∇p is Darcy's rate in porous media.
Analogously, the heat flux coupling relationship of the fluid between the domain and the fracture satisfies 53 where u f = −k f ∕ ∇ T p is the flow rate in fractures; the heat flux n ⋅ q l = n ⋅ ( − m l ∇T l ) denotes the heat flux exchange of the fluid between the porous media and the fracture.
Meanwhile, the temperature-dependent fluid thermodynamic properties are incorporated into the aforementioned mass and energy conservation equations to account for the bidirectional HT coupling. The changes in fluid density, viscosity, heat conductivity, and specific heat capacity with the temperature are described by some empirical formulas. 20

| Governing equation for stress equilibrium
Based on the linear elasticity, effective stress, and thermal stress theory, the thermo-poroelastic equation describing the stress-strain relationship of porous media is described as follows where ij is the total stress; G and λ are Lame's constants; tr is the trace operator; the third and fourth terms on the right hand of Equation (13) represent the poroelastic effect and the thermoelastic effect, respectively; is the bulk modulus of the drained porous media; T = m l + (1 − m ) m is the volumetric thermal expansion coefficient of porous media; and ij is the Dirac dealt function; p ( = m = f ) represents Biot's coefficient.
ij eff = ij + p p ij is the effective stress. In the context, the convention for the stress specifies the tension with positive notation.
Combined the equilibrium equation, the deformation equation of porous media can be derived as follows where f i denotes the external body force.
The change in effective stress has a significant effect on the permeability of the fractured reservoir. During geothermal mining, the injection and extraction of the fluid will result in the variation of the pore pressure. In addition, the large temperature difference between the injected cold fluid and the hot rock mass can induce thermal tensile stress. The thermo-poroelastic effect will cause the fracture closure, opening, or shear slip, which can affect mass and heat transport during heat extraction.
A model describing the opening and closure processes of the stress-dependent fracture has been developed by Barton and Bandis. 54 The change in the initial aperture of the fracture under in situ stresses can be expressed as follows where e 0 is the initial aperture of the fracture; n eff is the effective normal stress acting on the fracture surface; and nref is the effective normal stress required to cause 90% reduction in fracture aperture.
When the shear stress acting on the rough fracture surface exceeds its shear strength, the shear dilation will be induced due to the shear slippage of the fracture. The shear dilation contributes to increasing the aperture of the fracture. The shear strength of the rough fracture can be estimated by Patton's saw-tooth fracture model 55 where basic is the basic friction angle of the fracture; eff dil is the effective shear dilation angle of the fracture; and dil is the experimentally measured dilation angle.
The increase in the fracture aperture due to shear dilation can be calculated based on the shear displacement 56 where U s is the shear displacement; n is the shear stress acting on the fracture surface; K s is the shear stiffness, which can be approximately calculated as K s = 27 24 G l , where G is the shear modulus of the surrounding rock; l is the half-length of the fracture.
Therefore, for the closed fracture, the stress-dependent fracture aperture can be evaluated as 56 where e res is the residual aperture of the fracture under high effective normal stress.
When the fluid pressure in the fracture is bigger than the normal stress acting on the fracture surface, the fracture will be induced to open and become the hydraulic fracture. The pressure sensitivity of the fracture is considered by its hydraulic aperture increment associated with the change in fluid pressure as following empirical equation 57,58 where f is a constant characterizing the compliance of the fracture with respect to the fluid pressure change. Consider that the fracture aperture is in several millimeters, the compliance-related coefficient f should be of the order of ~10 −7 -10 −9 Pa −1 .
By integrating Equation (19) and introducing the effective stress, the hydraulic aperture of the fracture can be expressed as follows Assume that the parallel fracture surfaces have dissimilar morphology to account for the effect of the surface roughness on flow behavior through the fracture system. Based on the Cubic law, the equivalent permeability of the rough fracture can be obtained as follows where f is the roughness factor of the fracture.
In addition, the pressure-dependent pore permeability is also taken into consideration in this model as follows 59 where m eff = −tr( ij eff )∕3 the mean effective pressure, p is also a constant characterizing the compliance of porous media with respect to the fluid pressure change. Generally, the compliancerelated coefficient p is smaller than f due to less pressure sensitivity. Note that compression is here considered positive for the mean effective pressure.

| Numerical implementation
As aforementioned theories, the heat extraction from fractured geothermal reservoirs involves solving an intercoupled THM process. The complex THM coupling model combined with the corresponding initial and boundary conditions can be solved by the finite element method (FEM). Unfortunately, the computational cost of a fully coupled approach to solving the THM problem is always huge, especially in geothermal reservoirs with widely random fractures. Alternatively, the segregated approach, which solves each physics sequentially until convergence, is adopted to solve the THM coupling problem. Although more iterations may be required for the same problem, each loop through the segregated solution approach can be much faster than the Newton-Raphson step required for the fully coupled approach. It is also noticeably profitable to model all fractures as interior boundaries, because it eliminates the need to create a geometry with a high aspect ratio. Otherwise, a very long and narrow fracture domain would require a dense mesh consisting of great number (15) Δe n = e 0 1 + 9 n eff ∕ nref of tiny elements. Of course, it is still inevitable that the finer mesh is required in some regions where the fracture density is large. Unstructured triangular meshes are taken in this paper to accommodate different size subdomains. The meshes in some narrow subdomains are refined to improve the accuracy of the simulation. The detailed solution process for the THM coupling problem is shown in Figure 2.

| Validation of heat convection in a fracture
Precisely capturing the heat extraction from fracture surfaces is extremely important for evaluating heat recovery from fractured reservoirs. To testify the accuracy of the LTNE model proposed in this paper, a two-dimensional geometry model embedded a smooth fracture with constant aperture was used to model the heat exchange between the fracture wall and the flowing fluid (as shown in Figure 3). The solid temperatures on the up and down outer boundaries are maintained constant (T 0 ), while the left and right boundaries are thermally isolated. The initial temperature of the fluid in the fracture is identical with the initial temperature of the solid (T f0 = T s0 = T 0 ). It is assumed that the fluid injected with constant rate flows in parallel through the impermeable fracture walls. The temperature-dependent fluid thermodynamic properties are ignored. A constant heat transfer coefficient is assigned along the fracture plane. The injection fluid with low temperature (T in ) is gradually heated to a high temperature (T out ) fluid by heat convection. Another reasonable assumption for the thin fracture is that a uniform fluid temperature crosses the fracture section. The parameters used in this simulation are listed in Table 1.
There exists an analytical solution for the heat convection in a single fracture. The fluid temperature in the fracture can be obtained as follows 26 where x is the distance from the inlet end to the outlet end; h is the fracture heat transfer coefficient; and H is the height of the model. Figure 4 compares the analytical solution and numerical solution of the fluid temperature distribution along the fracture. The presence of significant good agreements for different heat transfer coefficients demonstrates the accuracy of the numerical algorithm. As the heat transfer coefficient increases, the fluid temperature rises faster, indicating that the heat convection has a profound effect on heat extraction.

| Validation of fully THM coupling analysis
This section aims to verify the THM coupling model through a case that accounts for the effect of the change in wellbore temperature on the temperature, pressure, and stress fields of the reservoir. The model considers a borehole with radius r w = 0.1 m in a saturated infinite granite reservoir, as shown in Figure 5. The outer boundary conditions are prescribed as The wellbore is suddenly cooled by water and maintained at 80°C. The wellbore wall is a free boundary. The initial pressure and temperature of the reservoir are 0 MPa and 200°C, respectively. Considering the symmetry of the model, only a quarter of the wellbore is modeled. An infinite domain is adopted to avoid the boundary effect. All parameters needed in this simulation were collected from the literatures 51,60 and listed in Table 2. It is worth noting that these parameters should be first converted to the corresponding model parameters in Section 2. This can be easily implemented based on the definition of these parameters. 51 The coupled thermo-poroelastic problem has been extensively studied by analytical and numerical methods, which can be found in literatures. 51,60 Figure 6A-D illustrates the changes in temperature, induced pore pressure, induced radial stress, and induced tangential stress around the wellbore, respectively. It can be observed that the numerical results are in agreement with the analytic solutions. This indicates the validity of this numerical approach for the case of thermos-poroelastic coupling problem under consideration. Therefore, the model proposed in this paper can be used to reveal the thermo-poroelastic responses of geothermal reservoirs due to cold fluid injection during heat extraction.

| STUDY ON HEAT EXTRACTION FOR CONCEPTUAL EGS
A typical EGS project located at Rittershoffen in Alsace, France, was launched in 2010. After drilling, the hydraulic testing found that the injectivity of well GRT-1 was insufficient for industrial exploitation. Subsequently, thermalchemical-hydraulic (TCH) stimulation was carried out to increase the connectivity between the well and the reservoir. After the TCH stimulation, the injectivity and permeability of the reservoir were greatly improved. At the end of the injection, the seismicity monitoring illuminated a planar 300 m long structure with a vertical extent of 200 m, as shown in Figure 7A. Subsequent tracer testing also confirmed that good hydraulic connectivity between well GRT-1 and well GRT-2 has been established. 8,9 Two groups of fractures with different azimuth angles are specified to fit the microseismic events in the stimulation zone. The corresponding fracture parameters are listed in Table 3. Fracture density is a variable parameter associated with microseismic events, and its influence on heat extraction performance will be discussed later. Finally, the fractal fracture networks are generated by the Monte-Carlo method. Doublet well layout with a distance of 200 m is designed in the simulation zone, as shown in Figure 7B. The same mass flux Q(r w , t) = −0.12 kg/(m·s) is assigned to the production well to establish the flow circulation, which means that the circulation loss is not currently taken into account in the performance evaluation of heat extraction. All fractures within the domain are regarded as internal boundaries, implicitly considering the mass and energy exchange between porous media and fractures.

| Basic input parameters
Most of the parameters required in this simulation were collected from the literatures 8,9,20,25,33,[38][39][40][41]56 and listed in Table 4. The temperature-dependent fluid thermodynamic properties are described by some empirical formulas and are not counted in Table 4.

| Evolution of pore pressure
Temporal and spatial evolution of fluid pressure in pores and fractures is demonstrated in Figure 8. It can be seen that a pressurized zone around the injection well is formed as the  fluid builds up in the wellbore. This is because in the early stage, the low permeable reservoir is insufficient to accommodate continuous injection of fluid. Conversely, the heat fluid extraction from the production well results in a pressure depletion zone around the production well. The hot fluid occupying the pore space is gradually replaced by the injected cold fluid. The increase in fluid viscosity will result in a bigger flow resistance in the cooled zone. It can be seen that the injection pressure gradually increases with time. An increasingly large pressurized zone moves toward the production well. It is taken for granted that the overall reservoir will eventually reach a steady state with the circulation of the fluid. It is obvious that the fluid pressure in pores is higher than that in fractures at the same location. This can be attributed to the higher flow resistance in low permeability pores. Most fractures are only partially connected, and some isolated small-scale fractures are widely distributed inside the domain (see Figure 7B). Although the seepage process is dominated by these interconnected large-scale fractures, the partially connected small-scale fractures also play a significant role in fluid circulation. It is observed that the pressurized zone around the injection well propagates more uniformly, rather than directionally, toward the production well. This result reveals that no preferential flow channels are formed between the injection well and the production well ( Figure 8D). Figure 9 illustrates the distribution of the fluid velocity vector after 40 years. It can be seen that the fluid velocity in fractures is much higher than that in porous media, especially in interconnected fractures. Heat convection will be the dominant heat transfer process in these large-scale fractures. It is also evident that the fluid velocity in some isolated small-scale fractures is of the same order of magnitude as that in high-permeability pores. Some regions that are free of fractures or enclosed by interconnected fractures are almost nonflowing (called as dead zones, u < 10 −10 m/s), where the heat convection process can be ignored and the energy transport is mainly dominated by heat conduction.

| Evolution of temperature
In the LTNE theory, there is a temperature difference between the solid phase and the fluid phase. The equivalent temperature T eq in pores and fractures is defined to account for the integrated effect of the fluid and the solid as follows Temporal and spatial evolution of the equivalent temperature is shown in Figure 10. It can be found that in the early stage, heat extraction from the hot rock mass occurs (24) Compliance-related coefficient Compliance-related coefficient primarily in these fractures connecting with the injection well. In this case, the matrix temperature drawdown is dominated by heat convection. With continuous injection of cold fluid, the rock matrix around these interconnected fractures will be further cooled down by heat convection and heat conduction. It is observed that the low-temperature zone is progressively expanding toward the production well with heat extraction time. There are no obvious preferential channels for heat transport between the injection well and the production well. The partially connected small-scale fractures within the reservoir dominate local mass and heat transport. It can be seen that the cooling of the hot rock mass toward the production well occurs along many flow paths in parallel except for Figure 10A. More pathways for heat transport contribute to avoiding the thermal short-circuiting and delaying thermal breakthrough. It can be seen that after 40a, a significant portion of the reservoir has been cooled by the injected fluid. Figure 11 illustrates the temperature evolution in the production well. The thermal breakthrough takes place at about 12.5 a, before which the matrix temperature around the production well remains constant. Although the thermal breakthrough takes place along some main fractures, the low-temperature profile is relatively uniformly propelled throughout the heat extraction process (see Figure 10C,D). The reduction in fluid temperature by only 21 K from 10a to 40a confirms the sufficient heat supply from other directions. This result indicates that most of the heat can be extracted from the geothermal reservoir. Therefore, generating a partially connected but more scattered fracture network is conducive to obtaining more heat production and prolonging the life service of the EGS. It is also evident from Figure 11 that in the production well, the temperature of the fluid is different from that of the matrix. To further reveal the LNTE effect, it can be seen from Figure 12 that the temperature difference between the solid phase and the fluid phase (up to 6 K) is relatively obvious in the fracture (see Figure 7B: fracture AB). This can be attributed to the higher flow velocity in the fracture. The temperature difference can be further enhanced in some fractures with higher flow rates. In porous media (see Figure 7B: line CD), the rock matrix and pore fluid also show a small temperature difference of about 1-3 K. This result is consistent with the result of Ref. 33 Therefore, it is more reasonable to use the LNTE theory to analyze the heat extraction from the HDR reservoir. The hot rock mass is gradually cooled by the injected fluid. The temperature difference between the rock matrix and the flowing fluid is negligible with the increase in time. This result reveals that the nonequilibrium heat transfer eventually converges to the thermal equilibrium.

| Evolution of effective stress
The effect of thermal stress on effective stress is investigated by comparing the results of the thermo-poroelastic model against those of poroelasticity. The effective stresses of thermo-poroelasticity and poroelasticity in the x and y directions are demonstrated in Figures 13 and 14, respectively. Note that tensile stress is considered positive. It can be clearly seen that the effective stress gradually increases as the cold fluid is continuously injected. The change in effective stress is more pronounced when considering the thermoelastic effect. This can be attributed to the pore pressure and thermal stress effects generated by the injected fluid. When the thermoelastic effect is considered, the thermal tensile stress is induced due to the temperature difference between the hot rock mass and the injected cold fluid. The shrinkage of the rock matrix in the cooled zone will enhance the reservoir permeability, which can further intensify the effect of poroelasticity on the effective stress. The influence of the poroelastic effect F I G U R E 1 1 Change in the production temperature at different times. Note that T l is the temperature of the production fluid, T m is the matrix temperature in the production wellbore, and T eq is the equivalent temperature on the effective stress is limited, while the thermoelastic effect can affect the overall reservoir with the circulation of the cold fluid. It is noteworthy that some outside regions (yellow ellipse), where the matrix temperature drop is insignificant, are subjected to compressive stress due to the effect of strain compatibility. 38 It can also be seen that in the early stage, the effective stress around the production well decreases slightly due to fluid extraction. The change in the magnitude of effective stress has a profound effect on reservoir permeability, which dominates the evolution of the mass and energy transport. Figure 15 illustrates the change in the direction of the maximum horizontal effective principal stress. With the injection of the cold fluid, the direction of the maximum horizontal effective principal stress rotates in some zones. The difference between the maximum and minimum horizontal effective principal stresses is relatively insignificant, only 4 MPa. Therefore, it can be seen that the reorientation of the maximum horizontal effective principal stress always takes place. This can be attributed to the heterogeneous induced stress generated by the thermo-poroelastic effect, which can easily alter the relative magnitude of the maximum and minimum horizontal effective principal stresses. The direction of the maximum horizontal effective principal stress is time-dependent, which evolves along with the contour of the lowtemperature distribution. It is noteworthy that the change in the direction of the horizontal effective principal stress can result in a strike-slip regime. 41 On the one hand, this will contribute to increasing the aperture of fractures to improve the reservoir permeability; on the other hand, the induced shear slip may trigger the seismicity. The equivalent temperature distribution along the lines XY and EF (see Figure 7B) for the poroelastic model and the thermo-poroelastic model is demonstrated in Figure  16. It can be seen from Figure 16A that at the early state, there is almost no difference in the equivalent temperature of the two models. This result reveals that heat extraction in this stage is dominated by the poroelastic effect. With the injection of the cold fluid, the thermoelastic effect gradually plays an important role in heat extraction. It can be found that in some regions significantly affected by the injected cold fluid, the equivalent temperature of the thermo-poroelastic model is lower than that of the poroelasticity (eg, up to 11.6 K at 10a). However, in other regions, the equivalent temperature of thermo-poroelasticity is higher than that of the poroelasticity (eg, up to 3.6 K at 10a). The phenomenon can be further clarified by Figure  16B. The equivalent temperature of the thermo-poroelastic model is higher before the heat reaches the mid-perpendicular line EF between XY, after which the equivalent temperature of the thermo-poroelastic model is lower. The reason is that when considering induced thermal stress, the reservoir permeability in the cooled zone is further improved, which contributes to increasing the flow rate to intensify the heat convection. Therefore, more heat can be extracted from the surrounding hot rock mass, resulting in a lower matrix temperature and bigger cooled zone. In addition, the flowing fluid in the thermo-poroelastic model is heated faster by the hot rock mass, and to a certain extent, its temperature will be equal to the rock temperature. There is no temperature difference between the rock matrix and the pore fluid to induce thermal stress. Therefore, the equivalent temperature in farther zones is higher when considering the thermo-poroelastic effect. In the later stage, the thermal breakthrough has taken place. It can be seen that there exists at most 3-4 K difference between poroelasticity and thermo-poroelasticity after 30a ( Figure  16A). These results further confirm that induced thermal stress has a significant effect on the heat extraction of the geothermal reservoir.

| Evolution of permeability
Temporal and spatial evolution of the stress-dependent permeability is demonstrated in Figure 17. It is evident that the injection of the cold fluid results in permeability enhancement around the injection well. This makes sense because the injected low-temperature fluid can induce thermal stress and pore pressure, resulting in an increase in the effective stress (note that tension is positive; see Figure 13). Consequently, pre-existing fractures may be induced to slip or even open, which brings about an increase in the fracture aperture. Meanwhile, the decrease in the mean effective pressure (note that compression is positive) will result in the permeability Conversely, it can also be seen that at the early stage, a slight decrease in permeability occurs near the production well due to fluid extraction. In the long term, the permeability of the overall reservoir is gradually improved with the supply of the injected fluid. Therefore, the stress-dependent fracture and pore permeability can be reasonably captured by the proposed model in this paper. The permeability change associated with fluid injection has a direct effect on the mass and energy transport. Especially, many improved secondary fractures contribute to extracting more heat from the hot rock mass. Figure 18 illustrates the changes in fluid pressure and fracture aperture in injection and production wells for poroelastic and thermo-poroelastic models. It can be seen that the injection pressures of the two models have a similar trend, but the increase in the injection pressure for poroelasticity is faster than the thermo-poroelasticity. At the beginning, the injection pressure raises sharply due to insufficient accommodation capacity of the reservoir. As the cold fluid is continuously injected, the increase in the fracture aperture contributes to lowering the flow resistance. Meanwhile, the increase in fluid viscosity will result in a bigger flow impedance. The combined effect of the two causes the injection pressure to maintain a relatively stable state. Since the fracture network is not fully connected, the injection pressure will continuously increase with the injection of fluid. It can be inferred that the injection pressure should maintain at a high level, to achieve the requirement of constant injection rate at 0.12 kg/(m·s). The change in injection pressure is similar to the result of Yao et al. 40 When considering the thermo-poroelastic effect, the change in the fracture aperture at the injection point is more significant. The closed fracture undergoes two stages with the change in the effective normal stress (I: shear dilation; II: induced opening). However, when ignoring the thermoelastic effect, the fracture aperture is only increased by 0.075 mm from 0a to 40a. There is little difference in production pressure between the two models. This is because the fracture aperture of the thermo-poroelastic model is at most 0.1 mm larger than that of the poroelastic model. In general, the permeability enhancement is primarily dependent on the induced stress generated by the injected cold fluid. Therefore, the fully coupled thermo-poroelastic effect must be taken into account when investigating the heat extraction of fractured geothermal reservoirs.

| Effect of fracture network morphology on heat extraction performance
In this section, we discuss the effect of fracture network morphology on heat extraction. The parameters for generating the fracture network are consistent with the basic case ( Figure 7B) except for the fracture density. Five cases with different fracture network densities and connectivity are adopted in this study, as shown in Table 5. The fracture network in Case 3-1 is a derivative of the Case 3 by slightly adjusting the connectivity. The resulting fracture network for the five cases is shown in Figure 19. The required simulation parameters are also consistent with those in Section 4.2. Some evaluation parameters are defined to characterize the long-term heat extraction performance of the EGS. The total heat extraction ratio η is defined as the extracted thermal energy from the EGS divided by the heat stored in the EGS where T m,init is the initial temperature of the reservoir; T m (t) is the reservoir temperature at time t; T in is the injected fluid temperature; and S represents the area of the reservoir. In this equation, (T m,init − T m (t))∕(T m,init − T in ) represents the local heat extraction ratio.
Thermal breakthrough time (t b ) is defined as the critical point of the stable production, after which the temperature of the production fluid begins to decline obviously.
The service time or lifetime of the EGS (t s ) is defined as the operation time, at which the production fluid temperature reaches the abandonment temperature. Currently, different criteria of abandonment temperature have been used in some literatures. These criteria generally range from 5% to 15% decline of the production temperature compared with its maximum value. 61 According to Ref., 33 the abandonment temperature below the maximum production temperature of 10 K is adopted in this study to determine the service time of the EGS.  Figure 20 illustrates the evolution of the local heat extraction ratio of different fractured geothermal reservoirs. It can be seen that the morphology of the fracture network has an extremely significant effect on the heat extraction process. For Case 1, at the commencement of the production phase, the heat extraction along these fractures connected with the injection well is faster than the other cases. The heat extraction profile tends to propagate more evenly in later stages. This is because the same injection rate is distributed into fewer fractures such that each fracture carries a greater flow rate. In addition, the lower fracture density (inadequate SRV) will result in a bigger hydraulic impedance. A direct result is that an extremely high-pressure gradient is required to drive the injection fluid flow from the injection well to the production well as shown in Figure 21. At the high injection pressure, the injected fluid diffuses faster and more evenly to the production well due to the poroelastic effect. By increasing the fracture density (as shown in Case 2 and Basic case), the heat extraction process will occur in parallel along widely distributed fractures. Compared with Case 1, although the initial heat extraction ratio is relatively low, it is even higher in the later stage. In addition, the injection-production pressure difference decreases significantly with the increase in the fracture density (see Figure 21). Therefore, an undeveloped fracture network may result in the inadequate injectivity. For Case 3, it can be seen that some fractures in the reservoir are fully interconnected, the preferential pathways formed dominate the mass and heat transport of the injection fluid. The fluid flows   primarily along interconnected fractures, which results in a smaller low-temperature sweep zone and the earlier thermal breakthrough. By reducing the connectivity of the fracture network in Case 3, it is obvious that the heat extraction scope of Case 3-1 is obviously larger than that of Case 3. Although thermal breakthrough has taken place along some main (but nonpreferential) channels, those partially connected fractures still contribute to the later heat supply. The result further reveals that generating a complex and scattered but not fully interconnected fracture network is conducive to extracting more heat from geothermal reservoirs. Figure 22 illustrates the effect of fracture network morphology on the production fluid temperature and the total heat extraction ratio. Compared with the Basic case, when the fracture density is relatively lower, the temporal evolution of the production fluid temperature shows three distinct stages. The first stage associated with the thermal breakthrough time shows a stable production temperature. After this stage, the production temperature decreases at a fast rate, indicating insufficient thermal compensation from the surrounding rock mass. In the third stage, the production temperature gradually becomes relatively stable with the supply of the heat, especially for Case 1. When the fracture density is larger compared with the basic case, it is evident that the temperature drops sharply after thermal breakthrough (Case 3). The negative effect can be improved by changing the connectivity of the fracture network (Case 3-1). Properly increasing the fracture density can delay the thermal breakthrough time, but an excessive fracture intensity will result in a significant advance in thermal breakthrough, especially for the fully connected fracture network (ie, Case 3 about 1.5a). The density and connectivity of the fracture network also have a profound effect on the lifetime of the EGS. For the five cases (see Table 5), the lifetimes are about 26.5a, 14a, 17a, 5.5a and 11a, respectively, and the corresponding heat extraction ratios are 0.57, 0.63, 0.59, 0.05 and 0.17, respectively. By properly changing the density and connectivity of the fracture network, it can be seen that the production performance of the EGS including the thermal breakthrough time, service time, and heat extraction ratio can be obviously improved. Therefore, when stimulating the geothermal reservoirs, it is necessary to avoid generating  preferential channels between the injection well and the production well that may result in the thermal short-circuiting. Of course, it is noteworthy that the undeveloped fracture networks or inadequately stimulated geothermal reservoirs will also be detrimental to the heat extraction.

| CONCLUSIONS
A fully coupled THM model with embedded random fractures was developed to investigate the thermo-poroelastic responses during heat extraction of geothermal reservoirs. The effect of fracture network morphology on heat extraction efficiency was analyzed based on production temperature, thermal breakthrough time, service time, local heat extraction ratio, and total heat extraction ratio. Some remarkable conclusions are summarized as follows.

1.
A pressurized zone located around the injection well is generated as the injection fluid builds up in the wellbore. Conversely, a pressure depletion zone around the production well is induced due to fluid extraction. The fluid velocity in interconnected large-scale fractures is much higher than that in porous media. However, in some isolated small-scale fractures, the flow rate is comparable to the flow rate in porous media. In the early stage, interconnected large-scale fractures are the main seepage pathways, but the partially connected small-scale fractures will play a dominant role in the later fluid circulation. 2. The generated low-temperature zone progressively moves toward the production well as the cold fluid is continuously injected. Heat extraction from the hot rock mass is dominated by different heat transfer regimes. In interconnected fractures, the matrix temperature drawdown is dominated by heat convection. In some isolated matrix regions, the thermal transport by conduction is the dominant process. There is a temperature difference between the rock matrix and the flowing fluid, especially in fractures. It is more practical to use the LNTE theory to analyze the heat extraction process in the fractured geothermal reservoir. 3. Comparing the thermo-poroelastic model with the poroelastic model, it can be observed that the thermoelastic effect will intensify the effect of poroelasticity on the effective stress. The thermo-poroelastic effect can significantly alter the direction and magnitude of the effective stress. On the one hand, the increase in the effective stress ("positive in tension'') will contribute to improving the permeability of the overall reservoir. Consequently, the evolution of the reservoir equivalent temperature and injection pressure in the thermo-poroelasticity is obviously different from the poroelasticity. On the other hand, the change in the effective stress can induce fracture shear slip, which can trigger microseismic events. The fully coupled thermo-poroelastic effect should be taken into account to reveal the underlying mechanisms during the heat extraction process of fractured geothermal reservoirs. 4. The stress-dependent fracture aperture and pressure-dependent pore permeability can be reasonably captured by the proposed model in this paper. It is evident that the injection of the cold fluid results in an enhancement of the overall reservoir permeability. This is mainly attributed to the induced thermal stress and pore pressure. It is also apparent that fractures are much more sensitive to effective stress than porous media. It is of great importance to evaluate the evolution of the permeability under in situ stresses, since it has a significant effect on the mass and heat transport in fractured geothermal reservoirs. 5. The morphology of the fracture network has a profound effect on heat extraction performance and injectivity.
Undeveloped fracture networks will result in a much higher injection pressure, earlier thermal breakthrough, and shorter service time, but higher heat extraction ratio. Conversely, adequate fracture networks can improve the injectivity. A fatal disadvantage is that fully interconnected fracture networks may bring about an earlier thermal breakthrough, shorter service time, and lower heat extraction ratio. Properly increasing the fracture density can delay the thermal breakthrough time, prolong the service life, and improve the heat extraction performance. Generating a complex and partially connected fracture networks without preferential channels is conducive to extracting more heat from geothermal reservoirs.