Induced seismicity closed-form traffic light system for actuarial decision-making during deep fluid injections

The rise in the frequency of anthropogenic earthquakes due to deep fluid injections is posing serious economic, societal, and legal challenges to many geo-energy and waste-disposal projects. Existing tools to assess such problems are still inherently heuristic and mostly based on expert elicitation (so-called clinical judgment). We propose, as a complementary approach, an adaptive traffic light system (ATLS) that is function of a statistical model of induced seismicity. It offers an actuarial judgement of the risk, which is based on a mapping between earthquake magnitude and risk. Using data from six underground reservoir stimulation experiments, mostly from Enhanced Geothermal Systems, we illustrate how such a data-driven adaptive forecasting system could guarantee a risk-based safety target. The proposed model, which includes a linear relationship between seismicity rate and flow rate, as well as a normal diffusion process for post-injection, is first confirmed to be representative of the data. Being integrable, the model yields a closed-form ATLS solution that is both transparent and robust. Although simulations verify that the safety target is consistently ensured when the ATLS is applied, the model from which simulations are generated is validated on a limited dataset, hence still requiring further tests in additional fluid injection environments.


Results
Predictive hazard model. A predictive model lies at the heart of any risk assessment. In the case of induced seismicity, a wide range of statistical and physics-based models exists 16,17 . Undoubtedly, more work is needed to develop, calibrate and validate new models; however, we believe that the missing link is the use of such models for deriving and monitoring quantitative risk thresholds. Here, we use as example a simple and yet robust model that forecasts the piecewise induced seismicity temporal, here daily, rate λ(t, m ≥ m 0 ; θ) as: where V  (t) is the injection flow rate as a function of time t in m 3 /day, θ = [b, a fb , τ] a set of model parameters describing the underground characteristics (earthquake size ratio, activation feedback in m −3 and mean relaxation time in days, respectively), m 0 the minimum magnitude cutoff, and t shut-in the shut-in time, also in days. Both a fb and b can also be functions of time t (see Decision variable section). In this model, the injection or operation phase is described by a linear relationship between λ(t, m ≥ m 0 ) and  V (t) in line with previous studies [17][18][19] . It derives directly from the linear relationship between  V and overpressure 17 , hence assuming no change of injectivity during any given stimulation. The post-injection phase is described by a pure exponential decay representative of a normal diffusion process 17 . Although the Modified Omori Law is sometimes used to describe post-injection seismicity 20 , reasons remain mostly historical 21,22 . The proposed alternative is verified to be consistent with the tested data (see Methods section) and preferred on analytical grounds, being directly integrable in contrast with the Modified Omori Law, which is conditional on parameter values and may require the definition of an ad hoc upper bound 22 . Finally, no maximum magnitude M max is imposed. It follows from Eq. (1) that induced seismicity is characterized by both the injection profile  V t ( ), and the underground feedback described by the three-parameter set θ. The main limitations of the proposed model are discussed in detail later on.
In contrast with complex fluid modelling 16 , the closed-form Eq. (1) can be computed on-the-fly; moreover, it includes the mean relaxation time, τ, hence taking into account the long-term underground feedback after shut-in. Finally, being integrable, it leads in turn to a closed-form ATLS, as demonstrated in the next section. Although the physical process governing the rate of induced seismicity is more complex than what is represented by Eq. (1), this rate model is proven to be valid in a Poissonian probabilistic setting (see Methods section). Moreover, the physical processes are either not still clear (in fact, there is not an unanimous consensus among scientists), or computationally expensive. Therefore, pragmatism imposes the use of statistical models until both an agreement is found on the physics of induced seismicity and computational time of complex physical modelling is reduced.
The model (Eq. 1) was fitted to six induced seismicity sequences observed in fluid injection experiments from enhanced geothermal systems (EGS) 13,[23][24][25][26] , the initial stage of a long-term brine sequestration 27 , and one fracking at an oil field 28 (Table 1). The model succeeds to describe most of the data as shown in Figure 1 (see results of statistical tests in the Methods section -Note that the rare outliers above 3σ may be due to missing on-site data that may affect seismicity, such as unknown technical operations on wells, or to second-order physical processes missing in Eq. (1), as discussed below). The parameters are found to range over 0.8 ≤ b ≤ 1.6, −2.8 ≤ a fb ≤ 0.1 m −3 and 0.2 ≤ τ ≤ 20 days and show a relatively large scattering between sites and between different stimulations at a same site (Table 2). It is important to note that we specifically chose those six datasets, as they had been made publicly available. Those cases are characterized by high pressures and flow rates, and are rich in induced earthquakes. This represents a selection bias and the a fb parameter could decrease to much lower values elsewhere (Fig. 2). For instance, most injection wells in the U.S. do not cause felt earthquakes 29 ; large a fb variations between regions and sites might be explained by different regional crustal stresses 30 . The present sites provide however natural laboratories to test our model and the associated ATLS, without generalizing or inferring any high level of risk for all existing deep fluid injections. Figure 2 illustrates the parameters' scattering, including results from past studies 18 for EGS and other hard rock settings enlarging the range of values to 0.7 ≤ b ≤ 2.2 and −4.2 ≤ a fb ≤ 0.4 m −3 . It has been shown that in fracking environments, the activation feedback can be as low as a fb = −9.25 m −3 18 . Safety criterion. A safety criterion is a probability of exceedance that can be fixed with respect to different safety metrics, such as fatalities, economic loss, building damage or level of nuisance 14 . Given the selected metric, the corresponding safety criterion can be converted in the magnitude space into the probability of exceedance Pr(m ≥ m saf ) = Y, which will ensure that the acceptable level of risk is preserved at all time: with m saf the magnitude at which the given safety limit (e.g., damage, fatality) is reached. Note that Eq. (2) is a closed-form expression where V(t) is the cumulative injected fluid volume, and m saf and Y are derived from the safety criterion (see Methods section). Moreover, the set of parameters θ is updatable at any given time. The mapping from risk to earthquake magnitude is required to control injection operations based on short-term observations (see Decision Variable section). While peak ground velocity (PGV) is a more direct measure 12 , conversion to magnitude is in any case inevitable to estimate the risk potential of larger earthquakes from the b-value. It is also a unique measure, while PGV requires a location that is not trivial to assign. It should be added that no maximum magnitude M max is imposed in Eqs (1-2). This implicitly assumes that both small to medium-size induced events and large triggered earthquakes on existing faults are treated the same way. This remains debated 31 although a recent study 19 demonstrated that the observed M max in fluid injections is compatible with the null-hypothesis of the Gutenberg-Richter law with no upper limit. The role of M max (and therefore of triggered earthquakes) is however only critical when the risk of fatalities (e.g., individual risk IR) is evaluated. For nuisance or minor damage thresholds, risk is more likely dominated by medium-size induced events. This important discussion has no significant impact on the method proposed, as proved in the Methods section. Before an ATLS is set, the likelihood of failure of the planned project with respect to a specified limit state function defined by the safety criterion in magnitude space (Eq. 2) can approximately be determined. In this study, we select as main metric the annual individual risk (IR) over the entire project period, and as safety criterion IR ≤ 10 −6 (i.e. the probability that a statistically representative individual dies for the introduced hazard), which is a threshold commonly enforced for hazardous installations 14 . Figure 2 shows the acceptable domain for a fixed limit state function, tested for different injection scenarios in which a hypothetical project plan is to inject a total volume V = 10,000 m 3 of fluids at a depth of 4 km 32 . The injection profile is assumed to be flat with a constant flow rate V  = 1 or 10 m 3 /min, having an impact on injection duration and tail behaviour (Eq. 2; Fig. 2). The project is considered to be located at a distance d = 0 km or 50 km from the nearest building. For a given site, with no knowledge of the underground feedback to fluid injection, project operators and regulators in an EGS setting (where high seismicity rates are common during stimulation) could use the known θ scattering for an a priori parameterization. This preliminary assessment shows the likelihood of the project to pass or fail the safety threshold. As shown in Figure 2, the results can be ambiguous, due to the large uncertainties associated with subsurface characteristics. Nevertheless, it provides a preliminary assessment of the risk reflecting the limited knowledge of the induced seismicity process. Future estimations of θ at additional injection sites will likely refine the results and improve the decision process. In addition, rules of decision-making under uncertainty can account for that ambiguity 9,33 . Decisions become more obvious in cases in which the diagram would be entirely green (clear go) or red (clear no-go). Underground stimulation activities in areas with low exposure (e.g. remote EGS plant locations with large distance d from the nearest habitations) would evidently have a lower induced-seismic risk and, thus, shrink the red area. The termination of the 2006 Basel EGS project was due to the high induced-seismic risk emerged from the high exposure of the urban built environment 7,9 . Note that Eq. (2) can be used to predetermine a distance d for which the induced-seismic risk would become acceptable-conditional to a given injection profile V(t) and parameter set θ (since m saf is a function of d; see Methods section).
Decision variable. The ATLS decision variable must be selected and updated as new data allows estimating θ more accurately, or if the planned injection scheme is changed. Here, a threshold earthquake magnitude m th is used as decision variable. In particular, m th is defined as the magnitude value for which mitigating actions must be taken, here corresponding to stopping injection, i.e.
th a bm shut in saf 10 fb saf (see Methods section). If m th is updated "on-the-fly", the project is guaranteed to meet the defined safety criterion.
To avoid reaching m th before the planned stop of the fluid injection one may be inclined to reduce the flow rate V  , but here it would only delay the time at which the injection must stop, as the risk is mostly controlled by total volume injected V (Eq. 2; the secondary role of  V is highlighted in Figure 2 by the change of width of the yellow band, for different  V and τ). It is common practice to reduce V  however 13 , but results of such action remain unclear 34 since verification of a safety threshold requests a large number of experiments (see simulations below); indeed: (i) It is plausible that such action has no overall effect, the risk remaining the same in average over a fixed V; (ii) If such action has an effect, it would indicate that the model of Eq. (1) does not properly describe the role of different injection strategies. Despite the proposed model being classified as a pure statistical method, it is based on physical considerations. Its first term, for example, builds on the linear relationship between volume change and pressure 17 . Other relationships can be envisioned such as a bilinear relationship indicative of a change of injectivity 26 , which may explain some second-order relationship observed between flow rate and M max 35 . The linear relationship could also be shifted in time by including a minimum pressure threshold 17 below which no induced seismicity is triggered. Adding such processes would likely allow for smarter mitigation strategies in which the shape of the injection profile would play a role. Since any model change would require the inclusion of additional parameters (which have yet to be constrained), and since Eq. (1) is verified to be consistent with most of the data tested, we consider Eq. (1) to be a reasonable first-order model for the proposed ATLS.
To validate the ATLS in a realistic time-dependent setting, we simulate the Basel induced seismicity sequence using Eq. (1) with  V Basel (t) (Figs 1; 3a) and θ Basel (t) = [b(t), a fb (t), τ = 1.12 day] (Figs 2; 3b; Table 2). We use the safety threshold IR ≤ 10 −6 with the nearest building above the borehole (d = 0 km), yielding m saf = 5.8 (see Methods section). The value m th as a function of time is shown in Figure 3c with two examples of induced seismicity time series. The first one, in grey, is a reproduction of the 2006 Basel experiment; the second one, coloured, is the case where the ATLS is used. Figure 3d finally shows that the safety target, while not respected for synthetic versions of the Basel experiment, is reached once the newly proposed ATLS (Eq. 3) is considered. Due to the stochastic nature of the earthquake process, the operational safety target is only reached on average over multiple sequences. It is plausible that an improved model with a non-linear relationship between flow rate and overpressure would open the possibility for mitigation strategies based on different injection profiles, hence potentially avoiding prematurely stopping the stimulation and the project.

Discussion
The main purpose of this study was to present a statistical-based first order ATLS, which verifies that a quantitative safety criterion is ensured. An important benefit of the outlined ATLS approach to both operators and regulators is its transparency and execution speed (being a suite of simple closed-form expressions). By principle, its adoption would make any project in compliance with the safety threshold whatever the response of the underground. This does not ensure that a project is financially successful, but it gives the operator the maximum allowable chance to reach success, based on a quantitative risk-based method.
It is important to note that existing traffic light systems based on heuristics can already provide reasonable results (for instance, the magnitude threshold m th = 2.9 fixed during the Basel experiment 13 is not far off the m th (t) computed with the new ATLS; Fig. 3c). However, magnitude thresholds enforced by different jurisdictions vary significantly (with 0.5 ≤ m th ≤ 4 10 ) with no clear link with the standard risk-based safety criteria used in other hazardous industries 14 . Although the application of the ATLS in the decision process might appear at first complex, there is substantial evidence that the algorithmic (or actuarial) approach is superior to the so-called clinical approach of informal judgement 36,37 . This must apply too to induced seismicity prognostics, experts necessarily basing their judgements, consciously or intuitively, on past observations shown here to be reasonably well described by Eq. (1) (Fig. 1 and Methods section).
Finally, we suggest the following future directions: (1) improve the model (Eq. 1) by relating directly overpressure instead of flow rate to induced seismicity. Due to potential changes in injectivity, gas kicks and other processes, overpressure is likely to provide a better proxy than injected fluid volume; (2) test the statistical model in other fluid injection environments; (3) improve the updating of the parameter space by using a hierarchical Bayesian framework where the uncertainties of the model parameters are taken into account; and (4) modify the mapping from risk to magnitude space for site-specific conditions which likely vary between fluid injection locations. The present study demonstrated the power of the actuarial approach and should be considered as a proof-of-concept for future physics-based induced seismicity models, more sophisticated engineer-based risk assessments, and improved mitigation strategies.

Methods
Time series analysis. The induced seismicity temporal rate model λ(t) of Eq. (1) is fitted to the induced seismicity data by using the maximum likelihood estimation (MLE) method 38 . The probability p i that the observed number n i of induced earthquakes results from a Poisson process 39   (1). The maximum likelihood estimate of θ is finally θ MLE = arg max θ LL(θ, X). m 0 is fixed to M c , the completeness magnitude defined as the magnitude bin with the highest number of events 40 . b is estimated independently of λ, also based on the MLE method 41 . The model is fitted to 8 datasets (from 6 stimulations in various injection settings; Table 1); the resulting maximum likelihood estimates are listed in Table 2.
Sensitivity analysis. Temporal changes in a fb and b are evaluated for induced seismicity sequences that are large enough (i.e., made of hundreds of events, such as in 1994 Paradox Valley, 2006 Basel and 2014 Newberry). The parameters are estimated using a moving window with constant event number n = 100. Before n is reached, MLE estimates obtained in retrospect (Table 2) are used, as shown in Figure 3b. In a prospective case, generic values should be used, e.g. the median or mean of values taken by the parameters in previous experiments. Since the post-injection phase is not considered in the sensitivity analysis, a fb is directly obtained from Eq. (1) so that Post-injection data analysis. While the linear relationship between λ and V  is well established [17][18][19] , the pure exponential behaviour of the induced seismicity post-injection tail has only been demonstrated for the 2006 Basel case 17 . Here, we compare three relaxation models: pure power law t t by using the Akaike Information Criterion 44,45 . We find that, of the 7 datasets (discarding KTB94a, which tail is likely cut), 5 are best described by the pure exponential function. The 2 others are best described by a stretched exponential function with stretching parameter β = 0.9 and 0.7, for PV94 and NB14b, respectively (β = 1 representing the pure exponential). This justifies the use of a pure exponential in Eq. (1) hence limiting θ to a simple three-parameter set.
Model goodness-of-fit. Figure 1 shows the general agreement between model and data by visual inspection. Additionally we test Eq. (1) against our datasets using the Kolmogorov-Smirnov (KS) confidence bounds 46 .
We first convert the dataset  t ( ) into a transformed dataset ( ) ∼ Τ  as follows:  Figure 4 with the dashed lines representing the two-sided 95% and 99% confidence intervals. The following conclusions can be drawn: The model performs well for the datasets KTB94b, B06 and G11; fairly well for the datasets KTB94a, PV94, NB12 and NB14b; and poorly for dataset NB14a. Given the range of different datasets, we conclude that the rate model expressed in Eq. (1) describes fairly well the relationship between injection flow rate and fluid-induced seismicity rate.
Time series simulation. We simulate induced seismicity time series with time-dependent flow rates V  (t) and θ(t) = [a fb (t), b(t), τ] by using the thinning method for the injection phase 47 . For the post-injection phase for which τ is assumed constant, as no control over post-injection seismicity is possible, the inversion method is used instead to simulate event occurrence times t 48 . Magnitudes m are also simulated using the inversion method (for both injection and post-injection phases with b time-dependent in the first phase).
Risk-to-magnitude mapping. We translate the safety threshold IR = 10 −6 as Pr(I = 9) = 10 −5 assuming that a statistically average "poor" building collapse is reached for a ground intensity I = 9 (i.e., most buildings of class A in the EMS98 code suffer collapse) 49 and that once a building collapses, there is a 10% chance of individual fatality 50   The correction is based on the observation that induced earthquakes would seem less severe in average than tectonic ones. We assume that the Modified Mercalli Intensity (MMI) 49 and the USGS Community Internet Intensity (CII) 51,52 are equivalent for sake of simplicity. For a specific project, a fully probabilistic risk approach 9 is recommended to derive the parameters Y and m saf of the ATLS closed-form expressions (Eqs 2-3). Here, the safety thresholds shown in Figure 2 use Y = 10 −5 with m saf (d = 0 km) = 5.8 and m saf (d = 50 km) = 7.9 (equivalent in the tectonic case to 5.0 and 7.1, respectively). This simplified approach is used to illustrate our ATLS proof-of-concept for a general case with no site-specific conditions. A detailed risk analysis, which would integrate risk over the full magnitude range 9 , remains outside the scope of the present study. The high values of m saf are due to using IR as the safety metric. Some experts may disagree that such high magnitudes can be reached in the induced seismicity context 31 , although it is statistically plausible 19 . The choice of the safety metric has however no impact on the validity of the proposed method. For example, using the minor (1). Using the Kolmogorov-Smirnov (KS) confidence bounds, we find that the model performs well for the datasets KTB94b, B06 and G11; fairly well for the datasets KTB94a, PV94, NB12 and NB14b; and poorly for dataset NB14a (see the Methods section for details).
SCIEnTIFIC RePoRTs | 7: 13607 | DOI:10.1038/s41598-017-13585-9 damage threshold Pr(I = 6) = Y 2 (i.e., most buildings of class A in the EMS98 code suffer negligible to slight damage) would yield m saf2 (d = 0 km) = 4.0. From Eq. (3), the same ATLS m th would be computed for both Pr(I = 9) = 10 −5 and Pr(I = 6) ≈ 10 −3.1 assuming e.g. b = 1 in log 10 (Y 2 ) = log 10 (2) with an impact on the results only if Y was tending to 1, which is unlikely in safety norms). We then define the ATLS as the operational magnitude threshold m th at which the injection is stopped in order to meet the safety target. Note that m th also provides the completeness magnitude to attain in the region and thus a basis for seismic network monitoring planning 42,53 . We thus obtain the following system of equations (for Y << 1): The second equation is always true since the expected number of events with m ≥ m th is one, given the assumption that injection is stopped for t = t shut-in as soon as m ≥ m th is first observed. Substituting V(t shut-in ) in the first equation of this system yields Eq. (3). It is worth noting that omitting the post-injection tail effect with τ = 0 yields the basic frequency-magnitude distribution threshold m th = log 10 (Y)/b + m saf .
Data availability. All the data used in this study are publicly available. For more information, please contact arnaud.mignan@sed.ethz.ch.