Statistical bounds on how induced seismicity stops.

Earthquakes caused by human activities receive scrutiny due to the risks and hazards they pose. Seismicity that occurs after the causative anthropogenic operation stops has been particularly problematic-both because of high-profile cases of damage caused by this trailing seismicity and due to the loss of control for risk management. With this motivation, we undertake a statistical examination of how induced seismicity stops. We borrow the concept of Båth's law from tectonic aftershock sequences. Båth's law anticipates the difference between magnitudes in two subsets of seismicity as dependent on their population count ratio. We test this concept for its applicability to induced seismicity, including ~ 80 cases of earthquakes caused by hydraulic fracturing, enhanced geothermal systems, and other fluid-injections with clear operational end points. We find that induced seismicity obeys Båth's law: both in terms of the magnitude-count-ratio relationship and the power law distribution of residuals. Furthermore, the distribution of count ratios is skewed and heavy-tailed, with most earthquakes occurring during stimulation/injection. We discuss potential models to improve the characterization of these count ratios and propose a Seismogenic Fault Injection Test to measure their parameters in situ. We conclude that Båth's law quantifies the occurrence of earthquake magnitudes trailing anthropogenic operations.

Introduction

The injection of fluids in the subsurface has the potential to cause faults to reactivate with earthquake slip[1,2]. Numerous high-magnitude cases of induced seismicity causing significant impacts[3-6] have led to increased attention. Study topics have covered a broad range of aspects including triggering mechanisms, maximum magnitude, potential risks/hazards, and managements strategies[7,8]. Despite the proliferation of scientific studies, relatively little attention has been devoted to understanding how induced earthquake sequences stop. This is surprising because earthquakes trailing well shut-in have been identified as one of the most significant sources of risk for induced seismicity management[9]. In practice, there are many cases where the largest magnitude earthquake trails the shut-in of the responsible well[5,10-12] or where abrupt ‘jumps’ in magnitude have occurred[13].

In this study we build a conceptual framework for statistical bounds on how induced earthquake sequences stop. We start by following the development of Båth’s law for aftershocks of natural earthquakes and adapt it through a set of assumptions and justifications for induced earthquakes. We test the validity of this framework on sequences of induced seismicity and find that in most cases the majority of earthquakes occur during the stimulation/injection periods, rather than trailing well shut-in. As well, the largest events occurring during the stimulation/trailing periods follow Båth’s law. We also discuss the implications of our results and suggest recommendations for in situ tests that could be useful for constraining induced earthquake risks. Finally, we speculate on how improved physical models could be used to build upon our empirical analysis of how induced earthquakes stop.

Background information on tectonic aftershocks and Båth’s law

Båth’s law is an empirical observation that the average magnitude difference between a mainshock earthquake and its largest aftershock is approximately 1.2 magnitude units[14]. Statistical models of earthquakes[15,16] have come to interpret Båth’s law as a probabilistic consequence of two more fundamental observations[17-20]. The first being an aftershock law (i.e., Omori’s law), which describes the temporal decay of earthquake rate following a main shock[21-23].In Omori’s law is the time-dependent rate of aftershocks, K describes the productivity of the aftershock sequence and strongly influences the total number of aftershocks, p is the decay exponent, and c is an ad hoc parameter added to ensure the equation doesn’t diverge at the time of the mainshock (i.e., t = 0). In this formulation, the value of p is bounded such that for a finite number of earthquakes to occur[24]. The second fundamental observation is the Gutenberg Richter magnitude frequency distribution (GR-MFD), which describes the relative number of small and large earthquakes in a sample[25,26].In the GR-MFD N is the sample size of the catalogue, a is the a-value which is related to the productivity of the sequence, and b is the b-value, which controls the proportionality of big-to-small magnitude earthquakes.

Omori’s law governs the total number of earthquakes that follow a mainshock N and the GR-MFD controls the distribution of earthquake magnitudes ‘randomly’ sampled in that population. The combination of these two observations gives the expression for the average magnitude difference (where N > N > 1).where the total number of earthquakes in a catalogue N is an observable and the number of aftershocks N can be forecast/estimated from the Omori’s law parameters. Typical estimates of Omori and GR-MFD parameters yields estimates of the average magnitude difference that are similar to the empirical observation of ~ 1.2 M[27,28].

Deriving a trailing earthquake model for induced seismicity

Following a similar rationale, we derive a trailing earthquake model from the statistical viewpoint of induced earthquakes as a point-process. We first assume that induced earthquake magnitudes are independent and identically distributed random variables that follow the GR-MFD. This assumption is supported by numerous cases that obeying this distribution[2,29,30], similar to naturally occurring tectonic sequences. From this assumption an expression for the distribution of the maximum magnitude earthquake follows (in the limit of large N).This is often referred to as the sample size hypothesis[29] for the distribution of maximum magnitudes Mmax, where the newly introduced parameters are the magnitude of completeness MC[31] and a confidence level that can be utilized to sample the Mmax distribution. The distribution of expected maximum magnitudes for different population sizes is shown in Fig. 1a.

Figure 1

Concepts of Båth’s law and a simple adaptation for trailing seismicity. (a) Probability density functions (solid lines) and cumulative density functions (dashed lines) for the expected maximum magnitude for various population sizes N, based on numerically sampling Eq. (4). (b) Probability density functions (solid lines) and cumulative density functions (dashed lines) for the magnitude difference ΔM for various population ratios R, based on numerically sampling Eq. (5). (c) Simple model of induced seismicity, where earthquakes lag injection/shut-in by time ΔT and trailing seismicity decays with time.

Based on this sample size hypothesis, we can consider an archetype of an induced earthquake sequence that is broken into two segments: the injection/stimulation phase (denoted by a subscript S) and the trailing shut-in phase (denoted by a subscript T) that follows (Fig. 1c). The stimulation phase would entail a population count of N stimulation earthquakes that occur until the anthropogenic operation is halted at time T, lagged by stimulation events for time ΔT, and followed by aftershock-like events (starting at an earthquake rate approximately equal to the stimulation rate). We consider trailing seismicity to be those N events that occur after time T i.e., the sum of lagged stimulation events and the aftershock-like events. Taking the difference of the maximum magnitudes (Eq. 4) for both segments and the stimulation segment yields the relationship for the trailing magnitude difference below.This produces a form for the differential maximum magnitude that is reminiscent of Båth’s law (Eq. 3), where 1/R = N/N = 1 + N/N = 1 + R. Here we note that ΔM and are physically observable quantities that may be subjected to hypothesis testing.

In continuing to develop our model, we assume a further constraint on the stimulated and trailing earthquake counts N and N. Similar to the prior logic in Båth’s law, we bring additional quantification to the value of . First, the stimulated earthquake count N has been suggested to have a linear relationship with injection volume[32].This expression effectively modifies the a-value of the GR-MFD to have an explicit dependence on the volume injected V(t) (in m3) and the background potential for earthquakes i.e., the seismogenic index Σ. This relationship assumes a non-homogenous Poisson process for the stimulation events[33] and has been well fit to numerous induced earthquake sequences[29,34,35]. Second, the total aftershock count is the integration of an aftershock rate decay function from the time of the mainshock to infinite time. In practice, induced sequences have often been observed to fit aftershock-like rate decay sequences following shut-in[36]. While there is evidence to suggest a linear dependence between injected volume and cumulative earthquake counts[29], we note that alternative models[7,37] could be utilized without changing the results of our study. Because there is an abundance in the choice of aftershock rate decay model, we summarize the prominent ones in Table 1. For each of the aftershock rate decay models we provide the rate decay function , the cumulative count function , and the total aftershock count . We also refer the reader to prior summaries of tectonic aftershock models[38]. In this study we included the most general version of aftershock models that (with parameter constraints) lead to finite cumulative event counts and are widely adopted by the seismological community.

Table 1

Aftershock models proposed for tectonic earthquakes.

Model name [References]	Aftershock rate decay	Cumulative aftershock count	Total aftershock count
Modified Omori[39]	n^ot=Kt+c^-p	N^ot=Kp-1c^1-p-t+c^1-p	N^o∞=Kc^1-pp-1
Exponential[40]	n^et=n^e0e^-t/τ	N^et=τn^e01-e^-t/τ	N^e∞=τn^e0
Stretched exponential[41]	n^st=qN^⋆e^dt_0^q1t+dt+dt_0^qe^-t+dt_0^q	N^st=N^⋆1-e^dt_0^qe^-t+dt_0^q	N^s∞=N^⋆
Cut-off Power Law[42]	n^ct=Kt+c^-pe^-t/τ	N^ct=Kτ^1-pe^c/τΓ1-p,cτ-Γ1-p,t+cτ	N^c∞=Kτ^1-pe^c/τΓ1-p,cτ
Gamma[43]	n^γt=At^-qγq,λ_bt-γq,λ_at	N^γt=Aq-1λbq-1-λaq-1--λaq-1e^-λ_at+t^1-pΓq,λ_at+λbq-1e^-λ_bt-t^1-pΓq,λ_bt	N^γ∞=Aq-1λbq-1-λaq-1
Epidemic type[44]	n^Et=∑_iK_it-t_i+c^-p	N^Et=∑_iK_ip-1c^1-p-t-t_i+c^1-p	N^E∞=∑_iK_ic^1-pp-1

For each model we provide an expression for the aftershock rate function , the cumulative aftershock count function , and the total number of aftershocks . Corresponding references are also provided for each model. Note that and are the lower and upper incomplete gamma functions, respectively. Superscripts are used to denote the aftershock decay model used.

We can expand upon Eq. (5) by imposing our archetypal induced earthquake sequence (Fig. 1c) in combination with the aftershock and stimulation models. Combining these models with the archetypal sequence allows for closed form expression of , and thus the distribution of the expected maximum earthquake magnitude.Here, has been described for both the Modified Omori (, Eq. 7) and Exponential (, Eq. 8) aftershock rate decay models. Note that we’ve assumed the initial aftershock rate is some fraction f of the ending stimulation rate—to account for slight increases or decreases in earthquake rate after the lagged shut-in response.

Testing the model’s validity

In adapting Båth’s law to induced seismicity, we’ve made several assumptions (Fig. 1): that the sequence is described by a nonhomogenous Poisson process, that the extended fault system can effectively be modelled as a point-process, that there’s a single lag time between operations and seismogenic response, that statistical parameters are stationary during the observation interval, and that there are no injection-volume or tectonic limits on the maximum magnitude. To validate (or refute) these assumptions, we begin by examining induced seismicity cases.

We focus only on cases where there has been a clearly delineated end of operations (Fig. S1). This criterion dominantly encompasses cases related to hydraulic fracturing[2] or enhanced geothermal systems[45]. We select cases of induced seismicity where both catalogue and injection time series information is available. Cases include the Duvernay[11,13,46,47] and Montney[48] in Canada, the Bowland Shale in the UK[49,50], the Appalachian Basin[51,52], SCOOP & STACK play[53-55], Eagle Ford[56], and Arkoma Basin[57] in the United States, the Sichuan Basin in China[6,58,59], Soultz-sous-Forêts in France[60,61], Cooper Basin[62] and Paralana[63] in Australia, KTB[64] and Groß Schönebeck[65] in Germany, Basel[66] and St Gallen[67] in Switzerland, Helsinki in Finland[68], Pohang in South Korea[5,69], and Berlín in El Salvador[70]. Many of these cases have data available in online repositories[71], within aggregate comparison studies[12,29,72], or as supplements to the individual case studies.

Briefly, individual cases are truncated at their magnitude of completeness (MC) and then subdivided into two populations: those occurring during stimulation or those trailing stimulation. We report statistics such as earthquake counts, largest magnitudes, b-value, and MC, which we summarize in Table S1. For additional details on the data acquisition and preprocessing, we direct the reader to the Supplementary Information (Section S1).

Examining the Statistics of RS

From our compiled database (Table S1), we examine the parameter observed during numerous induced seismicity sequences (Fig. 2a). shows a highly skewed distribution—tending towards values of 77% on average and 86% as a median (Fig. S5). This is comparable to observations in prior studies on hydraulic fracturing that found ~ 90% of events occurred during stage stimulation in the Duvernay[35] and ~ 87% of events in the Sichuan Basin were ‘forced’ under a modified Epidemic Type Aftershock Model[6,58]. While the distribution is skewed towards high values, we also note that it is heavy-tailed with a handful of cases with low values. For example, the enhanced geothermal stimulation at Pohang shows a large trailing aftershock sequence[5,69].

Figure 2

The empirical distribution of and bootstrapped beta distribution fit. (a) Histogram of for all the induced seismicity cases (blue bars) with its median and mean values (vertical lines), best fit to the beta distribution (crimson lines), compared against a random distribution (grey bars). (b) Scatterplot of bootstrapped beta distribution parameters (blue circles) with the best fit parameter (red circle). Inset histograms show the univariate distribution of beta parameters.

To further examine the skewed nature of the distribution, we compare our empirical observations to the distribution expected from a random process (Fig. 2a). To model the random process, we consider a uniform random variable for both N and N. We then perform the two-sample Kolmogorov–Smirnov test between our empirically observed and the random process model. The two-sample Kolmogorov–Smirnov test is a non-parametric hypothesis test that assesses the probability that two samples are drawn from the same (arbitrary) distribution[73]. In performing this test, we find extremely small p values that are orders of magnitude less than 5% (Fig. S6)—which rejects the null hypothesis that the two samples are from the same distribution. The failure of to be modelled as a random process suggests that there is a physical process underlying this distribution (e.g., Eqs. 7 and 8).

A better representation of the statistics of the would be useful for estimation of differential magnitudes ΔM (e.g., Eq. 5). Towards this end, we fit our histogram results to a beta distribution—a two-parameter (α and β) family of distributions defined over the interval 0–1[74]. We select the beta distribution due to its flexibility in fitting probability density functions over this interval. Further, to capture the variability in the fitted parameters we performed a 1000-trial weighted bootstrap analysis (see Supplementary Information Section S3 for more details). We find that the mean best-fit parameters for α and β are 1.016 and 0.271, respectively. Similarly, the median best-fit parameters for α and β are 0.940 and 0.238, respectively. The covariances (σαα, σββ, σαβ) for these parameter fits are 0.1151, 0.0130, and 0.0317 (Fig. 2b). In performing sensitivity tests, we do not find that fitted parameters differ appreciably, when comparing hydraulic fracturing or enhanced geothermal system subsets. However, larger datasets could be useful to discern systematic differences in according to operation and setting. These fitted distribution results could constrain trailing seismicity differential magnitudes ΔM, provided the validity of using in Båth’s law for induced seismicity was established.

Testing Båth’s law for induced seismicity

With constraints on the empirical distribution of , we examine the statistical validity of Båth’s law for induced seismicity. Equation (5) predicts that there is a linear relationship between the trailing-stimulation maximum magnitude difference ΔM (left hand side of equation) and the b-value scaled logarithm of the population ratio (first term of the right-hand side), where the remaining term describes the statistical distribution of the residuals (second term of right-hand side).

Based on this, we first consider a weighted linear regression applied to our dataset (Fig. 3a). We observe a reasonable agreement between observation and expectation. The fitted parameters for the linear regression are 1.06 ± 0.16 for the slope and 0.58 ± 0.18 for the y-intercept. The theoretically expected slope of 1.0 agrees within error of our regression-fitted slope (see Supplementary Information Section S4 for more details). The observation of a y-intercept term is the result of a correction term that conditions on the sample count, when there is an unknown true sample count (see van der Elst et al.[29] Sections 2.2 & A4 for additional details). The value of this correction term approaches (i.e., ~ 0.43) as the population count N becomes large, which agrees within error of our regression fitted parameter. While both the empirically derived slope and y-intercept agree within error to their expected values, we note that the best estimate of both values are higher than those expected; this larger value bias could be the result of a publication bias where scientific studies preferentially describe cases with ‘exceptional behavior’, like large magnitudes, red-light events, or lingering sequences of trailing seismicity.

Figure 3

Testing Båth’s law for induced seismicity. (a) Linear regression of magnitude differences ΔM versus logarithmic population ratios : data is separated by anthropogenic type (circles) and shown alongside the 1:1 line (solid gray line), the relationship expected by Båth’s law (solid black line), the regression fit to data (dashed line), and the corresponding 5/95 percentiles (dotted line). (b) The residuals to the expected Båth’s law relationship are shown as probability density function for the expected (solid black line) and empirically derived from data (blue bars). Inset panel shows the cumulative density function for the expected (solid black line) and empirically derived residuals (blue circles/line).

The second test we consider is an examination of the distribution of the Båth’s law residuals (Fig. 3b). In this test, we rearrange the terms in Eq. (5) such that just the power law distribution remains: i.e., . Note that individual correction terms are accounted for using the true population size, rather than the large N limit. We find a correspondence between the histogram of data residuals and the expected probability density function for Båth’s law. More quantitatively, we compare the expected and observed residual samples using the two-sample Kolmogorov–Smirnov test. We find a p value (~ 45%) that is larger than 5%, meaning we are unable to reject the null hypothesis that the two samples are from the same distribution. Most bootstrapped perturbations (~ 61%) support this result (Fig. S7). Of course, additional high-resolution cases of induced seismicity would be helpful for further scrutiny. The observation that induced earthquakes trailing the stimulation follow Båth’s law, both for expected magnitudes versus population counts, and for their residual distribution suggests that this approach should be useful for evaluating trailing seismicity.

Discussion

In this section we discuss the interpretations and implications of our results. We divide the discussion into three parts: (1) the applicability of Båth’s law for induced seismicity, (2) the relative importance of parameters in Båth’s law, and (3) a diagnostic test for operators to use in situ when encountering seismogenic faults.

The applicability of Båth’s law for induced seismicity

To date, many models have been introduced to capture the anticipated maximum magnitude of induced earthquakes from either physical or statistical perspectives[7,29,37]. Typically, these models have focused on the earthquakes occurring during the stimulation process. In this paper, we have described a statistical framework to extend these considerations to those earthquakes that trail after well shut-in based on Båth’s law for tectonic earthquakes. This conceptual extension is partly driven by the observation of large earthquake magnitude jumps[13], high-profile cases of large trailing events[5,11,12,48], and the contribution of trailing events to earthquake risk management[9].

Our results suggest that Båth’s law is a viable foundation for anticipating the magnitude of earthquakes trailing an operation that produces anthropogenic seismicity. This assertion is supported by the observation that: (1) induced seismicity cases are well fit by the linear relationship between magnitude differences ΔM and logarithmic population ratios that is predicted by Båth’s law (Fig. 3a), and (2) the residuals from this relationship are approximately consistent with the power law distribution anticipated by Båth’s law (Fig. 3b). Based on this, we further suggest that the description of Båth’s law (Eq. 5) will have predictive power for quantitative forecasts of trailing seismicity. Of course, this is predicated on knowing/estimating the population ratios —either a priori or updating in real-time during an operation. In this sense, the empirical distribution of and its best fit to a beta distribution (Fig. 2) provides some utility—possibly as prior information in a Bayesian updating scheme. Much of the literature that has been described for tectonic aftershock forecasts through extreme value statistics and Bayesian approaches[22,75], will likely be adaptable for the purpose of constraining N. Additionally, models of stimulation seismicity that consider the full injection history should provide better predictability to the anticipated stimulation events N[50].

This interpretation is in line with recent work that has found that induced earthquake sequences have statistical properties that are similar to (or subtle modifications of) their tectonic counterparts[29,36]. These results could be indicative of the importance of fluid-flow along faults for tectonic earthquakes, which has been suggested in some studies as important factor for the aftershock process[76,77]. General links between statistical models and physical processes have been made for earthquakes and aftershocks[78]. Thus, incorporating and testing physical models of aftershocks that consider realistic fault geometries, fault rheology, stress, material heterogeneity, and fluid flow will likely be important for better understanding post shut-in induced earthquakes.

Relative importance of Båth’s law parameters

A better understanding the physical processes that underlie the observed statistics would allow for better understanding of trailing seismicity. While the collected dataset (Table S1) is too small to meaningfully discern systematic trends in that vary with anthropogenic type, operational parameters, fault architecture, or geological setting, we may still apply our simple model to better understand the relative importance of parameters that contribute to a skewed distribution of . We focus first on the exponential decay model (Eq. 8) both for its simplicity and the breadth of induced seismicity studies it has been applied to[79,80]. If we assume representative values for a hydraulic fracturing operation (ΔT = 1 h,  = 3 days,  = 1.0, T = 2 weeks), we find that  = 82%; a value between those reported for the mean/median of the empirical . To build on this further, we assume that the parameters ΔT and follow a lognormal distribution while follows a Gaussian distribution so that we may compute a synthetic distribution. From this, we perform a Markov Chain Monte Carlo (MCMC) optimization to maximize the Kolomogorov–Smirnov p value between the empirical and synthetic . The MCMC approach randomly samples the parameter space to find optimal mean and standard deviations of these parameters (ΔT, , ) that are bounded within those reasonably observed in the literature[2,79].

The MCMC approach produces synthetic distributions that fit the empirical distribution well (likelihoods of up to 90%). The input parameters distributions needed to fit well (Fig. 4) suggests that a skewed distribution of input parameters (ΔT or ) is required to match the skewed distribution. Such a skewed distribution includes chances for black swan occurrences of extreme value events. In practice, this is what is observed in cases of induced seismicity: relatively rare events of large consequence. For example, the earthquakes at Pohang were exceptional in that a dormant fault system was reactivated by the geothermal stimulation program[5]. Similarly, most events in a 2018 episode in British Columbia occurred as part of a delayed trigger after hydraulic fracturing stage completion[48]. The largest event during the second well stimulation at Preston New Road occurred 60 h after the end of the last stimulation[50]. Earthquakes were apparently induced at distances of up to ~ 5 km away from the enhanced geothermal stimulation at Strasbourg-Vendenheim[81]. These consequential trailing seismicity events suggest the importance of better modelling and constraining parameters relevant to Båth’s law. We note, however, that there are significant limitations to our MCMC approach, as the distributions of ΔT and trade-off with one another. Thus, empirical observations of aftershock parameters ( or c and p) and injection responses (Σ, b, ΔT) would be needed to better constrain these values/distributions.

Figure 4

Synthetically fitting the distribution to model parameters. (d) The archetypal model for trailing seismicity ratios (Eq. 8) is fit (grey bars) to the empirical distribution of (blue bars). The distribution of model parameters (a) ΔT, (b) f, and (c) are also shown.

The Seismogenic fault injection test (SFIT)

Previously we showed that skewed values of ΔT or are required for a skewed —and represent an increased potential for large trailing magnitude differences ΔM. Given the possibility of consequential impacts from large magnitude trailing events, it would be prudent to measure these parameters in situ to reduce the uncertainties associated with fault reactivation.

Similar to Diagnostic Fracture Injection Tests[82] where fluids are intentionally injected into a reservoir to constrain its hydrological properties and stress state, we suggest an analogous Seismogenic Fault Injection Test (SFIT) to calibrate the seismic response to fluid injection (Fig. 5). Ideally, an inventory of faults would be available from a pre-assessment to injection. Injection would be performed with the intention of measuring both the time-delayed response to injection (Σ, b, ΔT) and the aftershock-like decay of seismicity ( or c and p) for each of the known faults[79]. The SFIT should be designed to constrain parameters of greatest consequence (ΔT and ). In the context of our study, unknown faults activated at distance from the stimulation interval (e.g., large ΔT) or with prolonged trailing sequences (e.g., large ) pose the greatest risks. With this in mind, an individual injection phase would span a range of injection rates, potentially sweeping through them cyclically in a ‘chirped’ manner (Fig. 5); the intent would be to both discern any rate-dependent features and better highlight potential time lags ΔT. As well, individual injection phases should be repeated to bolster confidence in the measurements, in particular the decay rate which only is well constrained after shut-in[79]. Ideally, these repeated injections phases would be spaced at non-periodic intervals to unambiguously measure delay times ΔT[83,84]. Finally, if we are making suggestions outside of what the simple models predict, varying the shut-in procedure for individual injection phases would be useful to discern safer responses to unwanted seismicity.

Figure 5

The Seismogenic Fault Injection Test (SFIT). Schematic diagram of a diagnostic test to measure a fault’s seismogneic response through parameters such as Σ, b, ΔT, and (or c and p). Repeated injection tests (at non-periodic intervals) are performed to verify the accuracy of measurements. Inset diagram shows how individual tests could be ‘chirped’ in their injection rate to better constrain parameters.

Conclusions

In summary, we have tested the applicability of Båth’s law to cases of induced seismicity with clear ends to fluid injection. Båth’s law fits the expected relationship, as confirmed by a linear regression of the collected data. Additionally, the residuals to the data follow the expected power law distribution, as confirmed by Kolomogorov–Smirnov tests. Thus, our compilation of population ratio statistics, including a fit to the beta distribution, could be useful for forecasting trailing seismicity. Finally, simple point-like models describing highlighted parameters of interest (ΔT and ); we propose in situ tests to help quantify the risks associated with seismicity trailing an anthropogenic operation. Variables used throughout the paper are summarized in Table 2.

Table 2

Summary of equation parameters and variables.

Variable	Description
a	GR-MFD productivity parameter, a-value
b	GR-MFD slope parameter, b-value
M, ΔM	Earthquake magnitude and magnitude difference, respectively
M_C	Magnitude of completeness
Σ	Seismogenic Index
Vt, V˙t	Cumulative injection volume and injection rate, respectively
nt	Earthquake rate
N, N_A,N_T,N_S	Cumulative earthquake counts: total, aftershocks, trailing, and stimulation
R_S, R_TS	Earthquake count ratios: related by 1/RS = N/NS = 1 + NT/NS = 1 + RTS
u, u_T,u_S	Confidence level variables
ΔT	Simple archetype model, stimulation-earthquake response delay
T_S	Stimulation/injection time interval
f	Stimulation-trailing rate change factor
K	Omori aftershock productivity parameter
p	Omori aftershock decay exponent
c	Omori aftershock singularity parameter
τ	Exponential aftershock mean decay time

Here we list the relevant parameters used in equations throughout our study, for convenience to the reader.

Supplementary Information 1.

Supplementary Information 2.