Gas extraction from the Groningen gas field resulted in significant induced seismicity. We analyze the magnitude‐frequency distribution of these earthquakes in space, time and in view of stress changes calculated based on gas production and reservoir properties. Previous studies suggested variations related to reservoir geometry and stress. While we confirm the spatial variations, we do not detect a clear sensitivity of b‐value to Coulomb stress changes. However, we find that b‐value correlates positively with the rate of Coulomb stress changes. This correlation is statistically significant and robust to uncertainties related to stress change calculation. This study thus points to a possible influence of stress change rate on the probability of the magnitude of induced earthquakes.
The factors influencing earthquake magnitude is a subject of active research. While fault geometry (Stirling et al., 1996;T. H. Goebel et al., 2017) plays a role, it is commonly admitted that higher stresses promote larger magnitudes (Scholz, 1968(Scholz, , 2015;;Schorlemmer et al., 2005;Rivière et al., 2018). The effect seems visible in examples of injectioninduced seismicity (Mukuhira et al., 2021(Mukuhira et al., , 2024) ) and seismicity induced by gas extraction at Groningen (Bourne & Oates, 2020;Kraaijpoel et al., 2022;Muntendam-Bos & Grobbe, 2022). Here, we reanalyze the case of Groningen (Fig. 1a, b). The motivation for this reanalysis is that previous studies have only considered annually averaged stress changes, while in reality gas production is highly seasonal (Fig. S1). We use a geomechanical model calibrated with measurements of surface deformation to calculate stress changes within and outside the reservoir with subannual resolution (Meyer et al., 2023;Acosta et al., 2023) (Fig. S1). Hereafter, we present the setting and review previous studies of induced seismicity at Groningen. We then present our data and methodology and discuss our results.
The Groningen gas reservoir (Fig. 1) consists of a Permian porous (∼20%) sandstone at a depth of 3 km (de Jager & Visser, 2017; Bourne et al., 2014). Reservoir thickness increases from ∼90 m to the southeast to ∼300 m in the northwest (de Jager & Visser, 2017). Earthquakes were first detected in 1991. Seismicity rate increased until 2013, and then decreased as extraction slowed (Muntendam-Bos & Grobbe, 2022) (Fig. 1d). Seismicity rate in space and time is related to reservoir compaction and can be reasonably well forecasted (Bourne & Oates, 2017;Bourne et al., 2018;Candela et al., 2019;Richter et al., 2020;Smith et al., 2022;Dempsey & Suckale, 2023;Kaveh et al., 2023). Several studies have also examined the magnitude-frequency distributions (MFD) of these earthquakes (Kaveh et al., 2023;Zöller & Hainzl, 2022;Zöller & Holschneider, 2016;Bourne et al., 2018). Estimating magnitudes is of foremost importance (Zöller & Holschneider, 2016;Bommer & van Elk, 2017), as they feed directly into risk analysis.
The MFD of earthquakes is commonly quantified using the b-value (Gutenberg & Richter, 1944). The b-value gauges the number of smaller earthquakes relative to larger ones, with high b-values indicating higher relative frequency of smaller earthquakes (Fig. 1c). In seismic hazard assessment, it is common to assume a stationary magnitude probability in time but with spatial variations. Variations in space have been identified at Groningen (Kraaijpoel et al., 2022;Gulia, 2023;Boitz et al., 2024;Muntendam-Bos & Grobbe, 2022), but other studies reveal possible variations in time. Bourne and Oates (2020) and Kraaijpoel et al. (2022) also report lower b-values at higher stresses. Their findings are consistent with observations of b-value in the laboratory (T. Goebel et al., 2013;Rivière et al., 2018) and the negative correlation between b-value and differential stress observed for natural earthquakes (Scholz, 2015).
However, three lines of evidence potentially undermine the correlation between bvalue and stress at Groningen. First, the earthquakes appear to follow a spatial dependence, as opposed to a stress change dependence. For earthquakes recorded within the reservoir outline between 1991 and 2023, those associated with higher Coulomb stress changes (CSC) occupy the center of the reservoir, where there is greater compaction, while those with lower CSC occur towards the edges (Fig. 1b). While CSC follows this radial pattern, the b-values do not. Muntendam-Bos and Grobbe (2022), Gulia (2023), andKraaijpoel et al. (2022) observe b-value increasing from the northwest to southeast, following the pattern of reservoir thickness. The second line of evidence relates to temporal resolution. Bourne and Oates (2020) and Kraaijpoel et al. (2022) use the Elastic Thin Sheet model of Bourne and Oates (2017), calculating CSC on an annual resolution. However, gas extraction and CSC follow a seasonal cycle (Acosta et al., 2023) that is not accounted in either study. Third, b-value may instead be controlled by the rate of CSC. Gulia (2023) observed a negative correlation between b-value and compaction rate, a quantity that scales with stress change rate. We are thus motivated to revisit the relationships between b-value and stress in the Groningen reservoir. To benchmark our findings with the existing literature, we also compute variations with time and space.
3 Data and methods
We use the earthquake catalog of the Royal Netherlands Meteorological Institute (KNMI, 2023) and select earthquakes recorded between December 1991 and July 2023 within the reservoir outline. This catalog contains 1487 events with local magnitudes between -0.18 and 3.60. The magnitude of completeness (M C ) decreased from ∼1.50 to ∼0.50 between 1991 and 2014 as seismic monitoring improved (Dost et al., 2017;Smith et al., 2022).
Stress changes throughout the reservoir are calculated based on monthly pressure changes, the maximum temporal resolution of the available information on extraction data (Oates et al., 2022;Acosta et al., 2023). We first use a vertical flow equilibrium model, benchmarked with fluid pressure measurements at wells (Meyer et al., 2023), to calculate pore pressure changes in the reservoir. From pressure changes, compaction and stress changes are calculated using the geomechanical model of Smith et al. (2022) (see Text S1). The reservoir is modeled as poroelastically deforming cuboids. Each cuboid is ∼500 m in the horizontal directions and has a thickness equal to the local thickness of the reservoir. Probabilistic assessment of earthquake location by Smith et al. (2020) show that the earthquakes' depth distribution peaks just above the reservoir. We therefore sample the stress field at a depth of 5 m above the reservoir and at the center of each cuboid in the horizontal plane. We hereafter refer to this sampling scheme as the reference model.
The reference model provides a good basis to predict the spatial and temporal varia--3-manuscript submitted to Geophysical Research Letters tion of seismicity rate (Kaveh et al., 2023). To verify our results are robust to different sampling locations in our model, we test additional sampling schemes (see Table S1).
For each earthquake, we estimate the Coulomb stress change (CSC) using the normal and shear stress change calculated at the nearest time-step and sampling location.
We assume a friction coefficient of 0.6 and a fault orientation optimally oriented for CSC.
We also estimate the time derivative of CSC using backward differencing with the previous month. This derivative is hereafter referred to as the CSC rate. CSC rate follows a seasonal cycle, with the winters characterized by a higher rate due to greater extraction. However, CSC rate experiences a long-term decline from September 2015 (Fig. S1) due to a year-round decrease (Smith et al., 2022) and an important reduction of the seasonal swings in extraction activity (Muntendam-Bos & Grobbe, 2022). We therefore only use earthquakes occurring before September 2015, when seasonal variations are strong and probably well resolved.
We use three methods to calculate b-value:
1. Maximum likelihood method (Aki, 1965), assuming a M C of 1.20. Kijko and Smit (2012), accounting for time-varying catalog completeness. We estimate the time-evolution of M C from Smith et al. ( 2022) (see Fig. S2).
3. A method proposed by van der Elst (2021) that uses the absolute values of magnitude differences between sequential earthquakes, hereafter referred to as the "babsolute" method. Because method is insensitive to the M C of the catalog, the entire catalog can be used to better resolve b-value variations. We opt to use babsolute rather than b-positive, which uses only the positive magnitude differences and therefore only half of the data. The b-positive method alleviates the fact that, during an aftershock sequence, smaller events are obscured by the temporary increase in seismic noise due to coda waves (van der Elst, 2021). In our case, the proportion of aftershocks is small, making the possible bias negligible. We verify this claim by checking that the b-positive and b-absolute methods yield values that are generally consistent (Fig. S3).
These methods differ by how they address catalog incompleteness at low magnitudes. This effect must be taken into account, as variation in the M C with time could be a source of bias.
For an ideal earthquake catalogue following a non-truncated Gutenberg-Richter distribution, all methods should yield the same b-value within uncertainties. In practice, this might not be the case due to departures from the ideal log-linear MFD.
We compute variations in b-value with time, space, CSC, and CSC rate. We group the earthquakes in overlapping bins with respect to each variable, in ascending order, and calculate the b-value of each bin. For each bin, we randomly sample, without replacement, a subset of earthquakes and calculate the b-value of each subset. We remove any earthquakes below the M C or the magnitude differences between sequential earthquakes that are smaller than the catalog resolution. Even after sample removal, most bins and subsets end up at sizes of ∼150 and ∼100, respectively. We conduct 50 iterations of sampling, after which we calculate the median, 16th percentile, and 84th percentile of b-value across all samples. We designate the median as the b-value of that bin.
-4-manuscript submitted to Geophysical Research Letters When using the modified maximum likelihood method of Kijko and Smit (2012), we divide each subset into groups of differing M C . We calculate the b-value of each group if it contains at least 10 earthquakes. We then use the average of these b-values, weighted by the number of earthquakes per group, as the b-value of that subset.
All three methods feature a gradient of b-value increasing from northwest to southeast (Fig. S5). These observations are consistent with Muntendam-Bos and Grobbe (2022), Gulia (2023) and Kraaijpoel et al. (2022). This gradient may be attributed to the decrease in reservoir thickness from northwest to southeast (Kraaijpoel et al., 2022), and indeed, we also uncover a negative correlation between b-value and reservoir thickness (Fig. S5d). This correlation between b-value and reservoir thickness is physically intuitive. Assuming seismogenic faults "cut through" the height of the reservoir, larger faults would exist where the reservoir is thicker, allowing for larger events and hence a lower b-value. However, the explanation for the higher b-values to the east and west is less clear.
We also find, as observed by Kaveh et al. (2023) and Gulia (2023), b-value increasing and decreasing with time (Fig. 2). This pattern roughly matches that of seismicity rate, which peaked in 2012-2015 (Fig. 1d). This faint agreement suggests that b-value might depend on the stress rate, rather than stress, as seismicity rate correlates with stress rate to the first order.
We do not see a clear variation of b-value with CSC (Fig. 2). However, when we calculate CSC using the Elastic Thin Sheet model (Bourne & Oates, 2017), updated with a monthly resolution, we recover the negative correlation between b-value and CSC observed by Bourne and Oates (2020) (Fig. 2c). This negative correlation might originate from the difference in spatial distribution between earthquakes at low and high CSC.
We do observe b-value increasing with CSC rate (Fig. 3a). From 10 Pa/day to 25 Pa/day, b-value increases from 0.6-0.8 to 0.9-1.1. Above 25 Pa/day, the b-value from the maximum likelihood methods plateaus around those values, whereas the b-value from the b-absolute method steadily declines to ∼0.8 (Fig. 3a).
We test whether the correlation between b-value and CSC rate is robust to the sampling locations within the mechanical model. Following Kaveh et al. (2023), we calculate CSC rate at the edge of each cuboid and at depths of 1, 10, and 50 m above the reservoir (see Table S1). We also test whether our correlation holds when using the Elastic Thin Sheet model of Bourne and Oates (2017) For all cases, we find a positive correlation between b-value and CSC rate (Fig. 3b, S6), underscoring the robustness of our observation.
We assess the statistical significance of the positive correlation between b-value and CSC rate. We divide the catalog into two groups: a "low CSC rate" group below 15.7 Pa/day and a "high CSC rate" group above 15.7 Pa/day. A two-sample, two-tailed Kolmogorov-Smirnov test between these groups yields a K-S statistic of 0.078, corresponding to a p-value of 0.154. The probability that the high and low CSC rate groups originate from the same distribution is thus 15.4 percent. To evaluate the probability of having such dissimilar distributions by chance, if they are characterized by the b-value, we also conduct a jacknife test. We calculate the b-value of each group and compare these values to those calculated from random samples of half the catalog, taken without replacement. A total of 800 samples are taken. Figure 4a plots the distribution of b-values calculated from each sample, alongside b-values calculated from the high and low CSC -5-manuscript submitted to Geophysical Research Letters rate groups, both of which fall on the edges of that distribution. The probability that a randomly-selected sample of earthquakes has a b-value as low as that of the low CSC rate group is 2.9 percent. Likewise, the probability for a b-value as high as that of the high CSC rate group is 1.0 percent. All b-values are calculated using the maximum likelihood method with a M C of 1.2.
We also validate that the b-values of each group reflect their respective MFDs. Figure 4b shows the MFD of the randomly selected samples, as well as the low and high CSC rate groups. These distributions are consistent with the trend in our estimated b-values.
The distribution of the high CSC rate group is steeper than that of the low CSC rate group. Furthermore, the distribution of either group falls at the edges of the distributions of the randomly selected groups. We repeat the same analyses, but for the CSC rates calculated using the Elastic Thin Sheet model (Bourne & Oates, 2017), and obtain similarly convincing results (Fig. S7).
We acknowledge that the b-value might not be the best quantity to characterize the MFD. Bourne and Oates (2020) argue that these distributions are better quantified using a "taper" at higher magnitudes. However, in the case of the groups with high and low CSC, the b-value is a better metric to characterize the MFD. This is evident upon visual inspection of the MFD of each group. Earthquakes at low CSC rate taper from M =∼2.65, while those at high CSC rate taper from M =∼2.95. These tapers alone, however, are insufficient to capture the difference in MFD. The high CSC rate group is characterized by larger b-values, which also reflect the steeper slope of its MFD (Fig. 4c).
We further verify the appropriateness of the b-value using K-S tests. We randomly generate 1000 synthetic catalogs obeying the Gutenberg-Richter distribution (Li et al., 2023), with b-values equal to that of the low and high CSC rate groups. We then conduct two-sample, two-tailed K-S tests between the synthetic and observed earthquakes, to assess the probability that they originate from the same distribution. For both high and low CSC rate groups, the vast majority of resulting p-values exceed 0.10, meaning we cannot discount the possibility that both groups follow a Gutenberg-Richter distribution (Fig. S8).
In any case, the difference in MFD calculated at high and low CSC rate cannot be interpreted as a result of a different detection level. The shape of both distributions roll over similarly at low magnitudes, indicating a similar detection level.
We find a positive correlation between b-value and CSC rate, but no significant correlation between b-value and CSC when seasonal variations are accounted. Our observations therefore seem to contradict the numerous studies which show the influence of stress on the MFD (T. Goebel et al., 2013;Scholz, 1968Scholz, , 2015;;Rivière et al., 2018;Tan et al., 2019;Dublanchet, 2022;Ito & Kaneko, 2023;Bourne & Oates, 2020;Gulia, 2023;Muntendam-Bos & Grobbe, 2022;Kraaijpoel et al., 2022;Mukuhira et al., 2024). The effect of stress on earthquake magnitude might not be visible in our study because the influence of CSC rate is dominant. There is experimental evidence that stress rate can influence nucleation size, which in turn affects magnitude. Guérin-Marthe et al. (2019) showed that the critical nucleation length of laboratory earthquakes decreases approximately with the logarithm of shear stress rate. Numerical experiments by Dublanchet (2020) show that decreasing nucleation length likely inhibits earthquake growth, resulting in larger b-values. Assum--6-manuscript submitted to Geophysical Research Letters ing that increasing shear stress rate correlates with increasing CSC rate, our results may indicate a decrease of nucleation length with increasing shear stress rate. We caution, however, that this is merely speculation and that experiments involving both normal and shear stress rate are necessary before definitive conclusions can be drawn.
We demonstrate that the MFD of induced earthquakes at Groningen varies in space, possibly in relation to the reservoir thickness. Because the b-value seems to increase with the CSC rate and plateaus, a conservative hypothesis for seismic hazard assessment would be to assume a b-value at the lower end of the distribution (0.6-0.8). The prospective forecast for the current plan to shut down production at Groningen by 2025 would, then, not be very different from the forecast obtained by Bourne and Oates (2020) based on the hypothesis of a correlation between annually averaged stress and the MFD. Since the CSC will remain high in the future while the CSC rate will decrease (because of the decision to shut down production), the two hypotheses would yield similar forecasts. They would, however, diverge if variations in CSC and CSC rates were not anti-correlated. In any case, the dependence of b-value on the CSC rate could be implemented in a stressbased seismicity forecast like Kaveh et al. (2023).
We note that the b-absolute and b-positive values are smaller than the values obtained with the other methods. This is probably because, to maximize the amount of data used, we assumed that magnitude differences as small as (∼0.01) can be resolved. -11-manuscript submitted to Geophysical Research Letters Figure 2. (a) Time series of the number of earthquakes per year. (b) Variation of b-value with time. Variation of b-value with CSC using the (c) reference model and (d) Elastic Thin Sheet model. b-values are calculated using maximum likelihood (blue), modified maximum likelihood with varying MC (pink), and b-absolute methods (orange). Shaded regions represent the 16th through 84th percentiles of each bin. Note that the b-absolute method yields b-values at extreme values of time and stress, as the lack of limitation by catalog incompleteness allows the use of more data.