Effects of injection rates and protocols
Observations derived from field applications and laboratory studies indicate a correlation between induced-seismicity rates and fluid-pressurisation rates17,18. In this paragraph we show that rates of (mathrm{delta FCS}) map directly onto seismicity rates (Eq. 6), thereby enabling to control for the probability of the occurrence of a seismic event at a given time. More information about the model set up, boundary and initial conditions are provided in Supplementary Information 2.
Figure 2 reveals how injection protocols can affect the induced seismicity, both in terms of rates and magnitudes. In line with existing observations, we find that seismicity rate is proportional to injected pore pressure and—most importantly—to imposed pressurisation rates. In our formulation, the dynamic triggering of seismic events in space and time is dictated by the ability of the reservoir to relax induced variations in FCS, and these variations also control post-injection seismicity rates. Accordingly, modelled decay in seismicity rates at shut-in can be fitted by a modified Omori’s law32 in the form of (see also Supplementary Information 4):
$${N}_{left(ge Mright)}left(t>{t}_{0}right)= frac{1-({frac{t}{{t}_{0}})}^{left(1-pright)}}{p-1}{N}_{left(ge Mright)}left(t={t}_{0}right),$$ (8)
where ({N}_{(ge M)}(t>{t}_{0})) is the number of cumulative events after shut-in, ({N}_{(ge M)}(t={t}_{0})) describes the number of cumulative events at shut-in and p denotes the exponent of the power law distribution.
Figure 2 Effects of injection protocols on induced seismicity rates and magnitudes. (a) Injection rates (red curves) and computed wellhead overpressure (blue curves) for the tested injection protocols. (b) Resulting seismicity rates and the power-law exponents (Eq. 8) that best fit the post-injection decay curves (empty circles). (c) Cumulative exceedance probability of a given magnitude Mw (Eq. 9). Each scenario lasted 6.5 days and used the same amount of net injected fluid volume (10,800 m3) and hydraulic energy (Figure S1). Different cases (upper and middle panels) have been shifted in time to ease visual comparison. Full size image
When considering monotonic injection, higher rates increase the overall seismic potential and lead to a steeper seismicity-rate curve during stimulation—followed by a more-gradual decay at shut-in (smaller power-law exponent, p)—and therefore to a higher post-injection seismic risk.
The magnitude distribution also depends on imposed injection rates, with higher rates being more likely to produce higher-magnitude events10. This is shown in the lower panel of Fig. 2 where we plot the obtained induced earthquake probability distributions for each protocol. The latter probabilities are computed by assuming that the occurrence of a given-magnitude earthquake can be described by a Poisson process. Therefore, the probability of exceeding a given magnitude (Mw) is given by:
$$mathrm{w}left({mathrm{M}}_{mathrm{w}}right)= 1-mathrm{ w}left(0,{mathrm{M}}_{mathrm{w}},{mathrm{N}}_{ge mathrm{M}}right)= 1 -mathrm{ exp}left(-{mathrm{N}}_{>mathrm{M}}right),$$ (9)
where (mathrm{w}(0,{mathrm{M}}_{mathrm{w}},{mathrm{N}}_{ge mathrm{M}})) is the probability of having no earthquakes of a magnitude greater than or equal to ({mathrm{M}}_{mathrm{w}}) within the specified time given the expectation of having ({mathrm{N}}_{ge mathrm{M}}) earthquakes whose magnitudes are greater than or equal to ({mathrm{M}}_{mathrm{w}}).
Our results demonstrate that it is possible to reduce the seismic risk during both active stimulation and shut-in by controlling rates and pressure magnitudes via stepwise-controlled injection or a cyclic-stimulation strategy. The role of cyclic-injection strategies or step-rate tests in influencing seismic-hazard potential depends on (1) the duration of each cycle, (2) imposed injection- and pressurisation rates, and (3) the time needed for the induced overpressure to equilibrate between each stage (Fig. 2).
With respect to the post-injection phase, cyclic rock stimulation is characterised by significantly faster relaxation upon shut-in and therefore by higher decay rates of induced seismicity (i.e. higher p values) than is monotonic or step-rate fluid injection. Thus, given equivalent tectonic conditions (the same ({Sigma }_{0})) and the same injected volume and hydraulic energy (Figure S1), the seismic hazard of a cyclic stimulation is significantly smaller (especially during the post-injection phase) than the hazard of monotonic injection strategies.
Deriving generic conclusions on the influence of the discussed stimulation protocols (which should be evaluated for individual cases) lies beyond the scope of this study. However, our results—which have been corroborated by novel laboratory17,18 and field evidence33—demonstrate the importance of systematically controlling the injection protocol to minimise the seismic hazard induced by subsurface stimulation.
Role of thermo-poroelastic stress transfer during long-term circulation of geothermal fluids
The dynamics of fluid-induced seismicity have been traditionally limited to the direct effects of pore pressure increase within the stimulated rock volume, and the effects of thermo-poroelastic stress transfer have only recently grown in importance14,34,35,36,37,38. These additional mechanisms are particularly relevant for geothermal applications over relatively long time scales of concurrent cold-water injection and hot-water production.
In our approach we can handle both long-range (quasi-)static thermo-poroelastic and direct pore-pressure-effect interactions over time and we can effectively translate these interactions into seismicity-rate changes and into the probability of the occurrence of a seismic event of a given magnitude at a given time. To further understand how thermo-poroelastic stress transfer can affect seismic hazard, we simulated a long-term circulation scenario (via simultaneous injection and production) and analysed the seismogenic stability of the system. We targeted a geothermal doublet lifetime of 30 years and continuously re-injected geothermal fluids (at 70 °C at 50 l/s) into the reservoir (at 150 °C).
We identified two main stages in the evolution of the system: (1) a short-term regime (before pore pressure equilibrates), during which instability is driven by the induced pore pressure increase near the injection well and variations in FCS are limited by the fluid-pressurised front, and (2) a long-term regime (after pore pressure equilibration), which shows a dominance of thermo-poroelastic effects (Figs. 3, 4). During this stage, thermo-poroelastic stress variations control the geomechanical response of the reservoir and result in a larger spatial footprint of positive variations in FCS that extend beyond the pressurised front (Fig. 3).
Figure 3 source and production-sink. During the early stages of the system evolution, variations in FCS are bounded by the pore pressure front and show an abrupt decay away from the injection boreole that is driven by pore pressure diffusion only. Once pressure equilibrium is reached in the reservoir, variations in FCS are driven by changes in thermo-poroelastic stresses, thereby resulting in a larger seismic footprint and a more-gradual spatial decay away from the injection source. Thermo-poroelastic stress transfer and induced variations in FCS during long-term geothermal operations. Upper panels: computed variations in FCS for a thermo-poroelastic simulation of a 30-year continuous injection and production. Left: model results at the time when pore pressure equilibrated in the system (maximum extent of the pore pressure front). Right: model results at the final stage of the operations. Background colours indicate computed values of (mathrm{delta FCS}), and isocountours of computed pore pressures (left panel) and thermal stresses (right panel). Lower panels: variations in computed (mathrm{rm M}[mathrm{delta FCS}]) extracted from the 3D model along a line passing through the injection- Full size image
Figure 4 Cumulative exceedance probability of a given magnitude Mw (upper panel) and its evolution as a function of the total circulated volume (lower panel) in a geothermal well doublet during 30 years of continuous injection and production. Continuous curves indicate cases for which thermo-poroelastic effects are considered, and dashed curves indicate cases for which these additional effects are disregarded. The difference between the curves increases with time and peaks after pore pressure equilibrates in the system. Full size image
Variations in the stress state due to coupled thermo-poroelastic stress also affect the seismicity distribution (Fig. 4). By disregarding thermo-poroelastic effects, the forecasted annual probability of exceedance (Eq. 9) is biased towards lower magnitudes for all times and/or circulated fluid volumes. The difference between the two models systematically increases with the circulated volume, thereby suggesting that cold-fluid re-injection can induce larger events than pore pressure diffusion alone, especially at later times in the system evolution. Computed differences concern mostly the highest magnitudes in the explored range (Mw > 2), thereby illustrating the significant risk of the long-term circulation of geothermal fluid, which must be accounted for in any seismic-risk-assessment- and mitigation strategy.
Modelling of the 2007 hydraulic stimulation at the geothermal site of Groß Schönebeck, Germany
In this section, we apply our model to a real-world case study and address a relatively complex stimulation experiment performed at the scientific geothermal facility of Groß Schönebeck in Northern Germany22. This 5-day water-frac treatment targeted the volcanic-rock section of the reservoir and consisted of five distinct cyclic-injection phases, with maximum flow rates of up to 9 m3/min (Fig. 5). In total, 13,170 m3 of slickwater and 24.4 tons of meshed-quartz sand were injected during the treatment while maintaining a wellhead pressure below 58.6 MPa (Figure S2). During the last two cycling stages, micro-seismic events with moment magnitudes ranging between − 1.8 and − 1.0 were recorded by a downhole seismometer located close to the injection section of the well. The hypocentre distribution exhibited a degree of directional ordering, with seismic activity migrating away from the injection well with time. The relocated seismic catalogue was found to be consistent with the presence of a seismic fault plane that was favourably oriented for re-activation under the in-situ stress state39 (Figure S3).
Figure 5 Modelling the 2007 reservoir treatment at Groß Schönebeck22. Upper Panel: imposed injection rates (blue curve) as well as measured (black curve) and computed (grey curve) bottom-hole overpressure ((mathrm{Delta BHP})) for the 2007 stimulation experiment carried out in Groß Schönebeck. Lower panel: predicted (continuous curves) maximum magnitudes and seismic moment as a function of the total injected volume. Dashed curves are predictions derived from Eq. (14) in10 for different values of the tectonic seismogenic index (({Sigma }_{0})) that are considered reasonable for the Groß Schönebeck site26. The dashed rectangle encloses the recorded micro-seismicity. Full size image
The aim of this modelling example is to showcase the forecast capability of our approach when applied to a complex geological set-up and injection protocol by also considering additional effects from the nonlinear physics on the resulting seismic response of the reservoir. For this purpose, we set up a 3D thermo-hydro-mechanical model of the stimulation treatment, which integrates details of the geological architecture of the reservoir, including mapped major fault zones, the open-hole section of the well, as well as the stimulated fracture. The hydraulic stimulation was considered implicitly by imposing variations in the fracture aperture as computed by a frac model22, which were also used to compute the evolution through time of the fracture permeability (Figure S2). Supplementary Information 3 provides a detailed description of the model set-up, including rock-properties and imposed boundary and initial conditions. The overall goal of the simulation was to match the spatio-temporal evolution of monitored micro-seismicity as well as their magnitude distribution. We identify a seismic event at a specific location following a Mohr–Coulomb-type rupture model, where rupture occurs if resolved stress magnitudes are higher than the critical yield strength of the rock. By assuming a constant friction coefficient, such critical value can be recast in term of an adimensional parameter, the slip tendency given by the ratio of shear to normal loading39. It follows, that an induced event will occur if the value of the slip tendency exceeds the value of the imposed static friction coefficient (considered equal to 0.6 in this study).
Model validation was done by a history matching against the bottom-hole pressure (BHP) measured during the treatment (Fig. 5, upper panel). We then used the results from the best-fitting model to compute the spatio-temporal evolution of the micro-seismic events and expected variations in computed magnitudes (Fig. 5; Figure S3). In Figure S3 we present the results of the best fitting model, in terms of the hydro-mechanical state of the reservoir after the fourth stimulation stage. We computed a maximum fluid-pressure increase of approximately (Delta mathrm{p}=) 8.1 MPa along the seismic plane, higher than the level required to drive instability based on in-situ stress conditions ((Delta mathrm{p}=) 6.2 MPa). We also observed no appreciable temperature changes during the entire stimulation and therefore concluded that for this stimulation, a pore pressure increase that was enhanced by poroelastic stress transfer triggered the reactivation of the fault, thereby triggering the observed micro-seismicity. Indeed, computed values of slip-tendency at the locations of the recorded hypocentres are in line with the fault instability in these locations and also match the monitored temporal pattern, with events that show an upward propagation along the seismic plane (Figure S3; Movie S1).
The model is also able to reproduce the upper bound of observed magnitudes of the induced events (Fig. 5 lower panel). Relying on a classical log-linear scaling between the total injected volume and maximum expected magnitudes9,10, we systematically overestimate the magnitude limits of the induced microseismicity, which occurred during the last two injection stages in Groß Schönebeck (dashed curves in Fig. 5). On the contrary, our model can capture the observed breakdown of the log-linear scaling between the injected volume and thereby to match the limits of induced event magnitudes. We interpret this result as a stress-memory effect of the reservoir caused by the operational control on the maximum wellhead pressure (Figure S2), which is similar to the Kaiser effect often discussed in global geothermal systems28,33.
https://www.nature.com/articles/s41598-021-02857-0
