The exponential rise of induced seismicity with increasing stress levels in the Groningen gas field and its implications for controlling seismic risk

S U M M A R Y Induced seismicity typically arises from the progressive activation of recently inactive geological faults by anthropogenic activity. Faults are mechanically and geometrically heterogeneous, so their extremes of stress and strength govern the initial evolution of induced seismicity. We derive a statistical model of Coulomb stress failures and associated aftershocks within the tail of the distribution of fault stress and strength variations to show initial induced seismicity rates will increase as an exponential function of induced stress. Our model provides operational forecasts consistent with the observed space–time–magnitude distribution of earthquakes induced by gas production from the Groningen field in the Netherlands. These probabilistic forecasts also match the observed changes in seismicity following a significant and sustained decrease in gas production rates designed to reduce seismic hazard and risk. This forecast capability allows reliable assessment of alternative control options to better inform future induced seismic risk management decisions.


I N T RO D U C T I O N
Industrial activities transport significant fluid volumes into or out of geological formations for the purposes of hydrocarbon extraction, waste-water disposal, CO 2 storage and geothermal energy provision. The changed pore fluid pressures inside these formations create poroelastic deformations which induce surface elevation changes (Geertsma 1973;Segall 1989) and occasionally earthquakes (Majer et al. 2007;Suckale 2009;Davies et al. 2013;Ellsworth 2013). Regulation and control of these activities often relies on a probabilistic analysis of the associated future seismicity, ground motion hazards and building damage risks-much like existing best practice for managing natural seismic hazard and risk (Cornell 1968;McGuire 2008). Such probabilistic analysis aims at forecasting probability distributions for the number, location, origin time and magnitude of future earthquakes. Leading forecast models for natural seismicity (e.g. Helmstetter et al. 2006) rely on statistical stationarity so that forecasts are based on historic frequencies. Induced seismicity requires a fundamentally different approach because the operational activities are not stationary in time and may also be adapted to reduce the seismic hazard.
Existing poroelastic theories succeed in describing the occurrence of induced seismicity (Segall 1989;Segall et al. 1994;Zoback & Zinke 2002;Segall & Lu 2015;Shirzaei et al. 2016), but lack sufficient detail to quantitatively match the statistical space-timemagnitude history of observed seismicity or to assess seismic hazard and risk. The missing element seems to be heterogeneities that localize and amplify stress build-up on pre-existing faults. Allowing for small-scale heterogeneities in fault strength and treating these as a uniformly distributed stochastic quantity predicts a constant rate of induced seismicity relative to the incremental stresses found to be in good agreement with observed rates of induced seismicity (Shapiro et al. 2007;Langenbruch & Zoback 2016) for injected fluid volumes up to 10 7 m 3 . For larger induced deformations, an exponential-like rise of induced seismicity rates is witnessed in laboratory experiments (Scholz 1968), volcanic magma accumulations (Lengliné et al. 2008) and the Groningen gas field (Bourne et al. 2014) where the bulk reservoir volume change exceeds 2 × 10 8 m 3 (Fig. 1). We will show this initial exponential-like rise of induced seismicity is a general consequence of fault stress disorder where initial fault failures correspond to extreme stress values in the tail of the frequency distribution.
Our new theory for fault reactivation induced by poroelastic deformations provides an operational, stochastic-physics model to account for the early spatial-temporal evolution of induced seismicity rates and magnitudes as a function of reservoir pore pressure, strain and topographic gradients (Fig. 2) within a heterogeneous, thinsheet reservoir geometry (Fig. 3). Elastic and geometric reservoir The observed relative increase in annual earthquake rates, dN/dt, exceeds the relative changes in annual reservoir volume change, dV/dt. Simulated seismicity from a stochastic model with disordered fault strength matches these observation trends and variability within the 95 per cent prediction interval (grey bands).
heterogeneities govern the incremental Coulomb stresses induced by pore-pressure changes. Under these incremental Coulomb stress loads, the weakest parts of pre-existing geological fault surfaces may experience frictional failure resulting in an induced earthquake. Structural heterogeneities act to localize shear stress development and seismicity within regions of greatest structural gradients, often associated with faults that partially offset the reservoir. Elastic heterogeneities may also localize induced shear stresses and seismicity within regions of greatest reservoir compressibility, often associated with higher reservoir porosities. Fault friction heterogeneities, due to geometric irregularities of the slip surface and frictional slip history localize seismicity within the weakest fault segments. Friction heterogeneity also governs an exponential-like increase in the rate of induced seismicity and an increase in the expected magnitudes as the fraction of fault segments loaded close to failure increases.
The resolvable elastic and geometric heterogeneities are represented by a deterministic, smoothed, poroelastic, stress-strain tensor field. These heterogeneities are estimated using geodetic measurements of surface displacements, geophysical imaging of the reservoir geometry and in-well monitoring of reservoir pore pressures. In contrast, frictional fault strength and initial stress heterogeneities are not directly observable so they are represented by a single invariant probability distribution of initial fault stress and a transient stochastic function for stress changes due to previous earthquakes. This results in a stochastic model for earthquake occurrences and their aftershocks that may be estimated using the observed earthquakes. Application of this theory to a major European hydrocarbon gas field adjacent to the city of Groningen in the Netherlands demonstrates a statistical earthquake history match and correctly forecasts the response to an unprecedented 50 per cent reduction in fluid extraction rates in just 3 yrs.

Poroelastic thin-sheet model
Fault instability corresponds to a maximum Coulomb stress, C, on an optimally oriented fault plane of at least zero, where and C i is the initial Coulomb stress, τ m is the change in maximum shear stress, μ is the friction coefficient and σ m is the change in mean total normal stress. For a heterogeneous, faulted reservoir geometry that is laterally extensive relative to its thickness and depth below the Earth's surface (Fig. 3), poroelastic displacements and strains induced by porepressure changes are predominantly vertical and uniaxial, respectively. Within this poroelastic, thin-sheet model (Bourne & Oates 2017b), in the limit of large geometric heterogeneities the vertically averaged maximum Coulomb stress, C, is where x is the reservoir position vector, t is time, γ = (1 − 2ν)/(2 − 2ν) and ν is Poisson's ratio, zz (x, t) is the vertical strain, (x) is the magnitude of initial topographic gradients of the top reservoir surface and H (x) is a poroelastic modulus defined as 1/H = 1/H r + 1/H s where H r = P/ zz which is the ratio of porepressure change to vertical strain, and H s = 3K s (1 − ν)/(1 + ν) Figure 2. A network of physical processes defines the stochastic model used to forecast induced seismicity. The observed elastic and geometric reservoir heterogeneities constrain a smoothed deterministic model of incremental Coulomb stresses, C, whilst the unobserved frictional heterogeneities, C i , are modeled stochastically. Using Bayesian inference, the hidden quantities (grey circles) are constrained by the resolvable heterogeneities (white circles) as witnessed by the available reservoir monitoring systems. which is a constant dependent on K s and ν, the bulk modulus of the skeleton material and the Poisson's ratio, respectively. In the limits of large or small reservoir stiffness, H r , relative to a material constant H s , the quantity H zz equals pressure depletion, P, or vertical strain, zz , respectively. The spatio-temporal evolution of zz (x, t) and P(x, t) may be inferred from geodetic measurements of induced surface displacements, and downhole reservoir pressure gauges, respectively. Reservoir thickness and the top reservoir surface, z, may be mapped from a reflection seismic image. Observed vertical strains depend on the static elastic reservoir properties and the time-varying distribution of reservoir pore-pressure changes. Observed topographic gradients encode information about the static structural heterogeneities that localize shear strains induced by vertical strain. The material constant, H s , determines the extent to which Coulomb stress changes correlate with reservoir strains (H s H r ) or pore-pressure changes (H s H r ). Due to uncertainties in measured epicentres and reservoir strains, and the limited number of observed earthquakes, the covariate, is smoothed with a Gaussian kernel of length-scale, L s .

Extreme threshold failure model
Fault strength heterogeneity is not mappable and so is treated as a stochastic property by simply representing the initial Coulomb stress on fault segments, C i , as an independently and identically distributed random variable. Under increasing load, initial fault failures occur for the largest values of C i within the tail of the distribution. The probability of fault failure, P f , under the condition C ≥ 0, is equivalent to the probability of the initial Coulomb stress exceeding a reservoir deformation-dependent threshold, C i ≥ − C. Regardless of the particular distribution of C i , and according to Extreme Threshold Theory (Pickands 1975;Coles 2001) P f follows the Generalized Pareto (GP) distribution. Using the simplest form of the GP distribution, an exponential trend, the corresponding rate of independent failure events with time and per unit area within some infinitesimal region is described by the Poisson point process intensity function, λ m , such that: where θ 0 and θ 1 are scalar parameters describing the location and scale of the GP distribution, h(x) is the reservoir thickness map and˙ C is the temporal derivative of C(x, t). The parameter vector {θ 0 , θ 1 , H s , L s } links the resolved reservoir deformation covariate, C, to the observed activity rate of induced independent earthquakes, λ m . Simulated moment magnitudes are independently and identically distributed according to a truncated exponential distribution (Cornell & Van Marcke 1969) where β describes the decrease in probability with magnitude, and M max is the maximum possible magnitude. Reloading previously inactive faults induces initial seismicity with expected magnitudes that increase with load as observed experimentally (Zang et al. 1998;Ojala 2003;Heap et al. 2009), theoretically for fibre bundle failures (Hemmer & Hansen 1992), and induced by mining (Urbancic et al. 1992;Gibowicz & Lasocki 2001) and reservoir compaction (Bourne et al. 2014). We model variation in β as the inverse power of C: which becomes β = β 0 in the limit C 0 C or n = 0. The parameter vector {β 0 , 0 , n} models the covariance between the observed reservoir deformations, , and earthquake magnitudes, M.

Aftershock model
Aftershocks induced by previous earthquakes as described by the Epidemic Type Aftershock Sequence (ETAS) model (Ogata 2011) are included by extending the intensity function as: where λ m is the intensity function for independent main events, M j is the magnitude of the jth event, and t i , t j , x i , x j are the origin times and locations of the i th and j th events, respectively. The aftershock triggering function, f, is defined as where t is the interevent time, and r is the interevent distance. The ETAS parameter vector {K, a, c, p, d, q} defines the productivity and spatial-temporal clustering of aftershocks triggered by previous main event and aftershocks. The observed vector of origin times, epicentres and magnitudes for n events, X = (X 1 , . . . , X n ), are considered to result from a probability model p(X|θ) with an unknown model parameter vector, θ = {θ 0 , θ 1 , H s , L s , β 0 , 0 , n, K , a, c, p, d, q}. The joint posterior probability of model parameters after incorporating these observed events, p(θ |X), is obtained based on the Markov Chain Monte Carlo method of Bayesian inference using the Metropolis-Hastings algorithm (Patil et al. 2010).

G RO N I N G E N G A S F I E L D
Western Europe's largest natural gas field comprises a 40 by 50 km wide and 250 m thick accumulation within the Rotliegend reservoir at 3 km depth (Stauble & Milius 1970), located in the northeast of the Netherlands, adjacent to the city of Groningen (Fig. 4), and below a present population of at least 500 000 residents and 250 000 buildings. The recoverable volume of Groningen gas is about 2800 billion cubic metres of which about 70 per cent has been produced so far by some 300 wells since 1963, leaving a further 800 billion cubic metres producible over at least the next 40 yrs. Reservoir porepressures gauges measured uniform depletion from 25 to 12 MPa at a mean rate of 0.2 MPa yr −1 since 1995. Surface displacements witnessed by field-wide optical leveling surveys acquired 14 times since 1964 (Hejmanowski 2000), monthly Interferometric Synthetic Aperture radar (InSAR) images since 1995 (Ketelaar 2009) Geophysical reflection seismic surveys image a dense heterogeneous network of pre-existing normal faults that transect and partially offset the reservoir interval (Fig. 4b). The Royal Netherlands Meteorological Institute (KNMI) has operated a local near-surface seismic monitoring network in the northeast of the Netherlands reporting all M L ≥ 1.5 events within the field since 1995  with standard epicentral errors of 500 m, similar to the mapped fault spacing. Epicentral density is heterogeneous and localized primarily in proximity to the largest reservoir shear strains close to the centre of the field where reservoir compaction and fault offsets are also largest. Normal-faulting focal mechanisms consistent with these faults were found for the largest few events , including the largest event, M L = 3.6 in 2012 (Dost & Kraaijpoel 2013). Two downhole monitoring arrays located within and above the reservoir and operated since 2013 find almost all hypocentral depths to be inside the reservoir to within a 50 m depth error (Daniel et al. 2016) indicating these events initiate inside or close to the reservoir interval. Given the reservoir thickness is 100-300 m, larger magnitude events, M > 3, may propagate outside the reservoir into a different stress regime.
Between 1995 and 2014 annual produced gas volumes varied from 23 to 54 × 10 9 m 3 (bcm) according to market demand. Correspondingly, the rates of annual reservoir pressure depletion and volume decrease varied from 0.2 to 0.35 MPa and 3 to 5 × 10 6 m 3 , respectively. Following widespread concern about the increasing trend of seismicity, production rates from the centre of the field were minimized from 2014 January onward and total annual production rates reduced from 53.87 bcm (2013) to 42.41 bcm (2014), 28.10 bcm (2015) and 27.59 bcm (2016). Such large, sustained production rate decreases were unprecedented during seismic monitoring since 1995.

F O R E C A S T P E R F O R M A N C E R E S U LT S
Our model reliably forecasts the 2012-2017 space-time-magnitude density of earthquakes induced in response to this production control measure according to pseudo-prospective tests of ensemble simulations using the joint posterior distribution of parameter values estimated from the period 1995-2012 (Fig. 5).
Due to aftershocks, which represent 10-20 per cent of the 2012-2017 forecast events (Fig. 5c), the relative abundance of shorter interevent times (Fig. 5d) and distances (Fig. 5e) match observations within the 95 per cent prediction intervals. Equivalent Poisson models excluding aftershocks underpredict this clustering and the interyear rate variability. The observed 2012-2017 epicentral density is within the 95 per cent prediction interval (Fig. 5f) despite the minimization of production from the zone of greatest seismicity during this time.
Forecast performance tests (Schorlemmer et al. 2010;Marzocchi et al. 2012) for the event number (N-test, Fig. 6a) and the log likelihood of space-time-magnitude distribution given the forecast model (L-test, Figs 6b and c) are within their 95 per cent prediction intervals. The relative log-likelihood performance of alternative models with less heterogeneity are all significantly worse (R-test, Figs 6d and e).
Relative to the full model, which comprises of activity rates with an exponential strain trend and aftershocks (EST-ETAS) and b-values with an inverse power-law strain trend (IPST), forecast performance reduces without aftershocks (EST), elastic heterogeneity (EPT, IPPT), initial stress and geometric heterogeneity (CT, PT) and without any heterogeneity (UNI). The complete model distributes activity rates with an exponential shear strain trend and aftershocks (EST-ETAS) and magnitudes with an inverse power-law shear strain trend of β values (IPST). The null model has uniform activity rates (UNI) and constant β-value magnitudes (UNI). Excluding aftershocks (EST) degrades forecast activity rates. Excluding elastic heterogeneity (EPT, IPPT) further reduces performance, and without initial stress heterogeneity activity rates (PT, CT) are similar to the null model and magnitudes (IPPT) are worse still. Forecast performance critically depends on elastic and geometric reservoir heterogeneities that localized induced shear stresses, fault heterogeneities that govern initial fault stresses, and the transient local stress transfers that trigger aftershocks.

D I S C U S S I O N A N D C O N C L U S I O N S
Field observations within the Groningen gas field indicate that some 15 MPa mean reservoir pore-pressure depletion and up to 0.35 m compaction are associated with seismicity rates that increased as an exponential-like function of induced stress and b-values that decreased like an inverse power function of induced stress. We have shown this behavior is consistent with a simple physical model that assumes reservoir deformations are elastic, and fault reactivations are Coulomb frictional failures in the tail of a spatially invariant initial stress distribution. Aftershocks are included according to the empirical ETAS model. This model yields a statistically significant fit to, and a physical explanation of, the observed space-timemagnitude evolution of induced seismicity as a function of induced stress changes.
Pseudo-prospective test results indicate probabilistic seismicity forecasts based on an ensemble of history-matched model instances are consistent with the response observed to subsequent changes in the rate and spatial distribution of gas production that alter reservoir stress states. On this basis, the model appears suitable to support probabilistic seismic hazard and risk analysis. Nonetheless, the reliability of these seismicity forecasts remains limited by our model assumptions, the resolution of observed deformations and seismicity, and the presumption that future deformation and seismicity mechanisms will not significantly differ from the past. We seek to limit forecast uncertainties by restricting the forecast period, typically to 5 yrs, by regularly re-evaluating forecast performance and, if necessary, adapting the model.
An exponential-like trend in initial activity rates is a general consequence of fault strengths being independently and identically distributed such that initial frictional failures occur in the tail of a single unknown probability distribution as described by Extreme Threshold Theory. In the limit of small incremental Coulomb stresses, relative to the variation in initial stresses, C 1/θ 1 , this exponential trend is indistinguishable from a linear trend. However, in pseudo-prospective tests the exponential trend significantly out-performed the linear trend (Bourne & Oates 2017b) indicating the Groningen field exceeds this linear limit and so the number of events per unit strain increases with increasing reservoir strain (Fig. 1b).
Within this model, at constant reservoir strain rate the expected seismicity rate, λ, doubles for every increment of strain zz , where zz = −log (2)/(θ 1 γ H ). This doubling strain depends on geometric heterogeneity, , elastic heterogeneity, H, and initial fault stress heterogeneity θ 1 . An instantaneous decrease in reservoir strain rates causes an instantaneous, proportionate and temporary reduction of induced seismicity rates. Thereafter, seismicity rates continue increasing as the doubling strain is unchanged, although the corresponding doubling time increases due to the lower strain rate (Fig. 5a, 2017-2021).
Maintaining a constant activity rate, λ m , requires, according to eqs (3) and (4), exponentially declining strain rates, such thaṫ zz =˙ 0 exp(−θ 1 γ H ( 0 − zz )) where 0 and˙ 0 are the vertical strain and strain rate at the time of intervention. This exponential decline in strain rates requires a commensurate exponential decline in gas production rates. Should failures progress beyond the tail into the body of the initial fault stress distribution, we expect activity rates would transition from an exponential function of incremental Coulomb stress to a linear function of Coulomb stress rate. We do not observe any evidence within Groningen seismicity for any such transition yet. As induced seismicity is not a time-stationary process methods that forecast activity rates as equal to recently observed rates (Petersen et al. 2017) are likely biased and insensitive to planned control measures.
Stochastic variability in our forecast requires proper characterization as it may limit the effectiveness of a production control measure. Aftershocks significantly increase stochastic variability due to magnitude-dependent triggering (eq. 8). We include aftershocks in the forecast as they are numerous (Fig. 5c) and equally likely to generate events with magnitudes that contribute to seismic hazard (Bourne et al. 2015) and risk. Aftershock removal methods (Gardner & Knopoff 1974;Zaliapin et al. 2008) for conventional Probabilistic Seismic Hazard Assessment will underestimate hazard in the presence of aftershocks and create bias due to non-unique removal criteria. Forecasts based on Monte Carlo simulation of a stochastic process that includes aftershocks such as ETAS avoid these concerns. Despite the additional model complexity these parameters are sufficiently well constrained to improve forecast performance (Fig. 6d).
The existence of large and uncontrollable stochastic variability combined with an exponential-like dependence on induced strains means control measures to reduce seismic risk may be more effective and sustainable than reducing seismicity and the associated ground motion hazards. For instance, selective building strengthening during a temporary period of lower seismicity rates achieved through lower production rates will lower risks despite higher hazards returning over time. If designed properly, the risks at peak hazard should be acceptable, otherwise further production controls will be required to allow more time for additional building strengthening or to lower the peak hazard. Although applied here to the Groningen gas field, these physical principles and control measure design considerations are relevant to mitigating seismic risks induced by fluid production or injection elsewhere.

A C K N O W L E D G E M E N T S
The data are available in an open access data repository (Bourne & Oates 2017a). This study used SciPy (Jones et al. 2001) and Matplotlib (Hunter 2007). SJB: conceived and developed the model. SJO: contributed to the data analysis. JvE: takes operational responsibility for the results. The authors thank Nederlandse Aardolie Maatschappij and Shell Global Solutions International for permission to publish.