Original Article

Confirming b-value spatial variations in Groningen gas field using machine-learning-derived seismic data

Laura Gulia^1,2* and Stefan Wiemer¹

¹Swiss Seismological Service, ETH, Zurich, Switzerland; ²University of Pisa, Pisa, Italy

Abstract

The Groningen gas field is Europe’s largest onshore reservoir. Continuous production between 1963 and 2024 has driven reservoir compaction and subsidence, progressively loading pre-existing faults and inducing earthquakes. Recently examined spatio-temporal variations of the Gutenberg–Richter b-value, a proxy for differential stress, across the field, mapping persistent low-b regions in the north and documenting a decade-long trend towards spatial homogeneity closely related to production.

This study revisits Gulia’s results by integrating the original Royal Netherlands Meteorological Institute (RNMI) dataset with the enhanced earthquake catalogue of which newly identifies 660 events and reveals non-clustered and swarm behaviours amongst induced seismicity. The combined/improved catalogue allows to better capture recent trends coinciding with reduced extraction rates and a shift in production towards the south. The aim of this work is to evaluate whether the patterns and hypotheses proposed remain consistent considering the enhanced dataset and to update the analysis with the most recent data. We also evaluate the impact on the evolution of the b-value in the probability of larger events, finding that the continued decrease in b-value compensates for the decrease in seismic activity.

Keywords: Groningen gas-field; induced seismicity; b-value; hazard

 

 

Introduction

The Groningen gas field, located in the northern Netherlands, was discovered in 1959, and production started in 1963. Over subsequent decades, Groningen has been the backbone of the Dutch gas supply system, playing a unique role as a balancing producer capable of modulating production to balance seasonal demand in the Netherlands and much of north-western Europe. Approximately 1,500 faults have been identified to date through comprehensive 3D seismic surveys (NAM, 2016), and the analysis of focal mechanisms indicates that normal faulting is the dominant style of faulting of the field (Willacy et al., 2019). Production caused gradual pressure reduction, leading to an increase in effective vertical stress, progressively reactivating favourably oriented faults: Seismicity is, indeed, associated with the reactivation of northwest–southeast-trending normal faults, at the reservoir level (Bourne et al., 2014; Kortekaas & Jaarsma, 2017; Zbinden et al., 2017).

The region was considered aseismic till 1991, when the first reported earthquake (Ml 2.3) occurred, and following the first events, the rate of seismicity increased non-linearly. Up to 2003, seismic activity remained limited and was primarily concentrated in the central sector of the field, corresponding to the main extraction zone. From 2003 onwards, driven by growing market demand, production levels rose, accompanied by an increase in both the frequency and magnitude of earthquakes. This culminated in a Ml 3.6 event in 2012, which produced significant non-structural damage, increased public concern and drove substantial regulatory intervention (Dost & Kraaijpoel, 2013). Starting from 2014, production measures aimed at lowering the level of seismicity have been implemented, and the production was gradually reduced by 80% in the central part of the field and instead moved towards the south. Public concern prompted extensive research, regulatory intervention and progressively stricter production limits. In 2019, the Dutch government announced an accelerated phase-out of extraction, with production ceasing entirely on 1 October 2023, and its permanent termination legislated in April 2024 (Rijksoverheid, 2024).

One of the fundamental empirical laws of seismology is the Gutenberg and Richter relationship (Gutenberg & Richter, 1944; Ishimoto & Ida, 1939): it estimates the number of earthquakes N larger than or equal to magnitude M, via the formula log(N) = a − bM, whereby the a-value is a volume productivity measure, and the b-value quantifies the frequency-magnitude distribution slope. It has been widely recognised that the b-value is inversely correlated with differential stress on faults: lower values tend to occur in regions of high stress concentration, whereas higher values are generally associated with lower stress levels (e.g. Gulia et al., 2018; Scholz, 1968; Tormann et al., 2015; Wyss, 1973), and b-value is nowadays accepted as a key tool for not only characterising seismic hazard but also assessing real-time stress changes in seismically active areas (Gulia & Wiemer, 2019; Gulia et al., 2024). But seismic catalogues are often incomplete, particularly for smaller magnitude events, due to limitations in detection thresholds or network coverage (Woessner & Wiemer, 2005). Such incompleteness can bias estimates of the b-value. To mitigate biases derived from incompleteness, van der Elst (2021) introduced the b-positive method, which focuses on the positive magnitude differences between successive earthquakes. This approach enhances robustness against catalogue incompleteness and detection thresholds. Subsequent studies have validated its effectiveness: Tinti and Gasperini (2024) confirmed its reliability in estimating b-values from binned data, while Gulia et al. (2024) demonstrated its better performance over traditional methods in real-time seismic monitoring.

A related approach, b-more positive (Lippiello & Petrillo, 2024), further refines the method to account for specific statistical biases, but given the limited volume of seismicity in the present analysis, its use is not necessary here. Lombardi et al. (2025) conduct a systematic seismological experiment to quantitatively assess how much b-value and completeness-magnitude (Mc) estimates vary when different, widely adopted seismological methods are applied to the same earthquake catalogue. During periods of intense seismicity (e.g. aftershock sequences), the discrepancies average approximately 0.4 for Mc and 0.3 for b-value. Other source of inaccuracies can lead to misleading conclusions about seismic hazard and stress distribution: Gulia et al. (2022) demonstrated this in their comment on Taroni et al. (2021), highlighting that the observed b-value patterns in Italy were significantly influenced by quarry blasts, leading to biased interpretations of seismicity, stress variations and incorrect correlations to heat flow. Therefore, before proceeding with any interpretation, it is essential to make an accurate estimate of b-value: b-positive allows to successfully overcome the problems related to incompleteness, while statistical techniques can help identify and remove non-tectonic events (Wiemer & Baer, 2000).

In the 2023 study of the evolution of the Groningen gas field in terms of b-value, Gulia hypothesised that the evolution of seismicity in the field would exhibit systematic spatial and temporal variations in that parameter. Considering that differential stresses should change significantly with production, the following patterns are to be expected:

1. Low b-values would be observed in regions of highest differential stress, corresponding to the central area of the field where production was concentrated until 2014 and where the largest recorded event (Huizinge, Ml 3.6 earthquake) occurred.
2. High b-values would appear in the Southern and on the fringe parts of the field, where stress concentrations were lower.
3. These spatial variations were expected to evolve over time, reflecting both gradual stress loading and changes in production patterns. Specifically, after 2014, when production in the central field decreased and shifted southward, Gulia hypothesises that b-values in the central area would stabilise.

The results of Gulia (2023) largely confirmed these hypotheses, providing a comprehensive, production-aware assessment of Groningen b-values in space and time. Analysis of seismic data from the Groningen field revealed, indeed, a clear spatial variation of b-values, with persistent lower b-values in the north, particularly near the Huizinge epicentral area, contrasting with higher values in the south; a temporal evolution consistent with production changes: the central field showed a stabilisation of b-values after 2014, while regions surrounding the central area experienced a gradual decrease in b-values, reflecting the southward shift in production. Briefly, b-value becomes more spatially uniform over the past decade. These observations support the link between differential stress changes induced by production and the seismicity characteristics, demonstrating that both the activity rate and b-value are sensitive indicators of evolving stress conditions in a producing gas field. Overall, the 2023 analysis provides strong evidence that production-driven stress changes govern both the magnitude distribution and spatial pattern of seismicity in Groningen, aligning with the initial expectations, offering insights for hazard mitigation strategies in extraction-induced seismic contexts.

In 2024, Sirorattanakul et al. created an enhanced earthquake catalogue for the Groningen gas field covering the period 2015–2022, generated using a deep-learning-based workflow. This catalogue contains 1369 events, 660 of which are new compared to the KNMI catalogue, greatly improving completeness, especially for smaller magnitudes. It also characterises seismic clustering, identifying both swarm-like and aftershock-like patterns.

Building upon this advance, the present study leverages the deep-learning-enhanced seismicity catalogue from Sirorattanakul et al. (2024) and the b-positive method to validate the b-value patterns previously reported by Gulia (2023), particularly the spatial heterogeneity, testing the persistence of the original spatial patterns, temporal trends and hypothesis. However, rather than representing a simple update of Gulia (2023), this work, based on enhanced data, seeks to reinforce the use of b-values as a robust, operational parameter for characterising the seismic hazard assessment in a production-related setting.

Data

The dataset consists of three components (Figure 1):

1. The pre-2015 catalogue compiled and used in Gulia (2023), derived from the Royal Netherlands Meteorological Institute – KNMI – detections and locations with local magnitude, Ml, consistent with the configuration in Gulia (2023).
2. The 2015–2022 enhanced catalogue of Sirorattanakul et al. (2024; hereinafter referred to as SIR24), which spans the period January 2015–26 December 2022 and comprises 1,369 located earthquakes detected via an advanced deep-learning workflow. This dataset represents almost double the number of events recorded in the standard KNMI catalogue, with 660 newly identified events previously absent from the KNMI’s record – thereby significantly reducing the estimated Mc (Sirorattanakul et al., 2024). A nearest-neighbour distance analysis reveals that 72% of the recorded events are independent background events, whilst 28% form clusters; of these clustered events, roughly 55% are swarm-like, with the remainder being aftershock-like. A significant portion of new detections are the five bursts of small-magnitude (Ml 0.5–1.5) swarm-like sequences that double the earthquake rates between November 2016 and May 2017, around the time of the maximum compaction rate, suggesting that they might be related to the large strain rate. These swarms propagate along trajectories that do not align with mapped faults or known structures; notably, the swarms cluster temporally around the period of maximum reservoir compaction, estimated between November 2016 and May 2017, within the Zechstein layer above the anhydrite caprock, distinctly separated from the directly induced seismicity within the reservoir (Sirorattanakul et al., 2024).
3. The post enhanced catalogue, events from 28 December 2022 to 29 October 2025, derived from the Royal Netherlands Meteorological Institute – KNMI – detections and locations with local magnitude, Ml.

Figure 1. Magnitude–time (M–t) plot (top) for the composite Groningen catalogue and key milestones of the field history (bottom). Black circles: KNMI catalogue events; red circles: events from the deep-learning-enhanced catalogue of Sirorattanakul et al. (2024).

The total number of events in the composed catalogue is 2,298, and coordinates have been converted from latitude and longitude to the standard system for The Netherlands – Rijksdriehoeksmeting (RD), expressed in terms of easting and northing coordinates (meters), for consistency with the KNMI dataset.

Magnitude conversion

Local magnitudes in the SIR24 catalogue are computed using the same procedure adopted by KNMI (NORSAR, 2018), and the catalogue shows good agreement with KNMI data in terms of event locations. Magnitude differences are generally small: the average difference is less than 0.1 magnitude unit, and only 78 events (12%) exhibit differences larger than 0.2 units (Sirorattanakul et al., 2024).

The construction of a homogenised earthquake catalogue is a fundamental preliminary step. Before merging the two catalogues, it is therefore necessary to convert the SIR24 magnitudes, so they are consistent with the KNMI dataset. Several regression methods are available in literature for this purpose (see Gasperini et al., 2026, for an overview). In general, the regression problem can be framed as the search for the functional model that best reproduces a dependent (response, predicted) variable y as a function of an independent (explanatory, predictor) variable x (Gasperini et al., 2026), and an appropriate regression method should consider experimental errors in both variables. For decades, ordinary least-squares (OLS) regression was the standard method used for magnitude conversions (e.g. Giardini et al., 1997). This approach implicitly assumes that only the dependent variable y is affected by measurement errors, while x is treated as exact.

In the SIR24 catalogue, individual magnitude uncertainties are reported, whereas no corresponding uncertainty estimates are provided in the KNMI catalogue. Generalised orthogonal regression (GOR) was therefore applied under the simplifying assumption that the errors affecting both variables are equal in magnitude and set to unity. This choice ensures a symmetric treatment of the two catalogues and is consistent with the use of GOR as a conservative alternative to OLS when comparable, but not fully constrained, uncertainties affect both axes (e.g. Castellaro et al., 2006). By explicitly accounting for errors in all coordinates and minimising the sum of squared orthogonal distances to the regression line, GOR yields parameter estimates that are less biased and more statistically robust than those obtained with OLS in the presence of measurement errors on both axes.

Figure 2 displays the relationship between KNMI and SIR24 magnitudes, together with the regression lines obtained using OLS and GOR. Both methods return regression coefficients that are very close to the 1:1 line (OLS: ~0.99; GOR: ~1.03), with negligible intercepts, indicating an almost identity transformation between the two magnitude scales. The small deviation of the GOR slope above unity suggests only a marginal tendency for SIR24 magnitudes to be slightly larger than KNMI magnitudes at higher values, but this effect is well within the stated uncertainties. Overall, the regressions confirm an excellent linear correspondence between the two catalogues and do not reveal any substantial systematic bias in the reported magnitudes.

Figure 2. (a) Comparison between local magnitudes (Ml) reported in the KNMI and SIR24 catalogues, with ordinary least squares (OLS) and generalised orthogonal regression (GOR) fits. (b) Density representation of the KNMI–SIR24 magnitude pairs. The near-identity relationship supports combining the catalogues without introducing a significant magnitude scale bias.

We then implemented the Student’s t-test, a parametric procedure used to assess whether the difference between the means of two sets of observations is statistically significant, given the variability within each group. For independent samples, the test compares the sample means under the null hypothesis that the corresponding population means are equal, using a test statistic given by the difference between the sample means divided by an estimate of the standard error of that difference, and evaluates it against a t-distribution with the appropriate degrees of freedom. The method assumes approximate normality of the data, independence of observations and – when the pooled-variance formulation is used – homogeneity of variances.

Applied here to compare the local magnitudes reported in the KNMI and SIR24 catalogues, the t-test yields a p-value of 0.84, indicating that the observed difference between the mean magnitudes of the two catalogues is small relative to the within-group variability and is not statistically significant at the conventional 0.05 level. Consequently, there is no statistical evidence for a systematic bias between the KNMI and SIR24 magnitude scales, and the two catalogues can be reliably combined into a single, homogeneous dataset for subsequent analyses.

Method

b-value mapping

In Gulia (2023), spatial b-value maps were performed using two independent estimators to ensure robustness: the maximum-likelihood formula of Aki (1965), corrected by Utsu (1966) for binned magnitudes, with a fixed Mc, as previously defined by other authors (Dost et al., 2017; Muntendam-Bos, 2020; Paleja et al., 2016) and the b-positive estimator (van der Elst, 2021), which does not require a precise Mc and is robust to temporal changes in catalogue completeness, enabling calculation on datasets with incomplete or variable completeness. The gas field was discretised on a dense 1-km grid, selecting events within a constant 5-km radius for each node, and separate maps were generated for the two key production periods (1991–2013 and 2014–2022). These choices are only apparently arbitrary, as they are grounded in previous observations. In particular, Gulia et al. (2018) demonstrated that spatio-temporal variations in the b-value are closely associated with fault geometry. Consequently, the search radius should be selected in proportion to the dimensions of active faults or to the largest-magnitude event observed within the study area. According to the empirical scaling relationships of Wells and Coppersmith (1994), a normal-faulting earthquake with a magnitude of 3.6 would correspond to a surface rupture length of approximately 0.8 km. However, adopting a radius of this size would result in an insufficient number of events for reliable spatial mapping. In Gulia (2023), a radius of 5 km was identified as the minimum value capable of providing an adequate sample size for robust b-value estimation. It is important to emphasise that, given the 1 km grid spacing and the use of overlapping sampling volumes, the grid-based approach enables the detection of spatial variations whose location and scale are not known a priori. This reduces the dependence on preliminary assumptions regarding the geometry and extent of the underlying processes. Furthermore, overlapping volumes, joined to a dense spaced grid, introduce a degree of continuity between neighbouring observations while simultaneously reducing the likelihood that a localised anomaly falls precisely on a cell boundary, where it could otherwise be artificially fragmented or underestimated. This approach is conceptually analogous to the moving-window technique commonly employed in time-series analysis. Overall, a grid composed of overlapping volumes represents an effective compromise between spatial resolution, statistical robustness and the ability to identify emerging spatial patterns. As an alternative approach to map b-value, Kamer and Himer (2015) proposed to adopt a non-arbitrary partitioning scheme based on a Voronoi tessellation: this method cannot be applied in the Groningen field due to the limited amount of data and small volume (further details in the Supplementary Material).

To update the maps, we applied the b-positive method (van der Elst, 2021), following Gulia (2023): for each node of a 1-km-spaced grid, the b-value was computed in a 5-km radius, only if there were at least 50 positive differences in magnitude within the node’s sample.

Time-series analysis

The north–south comparison through b-value time-series follows Gulia (2023), as well as the sub-fields definition, that is north and south of 593,000. For each sub-field, b-values are estimated using a moving-window approach with one-event steps, assigning each estimate to the time of the last event in the window (Tormann et al., 2013). After selecting the events in each window, we processed the catalogue following the approach of Gulia (2023) but increased the window size to 200 events to exploit the abundant data and quantified the associated uncertainty by means of bootstrapping (Figure 4a–b). To compute the a-value, we considered 1.5 as Mc for the entire period, as this is the highest magnitude threshold for which the seismic network was originally designed (Dost et al., 2017). The rates at this magnitude can therefore be regarded as complete in time. The annual probability (Figure 4d) was subsequently estimated as the inverse of the recurrence time, as inferred directly from the Gutenberg–Richter relationship.

Results

b-value mapping – results

The overall pattern and the correlations with the field’s production history and with a priori expectations, due to changes in gas production, are fully preserved. Gas production declined from ~54 bcm/yr in 2013 (Van Thienen-Visser et al., 2015) to ~9 bcm/yr in 2020 (Boitz et al., 2024), and to ~2.7 bcm/yr in 2024 (https://www.tno.nl), culminating in complete shutdown in 2024. As production decreased, overall seismicity was strongly reduced; consequently, the b-value could be estimated only for a subset of grid nodes. In more detail, spatial maps in the time intervals covered by the Sirorattanakul et al. (2024) dataset confirm the persistent homogeneous b-values in the second interval, almost fully covered by the enhanced catalogue. Here as well as in the analysis by Gulia (2023), the southern sector exhibits a marked decrease, so that the whole field converges to a single, low and uniform value, with values typically in the range 0.8–1.0, without any re-establishment of a pronounced north–south gradient. The initial decline in production in the north, accompanied by a southward shift, tended to uniform stress conditions across the field, compatible with stress redistribution. Compared with the Gulia (2023) map in Figure 3b, the increased number of events now allows us to map almost the entire field (Figure 3d), and the previously inferred homogeneity is thereby confirmed; however, the lowest values are still concentrated in the northern part.

Figure 3. Comparison of spatial b-value maps from (a–b) Gulia (2023) and (c–d) this study, for the same two production periods (see text). Maps in this study are computed with the b-positive estimator on a 1-km grid using a 5-km sampling radius (minimum 50 positive magnitude differences per node). The red circle in panel (d) marks the epicentral location of the 14 November 2025 Ml 3.4 event (shown for reference).

Time-series analysis – results

The overall temporal evolution of the two curves (northern and southern sectors) remains consistent with Gulia (2023). From 2012 to approximately mid-2015, the northern part of the field shows only minor b-value fluctuations around 0.8, whereas higher values – ranging between 1.0 and 1.3 – are observed in the southern sector. Following the production reduction and the associated southward migration, the two time-series converge towards a common value until around mid-2017, when the b-value in the southern region increases up to 1.6. This rise coincides with the occurrence of five swarms present in the Sirorattanakul et al. (2024) catalogue between November 2016 and May 2017, not recorded in the KNMI dataset. Earthquake swarms are defined as groups of clustered events in time and space with similar magnitude, not including a major event (Mogi, 1963). They often show larger b-value when compared with background seismicity (e.g. Wiemer & Wyss, 2002); however, some swarms exhibit a low b-value (e.g. Hainzl & Fischer, 2002; Passarelli et al., 2015). In the Sirorattanakul et al. (2024) catalogue, the swarms are composed of events in the Ml range 0.5–1.5; that is, only the portion of small events is enriched, thus resulting in a higher b-value. Outside this period, no other noticeable swarms were found by the authors. The first four swarms (Figure 9 in Sirorattanakul et al., 2024) occurred in the southern part (as defined in Figure 3), and only the fifth, in 2017, in the northern part: the time series of the northern sector shows an increase only after 2017, in full temporal agreement.

The b-value remains approximately steady, with some fluctuations, in both sectors from 2018 until late 2021, when a decrease is noted in both sectors, due to the occurrence of a M 3.1 event (16 November 2021) in the southern part and to two M 3.1 events (24 September and 8 October 2022) in the northern part of the field. After this date, the maximum measured magnitude is in the dataset (Figure 1). The two curves evolve broadly in parallel. Kaveh et al. (2024) derive a b-value time series using the b-positive approach but applied to the full KNMI earthquake catalogue. Between 2012 and 2018, the b-value time series differs from both Gulia (2023) and the present study, primarily because Kaveh et al. (2024) estimated the b-value using the full catalogue. Seismic productivity is higher in the southern sector of the field, as reflected by the a-value time series (Figure 4c), and consequently, when the entire dataset is considered, the inferred b-value is disproportionately influenced by the southern population. After 2018, event rates decline, and the b-value becomes more spatially uniform; accordingly, the three b-value time series (Gulia, 2023; Kaveh et al., 2024; and this study) converge towards comparable values (0.9 around 2022).

Figure 4. (a) b-value time series from Gulia (2023) and (b) updated b-value time series from this study for the northern and southern sub-fields (division at y = 593,000 m), using a moving window (200 events, one-event step) and bootstrap uncertainty. (c) Time series of annual event rates (a-value) computed assuming Mc = 1.5. (d) Corresponding annual probability of an Ml 3.4 event inferred from the Gutenberg–Richter parameters.

Time-dependent earthquake rate and the Ml 3.4 event of November 2025

The overall rate of seismicity in the Groningen area is known to be steadily decreasing due to first the reduction and then the shut in of seismicity, and this is shown also in the a-value that we compute in Figure 4c. The a-value especially in the south decreased by about 1, indicating a decrease with a factor of 10. Assuming a constant b-value, this would translate into a tenfold decrease in the rate of the hazardous larger event of magnitudes above 3; however, the probability of larger events depends strongly on the b-value. The decrease in the b-value that we document here thus leads to a nearly constant chance of a larger event in the very last years, as shown in Figure 4d. This important effect is sometimes undepreciated in discussion on the hazard and risk in the Groningen region.

While we were finalising this manuscript, a Ml 3.4 event occurred on 14 November 2025 near Zeerijp, in the northern sector of the field (red circle in Figure 3d). Figure 4d shows the annual probability of a Ml 3.4 event in the northern and southern parts of the field. Despite the convergence towards a uniform b-value, the northern sector still exhibits a systematically higher annual probability than the southern one, due to the different rates in the two sectors (Figure 4c). The annual probability from our estimates of an event of about M 3.4 in the northern sector is ~0.28; this is broadly comparable in magnitude scale to the probability (~0.3) reported by Kaveh et al. (2024) for exceeding M 3.6 within their uncertainty-quantified forecasting framework.

The occurrence of the 14 November 2025 Ml 3.4 earthquake in the northern sub-field is therefore consistent with these probabilistic expectations and illustrates that b-value homogeneity does not necessarily imply uniform hazard: spatial variations in rates remain critical. This example underscores the need to interpret b-value maps jointly with productivity and to maintain high-resolution monitoring during the final stages of production and throughout the post-closure period.

Discussion and conclusions

The closure of a field in 2024 and the consequent declining in seismic rates do not exclude the potential for stronger event: case studies around the world show that maximum induced earthquake magnitudes sometimes occur after the termination of energy projects (Schultz et al., 2022). While the overall seismicity rate and consequently the seismic hazard are decreasing, the possibility of triggering large earthquakes remains a major concern (Boitz et al., 2024), which requires continued monitoring, essential for public safety.

In this context, the b-value assumes particular importance because it provides information that is complementary to seismicity-rate observations. Whereas earthquake rates describe how frequently seismic events occur, the b-value reflects the relative proportion of small and large earthquakes and, therefore, indirectly captures changes in the underlying stress conditions of the reservoir. Laboratory experiments, tectonic studies and investigations of induced seismicity have repeatedly shown that lower b-values are generally associated with higher differential stress, while higher b-values correspond to more heterogeneous and less critically stressed conditions. As a result, the b-value can be interpreted as a practical ‘stress-meter’, capable of tracking stress-field evolution through the statistical properties of earthquake populations.

Gulia (2023) shows that changes in b-value can trace the evolution of reservoir loading driven by subsidence and compaction: reductions and spatial shifts in production since 2014 modified the stress field, which is reflected in the b-value patterns before and after this operational change. The enhanced, machine-learning-derived earthquake catalogue by Sirorattanakul et al. (2024), integrated in the KNMI dataset, offered the opportunity to revisit and deepen the analysis by Gulia (2023) using a substantially larger number of events and leveraging the b-positive estimator (van der Elst, 2021). The updated analyses largely corroborate the expectations. The updated analyses largely corroborate the previous findings while extending them into the final years of production and the field-closure phase.

The results confirm the trend towards sustained homogeneity following production changes. Only part of the field is mapped, which not only is a consequence of the imposed resolution and data-availability criteria but also reflects the scarcity of events, which prevents stable b-value estimates in several grid nodes. The b-values observed in the time series (Figure 4b) are slightly lower than those obtained from the spatial maps (Figure 3c–d). This discrepancy arises because the maps display b-values only at grid nodes where the conditions required for a robust calculation are fulfilled, whereas the time series aggregates all events within the entire study region. In this sense, the time series captures the overall temporal evolution of b-value across the region, while the maps add complementary spatial detail by highlighting localised variations where the data density is sufficient.

In particular, spatial mapping confirms the trend towards field-wide homogeneity after production reduction and shift. Time-series comparisons between northern and southern sub-areas indicate a post-2014 stabilisation of the parameter, compatible with depletion-driven stress changes. The enhanced catalogue, comprising roughly twice as many post-2015 events as the KNMI catalogue and revealing five previously unrecognised swarms, enabled more extensive analysis of the last 4 years, when production was sharply curtailed: the b-value remains spatially homogeneous. These findings strengthen the interpretation of the b-value as an operational proxy for stress conditions and are particularly significant because they demonstrate that the b-value responds systematically to operationally induced stress changes and continues to provide meaningful information even when seismicity rates decrease. The principal outcome of this study is therefore not simply the confirmation of the patterns reported by Gulia (2023) but the independent validation of the b-value as an operational proxy for stress conditions in the Groningen reservoir. In other words, the b-value emerges as a measurable, continuously updateable indicator of stress redistribution, capable of translating complex processes into an observable seismological parameter.

Finally, the fact that even small changes in b-value impact substantially the probability of larger events is highlighted by our study. Time-dependent hazard and risk calculation in general and also in the Groningen area hinge critically on b-value extrapolation and also on b-value truncation or roll-over (e.g. Kaveh et al., 2024): the probabilistic seismic hazard analysis conducted by KNMI (Dost & Spetzler, 2015; Spetzler & Dost, 2017) incorporated zone-dependent magnitude-frequency parameters. In the 2017 update, the adopted values are 0.9 for both Central North and Central South and 1.1 for the ‘Active Area’, which groups the rest of the field (Figure 2 and Table 1 in the Report by Spetzler & Dost, 2017). Gulia (2023) reports low and homogeneous b-values around ~0.7–0.8 in the field following the shift in production in 2014; consistently, the enhanced deep-learning-based catalogue of Sirorattanakul et al. (2024) yields a relatively low field-wide estimate (b ≈ 0.86) for 2015–2022. Our study suggests that the hazard and risk in the Groningen area remains higher than one would assume in a model with a constant b-value.

From an operational perspective, routinely updated b-maps and time-series of selected regions can act as complementary, physics-informed indicators alongside rate-based models and geomechanical assessments. Because the parameter is sensitive to evolving stress conditions, it offers an additional layer of information: when combined with modern machine-learning-based earthquake detection, robust event association and precise relocation procedures, b-value monitoring provides a coherent, physics-based framework for tracking reservoir evolution. Such an approach can, in principle, support risk-informed decision-making during all stages of reservoir development and operation; in the Groningen case, this now includes the final stage of field closure and the post-closure period. Analysing b-values through space and time alongside seismicity rates hence offers a practical means of assessing whether stressing conditions are evolving towards an increased or decreased chance of larger ruptures.

Acknowledgments

The authors thank the Editors, the two anonymous reviewers, Emanuele Biondini, Paolo Gasperini, Gianfranco Vannucci and Krittanon Sirorattanakul.

L. G. was partially supported by the Forecasting and understanding the seIsmic cycle through b-value (FARIA) project, funded by the European Research Council Executive Agency under Grant Agreement Number 101200403.

References

Aki, K., 1965. Maximum likelihood estimate of b in the formula logN = a – bM and its confidence limits, Bulletin of the Earthquake Research Institute, The University of Tokyo 43: 237–239.

Boitz, N., Langenbruch, C. & Shapiro, S.A., 2024. Production-induced seismicity indicates a low risk of strong earthquakes in the Groningen gas field. Nature Communications 15: 329. DOI: 10.1038/s41467-023-44485-4

Bourne, S.J., Oates, S.J., van Elk, J. & Doornhof, D., 2014. A seismological model for earthquakes induced by fluid extraction from a subsurface reservoir. Journal of Geophysical Research: Solid Earth 119(12): 8991–9015. DOI: 10.1002/2014JB011663

Castellaro, S., Mulargia, F. & Kagan, Y.Y., 2006. Regression problems for magnitudes. Geophysical Journal International 165: 913–930. DOI: 10.1111/j.1365-246X.2006.02955.x

Dost, B. & Kraaijpoel, D., 2013. The August 16, 2012 earthquake near Huizinge (Groningen). KNMI report. De Bilt: Royal Netherlands Meteorological Institute (KNMI), 26 pp.

Dost, B., Ruigrok, E. & Spetzler, J., 2017. Development of probabilistic hazard assessment for the Groningen gas field. Netherlands Journal of Geosciences 96(S1): S235–S245. DOI: 10.1017/njg.2017.20

Dost, B. & Spetzler, J., 2015. Probabilistic seismic hazard analysis for induced earthquakes in Groningen. KNMI report. De Bilt: Royal Netherlands Meteorological Institute (KNMI), 13 pp.

Gasperini, P., Biondini, E. & Lolli, B., 2026. The role of ‘equation error’ in empirical regressions for seismic magnitude conversions. Bulletin of the Seismological Society of America 116(1): 225–237. DOI: 10.1785/0120250007

Giardini, D., Di Donato, M. & Boschi, E., 1997. Calibration of magnitude scales for earthquakes of the Mediterranean. Journal of Seismology 1: 161–180. DOI: 10.1023/A:1009722406153

Gulia, L., 2023. Time–space evolution of the Groningen gas field in terms of b-value: insights and implications for seismic hazard. Seismological Research Letters 94(4): 1807–1820. DOI: 10.1785/0220220396

Gulia, L., Gasperini, P. & Wiemer, S., 2022. Comment on ‘high-definition mapping of the Gutenberg–Richter b-value and its relevance: a case study in Italy’ by M. Taroni, J. Zhuang, and W. Marzocchi. Seismological Research Letters 93(2A): 1089–1094. DOI: 10.1785/0220210190

Gulia, L., Rinaldi, A.P., Tormann, T., Vannucci, G., Enescu, B. & Wiemer, S., 2018. The effect of a mainshock on the size distribution of the aftershocks. Geophysical Research Letters 45: 13277–13287. DOI: 10.1029/2018GL080619

Gulia, L. & Wiemer, S., 2019. Real-time discrimination of earthquake foreshocks and aftershocks. Nature 574: 193–199. DOI: 10.1038/s41586-019-1606-4

Gulia, L., Wiemer, S., Biondini, E., Enescu, B. & Vannucci, G., 2024. Improving the Foreshock Traffic Light Systems for real-time discrimination between foreshocks and aftershocks. Seismological Research Letters 95(6): 3579–3592. DOI: 10.1785/0220240163

Gutenberg, B. & Richter, C.F., 1944. Frequency of earthquakes in California. Bulletin of the Seismological Society of America 34: 185–188. DOI: 10.1785/BSSA0340040185

Hainzl, S. & Fischer, T., 2002. Indications for a successively triggered rupture growth underlying the 2000 earthquake swarm in Vogtland/NW Bohemia. Journal of Geophysical Research: Solid Earth 107(B12): 2338. DOI: 10.1029/2002JB001865

Ishimoto, M. & Iida, K., 1939. Observations of earthquakes registered with the microseismograph constructed recently. Bulletin of Earthquake Research Institute, the University Tokyo 17: 443–478.

Kamer, Y., & Hiemer S., 2015. Data-driven spatial b value estimation with applications to California seismicity: To b or not to b, Journal of Geophysical Research: Solid Earth 120: 5191–5214. DOI: 10.1002/2014JB011510.

Kaveh, H., Batlle, P., Acosta, M., Kulkarni, P., Bourne, S.J. & Avouac, J.-P., 2024. Induced seismicity forecasting with uncertainty quantification: application to the Groningen gas field. Seismological Research Letters 95(2A): 773–790. DOI: 10.1785/0220230179

Kortekaas, M. & Jaarsma, B., 2017. Improved definition of faults in the Groningen field using seismic attributes. Netherlands Journal of Geosciences 96(S1): S71–S83. DOI: 10.1017/njg.2017.24

Lippiello, E. & Petrillo, G., 2024. b-More-incomplete and b-more-positive: insights on a robust estimator of magnitude distribution. Journal of Geophysical Research: Solid Earth 129, e2023JB027849. DOI: 10.1029/2023JB027849

Lombardi, A.M., Faenza, L., Lucente, F.P., Taroni, M., Gulia, L., Tramelli, A., Biondini, E., Console, R., Convertito, V., Godano, C., De Gori, P., Falcone, G., Garcia, A., Gasperini, P., Lolli, B., Monuori, C., Murru, M., Sandri, L., Spassiani, I. & Vannucci, G., 2025. An experiment on earthquake size distribution estimations reveals unexpected large epistemic uncertainty across methods. Geophysical Journal International 242(2): ggaf193. DOI: 10.1093/gji/ggaf193

Mogi, K., 1963. Some discussions on aftershocks, foreshocks and earthquake swarms: the fracture of a semi-infinite body caused by an inner stress origin and its relation to the earthquake phenomena (Third paper). Bulletin of the Earthquake Research Institute, University of Tokyo 41: 615–658.

Muntendam-Bos, A.G., 2020. Clustering characteristics of gas-extraction induced seismicity in the Groningen gas field. Geophysical Journal International 221(2): 879–892. DOI: 10.1093/gji/ggaa038

NAM, 2016. Winningsplan Groningen Gasveld 2016, Nederlandse Aardolie Maatschappij, BV (Assen), Accessed September 2022. www.nam.nl/feiten-en-cijfers/gaswinning.html (in Dutch)

NORSAR, 2018. Review of the public KNMI induced earthquake catalogue from the Groningen gas field, Project Phase 1, WP1. Accessed March 2024. https://kemprogramma.nl/file/download/18f4a605-4bbc-401a-b68e8bd840b3d05b/1562832346kem11%20norsasodm_groningenreview_wp1.pdf

Paleja, R., Bierman, S. & Jones, M., 2016. Impact of production shut-in on interevent time in Groningen: a statistical perspective. Shell report. The Hague: Royal Dutch Shell, 35 pp. Accessed September 2022. www.nam.nl/feiten-en-cijfers (in Dutch)

Passarelli, L., Hainzl, S., Cesca, S., Maccaferri, F., Mucciarelli, M., Roessler, D., Corbi, F., Dahm, T. & Rivalta, E., 2015. Aseismic transient driving the swarm-like seismic sequence in the Pollino range, southern Italy. Geophysical Journal International 201(3): 1553–1567. DOI: 10.1093/gji/ggv111

Rijksoverheid, 2024. Besluit van 25 april 2024 tot vaststelling van het tijdstip van inwerkingtreding van de Wet van 17 april 2024, houdende wijziging van de Gaswet en Mijnbouwwet in verband met de beëindiging van de gaswinning uit het Groningenveld. Staatsblad 2024: 114. Accessed 4 January 2026. https://zoek.officielebekendmakingen.nl/stb-2024-114.html

Scholz, C.H., 1968. The frequency–magnitude relation of microfracturing in rock and its relation to earthquakes. Bulletin of the Seismological Society of America 58: 399–415. DOI: 10.1785/BSSA0580010399

Schultz, R., Ellsworth, W.L. & Beroza, G.C., 2022. Statistical bounds on how induced seismicity stops. Scientific Reports 12: 1184. DOI: 10.1038/s41598-022-05216-9

Sirorattanakul, K., Wilding, J.D., Acosta, M., Li, Y., Ross, Z.E., Bourne, S.J., van Elk, J. & Avouac, J.-P., 2024. Bursts of fast propagating swarms of induced earthquakes at the Groningen gas field. Seismological Research Letters 95(1): 130–146. DOI: 10.1785/0220240107

Spetzler, J. & Dost, B., 2017. Probabilistic seismic hazard analysis for induced earthquakes in Groningen. KNMI report. De Bilt: Royal Netherlands Meteorological Institute (KNMI).

Taroni, M., Zhuang, J. & Marzocchi, W., 2021. High-definition mapping of the Gutenberg–Richter b-value and its relevance: a case study in Italy. Seismological Research Letters 92: 3778–3784. DOI: 10.1785/0220210017

Tinti, S. & Gasperini, P., 2024. The estimation of b-value of the frequency–magnitude distribution and of its 1σ intervals from binned magnitude data. Geophysical Journal International 238(1): 433–458. DOI: 10.1093/gji/ggad091

Utsu, T., 1966. A statistical significance test of the difference in bvalue between two earthquake groups, Journal of Physics of the Earth 14: 34–40.

Van Thienen-Visser, K., Sijacic, D., Nepveu, M., van Wees, J.-D. & Hettelaar, J., 2015. Response of induced seismicity to production changes in the Groningen field. TNO Report TNO 2015 R11367: 56 pp.

Tormann, T., Enescu, B., Woessner, J. & Wiemer, S., 2015. Randomness of megathrust earthquakes implied by rapid stress recovery after the Japan earthquake. Nature Geoscience 8(2): 152–158. DOI: 10.1038/ngeo2343

Tormann, T., Wiemer, S., Metzger, S., Michael, A. & Hardebeck, J.L., 2013. Size distribution of Parkfield’s microearthquakes reflects changes in surface creep rate. Geophysical Journal International 193(3): 1474–1478. DOI: 10.1093/gji/ggt093

Utsu, T., 1966. A statistical significance test of the difference in bvalue between two earthquake groups, Journal of Physics of the Earth 14: 34–40.

van der Elst, N.J., 2021. B-positive: a robust estimator of aftershock magnitude distribution in transiently incomplete catalogs. Journal of Geophysical Research: Solid Earth 126(2): e2020JB021027. DOI: 10.1029/2020JB021027

Wells, D.L. & Coppersmith, K.J., 1994. New empirical relationships among magnitude, rupture length, rupture width, ruptur area, and surface displacement, Bulletin of the Seismological Society of America 84: 974–1002.

Wiemer, S. & Baer, M., 2000. Mapping and removing quarry blast events from seismicity catalogs. Bulletin of the Seismological Society of America 90(2): 525–530. DOI: 10.1785/0119990104

Wiemer, S. & Wyss, M., 2002. Mapping spatial variability of the frequency-magnitude distribution of earthquakes. Advances in Geophysics 45: 259–302. DOI: 10.1016/S0065-2687(02)80007-3

Willacy, C., Dedemvan, E., Minisini, S., Li, J., Blockland, J.W., Das, I. & Droujinine, A., 2019. Full-waveform event location and moment tensor inversion for induced seismicity. Geophysics 84(2): KS39–KS47. DOI: 10.1190/GEO2018-0212.1

Woessner, J. & Wiemer, S., 2005. Assessing the quality of earthquake catalogues and the magnitude of completeness. Bulletin of the Seismological Society of America 95: 684–698. DOI: 10.1785/0120040007

Wyss, M., 1973. Towards a Physical Understanding of the Earthquake Frequency Distribution, Geophysical Journal International 31(4): 341–359. DOI: 10.1111/j.1365-246X.1973.tb06506.x

Zbinden, D., Rinaldi, A.P., Urpi, L. & Wiemer, S., 2017. On the physics-based processes behind production-induced seismicity in natural gas fields. Journal of Geophysical Research: Solid Earth 122: 3792–3812. DOI: 10.1002/2017JB014003