The risks of long-term re-injection in supercritical geothermal systems

Supercritical geothermal systems are appealing sources of sustainable and carbon-free energy located in volcanic areas. Recent successes in drilling and exploration have opened new possibilities and spiked interest in this technology. Experimental and numerical studies have also confirmed the feasibility of creating fluid conducting fractures in sedimentary and crystalline rocks at high temperature, paving the road towards Enhanced Supercritical Geothermal Systems. Despite their attractiveness, several important questions regarding safe exploitation remain open. We dedicate this manuscript to the first thermo-hydro-mechanical numerical study of a doublet geothermal system in supercritical conditions. Here we show that thermally-induced stress and strain effects dominate the geomechanical response of supercritical systems compared to pore pressure-related instabilities, and greatly enhance seismicity during cold water re-injection. This finding has important consequences in the design of Supercritical Geothermal Systems.

E nhanced Supercritical Geothermal Systems (ESGS) are still widely unexplored and fundamental questions are not only to be answered, but to be formulated in the first place. Provided it is possible to stimulate permeability of the reservoirs, is it feasible to circulate fluids between two wells as in traditional Geothermal Systems (GS)? A positive answer to this question entails the second one: can fluid circulation be maintained safely and what is the potential for fluid injectioninduced seismicity in ESGS? Additionally, what could be the effect of reservoir cooling by cold water re-injection? Addressing these questions is crucial for the successful development of ESGS.
The last two decades have witnessed an increasing research and industrial interest towards geothermal energy projects in magmatic environments [1][2][3][4][5][6][7][8] : as a consequence of high temperature and depth, fluids can be found in supercritical (SC) conditions, which for pure water occurs at temperature T > 374 C and pressure p > 22:064 MPa (http://www.iapws.org/). The idea of producing geothermal power by drilling directly into magmatic environments dates back more than a decade with pioneering projects like the Icelandic Deep Drilling Projects (IDDP1 and IDDP2), at the Krafla 1 and Reykjanes 9 fields, in Iceland. Producing SC geothermal fluids can have enormous advantages: the high enthalpy per unit mass implies up to a tenfold increase in power generation 2 and low fluid viscosity implies greatly enhanced hydraulic conductivity. The latter favours advective heat transport and diminishes pressure buildup. Only few wells worldwide have encountered SC fluid conditions 6,7 and current research projects are located in active volcanic areas such as Kakkonda, Japan 10 , Taupo Volcanic Area, New Zealand 11 , Larderello, Italy 12,13 , Krafla and Reykjanes, Iceland 14,15 , The Geysers and Salton Sea, the USA [16][17][18] and Los Humeros, Mexico [19][20][21][22] .
Volcanic formations at depth are usually data-poor environments as a consequence of the technical difficulties associated with drilling exploration wells in high-temperature and highly corrosive environments. Indirect geophysical methods 5,21 in combination with advanced numerical simulations 4,8 can overcome the lack of information and provide preliminary assessments of the deep thermal structures and convective anomalies. Alongside, laboratory experiments are essential to build realistic geomechanical models based on insights from the rheological behaviour of rocks at high pressure and temperature and its effect on fluid circulation [23][24][25][26] . Close to the brittle-ductile transition conditions of pressure and temperature, new findings suggest that fractures are sufficiently permeable to allow fluid circulation 27 and, in case of insufficient fracture density, enhancement strategies are likely to be successful 28 , proving that ESGS development in low-permeable reservoirs close to the ductile crust is possible (Fig. 1). Some insights can be gained from Enhanced Geothermal Systems (EGS), in which fluid is injected in two stages 29 : a first stage in which permeability enhancement is achieved by hydrofracturing and/or hydro-shearing, chemical or thermal stimulation, and a second stage where fluids are continuously produced and re-injected to facilitate advective heat transport and reach fluid mass balance. The injection of fluids [30][31][32][33][34][35] and the associated effects of reservoir cooling 30,34,36 are linked with an increased likelihood of inducing seismic events, an occurrence that in the past has jeopardised several geothermal projects, sometimes leading to the suspension of all operations [37][38][39][40][41][42][43] . During stimulation injection, micro-seismicity is associated with hydrofracture propagation and larger events happen during prolonged injection 44 . Nevertheless, there are reported cases in which M > 5 earthquakes were likely triggered during the stimulation phase of deep geothermal systems 40,42,43 . In this study, we focus exclusively on the prolonged fluid re-injection and assume a successful stimulation of the wells' surroundings already took place. Although the connection between cooling and reservoir stability is known and numerical studies applied to traditional geothermal systems have highlighted the influence of thermal gradients 45 , fault and rock permeability 46 , stress regime 30 , competition between thermal and hydraulic processes 47 and tensile fracturing effects on permeability 48 , the issue has never been addressed in ESGS systems. Specifically, cooling in ESGS could take on a whole new dimension, as the potential temperature drop and associated thermal stress variation can be very high. The complexity and strong non-linearity of the processes involved in ESGS make them unique environments which require state-ofthe-art numerical analyses to understand the main physical mechanisms controlling their geomechanical behaviour.
Here, we investigate for the first time the feasibility and safety of prolonged fluid production and re-injection in ESGS by performing coupled thermo-hydro-mechanical non-linear analyses with innovative features such as poro-mechanical changes, porosity-dependent permeability, full range of water equation of state, non-linear fluid compressibility, fault stability and rate of seismic production. We have developed a full model of the reservoir ranging from surface to 7 km depth, applied an extreme cooling scenario considering changes in the water properties during cold water injection and provided interpretations in terms of seismicity increase and fault reactivation. We have found that seismicity increases mainly as a consequence of the thermal quenching effects and that pore pressure increments play a secondary role in controlling rock and fault stability. The thermal effects are delayed in time because of the slower heat advection processes compared to pressure diffusion. The increased viscosity during cold fluid re-injection plays a secondary role compared to thermal quenching. This work constitutes the first of a kind numerical study of the complex geomechanical processes induced in ESGS.

Modelling approach
Theoretical model. We have employed transient numerical finite element analyses of coupled thermo-hydro-mechanical (THM) processes to investigate an idealised ESGS systems formed by an injection and a production well (Fig. 2a). The non-linear system of partial differential equations describing the mass, energy and ð1Þ and is formulated in terms of primary variables describing the fields of pore pressure p, temperature T and displacement u. t is time and k is the permeability tensor. Subscripts 's' and 'w' stand for solid and for water, respectively, and I is the identity tensor. The fractured rock properties are set to values typical of a crystalline rock with porosity n ¼ 0:01, solid density ρ s ¼ 2700 kg m À3 , linear thermal expansion coefficient α s ¼ 1 Á 10 À5 K À1 , specific heat capacity c s ¼ 950 J kg À1 K À1 , thermal conductivity λ s ¼ 3I W m À1 K À1 , Young's modulus E ¼ 60 GPa, Poisson's ratio ν ¼ 0:25 and Biot's coefficient α ¼ 1 À K=K s ¼ 0:5, where K is the drained bulk modulus and K s is the intrinsic bulk modulus of the solid phase. The gravity acceleration vector is indicated by g and the source terms for the thermal and hydraulic problem are Q T and Q H , respectively. The fault has different mechanical properties with n ¼ 0:05, E ¼ 20 GPa and α ¼ 1:0. The equation of state (EOS) of water is taken after the IAPWS-IF97 standards using the freesteam application (http://freesteam.sourceforge.net/) and properties of water like density ρ w , dynamic viscosity μ w , specific heat capacity c w , thermal conductivity λ w , compressibility β w and thermal expansion α w are functions of its thermodynamic state (see Methods for the computation of fluid compressibility and thermal expansion). The relation between permeability evolution and porosity is still a matter of deep research and debate in literature. Employing the formulation of Kozeny-Carman 51,52 , in which permeability is proportional to the cube of porosity k / n 3 , is common practice 53 . The Kozeny-Carman relationship, which is formulated for flow in given shape channels 54 or assemblies of perfect spheres, takes into account tortuosity in the flow but becomes inaccurate for fractured media. Recent research empirically linked permeability of fully fractured k f and intact k i rock to porosity n, based on permeability measurements on 70 different samples of volcanic and magmatic rocks 50 . We have introduced a formulation that interpolates between intact k i and fractured k f rock permeability as log k ¼ 1 À ω ð Þlog k i þ ω log k f , in which ω is a newly defined parameter representing either fully fractured (ω ¼ 1) or intact (ω ¼ 0) rock (Fig. 2b).
The Coulomb Failure Stress CFS ¼ jτ n j þ μσ 0 n is suited to estimate potential failure 34 , where τ n and σ 0 n are, respectively, the shear and the effective normal stress acting on a given plane. We consider the orientation of the main fault (i.e., dipping 75 ) and a frictional coefficient μ ¼ 0:577, which corresponds to a friction angle ϕ ¼ 30 . ΔCFS > 0 indicates increased instability in fractures and faults with that given orientation. The mobilised friction coefficient, which indicates the proximity of the stress state to failure conditions, is μ fr ¼ τ n = Àσ 0 n À Á . The Drucker-Prager criterion is an alternative to estimate the failure potential in terms of invariants of the stress tensor (direction independence) where c 0 ¼ 30 MPa is the cohesion of the fractured rock, q dp is the Drucker-Prager deviatoric stress, σ 0 m ¼ tr σ 0 ð Þ=3 is the mean effective stress, σ 0 ¼ σ þ αpI is the effective stress tensor. The mobilised failure ratio for the Drucker-Prager model is is the acting deviatoric stress and s ¼ σ 0 À Iσ 0 m . The prediction model of time-dependent earthquake nucleation in injection-induced seismicity is based on rate and state theory of friction on faults 55,56 . The equation for the rate of seismic production relative to the background seismic noise, R, A is a constitutive parameter and equals 0:02 for hydrothermal systems 55 , _ τ 0 ¼ 1 Â 10 À3 MPa yr À1 is the background shear stressing rate, τ c ¼ CFS and the time derivative symbol is _ x ¼ dx=dt.  (Fig. 2a). We first simulate an initialisation phase to reach the conditions of a SC reservoir located above a magmatic intrusion. The bottom central part of the model is assumed to be at a constant temperature of 550 C over a width of 2 km. The hydraulic conditions are noflow at the bottom, hydrostatic pressure at the sides and a fixed atmospheric pressure at the top. The thermal conditions at the sides are fixed temperature with the average gradient of 30 C km À1 . The initialisation phase is performed in hydrothermal conditions. The model is solved until temperature and pore pressure profiles are reasonably close to what are believed to be representative conditions of a SC reservoir 8,17,[58][59][60][61] .
The injection phase, in which the coupled THM system of Eq. (1) is solved, uses the p and T distribution obtained at the end of the initialisation phase. The initial stress state is given by a lithostatic vertical stress σ V corresponding to a solid density of 2700 kg m À3 and a normal faulting stress regime, i.e., The wells, which are spaced 500 m horizontally, are located at 5.5 km depth and are represented as points (Fig.  2c). A permeable (k ¼ 1:0 Â 10 À13 I m À2 ) inclined fault with a dip of 75 and a thickness of 1 m is placed between the two wells. Rock permeability depends on the fracture density ω (Fig. 2b), which follows the distribution indicated in Fig. 2c. We hypothesise a radius of 100 m in which, through stimulation (not modelled here), a permeability increase of up to three orders of magnitude around the two wells was reached prior to injection.
To estimate the downhole injection temperature difference ΔT, because of data scarcity in SC reservoirs, we have extrapolated values from records on existing wells at lower temperature 62 , leading to ΔT ¼ 322 C for an initial reservoir temperature of 457 C (see Table 1). Thus, we have run 2 different cases of injection temperature (ΔT ¼ 300 C and ΔT ¼ 0°C) with constant injection and production volumetric rate q v ¼ 8:91 Â 10 À5 m 3 s À1 , in the hypothesis of an injection distributed over a 300 m section of open well, and for a total of 25 years of exploitation.

Results
Initial conditions of the reservoir. The simulated SC geothermal reservoir presents a density-driven instability in the form of two convective cells with upward flow at the centre of the model, where the heat source is located (Fig. 3a). SC water is located close to the deep magmatic source at a depth between 4 and 7 km. The maximum flow velocity, which is located within the SC region and drives the thermal plume upward, enhances the advective component of heat transport and is responsible for the pressure drop above the heat source, resulting in a sub-hydrostatic pore pressure at the centre of the model. The results are in agreement with previous computations of deep SC geothermal systems 8 . The temperature distribution along depth at the centre of the model compares well with the one of worldwide SC geothermal sites 9,17,58-61 (Fig. 3b). The geothermal systems used for comparison have differing geology, which may influence the temperature distribution at each site. Despite such site-specific features, the general trend is well captured by our model.
Injection phase. Cold water is re-injected during the exploitation phase of geothermal systems, which redirects locally the initial vertically density-driven flow to a horizontal flow from the injection to the production well (Fig. 4a). Stream lines velocity is three orders of magnitude higher than that in the natural convection system. Water accelerates along the more permeable fault, presenting an upward flow because of the heat gradient and the subsequent buoyant forces related to the lower water density in the lower part of the model. The equivalent permeability of the reservoir, accounting for the stimulated region around the wells, is computed with the Cooper-Jacob method and is k eq ¼ 1:74 Â 10 À15 I m 2 , slightly higher than initial. Pressure changes are limited to 10 MPa increase and 2 MPa decrease. The gradient is lower within a 100-m radius from the wells because of the higher permeability. The pressure distribution is symmetric between the two wells in the early stages of injection, with the neutral point located in the mid-point, i.e., in the fault. In the long-term, the pressure distribution becomes asymmetrical and the neutral point shifts towards the injection well because the injected cold water has a higher viscosity that increases the head losses around the injection well (see Fig. 4b). Water transitions from SC to liquid phase as the temperature drops (cyan line), but does not evaporate because the pressure remains above 22.064 MPa. After 25 years of exploitation, the liquid water front in the reservoir has reached the fault. While the fluid flow quickly establishes the pressure gradient between the wells, the thermal diffusion is a slower process which generates an additional time dependency in the problem.
Fault stability and induced seismicity. The injection of cold water cools down and contracts the rock around the well (Fig. 4a), causing negative volumetric strain ϵ vol ¼ tr ϵ ð Þ < 0 over an area that increases with time following thermal diffusion. This circular shape is not centred in the injection well because cold water tends  to sink as a result of its higher density. Deformation is negligible surrounding the production well, which is in isothermal conditions, hinting towards the dominant effect of thermal changes in the geomechanical response of the reservoir. After 25 years of cold water re-injection, ΔCFS < 0 (stability increase) occurs above the production well and above and below the cooled region around the injection well as a consequence of a more compressive normal stress and reduced shear stress (Fig. 5a). ΔCFS > 0 (stability decrease) occurs between the injection and production well: cooling contracts the rock and induces a normal stress reduction that drives the fault toward failure conditions. The yellow line in Fig. 5b represents the boundary of null pore pressure variation after 25 years of isothermal injection: pore pressure increases are located on the right of the yellow line (injection side), while pore pressure decrease are located to the left of the yellow line (production side). In isothermal injection, flow behaviour dominates instability, which in turn migrates downward as a consequence of the pore pressure M DP values indicate that the rock is failing in two potential ways: when M DP > 1:0, frictional-compressive failure is expected and for M DP < 0, the state of stress exceeds hydrostatic tension and tensile failure is most likely (Fig. 5c and d). Tensile failure is a consequence of the rock thermal contraction as it roughly follows the cooled down area. The inset in Fig. 5c provides a schematic representation of the stress paths that lead to tensile failure and compares the here assumed case of cohesive strength with an hypothetical cohesionless one. For isothermal injection, a much smaller portion of the reservoir exhibits increase in the mobilised shear, which remains far from failure M DP < 1:0. As expected, no tensile mobilisation is observed because the pressure increase does not drive the rock toward tensile failure, which would essentially mean hydro-fracturing the rock.
The effective stress changes caused by pressure and temperature variations alter the seismicity rate. Figure. 6 shows the rate of seismicity production in logarithmic scale at the fault and in the fractured rock for cold water injection and isothermal injection. The seismicity rate increases along the fault in both cases. However, the rate of seismic increase within the main fault is significantly enhanced when cooling occurs. Seismicity rate increase is observed in the lower part of the fault after a bit more than 5 years, but is quickly reverted as cooling effects are dominant. The reversion of seismicity from lower to higher parts of the fault is slower when no cooling occurs and only pore pressure changes act in the system. While the seismicity rate in the fault increases almost one order of magnitude for the isothermal injection, there is a considerable four orders of magnitude increase for cold water injection. It should be noted that the seismicity rate increases quickly between 5 and 10 years; as it will be shown in the next paragraph, this corresponds to a time when the mobilised strength in the fault is still not at peak and the cooling front has not reached the fault yet. The observation point in the fractured rock is closer to the injection well, which implies greater thermal stress changes and hence increased induced seismicity in earlier times compared to the further away fault. Overall, cooling has important effects on the rate of seismicity production.
Fault mobilisation is assumed whenever μ fr ! 0:6 63 and its evolution with time is shown in Fig. 7a for both scenarios. In isothermal conditions, the mobilised friction in the fault never exceeds critical values and it can be considered stable throughout the whole simulation. On the contrary, for cold water injection, the fault undergoes progressively larger mobilisation size. The size of the unstable patch L r increases with time as the thermal front approaches the fault. Because rock contraction acts as a long-range mechanism to destabilise the fault, the temperature change in the fault alone would not provide sufficient insight. To account for this effect, we computed the point at which the temperature decrease halfway, i.e., the point at which T ¼ 300°C. The point of average temperature migrates with time during thermal diffusion and its distance from the fault along the line connecting the two wells is D c (Fig. 7a). We have plot D c vs. L r and found that a linear relation exists between this two values, as the length of the stress redistribution induced by cooling is proportional to the radius of the cooled region (Fig. 7b). The size of the mobilised fault patch L r is instead proportional to the logarithm of time / log 10 t ð Þ, which is a consequence of the thermal diffusion process (Fig. 7c).

Discussion
SC geothermal systems are located near thermal instabilities (heat plumes) and convective cells establish upward fluid migration and a sub-hydrostatic pressure in the deeper part of the reservoir as a consequence of the reduced fluid density. Estimates of the temperature and pressure distribution are of great help for preliminary assessment of geothermal potential and for project design. Thus, through numerical modelling, concurring hypotheses can be tested with the aim of better locating potential SC resources. Proper assessment of resources and initial conditions is fundamental to design geothermal exploitation systems, from conception to decommissioning. Numerical modelling can be applied to specific geothermal sites, accounting for local geological structures and physical parameters of the rocks.
During cold water re-injection, two concomitant timedependent processes are active: a fast (/ 10 6 s) diffusion process that counteracts the pore pressure gradient generated by the forced flow, and a slow (/ 10 8 s) advection process that transports heat in the fractured rock. The mechanical behaviour of the reservoir is influenced by these two processes, the short-term time dependence is governed by fluid flow and the long-term one by thermal effects. Viscosity increases during cold water re-injection, which in turn increases the total water head and shifts the nullpressure point towards the injection well (Fig. 4). Sinking effects are evident and the liquid water front extends asymmetrically around the injection well. The liquid front reaches the fault after roughly 25 years of injection, implying that SC fluid will still be produced. We have assumed a circular shape of the stimulated area around the wells because stimulation with supercritical water forms a cloud of increased permeability instead of a single planar fracture 27 . Nevertheless, hydrofracturing and hydroshearing with supercritical water is still a poorly understood topic. Different permeability structures around the wells will affect timing and orientation of the thermal front through advective processes, but the general flow conditions will not change in comparison to the current assumptions. Cooling increases the failure potential and the number of recorded seismic events in both the fractured rock and the main fault. While instabilities in the fractured rock are expected to produce microseismicity, as the existing fractures have limited extent and hence likewise limited rupture area, greater events are associated with reactivation of the main fault that will start slipping after roughly 10 years. Tensile failure is a consequence of thermal contraction and might have a beneficial effect by increasing the permeability, but can also greatly enhance the rate of microseismicity by reducing the mean normal stress acting on the fractures. What is likely to happen is an increase of frequency in small earthquakes in connection with cold water injection during early time. Larger events are delayed and their timing depends essentially on the distance between the cooling front and the existing major fault/s, with important consequences in terms of reservoir management and lifetime assessment. This pattern has been observed at The Geysers, where the number of seismic event closely correlates with injection history 64 .
The magnitude of the seismic events depends on the critical nucleation length in quasi-static regime 65 , which is controlled by frictional weakening associated with slip and dynamic effects of melting and pressurisation 66,67 . For injection-induced seismicity, nucleation size and dynamic fault rupture arrest strongly depend on pore pressure increment in the fault and on the acting background stress state 35 . It is therefore difficult to compute the magnitude of the larger events that involve the slip into the main fault, due to the quasi-static approach adopted here. Moment magnitude estimates can be obtained by enhancing the modelling approach and including discontinuous fault displacement in the finite element formulation with frictional weakening, and eventually coupling the analyses with a dynamic solver to estimate runaway ruptures and propagating crack fronts in dynamic conditions 68,69  The consequences must be evaluated on individual cases and holistic risk assessment approaches will provide estimates of potential damage, infrastructure disruption and disturbance to the local population. State-of-the-art numerical analyses are a fundamental tool to address the complexities involved in ESGS: entangled with monitoring techniques, they provide precious insights into the physical behaviour of ESGS and will prove essential for long-term safe reservoir management.

Methods
Porosity evolution. We report here the detailed derivation of mass, energy and momentum balance equations of porous media leading to the complete description of thermo-hydro-mechanical (THM) coupled processes presented in Eq. (1). Subscripts w stand for water, s for solid and m for porous medium. The sign convention of stresses and strains is that of solid mechanics, i.e., positive is for tensile. The equations are formulated as a function of the primary variables of pore pressure, p, temperature, T and displacement vector, u. The mass balance of the solid skeleton writes where n is porosity, ρ s is solid density, v s is the solid's velocity, t is time and the operator d s =dt represents the material time derivative of the solid defined as Expanding Eq. (3) yields an expression for the rate of change in porosity where _ ϵ v is the volumetric strain rate. The solid density can be expressed as a function of temperature T, mean effective stress σ 0 m ¼ tr σ 0 ð Þ=3 and pore pressure p 71 such that applying chain derivation yields with linear thermal expansion coefficient and compressibility of the solid skeleton to pore pressure changes where K s is the intrinsic bulk modulus of the solid phase. The compressibility of the solid phase in response to effective stress changes writes 1 ρ s @ρ s @σ 0 so that the solid density constitutive equation becomes The effective stress acting between the grains due to volumetric strains is reduced by grain volume changes due to thermal effects and due to compression in response to fluid pressure where K is the drained bulk modulus and Biot's coefficient is defined as Eq. (6) can now be recombined into the constitutive equation of the solid grain density, yielding 71 1 Introducing Eq. (13) into Eq. (5), the final expression of porosity variation becomes Mass balance equation. Recalling Eq. (5), the mass balance of the solid skeleton writes The mass balance of the fluid flow writes where v w is the fluid velocity and Q H is the source term. Recalling the definition of material derivative in Eq. (4), we can write the transformation between material derivatives with respect to the fluid phase d w =dt to material derivatives with respect to the solid phase d s =dt, i.e., for a given quantity a The fluid balance equation can then be rewritten as Substituting the expression of porosity material derivative from Eq. (18) into solid mass balance Eq. (15), with the hypothesis of small gradients and re-arranging, yields The density of water depends on pressure p and temperature T, so that its constitutive equation writes where β w is water compressibility and α w is water volumetric thermal expansion coefficient. The seepage velocity is v ¼ n v w À v s ð Þand relates to pressure changes through Darcy's law where k is the intrinsic permeability tensor, μ w is the water dynamic viscosity and g is the body force vector (only gravity is present in this case). Recalling that volumetric strain rate is _ ϵ v ¼ ∇ Á v s and the expression of solid compressibility Eq. (13), the continuity equation of the porous medium Eq. (19) can be rearranged into The first coefficient of Eq. (22) represents the storage term S, while the second coefficient describes the average thermal expansion of the porous medium. Together with the volumetric strain term, fluid-solid differential thermal expansion effects are captured.
Energy balance equation. Advective and conductive heat transfer in the porous medium are derived from the heat balance equation where c is specific heat and λ is the thermal conductivity, which for porous media write and where I is the identity matrix. The first term of Eq. (23) represents the heat storage term, the second term is related to conductive heat transport via Fourier's law and the third term is the advective heat transport term. The right-hand side Q T is a thermal source term. Thermal effects due to deformation have been neglected.
Linear momentum balance equation. The linear momentum balance of the porous medium writes The first term is the divergence of Cauchy's second order total stress tensor σ and the second term represents body forces due to gravity. The effective medium density is expressed as and the stress acting exclusively on the solid skeleton is Biot's effective stress σ 0 which writes The constitutive equation of the thermo-elastic solid writes where K and G are the bulk and shear moduli of elasticity, respectively, the infinitesimal strain tensor ϵ is the symmetric part of the displacement gradient and ϵ D is the deviator of the strain tensor. Other elastic constants can be employed, such as Young's modulus, E, and Poisson's ratio, ν.
Fluid compressibility and thermal expansion. We have applied a finite difference scheme to compute the water's compressibility and thermal expansion Permeability distribution. The fracture density around the two wells follows a spline bell-shape distribution where R p ¼ 100 m and ω 1 is taken as a Weibull distribution with ω 0 ¼ 0:25 andk ¼ 10.
Numerical implementation. The spatial discretization is achieved with the Galerkin method and the temporal discretization is carried out with an implicit finite difference scheme 72,73 . The weak forms of the balance equations write, for the mass Z Ω S @p @t ψdΩ þ where q T ¼ À λ m ∇T þ ρ w c w vT groups together the advective and diffusive heat fluxes of the heat equation 73 . Linear test functions ψ 2 V 1 & H 1 Γ Ω ð Þ 1 are employed for the heat and mass transport equations (TH) and ψ 2 V n & H 1 Γ Ω ð Þ n is a quadratic test function in n-dimensional space applied to the weak form of the momentum balance equation, Ω denotes the whole problem domain and Γ its boundary. The primary variables are approximated by admissible interpolation functions belonging to the Taylor-Hood space of finite elements. The system is re-written as where M and K are mass and Laplace matrices, respectively, functions f denote the right-hand side and the indices indicate the process to which they are specific. The mechanical stiffness matrix writes K M ¼ R Ω B T DBdΩ, with the B containing the gradient of the shape functions N to calculate the strains ϵ. The matrix D is the constitutive tangent matrix and the shape functions approximate the sought solution u in terms of its nodal valuesû as The general form of the mass and Laplace matrices writes and whereM is a scalar material parameter,K is a material tensor and linear or higher order shape function are used based on the process. The solution is achieved through a staggered scheme that combines Newton-Raphson and Picard solution methods 72 .

Data availability
Data generated in this study can be requested to the corresponding author. The numerical code employed for the analyses is open source and can be freely downloaded at https://doi.org/10.5281/zenodo.3294169. All data supporting the findings of this study are available within the paper and the methods section. Further information on compiling the code, as well as tutorials and simulation examples about OpenGeoSys are freely accessible at https://www.opengeosys.org/.