Event statistics

In this study, we used the induced event file maintained by the KNMI (https://www.knmi.nl/kennis-en-datacentrum/dataset/aardbevingscatalogus) containing all induced earthquakes recorded since 1986 in the Netherlands. Not all induced events occurred within the Groningen gas field. We therefore restricted the catalogue to the Groningen events by imposing a simple latitude-longitude filter (53.0931-53.4909N and 6.5516-7.1048E). By using such a simple filter, we probably lost a few Groningen events and included a few from nearby fields, but this will not greatly influence our statistical arguments below. We also discarded all events before 2002 as the instrumentation was sparse before that date. Rather than looking at magnitudes, a more physically meaningful quantity is the seismic moment, which is proportional to seismic energy release. We therefore converted the local magnitudes m _l given by the KNMI into seismic moments using log₁₀ M ₀ = 1.5m + 9.105, where M ₀ is the seismic moment in [Nm] and m the moment magnitude (Hanks & Kanamori, Reference Hanks and Kanamori1979). Before doing so, we had to change local magnitudes m _l into moment magnitudes m using the expression in Dost et al. (Reference Dost, Edwards and Bommer2018). This expression is valid for local magnitudes between 0.5 ≤ m _l ≤ 3.6. Restricting our catalogue further to events with a magnitude greater than 0.5 leaves us with 1372 events. We evaluated the probability distribution of the seismic moment using the Gaussian kernel density estimator from the SciPy statistics package (https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gaussian_kde.html). We see (Fig. 1) that the probability density of the seismic moment presents a power law decay over four decades with a slope α = 1.69 ± 0.01 giving the power law p(M ₀) ∝ M ₀ ^−α . Using known relations from seismology (Ben-Zion, Reference Ben-Zion2008), we can relate the slope of the seismic moment probability density to the Gutenberg–Richter law, which is cumulative, to give a Gutenberg-Richter slope b = 1.04 ± 0.02. This is close to what has been reported elsewhere (Dost et al., Reference Dost, Ruigrok and Spetzler2017). It is worth noting that the reported b-values vary in different regions of the gas field, which is most likely due to changes in compaction (Bourne et al., Reference Bourne, Oates, Elk and Doornhof2014), or changes in magnitude of completeness used in various parts of the region.

Fig. 1. Probability density of the seismic moment p(M ₀) for the Groningen catalogue using a cut-off magnitude of 0.5. The slope is 1.69.

Another statistic easily evaluated is that of interevent times. We follow Corral (Reference Corral2004) and estimate the probability density for interevent times, normalized by the average interevent time < t _i >, which is 5.06 days, using the same events. We observe (Fig. 3) that for longer times, the normalized interevent times follow a Gamma distribution $\gamma \left(0.64,1.70\right)$ resulting in a negative slope of γ = 0.36 ± 0.02 giving a power law behaviour p(t _i /<t _i >) ∝ (t _i /<t _i >)^−γ over two decades. This is slightly steeper than the values reported by Post et al. (Reference Post, Michels, Ampuero, Candela, Fokker, van Wees, van der Hofstad and van den Heuvel2021), but close to the universal values of Corral (Reference Corral2004). At shorter times, the slope is much steeper and is thought to be due to the presence of aftershocks following an Omori-type behaviour (Corral, Reference Corral2004). We find an Omori slope p = 0.82 ± 0.06. Although on the low end, this is within the range of expected values (Utsu et al., Reference Utsu, Ogata and Matsu’ura1995).

Both seismic moment and interevent times show a plateau at very low values of the respective variables. This is most likely due to the fact that the existing network, although dense, does not record all small events. There is a magnitude of completeness, above which all events are captured by the monitoring network. Dost et al. (Reference Dost, Ruigrok and Spetzler2017) estimated this magnitude of completeness to be above 1.2 for the Groningen area. If we now repeat the statistical analysis for a magnitude cut-off of 1.3, we are left with 694 events. The new probability density functions can be seen in Figs. 2 and 4. We now find α = 1.66 ± 0.02 and γ = 0.31 ± 0.05. The uncertainty in the slope of the intervent time curve is significantly larger, due to the dip at 10⁻³, which is likely due to the significantly smaller number of events using this higher cut-off. Within the uncertainties though, and certainly within two standard deviations, the slopes are the same. Introducing a magnitude of completeness eliminates the plateau in Fig. 2, while the interevent time curve remains largely unchanged (the Omori slope is now p= 0.87 ± 0.05). The average interevent time < t _i > for this higher magnigtude cut-off is 8.87 days. We checked that for other magnitude cut-offs, the power law parts of p(M ₀) and p(t _i /<t _i >) remained the same within the uncertainties, indicating that they are indeed scale-invariant.

Fig. 2. Probability density of the seismic moment p(M ₀) for the Groningen catalogue using a cut-off magnitude 1.3. The slope is 1.66.

Fig. 3. Probability density of the normalised interevent time p(t _i /<t _i >) for the Groningen catalogue using a cut-off magnitude of 0.5. The Omori slope is 0.82, and the Gamma slope is 0.36.

Fig. 4. Probability density of the normalised interevent time p(t _i /<t _i >) for the Groningen catalogue using a cut-off magnitude of 1.3. The Omori slope is 0.87, and the Gamma slope is 0.31.

Implications of the scale-invariance

In Benzi et al. (Reference Benzi, Castaldi, Toschi and Trampert2022), we showed that if p(M ₀) and p(t _i /<t _i >) are scale-invariant, α and γ are related for two specific cases of energy loading in various systems. Let us generalise this result. Assume that an earthquake releases the elastic energy that has been stored in the system over time. After the earthquake released it, energy will start being stored again and we assume that it grows as some power of the normalized interevent time. The released energy is proportional to the seismic moment (Ben-Zion, Reference Ben-Zion2008). We can then consider a new random variable $X={E_{{\rm stored}} \over E_{{\rm released}}}={\left(t_{i}/\lt t_{i}\gt \right)^{\varepsilon } \over M_{0}}$ . This constant power ϵ can be interpreted as the average over the considered time span 2002–2021 and does not model possible changes due to variations in the volume of gas extraction over the same time span. Indeed Bierman et al. (Reference Bierman, Paleja and Jones2015) analysed the rates of earthquakes in Groningen and found a possible correlation with changes in gas extraction for events of magnitude lower than 1.5. Using the ratio distribution rule, p(t _i /<t _i >) ∝ (t _i /<t _i >)^−γ and p(M ₀) ∝ M ₀ ^−α , the probability density for X is

(1) $$p\left(X\right)=X^{{{{1 \over \varepsilon}}} -{\gamma \over \varepsilon} -1}\int M_{0}^{{{1 \over \varepsilon}} -{\gamma \over \varepsilon} -\alpha }dM_{0}, $$

where the integration is over the part where the seismic moment shows a power law. Now let us consider the scale transformation

(2) $$M_{0}\rightarrow \lambda M_{0} $$

We then find that

(3) $$ p\left(X\right)\rightarrow \lambda ^{{1 \over \varepsilon} -{\gamma \over \varepsilon} -\alpha +1}p\left(X\right). $$

Since both p(M ₀) and p(t _i /<t _i >) are scale-invariant, p(X) is too. This implies that the exponent of λ has to be 0 and

(4) $$ \gamma +\varepsilon \alpha =1+\varepsilon. $$

Using the measured slopes for the seismic moment and interevent time above, we find that for a magnitude cut-off of 0.5, ϵ = 0.93 ±0.03and ϵ = 1.05 ± 0.04 for events above the magnitude cut-off of 1.3. For all practical purposes, we infer that the elastic energy within the Groningen gas reservoir roughly increases linearly with time until a new event occurs, i.e. ϵ ≈ 1. It is remarkable that this loading is very similar to that of tectonic earthquakes (Benzi et al., Reference Benzi, Castaldi, Toschi and Trampert2022).

Modelling interevent time using gradient boosted trees

Except for very short interevent times, p(t _i /<t _i >) closely follows a Gamma distribution (Figs. 3 or 4). A Gamma distribution could be an indication of some correlation between magnitude, space and time (Corral, Reference Corral2005). Furthermore, a clear Omori slope is present at short interevent times. Another common assumption is that some events are truly independent and follow a Poisson process, resulting in an exponential distribution. An important question therefore is how many aftershocks are following Omori’s power law, how many events a Gamma distribution and how many independent events an exponential distribution in the Groningen catalogue. Without explicitly distinguishing between aftershocks (Omori) and otherwise correlated events, the amount of correlated events was estimated to be between a few percent up to 27% (Luginbuhl et al., Reference Luginbuhl, Rundle and Turcote2018; Bourne et al., Reference Bourne, Oates and Elk2018; Candela et al., Reference Candela, Osinga, Ampuero, Wassing, Pluymaekers, Fokker, van Wees, de Waal and Muntendam-Bos2019; Muntendam-Bos, Reference Muntendam-Bos2020; Post et al., Reference Post, Michels, Ampuero, Candela, Fokker, van Wees, van der Hofstad and van den Heuvel2021). Most of the studies based on statistics assume stationarity to perform some form of declustering, except for the latter. We propose a machine learning route to interpret the interevent time distribution directly without any assumptions, except that the observed distribution is a mixture of Omori distributed, Gamma distributed and exponentially distributed events.

We generated 10⁵ samples from a power law, 10⁵ samples from a Gamma distribution, whose coefficients have been taken from the interevent times above, and 10⁵ samples from an exponential distribution. We then uniformly draw a random number n _g ∈ [0,1] representing the fraction of samples from the Gamma distribution and n _e ∈ [0,1] representing the fraction of samples from the exponential distribution. If n _g + n _e ≤ 1, then n _p = 1 − n _g − n _e is the fraction of power law samples. If however, n _g + n _e > 1, we assume that the interevent distribution is dominated by either the Gamma or the exponential distribution, i.e. if n _g > n _e , n _e = 0 and n _p = 1 − n _g , else if n _g < n _e ,n _g = 0 and n _p = 1 − n _e . We then generate a set of in total 10⁴ samples randomly composed of a fraction n _g Gamma distributed samples, a fraction n _e exponentially distributed samples and a fraction n _p power law distributed samples. Similar to the Groningen event data, we estimate the probability density for this set using a Gaussian kernel density estimator. The aim is to infer the fractions n _g , n _e and n _p directly from the probability density curve. To solve this regression problem, we used gradient boosting trees (Friedman, Reference Friedman2001) implemented by the Python machine learning package (https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.GradientBoostingRegressor.html). We created a training set of 5000 such probability density curves, which serve as input to the gradient boosting regressor and the fractions n _g , n _e and n _p are the target.

Splitting the training set into 80% training data and 20% test data, we show the performance of the trained gradient boosting regressor on the test data not previously seen by the regressor model (Fig. 5). The output points nicely fall on the diagonal, indication that the model has successfully learned the sample fractions from the input data. Another indication that the learning was successful is that the R ²-score for the test data is 0.98, similar to the score for the training data. We experimented with increasing the training set and slightly changing the parameters of the input distributions, but results remained similar. As can be seen in Fig. 5, the model is not perfect and there is some spread around the diagonal, which will result in some uncertainty in the inference from the observed interevent time distribution. We can directly use this spread to infer the uncertainty on our inference by estimating the standard deviation of the points around the diagonal.

Fig. 5. Performance of the gradient boosting random forest for a data set of 1000 unseen distributions.

Regardless of the cut-off magnitude, we find n _p = 0.03 ± 0.02, n _g = 0.92 ± 0.03 and n _e = 0.03 ± 0.03. It is not surprising that this inference is independent of the magnitude cut-off as the interevent time distribution is scale-invariant. Figure 6 shows that indeed the inferred fractions closely match the observed interevent time distribution, except for the plateau at the shortest times, which has not been modelled.

Fig. 6. Hundred realizations of the interevent times using the inferred proportions of the Omori, Gamma and exponential distributions within their standard deviation. The observed interevent time distribution is shown by red dots and coresponds to a cut-off magnitude of 0.5. Note that the narrow plateau at very short times has not been included in the modelling.