Seismic Moment Evolution During Hydraulic Stimulations

Analysis of past and present stimulation projects reveals that the temporal evolution and growth of maximum observed moment magnitudes may be linked directly to the injected fluid volume and hydraulic energy. Overall evolution of seismic moment seems independent of the tectonic stress regime and is most likely governed by reservoir specific parameters, such as the preexisting structural inventory. Data suggest that magnitudes can grow either in a stable way, indicating the constant propagation of self‐arrested ruptures, or unbound, for which the maximum magnitude is only limited by the size of tectonic faults and fault connectivity. Transition between the two states may occur at any time during injection or not at all. Monitoring and traffic light systems used during stimulations need to account for the possibility of unstable rupture propagation from the very beginning of injection by observing the entire seismicity evolution in near‐real time and at high resolution for an immediate reaction in injection strategy.


Introduction
The effort to reduce production of greenhouse gases, especially in Europe, has fueled the search for clean and renewable energy sources. Geothermal energy has long been investigated as a potential complement and long-term replacement for traditional fossil fuels in electricity and heat production, with its significant baseload capacity. However, unlike in hydrothermal systems in regions of elevated heat flow, projects aiming to extract fluids hot enough to produce electricity need to reach large depths, typically >4,000 m in crystalline basement or sedimentary basins. In order to develop deep geothermal reservoirs in low-permeability rocks, the formation needs to be hydraulically stimulated. Creation of Enhanced Geothermal Systems (EGS) opens fluid flow paths by injection of large quantities of water, which is typically accompanied by seismicity (Ellsworth, 2013;Majer et al., 2012).
Large induced earthquakes have led to the termination or suspension of several EGS projects, such as the deep heat mining project Basel in Switzerland (Giardini, 2009) and St. Gallen (Diehl et al., 2017). The occurrence of a M W 5.5 earthquake in 2017 near Pohang, South Korea, has been linked to a nearby located EGS project (Ellsworth et al., 2019). In densely populated urban areas substantial public concern about EGS projects now exists and is mainly related to stimulation-induced seismicity. In contrast, two large-scale EGS projects in Australia have been operational for years, and their remote location has prevented any acceptance issue arising from induced seismicity (Albaric et al., 2014;Baisch et al., 2006). However, most recently Kwiatek et al. (2019) reported on the successful and safe stimulation of the world's deepest EGS project so far, located near Helsinki. The authors attributed their success to the high-precision, almost real-time monitoring and analysis of seismic data feeding into a traffic light system and guiding the adaptive stimulation strategy .
Nucleation and propagation of earthquake ruptures related to fluid injection pose fundamental questions related to understanding and controlling injection-induced seismicity. Early models such as McGarr (1976) assumed that volume changes induce local changes in deviatoric stresses, which are then relaxed by induced earthquakes. Shapiro et al. (2010) defined a seismogenic index characterizing the level of seismic activity expected from injecting fluids into the rock formation, expanding on earlier work on the occurrence probability of events above a certain magnitude threshold (Shapiro et al., 2007). McGarr (2014) defined a critical pore pressure change caused by fluid injection inducing seismic events. The predicted cumulative seismic moment and maximum event magnitude are expected to increase with total fluid volume injected. Shapiro et al. (2011) postulated that the expected maximum magnitude may be related to the minimum axis of an ellipsoid covering the spatial distribution of earthquake hypocenters representing the fluid-stimulated volume. However, Norbeck and Horne (2018) recently analyzed numerically potential factors influencing rupture growth beyond the stimulated volume affected by pressure perturbations. Also, Galis et al. (2017) used a fracture mechanics model to analyze the relation between injected fluid volume and stable rupture propagation. The authors differentiate between stable self-arrested ruptures and unstable "runaway" ruptures, which propagate out of the stimulated volume on prestressed faults. Furthermore, Galis et al. (2017) found a fracture mechanics-based scaling relationship between maximum expected magnitude of self-arrested events and injected volume. These models predicting maximum magnitudes of induced events implicitly assume a stable rupture propagation process (Galis et al., 2017). In contrast, van der Elst et al. (2016) argued that induced earthquakes occurring along tectonic faults favorably oriented with respect to the tectonic stress field are unstable. The maximum expected magnitude is then only limited by regional tectonics and fault connectivity, which also holds for induced earthquakes.
Here, we compare the evolution of different parameters derived from seismic catalogs with injection parameters from several EGS projects and one scientific drilling project (Table 1). Unlike what was previously done, we treat each project as a dynamic experiment, rather than constraining the analysis to the final arrested state. The seismicity evolution with progressive fluid injection from the stimulations is compared to current models of injection-induced seismicity during the entire stimulation period. Our analysis provides insight about how seismic moment and related parameters of induced events evolve during injection and how these parameters are affected by tectonic faulting regime.

Data and Methods
We analyzed data collected from different studies investigating a total of eight stimulation projects, all located in crystalline basement rock (Table 1). These include the most prominent European EGS projects in Basel, Switzerland (BAS) (Häring et al., 2008) and Soultz-sous-Forêts (STZ), France. In Soultz, three different stimulations over the course of 10 years were performed in different wells and different depths. Therefore, we differentiate between the injections in 1993 (STZ93), 2000 (STZ00), and in 2003 (STZ03) (EOST & GEIE EMC., 2017, 2018a. We also included the deepest EGS Project to date (St1), located in Helsinki, Finland (Kwiatek et al., 2019). Furthermore, we included the fluid injection experiment from the German super deep scientific drilling hole (KTB) (Baisch et al., 2002), two Australian EGS projects, located at Paralana (Para) (Albaric et al., 2014) and the 2003 Cooper Basin (CBN) injection (Baisch et al., 2006), and the EGS project near Pohang, South Korea (Ellsworth et al., 2019). Finally, we also considered a single well injection period at the Berlín geothermal field (BGF), El Salvador, representing the only hydrothermal site considered here (Bommer et al., 2006;Kwiatek et al., 2014). The projects are located at sites dominated by strike-slip (e.g., KTB), normal faulting (e.g., Soultz) and transpressive (e.g., St1) tectonic stress regimes in a depth range between 3.5 and 9.0 km ( Table 1).
The main parameters analyzed in this study are the maximum observed moment magnitude during and following the stimulation M max,obs and cumulative seismic moment M 0,cum with injected fluid volume and hydraulic energy E hyd , that is, the product of well head pressure and injected volume integrated over time. We correct M 0,cum for each sequence accounting for missing events caused by different magnitude of completeness using computed synthetic catalogs (Figure 1b; see supporting information for more details). Furthermore, we compare seismic injection efficiency I E , that is, the ratio of cumulative radiated energy and hydraulic energy for the different stimulation projects (Maxwell et al., 2008). The radiated seismic energy E 0 for each event is estimated from the seismic moment M 0 and stress drop Δσ following Hanks and Kanamori (1979) as where G is the shear modulus calculated as G ¼ ρV 2 s , with density ρ and shear wave velocity V s , and η R is the radiation efficiency. When available we used already calculated catalogs of radiated energy (Kwiatek et al., 2014;Kwiatek et al., 2019) or computed radiated energy using average stress drops available (Baisch et al., 2009;Charléty et al., 2007;Goertz-Allmann et al., 2011;Jost et al., 1998) and assumed a median radiation efficiency of 0.46 (McGarr, 1999). The KTB data were corrected for the finite bandwidth applied in previous studies (Ide & Beroza, 2001). For the Pohang and Paralana stimulations no stress drop or radiated energy calculations were available, and so we assumed an average stress drop of 3 MPa (Cocco et al., 2016).

Maximum Magnitude With Injected Volume
For most of the respective stimulation period, all analyzed data sets show an increase of maximum observed moment magnitude M max,obs with injected volume that roughly fits to a linear trend in a double logarithmic plot ( Figure 1a). However, the slopes of the respective evolutions differ significantly between sites. The increase of M max,obs with injected fluid volume follows a slope between 1 and 1.5 as predicted by McGarr (2014) and Galis et al. (2017), respectively, for most of the projects analyzed here. The trends observed for KTB, Paralana, and St1 projects are best predicted by the model of Galis et al. (2017). In contrast, Soultz93 and Soultz00 stimulations show a slope, which is better modeled by McGarr (2014) with an even smaller slope for BGF. In striking contrast, stimulations performed at the Cooper Basin and Pohang sites show a much steeper increase of M max,obs with injected fluid volume, which for Pohang finally led to the occurrence of the largest known event induced by stimulation of a geothermal reservoir (M W 5.5). We note that these two steep trends are not predicted by any of the existing models.  (2014) and Galis et al. (2017) models (Figure 1a). For Basel the occurrence of a M W 3.0 just after terminating the injection (shut-in) is a large deviation from the previous linear trend. Soultz03 displays a steep increase of M max,obs with injected volume at the beginning of the stimulation, which then decreases onto a slope close to unity after 1,000 m 3 of injected fluid.
Interestingly, even injection at the same site does not yield the same evolution of M max,obs , which might be relevant for predicting M max,obs for subsequent injections at the same site. The reservoir in Soultz-sous-Forêts was stimulated during three different periods, in 1993 at 3.5 km depth and in 2000 and 2003 at 4.5-5 km depth. In 1993, seismic activity was limited to very small magnitudes, with M max,obs not exceeding 0.3. Seismicity evolution was mostly linear and is well fitted by the model of McGarr (2014). This stands in contrast to the later injections in 2000 and 2003, resulting in seismic events of much larger magnitudes.
In 2000, within a few hours after stimulation started a M W 1.8 event occurred, which was followed by a linear increase in magnitudes to M W 2.5 at the end of injection. Likewise, the stimulation in 2003 is characterized by first a steep increase of maximum magnitude at the onset of injection followed by a more modest increase for most of the remaining injection period.

Seismic Moment Evolution
Comparing the cumulative seismic moment M 0,cum to the cumulative volume of fluid injected, we find a clear linear relationship (Figure 2a). Most of the seismic sequences display a stable long-term increase of M 0,cum with a slope of about 1, consistent with the model of McGarr (1976McGarr ( , 2014. This suggests that for the dominant part of the injection period there exists a linear response between seismic energy and injected fluid volume. Many sequences are characterized by an initial phase of rapid increase in released seismic moment followed by an almost linear slope. The steep initial increase is especially pronounced in the Soultz03 injection. In striking contrast to most EGS projects analyzed here, Cooper Basin and Pohang show a much steeper increase of M 0,cum with volume injected (Figure 2a). Both slopes are similar to the initial slopes of the other projects yet remain high for the entire stimulation, not switching to the linear slope observed for the other seismic sequences. Interestingly, the Paralana data set shows a jump in M 0,cum during an otherwise linear evolution, caused by the occurrence of several large magnitude earthquakes (2 < M W < 2.5) compared to the remaining seismicity during the same time period. Following the sudden increase, the trend reverts back to the previous linear slope until shut-in. For the Basel sequence the evolution of M 0,cum follows a linear increase with injected volume for the entire injection period, only to deviate shortly after shut-in due to the occurrence of a M W 3.0 earthquake.
We find a linear relationship between M 0,cum and hydraulic energy E hyd comparable to the relation between M 0,cum and injected volume. (Figure 2b) However, several of the projects collapse to a similar trend or are separated by a distinct offset in seismic moment. The slope of the relation of M 0,cum and E hyd is about unity for the most part. Different trends are found for the final stimulation phase of Cooper Basin and the post shut-in phase of Basel. Although punctuated by several large events, the Cooper Basin stimulation follows mostly a similar trend as the other projects. In general, most of E hyd is applied during the period of linear moment increase. Only during the Pohang stimulation M 0,cum increased much more rapidly than for all the other projects.
Note that about the same total amount of E hyd was applied in several of the projects but resulted in very different cumulative seismic moment release. The later injections into the deep part of the Soultz-sous-Forêts reservoir in 2000 and 2003 display very similar relations between M 0,cum and E hyd . In contrast, the 1993 stimulation at shallower depth, with comparable fluid volumes injected and hydraulic energy, exhibits a similar evolution trend yet more than 3 orders of magnitude lower total seismic moment release.

Seismic Injection Efficiency
We analyzed seismic injection efficiency I E for each of the stimulation projects ( Figure 3). Injections performed at Soultz93, KTB, and BGF display very low I E , on the order of 10 −5 . Soultz93 shows a steady

Geophysical Research Letters
BENTZ ET AL.
increase, while KTB exhibits a steady decrease of I E during the active stimulations. Efficiency at BGF first increases and then slowly decreases again. In Pohang and Cooper Basin the injection started with relatively low efficiencies, about 10 −6 . At Pohang the efficiency increased very rapidly during stimulation, interrupted by small periods of decreasing I E coinciding with pauses in the injection activity. The final value reached at Pohang after the occurrence of the M W 5.5 event is larger than 1, indicating that more energy was radiated seismically than hydraulically injected. Cooper Basin exhibits a very irregular growth of I E , increasing during stimulation in two major steps to a final value of 2 · 10 −2 . Basel and Soultz00 show an overall decreasing efficiency during the course of the stimulation, although the initial trend at Soultz00 is dominated by the early M W 1.8 event within hours of injection start. Paralana and Soultz03 first show an increase of I E up to 10 −2 , followed by a slow decrease for the main part of the injection activity. Similarly, for the majority of projects, I E converges toward relatively stable values between 10 −3 and 10 −2 . It is interesting, that a steady decrease of I E at Basel still resulted in the large project stopping M W 3.0 event.

Discussion and Conclusion
Recently, Kwiatek et al. (2019) reported on a successful attempt controlling the evolution of seismic magnitudes during an EGS operation by adjusting the stimulation protocol. But whether or not maximum magnitudes of injection-induced seismicity can indeed be controlled and if events have the same magnitude limit as tectonic earthquakes is currently still a matter of debate (e.g., Hofmann et al., 2018;van der Elst et al., 2016). It is generally assumed that small perturbations of effective stresses are sufficient to initiate failure to critically stressed faults (Zoback & Townend, 2001). The ensuing rupture is expected to propagate until stored elastic energy is consumed at the rupture tip or the rupture hits a sufficiently strong barrier. However, as for tectonic earthquakes, the conditions controlling rupture arrest and final magnitude of induced earthquakes are still poorly understood (Galis et al., 2017;Garagash & Germanovich, 2012;Ripperger et al., 2007). Controlling induced event magnitudes during a stimulation campaign amounts to successfully maintaining event size within the limits of stable rupture propagation (maximum arrestable magnitude of Ripperger et al., 2007;Galis et al., 2017;pressure-constrained ruptures of Norbeck & Horne, 2018). Existing models relating maximum expected event magnitude to injected fluid volume assume an initial phase of stable growth of cumulative seismic moment before runaway ruptures occur (Galis et al., 2017;McGarr, 2014). In this study we examined prominent examples of EGS stimulations with regard to maximum observed moment magnitude, seismic moment evolution, and injection efficiency in response to injection operations.
A striking observation is that the general properties and evolution of seismicity are largely independent of the tectonic stress regime. The analyzed projects cover a range of tectonic settings from dominantly strike-slip and normal faulting to transpressive regimes (Table 1), yet no relationship between M max,obs , seismic moment evolution, or injection efficiency and background stress regime could be observed. Also, we did not observe a substantial variation in b value between different projects where stimulations were performed in different stress regimes (Figure 1b) as has been observed previously (Schorlemmer et al., 2005). This is in agreement with the analysis of 41 case studies of fluid injection projects (Evans et al., 2012) and supported by observations of changes in the local stress state as a function of depth (e.g., in Soultz, Cuenot et al., 2006). In addition, it may be expected that reservoir stresses are considerably modified due to poroelastic stress changes caused by fluid injection (Martínez-Garzón et al., 2013;Schoenball et al., 2014). Instead, we suggest that local tectonic features of the reservoir, such as characteristic fault length, orientation, and frictional properties, govern the seismic response of the reservoir. This may be illustrated by comparing different stimulation phases in Soultz. Here the seismic activity level and M max,obs change substantially between  (Figure 1). During later stimulations, much larger event magnitudes are observed for similar fluid volumes injected, possibly related to a strongly perturbed stress state from previous injections and/or activation of a different structural inventory.
Existing models predicting maximum magnitudes with injected volume fit the data reasonably well. Accounting for site-specific parameters, such as b value, some data may better be approximated by the rupture physics-based model of Galis et al. (2017) or the model of van der Elst et al. (2016), while others follow the relationship of McGarr (2014). However, what is important to note is that for most of the sites the temporal evolution of maximum observed magnitudes with the injected volume generally follows a constant slope, at least for substantial parts of the injection. While the maximum expected magnitude is considered a key parameter for each injection site, we find that continuous tracking of temporal changes of cumulative seismic moment in relation to injection parameters provides insight into site-specific seismicity evolution and potentially provides a first-order seismic hazard information. The data show two distinct populations among the described projects. A larger group displays a stable evolution of cumulative seismic moment across the largest part of injection activity; this includes 8 out of 10 analyzed data sets. We interpret this behavior following Galis et al. (2017) as an indication for stable growth of self-arrested ruptures. In this scenario the size of the earthquakes is likely governed by the induced pore pressure perturbation (Norbeck & Horne, 2018). If the pore pressure, that is, the injected volume or applied hydraulic energy, grows, so does the rupture area. This would result in constantly increasing magnitudes and cumulative moment with injected volume, a behavior we can clearly observe ( Figure 2). Interestingly, the slope characterizing increase in cumulative seismic moment with injected fluid volume and with hydraulic energy input is close to unity for most stimulations. This is in excellent agreement with the prediction of McGarr (2014) suggesting a linear relation between volumetric strain in the reservoir and cumulative seismic moment. Most of the sequences start with a short phase of steeply increasing seismic moment with injected fluid, which stabilizes to a constant slope afterward (Figure 2a). This may be seen as an initial adaption phase of fracture opening and fault activation around the well with a subsequent phase of stable rupture propagation following the expanding pore pressure front.
The second group of projects, including Pohang and Cooper Basin, shows a steep increase of maximum magnitude and cumulative seismic moment, which appears to grow unbound and does not stabilize. These trends are not captured by any of the models and suggest activation of runaway ruptures, where rupture size is only limited by the size of tectonic faults. For Pohang, immediately after onset of injection, seismic moments start to grow seemingly unbound. This is in contrast to the Basel and Cooper Basin project. The substantial part of Basel injection remains stable only to produce a large M W 3.0 event after shut-in. This highlights the problem that exceeding the maximum size of arrestable ruptures may occur at any point during injection, resulting in unstable rupture propagation (Galis et al., 2017). In addition, typically a significant time delay may exist between seismic response and change of injection parameters due to hydraulic diffusion characteristics of the reservoir. The Cooper Basin injection displays time periods of stable growth, punctuated by substantially larger events that may represent unstable ruptures. This explanation of the behavioral change is also indirectly contained in the assumptions of the current models. Following McGarr (2014) and Kostrov (1974), the ratio between average stress drop and the volume affected by injection should scale linearly with cumulative seismic moment. Any observed change in this scaling would indicate a considerable change in either stress drop or stimulated volume, pointing at a fundamental property change of the injected system. Seismic injection efficiencies observed for the different sites vary by 7 orders of magnitude. KTB, Soultz93, and BGF show very low efficiencies, below 10 −4 , indicating a seismically very inefficient process, such as fracture creation (Maxwell et al., 2008) or thermal diffusion in a high-heat environment in the case of BGF (Kwiatek et al., 2014). Most of the other projects converge toward 10 −3 to 10 −2 , values typical for tectonic earthquakes (Kanamori, 2001). Suggesting that the same energy balance holds between tectonic and most injection-induced earthquakes in crystalline basement. The large seismic injection efficiency at Pohang, suggesting more energy radiated seismically than injected hydraulically, most likely reflects that the energy released through the large M W 5.5 event was dominantly elastically stored on an already critically prestressed fault. While we assumed a median radiation efficiency η R of 0.46 (McGarr, 1999), the percentage of energy associated with stress drop that is radiated could potentially vary from site to site. However, McGarr (1999) reported values of radiation efficiency between 0.08 and 0.92, which would lead only to very small changes in the computed injection efficiencies compared to the overall trends ( Figure 3).
In agreement with suggestions of Galis et al. (2017), our observations of a potential transition from stable to unstable growth of seismic moment of injection-induced seismicity has important implications for the design of traffic light systems aiming at controlling the seismic hazard and risk imposed by hydraulic stimulations. The physical model predictions appropriately describe the pressure-controlled injection phase in terms of seismic moment evolution. Any potential trends of seismic moment evolution observed during injection that are completely at odds with the models should be considered seriously, possibly leading to a reassessment of further injection operations. For example, a behavior pointing at unstable rupture growth from the start of injection as found for Pohang, needs to be identified by monitoring not only the exceedance of a critical maximum magnitude but also by near-real-time detection of spatiotemporal evolution of event magnitudes or seismic moment with injected volume. It would be crucial to differentiate between stimulations resulting in runaway ruptures right from the beginning of the injection activities and stimulations that may become stable after an initial adaption phase. So far, the only observed clue to differentiate the two cases is the injected fluid volume up until the steep seismic moment increase persists. In a multistage stimulation project these observations could be utilized by analyzing the seismic moment evolution during each injection stage and modify subsequent injection strategies based on the observed trends. Such a strategy was successfully tested in the recent St1 stimulation project in Helsinki (Kwiatek et al., 2019). Second, the case of a self-arrested rupture becoming potentially unstable with further fluid injection needs to be identified early enough for operators to change injection strategies or cease stimulation. In practice this may prove difficult, as the Basel stimulation showed. Seismic moment and maximum magnitude growth were stable, and injection efficiency was decreasing during main stimulation; nevertheless, the stimulation resulted in the project-stopping M W 3.0 event after shut-in. Thus, trends in the seismic efficiency are not necessarily a reliable test for safe operations. However, the recent success in managing maximum magnitudes in the St1 stimulation in Helsinki, through close interaction of operator and seismic monitoring, highlights the potential of a new generation of traffic light systems combining near-real-time seismic monitoring and subsequently adapted injection strategies. Optimizing such rapid-response traffic light systems may be key to safe development of geothermal reservoirs even in urban areas. This may assist in geothermal energy as a valuable and economic contribution to a mix of renewable energies in implementing the ongoing energy transition.