The Influence of Coupled Thermomechanical Processes on the Pressure and Temperature due to Cold Water Injection into Multiple Fracture Zones in Deep Rock Formation

A technique to produce geothermal energy from deep rock formations at elevated temperatures consists of drilling two parallel deep boreholes, the second of which is directed so as to intersect a series of fractures produced by hydraulic fracturing in the first borehole. Then, the first borehole is used for injection of cold water and the second used to produce water that has been heated by the deep rock formation. Some very useful analytical solutions have been applied for a quick estimate of the water outlet temperature and injection/production pressures in this enhanced geothermal system (EGS), but they do not take into account the influence of thermomechanical and hydromechanical effects on the time evolutions of the pressure and temperature. This paper provided help for the engineering design of the EGS based on these analytical solutions, by evaluating the separate influences of the thermal (T), hydromechanical (HM), thermo-hydro-mechanical (THM) effects on the fluid pore pressure and temperature. A thermo-hydro-mechanical (THM) model was developed to simulate the heat extraction from multiple preexisting fracture zones in the hot rock formation, by considering permeability changes due to the injection pressure as a function of changes in the mean effective stress. It was found that the thermal effects (without coupling with mechanical effects) led to a decrease of the transmissivity of the fracture zones and a consequent increase in the injection pressure, by a maximum factor of 2. When the temperature is constant, the influence of the hydromechanical effects on the fluid pore pressure was found to be negligible, because in such scenario, the variation of the mean effective stress was 3MPa, which was associated with a maximum increase in the initial permeability of the fracture zone only by a factor of 1.2. Thermo-hydro-mechanical effects led to a maximum increase in the permeability of the fracture zones of approximately 10 times the initial value, which was associated with a decrease in the fluid pore pressure by a maximum factor of 1.25 and 2, when hydrological and thermohydrological effects were considered, respectively. Changes in temperature were found not to be affected significantly by the thermomechanical and hydromechanical effects, but by the flow rate in the fracture zones. A sensitivity analysis was conducted to study the influence of the number, the initial permeability, the elastic modulus and the residual porosity of the fracture zones, and the elastic modulus of the confining intact rock, on the simulation results. The results were found to be the most sensitive to the number and the initial permeability of the fracture zones.


Introduction
An enhanced geothermal system (EGS) enables the extraction of geothermal energy under conditions where conventional production is not economic [1,2]. This is the case for a site where the reservoir rock at depth is hot but has insuffi-cient permeability for fluid production. The EGS technology consists of creating new tensile fractures or shear reactivation of preexisting fractures through water injection into a deep borehole and consequently increasing the permeability of the reservoir. The heat exchange surface for thermal extraction can be substantially increased by creating a simulation zone, in which a network of preexisting fractures is dilated through shear reactivation [3]. Then, a second borehole is drilled to intersect these fractures. Thus, water circulating down one borehole, through the fractures, and up the other would carry off heat from the hot rock to the surface. For the design of this type of geothermal systems, two main issues need to be addressed: (1) time evolution of the temperature of the water at the production borehole and (2) time evolution of the water pressure at the injection and production boreholes. The production of geothermal energy is not effective if the difference in the fluid pore pressure between the injection and production borehole is too high or if the temperature at the production borehole decreases too much during the expected lifetime of the geothermal power plant. In the first case, pumping cost to generate the geothermal power will be too high, and in the second case, the produced water will be cooled down too quickly for use as a source for geothermal energy.
Analytical solutions are available to address the two main issues identified above. An existing analytical solution for temperature is based on the use of a simplified model of parallel fractures created by hydraulic fracturing in the hot rock through which injected cold water can flow to the production borehole and is heated up along the way [4]. This solution has proved to be very useful for a first-order evaluation of the temperature evolution of the produced water with time. The classic Theis equation [5] enables to calculate the water pressure at the injection and production boreholes, due to water injection and production. Thus, these two analytical solutions can be used for a quick estimate of the water outlet temperature and injection/production pressures in the EGS. However, they do not take into account factors such as stress-induced permeability changes in the fractures during injection, temperature-dependent fluid, and rock properties, in situ stresses and elastic properties of the intact rock and fracture zones. Analytical solutions for the mechanical perturbations induced by nonisothermal injection into a permeable medium are available [6], but they are restricted to unidirectional and radial reservoirs, neglect coupling effects such as variations in viscosity due to changes in temperature or variations in pressure driven by thermal strains, and are difficult to apply when the reservoirs are surrounded by rock with different hydromechanical properties. Due to the limitation of the analytical solutions, the study of the coupled thermo-hydro-mechanical processes induced by cold water injection often requires the development of numerical models. Recent numerical modelling studies consisting of hydromechanical and thermo-hydro-mechanical simulations due to the injection of cold water/CO 2 in hot reservoirs through an injection borehole intersecting preexisting fractures or fault zones have been conducted [7][8][9][10]. To the authors' knowledge, no numerical study of separate coupled thermomechanical, hydromechanical, and thermohydromechanical processes due to injection of cold water in multiple fractures in hot rocks intersected by an injection and production borehole doublet has been done.
Because the time evolutions of the fluid pore pressure and temperature are essential for the design of this type of EGS and the very useful analytical solutions have the limitations described above, we conduct a numerical study, based on a simple thermo-hydro-mechanical (THM) numerical model, to understand the influence of separate thermal (T), hydromechanical (HM), and thermo-hydro-mechanical (THM) effects on the fluid pore pressure and temperature, due to injection of cold water in multiple fracture zones in hot rocks. Such an understanding would allow estimates of these coupled thermal and mechanical effects that can be used to complement or quantitatively qualify the results of the analytical solutions for a quick assessment of the pressure and temperature of water out of the EGS as a function of time. This paper is aimed at contributing for the design of this type of geothermal systems from an engineering perspective.
In the next section, the analytical solution [4] for the temperature of produced water is briefly presented. Then a coupled thermo-hydro-mechanical (THM) model is developed within the framework of TOUGH-FLAC [11,12]. Recent studies have been made on modelling of hydraulically induced fracture propagation [13,14]. In this paper, the fracture zones are assumed to be already created by hydraulic fracturing and hydroshearing, and hence, those mechanisms are not simulated. Changes in the fracture initial permeability are mainly due to changes in the stress normal to the fracture zones. The numerical model is first verified against the analytical solution for temperature, without considering the couplings. Then, in the following sections, results of the effects of H, TH, HM, and THM couplings are presented. A sensitivity study is conducted to study the influence of some key parameters, namely, the number, the initial permeability, the elastic modulus and the residual porosity of the fracture zones, and the elastic modulus of the confining intact rock, on the simulation results. The paper is concluded with some remarks.

Analytical Solutions
2.1. Analytical Solution for Temperature. In [4], an analytical solution for heat extraction from multiple fractures in hot rock formation is presented. This solution is based on a linear model involving an infinite series of parallel, equidistant, vertical fractures of uniform thickness, separated by blocks of homogeneous and isotropic, impermeable rock, the width of individual fracture being assumed to be negligible in comparison with spacing between the fractures (see Figure 1).
Owing to the spatial periodicity of the temperature field is the parallel fracture system; it is possible to replace the infinite system by a finite one consisting of a single vertical fracture between two matrix blocks with an insulating outer boundary at a distance from the midplane of the fracture equal to half the fracture spacing. As illustrated in Figure 1, a rectilinear coordinate system is placed such that the x = 0 plane coincides with the midplane of the fracture. In this model, water is injected at z = 0 and is flowing upward in the fracture (see Figure 1).
The result for the dimensionless outlet water temperature T WD at a distance z from the injection point can be expressed in a general form as a function of dimensionless parameters half-fracture spacing x ED and time t D ′ . These dimensionless parameters depend on the volumetric flow rate Q per fracture per unit thickness of the system, the distance z between the 2 Geofluids injection and production boreholes, the half-fracture spacing x E , the rock thermal conductivity K R , the rock-specific heat c R , the water density ρ w , the water-specific heat c w , the initial rock temperature T R0 , the temperature T W of the outlet water, and the temperature T w0 of the injected water. Figure 2 shows the dimensionless water outlet temperature T WD versus dimensionless time t D ′ , which can be used for a quick estimate of T WD [4].

Analytical
Solution for Pressure. In [5], an exact analytical solution for the transient drawdown in an infinite uniform confined aquifer is presented. The analytical solution of the drawdown as a function of time and distance is expressed by where h 0 is the constant initial hydraulic head, Q is the constant flow rate abstracted from the borehole, S is the aquifer storage coefficient, y, z is the distance (in the plane of the fracture) at any time after the start of pumping, T is the aquifer transmissivity, and t is the time. The storage coefficient of a confined aquifer is a dimensionless parameter defined as the volume of water released from storage per unit surface area of the aquifer per unit decline in a hydraulic head. This coefficient is calculated by where S s is the specific storage and b is the thickness of the aquifer. The specific storage is the volume of water that a unit volume of aquifer releases from storage under a unit decline in the head. This parameter is related to the compressibility of the water and the aquifer, according to where ρ is the density of the water, g is the gravitational acceleration (9.81 m/sec 2 ), α is the compressibility of the aquifer, n is the porosity, and β is the compressibility of the water (4:4 × 10 -10 m sec 2 /kg).

Coupled THM Numerical Model
3.1. TOUGH-FLAC Code. To study coupled thermohydromechanical (THM) processes in rock formation, the TOUGH-FLAC [11,12] is used. It is a numerical simulator linking a finite-volume multiphase flow code TOUGH2 [15] and a finite-difference geomechanical code FLAC3D [16]. In a TOUGH-FLAC analysis of coupled thermohydromechanical problems, TOUGH2 and FLAC are executed on compatible numerical grids and linked through external coupling modules, which serve to pass relevant information between the field equations that are solved in the respective codes. A TOUGH-to-FLAC link takes multiphase pressure, saturation, and temperature from the TOUGH2 simulation and provides the updated temperature and fluid pore pressure information to FLAC3D. After data transfer, FLAC internally calculates thermal expansion and effective stresses. Finally, a FLAC-to-TOUGH link takes the element stress and deformation from FLAC3D and updates the corresponding element porosity, permeability, and capillary pressure to be used by TOUGH2. A separate batch program controls the coupling and execution of TOUGH2 and FLAC3D for the linked TOUGH-FLAC simulator. The calculation is stepped forward in time with the transient thermohydraulic analysis initialized in TOUGH2, and at each time step or at a TOUGH2 Newton iteration level, a quasistatic mechanical analysis is conducted with FLAC3D to calculate stress-induced changes in porosity and intrinsic permeability. The resulting thermo-hydro-mechanical analysis may be explicitly sequential, meaning that the porosity and permeability are evaluated only at the beginning of each  3 Geofluids time step or implicitly sequential, with permeability and porosity updated on the Newton iteration level towards the end of the time step, using an iterative process. In this paper, because the thermo-hydro-mechanical changes are relatively slow, the sequentially explicit solution is used.

Model to Consider Stress-Induced Changes in
Permeability due to Thermo-Hydro-Mechanical (THM) Effects in the Fracture Zones. Two coupling approaches are usually used within the TOUGH-FLAC framework to consider stress-induced changes in porosity and permeability in fractured rocks [17]. The first approach enables to consider porosity and permeability changes in petroleum reservoirs as a function of changes in the volumetric strain [18]. This model relates first the porosity ϕ at a given stress to the isotropic volumetric strain variation ε v and then the permeability k at a given stress to changes in porosity, according to the following equations: where ϕ i is the initial porosity, k i is the initial permeability, and n is a power law exponent. Changes in the volumetric strain result from changes in the effective stress and temperature [19]. The empirical relation between permeability and porosity expressed in Equation (5) has been shown to be widely applicable to geological materials. The power law exponent n can vary between 3 and 25 for consolidated geological materials [20,21]. This model has been applied to model stress-induced changes in porosity and permeability in fault zones with a power law exponent n of 15 [17,[22][23][24].
In the second approach, changes in permeability and porosity are function of the changes in the effective mean stress. Changes in the effective mean stress result from changes in the volumetric strain, due to changes in fluid pore pressure and changes in the temperature. This model is based on laboratory experiments done on porous sedimentary rock [25]. According to the coupling model, the porosity, ϕ, is related to the mean effective stress σ M ′ (see Figure 3) according to the following equation: where ϕ 0 is the porosity at zero stress and ϕ r is the residual porosity at high stress. The permeability is correlated to the porosity according to the following exponential function [25]: where k 0 is the permeability at zero stress. This approach has been applied to model stress-induced changes in porosity and permeability in fault zones [17]. In fracture zones, an approach based on the effective normal stress is applied in [8]. The models described above are empirical, and they are developed for sedimentary rocks. The mean effective stress approach is sensitive to the initial and residual porosity, while the volumetric strain approach is sensitive to the initial porosity and the power law exponent n. Both approaches can be used to obtain an estimation of the pressure in Dimensionless water outlet temperature Figure 2: Dimensionless water outlet temperature T WD versus dimensionless time t D ′ (extracted from [4]).  Geofluids fractures or fault zones, but they should be calibrated against site-specific data for an accurate representation of the stressinduced changes in porosity and permeability. Simulations done with the volumetric strain approach for a time period of 30 years, by using a power law exponent n of 3 and 15, lead to a maximum ratio between the final and initial fracture permeability less than 2. This result is because (1) the displacement normal to the outer boundaries of the model used to represent the multiple fracture zones (see Figure 1) is restricted and (2) the fracture zones are surrounded by intact rock which is stiffer. These changes in the fracture permeability are very small, given that in the literature [8], cooling effects may lead to an increase in the fracture permeability up to two orders of magnitude. Hence, the mean effective stress approach is used to model the stress-induced changes in fracture permeability. We choose an initial and residual porosity of the fracture zones, as such we get an increase in the initial fracture permeability by one order of magnitude. Further, to evaluate the influence of a larger increase in the fracture permeability (approximately two orders of magnitude), a sensitivity study on the residual porosity of the fracture zones is made. Results of our simulations show that the effect of the changes in the permeability of the intact rock on the fluid pore pressure is negligible, and hence, they are neglected.
3.3. Thermo-Hydro-Mechanical (THM) Design. A thermohydromechanical (THM) model is implemented within the framework of the existing TOUGH-FLAC simulator [11,12] (see Figure 4). The model is based on the concept of extracting heat from multiple fracture zones in hot rock (see Figure 1). Hydraulic fracturing and hydroshearing are not simulated.
The model has 2000 m by 2000 m by 110 m and considers the terrain between the depths of 6000 m and 8000 m. A fracture zone, with a square region of 1000 m by 1000 m, is located in the centre of the model. The x-coordinate axis is perpendicular to the fracture zone, and the z-coordinate axis is vertical. An injection and production boreholes are located perpendicular to the fracture zone. They are placed symmetrical to the centre of the fracture zone, so that the distance between each borehole and the lateral boundaries is the same. By taking that as the starting point, then we evaluate the changes of EGS operational results due to various coupled analyses. The thickness of the fracture zone can be several tens of meters [26]. In our case, a thickness of 10 m for the fracture zone is assumed, which is of the same order of magnitude as in [27]. The fracture zone is surrounded by 50 m of intact rock on both sides, which implies that, according to Figure 1, the fracture spacing (calculated from fracture centres) is 110 m. A mesh with 19800 elements over the whole domain and an element of 20 m by 10 m at the location of the injection and production boreholes is used. Similar mesh refinement is used in [28].
Null displacements are set normal to the six surfaces of the model. The boundaries of the model are assumed to be closed to fluid flow and heat transfer. In the two boundaries perpendicular to the x-coordinate axis, the boundary conditions of zero normal displacement, fluid flow, and heat transfer are justified by the symmetry conditions shown in Figure 1. Regarding the other four boundaries, results of our simulations show that for larger model domains, the solution for stresses and fluid pore pressure is found not to change significantly. In particular, it is found that if the upper boundary is placed at the terrain surface, the results are very similar to those presented in this paper. Note that the simulated geothermal system is constituted by an injection and production borehole doublet, and in this scenario, the influence of the boundary conditions on the fluid pore pressure around the boreholes is much less than that for only one borehole. The boundary conditions for fluid flow and heat transfer are the same as those used in [28,29].
The input parameters are shown in Table 1. The hydromechanical properties (density, elastic modulus, Poisson's ratio, specific heat, conductivity, and thermal expansion  10 -5 10 -5 5 Geofluids coefficient) of the intact rock are typical from crystalline formations and are extracted from standard reference materials (http://www.engineeringtoolbox.com). The permeability of the intact rock of 10 -18 m 2 is characteristic from crystalline formations at the depth of 7000 m [30][31][32]. The porosity is assumed to be 0.02, which is typical from crystalline formations [33]. The fracture zone has different hydraulic and mechanical properties. It is assumed that the fracture zone is more fractured than the surrounding intact rock, and therefore, it has a lower elastic modulus and higher permeability. The elastic modulus of the fracture zone is assumed to be equal to 10 GPa, which is of the same order of magnitude of the value used by [27,34]. The initial permeability of the fracture zone is assumed to be 10 -14 m 2 , which is of the same order of magnitude of the value used by [27,30]. The correspondent hydraulic conductivity of the fracture zone is 4:69 × 10 -7 m/s, which is obtained by considering the values of the water density (951.35 kg/m 3 ) and the dynamic water viscosity (0.0001989134 Ns/m 2 ) for the initial fluid pore pressure and temperature at the centre of the fracture zone. Given that the thickness of the fracture zone is 10 m, the transmissivity of the fracture zone is equal to 4:69 × 10 -6 m 2 /s. The porosity, density, specific heat, and conductivity of the fracture zone are assumed to be equal to those of the intact rock. A residual porosity of 0.019 is chosen, which is found to lead to a maximum increase in the permeability of the fracture zone of approximately one order of magnitude.
The ratio between the horizontal and vertical stresses is assumed to be 1.0. The temperature ranges between 132°C and 168°C, and its variation with depth is linear. The geothermal gradient of 18°C/km, observed in several sites in Sweden, is used [35]. The fluid pore pressure corresponds to the hydrostatic pressure. At the depth of the injection and production boreholes (7000 m), the initial fluid pore pressure and temperature are 68.7 MPa and 150°C, respectively. The distance between the injection and production boreholes is assumed to be 750 m. Cold water with a temperature of 47°C is injected at a constant total flow rate of 120 l/s. Heat is extracted at the production borehole at a production rate of 120 l/s. This flow rate is divided into the number of considered fracture zones. Further, a sensitivity study is conducted to evaluate the influence of the number, the initial permeability, the elastic modulus and the residual porosity of the fracture zones, and the elastic modulus of the confining intact rock, on the simulation results. The numerical simulations are run for a maximum period of 30 years. The initial time step is very small (≈0.1 seconds), and then it increases with time during the numerical simulations, which is done automatically by the TOUGH-FLAC simulator.

Verification of the Numerical Model against Analytical Solutions
In this section, we verify the use of our THM model for temperature and pressure calculations. To verify our model for a meaningful decrease in temperature at the production borehole, a flow rate of 12 l/s is used, which is equivalent to consider a number n F of fracture zones of 10. The comparison of the analytical solution for temperature provided by [4] with the results provided by our THM model at the production borehole is presented in Figure 5. The figure shows that at 30 years, the temperature decrease is approximately 10°C. The difference between the analytical and numerical solutions is less than 0.5°C. This result shows that our modelling can be used for the temperature calculation, with good accuracy. The comparison between the analytical solution for fluid pore pressure provided by the Theis equation [5] with the results provided by our model shows that the difference between both solutions is approximately 20% at location of the boreholes and less than 5% away of those. This is expected because at the boreholes, the fluid pore pressure gradient is very high and the model gives a constant value of fluid pore pressure over an element with 20 m by 10 m size. A much more refined mesh, with elements of a few centimetre size around the boreholes, would be needed to improve the accuracy of the fluid pore pressure solution, but simulations done with a logarithmic mesh and such refinement lead to problems related to the time step convergence. However, the used mesh is good enough for the purpose of this paper, which is not focused on absolute values but on changes in fluid pore pressure due to various coupled analyses. We are interested in relating the results obtained with and without considering stress-induced changes in permeability. The ratio between those two results does not depend as strongly on the size of the elements at the location of the injection and production.

Study of the Fluid Pore Pressure and Temperature including T, HM, and THM Effects
The analytical solutions for pressure and temperature do not consider the influence of hydromechanical (HM) and thermomechanical (TM) effects on the fluid pore pressure and temperature. In addition, the analytical solutions consider that the confining intact rock is impermeable. In this section, we use our numerical model presented in Figure 4 to analyse the influence of the permeability of the intact rock (by considering only hydrological H effects), the effects of   6 Geofluids temperature (T), and the hydromechanical (HM) and thermo-hydro-mechanical (THM) effects on the fluid pore pressure and temperature. In this set of calculations, we consider a flow rate of 6 l/s in the fracture zone in the model presented in Figure 4. This is equivalent to consider a number n F of fracture zones of 20, with a total flow rate of 120 l/s. Note that, as demonstrated further, a flow rate of 12 l/s (10 fracture zones), as used for the verification of our model against analytical solution for temperature, would lead to an extreme high and totally unrealistic increase in the fluid pore pressure, and in such scenario, the geothermal power plant is not efficient.

Study of the Influence of the Permeability k IR of the Intact
Rock on the Fluid Pore Pressure. This section is aimed at analysing the influence of the permeability of the intact rock on the fluid pore pressure. This analysis is conducted by considering only hydrological (H) effects. In this way, the temperature is kept constant, and the stress-induced changes in permeability due to thermal or mechanical effects are neglected. The permeability k IR of the intact rock is set to 10 -14 m 2 , 10 -16 m 2 , 10 -18 m 2 , and 10 -21 m 2 . The variation with time of the difference in the fluid pore pressure between the injection and production boreholes is presented in Figure 6.
The results show that the difference in the fluid pore pressure is constant after a short period of injection (approximately 10 days when k IR is 10 -18 m 2 or 10 -21 m 2 ). This result is also obtained with the Theis equation [5], which gives a sharp increase in the fluid pore pressure difference during the initial period of injection and production. When the permeability of the intact rock is equal to the permeability of the fracture zone, the difference in the fluid pore pressure at the injection and production boreholes is only 5 MPa, because in such scenario, the fluid flow is tridimensional. When the permeability of the intact rock decreases 2 orders of magnitude, from 10 -14 to 10 -16 m 2 , as expected, the difference in the fluid pore pressure at the injection and production boreholes increases substantially (approximately 20 MPa). The results obtained with a permeability of the intact rock of 10 -18 m 2 and 10 -21 m 2 are found to be practically identical.

Study of the Influence of the Thermal (T) Effects on the
Fluid Pore Pressure and Temperature. This section aims to analyse the influence of the thermal (T) effects on the fluid pore pressure and temperature due to thermally induced changes in fluid properties, such as viscosity and density. Stress-induced changes in permeability, due to thermomechanical (TM) and hydromechanical (HM) effects, are neglected. Figure 7 shows the variation with time of the difference in the fluid pore pressure between the injection and production boreholes when the hydrological (H) and thermohydrological (TH) effects are considered. When the changes in the temperature are considered, at 1 and 30 years of injection, the difference in the fluid pore pressure is approximately 1.4 and 1.6 times that difference obtained by considering the H effects only, respectively. This is because when changes in the temperature are allowed to occur, the decrease in the temperature around the injection borehole, caused by the cooling of the rock, leads to an expected increase in the water viscosity, which in turn leads to a decrease in the transmissivity of the fracture zone and a consequent increase in the fluid pore pressure. At the production borehole, it was found that over the 30 years of injection, the temperature remains practically unchanged.

Study of the Influence of the Hydromechanical (HM)
Effects on the Fluid Pore Pressure. This section is aimed at analysing the influence of the hydromechanical (HM) effects on the fluid pore pressure. The temperature is kept constant, and the stress-induced changes in permeability due to thermomechanical effects are neglected. Changes in permeability are due only to hydromechanical effects. Figure 8 shows the difference in the fluid pore pressure between the injection and production boreholes obtained with consideration of hydrological (H) effects and hydromechanical (HM) effects. The figure shows that the results are very similar. This result   7 Geofluids is because when the fluid pore pressure increases, the total stress also increases, because the displacement normal to the outer boundaries of the model is restricted, which results in a small maximum variation in the mean effective stress (approximately 3 MPa). In such scenario, the maximum increase in the initial permeability is approximately 1.2 times, which is not enough to cause significant changes in the fluid pore pressure.

Study of the Influence of the Thermo-Hydro-Mechanical (THM) Effects on the Fluid Pore Pressure and Temperature.
This section is aimed at analysing the influence of the thermo-hydro-mechanical (THM) effects on the fluid pore pressure and temperature. Figure 9 shows the variation with time of the difference in the fluid pore pressure between the injection and production boreholes obtained with consideration of the hydromechanical (HM) and thermo-hydromechanical (THM) effects. The figure shows that when the temperature is allowed to change (THM case), the difference in the fluid pore pressure starts to decrease approximately at 1 year of injection. At 30 years, this decrease is approximately 5 MPa, comparative with the case where the temperature is constant (HM case). This is because when the temperature around the injection borehole decreases, caused by the cold water injection, the effective stress decreases significantly (approximately 22 MPa). This decrease results in a maximum increase in the permeability of approximately 10 times the initial value. The comparison of the results presented in Figures 7 and 9 enables us to conclude that the mechanical deformation affects the fluid pore pressure by a factor of about 2. When THM effects are considered, it is found that at the production borehole, at 30 years of injection, the temperature is practically equal to the initial temperature (150°C). The temperature is not much affected by the stress-induced changes in permeability, because it depends mainly on the flow rate rather than the permeability of the fracture zone.

Sensitivity Study
This section presents the results of a sensitivity study on the influence of several key parameters on the simulation results (pressure, temperature). Those parameters are the number n F of fracture zones (in the context of the temperature model [4], see Figure 1), the initial permeability k F of the fracture zones, the elastic modulus E F of the fracture zones, the residual porosity ϕ F of the fracture zones, and the elastic modulus E R of the confining intact rock. The key parameters considered in the base case and sensitivity study are shown in Table 2. The remaining parameter values are indicated in Table 1. In this sensitivity study, the thermo-hydro-mechanical (THM) effects are included. If only TH effects are considered (no changes in permeability due to mechanical effects), the results do not depend on the elastic modulus of the confining intact rock and fracture zones, and for a flow rate of 6 l/s and an initial permeability of the fracture zones of 10 -14 m 2 , they are the same as shown in Figure 7.
6.1. Effect of the Number n F of Fracture Zones in the Context of Temperature Model [4]. The number n F of fracture zones is directly related with the flow rate in each fracture zone. For a total flow rate of 120 l/s, when n F is equal to 10, 20, and 30,     8 Geofluids the flow rate is equal to 12, 6, and 4 l/s, respectively. Figure 10 shows the variation with time of the difference in the fluid pore pressure between the injection and production boreholes and the temperature at the production borehole obtained for three values of the number n F of fracture zones: 10, 20, and 30. The results show that at 30 years of injection, the difference in the fluid pore pressure is approximately 15, 20, and 40 MPa, for n F equal to 30, 20, and 10. At 30 years of injection, the permeability at the injection borehole is found to be 10 times the initial value. The temperature at the production borehole decreases more when the number of fracture zones is smaller or the flow rate in each fracture zone is larger. For a flow rate of 12 l/s, at 30 years of injection, the decrease in the temperature is approximately 10°C.
6.2. Effect of the Initial Permeability k F of the Fracture Zones. Figure 11 shows the variation with time of the difference in the fluid pore pressure between the injection and production boreholes and the temperature at the production borehole, obtained for two values of initial permeability k F of the fracture zones: 10 -13 m 2 and 10 -14 m 2 . Additional simulations are done for permeability values smaller than 10 -14 m 2 and larger than 10 -13 m 2 . Results are not presented because in the former case, a very large difference in the fluid pore pressure between the injection and production boreholes is found, and in such scenario, the geothermal power plant is not efficient. In the latter case, the increase in the fluid pore pressure is very small, and hence, the results obtained with and without stress-induced changes in permeability are very similar.
The results show that, as expected, the maximum difference in the fluid pore pressure decreases approximately 10 times when the initial permeability of the fracture zones increases one order of magnitude from 10 -14 m 2 to 10 -13 m 2 . For an initial permeability of 10 -14 m 2 , the difference in the fluid pore pressure decreases from approximately 27 MPa, at approximately 100 days of injection, to approximately 20 MPa, at 30 years of injection. For an initial permeability of 10 -13 m 2 , after approximately 1 day of injection, the difference in the fluid pore pressure is approximately constant. At 30 years of injection, the permeability at the injection borehole is approximately 10 and 8 times the initial values of 10 -14 m 2 and 10 -13 m 2 , respectively. The temperature at the production borehole is very low sensitive to the initial value of the permeability of the fracture zones. Results of our simulations show that when the stress-induced changes in permeability due to thermomechanical (TM) effects are not considered, the temperature at the production borehole does not depend significantly on the initial permeability of the fracture zones, but on the flow rate. When thermomechanical effects are coupled, the thermal expansion coefficient changes, which results in slight variations in the temperature.
6.3. Effect of the Elastic Modulus E F of the Fracture Zones. Figure 12 shows the variation with time of the difference in the fluid pore pressure between the injection and production boreholes and the temperature at the production borehole, obtained for three values of elastic modulus E F of the fracture zones: 5, 10, and 20 GPa. The results show that when the 9 Geofluids fracture zones are stiffer, the decrease in the fluid pore pressure is more significant. In this scenario, the injection pressure induced by a constant flow rate is higher, which results in a larger decrease in the effective stress and consequently larger increase in the permeability of the fracture zones. At 30 years of injection, the difference in the fluid pore   pressure is 25, 20, and 17 MPa, for E F equal to 5, 10, and 20 GPa, respectively. The influence of E F on the temperature at the production borehole is found not to be significant. The elastic modulus E F affects the permeability of the fracture zones, which does not influence significantly the temperature at the production borehole.
6.4. Effect of the Residual Porosity ϕ R of the Fracture Zones. Figure 13 shows the variation with time of the difference in the fluid pore pressure between the injection and production boreholes and the temperature at the production borehole, obtained for two values of residual porosity ϕ F of the fracture zones: 0.019 and 0.018. The figure shows that when the porosity decreases from 0.019 to 0.018, at 30 years of injection, the difference in the fluid pore pressure decreases approximately at 4.2 MPa. This result is explained by an increase in the initial fracture permeability of approximately two orders of magnitude, for ϕ F equal to 0.018, against an increase in one order of magnitude, for ϕ F equal to 0.019. The fluid pore pressure curves obtained for the two analysed cases of ϕ F start to differ approximately at 1 day of injection, which coincides with the onset of permeability changes. The influence of the residual porosity ϕ F on the temperature at the production borehole is found to be negligible, because changes in the initial fracture permeability do not induce significant changes in the coefficient of thermal expansion.

Effect of the Elastic Modulus E R of the Confining Intact
Rock. Figure 14 shows the variation with time of the difference in the fluid pore pressure between the injection and production boreholes and the temperature at the production borehole, obtained for three values of elastic modulus E R of the confining intact rock: 20, 50, and 80 GPa. The figure shows that when the intact rock is stiffer, the difference in the fluid pore pressure between the injection and production boreholes decreases more. This is because in such scenario, the tensile stress in the fracture zones caused by the decrease of the temperature around the injection borehole leads to a decrease in the mean effective stress and consequent increase in the permeability of the fracture zones. At 30 years of injection, the difference in the fluid pore pressure ranges between approximately 18 MPa and 26 MPa, when the elastic modulus E R ranges between 80 and 20 GPa. Similar to the sensitivity studies presented above, the temperature at the production borehole is not significantly affected by the elastic modulus E R of the intact rock.

Conclusions
A thermo-hydro-mechanical (THM) model is proposed to provide help for the engineering design of this EGS by evaluating the influence of the thermal (T), coupled hydromechanical (HM), and coupled thermo-hydro-mechanical (THM) effects on the fluid pore pressure and temperature changes due to cold water injection in hot crystalline formation with fracture zones.
The results of the coupled T, HM, and THM analyses show that due to the hydrological effects only, the increase in the initial fluid pore pressure is practically constant when the permeability of the intact rock is equal or smaller than 10 -18 m 2 . As expected, temperature effects (without being coupled with the mechanical deformation) lead to a change in the fluid properties, such as viscosity and density, which results in a decrease in the transmissivity of the fracture zones and an increase in the injection pressure by a maximum factor of 2. Contrary to expectation, when the temperature is constant, the influence of the hydromechanical effects on the fluid pore pressure is negligible. This is because the displacement normal to the outer boundaries of the numerical model is restricted, and due to a constant injection rate, the compressive stress normal to the fracture zones increases. This results in a small maximum variation of the mean effective stress (approximately 3 MPa), which in turn leads to a maximum increase in the permeability of the fracture zones approximately only 1.2 times the initial value. When thermo-hydro-mechanical effects are considered, changes in the temperature due to the rock cooling result in larger changes in the mean effective stress. This result leads to a maximum increase in the permeability of the fracture zones of approximately 10 times the initial value and a consequent decrease in the fluid pore pressure by an approximate factor of 1.25 and 2, when hydrological and thermohydrological effects are considered, respectively. For all the coupled analysis, when the flow rate in the fracture zones is 6 l/s, the temperature at the production borehole remains practically unchanged over a period of 30 years.
A sensitivity analysis is conducted to study the influence of the number, the initial permeability, the elastic modulus and the residual porosity of the fracture zones, and the elastic modulus of the confining intact rock, on the simulation results. It was found that, as expected, when the number of the fracture zones decreases, for the same total flow rate, the flow rate per fracture zone increases, and as a result, the difference in the fluid pore pressure between the injection and production boreholes increases and the temperature at the production borehole decreases more. For a 12 l/s in each fracture zone, at 30 years of injection, the decrease in the temperature is 10°C. As expected, when the initial permeability of the fracture zones increases by one order of magnitude, the difference in the fluid pore pressure at the injection and production boreholes decreases by one order of magnitude. When the elastic modulus of the confining intact rock or of the fracture zones decreases, the changes in the effective stress are less, which results in a less decrease in the difference of the fluid pore pressure between the injection and production boreholes. At 30 years of injection, the difference in the fluid pore pressure increases approximately 1.5 times when the elastic modulus of the confining intact rock or of the fracture zones decreases by 4 times. The chosen values of the residual porosity of the fracture zones are directly related to the changes in their initial permeability. For a flow rate of 6 l/s, an increase in the fracture permeability from one to two orders of magnitude leads to a maximum decrease in the fluid pore pressure of approximately 4 MPa, at 30 years of injection. As expected, the temperature at the production borehole is practically unaffected by the initial permeability of the fracture zones, the residual porosity of the fracture zones, and the elastic modulus of the intact rock or the fracture zones. The time evolution of the temperature depends mainly on the flow rate in the fracture zones rather than the permeability.  Figure 14: Effect of the elastic modulus E R of the confining intact rock on the time evolution of the difference in the fluid pore pressure between the injection and production boreholes (a) and the temperature at the production borehole (b).

Data Availability
No data were used to support this study.

Conflicts of Interest
The authors declare that they have no conflicts of interest.