The increasing number of industrial operations related to hydrocarbon extraction, geothermal power production, hydraulic fracturing for shale gas exploitation, wastewater injection, water impoundment, hydrocarbon storage, and mining operations in recent years, and the potential for large-scale subsurface CO₂ storage in future, has increased the importance of understanding and de-risking induced seismicity both to the scientific community and to the public who live near such operations (Grigoli et al., 2017). The potential to induce seismicity by human activities is well-known (McGarr et al., 2002; Elsworth et al., 2016; Foulger et al., 2018; Keranen and Weingarten, 2018; Schultz et al., 2020). Military waste fluid injected in the Rocky Mountain Arsenal in the 1960's near Denver, Colorado (Healy et al., 1968), induced the so-called “Denver earthquakes.” Since then induced earthquakes related to mining (Arabasz et al., 2005; Fritschen, 2010), oil and gas field depletion (Bardainne et al., 2008; Van Thienen-Visser and Breunese, 2015) shale gas exploitation (Bao and Eaton, 2016; Clarke et al., 2019; Lei et al., 2019), geothermal exploitation (Häring et al., 2008; Deichmann and Giardini, 2009), and waste water disposal (Ellsworth, 2013) have been documented around the world (Baisch et al., 2019). In the UK, induced earthquakes related to hydraulic fracturing at Preese Hall (Clarke et al., 2014), and Preston New Road (Clarke et al., 2019) have been observed, and the latter led to an indefinitely imposed UK government moratorium on fracking.

Traffic light systems (TLS; Bommer et al., 2006; Majer et al., 2012; Mignan et al., 2017; Baisch et al., 2019) are used widely to manage hazard and risk due to induced seismicity in geothermal and hydrocarbon industries, whereby operations are continued (“green”), amended (“amber”), or stopped (“red”) based on the local event magnitude. In the original TLS developed by Bommer et al. (2006), the TLS thresholds are based on peak ground velocity, but other TLSs have been implemented based on earthquake magnitude or other ground motion parameters, such as peak ground acceleration (Ader et al., 2020). Depending on the industrial activities, criteria for a TLS may be very different. Baisch et al. (2019) and He et al. (2020) summarized some examples of existing TLSs that correspond to different industrial activities. In the UK the “amber” and “red” thresholds for induced seismicity related to unconventional oil and gas operations are set to local earthquake magnitudes M_L = 0 and M_L = 0.5, respectively, and this has led to multiple halts of hydraulic fracturing operations during the past few years (Clarke et al., 2019) and finally to an immediate moratorium of operations in November 2019.

The thresholds between zones in TLS are often defined based on limited case studies and on a priori assumptions in a best effort to provide simple schemes (Grigoli et al., 2017; Baisch et al., 2019). Consequently, they do not necessarily take into account the range of possible scenarios, nor uncertainties in event magnitudes, and hence some operations will incorrectly continue, increasing the risk of larger triggered earthquakes, while others will be wrongly halted. To ensure actions taken are robust, it is therefore necessary to estimate local magnitudes with uncertainties, and to consider them in the choice of M_L thresholds in TLSs.

Assessing accurate magnitudes for human-induced earthquakes such as shale gas stimulation, waste water storage, or enhanced geothermal systems is difficult, because they are affected by lack of knowledge about the Earth's subsurface between the source and receivers and by the magnitude scale used (Kendall et al., 2019). A standard approach to determine M_L is to first locate the earthquake and then apply an empirical scaling relation to the source-to-receiver distance (Gutenberg and Richter, 1942; Gutenberg, 2013). Source location-related uncertainties in M_L can then be evaluated using the location confidence ellipses. Unfortunately, estimating errors on M_L due to velocity model uncertainties, energy attenuation during propagation or site effects such as wavefield focusing is difficult.

It is well-known that the accuracy of hypocenter locations depends largely on the velocity model accuracy (Husen and Hardebeck, 2010). Various efforts have been made to estimate velocity model uncertainties, by including a correction term to traveltime curve predictions (Myers et al., 2007), making random perturbations around a given velocity model (Poliannikov et al., 2013) and by locating seismic events in an ensemble of velocity models obtained by a Bayesian analysis of independent data (Gesret et al., 2011; Hauser et al., 2011). Recently, Garcia-Aristizabal et al. (2020) analyzed different sources of uncertainty that can be relevant for the determination of earthquake source locations, and introduced a logic-tree-based ensemble modeling approach for framing the problem in a decision-making context. Their approach, however, is not fully probabilistic, but limited to a finite set of explored models.

Here we propose a way to calculate local magnitudes with uncertainties for microseismic events, and to include the uncertainties in the design of TLS. We use a 3D Monte Carlo non-linear traveltime tomography method to jointly invert for hypocenter locations and velocity model. This allows us to obtain posterior distributions for local magnitude M_L, which cover both velocity and source location uncertainties. Results clearly show that velocity uncertainties and station site effects are significant and change the zones of the TLS to which events are assigned, hence they directly affect safety related decisions. We then apply our method to the hydraulic fracturing induced seismicity at Preston New Road, UK and a mining site, and demonstrate that a red-light would have been encountered earlier if uncertainties would have been accounted for in the TLS thresholds.

Usually, local magnitudes M_L are calculated by first locating the earthquake using standard linearized earthquake location methods (e.g., Klein, 2002), which require simple assumptions about the unknown underlying subsurface seismic velocity structure, and then applying an empirical scaling relation to the source-receiver distance to determine M_L (Gutenberg and Richter, 1942; Gutenberg, 2013). The solution found by such location methods depends on the a priori best guess velocity model, and so it is not guaranteed to find a location near that of the true earthquake. They also cannot represent uncertainties on M_L related to velocity model uncertainties, energy attenuation during propagation, or site effects such as wavefield focusing.

The McMC joint inversion provides us 3D posterior distributions of seismic velocities (V_p and V_s), and of the earthquake hypocenter locations (Figures 4A,B). Therefore, we can calculate hypocentral distance posterior distributions (Figure 4C), which in turn allows us to estimate station-average local magnitudes M_L posterior distributions (Figure 4D) using a scaling relation (e.g., one of Equations 6–8). These distributions include the effects of velocity and source location uncertainties as well as the source radiation pattern on the pdf for event magnitudes. The station-averaged M_L posterior distribution for one source may have a width that spans more than one zone of the traffic light system (e.g., cyan distribution in Figure 4D) which indicates that velocity model uncertainties alone can change the TLS zone to which the earthquake is attributed. Thus, uncertainties affect real operational decisions.

We can now include the velocity and station site effect uncertainties in M_L in a traffic light system (TLS). To do this, we first calculate M_L threshold probability curves using the station-averaged M_L posterior distributions of the microseismic events at Preston New Road. Threshold probability curves describe the probability that an earthquake of a given magnitude is in any one zone of the TLS. They take into account velocity and station-site effect uncertainties, as well as attenuation and geometrical spreading in M_L. Furthermore, the threshold probability curves allow us to draw conclusions about the range of observed M_L for which the probability of any earthquake being in a zone drops below a given confidence level α. The last point is particularly interesting for regulators and operators because it enables them to define the thresholds between zones in such a way that the probability of an earthquake being in each zone of the TLS is always above a chosen confidence level α.

M_L threshold probability curves are obtained by shifting each of the 172 station-averaged event M_L pdfs along the local magnitude axis and estimating the percentage of the distribution lying in each of the three zones of the TLS (Supplementary Figure 3). By averaging over all threshold probabilities curves, we obtain one curve that describes the probability of an earthquake with a given magnitude being in any zone of the TLS. This can then be used to draw conclusions about (1) probabilities of earthquakes of a certain event magnitude M_L being in any one of the TLS zones, and (2) the M_L range for which the probability of any earthquake being in a zone drops below a given confidence level α.

Our approach here is approximate, but once the first inversion if performed, subsequent assignment of a new event's M_L to the correct zone of an adjusted TLS is trivial and can be done in real time, as explained below. In theory, however, the most rigorous approach to incorporating uncertainty into the calculation of M_L for any one event is to retrieve its full posterior M_L distribution, which is the averaged pdf across all stations for that one event. Although at this point we can do this for any existing event in our dataset, in practice we want to be able to do this for each new event that occurs, in real time. This presents a large challenge, however, since formally we must add the travel times from this new event to our dataset and re-run the whole sampling procedure again. We have added one new earthquake to the dataset and sampled 140,000 new models. This took 14,880 CPU-hours on the ARCHER HPC machine, and so remains practically impossible for real time monitoring. Furthermore, the M_L pdf of the new event is still sparsely sampled, and hence does not allow for a robust M_L uncertainty quantification.

We can use the three TLSs (Figure 8) to analyze retrospectively how decisions would have changed at Preston New Road under a TLS⁺ or TLS^– (Figure 9). We compare here the classification in the UK-TLS, a TLS^– and TLS⁺ for a 80% confidence level (Figure 5). That means, the risk of exceeding a TLS threshold in TLS^– is 20%, while in the TLS⁺, the certainty that a threshold was exceeded is 80%. The earthquake on October 19th would have been classified as “amber” in all three TLSs using the maximum probability magnitude, whereas it was classified as “green” by the operator. Hence, action would have been taken earlier and the probability of subsequent larger events would have been reduced. The same is true for the seismicity on October 24th (Figure 9B), where operations would have stopped immediately with a safety prioritizing system (TLS^–), and also in the UK TLS with a M_L that accounts for uncertainties. This demonstrates the importance of accounting for uncertainties in local magnitudes M_L in the decision-making process.

We acknowledge that the occurrence of induced seismicity is a multi-parameter phenomenon, depending on details of subsurface structures as well as on the complete history of operational measures and therefore cannot be predicted. Deformation processes may continue and can still induce seismicity after injections stop. The delay time between hydraulic fracturing completion and the cessation of the observed seismicity can be up to several years (Baisch et al., 2019). It is therefore speculative that an earlier stop would have prevented large magnitude post-injection seismicity at PNR. Nevertheless, it has been shown that lower M_L threshold values in the TLS used for the geothermal stimulation in Basel, Switzerland could have prevented larger magnitude post-injection seismicity (Baisch et al., 2019). We therefore argue that it is critical to establish systems which permit regulators to account for uncertainties while managing risk, as we propose here.

We implemented a fully Bayesian approach for analysing uncertainties, such as velocity model and source location uncertainties in local earthquake magnitudes and evaluate their influence on decision-making for induced seismicity. We conclude that these uncertainties are important, as they can make a difference of up to one or two magnitude unit, and hence directly affect operational decisions by potentially moving an earthquake two zones in a traffic light system (TLS) leading to radically different operational outcomes.

To build a site-specific probabilistic TLS that accounts for uncertainties, the following three steps are necessary: (1) run one fully non-linear hypocenter-velocity tomography for the site and calculate M_L posterior distributions for each earthquake. (2) calculate threshold probability curves, choose a desired confidence level α, and determine the M_L zones A and B below the desired confidence level. (3) attribute zone A and B to “green/amber” or “amber/red” according to the desired safety system (reduce the risk of larger magnitude events (TLS^–) or reduce the risk of halting operations unnecessarily (TLS⁺)). From this point on, real time assignment of any new event's M_L to the correct TLS zone is trival, yet incorporates all the uncertainty in the measurements.

We applied our method to anthropogenic seismicity at a hydraulic fracturing site in the UK, and demonstrate that a red-light threshold would have been encountered earlier in a TLS^–, which possibly could have prevented the UK-wide shut-down. We also applied our methods to mining-related seismicity at Thoresby Colliery, UK and find they apply equally well in this different setting. Hence, our approach can be applied to any site and any form of seismicity.

The original contributions presented in the study are included in the article/Supplementary Materials, further inquiries can be directed to the corresponding author/s.

CR processed the data, performed the analysis, prepared the figures, and wrote the paper. XZ developed the inversion code. BB provided the data of both case studies. All authors contributed to the interpretation of the results and the writing of the paper.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The handling editor declared a shared affiliation with one of the authors AC at time of review.