# Poroelastic stress relaxation, slip stress transfer and friction weakening controlled post-injection seismicity at the Basel Enhanced Geothermal System

Auregan Boyet, Silvia de Simone, Shemin Ge, Víctor Vilarrasa

► To cite this version:

Auregan Boyet, Silvia de Simone, Shemin Ge, Víctor Vilarrasa. Poroelastic stress relaxation, slip stress transfer and friction weakening controlled post-injection seismicity at the Basel Enhanced Geothermal System. Communications Earth & Environment, 2023, 4 (1), pp.104. 10.1038/s43247-023-00764-y. hal-04090065

HAL Id: hal-04090065

https://hal.science/hal-04090065

Submitted on 5 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Poroelastic stress relaxation, slip stress transfer and friction weakening controlled post-injection seismicity at the Basel Enhanced Geothermal System

Auregan Boyet ID 1,2,3✉, Silvia De Simone⁴, Shemin Ge ID ⁵ & Víctor Vilarrasa ID <fnref n="1" />,²,³

Induced seismicity is a limiting factor for the development of Enhanced Geothermal Systems (EGS). Its causal mechanisms are not fully understood, especially those of post-injection seismicity. To better understand the mechanisms that induced seismicity in the controversial case of the Basel EGS (Switzerland), we perform coupled hydro-mechanical simulation of the plastic response of a discrete pre-existing fault network built on the basis of the monitored seismicity. Simulation results show that the faults located in the vicinity of the injection well fail during injection mainly triggered by pore pressure buildup. Poroelastic stressing, which may be stabilizing or destabilizing depending on the fault orientation, reaches further than pressure diffusion, having a greater effect on distant faults. After injection stops, poroelastic stress relaxation leads to the immediate rupture of previously stabilized faults. Shear-slip stress transfer, which also contributes to post-injection reactivation of distant faults, is enhanced in faults with slip-induced friction weakening.

<footnote n="1">¹Institute of Environmental Assessment and Water Research, Spanish National Research Council (IDAEA-CSIC), Barcelona, Spain. ²Associated Unit: Hydrogeology Group UPC - CSIC, Barcelona, Spain. ³Global Change Research Group (GCRG), IMEDEA, CSIC-UIB, Esporles, Spain. ⁴Univ Rennes, CNRS, Géosciences Rennes, UMR 6118, 35000 Rennes, France. ⁵University of Colorado Boulder, Boulder, CO, USA. ✉email: auregan.boyet@csic.es</footnote>
Induced seismicity represents one of the main obstacles to the development of geothermal energy, which is a key low-carbon technology to reach mid-century net-zero carbon emission targets¹. In Enhanced Geothermal Systems (EGS), a fluid is circulated through a newly created and/or stimulated fracture network carrying heat through fluid circulation to the production well, increasing the generation of electricity. Seismicity of small magnitude ($M<2$) is in general observed, especially during the stimulation phase, in which a massive fluid injection is performed to enhance the permeability of pre-existing fractures such that the flow rates are sufficient for geothermal power production. This injection-induced seismicity has occasionally reached magnitudes large enough to be felt on the surface². Felt induced earthquakes are undesirable not only because they may injure people and damage buildings and infrastructure, but also because they cause a negative effect on public perception that may lead to project cancellation, as occurred at the EGS projects at Basel, Switzerland³, and Pohang, South Korea⁴. For these two cases, an intriguing common characteristic of induced seismicity by EGS stimulation is that the largest earthquakes take place after the stop of injection, when the induced seismicity potential is supposed to decrease because pore pressure drops.

Following these and other cases of poorly-understood induced seismicity⁵,⁶, the last two decades of research activity have extensively discussed the existence of multiple triggering mechanisms, including pore pressure buildup, poromechanical stress changes, and aseismic or seismic slip stress transfer⁷⁻¹⁰. By increasing the pore pressure, the effective normal stress acting on pre-existing fault surfaces, and consequently the shear resistance, reduces and may induce slip on these surfaces¹¹. Despite being the most common causal mechanism, pore pressure buildup may not be the only cause of induced seismicity, and in some cases large increases in pore pressure would be required to reach failure conditions. During hydraulic stimulation, fluid injection alters the pore pressure and temperature of the rock. The low temperature of the injected fluid (compared to the in-situ temperature) progressively cools down the vicinity of the injection well, becoming particularly significant during long-term fluid circulation¹². Poroelastic stresses propagate much ahead of the pressurized region and can trigger seismicity at large distance from the injection well¹³,¹⁴. Because these stresses are anisotropic, they can improve or worsen the mechanical stability of pre-existing faults depending on fault orientation¹⁵⁻¹⁷. Moreover, fluids flowing along preferential pathways have an anisotropic impact on the local stress tensor, leading to more pronounced anisotropy in the poroelastic stress redistribution¹⁸.

Pore pressure diffusion is also considered as a mechanism of post-injection induced seismicity¹⁷,¹⁹, as pore pressure continues to propagate in the reservoir after the stop of injection. Earthquake interaction is another potential mechanism²⁰. The stress variation caused by shear slip activation near the injection well may promote failure on nearby faults. This shear-slip stress transfer may be the result of both seismic and aseismic slip. Aseismic slip generated by pore pressure diffusion may precede seismic slip²¹,²² and its associated shear-slip stress transfer may seismically reactivate nearby or distant faults, including faults placed outside the pressurized area after injection stops²³, like in the case of the underground gas storage of Castor, Spain¹⁰. Yet, the cumulative stress transfer due to small events near the injection well is often insufficient to explain the reactivation of nearby or distant faults²⁴⁻²⁶. In some cases, only the combination of slip-induced stress redistribution, poro-thermo-elastic effects, pressure diffusion and fault weakening may explain the seismicity observed in the post-injection period¹⁸,²⁷. Additionally, fluid injection impacts the geomechanical properties of the rock, especially in the fault zones which concentrate preferential flow

paths. The activation of both seismic and aseismic shear slips degrades the fault frictional properties because part of the asperities at the fracture walls are deteriorated²⁸. This slip weakening further reduces the frictional resistance, promoting additional fault slip. Models taking into account this slip weakening allow to reproduce seismicity with larger events than models without friction weakening²⁹. The above-mentioned processes are highly influenced by the local stress field, fracture distribution and connection, rock mechanical properties, hydrologic factors and historical natural seismicity, which makes it a complex phenomenon to understand in detail².

Despite the progress in deepening the understanding of these processes, the ultimate causes of the high-magnitude post-injection seismicity remain not fully understood. In particular, the causes of the post-injection seismicity at the Deep Heat Mining Project at Basel (Switzerland) are not clear. During hydraulic stimulation in December 2006, event magnitudes up to $M_L2.6$ were recorded. Operations were then stopped, but an event of $M_L3.4$ ($M_W2.95$) occurred 5 hours after shut-in in the stimulation well³. Subsequently, the well was bled off, i.e., the wellhead was opened and hydrostatic pressure was imposed along the well. This felt post-injection induced seismicity led to the abandonment of the project. Numerous conjectures and studies have been developed since then to explain the observed seismic response, most of them focusing on pressure and stress redistribution consequent to fluid injection. Pressure diffusion has been shown to be the causal mechanism for part of the seismicity occurred during injection at Basel. Mukuhira et al.³⁰ compared the injection-induced pressure build-up, estimated by considering a homogeneous domain, with the critical pressure required for fault failure³¹. Similarly, Terakawa et al.³² and Terakawa³³ used the observed seismicity and critical pressure considerations to map the overpressure distribution in three dimensions, which was obtained exclusively by invoking the existence of preferential diffusion through pre-existing faults. They found that the inferred pressure is consistent with the wellhead pressure history. However, the values of overpressure necessary to explain seismicity are often unrealistically high and they only partially explain the co-injection seismicity. Slip-induced stress redistribution as a triggering mechanism at Basel EGS has been highlighted by showing that the seismicity rate is correlated with the interactions between seismic events²⁴,³⁴,³⁵. It has been suggested that co-injection induced seismicity was triggered by pore pressure diffusion¹⁹, while stress redistribution dominates in the post-injection induced seismicity³⁴. Andrés et al.³⁶ proposed a conceptual poroelastic 3D-model of fault reactivation to evaluate the potential causal mechanisms of induced seismicity at Basel. Although this was the first attempt to reproduce the seismicity at Basel acknowledging coupled hydro-mechanical (HM) effects, Andrés et al.³⁶ were only able to reproduce the temporal evolution of reactivation of a conceptual single fault plane. Moreover, the respective role of the different mechanisms and their combination were not deeply analyzed. After more than 15 years from the halting of the Basel EGS operations, a description of the main triggering mechanisms at Basel to reproduce the spatio-temporal observation of seismicity is still missing.

In this study, we simultaneously simulate pore pressure diffusion, poroelastic stress redistribution and shear slip-induced stress interactions, to understand their impacts on induced seismicity during and after injection at Basel. We build an explicit faulting model which is based on the seismic observations at Basel, and we solve the coupled HM problem associated with the hydraulic stimulation. We aim at identifying the mechanisms responsible of failure at different reservoir locations. We analyze the role of pressure diffusion, poromechanical stressing, stress variation due to shear slip activation and friction weakening
<figure><img src="image_3.png" /><figcaption>**Fig. 1 Geometry and conditions of the model.** a Plan view of the location of seismic events at Basel, sorted by clusters according to Deichmann et al.³⁷. Colors correspond to different clusters as indicated by the numbers. b Model setup showing geometry, fault network and boundary conditions. The central black dot represents the injection well and the gray area the damaged zone. Initial values of principal stresses and pore pressure are indicated, along with the boundary conditions of constant pressure prescribed at the outer boundaries (represented by the crossed circle) and no displacement perpendicular to the boundaries (represented by the triangles). Note the b represents the central region of the model containing the faults, but that boundaries are actually at 950 m away from the well. The direction of shear slip is also indicated for each fault. c Pressure evolution in time at the well location, according to Häring et al.³ (see inset in Figure b).</figcaption></figure>

consequent to shear slip activation, on the activation of post-injection seismic activity.

## Modeling approach

**Conceptual model.** We reproduce the setting of Basel EGS by building a numerical model of fluid injection into a fault network embedded in the crystalline basement. The network geometry is based on the interpretation of the seismic event clustering performed by Deichmann et al.³⁷. The faults are explicitly represented in our model, reproducing the pre-existing faults in the reservoir that were seismically activated (Fig. 1b). The faults may reactivate in shear mode as a result of water injection, but we disregard fracture creation or propagation during the hydraulic stimulation. The model domain consists of a plane strain 2D horizontal section intersected by a vertical injection well. This simplification of the reservoir is supported by the long open-hole section of the injection well (more than 300 m) and by the fact that most of the monitored events exhibits a focal mechanism with strike-slip movement and vertical dip³,³⁷,³⁸. These conditions imply that poroelastic stress and deformation are constant along the depth and that plane-strain conditions can be assumed on a horizontal section. We solve the fully coupled HM problem to estimate the pore pressure and stress variations, as well as fault reactivation, both during injection and after the stop of injection, including the bleed-off period after shut-in.

Pore pressure diffusion, poroelastic stressing and shear-slip stress transfer are analyzed as potential causal mechanisms of induced seismicity. Initially, we assess fault stability with purely hydraulic considerations, i.e., neglecting poromechanical effects. The direct effects of pore pressure diffusion on the stability of the faults are determined by applying the critical pressure theory³¹. For each fault, the critical pressure, $P_c$, corresponds to the maximum pressure value that the fault can sustain before reaching failure conditions, as expressed by

$$P_c = \sigma_n - \frac{\tau}{\mu}, \qquad (1)$$

where $\sigma_n$ and $\tau$ are, respectively, the normal and shear stresses acting on the fault and $\mu$ is the friction coefficient. The normal

and shear stresses define the fault slip tendency, and are calculated from the orientation of the fault and the regional stress field, which is assumed to remain constant during injection. The friction coefficient $\mu$ is an intrinsic property of the fault, which typically takes values around 0.6 in crystalline rocks³⁹. At Basel, since there are data of the onset of induced seismicity of each fault, i.e., fault activation, we can calibrate which is the actual friction coefficient $\mu_\varphi$ (the friction coefficient being the tangent of the frictional angle $\varphi$) for each fault. The obtained critical pressure $P_{c,\varphi}$ from the calibrated friction angles results in values that are more coherent with the pressure buildup and observed induced seismicity (Table 1). We compare the numerically simulated pressure variations in the vicinity of the faults with their critical pressure values to identify the direct impact of pore pressure diffusion on fault reactivation.

Following this simplified analysis, we analyze the fully coupled HM stress variation and consequent failure conditions in two modeling scenarios with elastic and visco-plastic fault mechanical behavior, respectively (for additional details see the Methods section). In the elastic scenario, fault deformation is reversible, and the shear displacement is small with no permanent static stress transfer and inversely proportional to the fault stiffness. In the visco-plastic scenario, faults respond to a Mohr-Coulomb failure criterion. When failure conditions are met, irreversible and abrupt shear slip occurs, with consequent irreversible stress redistribution, i.e., shear-slip stress transfer. The comparison of the two scenarios allows us to distinguish the effects of poroelastic stress from those of shear-slip stress transfer.

**Numerical model setup.** The 2D geometry represents a horizontal plane of 3.61 km² located at approximately 4630-m deep, coinciding with the injection depth in the crystalline basement at Basel. A set of faults is embedded in the rock matrix, with the fault network derived from the induced seismicity registered in the range from 3750 to 4750-m deep. We make use of the open-data seismic catalog by Deichmann et al.³⁷, which proposes a clustering of events occurring during injection and short-term post-injection stages that is based on their locations and focal
<table><caption>Table 1 Characteristics and estimated parameters for each fault in the numerical model of the Basel EGS project.</caption><thead><tr><th>Fault</th><th>Strike (°)</th><th>Dip (°)</th><th>Rake (°)</th><th>M<sub>w</sub></th><th>φ (°)</th><th>P<sub>c</sub> (MPa)</th><th>P<sub>c</sub> (MPa)</th><th>Net slip (m)</th><th>Slip area (m²)</th></tr></thead><tbody><tr><td>A</td><td>160</td><td>88</td><td>-165</td><td>1.75</td><td>27.5</td><td>54</td><td>48</td><td>0.0087</td><td>1 946</td></tr><tr><td>B</td><td>170</td><td>50</td><td>-50</td><td>2.10</td><td>19.5</td><td>66</td><td>53</td><td>0.0212</td><td>14 236</td></tr><tr><td>C</td><td>160</td><td>88</td><td>-155</td><td>2.95</td><td>24.5</td><td>59</td><td>49</td><td>0.0090</td><td>8 019</td></tr><tr><td>D</td><td>198</td><td>88</td><td>10</td><td>2.24</td><td>29</td><td>50</td><td>46</td><td>0.0124</td><td>846</td></tr><tr><td>E</td><td>206</td><td>71</td><td>-43</td><td>1.91</td><td>29.5</td><td>50</td><td>48</td><td>0.0101</td><td>1 403</td></tr><tr><td>F</td><td>200</td><td>85</td><td>14</td><td>1.75</td><td>28</td><td>52</td><td>47</td><td>0.0035</td><td>5 308</td></tr><tr><td>G</td><td>204</td><td>48</td><td>-22</td><td>1.99</td><td>30.5</td><td>50</td><td>49</td><td>0.0036</td><td>8 878</td></tr></tbody></table>

$M_w$ is the seismic event magnitude of the largest event of the cluster group associated with the corresponding fault. Friction angle $φ$ is calibrated to reproduce fault reactivation occurrence. Critical pressure $P_c$ is calculated with Eq. (1), using a friction coefficient of $μ = 0.6$ (corresponding to $φ = 30°$), while critical pressure $P_{c,v}$ using the friction coefficient that leads to failure at the time of the onset of induced seismicity for each fault, i.e., $μ_v = tan(φ)$. Net slip is calculated with Eq. (3).

<table><caption>Table 2 Hydro-mechanical properties of the materials composing the numerical model.</caption><thead><tr><th></th><th>Rock matrix</th><th>Damaged zone</th><th>Faults</th></tr></thead><tbody><tr><td>Young's modulus, E (GPa)</td><td>52</td><td>50</td><td>43</td></tr><tr><td>Poisson's ratio, ν (-)</td><td>0.25</td><td>0.25</td><td>0.25</td></tr><tr><td>Permeability, k (m²)</td><td>7.50·10<sup>-17</sup></td><td>1.00·10<sup>-14</sup></td><td>2.30·10<sup>-13</sup></td></tr><tr><td>Porosity, φ (-)</td><td>0.01</td><td>0.01</td><td>0.1</td></tr></tbody></table>

mechanisms (Fig. 1a). We translate these clusters into fault planes that constitute the fault network of the domain. From the 11 clusters proposed by Deichmann et al.³⁷, we simplify the fault network into an equivalent network composed of 7 pre-existing faults (Fig. 1b, Table 1 and Supplementary Table S1) by considering the microseismicity focal mechanisms, location, timing and magnitude (see Methods section).

The granitic rock is described as a porous isotropic material with linear elastic mechanical behavior. A damaged zone, with diameter of 20 meters, with a higher permeability and a lower stiffness surrounds the well to mimic an intensely altered region (Table 2). Faults are represented by continuous material elements with a thickness of 10 centimeters. For all materials, the values assigned to the hydraulic and mechanical parameters are the results of literature reviews and model calibration against field data of injection rate, wellhead pressure and seismicity activation³,³⁷. In particular, we calibrate the permeability of the rock matrix for a range from $10^{-18}$ to $10^{-16}$ m² and choose the value with the best fit between the monitored induced seismicity and pore pressure diffusion. The fault intrinsic permeability varies with fracture aperture following the cubic law as faults are deformed, with an initial value that is several orders of magnitude larger than the rock matrix permeability. Large permeability increase at the fault reactivation is important to model pressure diffusion through the domain⁴⁰. The specific friction coefficient for each fault is calibrated according to the slip tendency analysis, which is performed using the regional stress field³ and the focal mechanisms (Table 1) to match the fault reactivation with the monitored seismicity. In addition, a slip-weakening of the friction angle of 5° is assigned to fault C.

We set the initial conditions as defined by Häring et al.³, with the maximum principal stress aligned with the y-axis after a rotation of 36°, $S_{Hmax}$ having a strike of 144°N, and values of stresses and hydrostatic pressure as indicated in Fig. 1b. Constant pressure boundaries are set at 45 MPa, to match the hydrostatic pressure. Yet, the nature of the hydraulic boundary condition does not affect the results if pore pressure diffusion does not approach the domain boundaries during the simulated time. Normal displacement perpendicular to the boundary is

constrained to zero on all boundaries of the domain. The temperature at the depth of the reservoir is 190 °C and the injection fluid is water. We assume iso-thermal conditions, i.e., the injected fluid is in thermal equilibrium with the reservoir. To correctly reproduce the injection-induced overpressure of the real 3D domain into a 2D domain, we do not impose a fluid injection rate at the injection well, but we directly assign pressure variations as reported by Häring et al.³, after smoothing out the oscillations for computing purposes (see Fig. 1c). At day 6 after the start of the hydraulic stimulation, the well pressure is set to 45 MPa to reproduce the bleed-off.

## Results

**Comparison of modeled slip with Basel microseismicity.** We examine the spatial and temporal reactivation of the faults by observing their visco-plastic deviatoric strain ($ε_p$), which is a measure of irreversible shear deformations. The temporal evolution of $ε_p$ at three locations, at the center and close to the extremities of each fault, is analyzed together with the cumulative seismic moment, calculated (Eq. (4)) from the registered events³⁷ at the clusters associated with each fault (Fig. 2). $ε_p$ variations are observed during both co- and post-injection. For most faults, there is a remarkable correspondence between the numerically estimated values of $ε_p$ and the cumulative seismic moment of the cluster related to the fault, which highlights the ability of our model to capture the seismicity observed at Basel EGS.

Faults A and C are activated during the injection period. Fault A undergoes a progressive failure during co- and post-injection, while fault C undergoes two failures during injection. To better understand the role of friction weakening assigned to fault C, we repeat the simulations disabling this option. We observe that fault C fails only once if friction weakening is disabled (dashed lines in Fig. 2), i.e., the second failure does not occur. The effect of not considering friction weakening in fault C on the behavior of the other faults is significant for the multiple subsequent reactivations of faults A and E, located next to fault C (Supplementary Fig. S1). Ruptures occur both during and after injection for faults B and E, while fault F exclusively fails during the post-injection period. Values of visco-plastic strains at faults D and G are not large enough to interpret them as fault failures. We compare the scales of the numerically simulated displacements with the net slip estimated from the recorded magnitude of seismic events (Table 1, see also Methods) to verify the coherence of our model with the seismicity of the Basel EGS. Although it is not straightforward to differentiate between aseismic and seismic slip, the model results are overall temporally and spatially consistent with the monitored seismicity at Basel EGS³⁷. In Section 3.3, we analyze in detail the different patterns of failure to identify the triggering mechanisms.
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Temporal evolution of modeled plastic strain, $\epsilon_p$ (solid lines), at different locations on each fault (represented by different colors as indicated in the top-left plot), and of the cumulative seismic moment $M_o$ (gray solid lines) of the observed seismic events as reported by Deichmann et al.³⁷. At fault C, $\epsilon_p$ for the case without friction coefficient weakening is represented by the dashed lines. The vertical gray line represents the stop of injection at day 6. Note that the scale of the plastic strain is different for each fault.</figcaption></figure>

**Pore pressure diffusion.** Injection-induced pressure build-up spreads radially in the reservoir until reaching faults A and B, which are located around the injection well and alter the radial pressure propagation because they represent preferential pathways for pressure along their directions (Fig. 3a – observe the change of pressure gradient corresponding to fault A on day 2 in Fig. 3c). Pore pressure propagation is affected by the strain-dependent permeability of the faults, which causes an irreversible permeability enhancement up to four orders of magnitude upon failure (see Supplementary Fig. S2). After the stop of injection and subsequent bleed-off on day 6, pore pressure decreases drastically in the reservoir tending to recover the initial hydrostatic value of 45 MPa (Fig. 3b). However, residual fluid overpressure continues diffusing and it may trigger failure during the post-injection phase, as shown by the pressure variations along a cross-section at different times (Fig. 3c). Indeed, the peak of pore pressure propagates farther in the reservoir, as shown at fault C, where pore pressure is higher on day 8 than at the stop of injection on day 6.

To determine the triggering role of pore pressure diffusion, we compare the estimated critical pressure (Table 1) with the pressure reached for each fault (Fig. 4). We observe that, in some cases, the pressure increase is not sufficient to initiate rupture at the time of the observed seismicity. For example, at fault A, the critical pressure is only reached after 2 days of injection, whereas the fault is activated at the start of injection. The extreme situation of this case occurs in faults that are reactivated while the critical pressure is not reached, e.g., faults B, C and E (Fig. 4). In contrast, other faults are reactivated after the critical pressure is exceeded. For instance, the critical pressure is reached at fault F after 4 days of injection, but the fault is reactivated after the stop of injection, at day 7. These observations highlight that triggering mechanisms other than pore pressure buildup control induced seismicity.

**Stress variations.** We focus on five faults and we analyze the evolution of the effective normal stress $\sigma'_n$ and shear stress $\tau$
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Pore pressure evolution in the domain. **a** Contour plot of pressure at the end of injection (day 6), **b** contour plot of pressure in the post-injection stage (day 7), **c** spatial distribution of pressure at different days from day 2 to day 10 along a cross-section extending from the injection well to the bottom-left corner of the domain, which crosses faults A and C, as represented by the dashed white line in **a**. The damaged zone is represented by the gray band in Fig. 3c.</figcaption></figure>

acting on these faults, which allow us to illustrate different rupture patterns (Fig. 4). To facilitate the analysis, we show the variations of stresses with respect to their initial values. The failure potential along a given fault is expressed by the Coulomb Failure Stress (CFS), calculated as⁴¹

$$CFS = |\tau| - \sigma'_{n}\tan(\varphi), \quad (2)$$

where $\varphi$ being the initial friction angle. CFS depends on the stresses acting on the fault plane, and thus, on its orientation. A positive CFS indicates that failure conditions are reached along the fault. Here, we analyze the CFS variations with respect to the initial conditions, $\Delta$CFS, which illustrate the evolution of fault stability, i.e., a fault becomes less stable when $\Delta$CFS is positive. The elastic model indicates the effect of the poroelastic stressing, while the differences between the visco-plastic and the elastic models quantify the effect of the shear-slip stress transfer.

Both faults A and B undergo a significant increase of pore pressure (5–10 MPa) due to their proximity to the injection well. The pore pressure buildup leads to the decrease in the effective normal stress in both the elastic and visco-plastic scenarios. Concerning the shear stress variation, it is interesting to note that in the elastic scenario, $\Delta\tau$ is reduced during injection in fault A, whereas it increases in fault B. However, the opposite occurs when considering visco-plasticity of the faults as a result of the shear stress drop that takes places during slip. Before failure occurs, fluid injection causes a poromechanical response of the rock in which the rock expands and, as a result, total stress

components increase accordingly. The lower portion of the block of rock comprised between faults A and B, i.e., the half space with negative values of the y-coordinate, which contains the central point of these two faults whose results are plotted in Fig. 4, is compressed against fault A. Such deformation promotes a right-lateral movement of fault A (negative values of the shear stress) and a left-lateral movement of fault B (positive values of the shear stress). Since faults A and B undergo a right-lateral slip because of their orientation with respect to the principal stresses, fault A is destabilized by the poromechanical effect, while the upper part of fault B is destabilized and its lower part stabilized (in which lower and upper refer to the reference axis Y). This poromechanically-induced destabilization of fault A explains why failure is reached before its critical pressure is reached. The activation of shear slip, and the consequent stress redistribution, amplify failure on fault A (Fig. 2) and mitigate failure on fault B (compare the dashed and solid green lines in Fig. 4). These opposite effects, despite the similar orientation of the two faults, are due to the different location of the faults with respect to the antisymmetric variation of shear stress caused by poroelastic expansion (Fig. 5). Although critical pressure is not reached on faults A and B at the time of their activation, their rupture is mostly initiated by direct pressure effects and partially by the poroelastic effects. Once failure is initiated, shear-slip stress transfer plays a major role on fault stability, as shown by the lower reduction in the effective normal stress in fault B during the first two days of injection that stabilizes it. Furthermore, slip causes the fault to open up because
<figure><img src="image_6.png" /><figcaption>Fig. 4 Temporal evolution of effective normal and shear stresses changes (blue and yellow lines, respectively), Coulomb Failure Stress changes (green lines) and pore pressure (red line) under elastic (solid lines) and visco-plastic (dashed lines) behavior at the central location of the selected faults.</figcaption></figure>

The critical pressure for failure is represented by the horizontal dark red line. The gray shadowed regions correspond to times with increasing $ε_p$ for each fault according to Fig. 2. Note that the scale of stress variation is different in each plot. The behavior of the rest of the faults is shown in Supplementary Fig. S5.

<figure><img src="image_7.png" /><figcaption>Fig. 5 Shear stress variation $\Delta\tau_{xy}$ in the elastic scenario before and after the stop of injection. Observe the spatial distribution of the poroelastic stressing due to the injection-driven volume expansion during injection (a) and volume contraction during bleed-off (b and c), evolving as four antisymmetric lobes with respect to the injection well (black dot). Positive shear stress indicates left-lateral movement.</figcaption></figure>

of dilatancy, which causes undrained pore pressure drops that are subsequently recovered by diffusion. These dilatancy-induced pore pressure changes perturb the smooth evolution of pore pressure, stresses and fault stability observed when only elasticity is considered. Right after the stop of injection and subsequent bleed-off, poromechanical effects vanish. Therefore, fault A temporarily improves its stability and fault B fails for the second time coinciding with the immediate stop of poroelastic volume

expansion (see also Fig. 2). By bleeding off the well, i.e., opening the wellhead and achieving hydrostatic conditions within the well, poromechanical effects disappear almost completely, which may be the cause of the sudden second reactivation of fault B.

Fault C reactivates twice at day 4 and day 5. Both reactivations of fault C are due to pore pressure diffusion, poroelasticity and stress transfer from faults A and E, as shown by the difference between elastic and plastic Coulomb Failure Stress changes
<figure><img src="image_8.png" /><figcaption>**Fig. 6** Shear stress variation Δτ<sub>xy</sub> in the visco-plastic scenario at different times. **a** At day 5.8, fault B is affected by the injection. **b** At day 6.4, the rupture of fault B induces stress transfer, and amplifies the stress redistribution caused by poroelastic contraction due to bleed-off. **c** At day 6.7, fault F yields due to the combination of the shear-slip stress transfer and poroelastic mechanisms. Observe the evolution of the positive small lobe from poroelastic stressing after the bleed-off at day 6 that affects fault B. Note that the figure shows a close-up of the domain focusing on the area of interaction between faults B and F.</figcaption></figure>

(compare the solid and dashed green lines in the Fig. 4). Fault E fails for the first time during the injection and then twice during post-injection. Pore pressure variations are not sufficient to activate the fault, which is located outside of the pressurized region (Fig. 3). Additionally, poroelastic expansion has a stabilizing effect during injection because the induced right-lateral shear stress (Fig. 4) opposes to the left-lateral slip originated from the regional stress. On the other hand, elastic Δσ'ₙ slightly reduces and it tends to return to initial values after bleed-off, without completely recovering because residual overpressure continues propagating. Resulting from these variations, ΔCFS increases even after bleed-off and slightly starts decreasing at day 9 in the elastic case. Stresses in the visco-plastic scenario follow the same trends with a more complex behavior because each failure on the fault itself causes a drop in shear stress, and a drop in pore pressure that increases the effective normal stress, both contributing to improve stability. To better understand the role of shear-slip stress transfer, we run additional simulations in which alternatively each fault is the only one following a plastic-behavior, while the others follow an elastic behavior. We observe that without slip activation on fault C, fault E does not reach failure conditions. Overall, fault E is triggered by the combination of poroelastic stressing and shear-slip stress transfer during injection, and by the combination of stress transfer from the distant fault C and the poroelastic stressing relaxation due to the stop of injection and bleed-off in the post-injection stage.

Fault F is of special interest because it fails solely after the stop
of injection. Pore pressure reaches the critical pressure at day 4.
But similar to fault E, during the injection, poroelastically-
induced right-lateral movement opposes to the left-lateral slip
caused by the regional stress state and prevents the failure due to
the pore pressure increase. The first rupture of fault F, which
occurs 15 hours after the stop of injection, is caused by the
combination of poroelastic stress relaxation and shear-slip stress
transfer caused by the reactivation of fault B with the pore
pressure diffusion, as shown by the comparison with the
simulation in which only fault F has a visco-plastic behavior
(fault F never reaches failure conditions if fault B is not activated,
see Supplementary Fig. S3). Figure 6 shows the spatial evolution

of the shear stress τ<sub>xy</sub> around faults B and F. The reactivation of
fault B creates a local shear stress variation at the tips of fault F
(Fig. 6b) that destabilizes it (Fig. 6c). The rupture of fault F
stabilizes the fault due to the shear stress drop and an increase in
the effective normal stress caused by the dilatancy-induced pore
pressure drop. The second failure of fault F is caused by stress
transfer from the reactivation of a portion of fault B.

Discussion

We have revisited the intriguing case of co- and post-injection
induced seismicity at Basel EGS and we have identified the trig-
gering mechanisms by making use of coupled fluid flow and
geomechanics numerical simulations in a model that explicitly
includes a set of pre-existing faults based on in-situ observed
seismicity. We simulate fault reactivation, which is a necessary
condition for induced seismic events to occur, by considering the
non-elastic response of the faults. Despite a few simplifying
assumptions, simulation results are remarkably coherent with the
observed induced seismicity, both temporally and spatially. Our
results illustrate that accounting for poroelastic stressing and
non-elastic behavior, i.e., shear-slip stress transfer, is crucial to
reproduce the reactivation of certain faults.

In general, pore pressure diffusion is accepted as the main
triggering mechanism in EGS. Yet, this vision may be over-
simplifying, leading to inaccurate forecasting of induced seismi-
city, especially of post-injection seismicity. We estimate values of
critical pressures which are similar to the ones proposed by other
studies30,32,33. However, the simulation of the hydraulic stimu-
lation shows that pore pressure does not reach the critical pres-
sure for most of the faults. Thus, even though pore pressure
buildup has a direct effect on faults, especially in those in the
vicinity of the well, other triggering mechanisms are relevant and
should be taken into account to enable reliable estimates of
induced seismicity42.

Poromechanical volume expansion exhibits a wider and faster
front than pore pressure diffusion during co-injection, in
accordance with previous works by Duboeuf et al.⁴³ and
Krietsch et al.⁴⁴. This poroelastic effect, driven by fluid injection
that acts as a compressive loading, can improve or worsen the
<figure><img src="image_9.png" /><figcaption>**Fig. 7** Temporal evolution of the portion of each fault reactivated by the different combinations of triggering mechanisms. Faults G is not reactivated by any mechanisms.</figcaption></figure>

stability of faults depending on their orientation<sup>17,18</sup>. After the stop of injection, and more pronounced if bleed-off is applied, pressure gradients dissipate fast and poroelastic stress vanishes. Therefore, faults on which stability is enhanced by poroelastic stressing during injection are destabilized by volume contraction caused by bleed-off, e.g., faults B, E and F. The abrupt decrease of pore pressure and poroelastic stress are responsible for the immediate post-injection induced seismicity in certain zones of the reservoir, e.g., fault F. Lastly, shear-slip stress transfer affects the stability of nearby-faults and amplifies the induced seismicity (Fig. 6). Catalli et al.<sup>34</sup> suggested that the seismicity rate in the post-injection period is higher when stress transfer is taken into account than with pressure-induced seismicity only, enabling a better fitting of the observed seismicity. Our study confirms this finding, as we observe that only two faults in the vicinity of the well (faults A and B) fail when stress transfer is neglected while shear-slip stress transfer is required to reach shear failure conditions on the rest of the faults. The combination of pore pressure diffusion, poroelastic stress and shear-slip stress transfer<sup>18,27</sup> explains the fault rupture patterns underpinning the co-injection and post-injection induced seismicity at Basel. A quantification of the effects of the

reactivation capacity of each mechanism, or combination of mechanisms, is presented in Fig. 7, where we estimate the portion of each fault reactivated by direct pore pressure effects, poroelastic effects, i.e., pore pressure changes and induced poromechanical stresses, or the combination of poroelastic effects and stress transfer during injection and post-injection stages. To do so, we analyze the behavior of each fault mesh element. If pore pressure, estimated by the hydromechanical elastic model, reaches the critical pressure, we consider that fault portion as reactivated by direct pore pressure effects. If the Coulomb Failure Stress, estimated by the hydromechanical elastic model, is positive, we consider that fault portion as reactivated by the combination of pore pressure and poroelasticity. Finally, if a positive increment of the plastic strain is estimated by the hydromechanical visco-plastic model (and if the plastic strain has a value higher than $1e-3$), we consider that fault portion as reactivated by the combination of all the mechanisms. Pore pressure has a major impact on the stability of the faults located near the injection well before shut-in, and a minor impact after shut-in. Poroelasticity affects a wider region, but depending on the orientation and location of the faults, it can promote (faults B and C) or hinder (faults A and F) fault
reactivation. Note that fault F is reactivated if pore pressure is considered as the only triggering mechanism. However, poroelasticity hinders failure because it stabilizes the fault. Static stress transfer is responsible for the reactivation of faults that are not destabilized by the effects of pore pressure and poromechanical stresses (faults C and E). Fault D does not have a sufficient plastic strain increase and is not considered as reactivated, even if pore pressure alone, and its combination with poroelasticity could trigger fault reactivation. Post-injection fault reactivations are due to the fact that poroelasticity vanishes abruptly, and stabilization effects are quickly reversed; while pore pressure still diffuses after shut-in (fault F). These mechanisms are combined to the continuous stress redistribution from fault reactivation during and after the stimulation, which trigger the post-injection induced seismicity.

Despite our modeling approach permitting to identify the triggering mechanisms of the induced seismicity at Basel, simulation results cannot fully explain the reactivation of all faults. In particular, faults D and G, even though they reach failure conditions in our model, present a plastic strain that is too small to explain the observed cumulative seismic moment at these faults (Fig. 2). For the rest of the faults, the spatio-temporal evolution of faults reactivation correlates well with the observed seismic events. Yet, the largest event, which occurred at fault C shortly after the stop of injection, is not captured by our numerical model that reproduce two reactivations at times earlier than the reported ones in the injection stage (Fig. 2). This is probably due to our simplification of the fracture network into a few faults. It is possible that aseismic slip occurred on smaller fracture connecting the seismogenic faults, but no data is available to confirm this hypothesis. Our model identifies that multiple reactivations of fault C only occur with a slip-weakening friction, but the weakening, and reactivation of the real fault may have obeyed to a behavior different than just the linear strain-weakening assumed in our model.

For faults A, B, E and F, our 2D hydromechanical model of the Basel EGS is able to qualitatively reproduce the timing of fault reactivation. Although the adopted constitutive visco-plastic model does not allow to quantify the seismic magnitude, we compare the observed cumulative seismic moment with the numerically estimated plastic strain, which is proportional to slip and, therefore, to the moment magnitude. Interestingly, the temporal evolution of the numerically estimated plastic strains qualitatively corresponds to the one of the observed seismicity for most of the faults. Therefore, our results are relevant to identify the triggering effects of the multiple processes represented in the model.

Another simplification of our modeling resides in the adoption of a 2D domain. Although the monitored seismic events – exhibiting focal mechanisms with strike-slip movement and vertical dip – combined with the long open-hole section of the injection well – ensuring that overpressure field and poroelastic deformation are constant along the depth – suggest that a 2D horizontal section that crosses the vertical faults is reasonable to represent fault reactivation, static stress transfer is limited by this assumption. Indeed, stress redistribution from fault slip and earthquake-interactions are more complex in a 3D domain than in a 2D model. Improvements of our model could be achieved by modeling the 3D geometry of the fractured network, but such representation is extremely computationally challenging.

This novel analysis of the induced seismicity at Basel provides a substantial step forward in the general understanding of the physical processes that induce seismicity in the context of EGS hydraulic stimulation. Fault failures occurring during and after injection are located close to the well and farther away in the reservoir, respectively. Pore pressure diffusion and poromechanical

stress combined are the main triggering mechanisms during injection. Poromechanical effects extend farther and faster than pressure diffusion during injection (Fig. 8b). Furthermore, poroelastic stresses, depending on fault orientation with respect to the injection well, stabilize or destabilize faults during injection, and cause the opposite effect after the stop of injection as they rapidly diminish with pore pressure drop. After the stop of injection, pore pressure continues to advance further, leading to pore pressure increase far away from the well, which may induce some seismic events. Shear-slip stress transfer becomes dominant after the stop of injection, especially in faults far away from the injection well (Fig. 8c). In brittle rock, like the crystalline rock that is the target of most EGS projects, faults may present a slip-weakening friction, which may enhance the magnitude and frequency of induced seismic events$^{45}$. An analysis of the real-time monitored seismicity could allow calibrating the actual friction angle of the reactivated faults to improve the accuracy of the predictions of the reservoir stability. This improved understanding of the causal mechanisms of induced seismicity in EGS will contribute to have a better forecasting capability of induced seismicity and to come up with stimulation protocols that mitigate induced earthquakes, which are key points for the widespread development and management of geothermal projects.

## Methods

**Fault network.** To construct the geometry of the fault system, we use the Basel seismic database provided by Deichmann et al.$^{37}$, in which different seismic events are grouped on the basis of their focal mechanism, location and timing of occurrence. We assume that each cluster corresponds to a fault plane obtaining a network of 11 faults with different orientations and locations. To simplify the geometry, we combine a few clusters on the basis of similarities of focal mechanism, location and timing. This operation is applied to two sets of three clusters each, which allows us to reduce the original network of 11 faults to a simplified network of 7 faults (Supplementary Fig. S4 and Supplementary Table S1). To the representativeness of this simplification, we compare the failure-induced static stress redistribution in the simplified fault system and in the original fault system, i.e., in which each cluster is represented by a fault. We use Coulomb$^{46,47}$ to calculate the static stresses induced by fault slip according to linear elastic behavior of rock. A net slip *d* is imposed on each fault plane, calculated as$^{48}$

$$d = \frac{M_o}{G \times A}, \qquad (3)$$

where $G = E/(2(1 + \nu))$ is the shear modulus, which we assume to equal 21 GPa, E is Young's modulus, $\nu$ is the Poisson ratio, and A is the area of the slipping surface of the fault plane, which we estimate from the seismic cloud of the cluster. $M_o$ is the moment magnitude of the largest event of the cluster, derived from the magnitude of the seismic event $M_w$, and calculated by$^{49}$

$$M_o = 10^{1.5M_w+6.07}. \qquad (4)$$

Results of the slip-induced stress variation in the simplified and original systems are comparable, which confirms the effectiveness of the simplified fault network to model the case of Basel EGS (see Supplementary Fig. S4). Next, we simplify the three-dimensional model to a two-dimensional one by considering the projection of the cluster surfaces on the two-dimensional plane. It should be noted that Coulomb3 considers three dimensions, thus comparisons with the shear-slip stress transfer estimated by the two-dimensional HM numerical model have to be made carefully. Similarly, we can only qualitatively compare the net slip estimated through Eq. (3) (Table 1) with the shear slip displacements calculated by the HM numerical simulations to verify the consistency of the HM model.

**Hydro-mechanical problem.** We calculate the stress and pore pressure variations consequent to fluid injection by means of CODE_BRIGHT, a Finite Element Method (FEM) simulator that solves the fully coupled hydro-mechanical problem$^{50}$. Faults are modeled by finite-thickness elements, yet governed by a Mohr-Coulomb failure criterion together with a friction law. To reduce computational effort, we assume this continuum approach instead of an interface model with discontinuous displacement. The two methods may lead to equivalent results if a correct parametrization is performed$^{42,51}$. The mechanical governing equation to be solved is the momentum balance for the medium, expressed by

$$\nabla \cdot \sigma + \mathbf{b} = 0, \qquad (5)$$

where $\sigma$ is the stress tensor and $\mathbf{b}$ is the vector of body forces. The hydraulic governing equation is the mass balance of water, which for a fully saturated porous
<figure><img src="image_10.png" /><figcaption>**Fig. 8** Schematic illustration of the triggering mechanisms of induced seismicity. **a–c** Evolution of pore pressure, poroelastic stress and shear-slip stress transfer, respectively, before, during and after injection in a simplified generic fault network. The black dot in the center represents the injection well. **d–f** Mohr circles illustrating the stress state of the faults at each stage. **a** At initial conditions, faults F1 and F2 undergo right-lateral slip while fault F3 undergoes left-lateral slip (dashed arrows), in accordance with the maximum ($\sigma_H$) and minimum ($\sigma_H$) horizontal stresses. **b** During injection, pore pressure diffuses in the vicinity of the well. Poroelastic stressing extends farther and faster, and it exerts an inversed stress than the initial shear stress on F2 and F3, which are thus stabilized during injection. Combined with pore pressure, poroelastic stressing triggers the reactivation of F1, with the decrease of normal stress and the increase of the shear stress (**e**). Subsequently, F1 is stabilized by the shear stress drop (**c** and **f**). After the stop of injection, the pore pressure front continues to diffuse, while poroelastic stress relaxes (**c**). This change of direction leads to the increase of shear stress at the previously stabilized F2 and F3, which reach the failure envelope (**f**). The shear-slip stress transfer due to the reactivation of F2 affects F3, emphasizing the poroelastic effects until reaching failure.</figcaption></figure>

medium is expressed by

$$ \phi \frac{\partial P}{\partial t} + \frac{\partial}{\partial t} \varepsilon_v + \nabla \cdot \mathbf{q} = f^w, \qquad (6) $$

where $\phi$ is the rock porosity, $\beta$ is the fluid compressibility, $\varepsilon_v$ represents the volumetric strain, and $f^w$ is an external supply of water. $\mathbf{q}$ is the water mass flux and is expressed by Darcy's law

$$ \mathbf{q} = -\frac{k}{\gamma}(\nabla P - \rho \mathbf{g}), \qquad (7) $$

where $\gamma$ and $\rho$ are the fluid viscosity and density, respectively, $\mathbf{g}$ is the gravity vector and $k$ is the rock intrinsic permeability, considered as isotropic. The intrinsic permeability of the intact rock is considered to be a function of porosity by means of Kozeny's model

$$ k = k_o \frac{\phi^3}{(1-\phi)^2} \frac{(1-\phi_o)^2}{\phi_o^3}, \qquad (8) $$

with $k_o$ and $\phi_o$ being the reference values for intrinsic permeability and porosity of the rock matrix, respectively. Fault permeability variations are instead calculated as a function of the undergone deformations, adopting the “embedded model” proposed by Olivella and Alonso⁵², with fractures defined by their aperture embedded in a continuous finite element composed of rock matrix. The model assumes that variations of permeability are proportional to the square of the fault aperture variation, in agreement with the cubic law.

The coupling between Eqs. (5) and (6) is built through the elastic constitutive law, which relates stress tensor, $\boldsymbol{\sigma}$, strain tensor, $\boldsymbol{\varepsilon}$, and pressure as

$$ \Delta\sigma = K\varepsilon_v I + 2G\left(\varepsilon - \frac{\varepsilon_v}{3}I\right) - \Delta PI, \qquad (9) $$

where $K = E/[3(1 - 2\nu)]$ is the rock bulk modulus and $I$ is the first invariant of the stress tensor. While the intact rock is assumed to follow the standard linear elasticity model, the fault elements may present a non-elastic behavior according to a visco-plastic constitutive law, following a Mohr–Coulomb failure criterion, as

expressed by⁵³

$$ \frac{d\varepsilon_p}{dt} = \Gamma\langle\Phi(F)\rangle \frac{\partial\xi}{\partial\sigma}, \qquad (10) $$

where $\varepsilon_p$ is the visco-plastic strain, $\Gamma$ is the fluidity set at $1.00E-4 \text{ s}^{-1}\text{MPa}^{-m}$, $F$ is the yield function, $\xi$ is the flow rule and $\Phi(F)$ is the overstress function, which are defined as

$$ F = \sigma_m \cdot \sin\varphi(\eta) + \left[ \cos\theta - \frac{1}{\sqrt{3}} \sin\theta \cdot \sin\varphi(\eta) \right] \cdot \sqrt{J_2} - c(\eta) \cdot \cos\varphi(\eta), \qquad (11) $$

$$ \xi = \alpha \cdot \sigma_m \cdot \sin\psi + \left( \cos\theta - \frac{1}{\sqrt{3}} \sin\theta \cdot \sin\psi \right) \cdot \sqrt{J_2} - c(\eta) \cdot \cos\psi(\eta), \qquad (12) $$

$$ \Phi(F) = \begin{cases} 0, & \text{if } F \le 0 \\ F^m, & \text{if } F > 0 \end{cases}, \qquad (13) $$

where $c$ is cohesion equal to 2 MPa, $\eta$ is the weakening parameter of 0.01, $\alpha$ is a parameter for the plastic potential set at 1, $J_2$ is the second invariant of the deviatoric stress tensor and $\psi$ is the dilatancy angle, set at 3° for the faults. The invariants of the equations are $\sigma_m$, the effective mean stress, and $\theta$, Lode's angle. $m$ is a constant power to define the overstress function, set as 3 in the model. Fault elements are deformed elastically until stresses reach the shear yield surface ($F = 0$). When the yield surface is exceeded the fault begins to slip irreversibly, but stresses are allowed to remain beyond the yield surface for a range determined by the overstress function. To reproduce the degradation of frictional resistance consequent to the shear slip, we apply a friction weakening on fault C. The friction coefficient decreases linearly from the initial value, $\varphi^{\text{init}}$, to the residual one, $\varphi^{\text{res}}$, over a critical shear strain $\eta^*$
$$ \varphi(\eta) = \begin{cases} \varphi^{\text{init}}, & \eta \le 0 \\ \varphi^{\text{init}} + \frac{\varphi^{\text{rm}} - \varphi^{\text{out}}}{\eta^*}, & \eta \le \eta^* \\ \varphi^{\text{res}}, & \eta^* \le \eta \end{cases}. \quad (14) $$

The described model is able to reproduce the HM response to fluid injection and the activation of shear slip along the pre-existing fractures in the network.

## Data availability

The associated data is available at the repository DIGITAL.CSIC (https://digital.csic.es/handle/10261/275833). All numerical simulations, the presented results and findings of this study can be reproduced by the provided data.

Received: 30 September 2022; Accepted: 17 March 2023; Published online: 05 April 2023

## References

1. IEA. World energy outlook 2017 (IEA, 2017).

2. Majer, E. L. et al. Induced seismicity associated with Enhanced Geothermal Systems. Geothermics 36, 185-222 (2007).

3. Häring, M. O., Schanz, U., Ladner, F. & Dyer, B. C. Characterisation of the Basel 1 enhanced geothermal system. Geothermics 37, 469-495 (2008).

4. Ellsworth, W. L., Giardini, D., Townend, J., Ge, S. & Shimamoto, T. Triggering of the Pohang, Korea, Earthquake (Mw 5.5) by Enhanced Geothermal System Stimulation. Seismol. Res. Lett. https://doi.org/10.1785/0220190102 (2019).

5. Evans, K. F., Zappone, A., Kraft, T., Deichmann, N. & Moia, F. A survey of the induced seismic responses to fluid injection in geothermal and CO2 reservoirs in Europe. Geothermics 41, 30-54 (2012).

6. Grigoli, F. et al. Current challenges in monitoring, discrimination, and management of induced seismicity related to underground industrial activities: A European perspective: Challenges in induced seismicity. Rev. Geophys. 55, 310-340 (2017).

7. Ellsworth, W. L. Injection-induced earthquakes. Science 341, 1225942 (2013).

8. Ge, S. & Saar, M. O. Review: induced seismicity during geoenergy development—a hydromechanical perspective. J. Geophys. Res. Solid Earth 127, e2021JB023141 (2022).

9. Keranen, K. M. & Weingarten, M. Induced seismicity. Annu. Rev. Earth Planet. Sci. 46, 149-174 (2018).

10. Vilarrasa, V., De Simone, S., Carrera, J. & Villaseñor, A. Unraveling the causes of the seismicity induced by underground gas storage at castor, Spain. Geophys. Res. Lett. 48, e2020GL092038 (2021).

11. Raleigh, C. B., Healy, J. H. & Bredehoeft, J. D. An experiment in earthquake control at Rangely, Colorado. Science 191, 1230-1237 (1976).

12. Parisio, F., Vilarrasa, V., Wang, W., Kolditz, O. & Nagel, T. The risks of long-term re-injection in supercritical geothermal systems. Nat. Commun. 10, 4391 (2019).

13. Goebel, T. H. W., Weingarten, M., Chen, X., Haffener, J. & Brodsky, E. E. The 2016 Mw5.1 Fairview, Oklahoma earthquakes: evidence for long-range poroelastic triggering at >40 km from fluid disposal wells. Earth Planet. Sci. Lett. 472, 50-61 (2017).

14. Goebel, T. H. W. & Brodsky, E. E. The spatial footprint of injection wells in a global compilation of induced earthquake sequences. Science 361, 899-904 (2018).

15. Andrés, S., Santillán, D., Mosquera, J. C. & Cueto-Felgueroso, L. Delayed weakening and reactivation of rate-and-state faults driven by pressure changes due to fluid injection. J. Geophys. Res. Solid Earth 124, 11917-11937 (2019).

16. De Simone, S., Vilarrasa, V., Carrera, J., Alcolea, A. & Meier, P. Thermal coupling may control mechanical stability of geothermal reservoirs during cold water injection. Phys. Chem. Earth Parts ABC 64, 117-126 (2013).

17. Segall, P. & Lu, S. Injection-induced seismicity: poroelastic and earthquake nucleation effects: injection induced seismicity. J. Geophys. Res. Solid Earth 120, 5082-5103 (2015).

18. De Simone, S., Carrera, J. & Vilarrasa, V. Superposition approach to understand triggering mechanisms of post-injection induced seismicity. Geothermics 70, 85-97 (2017).

19. Bachmann, C. E., Wiemer, S., Goertz-Allmann, B. P. & Woessner, J. Influence of pore-pressure on the event-size distribution of induced earthquakes: pore pressure and earthquake. Geophys. Res. Lett. 39, https://doi.org/10.1029/2012GL051480 (2012).

20. King, G. C. P., Stein, R. S. & Lin, J. Static stress changes and the triggering of earthquakes. Bull. Seismol. Soc. Am. 84, 935-953 (1994).

21. Bhattacharya, P. & Viesca, R. C. Fluid-induced aseismic fault pore migration. Science 364, 464-468 (2019).

22. Cappa, F., Scuderi, M. M., Collettini, C., Guglielmi, Y. & Avouac, J.-P. Stabilization of fault slip by fluid injection in the laboratory and in situ. Sci. Adv. 5, eaau4065 (2019).

23. De Barros, L., Wynants-Morel, N., Cappa, F. & Danré, P. Migration of fluid-induced seismicity reveals the seismogenic state of faults. J. Geophys. Res. Solid Earth 126, e2021JB022767 (2021).

24. Catalli, F., Meier, M.-A. & Wiemer, S. The role of Coulomb stress changes for injection-induced seismicity: the Basel enhanced geothermal system: the Basel enhanced geothermal system. Geophys. Res. Lett. 40, 72-77 (2013).

25. Kettlety, T., Verdon, J. P., Werner, M. J., Kendall, J. M. & Budge, J. Investigating the role of elastostatic stress transfer during hydraulic fracturing-induced fault activation. Geophys. J. Int. https://doi.org/10.1093/gji/ggz080 (2019).

26. Schoenball, M., Baujard, C., Kohl, T. & Dorbath, L. The role of triggering by static stress transfer during geothermal reservoir stimulation: stress transfer during stimulation. J. Geophys. Res. Solid Earth 117, https://doi.org/10.1029/2012JB009304 (2012).

27. Yeo, I. W., Brown, M. R. M., Ge, S. & Lee, K. K. Causal mechanism of injection-induced earthquakes through the Mw 5.5 Pohang earthquake case study. Nat. Commun. 11, 2614 (2020).

28. Wibberley, C. A. J. & Shimamoto, T. Earthquake slip weakening and asperities explained by thermal pressurization. Nature 436, 689-692 (2005).

29. Eyre, T. S. et al. The role of aseismic slip in hydraulic fracturing-induced seismicity. Sci. Adv. 5, eaav7172 (2019).

30. Mukuhira, Y., Dinske, C., Asanuma, H., Ito, T. & Häring, M. O. Pore pressure behavior at the shut-in phase and causality of large induced seismicity at Basel, Switzerland. J. Geophys. Res. Solid Earth 122, 411-435 (2017).

31. Shapiro, S. A. Fluid-induced seismicity (Cambridge University Press, 2015).

32. Terakawa, T., Miller, S. A. & Deichmann, N. High fluid pressure and triggered earthquakes in the enhanced geothermal system in Basel, Switzerland: High fluid pressure in Basel. J. Geophys. Res. Solid Earth 117, https://doi.org/10.1029/2011JB008980 (2012).

33. Terakawa, T. Evolution of pore fluid pressures in a stimulated geothermal reservoir inferred from earthquake focal mechanisms. Geophys. Res. Lett. 41, 7468-7476 (2014).

34. Catalli, F., Rinaldi, A. P., Gischig, V., Nespoli, M. & Wiemer, S. The importance of earthquake interactions for injection-induced seismicity: Retrospective modeling of the Basel Enhanced Geothermal System. Geophys. Res. Lett. 43, 4992-4999 (2016).

35. Goertz-Allmann, B. P., Goertz, A. & Wiemer, S. Stress drop variations of induced earthquakes at the Basel geothermal site. Geophys. Res. Lett. 38, 2011GL047498 (2011).

36. Andrés, S., Santillán, D., Mosquera, J. C. & Cueto-Felgueroso, L. Thermoporoelastic analysis of induced seismicity at the Basel enhanced geothermal system. Sustainability 11, 6904 (2019).

37. Deichmann, N., Kraft, T. & Evans, K. F. Identification of faults activated during the stimulation of the Basel geothermal project from cluster analysis and focal mechanisms of the larger magnitude events. Geothermics 52, 84-97 (2014).

38. Kraft, T. & Deichmann, N. High-precision relocation and focal mechanism of the injection-induced seismicity at the Basel EGS. Geothermics 52, 59-73 (2014).

39. Ohnaka, M. Frictional characteristics of typical rocks. J. Phys. Earth 23, 87-112 (1975).

40. Miller, S. A. Modeling enhanced geothermal systems and the essential nature of large-scale changes in permeability at the onset of slip. In Crustal Permeability (eds. Gleeson, T. & Ingebritsen, S. E.) 353-362 (John Wiley & Sons, Ltd, 2016). https://doi.org/10.1002/9781119166573.ch27

41. Jaeger, J. C., Cook, N. G. & Zimmerman, R. Fundamentals of rock mechanics (John Wiley & Sons, 2009).

42. Zareidarmiyan, A., Parisio, F., Makhnenko, R. Y., Salarirad, H. & Vilarrasa, V. How equivalent are equivalent porous media? Geophys. Res. Lett. 48, e2020GL089163 (2021).

43. Duboeuf, L. et al. Aseismic motions drive a sparse seismicity during fluid injections into a fractured zone in a carbonate reservoir: injection-induced (A) seismic motions. J. Geophys. Res. Solid Earth 122, 8285-8304 (2017).

44. Krietsch, H. et al. Hydromechanical processes and their influence on the stimulation effected volume: observations from a decameter-scale hydraulic stimulation project. Solid Earth 11, 1699-1729 (2020).

45. Kirkpatrick, J. D. & Shipton, Z. K. Geologic evidence for multiple slip weakening mechanisms during seismic slip in crystalline rock. J. Geophys. Res. 114, B12401 (2009).

46. Lin, J. & Stein, R. S. Stress triggering in thrust and subduction earthquakes and stress interaction between the southern San Andreas and nearby thrust and strike-slip faults: Stress triggering and fault interaction. J. Geophys. Res. Solid Earth 109, (2004).

47. Toda, S., Stein, R. S., Richards-Dinger, K. & Bozkurt, S. Forecasting the evolution of seismicity in southern California: animations built on earthquake stress transfer. J. Geophys. Res. 110, B05S16 (2005).
48. Steketee, J. A. Some geophysical applications of the elasticity theory of dislocations. *Can. J. Phys.* **36**, 1168–1198 (1958).

49. Kanamori, H. & Brodsky, E. E. The physics of earthquakes. *Rep. Prog. Phys.* **67**, 1429–1496 (2004).

50. Olivella, S., Gens, A., Carrera, J. & Alonso, E. E. Numerical formulation for a simulator (CODE_BRIGHT) for the coupled analysis of saline media. *Eng. Comput* **13**, 87–112 (1996).

51. Cappa, F. & Rutqvist, J. Modeling of coupled deformation and permeability evolution during fault reactivation induced by deep underground injection of CO2. *Int. J. Greenh. Gas Control* **5**, 336–346 (2011).

52. Olivella, S. & Alonso, E. E. Gas flow through clay barriers. *Géotechnique* **58**, 157–176 (2008).

53. Vilarrasa, V., Bolster, D., Olivella, S. & Carrera, J. Coupled hydromechanical modeling of CO2 sequestration in deep saline aquifers. *Int. J. Greenh. Gas Control* **4**, 910–919 (2010).

## Acknowledgements

A.B. and V.V. acknowledge funding from the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation Program through the Starting Grant GEoREST (www.georest.eu) under Grant agreement No. 801809. IDAEA-CSIC is a Centre of Excellence Severo Ochoa (Spanish Ministry of Science and Innovation, Grant CEX2018-000794-S funded by MCIN/AEI/ https://doi.org/10.13039/501100011033). This research has been carried out within the framework of the activities of the Spanish Government through the "Maria de Maeztu Centre of Excellence" accreditation to IMEDEA (CSIC-UIB) (CEX2021-001198).

## Author contributions

A.B., S.D. and V.V. designed the study. A.B. conducted the numerical simulations and prepared the figures. A.B., S.D., V.V. and S.G. contributed to the interpretation of results, writing, and editing the paper.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s43247-023-00764-y.

**Correspondence** and requests for materials should be addressed to Auregan Boyet.

**Peer review information** *Communications Earth & Environment* thanks Vanille Ritz, Stephen A. Miller and Antoine B. Jacquey for their contribution to the peer review of this work. Primary Handling Editors: Luca Dal Zilio and Joe Aslin.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_11.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023