Production-induced seismicity indicates a low risk of strong earthquakes in the Groningen gas field

It has been known for decades that both subsurface fluid extractions and injections can cause earthquakes. The question of the maximum possible earthquake magnitude, Mmax, in connection with such geo-technical operations is important for understanding, evaluating and controlling their seismic hazard. This question is under-researched and controversial, especially when it comes to the long-term production of hydrocarbons such as in the Groningen gas field. The Groningen gas field is the largest in western Europe with a total gas volume of approximately 2900 bcm located in the Upper Rotliegend at approximately 3 km depth. The reservoir thickness increases from about 150 m in the south-east to close to 300 m in the north-west. The field was discovered in 1959; production started in 1965, peaked at more than 80 bcm/year in the 1970s, and then declined to 30–50 bcm/year. After the largest earthquake (Mw 3.5, ML 3.6) in 2012, lower production rates were mandated. Pore pressure decreased up to 1 MPa/year with a total depletion of approximately 28 MPa until 2022. The end of gas production is planned, with partial backup operation and full closure planned for 2030. Seismicity caused by gas production is the reason for closure. The first earthquake was recorded in 1991, followed by a relatively constant annual number (around 10–20 above ML = 1.2) until 2002, after which the rate increased. Until February 2022, 1474 earthquakes with magnitudes between ML = −0.8 and ML = 3.6 were registered. Because seismicity may continue after extraction stops and maximum induced magnitudes can occur post-operations, the maximum possible magnitude remains intensively discussed. To understand the nature of a maximum possible earthquake, a distinction must be made between events predominantly caused by geotechnical interventions and tectonically prepared earthquakes triggered by anthropogenic influences. Different definitions exist; here, the authors use the definition that induced earthquakes rupture entirely within the stimulated rock volume, while triggered earthquakes have ruptures that run outside the stimulated volume. This approach relates directly to observable features and is close to the concept of a runaway earthquake. Induced seismicity is mainly controlled by technical operations and is easier to manage (e.g., via traffic-light systems) than triggered seismicity. Understanding both the maximum induced and maximum triggered earthquake is vital for hazard assessment, particularly at Groningen. The study investigates whether induced seismicity can indicate the possibility of triggering large tectonic earthquakes. The authors develop an approach to estimate the maximum possible induced magnitude MI for fluid extraction from a flat underground reservoir, applying the Lower-Bound (LB) frequency-magnitude statistic to Groningen. They then assess the worst-case probability of triggering larger tectonic earthquakes using the Seismogenic Index (SI) framework extended to production.

The paper reviews differing definitions of induced versus triggered earthquakes: (1) comparison of earthquake shear stress release to the stress impact of stimulation; (2) events that only occur due to geotechnical operations and are fully controlled by man-made stress changes versus triggered events on faults favorably oriented to tectonic stresses; and (3) a geometric definition where induced earthquakes rupture entirely within the stimulated rock volume, while triggered earthquakes rupture beyond it. The study adopts the third, observation-linked definition, aligned with the concept of runaway earthquakes. The authors discuss prior difficulties in constraining upper magnitude bounds or tapers of the Gutenberg-Richter (GR) statistic for Groningen using purely statistical methods, citing studies that highlight controversy and uncertainty regarding Mmax in long-term hydrocarbon production settings. They note that induced seismicity often deviates from GR behavior, with underrepresentation of large magnitudes and variable b-values, motivating physically based alternatives like the LB-statistic that incorporate finite stimulated volumes. They also build on prior work developing the Seismogenic Index (SI) to relate seismic response to stress changes in injection/production operations and recent formulations of worst-case exceedance probabilities (WCEP) for triggered events.

Assumptions and theoretical estimates: - Initially assume all Groningen earthquakes are induced and ruptures are limited to the production-stimulated reservoir layer. Ruptures are approximately isometric and occur on normal-faulting planes dipping ~60°. - Compute maximum reservoir-limited magnitudes using rupture geometry and stress drop: Mw = log10(XY) + log10(Δσ)/1.5 – 6.07, where XY is effective rupture length; largest rupture length XY = h/sin(π/3), with h the local reservoir thickness. Stress drops of 1 and 10 MPa bound observed values in Groningen. Maximum rupture lengths are ~200 m (SE) to ~500 m (NW), giving local Mmax ranges of 2.6–3.3 (SE) and 3.5–4.1 (NW). Lower-Bound (LB) frequency-magnitude statistic for induced seismicity: - Model the stimulated volume as a thin, laterally extensive reservoir layer; earthquakes occur on critically oriented faults; rupture length is approximately proportional to rupture height and not exceeding it significantly. - Derive the LB FM-distribution by integrating the GR probability density over the probability that ruptures fall entirely within the stimulated volume, yielding an upper magnitude limit MY (maximum possible induced magnitude). Exact and approximate forms: log10 NM = a – bM + log10(1 + 2 b⁻1 10^{-b(M–MY)} – 2 b⁻1 10^{-2b(M–MY)}). An accurate approximation decoupling b and MY is: log10 NM = a – bM + 2 log10(1 – 10^{-b(M–MY)}). - Fit LB and classical GR models to Groningen’s catalog (nonlinear fits), estimate parameters a, b, and MY, and compare models using Akaike Information Criterion (AIC). Also estimate temporal evolution of parameters by moving-window fits. Temporal and spatial analyses: - Temporal fits of GR and LB to estimate time evolution of a, b, and MY, with uncertainties from covariance matrices and additional statistical assessment. - Spatial mapping: For each map node, select events within 15 km radius and fit the LB approximation to obtain local a, b, and MY; report MY only where its uncertainty is ≤0.2 to ensure robustness. Seismogenic Index (SI) estimation and scenarios: - Estimate time-dependent SI using Σo(t) = log10(Nt^c) + Mc b(t) – δΣ(t), where Nt is number of events above completeness Mc. Two scenarios account for catalog completeness and potential triggering thresholds: (1) seismicity controlled by pore-pressure changes since start of production; (2) seismicity controlled by pore-pressure changes only after first observed earthquake in 1992. This brackets the SI range. Worst-Case Exceedance Probability (WCEP) for triggering large tectonic earthquakes: - Extend prior WCEP formulation (developed for injections) to fluid production. Worst case uses upper bound of SI and lower bound of b: W_M^{Mmax}(t) = 1 – exp[ – 10^{δΣ(t)} 10^{sup Σo(t) – inf b(t) Mmax} ]. - Compute δΣ(t) from poroelastic stress changes tied to pore-pressure depletion: δFCS(t) ∝ δPp(t). For a thin-layer reservoir approximation: δΣ(t) = log10[A H S (1–ns)/(ns sin φ)] + log10[δPp(t)], with ns = 0.375 and S = 5×10^-10 as Groningen parameters; A is reservoir area, H thickness, φ friction angle. Use observed pore-pressure change time series and fitted a, b to evaluate WCEP for target magnitudes. Data: - Use published Groningen production volumes, pore pressure, and earthquake catalog through February 2022; production peaked >80 bcm/yr in 1970s; cumulative depletion ~28 MPa by 2022; catalog includes 1474 events (ML −0.8 to 3.6).

- Maximum possible induced magnitude: LB-statistic applied to the complete catalog yields an upper bound MY ≈ Mw 3.97 (~4). Theoretical geometry with stress drops 1–10 MPa gives local induced Mmax of ~2.6–3.3 (SE) and ~3.5–4.1 (NW), consistent with observed spatial trends and maximum observed magnitude ML 3.6 (Mw 3.5). - Model comparison: Classical GR fit (b = 0.94, a = 3.96) statistically expects Mmax = a + b/2 ≈ 4.21 and implies only a 2.34% chance of observing no event > M 3.6 in the complete catalog, contradicting observations. LB provides a superior fit (lower AIC) and explains the deficit of large events via a physically based upper bound. - Temporal behavior: GR b-values are relatively stable near ~0.95; LB b-values start near ~0.85, fluctuate early due to limited statistics, and stabilize. Temporal changes in b and MY are mainly statistical; estimates from the full catalog (GR b ≈ 0.94; LB b ≈ 0.76; MY ≈ 3.97) are preferred. - Spatial segmentation: NW part exhibits lower b (~0.8) and higher MY (~4.1); SE part shows higher b (>1) and lower MY (~3.7). a-values are relatively uniform (central high, decreasing to rim). The Seismogenic Index shows no significant spatial variation. - Seismogenic Index (SI): Considering two scenarios (since start of production vs. since first observed event in 1992) constrains SI within a range; scenario 1 shows SI increase from 2000–2012 then leveling; scenario 2 remains nearly constant, consistent with stationary conditions and possible gradual involvement of cohesive faults. - Worst-Case Exceedance Probability (WCEP): Using upper SI and lower b bounds, the worst-case probability to trigger tectonic events up to Mw 5.5 becomes significant (up to ~30%). If triggering larger tectonic earthquakes were feasible under Groningen conditions, events of Mw ≥ 4 would be expected to have occurred several times given the long production history; they have not. - Implication: The agreement of observed maximum magnitudes with reservoir-limited predictions and the absence of larger events despite long-term stressing indicate that Groningen seismicity is predominantly induced and reservoir-limited. This suggests an inherently stable field where physical conditions for triggering large tectonic earthquakes likely do not exist, implying a low risk of strong earthquakes.

The study addresses the central question of the maximum possible earthquake magnitude associated with Groningen gas production by integrating physical constraints from reservoir geometry with statistics that account for finite stimulated volumes. The LB-statistic captures the observed deficit of large events and provides a physically meaningful upper magnitude bound near Mw 4, unlike the classical GR model which overpredicts the likelihood of larger events. The spatial and temporal analyses reinforce the interpretation that observed earthquakes are induced and limited by the reservoir. Spatial correlation between reservoir thickness and both observed maxima and LB-derived MY supports the reservoir-limited rupture model. Temporal variability in b and MY appears largely due to statistical fluctuations, with stable estimates from the full catalog providing reliable parameters for hazard assessment. By extending the Seismogenic Index framework and WCEP to production, the authors evaluate the worst-case probability of triggering larger tectonic earthquakes in the surrounding continuum. Given Groningen’s long production history, the worst-case analysis implies that if conditions for triggering larger tectonic events existed, several Mw ≥ 4 events would likely have occurred, which contradicts observations. This discrepancy suggests that the necessary physical conditions for runaway ruptures extending beyond the reservoir volume are absent, and the field is inherently stable. These findings support the use of induced-seismicity-based models to inform operational controls and hazard assessments and indicate a low likelihood of strong earthquakes in Groningen, even considering post-production seismicity.

The paper introduces a physically grounded framework for assessing maximum induced magnitudes and worst-case triggering probabilities in a long-produced gas field. By applying the Lower-Bound frequency-magnitude statistic with reservoir geometry constraints, the maximum induced magnitude in Groningen is estimated at approximately Mw 4, consistent with observed magnitudes and spatial trends tied to reservoir thickness. Extending the Seismogenic Index and Worst-Case Exceedance Probability formulations to fluid production shows that, under worst-case assumptions, significant probabilities for Mw up to 5.5 would have manifested as multiple strong events over Groningen’s long production history if triggering were feasible. Their absence indicates that Groningen is inherently stable and that conditions for triggering large tectonic earthquakes likely do not exist, implying a low risk of strong earthquakes. The approach provides a transferable methodology for other geo-energy projects to distinguish induced from triggered seismicity and to bound hazard using physically motivated statistics and stress-change metrics. Continued monitoring and application of LB and SI/WCEP analyses with improved catalogs can further refine hazard estimates.

- Catalog completeness: Early observations (pre-1992) are likely incomplete due to a sparse network; SI is therefore evaluated under two scenarios to bracket uncertainties. - Statistical uncertainties: Covariance-based uncertainties may underestimate true errors; early temporal variations of parameters (b, MY) are dominated by limited statistics. LB parameter coupling (b and MY) in the exact formulation can lead to inaccuracies in poor statistics; the approximation mitigates but does not eliminate this. - Model assumptions: Ruptures are assumed approximately isometric, limited to the reservoir layer, on critically oriented faults with ~60° dip. Reservoir is approximated as a thin, laterally extensive, planar layer, neglecting non-planar boundaries and fault complexity. - Stress-drop bounds: Stress drops are assumed in the 1–10 MPa range; local deviations could affect magnitude estimates. - GR and Poisson assumptions: Comparisons rely on GR fits and Poissonian independence for expected maximum magnitude calculations; deviations from these assumptions could alter expectations. - Spatial mapping: Use of a 15 km smoothing radius may obscure finer-scale heterogeneity; MY is mapped only where uncertainty ≤0.2, leaving gaps in data-sparse areas. - WCEP worst-case framing: Triggering probabilities use conservative bounds (upper SI, lower b). Actual probabilities may be lower; WCEP does not model specific fault geometries or dynamic triggering processes. - Parameterization of δΣ: Uses reservoir-scale poroelastic parameters (ns, S) and friction angle; uncertainties in these inputs and pore-pressure change estimates propagate into WCEP and SI.