Investigating seismicity rates with Coulomb failure stress models caused by pore pressure and thermal stress from operating a well doublet in a generic geothermal reservoir in the Netherlands

Abstract

The utilisation of geothermal energy in the Netherlands is primarily focused on deep sedimentary aquifers, which are often intersected by major faults. Geothermal operations (i.e. fluid production and injection) may alter the effective stress state along these faults and trigger induced seismic events. Pore pressure perturbations have been generally considered the main driver of injection-induced seismicity. However, thermal stresses caused by temperature gradients between the re-injected cold fluid and the reservoir rock may also contribute to the triggering of earthquakes in geothermal reservoirs. While existing geothermal power plants operating in sandstone reservoirs did not produce any major induced seismicity, it is a matter of debate whether a reduction in the temperature of the re-injected fluid could increase the seismic hazard potential. In this study, we applied modified Gutenberg–Richter statistics based on frictional Coulomb stress variations implemented in a coupled thermo-hydro-mechanical model to estimate the seismic hazard caused by the operation of a geothermal doublet. We conducted a systematic parametric study to assess and rank the impact of different intrinsic (geological) and extrinsic (operational) parameters on the induced seismic hazard potential. We identified a competing mechanism between induced variations in pore pressure and thermal stress within the reservoir in controlling induced seismicity. We found that stress changes induced by pore pressure variations are the main cause of seismic hazard, although thermally induced stresses also contribute significantly. The results indicate that by optimising the operational parameters it is possible to increase production efficiency while maintaining a long-term control over the fluid injection-induced seismicity.

Introduction

In the Netherlands, hydrocarbon exploitation is widespread, and the existing oil and gas wells combined with the knowledge of the subsurface is a suitable base for geothermal project development (Vinci et al., 2021). Geothermal energy is a local base load energy source for heating, cooling and electricity generation. At present, geothermal energy is extracted from deep sedimentary aquifers (Willems et al., 2017). These reservoirs are commonly intersected by fractures and faults. Permeable faults can improve production and re-injection as they act as preferential fluid pathways. On the other hand, faults also carry the risk of hosting seismic events induced by geothermal operations.

A series of fluid injection-induced earthquakes led to the cessation of geothermal projects in Basel, Switzerland (magnitudes 2.6–3.4) in 2006 (Deichmann & Giardini, 2009) and in Strasbourg, France (magnitudes 3–3.6) occurring between 2019 and 2021 (Schmittbuhl et al., 2021) as they were felt by the population, even though they resulted in only non-structural damage of residential buildings. In Pohang, South Korea, a magnitude 5.5 earthquake occurred in 2017, two months after a series of hydraulic stimulations related to an Enhanced Geothermal System (EGS) project, which was the most damaging earthquake since recording of seismicity started in Korea in 1905 (Kim et al., 2018; Yeo et al., 2020). We note that the events above were hosted within the crystalline basement, considered to be seismically more active (Buijze et al., 2019). Nevertheless, these observations highlight the importance of understanding the physical processes at play, when the in-situ stress conditions of a reservoir are modified in response to pressure and temperature variations caused by fluid injection (Zang et al., 2014).

Induced seismicity is an important point of discussion in the Netherlands where induced events have been recorded mainly related to gas production since as early as 1986 (Elk et al., 2017). The largest induced events in the country are related to the Groningen gas field with magnitudes of up to 3.6, which based on population concern ultimately led to the government’s decision to stop production before the depletion of the reservoir (Muntendam-Bos et al., 2022). It is worth noting that the seismic activity observed in the Groningen gas field is related to reservoir depletion rather than geothermal operations, while our study focuses on the latter. However, it is important to note that all geothermal wells in the Netherlands exploit the same sandstone formation (Rotliegend) as the gas field for production (Buijze et al., 2019). Furthermore, the sensitivity of the population to induced seismicity in connection to the Groningen gas field can have an impact on the perception of existing and future geothermal projects (Schultz et al., 2022). The placement of geothermal power plants in close proximity to consumers poses a major challenge, especially in the case of heat generation plants. Although it is not feasible to exclude the occurrence of minor induced seismic events that may be perceived by the population, measures are available to minimise their frequency and magnitude (Knoblauch & Trutnevyte, 2018).

To this end, it is common practice to conduct a preliminary risk assessment to support the optimisation of the operational strategy at a later stage of the project (Baisch et al., 2016). It is crucial to understand how operational parameters such as injection temperature, injection pressure, injection rate, injected fluid volume and distance to existing mapped faults, as well as the local geological conditions, such as rock properties, fault properties and the in-situ stress perturbations, influence the induced seismic hazard potential (Zang et al., 2014). In addition, the impact of temperature variations is becoming more relevant as the cascade use of geothermal resources to optimise energy efficiency becomes more important, leading to lower re-injection temperatures (Rubio-Maya et al., 2015).

The aim of this study is to systematically assess the impact of geological properties and operational parameters and to identify the most influential variables in terms of induced seismic hazard associated with geothermal operations in the Netherlands, with particular attention to the impact of re-injection temperature on dormant faults.

We apply a recently published approach relying on modified Gutenberg–Richter (GR) statistics (Cacace et al., 2021) to estimate the seismic hazard potential associated with the operation of a typical geothermal well doublet in the Netherlands. Our method is based on the simulated perturbations of the frictional Coulomb stress (FCS hereafter) caused by subsurface operations, combined with a thermo-hydro-mechanically coupled reservoir model (Cacace & Jacquey, 2017).

Methods and materials

Induced seismic hazard potential based on FCS variations

In this study, we apply a hybrid model that combines statistical and physics-based elements. We use modified GR statistics for induced seismicity as formulated by Cacace et al. (2021), which builds on the seismogenic index (SI) approach introduced by Shapiro et al. (2007, 2010).

The modified model considers a statistically homogeneous (Poissonian) distribution of pre-existing defects (i.e. faults and fractures) with a bulk concentration N in the porous medium (i.e., the reservoir). They are modelled as non-interacting, point-like defects, and each of them is characterised by a critical value (C) of FCS variation (δFCS). The occurrence of an earthquake at a particular fault depends on whether or not this critical value has been exceeded at this point according to a Mohr–Coulomb failure criterion (S. A. Shapiro, 2018). Furthermore, we assume that the value of C is randomly assigned to each defect in the population from a uniform distribution with a minimum (C _min) and a maximum (C _max) value. In our model, we assume a critically stressed crust by setting the value of C _min to 0.01 MPa, to be higher than the FCS variation caused by the solid Earth tide. The maximum value, C _max, is set to 10 MPa, an estimation of the upper limit of stress drop for fluid injection-induced seismicity based on several case studies (Dinske & Shapiro, 2013).

The stability of a pre-existing fracture or fault is determined by the resolved variations of FCS with respect to an undisturbed tectonic stressing state identified in our model by the tectonic SI (Σ ₀). Positive δFCS values promote instabilities, while negative δFCS values lead to fault stabilisation (Shapiro, 2018; Cacace et al., 2021).

The classical GR equation for injection-induced seismicity is expressed as (Shapiro et al., 2010; Shapiro, 2018):

(1) $$\log _{10} [N_{\geq M}(t)]=\left[{\rm{\Sigma}} _{0}+\delta {\rm{\Sigma}} \left(t\right)\right]-bM$$

where Σ ₀ is the tectonic SI quantifying the seismic reaction of a certain reservoir to a unit volume of injected fluid at a given location, independent of any other operational parameters (e.g. fluid injection rate, temperature or different injection protocols). The term δΣ(t) is a correction term to the SI value accounting for the effects of operational parameters. The b-value is a regional seismicity constant, which provides information about the ratio of large to small earthquakes in a certain region (Dinske & Shapiro, 2013).

Assuming monotonic fluid injection and that only the pore pressure increase is responsible for induced seismic instabilities, δΣ(t) can be expressed as (Shapiro et al., 2010; McGarr, 2014; Van der Elst et al., 2016; Shapiro, 2018):

(2) $$\delta {\rm{\Sigma}} \left(t\right)=\log _{10} [{V_{fluid}}(t)]$$

Cacace et al. (2021) modified Eqs. (1) and (2) by applying a Mohr–Coulomb failure criterion (Labuz & Zang, 2012) and considering variations of FCS, resolved on homogeneously distributed cracks and faults as the triggering mechanism of induced seismicity:

(3) $$\delta {\rm FCS}({{\textbf{x}}},t)=\delta \tau ({{\textbf{x}}},t)-\mu (\delta \sigma ({{\textbf{x}}},t)-\delta p({{\textbf{x}}},t))$$

with x denoting the three-dimensional position vector, δτ, δσ and δp are variations in shear stress, normal stress and pore pressure, respectively and µ is the friction coefficient. The term δσ accounts for poroelastic and thermally induced normal stress changes. Following this approach, δΣ(t) is expressed by the volume integral of resolved δFCS over the model domain:

(4) $$\mathit\Sigma \left(t\right)=\log _{10} \left\lceil \int _{V}{SM[\delta {\rm FCS}(\textbf{x},t)] \over \sin (\psi )}dV\right\rceil$$

where S is the uniaxial storage coefficient, ψ is the friction angle and M[δFCS( x , t)] is the minimum positive monotonic majorant of δFCS( x , t) (Parotidis & Shapiro, 2004).

A majorant is a function that is greater than or equal to a given function, at all points in its domain. The minimum positive monotonic majorant is a special type of monotonically increasing majorant that is defined as the smallest possible majorant that is greater than or equal to the given function and is positive throughout its domain. The role of the minimum positive monotonic majorant function is to ensure that an instability is not induced on the same defect twice unless the previous $\delta {\rm FCS}$ value is exceeded. This approach encompasses the Kaiser effect, which states that seismic events occur only if the previous maximum stress has been exceeded and is a well-known phenomenon in materials science and rock mechanics (Zang & Stephansson, 2009).

Based on Eqs. (1) and (4), the modified GR statistics can be expressed as:

(5) $$N_{\geq M}\left(t\right)=10^{{(\Sigma _{0}}-bM)}\int _{V}{SM\left[\delta {\rm FCS}\left(\textbf{x},t\right)\right] \over \sin \left(\psi \right)}dV$$

Equation (5) generalises the classical GR statistics for induced seismicity to account for any temporal and spatial variation in the effective stresses caused by arbitrary processes (e.g. by the non-linear coupled thermal, hydraulic and mechanical response of the reservoir to complex injection/production protocols). This is critical for studying the long-term-induced seismic hazard potential of a geothermal well doublet where the net injected volume of fluid is conceptually zero. For a more detailed description of the δFCS method, we refer to Cacace et al. (2021).

Numerical model setup

In this study, the finite element code GOLEM (Cacace & Jacquey, 2017) is used to carry out all numerical investigations presented later on in the manuscript. GOLEM is a numerical modelling software designed for simulating thermo-hydro-mechanical processes in fractured porous rocks. It is an object-oriented, scalable, and implicit finite element simulator based on the MOOSE framework (Gaston et al., 2009) that solves the coupled non-linear systems of equations expressing the conservation of fluid mass, solid mass, energy, and momentum.

Our model is based on a previously published calibrated and history matched model of the Groß Schönebeck EGS (Germany), which also targets the Rotliegend reservoir formation (Blöcher et al., 2018). This model was adapted to the tectonic conditions and geological properties of the Slochteren sandstone reservoir in the Netherlands.

The Slochteren Sandstone Formation (SSF), which is part of the Upper Rotliegend Group, was chosen as reservoir rock in our modelling study because it is the most heavily exploited formation by geothermal projects in the Netherlands. It has a large lateral extent, being present in most parts of the country. Currently, there are a total of eight doublets in operation at various locations targeting depths ranging from 1.9 to 2.3 km and temperatures between 70 and 100°C.

The SSF is overlain by the Zechstein Group, which is a relatively complex group composed by a succession of evaporites and carbonates with thin intercalations of claystone. The underlying formation of the SSF is the Carboniferous Limburg Group with a predominantly fine-grained lithology (Buijze et al., 2019; Mijnlieff, 2020).

For the purposes of our study, we have simplified the Zechstein Group to a single, homogeneous rock salt formation which serves as the top cap layer of our model. Likewise, the Limburg Group has been simplified to a homogeneous claystone formation which is used as the bottom cap layer.

All three units are intersected by a single inclined fault, approximated as a planar discontinuity embedded in the three-dimensional rock matrix. A damage zone is defined around the fault with a thickness of 50 m on both sides of the fault plane. In the base model, the permeability of the damage zone is equal to the surrounding matrix. A well doublet is implemented via one-dimensional finite elements, representing the open-hole section of the wells (i.e. from the top to the bottom of the reservoir in our model). The governing equations are homogenised by considering the surface area of the well as a scaling parameter (Cacace & Jacquey, 2017). In the base model configuration, the two wells are located 1 km apart and 250 m from the fault plane. The horizontal extent of the model is 4 × 4 km, which was chosen to avoid edge effects close to the boundaries (Figure S2). The thickness of the reservoir unit is 200 m, while the top and bottom units are 1 km thick to avoid numerical boundary effects inside the reservoir. The top view and side view of our 3D model geometry are schematically shown in Fig. 1, while a snapshot of the mesh is shown in Figure S1. The model geometry together with the mechanical, hydraulic and thermal parameters of the three geological layers is listed in Table 1. The parameter values used in our model are mean values representative of the SSF.

Figure 1. Conceptual geometry of the base model from top view (left) and side view (right).

Table 1. Geometrical and physical properties of the geological units and the fault in the base model.

The boundary and initial conditions are listed in Table 2. A constant temperature boundary condition is applied along the top (47.2°C) and bottom surfaces (115.4°C), while a constant hydrostatic pressure gradient is imposed along the four sides of the model (open boundaries).

Table 2. Summary of boundary and initial conditions used in the base model.

The initial stress state of our model was calibrated by matching the vertical gradients of the principal stresses σ ₁, σ ₂, and σ ₃ characteristic to the SSF (where σ ₁ = σ _V > σ ₂ = σ _H > σ ₃ = σ_h, resulting in a normal faulting regime) at a depth of −2300 m at the centre of the model (x = y = 0 and z = −2300 m). The principal stress gradients in our model are mean values based on available data for the Groningen gas field (Guises et al., 2015; Van Eijs, 2015; TNO, 2015; Osinga & Buik, 2019) and were chosen as ∂σ ₁/∂z= 22 MPa/km, ∂σ ₂/∂z= 15 MPa/km and ∂σ ₃/∂z= 14 MPa/km (Figure S3).

We initialised the horizontal stresses on the western and southern faces of the model with constant displacement boundary conditions, while the northern and eastern faces were fixed (zero-displacement). The constant displacement values were adjusted iteratively until the resulting horizontal stresses at the centre of the model matched the targeted values. The linear vertical stress gradient of 22 MPa/km was imposed across the entire model domain by a constant stress boundary condition applied on the top of the model, while a zero-displacement boundary condition was applied on the bottom.

As we already mentioned in the previous section, we set C _min (the lowest possible critical δFCS value that can be assigned to a fault in our model) to 0.01 MPa. This value was chosen to ensure that any perturbation beyond the effect of the solid Earth tide has the potential to induce a seismic event. This is a conservative approach, as it assumes the presence of critically stressed faults in our model.

The initial and boundary temperature and pore pressure gradients were chosen to be representative of the Dutch geothermal sites in the SSF. All fluid parameters were determined and adjusted based on the fluid properties measured at the Groß Schönebeck site (Blöcher et al., 2018).

After the boundary and initial conditions are established, we simulate the continuous injection and production of a geothermal well doublet over a period of 30 years, with the temperature of the injected fluid set at 30°C and the injection and production rates both fixed at 70 l/s.

In order to evaluate the induced seismic hazard potential caused by the FCS variations resulting from the operation of the well doublet, it is crucial to select a SI characteristic for the SSF. Shapiro (2018) conducted a study to determine the tectonic SI (Σ ₀) for the Groningen gas field. Based on recorded induced seismicity during the production period spanning from 1993 to 2015, the study found that the SI for the Groningen reservoir ranges from −5 to −4. Based on this information, we selected an SI value of −4.5, which ensured that the seismicity rates generated by our base model align with the findings of the case study review conducted by Buijze et al. (2019). This study indicates that geothermal operations in sandstone reservoirs in the Netherlands are unlikely to induce seismic events with magnitudes larger than 2.0.

Finally, we note that Muntendam-Bos and Grobbe (2022) found a spatial variation in the b-value for the Groningen gas field by analysing the seismicity catalogue for the period from 1991 to 2021, with median values ranging from 0.77 to 1.52. However, since a b-value of 0.95 ± 0.04 was determined for the full Groningen catalogue, we adopted a value of b = 1 for our base model, which is a common assumption for tectonic earthquakes (Shapiro, 2018).

Modelling scenarios

A systematic sensitivity analysis was performed to quantify the relative effects of operational and geological parameters in a setup relevant for seismic risk analysis for geothermal projects targeting the SSF. The reference of our sensitivity analysis is the base case model as described in the previous paragraph. All modelling scenarios simulated and analysed in this study are listed in Table 3. In each scenario, the value of only one parameter at a time is changed to the value indicated in Figs. 5 and 6 (to be shown later). All other parameters remain fixed and correspond to those of the base model. The end member values selected for the parameters studied are based on the available literature and were chosen to encompass the range of variations in the SSF across the Netherlands.

Table 3. Summary of the modelling scenarios in the sensitivity analysis performed in this study.

We performed additional analyses to investigate the influence of the SI and the b-value on the induced seismic hazard potential. To determine the sensitivity of the seismic hazard potential to the SI, we varied the SI values between −6 and −3, while leaving the b-value at 1, based on the study of Shapiro (2018). To assess the impact of the b-value on the seismic hazard potential, we varied the b-value between 0.6 and 1.6, while fixing the SI value at −4.5, considering the spatial variations of the b-value in the Groningen gas field discussed by Muntendam-Bos and Grobbe (2022).

Results

Results of the base case model

Figure 2 displays the results from the base model in terms of pore pressure and thermal stress (contour lines) as well as modelled FCS variations (background colours) as extracted from the 3D model along a horizontal plane at 2300 m depth. The model snapshots are taken after 25 days (Fig. 2a, b), when the pore pressure between the two wells reaches steady-state equilibrium (in other parts of the 3D model the pressure is still transient) and the thermal front starts to expand, and at the end of the life cycle of the doublet after 30 years (Fig. 2c, d). In panels (e) and (f), the values for pore pressure, temperature and M[δFCS] are shown as profiles along a line running through the well doublet after 25 days and 30 years of operation, respectively.

Figure 2. Pore pressure, thermal stress and FCS variation in the base model. The panels (a)–(d) show a horizontal plane inside the reservoir at 2300 m depth extracted from the 3D model. The pore pressure and thermal stress are shown by coloured isocontours while the background colour represents the FCS variation. Snapshots are taken at 25 days, when the pore pressure reaches equilibrium between the two wells (in other parts of the 3D model the pressure is still transient) and the thermal front starts to expand, and at the end of the operation after 30 years. In panels (e) and (f), pore pressure, temperature and M[δFCS] values are shown along a line passing through the well doublet at 25 days and 30 years, respectively. Blue dots: injection well; Red dots: production well.

Within the framework of the δFCS model, we compute the temporal evolution of the SI from variations in stress and fluid pressure induced by geothermal processes. Figure 3 shows the magnitude–frequency distribution of the simulated seismicity in the base model. The cumulative number of seismic events above a certain magnitude is represented on a logarithmic scale by different colours for 5-year intervals over the whole 30-year period of operation. The majority of the simulated seismic events are below M_w = 1. The number of induced events increases until the end of the simulated circulation period of 30 years.

Figure 3. Magnitude–frequency distribution of the simulated seismicity in the base model.

Figure 4 displays the computed cumulative magnitude exceedance probabilities. It indicates that there is a 90% likelihood of the predicted induced seismicity being less than magnitude 0.6 after 30 years of operation, and the probability of a magnitude 1.9 event decreases to 10%. However, we must keep in mind that the magnitudes of simulated induced seismic events are determined primarily by the chosen value of the SI and to a lesser extent by the b-value. The absolute magnitudes discussed in this study should be considered in conjunction with the uncertainties shown later in Fig. 8.

Figure 4. The evolution of cumulative magnitude exceedance probability in the base model. The magnitude of the predicted seismic events with a probability of 90% and 10% is shown at the end of year 30.

In the subsequent sections, we present the relative sensitivity and impact of the intrinsic (geological) and extrinsic (operational) parameters examined.

Results of the sensitivity analysis

The main objective of the sensitivity analysis was to investigate the effects of the temperature of the re-injected fluid and the resulting thermal stresses on the induced seismicity, but we also tested a wider range of intrinsic (Fig. 5) and extrinsic (operational) parameters (Fig. 6). Since the thermal effects associated with cold-water injection on induced seismicity are the main focus of our study, we investigated the effects of thermoelasticity in a separate scenario by applying a volumetric thermal expansion coefficient of zero (Fig. 6, S032; Fig. 7), so that no thermal stresses can be induced in the model.

Figure 5. Sensitivity analysis of earthquake occurrence probability for individual model parameters. P90 and P10 denote the probability of 90% and 10%, respectively, for inducing a seismic event of the given magnitude (scenarios S001-S017). The black and red dashed lines indicate the P90 and P10 values for the base case model.

Figure 6. Sensitivity analysis of earthquake occurrence probability for individual model parameters. P90 and P10 denote the probability of 90% and 10%, respectively, for inducing a seismic event of the given magnitude (scenarios S022-S032). The black and red dashed lines indicate the P90 and P10 values for the base case model.

Figure 7. Comparison of cumulative magnitude exceedance probability with and without thermo-poroelastic effects. Dashed curves represent the scenario when thermo-poroelastic effects are neglected, while solid curves show the results of the base model where these effects are considered.

Figures 5 and 6 show the results of the sensitivity analysis, in terms of the maximum moment magnitudes occurring at the end of year 30 with a probability of 90% (P90; left end of the bars) and with a probability of 10% (P10; right end of the bars) compared to the P90 and P10 values of the base model (black and red dashed lines, respectively, in Figs. 5 and 6). For comparison, the P90 (magnitude 0.6) and P10 (magnitude 1.9) values of the base model are marked in Fig. 4. The P90 values of all investigated scenarios range from magnitude 0.1 to magnitude 0.8, while for the P10 value the smallest and largest magnitudes are 1.4 and 2.1, respectively.

The influence of the investigated parameters (within the selected value ranges) is not strong enough to significantly increase the seismic hazard potential associated with the operation of the geothermal doublet. Among the intrinsic properties, the induced seismic hazard potential is more sensitive to the Young’s modulus, the horizontal permeability and the porosity of the reservoir layer, as well as the permeability of the fault. The operational parameters do not have a significantly stronger influence on the induced seismic hazard. However, the P90 and P10 values of all scenarios in this parameter group deviate from the base model, in contrast to the intrinsic properties, where the majority of the tested parameters have no or only a minor influence on the seismic hazard potential

Figure 7 shows the comparison between scenario S032, where thermoelastic effects are neglected (dashed curves), and the base model, where thermoelasticity is considered (solid curves). Initially, in year 1, the two scenarios share similar seismic responses, but as the operational phase progresses, the results of the two models start to deviate, with the base model showing a higher seismic risk. This is indicative of how the induced seismicity is controlled by pore pressure changes on a short time scale (months; Fig. 2a, b, e), while cooling-induced thermal stresses contribute to the system evolution only on longer time scales (years; Fig. 2c, d, f). When thermoelastic effects are neglected, the P90 value is magnitude 0.3, while the P10 value is magnitude 1.6, which is ∼0.3 magnitude lower compared to the base model.

The results presented in Fig. 8 show that higher SI values are associated with higher induced seismic hazard potential, which is reflected in both P90 and P10 values. In contrast, lower b-values are associated with a higher occurrence of large magnitude events in the seismic catalogue, resulting in higher P90 and P10 values, while higher b-values are associated with a higher occurrence of smaller magnitude events and a lower occurrence of larger magnitude events. These results are consistent with the observations of Dinske and Shapiro (2013), who suggest that local tectonic conditions have a significant influence on induced seismicity and strongly constrain the seismic hazard potential.

Figure 8. Influence of seismogenic index (SI) on the probability of earthquake occurrence in the base model. P90 and P10 denote the probability of 90% and 10%, respectively, for inducing a seismic event of the given magnitude.

We note that, to interpret absolute magnitude values of induced seismicity in an area of interest, it is of utmost importance to calibrate the hazard assessment curves against site-specific observed magnitude–frequency distributions. Such calibration necessitates the acquisition of seismic monitoring data during the stimulation phase of the particular geothermal project under investigation.

Discussion

A sensitivity analysis was performed on a generic model of a geothermal reservoir in the SSF with a single regional fault and exploited by a well doublet. We applied modified GR statistics on the simulated frictional Coulomb-stress perturbations with a fixed SI of −4.5.

As shown in Fig. 7, our models indicate that the seismic hazard is mainly controlled by the build-up of pore pressure due to fluid circulation. Even when thermal stresses are disregarded, the magnitude exceedance probabilities continue to increase throughout the operation. The spacing between the curves indicates that the pressure field is gradually approaching a quasi-stationary state (Fig. 7, dashed curves). However, when thermal stresses are also considered (Fig. 7, solid curves), the seismic hazard increases due to the continued temperature changes after the pressure equilibrium is reached (after ∼1 month; Fig. 2). This results in higher induced seismicity in terms of P90 and P10 values, with an increase of approximately 0.3 magnitude units.

In the following, we discuss which of the parameters have the greatest influence on seismic hazard. We interpret the effect of these parameters by associating them with changes in the pore pressure field and the thermal stresses.

The geometrical parameters have little or no effect on the induced seismicity. Out of the geomechanical parameters, only the Young’s modulus influences the P10 and P90 values, which is expected as the thermal stress magnitude is linearly dependent on this parameter. As such, a higher Young’s modulus corresponds to a greater potential for induced seismic hazard (Segall & Fitzgerald, 1998).

The hydraulic properties that have a significant impact on induced seismicity are the horizontal permeability and porosity of the reservoir, as well as the permeability of the fault damage zone. The porosity of the reservoir is directly linked to the uniaxial storage coefficient, which is shown in Equation (4) to have a direct effect on induced seismicity. As the porosity of the reservoir increases, so does the uniaxial storage coefficient, resulting in a higher potential for induced seismic hazard.

While a lower horizontal reservoir permeability, as shown in Fig. 5 (S009a), may also increase the seismic hazard potential, this assumption is only valid if the injection rate is fixed. In practice, the operator may have to reduce the injection rate if the reservoir has low permeability, which could potentially lower the seismic hazard. Conversely, a high permeability damage zone (as shown in Fig. 5, S015b) can help to connect the reservoir formation to the bounding cap rock formations and maintain hydrostatic pressure throughout the model. This prevents significant pressure build-up within the reservoir, reducing the seismic hazard. The high permeability damage zone also enables the cold-water plume to spread preferentially along the fault (Fig. 9), which limits fluid infiltration into the upper and lower cap layers and prevents the build-up of thermal stresses along the reservoir boundaries. This further contributes to reduced seismicity.

Among the thermophysical parameters, the coefficient of thermal expansion has the greatest impact on induced seismicity, as it determines how much thermal stress can accumulate for a given temperature difference. The heat capacity of the reservoir rock has a minor impact on induced seismicity, while the effect of thermal conductivity is negligible. This is because the Slochteren base case system is convection dominated due to the high permeability of the sandstone.

Operational parameters have a more significant influence on induced seismicity than intrinsic properties, but their variations do not lead to a significant increase in the seismic hazard. The temperature of the injected fluid affects induced seismicity as expected; lower temperatures result in higher thermal stresses and increased seismic hazard potential. However, reducing the injection temperature from 30°C to 15°C does not pose a significantly higher seismic hazard. The distance between the two wells (borehole spacing) affects the amount of cold fluid infiltrating the caprock and bedrock and therefore affects the build-up of thermal stresses. When the two wells are close enough to each other, the cold-water plume propagates preferentially towards the production well rather than vertically towards adjacent strata, which results in lower seismic hazard potential. Conversely, when the boreholes are farther apart, diffusion has a stronger influence and a larger amount of cold water can migrate into the caprock and bedrock, causing thermal stresses to build up, and higher seismic hazard potential.

The induced seismic hazard potential is proportional to the injection rate, as a higher injection rate leads to higher pore pressure around the injection well and a larger cold-water plume. The reservoir is primarily stimulated by mode I tensile fractures when higher injection rates are used. Designing low injection rates in the field can be an appropriate measure to increase the overall size of the fluid infiltration zone. This will allow significant parts of the reservoir to hydro-shear, in particular, at pre-existing fracture networks. Hydraulic shearing will allow to enhance permeability in a sustainable, permanent way since relocated asperities along opposite fracture walls will create residual pore space beneficial for fluid pathways (Rinaldi & Rutqvist, 2019).

Some of the parameters do not affect induced seismic hazard potential in this study, but this conclusion holds for the generic model presented, only. We present a case study specific to the SSF, and these intrinsic parameters are also characteristic of the same formation. The parameter values chosen in our sensitivity study represent the uncertainties in the available data, so it is possible that their influence is relatively small for our study area but may still be significant for a different reservoir.

It should be noted that the tectonic SI value determined by Shapiro (2018) for the Groningen gas field is only representative of the reservoir volume that hosted the seismic events used in the analysis to obtain it. It is possible that the induced seismic hazard potential is underestimated if basement faults are also activated by geothermal operation, as these are typically characterised by a higher tectonic SI (Dinske & Shapiro, 2013). Assuming a constant SI = −4.5 for all three geological units may not be appropriate in this scenario. According to Buijze et al. (2019), clay and shale layers interbedded with sandstone formations can serve as a hydraulic barrier, effectively isolating the formation from the seismogenic basement. In our study, we make the assumption that the Limburg Group, which underlies the SSF, plays this role as a hydraulic barrier. Therefore, we exclude the possibility of basement fault reactivation resulting from fluid injection in our analysis.

Another limitation of the present study is that the stress-dependent permeability variations of the faults and fractures were not considered. Cao et al. (2022) have shown that the cooling-induced permeability enhancement is the dominant mechanism for induced seismicity in the Hellisheiði geothermal field in Iceland. It is important to note, however, that the Hellisheiði geothermal field is a high enthalpy geothermal reservoir in fractured volcanic rock. It is therefore not appropriate to transfer this observation to the low enthalpy, porous Slochteren sandstone reservoir investigated in this study, where permeability variations are expected to be minimal. Nevertheless, the cooling-induced permeability enhancement is considered as subject for a future study.

Based on our study, it can be concluded that if the geology of a certain site is well characterised by an exhaustive exploration campaign and the area is not tectonically (i.e. seismically) active, then minor changes in the intrinsic parameter values would not cause a significant increase in the background seismic hazard.

On the other hand, by tuning the operational parameters it is possible to further control the induced seismicity on the long term. In summary, our study shows that while the build-up of pore pressure is the dominant factor in induced seismic hazard during the early stages of a geothermal well doublet’s operation, thermal stress must also be considered in the risk assessment for longer production periods.

Conclusion

We built a finite element reservoir model of the SSF with a single planar fault discontinuity and simulated cold-water injection and hot water production through a geothermal well doublet over 30 years. We performed a sensitivity analysis of different intrinsic (local geology) and extrinsic (operation of the well doublet) parameters on the induced seismic hazard using a modification of the GR statistics for induced seismicity, with a special emphasis on the connection of pore pressure vs. thermally induced stresses. The main findings are as follows:

1. 1. Our results suggest a competing mechanism between the pore pressure front and the thermal front propagating due to cold-water injection. We found that on short time scales (months), thermally induced stresses have no influence on the induced seismicity, as the thermally disturbed zone is restricted to the close vicinity of the injection well at this early stage of fluid injection. This observation is consistent with the traditional practice of relating fluid injection-induced seismicity only to the total injected fluid volume. After years of fluid circulation, the thermal front continues to expand, resulting in a significant contribution of thermally induced stresses to the induced seismicity rate (approximately +0.3 M over 30 years). Therefore, we argue that thermally induced stress changes should be considered in the risk assessment of long-term geothermal operations.

2. 2. We found that the pore pressure reaches a quasi-stationary state between the two wells after ∼1 month. The pressure front slowly continues to expand until the end of the operation, elevating the pore pressure inside the reservoir and partly in the more rigid caprock and base rock, too. This makes the pore pressure increase the dominant contributor to the seismic hazard in our model. The thermally induced stresses, however, further increase the potential induced seismic hazard.

3. 3. Our sensitivity analysis has shown that the most critical geological parameters in the risk assessment are ranked as follows:

   1. a. Young’s modulus of the reservoir and the adjacent formations

   2. b. Horizontal permeability of the reservoir

   3. c. Porosity of the reservoir and the adjacent formations

   4. d. Fault permeability and porosity

   5. e. Bulk modulus of the matrix

4. 4. Apart from the SI, we found that lowering the injection rate (and the production rate) can significantly reduce the induced seismic hazard potential, even though the total injected volume in our model should be conceptually zero (as injection rate = production rate). Decreasing the temperature of the re-injected fluid from 30°C to 15°C does not increase the seismic hazard significantly over 30 years of fluid circulation. This can be a relevant observation for a more efficient, cascade use of geothermal heat.

Figure 9. A cross-section, perpendicular to the fault, showing the temperature and pore pressure distribution around the injection well after 30 years of operation (a) in the base case scenario; (b) in scenario S15b, with a high permeability damage zone around the fault. The pore pressure is displayed by isobars while the background colouring shows the temperature distribution. Red line: fault; magenta lines: reservoir top and bottom; blue line: injection well.