Thermo-Poroelastic Analysis of Induced Seismicity at the Basel Enhanced Geothermal System

Geothermal energy has emerged as an alternative to ensure a green energy supply while tackling climate change. Geothermal systems extract the heat stored in the Earth’s crust by warming up water, but the low rock permeability at exploitation depths may require the hydraulic stimulation of the rock fracture network. Enhanced Geothermal Systems (EGS) employ techniques such as hydro-shearing and hydro-fracturing for that purpose, but their use promotes anthropogenic earthquakes induced by the injection or extraction of fluids. This work addresses this problem through developing a computational 3D model to explore fault reactivation and evaluating the potential for earthquake triggering at preexisting geological faults. These are included in the model as frictional contacts that allow the relative displacement between both of its sides, governed by rate-and-state friction laws and fully coupled with thermo-hydro-mechanical equations. We apply our methodology to the Basel project, employing the on-site parameters and conditions. Our results demonstrate that earthquakes which occurred in December 2006 in Basel (Switzerland) are compatible with the geomechanical and frictional consequences of the hydraulic stimulation of the rock mass. The application of our model also shows that it can be useful for predicting fault reactivation and engineering injection protocols for managing the safe and sustainable operation of EGS.


Introduction
The 2030 Agenda for Sustainable Development, adopted by world leaders in September 2015, was established by the United Nations and comprises 17 goals and 169 targets to be fulfilled by the ratifying countries. The Agenda aims to overcome all forms of poverty, while tackling climate change and environmental protection. Goal 7-"Affordable and Clean Energy"-aims to ensure access to affordable, reliable, sustainable, and modern energy for all [1]. Goal 7 goes hand-in-hand with Goal 13, "Climate Action", which aims to take urgent action to combat climate change and its impacts. Geothermal energy emerges as an alternative renewable energy to reach both goals [2,3], as it is affordable and clean. Geothermal systems extract the heat stored in the Earth's crust by warming up water or a mixture of water and gas. The fluid is circulated down through injection wells, heated by the contact with rocks, and returned to the surface through production wells to form a closed loop [4,5]. Hot water or steam is then transformed into a marketable product, such as electricity. Nevertheless, in most geothermal reservoirs, rock permeability at exploitation depths is very low, rendering geothermal projects economically infeasible. The challenge of permeability enhancement has been addressed by the so-called Enhanced Geothermal Systems (EGS) [2,6,7]. tendency to slip has proven very useful to understand some of the geomechanical challenges posed by subsurface energy technologies (e.g., [46,[57][58][59][60][61][62][63][64]). We aim to analyze the hydro-shearing effects during the stimulation at the Basel-1 well using our three-dimensional (3D) model. We use, as input, the same stress-field and rock properties estimated for the main hydrogeological units at Basel-1, and simulate the stimulation using a transient well-head pressure and flow evolution similar to the values reported at the study site. The proposed model is a useful tool for engineers and practitioners to answer the essential question for the development of the geothermal energy as an affordable, clean, and safe renewable energy source: what is the optimal exploitation protocol that minimizes seismic risk and maximizes economic and energy performance?

Materials and Methods
In this section, we describe our thermo-hydro-geomechanical model, which fully couples heat transport, rock deformation, and fluid flow; faults are described as interfaces with friction governed by the rate-and-state law. We perform numerical simulations of the hydraulic stimulation operations conducted at the DHM project in December 2006 to show the ability of our model to characterize fault reactivation. We adopt as inputs the recorded well-head pressure and flow evolutions, as well as on-site material parameters.

Frictional Strength and Resistance of Faults
We employ the Amontons-Coulomb theory as a constitutive model of rock friction. The frictional strength, τ * , that impedes sliding of a static, cohesionless contact interface is given by τ * =µ|σ| [65][66][67]. In general, the relation between these Coulomb magnitudes and the shear stress acting on the contact plane, τ, depends on the sliding regime. A static interface satisfies τ≤τ * and for sliding under quasi-static conditions, the relationship τ≈τ * holds.
In the fluid-saturated media, frictional strength is defined using the effective normal stress, σ =σ+p, where σ is the total normal stress acting on the contact, and p is the pore pressure of the fluid. In the above and following expressions, tensile stresses are positive, and pore pressure is positive when above the atmospheric value. We assume that effective normal stresses remain compressive on contact surfaces.

Rate-and-State Models for Interfaces
Faults are often assumed to be well-oriented for failure but locked prior to injection, in such a way that the onset of slip on an inactive fault-reactivation-is essential to understand the geomechanics of induced seismicity.
Rate-and-state formulations gather the traditional concepts of static and dynamic friction by including the dependence of µ on the slip velocity and history of sliding [66,[68][69][70]. These models were derived from laboratory experiments of unidirectional slip in the double-direct shear configuration, and account the response of µ after step changes in slip velocity or normal stress [71]. For a frictional interface that is sliding at velocity V, the definition of µ reads: where µ * is the steady-state coefficient at the reference slip velocity V * , a is the direct-effect parameter, and b is the friction evolution parameter. θ is the state variable, and θ * =D c /V * is its steady-state value, where D c is the characteristic slip memory distance over which τ * evolves once the system is perturbed [72]. Several definitions for θ have been proposed according to the rate-and-state friction models. Deep physical understandings and theoretical analyses [70,73,74], as well as comparisons with experiments of velocity steps [75,76], shear stress steps [77], and normal stress steps [78][79][80][81] have allowed us to elucidate the relative advantages and disadvantages of the proposed formulations. Our study focuses on those models that incorporate a dependence on the effective normal stress rate. In such a sense, the Linker-Dieterich model [49] generalizes the basic rate-and-state laws by including a term in the state evolution equation, applicable to both the "slip" [69] and "aging" [68] laws [73]: In the above equations, the empirical parameter α controls the stressing-rate effect on the state variable, ranging from 0 to µ 0 [82]. Since we focus on the reactivation of faults which are initially at rest, we adopt a regularization of the rate-and-state models in the limit of zero slip speed, V=0, proposed by Yang et al. [83]. µ is then defined as: with µ 0 being the initial friction coefficient. We implement in our model the equation for the aging law: which is equivalent to the regularized model used by Tal et al. [84] with a threshold velocity V th = V * .

Thermo-Hydro-Mechanical 3D Model of Fault Reactivation
The rock is modeled as a thermo-poroelastic saturated material with single-phase flow. We adopt the classical theory of linear poroelasticity [85,86] and solve for the combination of fluid pressure, rock deformation, temperature, and frictional contact on the fault [31,32]. The solid and mass conservation, as well as mechanical equilibrium are coupled using the effective stress. The quasi-static Biot equations for a porous medium read: where α B is the Biot coefficient, ε v =tr (ε) (with ε being the infinitesimal strain tensor) is the volumetric strain, k is the intrinsic permeability of the porous medium, η f is the fluid dynamic viscosity, ρ f is the fluid density, p is the pressure field, and σ is the total stress tensor. We consider a linear poroelastic material under small deformations, as well as plane strain conditions. Then, the effective stress tensor, σ =σ + α B pI, is a linear function of strains, σ =2Gε el + λtr(ε el )I, where λ and G are the Lamé constants, ε el =ε − ε th = 1 2 ∇u + ∇u T − ε th is the elastic strain tensor, the result of subtracting the thermal strains to the total strain tensor ε, with u being the displacement field. The storage coefficient, S=φχ f + (α B − φ)χ s , depends on the rock porosity, φ, and on the fluid and solid matrix compressibilities, χ f and χ s , where χ s =(1 − α B )/K, and K=λ + 2 3 G is the bulk modulus of the porous matrix.
Conservation of energy reduces to the heat equation [87]: where c s is the heat capacity of the rock, c f is the fluid heat capacity, κ s is the thermal conductivity of the rock, and κ f is the fluid thermal conductivity. (ρc) sat y κ sat are the saturated values of the product of density by heat capacity and thermal conductivity. Q includes the source or sink terms. Temperature changes propagate by diffusion and convection throughout the rock mass and the fluid. This changes produce thermal strains controlled by the thermal expansion coefficient α T : where T is the temperature field, T 0 is the reference initial temperature, and α T is the solid thermal expansion coefficient. We consider that the fluid keeps in a liquid state and its properties do not change with temperature or pressure, due to the high temperature and pressure conditions at the usual depths of the EGS projects. This couples the diffusive part of heat transport, while the convective part is modeled by introducing the Darcy velocity field v in Equation (7).

Case Study: The Deep Heat Mining Project in Basel, Switzerland
The DHM project in Basel was a milestone in geothermal energy. The know-how gained after the DHM project and the seismic events of 2006 boosted the development of EGS systems [52,53]. The first stage of the project was drilling the 5000 m depth Basel-1 well. The well reached a crystalline granitic rock basement at a temperature of 200 • C. The well was also employed for the hydraulic stimulation of the reservoir and the field characterization [56].
The number of fractures between 4629 m and 5000 m depth was between 0.2 and 0.3 fractures per meter [56]. The preferred fracture direction and orientation were NW-SE to NNW-SSE with dips greater than 60 • , although the measures of the hypocenter locations during seismic events detected new families of fractures [88]. The events with greatest magnitude which rolled around in 2006 took place on a family with a deviation of 10 • with respect of the maximum principal stress direction and a dip of 80 • . The orientation of the principal stresses was deduced from acoustic geophysical studies within the Basel-1 well. The minimum principal stress, σ h,min , had an orientation of 54 ± 14 • and the maximum one, σ h,max , 144 ± 14 • [88]. These orientations are consistent with the in situ stress state in the upper Rhine Basin deduced from previous seismic events in the crystalline rocky massif [56].
The magnitude of the principal stresses were also quantified. The tectonic ratio of the minimum principal stress, σ h,min , to the vertical one is 0.7, and the ratio of the maximum principal stress, σ h,max , to the vertical one is 1.6, in such a way that the vertical stress is the intermediate principal stress [56]. Mechanic boundary conditions are defined by the expressions: where d is the depth. In these expressions, compressive stresses are assumed positive. The imposed stresses at the boundaries are assumed constant in time, which is a feasible hypothesis even if simulation time exceeds one decade [89]. The temperature at the bottom of the well is between 190 and 200 • C, with a thermal gradient of 40 • C/km [55].

The Basel 3D Model
The geometry of the 3D model domain is a 1.5 km 3 cube that is supposed to be homogeneous porous rock mass. In the center of the cube, we situated a sphere divided into two hemispheres with their intersection plane representing the main fault. The injection well is modeled as a cylinder with a diameter of 1 m and a height of 380 m whose geometry is subtracted from the solid cube. The well is vertical and located at a distance 50 m away from the center of the fault (Figure 1). The diameter of the simulated well is higher than the real one due to mathematical issues, where the diffusion of the injected fluid should not be simulated with elements which are too fine. The geometry of the model is described in Figure 1.  We focus on the depth range between 4050 and 5550 m, and the injection takes place through the well between 4620 to 5000 m depth. The x axis is parallel to the maximum principal stress (σ x =σ h,max ), the y axis to the minimum principal stress (σ y =σ h,min ), and the z axis to the intermediate principal stress. We apply vertical stress to simulate overburden strata according to where the z coordinate ranges between 0 and 1500 m in our model.
The complexity of the fracture network requires the adoption of some simplifications to define the orientation of the fault plane, based on the hypocenter locations estimated from the motion data records of the historical seismic activity [88]. We adopt the most unfavorable fracture family as the preferred sliding plane, which has a direction deviated approximately 10 • with respect to the direction of the maximum principal stress, which dips 80 • towards the SW. Most of the detected hypocenters do not separate more than 50 m with respect to the upper injection section of the well [55]. For that reason, the fault plane of the model is oriented 10 • with respect to the x-axis and dipping 80 • to the SW, with the injection well-located at a distance 50 m away from its center. We refine the mesh of the model in the area of the central sphere, as shown in Figure 2, and force mesh conformity throughout the fault plane. We use tetrahedral elements-quadratic for the solid mechanics equations and linear for fluid flow and heat transport. Mechanic boundary conditions are equal to tectonic stresses and are imposed on three faces of the rectangular domain (Figure 2). We impose a reference hydrostatic pressure, p = ρ f gd = ρ f g(5550 − z), along all external boundaries, where d is the depth, and the temperature is set at the geothermal gradient T = 0.04 [C • /m]·d = 0.04 [C • /m]·(5550 − z).
We impede the displacements in the normal direction and impose no-flow (thermal and hydraulic) conditions on the other ones. At the fault plane, we impose the contact condition, allowing relative tangential displacements between its edges. We consider the fault is almost impermeable and has the same thermal properties as the rest of the domain, resulting in thermal continuity. We simplify the fluid flow around the injection well by imposing a volumetric flux along the boundary of an effective injection region (the lateral surface of the cylinder). During the simulations, temperature is also fixed in the injection well.
Based on the in-site measurements, we define feasible parameters for each physical process of the model.
The Young modulus, the Poisson ratio, and the density of the solid skeleton are usual values for crystalline rock formations of granitic type. The parameters related to the fluid, such as density or compressibility of water, are also usual values, unlike viscosity, whose value of 2.4 × 10 −4 Pa·s, is lower than the viscosity of water at ambient temperature (η f = 10 −3 Pa·s) due to the high temperatures at such depths [23]. Permeability and porosity have been chosen according to the characteristics of the rock mass. Lastly, thermal parameters have been taken from [90] that collects data of thermal conductivities and heat capacities from different materials. Coupling between the flow and mechanical problems have been included using the Biot coefficient α B = 1 and between the thermal and mechanical physics with the thermal expansion coefficient α T = 8 × 10 −6 .
The parameters of the "rate-and-state" model have been chosen within a feasible range to emulate the fault reactivation at a similar time scale. We consider the fault is arrested before reactivation, being the slip speed V = 0. Hence, Equations (1) and (2) can be simplified, so it is only necessary to define the parameters α, b, D c , and V * . We adopt the aging law and the rate-and-state parameters, b = 0.03, D c = 700 µm, V * = 10 −9 m/s, and µ 0 = 0.55. Table 1 lists the parameters of the model. The constitutive laws for fault strength are given by Equations (1) and (2), and the frictional contact on the fault is modeled using an Augmented Lagrangian formulation [91]. We solve, in a monolithically-coupled fashion, the field Equations (5)-(8) and the rate-and-state aging law (1) and (2) with the frictional contact variables [31,32].

Calibration
The hydraulic stimulation of the DHM geothermal project required the injection of 11, 600 m 3 of water, following the protocol depicted in Figure 1b [56]. This volume was injected prior to the 8 December 2006 earthquake sequence. Since we simulated the on-site conditions during the injection operations, we imposed an inflow velocity at the injection well of the model that is the result of dividing the flow rate by the lateral surface of the cylinder q iny =Q 0 /(2πrh), where r is the radius of the cylinder and h is the height where injection takes place. We assumed that the injection temperature remained constant and equal to ambient rock temperature. Figure 3a,b shows a plot of the simulated evolutions of pressure and injected flow rate computed with our model, and the registered data at the Basel-1 injection well in 2006. We adopted as the flow boundary condition the measured injected flow rate at the injection well (see Figure 3a). Our computed injected pressure evolution initially differs from the measured values on-site, Figure 3b. The difference between pressure observed in real data and model results arises from the weakening and fracturing of the rock in the vicinity of the well, allowing water to flow through the rock matrix and its fractures. After 5 days of injection, when the largest earthquakes occurred, the difference between our simulated results and the measured data were drastically reduced, showing that our model is properly adjusted.  The results of pore pressure and temperature fields are included in Figure 4. We define a reference horizontal plane at 4800 m depth (z = 750 m), where the changes of pore pressure and temperature cause by the injection of cold water are displayed. We show in Figure 4a the locations of the reference plane, the fault plane, and the injection well. We also include the buildup of pore pressure around the injection well and near the fault plane at the time that reactivation occurs (t = 5.5 days). The injection cools down a small area around the well which does not reach the fault (Figure 4b). Since the time for heat diffusion is higher than the one for pressure diffusion, pore pressure changes around the well are much faster than temperature changes. As fault temperature remains constant during the stimulation phase (Figure 4b), pore pressure increases on the fault plane (Figure 4c). We explain in the next section how pore pressure on the fault plane controls the fault reactivation and the frictional properties of the contact.  . In (b) we display the increment of temperature due to the injection at the reference plane, with an inset that zooms the surroundings of the injection well and shows that the temperature diffusion is much slower than pressure propagation. In (c) we show the results of the pore pressure increase at the fault plane. The vertical axis of the image corresponds to the maximum slope line of the fault plane and the horizontal axis corresponds to a horizontal direction in the 3D model deviated 10 • with respect to the x-axis.

Fault Reactivation
Fault reactivation is the onset of fault slip. It depends on the variables involved in fault stability, such as frictional strength µ|σ | or shear stress τ. We quantify fault stability through the Coulomb Failure Function, defined as CFF=µ|σ |-τ. The CFF equals to zero when the fault is at rest, and when CFF is less than zero the fault reactivates, given that the shear stress τ exceeds the frictional strength of the contact µ|σ |. Changes in Coulomb Failure Function can be used as a proxy for fault weakening (∆CFF<0) or strengthening (∆CFF>0). In that sense, we show in Figure 5 the increase in the failure function ∆CFF=CFF(t) − CFF(t = 0) on the fault plane, which indicates how the fault weakens due to the effect of fluid injection. We include in Figure 6 our computed results with the 3D model for the variables involved in the frictional stability on the fault plane. The distribution of effective normal stress, |σ |, acting on the fault plane (Figure 6a) at the onset of the slip shows that there is a decrease in fault effective compressions. Moreover, the spatial distribution is symmetric with respect to a vertical axis. The decrease in effective normal stress is almost the same as the increase in pore pressure (Figure 4c). Differences arise from poroelastic and thermal effects, and indicate that pore pressure changes dominate over thermal and poroelastic effects.
We plot the modulus of the shear stress τ on the fault at reactivation time in Figure 6b. Shear stress increases on the north side of the fault and decreases on the right side. This response is caused by the increase in pore pressure, as well as by the poroelastic effects accounted in our model [16][17][18][19]31,32]. The poroelastic effects are also coupled with the fault orientation and tectonic stresses, which contribute to the asymmetry of the results.
We show the value of the friction coefficient, µ, in Figure 6c, computed with the rate-and-state law (Equatons (3) and (4)). Since the fault is initially at rest (V = 0), the observed evolution of the friction coefficient is attributed to the decrease in the effective stress through the Linker-Dieterich term and the α-parameter in Equation (2) [79]. The Coulomb Failure Function CFF, Figure 6d, indicates that fault reactivation occurs after 5.5 days of injection. The asymmetry in the CFF distribution remarks the influence of shear stresses, in contrast with the symmetry of effective normal stress and friction coefficient that are directly pore-pressure dependent. This pattern of symmetry and asymmetry of stresses before and at reactivation also influence on the nucleation and rupture phases of the earthquake. We illustrate the evolution of stresses and frictional variables up to the fault reactivation at a control point on the fault. The point is located at the left-half part of the fault, and it is the first point where the ratio of the acting shear stress to the effective normal stress equals the frictional strength-that is, it is the first fault point at which slip occurs (blue area where CFF = 0 in Figure 6d). Figure 7 shows the evolution of the friction at the control point up to the reactivation. The mobilized friction (green line) is the ratio τ/|σ | that normalizes the shear stress acting on the fault with the effective normal stress, and the friction coefficient µ (blue line) is a dimensionless frictional strength. When the mobilized friction equals frictional strength at 5.5 days (τ/|σ | = µ), the fault reactivates.
We observe that due to the injection protocol that systematically increases the flow rate, and consequently, the pore pressure, the slopes of both the friction coefficient and the mobilized friction curves increase. Moreover, the friction coefficient changes from µ 0 = 0.55 to 0.61 at the reactivation instant due to the effect of the variation of the effective normal stress on the friction coefficient. It delays the fault reactivation, which changes from 4.6 days if effective normal stress rate is disregarded, and to 5.5 days if the rate is accounted for. Therefore, the effect of normal stress rate on friction coefficient needs to be taken into account, as it partially controls the reactivation time. These results verify that our models properly reproduce the on-site reactivation of the reservoir's representative fault after 5.5 days of injection, elapsed from December 3 to 8.
Our methodology may be a useful tool for engineers, practitioners, and stakeholders to assess fault reactivation under real conditions of natural stresses, temperatures, and injection protocols. The application of our model to the DHM project in Basel has shown its ability to predict fault reactivation and demonstrated that the earthquake sequence occurred in December 2006 may have been caused by the hydraulic stimulation of the rock mass. Our model can also be useful for assessing new injection protocols, as well as stimulating and managing the operation of EGS system in the short and long term.

Conclusions
Geothermal energy emerges as an alternative renewable energy to ensure access to affordable, reliable, sustainable, and modern energy for all the world. In most of the geothermal deposits, rock permeability at exploitation depths is very low, rendering geothermal projects economically infeasible. This drawback has been solved by the so-called Enhanced Geothermal Systems (EGS), where rock permeability is enhanced through the so-called rock stimulation. One of the most widely used techniques is hydro-shearing, which reactivates preexisting joints by initiating shear failure. Water is injected under high pressure, reducing normal stress across them and eventually triggering shear failure. A major environmental issue for these techniques is induced seismicity as a result of water injection.
Here, we have presented a finite element model for the simulation of fault reactivation in poroelastic media with rate-and-state friction law. Our model monolithically couples fluid flow, rock mechanics, heat transport, and rate-and-state friction. We have conducted three-dimensional simulations of fault reactivation with frictional strength governed by a Linker-Dieterich law, embedded in a poroelastic homogenous media, and driven by fluid injection. The Linker-Dieterich law accounts for the effect of effective stress rate on the friction evolution, and is able to explain the observed delays in fault reactivation. We applied our model to simulate the hydro-shearing effects during the stimulation at the Basel-1 well at the Deep Heat Mining geothermal project in Basilea (Switzerland). We adopted as input the stress field and rock properties estimated for the main hydrogeological units at Basel-1, and simulated the stimulation phase using a transient well-head pressure and flow evolution similar to the values reported at the study site. Our three-dimensional model satisfactorily reproduced the registered injection flow rate and pressure injection, as well as the time of the fault reactivation.
Our simulated results indicate that thermal effects are negligible during the stimulation phase. Temperature changes occur in a small area around the injection well and do not reach the fault plane. Nevertheless, pressure changes, together with poroelastic effects, weaken the fault, leading eventually to fault reactivation and the onset of slip. Our model was able to reproduce the instant at which fault reactivation occurred at the Basel-1 site, demonstrating that the earthquake sequence that occurred in December 2006 in Basel was caused by the hydraulic stimulation of the geothermal reservoir.
Our methodology emerges as a useful tool for engineers, practitioners, and stakeholders to assess fault reactivation under real conditions of natural stresses, temperatures, and injection protocols. The application of our model to the DHM project in Basel has shown its ability to predict fault reactivation in a real case. Our model can also be useful for assessing new injection protocols, as well as stimulating and managing the operation of an EGS system in the short and long term. Funding: This research was supported by the Universidad Politécnica de Madrid through the special programme for young scientists under grant VJIDOCUPM19DSS ("Programa Propio de I+D+I de la Universidad Politécnica de Madrid. Convocatoria de ayuda dirigida a jóvenes investigadores doctores para fortalecer sus planes de investigación"). DSS gratefully acknowledges funding from the Fundación BBVA though Becas Leonardo a Investigadores y Creadores Culturales 2019 (Proyecto realizado con la Beca Leonardo a Investigadores y Creadores Culturales 2019 de la Fundación BBVA).