The risk of fluid-injection-induced fault reactivation in carbonate reservoirs: an investigation of a geothermal system in the Ruhr region (Germany)

In this study, we investigate the potential for fluid-injection-induced fault reactivation and induced seismicity risk during simultaneous injection-extraction operation in a theoretical geothermal doublet system in a carbonate reservoir in the Ruhr region. Using a coupled three-dimensional thermo-hydro-mechanical approach, we investigate the probability of injection-induced fault rupture. We perform a sensitivity study assuming variability of the fault and matrix permeability, injection/production flow rates, well placement options, rock thermal properties, and evaluate the influence of thermally induced stresses. The ruptured fault areas were calculated based on a Coulomb friction law and a notion that the shear slip is controlled by the ratio of shear to effective normal stresses acting on a pre-existing plane of weakness in the in situ stress field configuration. Ruptured fault areas in the intrinsically not critically-stressed environment, using location-specific empirical correlations, were used to compute local moment magnitudes of potential earthquakes. Based on this study, we conclude that, in the long-term, thermally-induced stresses play a dominant role during fault reactivation and greatly increase the likelihood for induced seismicity. We, therefore, propose that a minimum safe distance between an injection well and a fault should be based primarily on the radius of a thermal plume generated during the expected lifetime of a geothermal system. Results from this study provide valuable insights for the development of future deep geothermal systems in the Ruhr region and other geothermal reservoirs worldwide. We perform numerical simulations of fluid-injection-induced fault reactivation in a theoretical geothermal doublet system in a carbonate reservoir. Ruptured fault areas, in the intrinsically not critically-stressed environment, using location-specific empirical correlations were used to compute local moment magnitudes of potential earthquakes. In the long term, thermally-induced stresses play a dominant role during fault reactivation and greatly increase the likelihood of induced seismicity. A minimum safe distance between an injection well and a fault should be based primarily on the radius of a thermal plume generated during the expected lifetime of a geothermal system. We perform numerical simulations of fluid-injection-induced fault reactivation in a theoretical geothermal doublet system in a carbonate reservoir. Ruptured fault areas, in the intrinsically not critically-stressed environment, using location-specific empirical correlations were used to compute local moment magnitudes of potential earthquakes. In the long term, thermally-induced stresses play a dominant role during fault reactivation and greatly increase the likelihood of induced seismicity. A minimum safe distance between an injection well and a fault should be based primarily on the radius of a thermal plume generated during the expected lifetime of a geothermal system.

Abstract In this study, we investigate the potential for fluid-injection-induced fault reactivation and induced seismicity risk during simultaneous injection-extraction operation in a theoretical geothermal doublet system in a carbonate reservoir in the Ruhr region. Using a coupled three-dimensional thermohydro-mechanical approach, we investigate the probability of injection-induced fault rupture. We perform a sensitivity study assuming variability of the fault and matrix permeability, injection/production flow rates, well placement options, rock thermal properties, and evaluate the influence of thermally induced stresses. The ruptured fault areas were calculated based on a Coulomb friction law and a notion that the shear slip is controlled by the ratio of shear to effective normal stresses acting on a pre-existing plane of weakness in the in situ stress field configuration. Ruptured fault areas in the intrinsically not critically-stressed environment, using location-specific empirical correlations, were used to compute local moment magnitudes of potential earthquakes. Based on this study, we conclude that, in the long-term, thermally-induced stresses play a dominant role during fault reactivation and greatly increase the likelihood for induced seismicity. We, therefore, propose that a minimum safe distance between an injection well and a fault should be based primarily on the radius of a thermal plume generated during the expected lifetime of a geothermal system. Results from this study provide valuable insights for the development of future deep geothermal systems in the Ruhr region and other geothermal reservoirs worldwide.

Introduction
Geothermal systems are believed to be promising sources of clean and renewable energy, being it either heat or electricity, especially for resource-poor countries. For the development of, especially lowpermeability, geothermal systems, reservoir enhancement, in a form of hydraulic, chemical, or thermal stimulation is often necessary (Huenges and Ledru 2011). Any geothermal system, either hydrothermal or enhanced have to sustain economic flow rates and provide sufficient economical outputs for a minimum period of 20 to 30 years (Tester et al. 2006;Breede et al. 2013). During the lifetime of a geothermal system, fluid circulation will alter the initial in situ stress state. In the presence of faults, these changes could lead to their shear reactivation and, as a consequence, induction of seismic events. One can distinguish between fault reactivation caused by either direct or indirect fluid-injection-or fluid-productioninduced stress changes. In the direct reactivation of faults, fluid pressure changes perturb the effective normal stress exerted on a given fault and induces its slip. In the indirect reactivation, change of the local thermo-poro-elastic stress regime due to fluid injection or production causes a fault to slip, even in the case if a fault is not or only weakly hydraulically connected to the injection area (Berre et al. 2020). Predicting fault reactivation during simultaneous fluid injection and production is especially challenging, as it inherently involves the complexity of fluid migration and thermo-poro-elastic disturbance of the in situ stress state. Additionally, chemical effects, related to the fluid-rock interaction, e.g., frictional weakening of faults (Hirono et al. 2013), could be of high significance in the case of carbonate reservoirs. A deeper understanding of the above-mentioned processes has a great value for understanding and predicting the induced seismicity hazard.
While many of the fluid-induced seismic events in geothermal areas remain moderate in size and could be considered relatively harmless, some have caused considerable structural and financial damages, injuries, people displacement, the resistance of the local population, redirection of drilling targets, and project abandonment. As an example, seismic event with a local magnitude, M L , of 2.4 in Landau from 2009 and 2.2 in Insheim from 2010 (both located in the Upper Rhine Graben, Germany), related to the geothermal system promoted changes to the operation practises (Groos et al. 2013;Buijze et al. 2019). Noticeable microseismicity in the aseismic Greater Munich area (southern Germany) related to the geothermal use in the city of Munich, lead to the abandonment of previous fault-related drilling targets and redirecting geothermal well paths into lithofacies of higher permeability (Böhm et al. 2021). The 5.5 M W earthquake in Pohang in South Korea from 2017, believed to be triggered by high overpressures created during drilling and following stimulation operations at, previously unknown, fault within an enhanced geothermal system system two months after stimulation operations terminated, caused injuries, significant structural damages, and forced people displacement. As concluded by Lee et al. (2019), seismic risk in the south-korean project, which violated the so-called 'volume hypothesis', could have been evaluated months before the earthquake, had the presence of the undetected fault, and its susceptibility to slip in the prevailing stress state, been recognised. It is, therefore, clear that a detailed understanding of the stress state change during fluid injection-production as well as its evolution in time and space, with consideration of hydraulic, thermal, and mechanical characteristics of the reservoir, is necessary. The forecasting simulation techniques, as attempted in this study, could a priori help to mitigate induced seismicity risk of deep geothermal systems, especially for densely populated areas, of which the Ruhr region is a primary example. This study focuses on simulating simultaneous fluid injection and depletion for a period of 20 years in a theoretical geothermal doublet system located in the southern Ruhr region to assess the possibility of fault reactivation and evaluate the induced seismicity risk. We aim at quantifying pressure and stress changes on faults, located close to a doublet system, with a physics-based model, assuming both isothermal and non-isothermal cases, focusing on a few endmember scenarios with varying (i) well placement, (ii) fault and (iii) reservoir rock permeability, (iv) production and injection flow rate, and (v) rock thermal properties. Due to the lack of available studies on fault-related properties, for the assessment of seismic hazard related to a potential fault slip, two end-member fault cohesion cases, imitating elastic clay-and brittle calcite-filled fault, were assumed. The undisturbed in situ stress conditions, prior to the fluid circulation, were extracted from a regional static 3D geomechanical model of the region. The results from this study contribute to a deeper understanding of the physical processes of fault reactivation in the Ruhr region and other geothermal reservoirs worldwide.

Geological setting and reservoir characteristics
The coal-bearing strata of the Ruhr region belong to the external fold and thrust belt of the latest stage of the Variscan orogen, located east of the Wales-Brabant Massif (Brix et al. 1988). The Variscan orogenic belt was formed during the Late Paleozoic convergence of the Euramerican and Grondwanaland continental masses (Ziegler 1990). The collision and convergence took place during the Devonian and Carboniferous geologic periods.
The Ruhr region is influenced by two main fault systems. The major faults strike in the NE-SW direction with dips ranging from steeply inclined to bed-parallel reaching lengths of up to 40 km. These reverse type faults have horizontal displacements in the order of tens to hundreds of meters with few reaching horizontal displacements of ∼ 2.5 km . The reverse faults are dissected by multiple NW-SE and N-S oriented normal faults resulting in a horst-andgraben structure with few strike-slip faults of various orientations also being observed in the region (Brix et al. 1988). Folds, mainly NE-SW oriented, vary significantly in shape and dimension with wavelengths of up to 10 km , increasing towards the north, having amplitudes in the order of several hundred meters.
Extensive hard coal mining activities exposed thick molasse-type clastic sediments of the Upper Carboniferous age ranging from shales, silt-to coarse-grained sandstones with 2 to 2.5 m thick coal seams of varied strength, all heavily deformed by thrusting and folding (Bachmann et al. 1971). North of the Ruhr region, thick, in places up to ∼ 1.8 km , Cretaceous strata overlays the Upper Carboniferous sequences (Drozdzewski 1993). Four deep exploratory boreholes, located north of the Ruhr region i.e., Münsterland 1 (Hesemann 1965), Versold 1, Isselburg 3 (Drozdzewski 1993), and Vingerhoets 93 (Eder et al. 1983) have reached Devonian carbonates. Carbonate layers of the Middle and Upper Devonian age, being part of the Devonian Reef Complex, are outcropping south of the Ruhr region close to the city of Iserlohn.
Based on laboratory studies carried out on rock samples from the Ruhr region, it can be concluded that low matrix permeability of < 10 −18 m 2 characterizes sandstones of the Upper Carboniferous age (Stöckhert 2015;Nehler et al. 2016;Brenne 2016;Alber and Meier 2020;Stoeckhert et al. 2020;Duda and Renner 2012) and permeability of < 10 −15 m 2 characterizes Devonian limestones and dolomites . Matrix porosity of both rock formations does not exceed, on average, 5 to 7 % (Thielemann et al. 2001;Jorand et al. 2015;Stoeckhert et al. 2020) with local increases of permeability being connected to the weaker coal layers (Alber and Meier 2020;Thielemann et al. 2001). Due to the low matrix porosity and permeability of the host rocks in the Ruhr region, it is expected that the future geothermal systems will rely on the secondary permeability, being it either lithofacies of increased hydraulic properties, pre-existing fracture networks, and/or faults within, primarily, Devonian carbonate units.

Numerical modeling approach
The numerical model used in this study simulates single-phase fluid flow, following Darcy's law, in faults, represented as 2D planes of weakness cutting the model volume, and the surrounding poro-elastic rock mass. Fluids are being injected into the reservoir (i.e., Devonian carbonate layer), through an injection and produced from a production well, within a so-called geothermal doublet system, under the mass balance principle. The slip on faults was simulated based on Coulomb friction law and a notion that shear slip is controlled by the ratio of shear to effective normal stresses acting on a pre-existing plane of weakness in the stress field configuration and coupled with thermo-poro-elastic response of the porous medium to the fluid circulation. The numerical simulations were performed using the COMSOL Multiphysics software (COMSOL AB 2021), which accounts for pore pressure, P p , and associated stress and strain changes, where the coupling strength is governed by the Biot's coefficient, . The applied fluid flow governing equations resulted from conservation of momentum and mass in a fully saturated porous medium, whereas the heat transport governing equations, were derived from the heat balance considering both advective and conductive heat transport. The full theory behind how fluid circulation affects the effective stress of pre-existing discontinuities in a geothermal reservoir using thermo-hydro-mechanical simulating approaches can be found in COMSOL AB (2021) or Kolditz et al. (2012). The effects of fault permeability enhancement due to the shear dilatation (e.g., Eyinla and Oladunjoye 2021), change of rock properties due to P p or temperature, T, influence (e.g., Ahrens et al. 2021), influence of fluid chemistry on rock mass and fault properties (e.g., Bächler and Kohl 2005), mechanisms of earthquake interactions (e.g., Catalli et al. 2016), and Kaiser effect (e.g., Li and Nordlund 1993) were neglected.
To fully evaluate the risk of fluid-induced fault slip and resultant seismicity in the Ruhr region and simultaneously account for a wide range of uncertainties related especially to the fault and matrix permeability, end-member scenarios were evaluated. These included (i) different well placement options (i.e., scenarios with injection well penetrating a fault and production well being located farther away from the same fault, and scenarios with production well penetrating a fault and an injection well located away from the same fault); (ii) different fault, k f , and (iii) reservoir, k res , permeability; iv) two extreme cases of injection and production rates based on common circulation rates found in geothermal systems worldwide (Evans et al. 2012; v) two extreme cases of thermal expansion coefficient of the reservoir rock assuming (a) a carbonate material with high amount of quarzite and/or clay and/or an influence of bedding and/or exposure to high reservoir temperature (b) pure and homogenous carbonate material. Additionally, an isothermal case, i.e., without thermal influence resulting from an injection of a colder fluid into the hotter reservoir, was explored to assess the significance of thermally induced stresses. To investigate the sensitivity of different input parameters on the fluid-induced fault reactivation potential, we analysed 14 cases in total, where each case shows simulation result with a change of one critical parameter leaving all other constant (Table 1). All case scenarios were performed for an assumed lifetime of 20 years with a constant injection fluid temperature of 65 • C and assuming a borehole skin zone radius of 1 m.

Slip tendency and Coulomb stresses
Once the in situ stress tensor and fault's positioning in the subsurface are known, acting shear, , and normal, n , stresses along any arbitrarily oriented planar feature can be resolved (Jaeger et al. 2007). In accordance with the Amonton's law, the slip tendency, T s , defined as a ratio between and a difference between n and pore pressure,P p , can be assessed (Morris et al. 1996) where, C 0 is the cohesion of a fault. The seismicity will be likely induced on a pre-existing zone of weakness in the subsurface, once the ratio of shear to effective normal stress, ′ n , exceeds its frictional resistance, (e.g., Healy et al. 1970) We assume that the model area is (i) tectonically quiescent and (ii) located in the so-called 'low-stress' environment, where the nucleation of a potential seismic event will be restricted to the size of the fault segment that is considered to be critically stressed (2) T s ≥ .

Page 5 of 16 38
Vol.: (0123456789) (i.e., where T s ≥ ). This is contrary to a case of a 'high-stress' environment, where once a seismic event nucleates, it could eventually propagate over the entirety of a fault plane (Wassing et al. 2021). Our assumption of a 'low-stress' environment was later confirmed by a low initial slip tendency of the simulated faults based on the regional 3D geomechanical model (Kruszewski et al. 2022). By summing up the critically stressed fault segments, the size of the reactivated (or ruptured) fault area, A r , can be resolved. It is worth mentioning, that due to the simplified friction model utilized in this study, it is not known when, during the lifetime of a geothermal system, given perturbed fault segment will actually slip and if it will be a coseismic or an aseismic slip.

Seismic event magnitude
Having computed the reactivated area on a fault plane allows estimating theoretical local moment magnitude, M L , due to shear slip from locationspecific empirical correlations. Firstly, we assume an approach by Abe (1975), based on which seismic moment, M 0 , can be straightforwardly determined from A r alone where, M 0 is expressed in dyne cm and A r in km 2 . Having estimated M 0 , correlations by Bischoff et al. (2010), developed from S-wave spectra from induced events with 0 ≤ M L ≤ 2.4 located in the eastern Ruhr region, were applied to compute M L of induced seismic events where, M 0 is in this case expressed in Nm. The applicability of empirical correlation by Bischoff et al. (2010), was assumed in this study as a boundary between acceptable (M L < 2.4) and unacceptable ( M L ≥ 2.4) seismic event magnitude connected to geothermal operations. This criterion can be considered as significantly less conservative than assumption of M L of 2.1 as an upper threshold for the stimulation operations in Finland (Ader et al. 2020) after which immediate stop of stimulation activities was necessary or M L of 1.7 which have already lead to a closure of the deep geothermal project in the eastern Netherlands (Vörös and Baisch 2019). In accordance with the German DIN 4150 standard, damages to the residential buildings can be ruled out for the maximum ground velocities of 5 mm/s, which accordingly to Bischoff et al. (2010) could be expected already for M L between 1.5 and 2.3.

Model geometry and discretization
For the investigation of fluid-induced fault reactivation and induced seismicity risk, and based on the initial slip, T s , and dilatation, T d (Ferrill et al. 1999), tendencies of major faults in the static regional 3D geomechanical model (Kruszewski et al. 2022), an area for a theoretical deep geothermal system located in the southern Ruhr region was selected. The model has a size of 10 km by 10 km by 10 km and consists of five geological layers and three fault zones. The horizon geometry was based on the CICM (Coal Inventory Calculation Model) (Juch et al. 1994;Juch 1997), whereas fault plane geometries were reconstructed from geological cross-sections and surface fault trace maps. The two wells, which are located 1.8 km from one another (measured from the bottom of open hole sections of both wells), were assumed to be drilled to a final depth of 5.5 km with a total length of an open hole section of 700 m, surrounded by Devonian carbonates, and a cased section between surface and a depth of 4.8 km . Such a conservative inter-well distance was selected to exclude a potential risk of a premature thermal breakthrough, as well as thermal short circuit. The directional well A is penetrating the N-S-striking and vertically-dipping normal fault at approximately 5.3 km depth. The vertical well B is located 1.6 km from the N-S-striking fault and 2.0 km from the closest NE-SW-striking reverse fault. All included faults extend from the model surface to a depth of 6 km (Fig. 1). The overburden layers consisting of Cretaceous and Carboniferous units are assumed as low-permeability caprocks with constant permeability, whereas the Devonian layer is assumed as a reservoir layer with an improved permeability. The mechanical, hydraulic, and thermal properties of rock layers and faults (all properties were held constant along faults) are listed in Table 2. As the actual lithological content of Devonian layers remains unknown, we assume two extreme cases of thermal expansion coefficient of the reservoir rock, . The highest value is based on the assumption of high clay and quarzite contents, influence of the temperature and bedding at the reservoir level (Johnson and Parsons 1944;Feng et al. 2020), whereas the lowest value assumes homogenous and pure carbonate rock formation (Robertson 1988). The latter is expected to be rather an unusual case in situ, whereas the real value of the carbonate rock should lie somewhere between these two extremes. A constant of 0.65 (Zimmerman 1991) and of 0.6 (Verberne et al. 2010;Zoback and Kohli 2019) were assumed across the model volume. We regard the of the reservoir of 0.6 as an optimistic value, which, in reality, could be as low as 0.3, if faults would have been filled with clay materials (Ikari et al. 2009).
The model volume was discretized into 271 100 tetrahedral elements with the average element quality

Boundary conditions
The initial stress state with average total principal stress magnitudes and total stress ratio, K-value, extracted from a regional 3D geomechanical model (Kruszewski et al. 2022) is shown in Fig. 1 and presents strike-slip regime at reservoir depth. Under the assumption of no cohesion, the average T s of the simulated segment of the N-S-striking fault, as shown in Fig. 1, is 0.22 with an average, T d of 0.92. The NE-SW-striking reverse faults have slightly higher average T s of 0.27 and much lower T d of 0.26. T s of both fault systems are much lower than the frictional threshold of 0.6 and, therefore, it can be assumed that all simulated faults are initially not critically stressed in the prevailing stress state configuration.
The bottom and side model boundaries are free to move in-plane and are fixed in the out-of-plane direction, whereas the upper surface is free to move in both directions. Faults, represented as 2D planes of weakness for computational efficiency, are assigned constant hydraulic aperture of 0.01 m and regarded as thin elastic layers with stiffness properties being lower than the surrounding rock matrix. We impose hydraulic gradient with linearly increasing P p at the sides of the model and a fixed atmospheric pressure of 0.1 MPa at the surface layer. The model bottom boundary at 10 km was defined as impermeable. The open hole sections of both wells, were represented by line sources, where fluid injection and production was defined with mass flow rates, ̇q . The vertical boundaries were defined as open boundaries for the heat transfer, where the heat source is caused by an injection well. Both production and injection wells are considered as 1D element for improved computational efficiency. The initial temperature is assumed to follow a geothermal gradient of 33 • C/km, which yields a reservoir temperature of approximately 175 • C at depth of 5 km (assuming the surface temperature of 10 • C). We account for the change of reinjection fluid, assumed here to be pure water, properties with temperature with an exclusion of water compressibility which we assume to be constant and amount to 4 × 10 −10 1∕Pa . Temperature-dependent re-injection fluid properties include density, dynamic viscosity, specific heat capacity, thermal conductivity, and thermal expansion coefficient.

Results and discussion
Numerical simulation results for different case scenarios against time, averaged over a quadratic area on a N-S-striking fault plane with a side length of 500 m and well A located in the middle of it, are presented in Figs. 2 and 3. Based on Fig. 2a, b, it can be observed that an injection into fault results in a significant -Poisson's ratio, k res -matrix permeability; -porosity;thermal conductivity; -coefficient of thermal expansion; C pheat capacity at constant pressure); A Hesemann (1965)   increase of Δ , in a worst case reaching 13 MPa, and a decrease of Δ � n , in a worst case reaching 25 MPa, after 20 years of geothermal operations. A higher Δ increase was observed for non-isothermal cases with high ̇q and . Both k res and k f proved to have insignificantly small effects on Δ , whereas only a minute Δ increase was observed in the isothermal scenario. For the case of Δ � n , its biggest decrease was observed for non-isothermal cases with higher ̇q and . As in the case of Δ , k f proved to have insignificantly small effect on Δ � n , whereas lower k res shown only slightly smaller Δ � n values than in the case of higher k res . The Δ � n for an isothermal case remains at approximately 1.7 MPa throughout the whole simulation time and is substantially lower that for the non-isothermal cases. As presented in Fig. 2c, the highest ΔT s increase is observed for non-isothermal cases with high ̇q and . Both k res and k f , proved to have negligible effect on T s . Only a minute and rather an irrelevant increase of ΔT s is seen for the isothermal case throughout 20 years of operations. Thermal reservoir properties proved, therefore, to have a substantial role in the estimation of an onset on a fault shear slip, where higher value provides higher Δ , Δ � n , and ΔT s .
Looking at Fig. 2d, it can be observed that the highest change of overpressure at the injection well, ΔP p , was observed for the case with low k res and amounts to 40 MPa after 20 years of geothermal operations. In comparison, significantly lower ΔP p of only 4 MPa was registered for the case with higher k res and is connected to the faster fluid diffusion into the model volume than in the case of low k res . No significant influence of k f on ΔP p was observed. For the case with ̇q of 50 l/s, ΔP p amounted to 1.3 MPa after 20 years from the start of operations, whereas in the isothermal case, increase of ΔP p , amounted to 12 MPa.
In cases where a production well is crossing the N-S-striking fault, as presented in Fig. 3, only a very small Δ increase, in a worst case reaching approximately 0.4 MPa, was registered (Fig. 3a). Conversely to the case with an injection well penetrating a N-S-striking fault, having a production well at a N-S-striking fault results in an increase of Δ � n , which in the most extreme case of low k res approaches 2.7 MPa (Fig. 3b). Such a perturbation of the initial stress condition leads to only a small and rather irrelevant ΔT s increase, which will not lead to any of  Bischoff et al. (2010). In cases 1.6, 1.7 as well as between 2.1 and 2.7, no A r and, therefore, no seismicity was recorded 38 Page 10 of 16 Vol:. (1234567890) the fault segments reaching their critical threshold (Fig. 3c). For the isothermal case, fault stabilization, i.e., ΔT s decrease, along the whole simulation time was observed. It can be, thus, said that in the case of a production well penetrating the N-S-striking fault, reactivation of the previously aseismic N-S-striking fault is not expected.
The reactivation of the N-S-striking fault, for cases with an injection well penetrating the fault, starts in the first months of geothermal operations (Fig. 4). For cases with ̇q of 150 l/s and C 0 of 0.2 MPa (e.g., Wyllie and Norrish 1996), which is assumed to be cohesion of a clay-filled fault, A r increases linearly and reaches approximately 0.2 km 2 at the end of fluid circulation. For the case with ̇q of 50 l/s, A r increases also linearly until the end of geothermal operations, but approximately factor 3 smaller than for the case with ̇q of 150 l/s. This indicates, that lower ̇q significantly reduces the amount of fault reactivation and, simultaneously, seismic hazard. In cases with C 0 of 20 MPa, which is assumed to be cohesion of a calcite-filled Change of Δ on a N-Sstriking fault plane between the start and 20 years of geothermal doublet operation for case 1.2 (non-isothermal); c Change of Δ on a N-S-striking fault plane between the start and 20 years of geothermal doublet operation for case 1.7 (nonisothermal); d Change of normal effective stresses, Δ � n , on a N-S-striking fault plane between the start and 20 years of geothermal doublet operation for case 1.6 (isothermal); e Change of Δ � n on a N-S-striking fault plane between the start and 20 years of geothermal doublet operation for case 1.2 (nonisothermal); f Change of Δ � n on a N-S-striking fault plane between the start and 20 years of geothermal doublet operation for case 1.7 (non-isothermal); g Change of slip potential, ΔT s , on a N-S-striking fault plane between the start and 20 years of geothermal doublet operation for case 1.6 (isothermal); h) Change of ΔT s on a N-S-striking fault plane between the start and 20 years of geothermal doublet operation for case 1.2 (non-isothermal); i) Change of ΔT s on a N-S-striking fault plane between the start and 20 years of geothermal doublet operation for case 1.7 (non-isothermal) Page 11 of 16 38 Vol.: (0123456789) fault (e.g., Li et al. 2020), A r does not exceed 0.03 km 2 and is the highest between the 5th and 13th year of geothermal production. A slight decrease of A r is observed in the last 5 years of geothermal production for cases with C 0 of 20 MPa. Based on estimated M L , it can be seen that for case scenarios with ̇q of 150 l/s and C 0 of 0.2 MPa, relatively quick, i.e., after 1.8 years, M L exceeds 2.4, whereas for the case with ̇q of 50 l/s this time amounts to 4 years. In cases with C 0 of 20 MPa, the time when M L exceeds 2.4 amounts to 2.9 years for ̇q of 150 l/s and 7 years for ̇q of 50 l/s. It is also observed that lower k res only slightly increases A r , whereas k f , seem to have no major influence for cases with C 0 of 0.2 MPa. For cases with C 0 of 20 MPa, on the other hand, lower k res indicates lower A r . The case with ̇q of 50 l/s presents more steady and mellow increase of A r , peaking around the 13th year of geothermal production, in comparison to the more sharp increase of A r for cases with ̇q of 150 l/s, peaking around the 5th year of operation. For the isothermal case, case with low , and for all cases with a production well penetrating a N-S-striking fault, stresses exerted on the N-S-striking fault stay below the frictional limit and, therefore, no fault reactivation is to be observed there.

Thermal effects
Comparing the non-isothermal cases with the isothermal case, as seen in Figs. 2 and 3, it is clear that thermally-induced stresses play a dominant role, in the long-term, during fault reactivation, as they greatly increase the likelihood of induced seismicity with pore pressure effects playing rather a marginal role during fault reactivation. Similar conclusions were drawn by Jacquey et al. (2015), Parisio et al. (2019), andWassing et al. (2021). This can be proven by comparing results between non-isothermal and isothermal cases, where stresses after 20 years of geothermal operations exerted on the N-S-striking fault plane never exceed the frictional threshold, and the fault could be considered stable throughout the simulation time of 20 years (Fig. 5). In the nonisothermal cases, the same fault undergoes a linear increase of reactivated area (Fig. 4). This indicates, that potentially large seismic events could be triggered by reservoir cooling during fluid circulation. Such a phenomenon can be explained by the propagation of the thermal plume (i.e., a rock volume where temperature decrease between an initial state and the end of geothermal operations is observed) from an injection well deeper into the rock volume. The extent of a thermal plume during different time intervals for case 1.2 is presented in Fig. 6a. As shown there, fluid re-injection cools down and contracts the rock volume around the injection well causing negative volumetric strains (i.e., thermal contraction). As, with time, thermal plume propagates deeper into the rock, it decreases Δ � n and increases Δ , which enables the favourably-oriented fault to exceed its frictional limit (Fig. 5h). Thermal effects during the early years of geothermal doublet operation are negligible, as temperature change of the rock matrix remains small. Additionally, an increase of Δ � n , and, therefore, fault stabilization is seen around the cooled region at the injection well causing a "halo" effect ( Fig. 5e,f).
The size of the perturbed fault segments for cases 1.2, 1.6, and 1.7 can be seen in Fig. 5. As it can be deduced from Fig. 5e and f, the thermal plume tends to sink in below the injection well due to the higher density of the colder re-injection fluid with the production well remaining in isothermal conditions (as shown in Fig. 6c). The radius of the thermal plume, which for the purpose of this study was assumed as a shortest horizontal distance between an injection well and the farthest reaching temperature perturbation within the model volume, after 20 years of geothermal operation, for the case with highest k res and k f , amounts to approximately 500 m ( Fig. 6a and c; shortest inter-well distance marked with L) with a pressure plume having approximately twice as large radius (Fig. 6c). It is worth mentioning, that the size of thermal plume will be a function of time, therefore, for cases with longer system lifecycles, larger thermal plume radii should be expected. Based on the radius of the generated thermal plume, a minimum safe distance between an injection well and a fault that would not lead to its reactivation could be proposed. This allowable distance for our case study, which when not obeyed could lead to an initiation of fault reactivation processes, will, therefore, amount to 500 m. The shape of the thermal plume is, however, highly influenced by the assumption of a homogeneous and isotropic reservoir extending from the bottom of the Lower Carboniferous layers to a depth of 10 km and is motivated by the lack of geological data from the study region. In the case of a reservoir being constrained by an impermeable sealing layer from below, the radius of the thermal plume may be different than simulated in this study, which will considerably affect the amount of fault reactivation and the extent of the thermal plume. An assumption of isotropic permeability is also a substantial simplification in this model. In reality, dense natural fracture networks could highly influence fluid flow in the study area and create significant permeability anisotropies.
Once stresses around the injection well turn tensile, hydrofracturing (i.e., propagation of smaller scale tensile fractures) may occur. This process can be considered positive, as it will result in an increase of near-borehole permeability, but will also lead to nucleation of small magnitude earthquakes, which could be deemed as undesirable. In contrary to the relatively small magnitude tensile earthquakes, larger, shear-dominant, seismic events will be delayed and their timing will depend on the extent of a thermal plume (Fig. 4). The thermal plume could potentially propagate also after fluid injection is halted, which will mainly dependent on the local hydro-geological conditions (i.e., low permeability reservoir will be more prone to post-injection effects than a conductive aquifer). Such delayed post-injection trailing effects have been observed in a few geothermal fields worldwide (e.g., Lee et al. 2019;Mignan et al. 2015;Vörös and Baisch 2019). It is worth mentioning, that as of today, and based on authors knowledge, none of the induced seismic events recorded in geothermal fields worldwide could be explained by to the influence of thermal contraction on the fault slip. Most of the seismic events, known in the literature, are believed to be caused by, mainly short-term, pressure effects. It remains, however, unclear if similar simulations, as carried out in this study, were undertaken and if thermal effects were actually regarded as important.
As proven above, thermally-induced perturbations of the initial stress state due to the reservoir quenching will play a dominant role in the long-term stability of deep geothermal systems in the Ruhr region. This could imply that (semi-)continuous and realtime temperature field monitoring, complementing seismic monitoring, could be considered as a crucial element of re-risking future geothermal projects. One of the methods, although limited to a narrow area in the close vicinity of a borehole, could include temperature monitoring utilizing fiber optic cable installed in a borehole (Reinsch 2012;Schölderle et al. 2021). Another method could include electrical resistivity tomography methods, which enable controlling the temperature distribution within the reservoir and are Page 13 of 16 38 Vol.: (0123456789) not limited to a near-borehole area. Such method was successfully field-tested e.g., in shallow geothermal wells in eastern Belgium (Hermans et al. 2014).

Model uncertainties and way forward
Uncertainties introduced to the model result from (i) model geometry, (ii) rock mass and (iii) fault zone properties, (iv) initial stress state, as well as (vi) basic model assumptions and boundary conditions.
One of the core model unknowns remains the uncertainty of the existence of Devonian layers underlying the Carboniferous strata across the Ruhr region as well as its thickness. The only indication of such comes from interpretations of the DEKORP 2-N seismic line (DEKORP 1990). Besides, fault zone geometry at greater depths remains uncertain. The presence and architecture of implemented fault zones may, in reality, differ from the ones constructed in the model, which is especially the case for the deeper subsurface. The assumptions made in this study shall be validated with 2-to 3-D seismic surveys and results from deep exploration drilling campaigns in the region. As it was proven in previous case studies, slip and resultant seismic activity may be induced on a previously unknown critically-stressed fault (e.g., Lee et al. 2019), which could be located below the depths of economical interests. To exclude such occurrences and de-risk future deep geothermal projects in the area, seismic campaigns prior to drilling and seismic monitoring during and after drilling operations could be implemented.
It is worth mentioning that uncertainty will be also connected to the properties of host rocks and faults. Especially, of 0.6 and of 0.65, which were assumed to be constant along the model volume, are considerable simplifications. Frictional properties may change from fault to fault and along a trace of one particular fault and may differ in different lithological units. To reduce the uncertainty of the aforementioned parameters, more detailed laboratory measurements in reservoir-specific conditions are sought.
Fault cohesion, as proved in Fig. 4, could be considered as one of the main parameters controlling the slip on fault planes, where higher cohesion may lead to potentially delayed and weaker earthquakes. Values assumed in this study are merely two extreme cases of more ductile clay-filled and more brittle calcitefilled faults assumed from literature sources. The actual fault properties for the region remain unknown. Therefore, future research in the induced seismicity risk in the Ruhr region, shall in detail investigate fault properties and their change along different fault orientation as well as along faults in situ.
The thermal properties of the reservoir rock can be regarded as one of the most crucial parameters in assessing thermal stress contribution and, therefore, the onset on fault slip. In this point in time, and until more laboratory data will be produced, the actual value of remains uncertain.
Although the developed model is location-specific, it can still provide insights into how a geothermal system in a low-permeability sedimentary setting could work and if there is a risk of fault reactivation connected to it. We recognise the simplicity of the model and simulation approach. Due to the limited availability of data and associated uncertainties, richer models with more assumptions will undoubtedly face a risk of overparameterization and will introduce additional sources of uncertainty. Once results from exploration drilling, seismic campaigns, and laboratory studies will be produced, model parametrization could be significantly improved by e.g., results from hydraulic tests, more complex friction models, or chemical coupling.

Conclusions
The results from this study highlight the complexity of fault reactivation processes and indicate that thermally-induced stresses, in the long-term, play a dominant role during fault reactivation and greatly increase the likelihood for induced seismicity. We observe that large seismic events, and their timing, will depend primarily on the quenching of the reservoir. To assure the safety of a future geothermal system in the Ruhr area, an injection well should not be placed on a fault, as local moment magnitudes exceeding 2.4, could be induced in the first years of geothermal doublet operation, even at relatively low operational flow rates. On the other hand, the placement of a production well will not cause fault reactivation. As thermally-induced stresses proved to have dominant effects on the long-term seismic risk, we propose to account for the radius of a thermal plume, generated during the expected lifetime of a geothermal system, as the minimum allowable distance between a fault and an injection well. This distance, in our location-specific study for an assumed system lifetime of 20 years, amounted to approximately 500 m Observations made in this study could allow de-risking of future geothermal projects in the region and enable better and more informed well placement and operational designs. We believe that similar forecasting models, accounting for subsurface 3D geometries, coupled physics, and process dynamics will allow for risk mitigation in future geothermal projects not only in the Ruhr region but also in other geothermal reservoirs worldwide. Additionally to the simulation efforts, (semi-) continuous temperature field monitoring, not exclusively limited to the borehole vicinity, could significantly help in de-risking future deep geothermal projects. We shall also emphasize that more field (primarily in situ stress) and laboratory (especially fault's frictional properties and thermal properties) data, under reservoir-specific conditions, as well as seismic monitoring prior, during, and after drilling operations, are of utmost importance to better understand governing mechanisms of fault reactivation and connected seismic risk. In accordance to the range of assumptions made in this analysis, a more detailed study incorporating results from seismic campaigns and exploratory drilling is required to improve the numerical simulation and reduce model uncertainties.