Transient evolution of the relative size distribution of earthquakes as a risk indicator for induced seismicity

ARTICLE https://doi.org/10.1038/s43247-022-00581-9 OPEN Transient evolution of the relative size distribution of earthquakes as a risk indicator for induced seismicity 1234567890():,; Vanille A. Ritz1 ✉, Antonio P. Rinaldi 1 & Stefan Wiemer 1 Induced earthquakes pose a substantial challenge to many geo-energy applications, and in particular to Enhanced Geothermal Systems. We demonstrate that the key factor controlling the seismic hazard is the relative size distribution of earthquakes, the b-value, because it is closely coupled to the stress conditions in the underground. By comparing high resolution observations from an Enhanced Geothermal System project in Basel with a loosely coupled hydro-mechanical-stochastic model, we establish a highly systematic behaviour of the bvalue and resulting hazard through the injection cycle. This time evolution is controlled not only by the specific site conditions and the proximity of nearby faults but also by the injection strategy followed. Our results open up new approaches to assess and mitigate seismic hazard and risk through careful site selection and adequate injection strategy, coupled to real-time monitoring and modelling during reservoir stimulation. 1 Swiss Seismological Service, ETH Zürich, Zurich, Switzerland. ✉email: COMMUNICATIONS EARTH & ENVIRONMENT | (2022)3:249 | https://doi.org/10.1038/s43247-022-00581-9 | www.nature.com/commsenv 1 ARTICLE T COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-022-00581-9 he frequency and detection of induced earthquakes has increased dramatically over the past few decades, leading the scientific community and industry to become more aware of the issue. Most sources of induced seismicity involve the underground injection or withdrawal of fluid1. Earthquakes have been associated with several human activities ranging from oil and gas extraction2–4 to underground wastewater storage5,6, and even artificial lakes for hydro-power7, tunnelling8,9 and nuclear waste disposal10. Among the various industrial operations, geothermal exploitation constitutes a promising approach for green energy production, but often requires fluid stimulation to increase the reservoir’s permeability to create Enhanced Geothermal Systems (EGS11,12). Such stimulation activities, however, are often linked to strong seismic events. Numerous cases around the world have led to much debate in the scientific community and industry as to whether this technology is really beneficial. Renowned examples are: M 2.9 in Soultz-sous-Forêts (2003, France)13, ML 3.4 in Basel (2006, Switzerland)14 and Mw 5.5 in Pohang (2017, Republic of Korea)15; but other geothermal exploitation schemes are not exempt from damaging induced earthquakes, as a recent sequence in Strasbourg culminating in an MLv 3.6 shows (2019–2021, France)16. The occurrence of induced earthquakes of potentially large magnitude requires the development of strategies and tools to mitigate the hazard posed by EGS stimulations17,18. Numerical models constitute the basis for the most advanced tools proposed to understand, forecast, and mitigate the hazard and risk accompanying geothermal stimulations. Multiple classes of models have been used to investigate fluid induced seismicity19, from the purely statistical20–22 to the fully coupled thermohydro-mechanical23–25. Hybrid models strive to combine the strength of both approaches, using the velocity of statistical methods, and the complexity of the principal physical processes26–31. A fundamental aspect of induced seismicity mitigation relies on linking modelling results and risk/hazard calculations32–34. Numerical models can also investigate the impact of injection scenarios on the risk posed by induced seismicity31,35–37. In this study we use the capabilities of a loosely coupled hydromechanical-stochastic model – TOUGH2-Seed38 – to investigate how seismicity is affected by the presence of a major fault zone and to understand if a different injection strategy may help with predicting the potential hazard. To assess an injection scenario, we follow a simple metric: the evolution of the Gutenberg–Richter b-value39, describing the relative size distribution of earthquakes, during both injection and post-shut-in phases. This approach allows us to go beyond statistical models that aim to simulate the rate/number of events (λ or the Gutenberg-Richter a-value)32,40,41 by modelling both the a-value and b-value. The approach is based on the assumption that the bvalue is a proxy for the state of stress in the subsurface: indeed, several studies have shown that the b-value is sensitive to the state of stress (inverse dependence observed both in laboratory42–44 and field45,46 studies). Modelling both the number of events and the b-value provides an additional proxy for seismic risk47,48. Indeed, the b-value is the true driver behind the level of risk/hazard as highlighted in Fig. 1a: for a comparable a-value, Pohang and Basel have vastly different probabilities of a magnitude 6 occurring. Usually, the probability of occurrence of a large event is used with a fixed bvalue calculated for a region, as the b-value shows variability depending on the tectonic context (Fig. 1a). However, more recent studies with high-resolution catalogues, both in natural contexts49,50 and in induced seismicity contexts5,51, have shown that the b-value changes through time and space within a given 2 region due to additional forcing (e.g. other earthquakes or fluid injection). These changes in the b-value (both spatial and temporal evolution) have been shown to help with the forecasting of larger events: the b-value tends to increase after a main shock and during the aftershock sequence49. However, it has been observed that in the case of doublets, the second event can be preceded by a decrease in the b-value (compared with the background) and a Traffic Light System based on b-value variation to understand when the hazard has passed has been suggested50. Thus, the bvalue can be used as a proxy for average stress conditions and as a tool for risk mitigation to discriminate between a typical aftershock sequence and a precursory sequence. In the case of induced seismicity in Basel, initial studies showed how the b-value dropped significantly after shut-in21, but with a high-resolution catalogue51, the b-value is actually observed to drop days before the shut-in and the two largest events (ML 2.9 during injection and ML 3.4 a few hours after shut-in14). Figure 1-b shows the evolution of the b-value in time in Basel using a high-resolution catalogue51 with a fixed time windows of 0.2 days. Supplementary Fig. S1 presents a more complete analysis of evolution of the bvalue, number of events above completeness and magnitude of completeness in time for the Basel case. The variability in the bvalue during the injection and a sharp drop around day 3 have been thought to be linked to pore-pressure perturbations and a change in the growth behaviou (...truncated)