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.

Introduction

Enhanced 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 injection-induced 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-8]: as a consequence of high temperature and depth, fluids can be found in supercritical (SC) conditions, which for pure water occurs at temperature and pressure  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 build-up. 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-18] and Los Humeros, Mexico[19-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-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 hydro-fracturing 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-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-43]. During stimulation injection, micro-seismicity is associated with hydro-fracture propagation and larger events happen during prolonged injection[44]. Nevertheless, there are reported cases in which 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-of-the-art numerical analyses to understand the main physical mechanisms controlling their geomechanical behaviour.

Fig. 1

Enhanced Supercritical Geothermal System. Schematic representation of a doublet system in an Enhanced Supercritical Geothermal System (ESGS), with pre-existing magmatic heat source that generates a convective hydrothermal circulation. A doublet system of injection/production alters the pressure and temperature fields, leading to potential fault instability

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 momentum balance of porous media is[49]and is formulated in terms of primary variables describing the fields of pore pressure , temperature and displacement . is time and is the permeability tensor. Subscripts ‘s’ and ‘w’ stand for solid and for water, respectively, and is the identity tensor. The fractured rock properties are set to values typical of a crystalline rock with porosity , solid density , linear thermal expansion coefficient , specific heat capacity , thermal conductivity , Young’s modulus  GPa, Poisson’s ratio and Biot’s coefficient , where K is the drained bulk modulus and is the intrinsic bulk modulus of the solid phase. The gravity acceleration vector is indicated by and the source terms for the thermal and hydraulic problem are and , respectively. The fault has different mechanical properties with ,  GPa and . 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 , dynamic viscosity , specific heat capacity , thermal conductivity , compressibility and thermal expansion are functions of its thermodynamic state (see Methods for the computation of fluid compressibility and thermal expansion).

Fig. 2

Numerical model and equivalent permeability. a Model setup of the complete reservoir for initialisation phase, in which a constant bottom temperature is applied to initialise the thermal plume. b Permeability curves of volcanic rocks are function of porosity[50] in a logarithmic relation with representing the degree of fracturing of the rock. c Detail of the geothermal doublet located at 5.5 km depth with the fault between the wells and the distribution of initial permeability for the injection phase

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 , 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 and intact rock to porosity , based on permeability measurements on 70 different samples of volcanic and magmatic rocks[50]. We have introduced a formulation that interpolates between intact and fractured rock permeability as , in which is a newly defined parameter representing either fully fractured () or intact () rock (Fig. 2b).

The Coulomb Failure Stress is suited to estimate potential failure[34], where and 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 ) and a frictional coefficient , which corresponds to a friction angle . 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 . The Drucker–Prager criterion is an alternative to estimate the failure potential in terms of invariants of the stress tensor (direction independence)where  MPa is the cohesion of the fractured rock, is the Drucker–Prager deviatoric stress, is the mean effective stress, is the effective stress tensor. The mobilised failure ratio for the Drucker–Prager model is , where is the acting deviatoric stress and .

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, , writes , where is a characteristic time, is a constitutive parameter and equals for hydrothermal systems[55], is the background shear stressing rate, and the time derivative symbol is .

Numerical model setup

The numerical solution is obtained via the finite element method (FEM) with a Lagrangian approach with respect to the solid. The model is implemented in the C++, object-oriented and open-source finite element solver OpenGeoSys[57]. We have built a representative FEM model of a typical deep SC geothermal system with a two-dimensional plane strain numerical model 7 km deep and 10 km wide (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 over a width of 2 km. The hydraulic conditions are no-flow 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 . The initialisation phase is performed in hydro-thermal 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-61].

The injection phase, in which the coupled THM system of Eq. (1) is solved, uses the and distribution obtained at the end of the initialisation phase. The initial stress state is given by a lithostatic vertical stress corresponding to a solid density of and a normal faulting stress regime, i.e., , with . The wells, which are spaced 500 m horizontally, are located at 5.5 km depth and are represented as points (Fig. 2c). A permeable () inclined fault with a dip of 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 , because of data scarcity in SC reservoirs, we have extrapolated values from records on existing wells at lower temperature[62], leading to for an initial reservoir temperature of (see Table 1). Thus, we have run 2 different cases of injection temperature ( and °C) with constant injection and production volumetric rate , in the hypothesis of an injection distributed over a 300 m section of open well, and for a total of 25 years of exploitation.

Table 1

Estimate of the injection temperature, obtained by extrapolation of values compiled in a literature review of several geothermal wells[62]. The linear extrapolation from literature values (with the estimated average temperature of Low Enthalpy systems of ) leads to the value of for the SC case

Category	T_min	T_max	<T>	ΔT
	^∘C	^∘C	^∘C	^∘C
Hot water	−	220	140	55
Low enthalpy	220	250	235	131
Medium enthalpy	250	300	275	186
High enthalpy	250	330	290	169
Supercritical	−	−	457	322

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.

Fig. 3

Initial conditions in the reservoir. a The solution of the hydro-thermal (TH) problem refers to a time of roughly 31’000 yr. Instability triggers two convective cells and supercritical resources are found in the deeper area and are delimited with a cyan line. b The temperature profile is compared with several geothermal wells in the world close to or at supercritical conditions, such as Reykjanes (RE)[59], Los Humeros (LH)[58], Larderello (LA)[61], Mofete (MO)[60], The Geysers (TG)[17] and IDDP-2 at Reykjanes[9]

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

Fig. 4

Pressure distribution between the wells and deformation. a The streamlines go from the injection well to the extraction well and the volumetric strain follows the thermal field after 25 years of injection. The cyan line indicates the liquid front, which is close to the fault. b The pressure distribution is asymmetrical in the long-term and the mid-point moves with time towards the injection well because of the higher viscosity of the liquid water around the injection well due to thermal quenching

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 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, (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). (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 increment and normal effective stress reduction. The lower magnitude of CFS in the isothermal injection implies that pore pressure effects play a smaller role than thermal effects in terms of fault stability.

Fig. 5

Fault reactivation and reservoir stability. Coulomb Failure Stress changes () after 25 years for the cold water (a) and isothermal injection (b) cases and Drucker-Prager mobilised failure ratio after 25 years for the cold water (c) and isothermal injection (d) cases. Negative values of imply : hydrostatic tension failure will be reached and tensile fractures will be likely to appear, as highlighted by the schematic representation in the inset in c

values indicate that the rock is failing in two potential ways: when , frictional-compressive failure is expected and for , 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 . 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.

Fig. 6

Rate of seismic production. The rate of seismicity production is computed at different locations in the main fault and in the fractured rock (curve colours in a refers to points which are indicated in the model in b), for isothermal and cold water re-injection. Injecting cold fluids generates significantly higher induced seismicity than injecting fluids in thermal equilibrium with the reservoir

Fault mobilisation is assumed whenever [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 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 °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 (Fig. 7a). We have plot vs. 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 is instead proportional to the logarithm of time , which is a consequence of the thermal diffusion process (Fig. 7c).

Fig. 7

Fault mobilisation. a Evolution with time of mobilised friction along the fault coordinate for cold water and isothermal injection. b, c The fault mobilised length () of the cold water injection has a linear dependence on the distance of the mean temperature () front (b) and a logarithmic dependence on time (c)

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 time-dependent processes are active: a fast ( s) diffusion process that counteracts the pore pressure gradient generated by the forced flow, and a slow ( 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 null-pressure 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]. Nevertheless, we have shown how the mobilised size of the fault increases logarithmically with time as the cooling front approaches it. The size of the fault mobilisation is quite large and, as an example, after 25 years of injection has reached  m in length and, considering the  m of open well section for injection, could have a surface in the order of . Thus, large events triggered by cooling of the reservoir are possible and their timing is controlled by the thermal diffusion process.

Our approach is a first attempt to model the complex geomechanical behaviour of ESGS. Several open questions remain at present to be addressed, and future studies should focus on the influence of inelastic behaviour (rock strength decays at high temperature[70]), on the role of geological structures and multiple faults, on the full three dimensional reservoir behaviour, on the presence of fluid mixtures (), on modelling the frictional behaviour of faults and dynamic rupture. Nevertheless, we have provided an important piece of evidence to support the role of cooling in connection to induced seismicity and fault destabilization. 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 stand for water, for solid and 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, , temperature, and displacement vector, . The mass balance of the solid skeleton writeswhere is porosity, is solid density, is the solid’s velocity, is time and the operator represents the material time derivative of the solid defined asExpanding Eq. (3) yields an expression for the rate of change in porositywhere is the volumetric strain rate. The solid density can be expressed as a function of temperature , mean effective stress and pore pressure [71] such that applying chain derivation yieldswith linear thermal expansion coefficientand compressibility of the solid skeleton to pore pressure changeswhere is the intrinsic bulk modulus of the solid phase. The compressibility of the solid phase in response to effective stress changes writesso that the solid density constitutive equation becomesThe 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 pressurewhere K is the drained bulk modulus and Biot’s coefficient is defined asEq. (6) can now be recombined into the constitutive equation of the solid grain density, yielding[71]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 writesThe mass balance of the fluid flow writeswhere is the fluid velocity and 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 to material derivatives with respect to the solid phase , i.e., for a given quantity The fluid balance equation can then be rewritten asSubstituting the expression of porosity material derivative from Eq. (18) into solid mass balance Eq. (15), with the hypothesis of small gradients and re-arranging, yieldsThe density of water depends on pressure and temperature , so that its constitutive equation writeswhere is water compressibility and is water volumetric thermal expansion coefficient. The seepage velocity is and relates to pressure changes through Darcy’s lawwhere is the intrinsic permeability tensor, is the water dynamic viscosity and is the body force vector (only gravity is present in this case). Recalling that volumetric strain rate is and the expression of solid compressibility Eq. (13), the continuity equation of the porous medium Eq. (19) can be rearranged intoThe first coefficient of Eq. (22) represents the storage term , 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 equationwhere is specific heat and is the thermal conductivity, which for porous media writeandwhere 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 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 writesThe 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 asand the stress acting exclusively on the solid skeleton is Biot’s effective stress which writesThe constitutive equation of the thermo-elastic solid writeswhere K and are the bulk and shear moduli of elasticity, respectively, the infinitesimal strain tensor is the symmetric part of the displacement gradientand is the deviator of the strain tensor. Other elastic constants can be employed, such as Young’s modulus, , and Poisson’s ratio, .

Fluid compressibility and thermal expansion

We have applied a finite difference scheme to compute the water's compressibilityand thermal expansion

Permeability distribution

The fracture density around the two wells follows a spline bell-shape distributionwhere  m and is taken as a Weibull distributionwith and .

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

heatand momentumwhere groups together the advective and diffusive heat fluxes of the heat equation[73]. Linear test functions are employed for the heat and mass transport equations (TH) and is a quadratic test function in -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 and are mass and Laplace matrices, respectively, functions denote the right-hand side and the indices indicate the process to which they are specific. The mechanical stiffness matrix writes , with the containing the gradient of the shape functions to calculate the strains . The matrix is the constitutive tangent matrix and the shape functions approximate the sought solution in terms of its nodal values asThe general form of the mass and Laplace matrices writesandwhere is a scalar material parameter, 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].