3D mechanical analysis of geothermal reservoir operations in faulted sedimentary aquifers using MACRIS

Accurate and efficient predictions of three‑dimensional subsurface stress changes are required for the assessment of geothermal operations with respect to fault stability and the potential risk for induced seismicity. This work extends the model capabili‑ ties of M echanical A nalysis of C omplex R eservoirs for I nduced S eismicity (MACRIS) to account for high‑resolution thermo‑elastic stress evaluations in structurally complex (i.e. faulted) and matrix permeability dominated geothermal systems. By adopting a mesh‑free approach suitable to industry standard flow simulation models, MACRIS is capable of preserving the complex 3D hydraulic development of the injected cold‑water volume and the 3D geometrical complexities of the reservoir model. The workflow has been applied to three‑dimensional models with clastic reservoir charac‑ teristics representative for low enthalpy geothermal exploitation in the Netherlands. The models are marked by a single fault, subject to no and normal offset. Comparison of simulated stress evolutions in MACRIS with alternative analytical solutions highlight the effects of stress arching involved in the poro‑ and thermo‑elastic stress devel‑ opments on complex faults intersected by or in direct contact with the cold‑water volume. Results are in agreement with previous studies and show the effect of thermal stressing to be dominant, arching of stresses to occur at the rim of the cold‑water vol‑ ume, and in cooling reservoirs, the intersection area of the cold‑water volume in direct contact with the fault plane to be the main driver for fault reactivation and subse‑ quent seismic potential. Moreover, results show the effects of stress arching (i) to be enhanced in the case of reservoir throw and flow compartmentalization, and (ii) to be reduced by a relative increase in conductive heat transfer between the reservoir and surrounding formations.


Introduction
Geothermal energy resources have great potential for district-and greenhouse heating, and electricity generation.The demand for geothermal power as a renewable and sustainable form of energy has therefore seen a significant growth over the past decades (Lund and Boyd 2016;Lund and Toth 2020).However, a limiting factor in the development of geothermal energy resources is the occurrence of induced seismic events which have been observed in numerous geothermal projects around the world.Operations in several geothermal projects have been suspended as a direct result of induced events with a magnitude that exceeds initial predictions (Häring et al. 2008;Evans et al. 2012;Zang et al. 2014;Buijze et al. 2019;Muntendam-Bos et al. 2021;Kinscher et al. 2023).This demonstrates the importance of improved seismic hazard assessment and enhancing the robustness of associated modelling approaches in the further development of geothermal power as a sustainable energy resource.
Induced seismicity is often associated with the development of Enhanced Geothermal Systems (or EGS) characterized by the short-term hydraulic stimulation of fractured and competent rock types to enhance reservoir permeability (Moeck 2014; Buijze et al. 2019).Over the past decades, induced seismicity in the context of the development of EGS has been extensively studied (Majer et al. 2007;Jeanne et al. 2014Jeanne et al. , 2017;;Wassing et al. 2014;Gaucher et al. 2015;Tomac and Sauter 2018;Zang et al. 2014Zang et al. , 2019;;Rathnaweera et al. 2020;Zhu et al. 2023).Amongst the most important findings is the correlation between the maximum seismic magnitude and the total volume of injected fluid for hydraulic stimulation, i.e. the maximum magnitude increases with time for continuous fluid injection (Shapiro et al. 2010;McGarr 2014).Although induced micro-seismic events ( M L < 1.5 ) are a known side effect of hydraulic stimulation (Zoback 2007), Zang et al. (2014) state that short-term stimulation in EGS has a higher tendency to produce larger induced events compared to equivalent oil and gas operations.This is largely related to rupture of nearby pre-existing fault zones which can generate isolated events of considerably larger magnitude (Zang et al. 2019), as observed in EGS projects in Soultz-sous-Forêts ( M L 2.7 ) (Baisch et al. 2010), Basel ( M L 3.4 ) (Häring et al. 2008) and Pohang ( M L 5.4 ) (Kim et al. 2018).As such, more recent studies have contradicted the assumption of induced ruptures to be confined to the injected volume and have proposed that the maximum seismic magnitude adheres to the limit of tectonic seismicity when induced events in the development of EGS occur on nearby pre-existing fault zones (Izadi and Elsworth 2014;Atkinson et al. 2016;Van der Elst et al. 2016).
Felt seismicity ( M L > 1.5 ) has also been reported outside EGS (Buijze et al. 2019).In particular in the operational phase, the injection of cold fluids creates a temperature contrast capable of generating significant thermal stresses, which, superimposed on the pressure induced stress changes, can lead to an increase in seismic hazard (Candela et al. 2018;Buijze et al. 2019;Parisio et al. 2019;Wassing et al. 2021;Kivi et al. 2022).In hydrothermal systems, relying predominantly on faults and fracture permeability for the long-term operational circulation of fluids (Moeck 2014), the Geysers high enthalpy geothermal field in California (Majer and Peterson 2007) is a prime example marked by an increase in induced seismicity ( M L > 4 ) as a result of increased fluid injection to support pressure maintenance in the reservoir.Geothermal projects targeting faulted and fractured limestones in Germany ( M L 2.1−3.5 ) (Seithel et al. 2019), the Nether- lands ( M L 1.7 ) (Muntendam-Bos et al. 2021) and Belgium ( M L 2.2 ) (Kinscher et al. 2023) have produced felt seismic events at different stages during fluid circulation, i.e. up to multiple years after the onset of cold-water injection.Candela et al. (2018) show that the injection of cold fluids in fractured geothermal systems can impose significant thermal stresses on pre-existing fault zones in the vicinity of the geothermal doublet, independent of injected volume.They argue that the destabilization of these fault zones as a result of long-term cooling of the reservoir rocks may form the primary control on induced seismicity in geothermal systems.A recent study by Wassing et al. (2021) has highlighted the spatial and temporal pattern of seismicity in a fractured and faulted carbonate reservoir.They emphasize the contribution of fracture spacing and the dominant effect of thermal cooling on stress alteration and subsequent seismic potential of a fault located close to an injection well.
Low enthalpy geothermal systems that are characterized by sedimentary reservoirs with high matrix permeability (i.e.mostly clastic reservoirs) are up till now not associated with felt seismicity (Moeck 2014; Buijze et al. 2019).For a clastic reservoir, Kivi et al. (2022) have investigated the effects of long-term cooling on the reactivation potential of faults located at long distances (~ 1 km) from the geothermal doublet.They show cooling-induced stress changes in the vicinity of the geothermal doublet to extend far away from the cold-water volume, especially in the long term, and increase the reactivation potential of distant faults.The previous studies highlight the control of thermal cooling and its spatial extent on potential stress instabilities of faults for geothermal reservoir operations.For sedimentary reservoir settings, enhanced geomechanical analysis techniques for seismic hazard assessment commonly include the development of a finite element model (FEM).A geomechanical FEM requires building a dedicated mesh with a very dense resolution to be able to capture sharp stress variations caused by reservoir compaction or dilation (Orlic and Wassing 2013).For structurally complex (i.e.faulted) reservoirs, building geomechanical FEM can be a major challenge as the required gridding system differs considerably from that applied in industry standard simulation models.Reservoir characterization and flow models are commonly constructed using a so-called corner-point grid geometry and thus require conversion to the tetrahedral gridding system typically applied in geomechanical FEM (Koutsabeloulis and Zhang 2009;Cappa and Rutqvist 2011;Orlic and Wassing 2013;Sanz et al. 2015;Lele et al. 2016).
Van Wees et al. (2019) have presented a novel mesh-free geomechanical modelling approach, known as Mechanical Analysis of Complex Reservoir for Induced Seismicity (MACRIS), capable of calculating induced poro-elastic stress changes and associated seismic moment response in structurally complex (i.e.faulted) hydrocarbon reservoirs.By adopting a mesh-free approach suitable to industry standard 3D reservoir characterization and flow simulation models, MACRIS specifically preserves (i) the complex 3D development of the pressure field and (ii) the 3D geometrical complexities of the reservoir, including reservoir throw and fault geometry, following high quality subsurface seismic data (Moeck 2014;Mijnlief 2020;Van Wees et al. 2018, 2020;Candela et al. 2022).Implementation of the Barnes-Hut algorithm (Barnes and Hut 1986) increases the spatial resolution of stress and reduces the computational intensity.MACRIS has been used successfully to analyse induced seismicity in the heavily faulted Groningen Gas Field in the Netherlands (Candela et al. 2018(Candela et al. , 2022)).
In this work, the existing model capabilities of MACRIS are extended to account for high-resolution thermo-elastic stress evaluations in structurally complex geothermal reservoirs.The novel implementation of thermo-elasticity in MACRIS is described in detail and benchmarked through comparison with a finite element (FEM) solution.
MACRIS is subsequently applied to models with clastic reservoir characteristics representative for low enthalpy geothermal exploitation in the Netherlands (Buijze et al. 2023).The models include over-and underburden rock and are marked by a single fault, subject to no and normal offset, which is located such that the hydraulically propagating cold-water volume originating from an injector well crosses the fault plane.Comparison of MACRIS with the analytical solution for uniaxial reservoir compaction enables a detailed analysis of the effects of stress arching involved in the poro-and thermo-elastic stress development in complex geothermal reservoirs.Results are in agreement with previous studies that show (i) thermal stressing to be the dominant mechanism in the reactivation of faults in the vicinity of a geothermal doublet (Candela et al. 2018;Wassing et al. 2021;Kivi et al. 2022); (ii) arching of stresses to occur at the thermal front of the cold-water volume (Wassing et al. 2021); and (iii) in cooling reservoirs, the intersection area of the injected cold-water volume in direct contact with the fault plane to be the main driver for fault reactivation and subsequent seismic hazard assessment (Wassing et al. 2021).Moreover, results show the effects of stress arching (i) to be enhanced in the case of reservoir throw and flow compartmentalization, and (ii) to be reduced by a relative increase in conductive heat transfer between the reservoir and surrounding formations.

Model geometry and parametrization
In the Netherlands, geothermal projects have mainly been developed in clastic reservoirs with high primary porosity and permeability, targeting the same formations in which the large hydrocarbon fields are located (Mijnlieff 2020;Buijze et al. 2023).A wealth of existing oil and gas data, including thousands of 2-to 6-km deep exploration and production wells, 2D and 3D seismic lines and hundreds of thousands of core plug measurements, has led to a well-constrained database of the Dutch subsurface (Van Wees et al. 2012, 2017).The occurrence of felt seismicity associated with the Groningen gas field (Muntendam-Bos et al. 2021) has led to an increase in the research effort on induced seismicity in the Netherlands (Van Wees et al. 2018;Candela et al. 2019Candela et al. , 2022;;Wassing et al. 2021;Buijze et al. 2023).The configuration of the three-dimensional geological model used in this work (Fig. 1) is characteristic of the Dutch subsurface and based on the development of geothermal doublet systems in the Netherlands (Willems et al. 2017;Buijze et al. 2019;Mijnlieff 2020;Van Wees et al. 2020).The model has a simplified geometry and consists of a highly permeable reservoir formation over-and underlain by impermeable seal and base formations, respectively.Each of the model formations are assumed homogeneous and isotropic in regard to both dynamic and mechanical properties.Linear elastic material behaviour is assumed with no contrast in elastic properties between the reservoir and surrounding formations.The lateral extent of the model is 5 km in both directions.A limited depth interval of 400 m is considered, with a reservoir thickness of 100 m.The grid cell dimensions of the discretized model are 50 m × 50 m × 10 m in the x , y and z directions, respectively.The geological model is placed in an extensional tectonic setting assuming a normal faulting regime.A single fault is considered and defined by the pillars of the model grid corresponding to the fault plane.Two geological scenarios are considered regarding the structural complexity of the model; (1) no offset along the fault plane, and (2) normal fault offset of half the reservoir thickness.In both geological scenarios, the mid-reservoir depth is 2.3 km and defined as halfway the top of the footwall and base of the hanging wall.The minimum horizontal stress is oriented perpendicular to the fault strike to ensure a favourable slip tendency of the fault.The model is operated by a geothermal doublet which is oriented perpendicular to fault strike and located such that doublet flow crosses the fault plane.Both wells are sub-vertical following the model grid and placed at equal distance from the fault.The wells have an 8.5-in.diameter and open-hole sections that are limited to the reservoir interval.Fault transmissibility is treated equal to the reservoir transmissibility, regardless of the offset of the fault, and is governed by the interconnected reservoir interval over the fault plane.Adopting an explicit coupling approach, the three-dimensional developments of pressure and temperature in the reservoir and over-and underburden formations are obtained from the Open Porous Media (OPM) Flow reservoir simulator (Rasmussen et al. 2021).Initially, the reservoir is assumed to be at hydrostatic pressure.A geothermal gradient of 31 °C/km is assumed with an average surface temperature of 10 °C.Cold water is injected into the reservoir at a temperature of 30 °C for a period of 50 years.The arrival of the cold-water volume (or thermal front) at an arbitrary position within the model is defined by a noticeable decrease in temperature (~ 1 °C) at the respective position.The arrival of the thermal front at the production well is often referred to as thermal breakthrough.The injection flowrate is chosen such that after 50 years thermal breakthrough occurs in the production well (Van Wees et al. 2012).Mechanical, thermal and hydrological model parameters and initial conditions are based on those relevant to geothermal doublet operations in the Netherlands (Buijze et al. 2023), and are summarized in Table 1.
Induced seismicity is a site-specific phenomenon depending on key geological and operational factors that need to be considered in seismic hazard assessment (Buijze et al. 2019).The model presented in Fig. 1 poses as a simplified representation of the subsurface.In reality, the sandstone reservoir formation is often intercalated with clay and shale layers resulting in hydraulically isolated reservoir zones of considerably lower thickness (Buijze et al. 2019;Mijnlieff 2020).Furthermore, many different doublet configurations are possible in regard to the distance and orientation relative to a fault.Such geological and operational variations directly affect the three-dimensional development of the subsurface pressure and temperature fields, and, in turn, affect the development of thermal stresses on a fault near the geothermal doublet.To this end, additional model scenarios are considered in which (1) the reservoir thickness is decreased; (2) the reservoir formation is intercalated with low permeability layers; (3) the injection well is located closer to the fault; and (4) the doublet is oriented parallel to a sealing fault.

Linear thermo-elasticity
The stress-strain relation for linear poro-and thermo-elasticity in a homogeneous and isotropic medium is given by (Fjaer et al. 2008): where σ ij , ε vol , ε ij and δ ij are the total stress tensor, volumetric strain (which is defined as the sum of the principal strains), strain tensor and Kronecker delta, respectively.Equation 1 is written in terms of the Young's modulus E and Poisson ratio ν , two elas- tic parameters describing the compaction of a material.The effect of temperature is accounted for in the final term on the right-hand side, with α ij being the linear thermal expansion coefficient and T describing the difference between initial and predicted reservoir temperatures in the simulation.The effect of pressure is captured by the effective stress σ ′ , which is defined as: in which β is the Biot coefficient and P the pressure differential in the reservoir.In order to obtain the maximum stress response the value of the Biot coefficient is equal to 1.Note that without definition of the effective stress Eq. 1 only describes the thermo-elastic stress response.The term holding the pressure differential in Eq. 2 is interchangeable (1) with the final term in Eq. 1 to obtain a relation describing solely the poro-elastic stress response (Fjaer et al. 2008).
Well-known geomechanical expressions of oil and gas production are reservoir compaction and associated surface subsidence (i.e.Mulders 2003;Fokker and Orlic 2006;Fjaer et al. 2008).In the case the reservoir is laterally uniform and not offset by the fault, the subsurface stress response can be approximated using an analytical approach based on uniaxial reservoir compaction.In this approach, two main assumptions are made; (1) the lateral strain is neglected, i.e. ε h = ε H = 0 , and (2) the total vertical stress is con- stant, i.e. �σ v = 0 .The latter results in any stress arching effects (i.e.shielding of stress by the overburden) being neglected (Fjaer et al. 2008).With the change in total vertical stress equal to zero, the vertical elastic strain follows from Eq. 1, and becomes for the effect of temperature and pressure, respectively.The effective vertical stress change follows directly from Eq. 2, i.e. �σ ′ v = −�P , and substitution of the vertical strain components in Eq. 1 yields the effective horizontal stress changes in both lateral directions as: which conform to Hooke's law and poro-elasticity (Fjaer et al. 2008), are non-zero despite the previously neglected lateral strains.The assumption of constant vertical stress in the uniaxial stress solution means that the total vertical stress does not change laterally, nor due to changing pressure or temperature (Fjaer et al. 2008;Van Wees et al. 2014).For reservoirs with high aspect ratios (i.e.lateral extent >> thickness) or charac- terized by lateral contrasts in pressure and temperature changes (i.e. at the rim of the injected cold-water volume or at offsetting faults), the effects of stress arching become increasingly important (Mulders 2003;Fokker and Orlic 2006;Van Wees et al. 2014).In such reservoirs, the uniaxial solution presented above is therefore likely to be inadequate in the approximation of the subsurface stress response and a Thermo-Hydro-Mechanical (THM) reservoir workflow beyond the uniaxial approach is required.The importance of stress arching effects in geothermal reservoirs and its implications for fault reactivation have been demonstrated by Wassing et al. (2021).Alternatively, inclusion theory and the inflation point source solution (i.e. the nucleus of strain concept) allow for evaluation of the full 3D subsurface stress field including the effects of stress arching (Van Wees et al. 2019;Jansen et al. 2019).These methods have been used in many subsidence studies for reservoir compaction over the past decades (Eshelby 1957;Geertsma 1973;Segall and Fitzgerald 1998;Fjaer et al. 2008;Fokker and Orlic 2006).
MACRIS employs a semi-analytical and mesh-free approach, based on the aforementioned inflation point source solution, in which the need to build a dedicated FEM model for the geomechanical analysis is eliminated.In this approach, the three-dimensional pressure and temperature fields serve directly as input and each grid cell is considered as a point source (i.e. a nucleus of strain) marked by a finite volume dV .The volumetric inflation V (positive for �T , �P > 0 ) of a point source is given by (3) in which the temperature and pressure-dependent vertical strains are obtained assuming a uniaxially compacting reservoir and are thus given by the expressions in Eq. 3.
Depending on the elastic properties of the medium, the volumetric inflation of the point source poses as its strength and allows for the computation of the subsurface displacement field (Okada 1992;Van Wees et al. 2019;Candela et al. 2022).The displacement on the fault plane is defined at regularly spaced (i.e.every 5 m) observation points (Fig. 2) in the dip direction along the fault pillars of the reservoir flow simulation grid (Van Wees et al. 2019;Candela et al. 2022).The poro-and thermo-elastic strain changes along the fault plane are obtained by integrating the contribution of each of the point sources, defined by the pressure and temperature changes in the reservoir flow model cells following Eq. 3.These elastic strain changes are subsequently used as sources for calculating induced stress changes on the fault's pillar (Fig. 2) using Eq. 1.Note that this approach eliminates the assumption of constant vertical stress and allows the effects of stress arching to be captured.A key aspect of MACRIS is the implementation of the Barnes-Hut algorithm (Barnes and Hut 1986).By re-discretizing the initial reservoir grid, the point source resolution further away and close to a fault are aggregated and refined, respectively, in order to increase the spatial resolution of predicted stress on the fault, while preserving computational performance (Van Wees et al. 2019).The MACRIS approach yields the full three-dimensional total stress response on the fault, including stress arching effects.The Coulomb stress CS on a fault is indicative of its stability and seismicity potential (Van Wees et al. 2019;Wassing et al. 2021) and is defined as (Zoback 2007) where σ s is the shear stress on the fault plane and µ is the friction coefficient of the fault.The normal stress value σ n is derived from the total stress tensor �σ ij obtained from either the uniaxial stress solution or the MACRIS approach.The effective normal stress σ ′ n on the fault plane is defined through Eq. 2 and accounts for the direct effect of the fluid pressure inside the fault on the total stress (Van Wees et al. 2019).In this work, (5)

Validation of thermo-elasticity in MACRIS
To validate the implementation of linear thermo-elasticity in MACRIS, the stress solution is compared with a 2D finite element model (FEM) solution extending from the benchmark study performed in Van Wees et al. ( 2019).For a laterally extensive reservoir and in the absence of conductive heat transfer, thermal cooling of the reservoir can be considered uniform.This allows use of the plane-strain assumption to (i) increase the efficiency of the benchmarking process whilst preserving a 3D representation of the reservoir model (Van Wees et al. 2019) and (ii) yield a solution that is in agreement with fully analytical solutions for laterally infinite reservoirs (Jansen et al. 2019).As such, a plane-strain elastic FEM of a geothermal reservoir subject to uniform cooling is constructed using DIANA, a commercial code widely used in engineering and subsurface applications (DIANA 10.1 User Manual 2016).The plane-strain finite-element mesh models a 2D section of 6 km width and 6 km depth and incorporates two reservoir compartments, separated by a fault with a throw of half the reservoir thickness (Fig. 3).The mesh is finest around the fault with a resolution of 1 m at reservoir depth and gets gradually coarser as you go away from the fault and the reservoir (for details of the FEM mesh and model setup see Van Wees et al. 2019).The FEM resolution has been chosen Fig. 3 Benchmark of MACRIS stress solutions for linear thermo-elasticity against the FEM solution from DIANA for a normal offsetting fault and a uniformly cooled reservoir.Note the stress peaks at the interfaces between formations as a result of the large temperature contrasts over the model depth interval.Model parameters and fault properties are as outlined in Table 1.The reservoir geometry is outlined in grey/black sufficiently high to capture high bending stresses close to the fault and at the interface of reservoir and adjacent rock.The model is laterally constrained on the sides (zero displacement in horizontal direction) and vertically constrained at the bottom (zero displacement in vertical direction), and subjected to gravity loading.These constraints effectively lead to a stress response on the sides of the model which corresponds to an infinite lateral extension of the reservoir.For MACRIS, a 3D representation of the faulted reservoir is used with a square geometry (similar to the model presented in Fig. 1) and with dimensions of 60 km laterally.The grid cell dimensions of the MACRIS model are 100 m × 100 m × 10 m in the x , y and z directions, respectively.The very large dimension of the 3D faulted reservoir is required to mimic the plane strain boundary condition of the 2D FEM model.The comparison between the 2D FEM and the central pillar of the 3D MACRIS model is presented in Fig. 3.The MACRIS results are in agreement with the high-resolution FEM solution.There exist, however, minor deviations at the transitions of the reservoir interval to the over-and underburden formations.The stress peaks observed at these transitions are a result of the sharp thermal boundary between the reservoir and the surrounding rocks, irrespective of fault offset.These stress peaks are likely to disappear as a result of conductive heat transfer from the surroundings into the reservoir, effectively smoothening the temperature contrast.MACRIS is thus shown to be capable of reproducing the high-resolution stress response with high accuracy.

The effect of temperature
The relative magnitude of thermo-and poro-elastic stresses during geothermal energy production can be quantified by rewriting Eq. 4 to obtain a ratio for the effective horizontal stress changes as a result of temperature and pressure changes (Segall and Fitzgerald 1998): Figure 4 (top) shows the horizontal stress ratio corresponding to the pressure and temperature solutions at mid-reservoir depth after 50 years of production in the case of no fault offset.The relative magnitude of thermo-elastic stresses is significantly larger (i.e.tenfold) compared to the poro-elastic effects within the cooled area.The effect of pressure is limited to the near-wellbore region evident by a local minimum of the horizontal stress ratio.As outlined by Van Wees et al. (2012), the coefficient of performance (COP) provides a measure of geothermal doublet efficiency in terms of the ratio of geothermal power and pumping power for fluid circulation in the injection and production wells and the reservoir.In the Netherlands, a COP of 15 is a representative target value (Van Wees et al. 2012).A restriction of 13.5 MPa km −1 (depth) on the allowed pressure changes in the injection well is simultaneously imposed by government (State Supervision of Mines and TNO-AGE, 2013).For the hydrostatic gradient listed in Table 1, this results in a targeted and maximum allowed pressure differential of ~ 3.6 MPa and ~ 7.3 MPa in the injection well, respectively.By adopting Peaceman's correction for the asymptotic decrease of the pressure differential with increasing radius from the well (White et al. 2013), the corresponding grid cell reservoir pressure differentials become ~ 1.8 and ~ 3.65 MPa at mid-reservoir depth (i.e.Fig. 4, bottom left), respectively.This (7) value provides a measure of the maximum magnitude of poro-elastic stresses in the near-wellbore region relative to static reservoir properties such as matrix permeability.
Figure 4 (bottom row) shows that the horizontal stress ratio remains positive as the relative magnitude of the poro-elastic stresses is increased at the injection well following a decrease in reservoir permeability.For high matrix permeability geothermal systems, thermo-elastic stresses and the effect of temperature are therefore shown to be dominant in the stress change induced within the cold-water volume.

Application of MACRIS to clastic reservoir models
For the two geological model scenarios defined in the previous section, the evolution of thermal and stress changes on the fault plane has been investigated in space and time.
The stress solutions obtained from MACRIS and the uniaxial solution are compared both in terms of absolute stress response as well as the degree to which MACRIS differs from the uniaxial solution.
Figure 5 presents the stress response on the fault plane in the case of no fault offset and focussed on the fault pillar positioned closest to the line connecting injector to producer well (see Fig. 1) marked by the earliest possible advent of the cold front on the fault plane.The previously observed stress peaks (Fig. 4) at the interface between reservoir and over-and underburden formations are absent in the case of non-uniform cooling.A gradual temperature contrast develops at these interfaces as a result of conductive heat transfer from the surroundings into the reservoir and yields a smooth stress signature over the model depth interval.On the lower fault segment (i.e.below mid-reservoir depth), MACRIS predicts a significantly higher shear stress (> 10%) at the arrival of the thermal front on the fault plane (i.e.20 years).The differential in thermal contraction that exists along the thermal front increments the shear stress; a phenomenon referred to as stress arching (Wassing et al. 2021).After passing of the thermal front, the difference in shear stress on the lower fault segment decreases as cooling on the fault plane becomes more uniform.The effect of stress arching is opposite on the upper fault segment as MACRIS predicts an increasingly lower shear stress over time (< − 10% to < − 18%).When thermal breakthrough occurs in the production well (i.e.50 years), thermal compaction of the reservoir has increased relative to the overburden formation.This causes the horizontal stresses in the overburden to increase, effectively reducing the shear stress (Geertsma 1973;Mulders 2003;Kivi et al. 2022).The temperature contrast at the interface between the formations becomes more gradual as conductive heat transfer progresses.Over time, this reduces the effect of stress arching directly at the interface.The effect of stress arching increments the Coulomb stress change in the reservoir at the arrival of the thermal front.After 50 years of production, both Coulomb stress solutions are shown to be in agreement on the reservoir interval.However, the MACRIS Coulomb stress solution does illustrate the continuous effect of stress arching in the overburden during progressive thermal compaction of the reservoir.
Figure 6 presents the stress response on the fault plane in the case of normal fault offset focussed on the pillar equal to Fig. 5.The difference between MACRIS and the uniaxial solution shows a similar pattern to the case of no fault offset at the arrival of the thermal front.However, in the case of normal fault offset, the difference in shear stress is further increased (> 15%) and more concentrated at mid-reservoir depth.Wassing et al. (2021) state that the effect of stress arching on the fault is increased when a sharp temperature front evolves along the fault plane.In those sections where the reservoir interval is not interconnected, steep temperature differentials develop over the fault plane.As a direct result, the increased effect of stress arching further increments the shear stress change within the interconnected reservoir interval.Note that the relative difference is most distinct at the top and base of the reservoir interval in the hanging-and foot-wall, respectively.The difference in shear stress on the reservoir interval does not decrease after passing of the thermal front as cooling on the fault plane remains nonuniform over time.At the time of thermal breakthrough, a considerable difference exists in the Coulomb stress solutions (~ 8%).When translated to a Coulomb Failure Function (Zoback 2007), this could lead to substantial misconceptions (i.e.fault failure or no fault failure) in the prediction of fault reactivation.Similar to the case of no fault offset, an increasingly lower shear stress is predicted by MACRIS in the overburden formation.Note that a sharp thermal boundary develops at the interface between the reservoir and basement formation in the hanging-wall.This causes the shear stress to decrease at the interface resembling the stress peaks observed in Fig. 4.
From the previously presented stress responses in Figs. 5 and 6, the effect of stress arching is most distinct in terms of shear stress.Therefore, focus is placed on the evolution of shear stress in order to further investigate the difference between MACRIS and the uniaxial solution.Figure 7 presents the shear stress solutions on the fault pillars at mid-reservoir depth.They are illustrated relative to the temperature solution in order to emphasize the dependency of the shear induced area on the intersection of the cold-water volume with the fault plane.In the case of no fault offset (Fig. 7, top), the shear stress on the fault pillars at thermal breakthrough shows, at first glance, a similar response between both approaches.The number of pillars on which a change in shear stress is induced coincides with the extent to which the cold-water volume laterally intersects the fault plane.However, subtle differences exist at the centre and boundaries of the intersected fault area.Following the temporal evolution of the thermal front, MACRIS consistently predicts a slightly higher (~ 5%) shear stress response as a result of stress arching effects.In the centre region, where cooling is more uniform, the stress response predicted by MACRIS is slightly lower (~ − 2%).This is a result of stress arching in the overburden, effectively reducing the shear stress at midreservoir depth.In the case of normal fault offset (Fig. 7, bottom), again the effect of stress arching becomes evident.The degree of cooling within the intersected area decreases, as compared to the case of no fault offset, and so does the number of shear induced pillars approximated by the uniaxial solution.Similar to the results presented in Fig. 6, considerable deviations in shear stress magnitudes (~ 18%) are observed at the central pillars of the fault as a result of stress arching effects.Moving towards the position of the thermal front, the deviation in the uniaxial solution reaches the point where no considerable shear stress change is induced on the fault pillars as compared to the MACRIS solution (~ 15%), e.g. the pillars at Y = 2100 m and/or Y = 2900 m.Including the effects of stress arching yields a significant shear induced stress change on these pillars even though they show less thermal cooling.As a result, the number Fig. 7 The shear stress solutions on the fault pillars at mid-reservoir depth for the case of no fault offset (top row) and normal fault offset (bottom row).(Right column) The relative difference in shear stress between both approaches over the production lifetime of the geothermal doublet of shear induced pillars increases relative to the lateral intersection of the cold-water volume as compared to the uniaxial solution and the case of no fault offset.
Figure 8 illustrates the difference in shear stress solution on the fault plane for the case of no (top row) and normal (bottom row) fault offset at the arrival of the thermal front on the fault plane (i.e.20 years) and when thermal breakthrough occurs in the production well (i.e.50 years).In both cases, the effect of stress arching clearly correlates to the intersection of the thermal front with the fault plane.After 20 years of production, the differences between MACRIS and the uniaxial solution are distributed similarly to the results presented in Figs. 5 and 6 for the centre fault pillar.In the case of no fault offset, the incremented shear stress on the lower fault segment follows the lateral movement of the thermal front as the intersected area increases over time.The effect of stress arching in the overburden that causes the difference in shear stress on the lower fault segment to reduce is limited to the interior of the cold-water volume where temperature is more uniform.A similar but more explicit pattern is observed in the case of normal fault offset.As a result of the increased effects of stress arching, the difference in shear stress change further increases (> 20%) within the interconnected reservoir interval and appears independent of the intersected area.The previously presented results show that the contribution of stress arching increases when a steep temperature gradient exists at the thermal front, and demonstrate that uniform cooling and conductive heat transfer into the reservoir reduce the effects of stress arching.To further elaborate on these findings and emphasize the three-dimensional effects of long-term cooling of the fault plane (1) the reservoir thickness is decreased (Fig. 9); (2) the reservoir formation is intercalated with low permeability layers (Fig. 10); (3) the injection well is located closer to the fault (Figs. 11 and 12); and (4) the doublet is oriented parallel to a sealing fault (Fig. 13).
Fig. 8 The relative difference in shear stress solutions against the temperature profile on the fault plane for the case of no fault offset (top row) and normal fault offset (bottom row) as seen from the position of the injection/production well.A positive difference indicates MACRIS to predict a larger stress response compared to the uniaxial solution.The reservoir geometry is outlined in grey/black where the solid and dashed lines mark the reservoir interval in the foot-and hanging-wall, respectively Figure 9 presents the stress response for a reservoir thickness of 40 m in the case of no (top row) and normal (bottom row) fault offset for the pillar that is in line with the orientation of the geothermal doublet (Fig. 2).The well flow rate has been adjusted to 115 m 3 /h in order to keep the timing of thermal breakthrough constant.This also ensures that, in terms of convective heat transport only, both the lateral extent to which the fault plane is intersected by the thermal front as well as the fluid velocity of the propagating cold-water volume remain equal to the originally presented scenario.Similarly, as an identical temperature differential is imposed on the reservoir, the rate of conductive heat transfer and the affected depth interval of the reservoir remain equal as well.Consequently, when reservoir thickness decreases, the relative contribution of conductive heat transfer increases creating an increasingly elongated shape of the thermal front to allow for a more gradual temperature gradient.In both cases, this not only reduces the magnitude of the absolute induced shear stress change on the fault pillar, it also reduces the effect of stress arching at the arrival of the thermal front on the fault plane (i.e.20 years).Over time (i.e.50 years), however, the increasingly elongated thermal front is shown to prolong the effect of stress arching along the fault pillar.In the case of normal fault offset, the effect of stress arching itself is increased similar to previously presented results.
Considering the original depth interval of 100 m, the reservoir is intercalated with 2 evenly spaced low permeability (i.e. 1 mD) layers of 20 m thickness resulting in 3 isolated reservoir zones of equal thickness.The well flow rate has not been adjusted in this case as the open-hole section of the wells covers the entire depth interval of the reservoir zone.Figure 10 presents the stress response for an intercalated reservoir in the case of no (top row) and normal (bottom row) fault offset.Similar to the results presented in Fig. 9, the diffusive shape of the thermal front decreases the effect of stress arching within the individual reservoir zones.The existence of isolated high permeability layers within the reservoir interval creates so-called thief-zones resulting Fig. 13 Results for the scenario in which the geothermal doublet is oriented parallel to a bounding fault.(Top row) the shear stress solutions on the fault pillars at mid-reservoir depth when thermal breakthrough occurs in the production well in the case the doublet is located at a distance of (left) 750 m and (right) 300 m from the fault.(Bottom row) the corresponding stress and temperature solutions on the fault plane, focused on the pillar that is in line with the injection well (i.e.Y = 1750 m ).The reservoir geometry is outlined in grey/ black in localized cooling on the fault plane.The stress signatures are observed to adhere strictly to the high permeability layers, especially at the arrival of the thermal front (i.e.20 years).In the case of no fault offset, the effect of stress arching is observed to increment the shear stress on the lower fault segment similar to the results presented in Fig. 5.In the case of normal fault offset, the spatial development of the cold-water volume is limited to the lower and interconnected reservoir zones.A steep temperature differential develops only at the base of the reservoir interval where the effects of stress arching are observed to increment the shear stress similar to the shear stress response presented in Fig. 6.At the top of the reservoir interval, conductive heat transport governs the stress change and, as a result, the effects of stress arching are decreased.In terms of the Coulomb stress solutions between both approaches, similar conclusions can be drawn as to the scenario of a uniform reservoir interval (Figs. 5  and 6).
Similar to Fig. 5, Fig. 11 presents the stress response on the fault plane in the case of no fault offset and focussed on the fault pillar in line with the geothermal doublet (Fig. 2).The location of the geothermal doublet is shifted such that the injection well is located at a distance of 300 m from the fault whilst maintaining the original doublet spacing.Although not addressed separately at this point, note that this configuration will increase the intersection area of the cold-water volume with the fault plane.As the relative contribution of conductive heat transfer is limited during the early stages of cold-water injection, a steep temperature gradient exists at the rim of the cold-water volume.The resulting increase in the degree and uniformity of cooling causes the magnitude of the induced stress change and the effects of stress arching to increase (Wassing et al. 2021).At the arrival of the thermal front on the fault plane (i.e. 5 years), the observed shear stress signatures of the uniaxial and MACRIS solutions show a similar pattern in which the effects of stress arching have increased significantly compared to the results presented in Fig. 5.The differential in thermal contraction increases not only along the thermal front, but also between the reservoir and overburden formation.This effectively decreases the shear stress response over the model depth interval.A consistently and progressively lower shear stress (~ − 10%) is predicted by MACRIS on the reservoir interval and lower fault segment, respectively, after passing of the thermal front (i.e. 15 years).Where the effects of stress arching were previously observed to increment the Coulomb stress change, the reduced shear stress causes a reduction in the Coulomb stress change on the reservoir interval.Initially, a similarly steep temperature contrast exists at the transition between reservoir and over-and underburden formations, evident by the sharp stress signatures over the model depth interval.Similar to the results presented in Fig. 5, the ongoing conductive heat transfer into the reservoir causes the stress signatures to smoothen over the transitions between the reservoir and over-and underburden formations.Moving towards the time of thermal breakthrough in the production well (i.e.50 years), the effects of stress arching on the reservoir interval decrease and both Coulomb stress solutions are shown to be in agreement.
Figure 12 presents the stress response on the fault plane in the case of normal fault offset focussed on the pillar equal to Fig. 11 with the injection well located at a distance of 300 m from the fault whilst maintaining the original doublet spacing.Similar to the case of no fault offset, the increased degree of cooling causes the magnitude of the induced stress change and the effects of stress arching to increase.The stress signatures of the uniaxial and MACRIS solutions are similar to those presented in Fig. 6, showing the shear stress to be incremented most within the interconnected reservoir interval.Similar to the results presented in Fig. 11, the increased effects of stress arching in the overburden are observed to consistently decrease the shear stress response in the surrounding formations.In the case of normal fault offset, however, the differential in thermal contraction originating from the juxtaposition of the reservoir is shown competent to counteract the effects of stress arching in the overburden and increment the induced shear stress change.As a result, the Coulomb stress change remains incremented during the interaction of the cold-water volume with the fault plane.At the time of thermal breakthrough in the production well (i.e.50 years), a considerable difference exists in the Coulomb stress solutions (> 18%) around mid-reservoir depth.
To corroborate the importance of the three-dimensional spatial development of the coldwater volume relative to the fault, the configuration of the geothermal doublet and fault offset are altered.Figure 13 presents the results for a doublet configuration in which the line connecting injector and producer wells is oriented parallel to and at a distance of 750 m (left) and 300 m (right) from a bounding fault.The reservoir flow is compartmentalized by an increase in fault offset equal to the reservoir thickness.With zero fault transmissibility, the injected cold-water volume is deflected by the fault plane and its development is increasingly directed towards the production well.In the case the doublet is located at 750 m from the fault, the time of thermal breakthrough in the production well is decreased to 40 years.At this time, MACRIS predicts a minor induced shear stress change regardless of the fact that the thermal front has not (yet) reached the fault.Thermal stresses within the cold-water volume appear to be transferred to the nearby fault, similar to the findings presented by Candela et al. (2018) and Kivi et al. (2022).In the case the doublet is located at 300 m from the fault, although similarly deflected, the cold-water volume is in direct contact with the fault plane without intersecting it.The number of shear induced fault pillars coincides with the contact area of the cold-water volume with the fault plane.When thermal breakthrough occurs in the production well (i.e. 30 years), a steep temperature gradient exists on the juxtaposed reservoir interval.The increased degree of cooling and the sharp temperature contrast are observed to similarly increase the effects of stress arching in the overburden formation (Figs. 11 and 12) and on the reservoir interval (Figs. 6 and 12), respectively.The induced shear stress change on the fault plane is, however, incremented even further (~ 40% at mid-reservoir depth) as a result of the increase in normal fault offset.Similar to the results presented in Fig. 7, the effects of stress arching causes both the number of shear induced fault pillars as well as the induced stress magnitude on these pillars to increase relative to the lateral contact area of the cold-water volume as compared to the uniaxial solution.For both cases, the incremented stress response predicted by MACRIS highlights the importance of accurately accounting for the effects of stress arching on bounding faults in seismic hazard assessment.

Discussion and conclusion
The model capabilities of MACRIS have been successfully extended to account for thermo-elastic stress changes in complex geothermal reservoirs.The semi-analytical stress solution for thermal cooling is shown to be in agreement with the FEM when reservoir throw is considered.The comparison of MACRIS with the analytical solution for uniaxial reservoir compaction highlights the need for accurate predictions of subsurface stress changes in complex geothermal reservoirs.MACRIS is shown able to predict the stress development on faults during long-term cooling of the subsurface by accurately capturing the effects of stress arching at the position of the thermal front (taking into account thermal diffusion) and when reservoir complexity increases.The implementation of thermo-elasticity adheres to the computational approach adopted in MACRIS making it suitable for the evaluation of thermal stresses in heavily faulted reservoirs (Van Wees et al. 2019;Candela et al. 2019Candela et al. , 2022)).To resolve the stress evolution on complex faults at the same accuracy as MACRIS, the computational cost of a dedicated FEM (i.e.100-500k cubic order elements required to capture high-resolution stress variations) is at least one order of magnitude higher (Van Wees et al. 2019).Using 8 logical processors, MACRIS is able to compute the stress development on all model fault pillars (Fig. 1) for ten timesteps in approximately 180 s.Furthermore, MACRIS enables parallel computing on individual fault pillars resulting in a reduction of the computational time that is proportional to the added number of logical processors.
MACRIS adopts an explicit coupling approach between the dynamic behaviour and mechanical response of the reservoir.As such, a static model is considered in which the mechanical parameters do not change with variations in reservoir pressure and temperature.Particularly in EGS, induced thermo-elastic stress changes affect reservoir performance by altering rock properties, including fracture aperture and permeability, and lead to long-range stress interactions with nearby pre-existing fault zones (Atkinson et al. 2016;Candela et al. 2018;Rathnaweera et al. 2020).In complex geothermal reservoirs characterized by high matrix permeability, the absence of fractures eliminates the effect of fracture dominated long-range stress interactions allowing for a static mechanical representation of the reservoir.However, short-range stress interactions do arise when the cold-water volume is deflected from the fault plane as illustrated in Fig. 13.Long-term cooling causes thermal contraction of the reservoir rock which can affect the matrix permeability.This, in turn, may affect the hydraulic propagation of the injected cold-water volume and the relative contribution of conductive heat transfer at the thermal front.Although not considered in this work, incorporating reservoir creep (or inelastic strain) as part of the thermal compaction strains could reduce the subsurface stress response and the associated potential for fault reactivation significantly.Since the stress response is driven primarily by the elastic component of strain (Van Wees et al. 2018), the predicted stress could be lowered by the ratio of creep to total strain.For gas depletion in similar clastic reservoirs creep contributes to approximately 50% of the total strain (Orlic and Wassing 2013;Pijnenburg et al. 2018), implying predicted stress changes should be half of those inferred from the fully elastic solution.The effects of creep can be effectively incorporated in MACRIS by treating nuclei of inelastic thermal strains as a convolution over time of the initially fully elastic strains as proposed for pressure depletion (Van Wees et al. 2018).In addition, creep could result in an apparent stiffness contrast (Van Wees et al. 2018) that, except for peaks at the interfaces between formations, would affect the predicted stress changes in a very moderate way.An indepth investigation of the potential effects of creep is considered beyond the scope of this paper.
In high matrix permeability clastic geothermal reservoirs, long-term cooling can cause a significant increase in Coulomb stresses on faults in the vicinity of a geothermal doublet.Similar to earlier findings on the importance of thermal stress changes associated with cold-water injection (Candela et al. 2018;Wassing et al. 2021;Kivi et al. 2022), the effect of temperature is shown to be dominant in the subsurface stress response.The magnitude of the thermo-elastic stress response is limited by the degree of cooling of the rock matrix.The potential for fault reactivation is therefore primarily dependent on the extent to which the cold-water volume intersects or is in direct contact with the fault plane and the injection temperature can provide operational control on the induced stress change and subsequent seismic hazard.Furthermore, as illustrated in Figs. 5 and  6, continuous conductive heat transfer (or thermal diffusion) into the reservoir increasingly smoothens the stress signatures at the interfaces between the reservoir and surrounding formations.
The effects of stress arching are shown to benefit from the development of a steep temperature gradient at the thermal front, similar to the findings presented by Wassing et al. ( 2021).An increase in conductive heat transfer into the reservoir is observed to counteract stress arching as the thermal front attains a more diffusive shape.Spatial deviations in the uniformity of the temperature gradient along the thermal front lead to localized shear stress increments (i.e.Fig. 10), highlighting the need for high-resolution stress solutions on the fault.Numerical diffusion is known to affect the temperature solution (Wesseling 2001) by introducing a certain amount of inaccuracy in the position and steepness of the hydraulically propagating thermal front.Although not presented in this work, the grid resolution of the geological model in Fig. 1 has been varied to investigate the effect of numerical diffusion on stress arching phenomena.As the horizontal grid resolution coarsens, a more gradual temperature gradient develops over the thermal front.Similar to an increase in conductive heat transfer, this reduces both the absolute magnitude of the induced stress change as well as the effects of stress arching on the fault plane.Coarsening of the vertical grid resolution introduces a minor deviation (< 2%) in the absolute magnitude of the stress change at mid-reservoir depth, i.e.where a peak in the stress signature exists.Both the effects of stress arching as well as the stress signature over the model depth interval are unaffected by coarsening of the vertical resolution.Following a mesh-free approach, the vertical resolution in MACRIS is defined independently from the reservoir flow model resolution.Although high-resolution stress predictions can still be obtained using MACRIS, the stress signature will adhere to the flow model resolution on which the temperature field is obtained.
Fault throw and fault related reservoir flow compartmentalization are shown to strongly affect the induced stress changes along the fault and enhance the effects of stress arching over the model depth interval.These findings are in agreement with previous studies of stress arching in depleting oil and gas reservoirs (Mulders 2003;Orlic and Wassing 2013;McGarr 2014;Van Wees et al. 2014, 2017;Candela et al. 2019Candela et al. , 2022)).Continuous conductive heat transfer along the fault plane is observed to smoothen the large temperature contrasts on the juxtaposed reservoir sections along the fault and subsequently reduce the effects of stress arching over time.In all cases, reservoir flow compartmentalization is shown to strongly affect the three-dimensional spatial development of the cold-water volume.Flow paths and the associated thermal front are deflected from the fault, effectively shielding faults from cooling.Consequently, as illustrated by Fig. 13, in most cases the degree of cooling and the induced stress changes on the fault plane will be considerably lower than for transmissive faults.Note, however, that despite the lesser degree of cooling the effects of stress arching increase as a result of the increase in fault throw.Furthermore, when the fault acts as a barrier to doublet flow, the assumption of radially symmetric flow used in many 2D analytical models incorporating fault throw becomes invalid (Mulders 2003;Fokker and Orlic 2006;Fjaer et al. 2008).
The assumption on fault transmissibility in the model setup provides maximum connectivity over the fault plane in order to highlight the effect of fault throw on the induced stress changes.Lowering the fault transmissibility by defining, i.e. a shale gauge ratio will result in a predominantly sealing character of the fault and promote flow compartmentalization. Therefore, similar to the results presented in Fig. 13, a reduction in the fault hydraulic properties will affect the spatial development of the cold water volume and enhance the effects of stress arching on the fault plane.
Uniform elastic properties of the subsurface and no elastic contrast between the reservoir and surrounding formations are assumed.Both the magnitudes as well as the patterns of stress changes can alter significantly as a result of vertical and lateral contrasts in elastic properties (Fokker and Osinga 2018;Van Wees et al. 2018).Considering an increase in the Young's modulus of the over-and underburden formations, effectively stiffening the reservoir surroundings, the change in normal and shear stresses at the interfaces between the reservoir and surrounding formations are expected to increase.The absolute magnitude of the induced stress change in the reservoir will remain similar as the reservoir stiffness is unchanged.However, as a result of the increased stress changes on the interfaces, the stress response is expected to become more uniform on the reservoir interval.Similar to the results presented in Figs.11 and 12, this will sharpen the stress signature over the model depth interval.A lateral contrast in elastic properties arises on those sections of the fault plane where the reservoir interval is not interconnected.The increased stiffness of the surrounding formations yields an increase in the horizontal stress that is expected to enhance the effects of stress arching on the fault plane, similar to previous studies of stress arching in depleting oil and gas reservoirs (Mulders 2003;Van Wees et al. 2018).
Stress arching at the thermal front is shown to increment the shear stress on the lower fault segment.Wassing et al. (2021) observe a similar effect, but on the upper fault segment.As illustrated by Kivi et al. (2022), this is due to the dip direction of the fault with respect to the location of the injection well and irrespective of fault throw.Fluctuations in fault throw along fault strike will affect the degree of cooling and the effects of stress arching on the fault plane locally.The results presented in this work consider a critically stressed fault.As the orientation of the minimum horizontal stress with respect to fault strike changes, or vice versa, the slip-tendency of the fault will decrease.Under similar operational conditions, the potential for fault reactivation and seismic hazard decreases regardless of a coincidental increase in the intersection area of the cold-water volume with the fault plane.A semi-analytical approach to estimate the potential cumulative seismic moment from the elastic stress solution is proposed by Van Wees et al. (2018), in which an average excess Coulomb stress is redistributed relative to the Mohr-Coulomb failure criterion to determine the slip length along the fault plane.Following this approach, the observed reduction in Coulomb stress in the surrounding formations as a result of stress arching may provide a buffer for the build-up of seismic moment on the fault by limiting the slip length along the fault plane.An increase in the effective stress ratio would shift the Mohr-Coulomb failure envelope to represent a lower slip-tendency of the fault, effectively decreasing the potential for fault reactivation and seismic hazard.As the effect of temperature is dominant, the value of the linear thermal compaction coefficient is expected to significantly affect the magnitude of the resulting stress change.A lower thermal compaction coefficient could additionally reduce the differential in thermal contraction and associated stress arching effects at the rim of the cold-water volume.Although beyond the scope of this work, addressing the sensitivities in fault geometry and parameter values would augment the robustness of the presented findings.
The results presented in this work apply predominantly to low enthalpy geothermal systems that are characterized by sedimentary reservoirs with high matrix permeability and are representative for geothermal exploitation in the Netherlands (Buijze et al. 2023).The development of geothermal doublets in these systems often involve well placement within singular fault blocks resulting in cold-water volume interactions with (i) sealing or bounding faults as addressed in Fig. 13 or (ii) (non-) offsetting faults located in between, and at various distances to, the wells as addressed in Figs. 5,6,11 and 12.In terms of stress arching effects, the presented results can be considered representative for fracture permeability dominated geothermal systems when focusing solely on the spatial and temporal development of the injected cold-water volume (Wassing et al. 2021).
In conclusion, the model capabilities of MACRIS have been successfully extended to capture the thermo-elastic stress changes in structurally complex (i.e.faulted) geothermal reservoirs.Similar to previous studies by Wassing et al. (2021) and Kivi et al. (2022), the reactivation potential of faults in geothermal reservoirs is predominantly governed by the intersection or contact area of the injected cold-water volume with the fault plane.Moreover, results show the temperature contrast along the thermal front to yield stress arching effects which are (i) enhanced as a result of reservoir throw and flow compartmentalization and (ii) reduced by a relative increase in conductive heat transfer between the reservoir and surrounding formations.

Fig. 1
Fig. 1 Schematic representation of the low enthalpy geothermal model with a high permeability clastic reservoir (marked in orange colour) for (left) the scenario with no fault offset and (right) the scenario with normal fault offset.The fault plane is outlined in red contour and dips in the direction of the minimum horizontal stress.The blue and red bars at the top of the model indicate the location of the injection and production well, respectively

Fig. 2
Fig.2Modelling workflow for geothermal reservoir studies indicating relative contributions of pressure and temperature on the subsurface stress response.The computational methodology of MACRIS is schematically outlined in the top right corner illustrating the use of sources and receivers for the calculation of induced stress changes on the fault's pillar (modified fromCandela et al. 2022)

Fig. 4 (
Fig. 4 (Top) the horizontal stress ratio over the line connecting injection to producer well relative to the absolute pressure and temperature solutions after 50 years of production in the case of no fault offset.(Bottom row) the effect of varying reservoir matrix permeabilities on the pressure solution (Van Wees et al. 2012) and the horizontal stress ratio along the orientation of the geothermal doublet.The horizontal stress ratio around the injection well is highlighted to emphasize the near wellbore region.The fault geometry is outlined in grey

Fig. 5 (
Fig. 5 (Top row) stress and temperature solutions on the fault plane in the case of no fault offset at (left) the arrival of the thermal front on the fault plane and (right) when thermal breakthrough occurs in the production well.(Bottom row) the relative difference between both approaches on the fault pillar is presented over the production lifetime of the geothermal doublet.A positive difference indicates MACRIS to predict a larger stress response compared to the uniaxial solution.The reservoir geometry is outlined in grey/ black

Fig. 6 (
Fig. 6 (Top row) stress and temperature solutions on the fault plane in the case of normal fault offset at (left) the arrival of the thermal front on the fault plane and (right) when thermal breakthrough occurs in the production well.(Bottom row) the relative difference between both approaches on the fault pillar is presented over the production lifetime of the geothermal doublet.A positive difference indicates MACRIS to predict a larger stress response compared to the uniaxial solution.The reservoir geometry is outlined in grey/ black

Fig. 9
Fig.9Stress and temperature solutions on the fault plane for a reduced reservoir thickness in the case of no fault offset (top row) and normal fault offset (bottom row).The relative difference in shear stress solutions between both approaches is presented over the production lifetime of the geothermal doublet.A positive difference indicates MACRIS to predict a larger stress response compared to the uniaxial solution.The reservoir geometry is outlined in grey/black

Table 1
Model parameters for the evaluation of the subsurface stress field using MACRIS and the uniaxial solution for reservoir compaction