Interevent-time distribution and aftershock frequency in non-stationary induced seismicity

The initial footprint of an earthquake can be extended considerably by triggering of clustered aftershocks. Such earthquake–earthquake interactions have been studied extensively for data-rich, stationary natural seismicity. Induced seismicity, however, is intrinsically inhomogeneous in time and space and may have a limited catalog of events; this may hamper the distinction between human-induced background events and triggered aftershocks. Here we introduce a novel Gamma Accelerated-Failure-Time model for efficiently analyzing interevent-time distributions in such cases. It addresses the spatiotemporal variation and quantifies, per event, the probability of each event to have been triggered. Distentangling the obscuring aftershocks from the background events is a crucial step to better understand the causal relationship between operational parameters and non-stationary induced seismicity. Applied to the Groningen gas field in the North of the Netherlands, our model elucidates geological and operational drivers of seismicity and has been used to test for aftershock triggering. We find that the hazard rate in Groningen is indeed enhanced after each event and conclude that aftershock triggering cannot be ignored. In particular we find that the non-stationary interevent-time distribution is well described by our Gamma model. This model suggests that 27.0(± 8.5)% of the recorded events in the Groningen field can be attributed to triggering.

1. Gamma IAFT model details and maximum-likelihood equations 1

.1. Notation and survival analysis.
In this Supplementary Information we use the standard notation in probability theory; random variables are represented by capital letters (e.g. X), while their realizations are written as lower-case letters (e.g. x). A probability is denoted as P(·), the expectation as E(·) and H0 represent a statistical null-hypothesis. Estimates of parameters will be indicated by a hat,·.
We will use (and shortly state) standard concepts from survival analysis. The mathematical proofs and definitions used are well described in several textbooks, e.g. (1), but will be mostly omitted here. In survival, or failure time, analysis the central concept is the hazard rate, the instantaneous event rate, equal to where S(u) = P(U > u) = 1 − F (u) is the survival function. In statistical seismology this instantaneous event rate is often referred to as the intensity function.
The following equivalence relations are of great importance. Let H(u) equal u s=0 h(s)ds, the integrated (or accumulated) hazard, and f (u) represent the probability-density function (pdf, f (u) = − dS(u) du ), then

Gamma
Instantaneous Accelerated-Failure-Time model. The pdf and survival function of a random variable U0 that follows a Gamma(τ0, k) distribution equal . (6) where Γ(k, s) = ∞ s x k−1 exp(−x)dx is the upper incomplete Gamma function. Equivalently, the hazard function equals Finally the integrated hazard function equals H0(u) = log (Γ(k)) − log Γ(k, u τ0 ) .
Let t equal the process time starting at t = 0. In this research we model the hazard by a Gamma hazard changing over time as a result of a time-varying scale parameter τ (t). More specifically, to model the interevent times in the Groningen field we model the Gamma scale parameter as a function of (time-varying) local characteristics. The scale parameter represents the inverse of the background rate and equals Since, assuming τ (t l + u) > 0, lim u→∞ h(u) − 1/τ (t l + u) = 0, the interpretation of 1/τ (t) as the background rate is indeed appropriate for this time-varying model. The IAFT model introduced can also be fit for more-parameter distributions of which the long-time hazard equals a Poisson hazard, i.e. converges to a constant. In the Groningen case study we have only access to the monthly compaction rate and therefore the scale parameter is constant per month. The background rate is thus a discontinuous function with jumps at the end of each month. Let n(u + t l ) equal the number of mutually exclusive months at time u after t l . In this context, let sj equal the time in days, after t = t l , at the end of month j, with the exceptions that s0 = 0 and s n(u+t l ) = u. Then, = n(u+t l ) j=1 log Γ(k, The local occurrence probabilities, P(event occurred in region x | event occurred at time t), equal The local (region specific) hazard at time t in region x is defined as hx(u, t l ) = w(x, t l + u)h (u). Finally, the log-likelihood of the observed interevent times and event locations in the Groningen field equals * ll((u1, x1), ..., (un, xn)) = log (f ((u1, x1), ..., (un, xn))) where n is the total number of events, ui, respectively xi, are the interevent time and the region of the i th event and t l,i = i−1 j=1 uj.

Maximum-likelihood estimates.
The parameters (k, τ0, β) of the Gamma IAFT model are estimated by maximizing the log-likelihood defined in Eq. (13). These maximum-likelihood estimates are hard to find numerically as a result of the existence of local maxima for our complex log-likelihood. If τ (t) is constant over the course of each interevent time, then the log-likelihood simplifies to Now, the joint distribution is a member of the exponential family. Therefore, the log-likelihood is concave and a unique maximum exists, see e.g. Appendix A in (3). Finding the parameters such that the gradient of the log-likelihood equals zero gives rise to this unique maximum. The components of the gradients can be derived and are equal to ∂ll((u1, x1), ..., (un, xn)) * In this research we ignore the ongoing (censored) interevent time at the end of the catalog since its influence was negligible. More generally, the log-likelihood for a process with censored times and ψ(x) = d dx log(Γ(x)) equals the Digamma function. Moreover, the Hessian of Eq. (14) is equal to ∂ll ((u1, x1), ..., (un, xn)) ∂k∂τ0 = n, where Our optimization strategy is as follows. We numerically find the roots of the gradient (here the Hessian is used). These parameter values are now used as the starting values to numerically (Broyden-Fletcher-Goldfarb-Shanno algorithm (4) is used) find the maximum-likelihood estimates for the complex model.

Determination of model covariates and parameter values for Groningen
In this study we have excluded events from the KNMI catalog based on a magnitude of completeness equal to 1.3 as suggested in (5).

Candidate covariates Groningen.
The covariates taken into account in the case-study analysis can be found in Table  S1. The stationary covariates are derived from the geological top Rotliegend surface model from the NAM provided via TNO (6). For the compaction covariates we use the compaction model provided by Shell (7) For stability of the numerical optimization of the likelihood the covariates are transformed. First the yF (x) is scaled by 10 −5 , yT (x) by 10 −1 , yĊ (x, t) by 10 5 and yC (x, t) by 10. Subsequently, all covariates are standardized by their mean (over the regions, see Table S1), e.g. yC (x, t) := 10yC (x, t) − 1.47. The numerical values of the stationary covariates obtained by this procedure can be found in Furthermore, we considered a critical total compaction by introducing a truncation c, mimicking the compaction level at which all faults in a region have reached their fault strength. The Gamma time-scale parameter used in the Groningen case study equals

Final model.
In statistics it is common to rely on a forward-or backward-selection procedure to select covariates that should be included in the model. However, the optimal model can easily be missed when covariates are highly correlated. We therefore fit all possible models including an effect of the total compaction (since this is the only covariate that can explain the increase in intensity over time). Each model is fit using the approach explained in Section 1.3. However, the derivative of ll ((u1, x1), ..., (un, xn)) with respect to the compaction-truncation variable c does not exist. The gradient can be derived given a value of c. For several values of the c variable we numerically find the roots of the gradient (here the Hessian is used). Subsequently, the different c fits are compared and the parameters giving rise to the highest log-likelihood are selected. In turn these estimates are used as the starting values of the numerical optimization. The final model is selected based on Akaike's information criteria (AIC) (8) to prevent overfitting. Below we present the (log-likelihood-based) best-fitting model with J covariates involved.
It is important to notice that the k estimate proves insensitive to the model choice, thus our estimate for the fraction of triggered events is insensitive to the model choice. Based on the AIC the model with covariates total compaction, compaction rate, total fault density, percentage of faults with specific strike angle, throw-thickness ratio and truncation is selected. The effect estimates and corresponding standard errors are presented in Table S4.
The final estimates define the development of τ (t), which evolution is displayed in Fig. S2. This background rate, 1/τ (t), is partitioned over the 15 regions using the weights of the Gamma IAFT model. For a subset of the regions, the evolution of the proportion of seismic events over time is displayed in Fig. S3.
The effect estimates and cap estimate presented in Table S4 can be transformed back such that they can be interpreted in the units of the original covariates introduced in Table S1, e.g. the cap variable in meters compaction is estimated as 0.94+1.47 10 with a corresponding standard error of 0. 044 10 , and are presented in Table 1 of the main text. The final model thus includes a   Table S3. Maximum-likelihood Gamma IAFT-model parameter estimates using J = 3, ..., 10 covariates respectively. truncation on the total compaction. This effect is shown for each region in Fig. S4. The Gamma IAFT model suggests that all faults in the most active region 10, the Loppersum area, were critically stressed at the end of 2007, when the local compaction was 0.24 meter.

Risk of bias due to successive-event overlaps in seismograms.
In real catalogs, short interevent times cannot be detected because waveforms of successive-event overlaps in seismograms. As a result it is impossible to observe short interevent times below some threshold. The observed interevent times are thus realization of a conditional distribution (random interevent time U > u0, for some threshold u0). If the probability that U ≤ u0 is reasonably large, then this conditional distribution seriously deviates from the full distribution. Parameter estimates obtained by fitting the non-conditional distribution will then be biased. Instead one can fit the conditional distribution by maximizing where fi, Si, Hi are the pdf, survival function and integrated hazard function of the i th interevent-time distribution. The potential bias is the result of the n i=1 Hi(u0) contribution that lowers minus the log-likelihood. It is thus important to evaluate n i=1 Hi(u0) after the model has been fitted. If the contribution is low, then the conditional distribution is close to the full distribution and the thresholding will not result in bias of the estimates.
In the Groningen catalog, which we used to fit the final model (Tables S3 and S4), the smallest interevent time was equal to 0.00069 days (1 minute). From Groningen seismograms, see e.g. (9,10), it becomes clear that the time between the first-P-wave arrival at a seismic station and the moment when the signal amplitude has decayed below the first-P-wave amplitude is in the order of seconds. Thus, it is reasonable to assume that all interevent times larger than or equal to 1 minute are observed. Using the model fitted in this study we find that 397 i=1 Hi(0.00069) = 0.2, which is extremely small compared to the final minus log-likelihood of 2335.8. In this case we thus do not suffer from time-threshold bias. In the case the thresholded distribution does deviate from the unconditional distribution one should maximize Eq. (27) instead of Eq. (13). As an example we present the thresholded-distribution fit using u0 = 5 minutes, such that we fit 395 interevent times, which estimates are presented in Table S5. Indeed the estimates are similar to those presented in Table S4.

Comparison with the memoryless IAFT model.
To test the presence of earthquake-earthquake interaction we compare the final Gamma IAFT-model fit with the fit of an Exponential (k = 1)-model. The maximum-likelihood estimates of the parameters of the Exponential model equal: Based on the likelihood-ratio test (p < 0.0001) we reject the hypothesis H0 : k = 1 and conclude that the improved fit cannot be explained by chance alone. It is important to note that both models are fitted to the Groningen dataset and thus roughly expected the same number of events in the period October 1 st 1995 -October 1 st 2018. The difference between the models is clearly illustrated by plotting the temporal hazards over time, see Fig. S5. The hazard of the Gamma IAFT model is far more spiky than the hazard of the Exponential IAFT model. The Gamma IAFT model results in a more clustered profile of events as illustrated in Fig. 3a.

Model validation with simulated Groningen data
The final model obtained in Section 2.2 for the Groningen field results in an estimate of the non-stationary driving function 1/τ (t) and an interevent-time distribution after each event time t l . Since in the Groningen case the time scales of variation of τ vs. t and h vs. u turn out to be well separated, the instantaneous distribution, following each t l , can be approximated by a Gamma distribution. As a result the time variation of τ can be scaled out, giving a single Gamma, Γ(k, 1), distribution vs. the variable U/τ (t l + U ) over the full time period. This function is shown in Fig. S6 together with the actual rescaled, byτ (t), Groningen interevent data. The pdf of (Wald-Wolfowitz runs test, p = 0.17), see Fig. S6. Thus, once we account for the change in intensity over the years the mixed distribution of background and triggered interevent times can be well described by this universal Gamma curve. However, since the data are scarce, the agreement between curve and data is in itself no sufficient validation of our model, and more formal validation tests will be presented below.

Simulation description.
To simulate past-event catalogs based on the final model we rely on the theory of survival analysis. In general a random variable with cumulative distribution function (cdf) F is generated by drawing a realization of a Uniform [0, 1] random variable V and evaluating the inverse cdf at this V , F −1 (V ). Equivalently, one could evaluate the survival function, S = 1 − F , at the random V . From survival analysis we know F (u) = 1 − exp (−H(u)). Our Gamma IAFT model defines the hazard rate, and thus the integrated hazard H(u), presented in Eq. (11), and the cdf.
The maximum-likelihood estimates of the final model theoretically follow an asymptotic multivariate normal distribution † , see e.g. (11). More precisely whereθn are the maximum-likelihood estimates for (k, β1, ..., βJ ) based on a sample of size n, θ = (k, β1, ..., βJ ), d → represents convergence in distribution and I θ equals the Fisher information matrix. Letθ represent the maximum-likelihood estimates of the Gamma IAFT-model parameters. We estimate the total varianceΣ, by computing (numerically) the inverse of the Hessian of the log-likelihood evaluated at the maximum-likelihood estimatesθ. The first step in the simulation of an earthquake catalog consists of drawing a realization θ of the multivariate normal N θ ,Σ .
Subsequently, a clock is set to t = 0 (October 1995 in the Groningen case study) and a interevent time U1 is realized. Furthermore, a multinomial random variable, with probabilities equal to the Gamma IAFT-model weights at time U1, is realized to determine the event location. Next, the clock is set to U1, after which a second interevent time, U2, is realized based on the integrated hazard after U1 and the location of this event is generated based on the weights at time U1 + U2. These steps continue until N events are realized such that Thousand generated catalogs based on the models fitted to the interevent times from the Groningen catalog (Gamma IAFT and Exponential IAFT) are shown in Fig. S7. † This theory holds for the model presented as Eq. (10) and Eq. (12). However, for model Eq. (26) used in the Groningen case, θ = (k, c, β 1 , ..., β J ) and the theory breaks down since c → min{c, y C (x, t)} is not differentiable w.r.t to c at y C (x, t). In practice the parameter estimates will not change if we use a smooth approximation of the min function. Using this approximation the theory does apply and thus allows us to rely on the discussed result by numerically deriving the inverse of the Hessian of the log-likelihood evaluated at the maximum-likelihood estimatesθ.  The empirical distribution of the Cox-Snell residuals can be compared to the unit Exponential distribution using the Kolmogorov-Smirnov (KS) test (13). For the final model the p-value of this test is equal to 0.681, from which we conclude that the proposed Gamma IAFT model seems suitable. The p-value of the KS test comparing the unit Exponential with the Cox-Snell residuals of the Exponential IAFT model instead equals 0.009, as a result of the underestimation of the clustering. This explains the superiority of our final IAFT model compared to the Exponential IAFT model once more.

Randomness.
The distribution of the Cox-Snell residuals can be used to test for overall performance. However, in non-stationary processes it is very important to investigate the randomness in the sequence of the residuals. E.g. in a setting where the Cox-Snell residuals are expected to be smaller at the beginning of the catalog and larger towards the end, the proposed model does not describe the underlying mechanism of the data. To test whether the Cox-Snell residuals (sorted on event date) are independent and identically distributed we apply Wald-Wolfowitz runs test (13). The resulting p-value equals 0.159, from which we conclude that there is no reason to doubt the appropriateness of the model.

Spatial fit.
We have stressed that appropriate modeling of the spatial heterogeneities is key to model the temporal clustering. To validate the spatial fit we use Pearson's chi-squared test (13). More specifically, the test statistic equals .
The expected counts, E Gamma IAFT [Xi], per region are estimated from the thousand simulations and presented in Fig. S9 (left). For the Groningen data the statistic equals χ 2 = 37.8 and the contributions per region are presented in Fig. S9 (right  The distribution of χ 2 assuming the Gamma IAFT model is approximated using the thousand simulations. Based on this empirical distribution we find P(χ 2 > χ 2 data ) = 0.048.
This probability is rather low and thus questions the spatial fit of our model. However, with a closer look at the expected and realized counts in chi-square statistic of χ 2 = 24.5, with a p-value of P( χ 2 > χ 2 data ) = 0.24. This shows that our model is capable of reproducing the spatial distribution over the field.

Hindcast.
3.3.1. Gamma IAFT model. Finally, we investigate how well we can forecast the number of events in the period 2014-2018 based on the information available before that period. This is of particular interest since the production of gas has drastically decreased from this period on, giving rise to compaction rates that have not been observed before. When we want to interpret our claims causally we should be able to predict event counts in a production scenario that has not been observed so far. The maximum-likelihood estimates of the Gamma IAFT model based on the period 1995-2013 are given in Table S7. Simulations for the 2014-2018 period of the Gamma IAFT model with the parameters from Table S7 are presented in Fig. S11.

3.3.2.
Ignoring spatial heterogeneity. The importance of taking care of the spatial heterogeneity can be illustrated by studying the hindcast with a model with no spatial (ns) heterogeneity and only describing the temporal process. Now, we model the scale parameter of the Gamma IAFT model, τns, as The maximum-likelihood estimates of this temporal model based on the period 1995 − 2013 can be found in Table S8. In this case, the compaction truncation did not improve the model fit. The simulations based on this model can be found in Thus The Gamma IAFT model can be used to stochastically decluster an earthquake catalog by labeling events as a triggered or background event based on random (Bernoulli) experiments.

Theorem 2. The probability of an event to originate from the background process, given the interevent time u, equals
where h k=1 (u) = 1 τ (t l +u) and equals the Gamma IAFT hazard at u when setting k equal to 1.
Proof. Let U = min(U0, U1), be the random interevent time starting at t l where U0 and U1 equal the times to the next background event and to the next triggered, by earthquake-earthquake interaction, event respectively. Here U0 is the interevent time, after t l of an inhomogeneous Poisson process with (background) rate 1/τ (t l + u). Using Bayes' theorem and the survival-analysis equivalence relation Eq. (4), given U = u, we find h 0 (u)+h 1 (u) and substituting h0(u) = 1/τ (t l + u) gives the result.
By drawing a Uniform random variable, V , for each event we label the event as a background event when V < p background (ui) and as a triggered event otherwise. To decluster the catalog of Groningen we start by drawing Gamma IAFT-model parameter realizations from the multivariate normal distribution discussed in Section 3.1 and start labeling each event (background or triggered) using p background (ui) given θj. Some examples of stochastically declustered catalogs of Groningen can be found in Fig. S13. At the end of the catalog (October 2018) the declustered catalogs contain on average 73% of the events.

Fraction of background events.
We would like to compute the % of background events in the Gamma IAFT model. To do so we are interested in the expected value of the background-event probability, the ratio between the background rate and the total hazard of the Gamma IAFT model, at interevent time U (after t l ). The complex interevent-time pdf can be derived using the survival-analysis equivalence relations whereh(u),f (u) andS(u) equal respectively the hazard function, pdf and survival function of the stationary Gamma distribution with scale parameter τ (t l + u). Since the pdf and cdf have no elegant closed form expression it is hard to derive Let us assume there are a finite number of intensity changes over the period of interest. Now, we can define the n intensity changes at time s1, s2, ..., sn after t l . During the interval [t l + si, t l + si+1) the background rate is constant and can be referred to as 1 τ i . After time t l + sn, the background rate stays constant at the level of 1 τ n+1 . Without loss of generality we assume s0 = 0. Then the expected proportion of background events can be computed as If the scale parameter τ (t) does not change during each interevent period, then each interevent time is Gamma distributed, with shape parameter k and a scale parameter that might differ per interevent time, and the expression simplifies to Γ(k+1) Γ(k) = k. If the scale parameter τ (t) does change, E [p background (U )] will differ from k, but can be derived via simulation. For each simulation 1 n j n j i=1 p background (ui) can be computed and by the Central Limit Theorem, the average over all simulations will converge to E [p background (U )]. For the Gamma IAFT model using the maximum-likelihood estimates, derived from the Groningen catalog, this fraction equals 0.73 and indeed equals k. This is a consequence of the differences between 1/τi and 1/τi+1 being small and thus the complex distributions are very comparable to (constant scale-parameter) Gamma distributions.

Confidence interval of the fraction of background events in Groningen.
In an induced-seismicity study, like Groningen, the stochastic process underlying the data can be modeled. The observed data can then be viewed as one specific realization of the stochastic model. Note that the fraction of background events can differ per realization of the stochastic model and can deviate from the expected fraction of the process (the latter has been discussed in the previous section). For the real data the event times are known and the fraction should thus be estimated conditioned on these event times. The i th event has a probability equal to p background (ui) to be a background event, thus the expected number of background events is equal to n i=1 p background (ui). For a Gamma IAFT model with the maximum-likelihood estimates, the expected number of background events is equal to the expected number of the stochastic process, E n i=1 p background (Ui) , discussed in the previous section, and thus equals 290 (73% of the total).
Each event label, Li = 1 event i is a background event , is Bernoulli distributed with the just-defined background probability. Therefore, the number of background events, N = n i=1 Li, follows a so-called Poisson binomial distribution with parameters p1, p2, ..., pn (14). Confidence limits of the number of background events in a realized catalog of a known Gamma IAFT model are thus equal to the quantiles of this Poisson binomial distribution. In practice the parameters of the Gamma IAFT model are unknown and therefore estimated, which introduces uncertainty. Let θ represent the real Gamma IAFT model parameters, then the distribution of N | θ is the Poisson binomial distribution just introduced. Now, let θ itself be a multivariate random variable with law f θ . Then, the distribution of N equals where Θ is the support of the θ variable, namely {k > 0, τ0 > 0, βC ∈ R, ...}. Let us assume that θ follows the asymptotic multivariate normal distribution of the maximum-likelihood estimate (of the Gamma IAFT model) discussed in Section 3.1.
Since it is analytically very difficult to compute the 2.5% and 97.5% quantiles we have used a Monte Carlo approach to derive the distribution of N . For j = 1 to j = 10000: 1. We have generated θj from N θ ,Σ .
2. We have generated Nj from the Poisson binomial distribution with p background (ui).
Finally, the 2.5% and 97.5% quantiles of the empirical distribution are computed. For the Groningen catalog (n = 397), we have concluded that the confidence interval (CI) for the number of background events equals (256, 324), which equals (64.5%, 81.6%) of all events. ‡ Note that the CI of k, based on the asymptotic normal distribution ofk, equals (66.8%, 79.2%) and is thus tighter. ‡ To verify that the empirical distribution represents the distribution of N we check that the CI based on only the first five thousand simulated catalogs, (255, 325), is similar to the CI based on all simulations.

Comparison with ETAS-based models and data
The Epidemic-Type Aftershock-Sequence (ETAS) model is often used in practice to analyze seismicity (15). In the ETAS model the occurrence rate of aftershocks at time t due to an event at time v decreases according to the modified Omori law K (c+t−v) 1+θ , where K, c and θ are model parameters (16). Each event of magnitude m triggers aftershocks with a rate proportional to 10 α(m−m min ) , where mmin is the magnitude of completeness. The hazard of the temporal ETAS model equals Here 1/τ (t) represents the background rate as in the Gamma IAFT model, n(t) equals the total number of past events before time t, vi and mi represent the event time and magnitude of the i th event. The ratio of background-to-total events can, only in steady-state, be directly related to the so-called branching ratio, g, the number of first-generation daughter events per background event (17). For θ > 0 and α < b, we obtain that g = K is the parameter of the Gutenberg-Richter distribution. The integrated hazard of the ETAS model equals The integrated hazard defines the survivor function S ETAS (t), which allows us to simulate with the ETAS model, following the same lines as explained in Section 3.1.

Groningen-based ETAS model.
In this study our focus is on the Gamma IAFT model rather than an ETAS model, since the small data base and the complexity of a non-stationary driving force asks for a parsimonious parameter set to avoid arbitrariness. We want to show that this parameter-poor Gamma model can nevertheless be effectively used on a Groningen catalog simulated with a temporal ETAS model. Based on the Cox-Snell test we assume that the τ0 and β estimates obtained from the Gamma IAFT fit to the Groningen catalogs are accurately describing the background rate. To find values for the parameters c, θ, α and g that give rise to synthetic catalogs that are representative of the real Groningen catalog we evaluate the log-likelihood  16)). In Table S9, per c value, the parameter sets that give rise to the two best fits and the corresponding log-likelihood values are presented. Accurate estimation of the four remaining ETAS parameters requires a sufficiently large number of observations. To illustrate that in the Groningen catalog this is not the case, we compared the two best fits, (c = 6h, α = 0.3, θ = 0.05, g = 0.8) and (c = 6h, α = 0.3, θ = 0.10, g = 0.4). Their log-likelihoods are insignificantly different and indicate a strong intercorrelation of θ and g. The resulting systems are however vastly different in their cross-over times for saturation of individual aftershock sequences (17), and hence in intersequence overlap. Consistent with the choice of the Gamma timescale τ (t) as the true input background rate, we selected the fit with g = 0.4, i.e., limited intersequence overlap (18). Synthetic ETAS-simulated catalogs mimicking the Groningen data are then obtained with c = 6 hours, θ = 0.10, α = 0.3 and g = 0.4 see Fig. S14. Based on 500 simulations our temporal ETAS model gives rise to a mean fraction of 76.6% background events over the course of the Groningen catalog, see Fig. S16a. (Indeed this does not equal 1 − g , as the system is not in steady state.) On a par with the Gamma IAFT model the temporal ETAS model did not consider individual event locations, although differences in compaction behavior over the Groningen field were accounted for via the used time-varying background rate.

Gamma IAFT-model fit to ETAS simulations.
To show the appropriateness of the Gamma IAFT model we first fitted it, following Section 1.3, to the union of all interevent times from all 500 simulations, by optimizing the log-likelihood: ll IAFT (u1,1, ..., u1,n 1 , u2,1, ..., u500,n 500 | k, β which gives in particular k = 0.746. Following the same procedure as in the paper and Section 3.2.1, we have compared the cdf of the Cox-Snell residuals from this Gamma IAFT-model fit, considering its maximum-likelihood estimates, to the cdf of a unit Exponential, see Fig. 3a (main text) and Fig. S15. We conclude that the interevent times generated by the temporal ETAS model are indeed well approximated by the Gamma IAFT-model distribution.

Background fraction estimation based on Gamma
IAFT model or ETAS model. As the fraction of aftershocks is a stochastic quantity (19), the 'true' expected background fraction, µBF (g, θ, c, β), differs per synthetic catalog and equals Fig. S16a. We would like to investigate the performance of the estimation of the expected background fraction by fitting the ETAS model or the Gamma IAFT model to the synthetic ETAS catalogs.

Efficiency.
For each synthetic catalog, we have compared the estimated background fractions to µBF (g, θ, c, β). A histogram of the difference is presented as Fig. S16b. The median of these errors equals 0.03 and −0.03 for the ETAS and the Gamma IAFT fit respectively. We conclude that, despite the finite sample size, the background rate and % of background events are appropriately estimated with the Gamma IAFT model. The variances of these background-fraction estimates are referred to as σ 2 ETAS and σ 2 IAFT . We findσ 2 ETAS = 0.05 and σ 2 IAFT = 0.04 respectively. The standard error is significantly higher when fitting an ETAS model (Brown-Forsythe test, H0 : σ 2 ETAS ≤ σ 2 IAFT , p < 10 −4 ). Hence, even for the synthetic Groningen-based ETAS catalogs, the fraction of triggered events can be estimated more efficiently by fitting a Gamma IAFT model.

5.4.
Power of the Cox-Snell test for validating Gamma behavior. Actual catalogs as well as synthetic ETAS catalogs do not always conform to Gamma interevent distributions (20). To show how sensitively our approach can validate the appropriateness of the Gamma distribution we simulated 1000 times 397 events for different stationary ETAS models. We used τ /c = 10 5 , setting the values of n and θ such that, depending on the branching ratio and the tail of the Omori function, both catalogs far from and close to the asymptotic steady state were obtained; the base case n = 0.9, θ = 0.03, which is still far from steady state, corresponds to Fig. 1 in (20). The settings give rise to distributions that deviate from Gamma behavior, see Fig. S17. Consequently, the (Gamma) maximum-likelihood estimators of the shape (k) and scale (τ ) parameters deviate from the real fraction of background events and the inverse background rate respectively. In these settings the maximum-likelihood estimator of the shape parameter,k, should thus not be interpreted as the fraction of background events.  Fig. S17. Empirical estimate of the pdf for U τ , using a logarithmic binning of all (1000 · 396) simulated interevent times with the ETAS model (black line), the Γ( % background, 1) pdf (cyan) and the Γ(k, 1) pdf (blue) for settings A-E described in Table S10.
Nevertheless, as shown in Table S11, the power (probability of rejecting the null-hypothesis) of the goodness-of-fit test (of the Gamma distribution), based on the Cox-Snell residuals discussed in 3.2.1, is very high for these scenarios. So the Cox-Snell test can sensitively verify the appropriateness of the Gamma model. The interpretation of the Gamma IAFT model presented in this research is only valid after such formal model validation.