The Groningen gas field has shown considerable compaction and subsidence since starting production in the early 1960s. The behaviour is understood from the geomechanical response of the reservoir pressure depletion. By integrating surface movement measurements and modelling, the model parameters can be constrained and understanding of the subsurface behaviour can be improved. Such a procedure has been employed to formulate new compaction and subsidence forecasts. The results are put into the context of an extensive review of the work performed in this field, both in Groningen and beyond. The review is used to formulate a way forward designed to integrate all knowledge in a stochastic manner.

Not long after discovery of the Groningen field, gas-production-induced compaction and consequent land subsidence was recognised to be a potential threat to groundwater management in the province of Groningen, in addition to the fact that parts of the province lie below sea level. More recently, NAM's seismological model also pointed to a correlation between reservoir compaction and the observed induced seismicity above the field. In addition to the already existing requirement for accurate subsidence predictions, this demanded a more accurate description of the expected spatial and temporal development of compaction.

Since the start of production in 1963, multiple levelling campaigns have gathered a unique set of deformation measurements used to calibrate geomechanical models. In this paper we present a methodology to model compaction and subsidence, combining results from rock mechanics experiments and surface deformation measurements. Besides the optical spirit-levelling data, InSAR data are also used for inversion to compaction and calibration of compaction models. Residual analysis, i.e. analysis of the difference between measurement and model output, provides confidence in the model results used for subsidence forecasting and as input to seismological models.

Understanding the mechanisms and key parameters controlling depletion-induced seismicity is key for seismic hazard analyses and the design of mitigation measures. In this paper a methodology is presented to model in 2D the static stress development on faults offsetting depleting reservoir compartments, reactivation of the fault, nucleation of seismic instability, and the subsequent fully dynamic rupture including seismic fault rupture and near-field wave propagation. Slip-dependent reduction of the fault's strength (cohesion and friction) was used to model the development of the instability and seismic rupture. The inclusion of the dynamic calculation allows for a closer comparison to field observables such as borehole seismic data compared to traditional static geomechanical models. We applied this model procedure to a fault and stratigraphy typical for the Groningen field, and compared the results for an offset fault to a fault without offset. A non-zero offset on the fault strongly influenced the stress distribution along the fault due to stress concentrations in the near-fault area close to the top of the hanging wall and the base of the footwall. The heterogeneous stress distribution not only controlled where nucleation occurred within the reservoir interval, but also influenced the subsequent propagation of seismic rupture with low stresses inhibiting the propagation of slip. In a reservoir without offset the stress distribution was relatively uniform throughout the reservoir depth interval. Reactivation occurred at a much larger pressure decrease, but the subsequent seismic event was much larger due to the more uniform state of stress within the reservoir. In both cases the models predicted a unidirectional downward-propagating rupture, with the largest wave amplitudes being radiated downwards into the hanging wall. This study showed how a realistic seismic event could be successfully modelled, including the depletion-induced stressing, nucleation, dynamic propagation, and wave propagation. The influence of fault offset on the depletion-induced stress is significant; the heterogeneous stress development along offset faults not only strongly controls the timing and location of a seismic slip, but also influences the subsequent rupture size of the dynamic event.

The Groningen field is the largest onshore gas field in Europe. Continuous production since 1963 has led to induced seismicity starting in the early 1990s. Production measures aimed at lowering the level of seismicity have been implemented since 2014. In this paper we start from an empirical relationship between the cumulative number of seismic events and cumulative gas production. We show that a better way to analyse the data is to relate the ratio of activity rate over production rate versus the cumulative production, such that the model parameters and their uncertainty can be determined. This also allows us to make predictions including the confidence intervals.

Using this model, we first performed regression analysis based on the larger seismic catalogue which includes all recorded events with a magnitude of 1.3 and larger, because we consider this value to be the magnitude of completeness since 1995. We have also performed regression analysis based on a smaller seismic catalogue consisting of all events with a magnitude of 1.5 and larger. This was done in order to be able to compare our forecast with forecasts performed by others. Our prediction for 2016, based on the seismic catalogue of all events with a magnitude of M≥1.5 (using only the events recorded before 2016), was 16±8 events. By the end of 2016, 13 such events had been recorded.

We discuss a number of factors which may influence the predictive power of the derived relationship and which require further study. For instance, we consider the delay between production and earthquakes which increases with decreasing reservoir pressure. In addition, the effect of seasonal fluctuation in Groningen production should be considered. Future work can be done to include these effects in the empirical model. We also investigated the challenges related to the applicability of the analysis to sub-regions of the Groningen field.

Previous locations of earthquakes induced by depletion of the Groningen gas field were not accurate enough to infer which faults in the reservoir are reactivated. A multiplet analysis is performed to identify clusters of earthquakes that have similar waveforms, representing repeating rupture on the same or nearby faults. The multiplet analysis is based on the cross-correlation of seismograms to assess the degree of similarity. Using data of a single station, six earthquake clusters within the limits of the Groningen field were identified for the period 2010 to mid-2014. Four of these clusters were suitable for a relocation method that is based on the difference in travel time between the P- and the S-wave. Events within a cluster can be relocated relative to a master event with improved accuracy by cross-correlating first arrivals. By choosing master events located with a new dense seismic network, the relocated events likely not only have better relative, but also improved absolute locations. For a few clusters with sufficient signal-to-noise detections, we show that the relocation method is successful in assigning clusters to specific faults at the reservoir level. Overall, about 90% of the events did not show clustering, despite choosing low correlation thresholds of 0.5 and 0.6. This suggests that different faults and/or fault segments with likely varying source mechanisms are active in reservoir sub-regions of a few square kilometres.

This paper reviews the evolution of a sequence of seismological models developed and implemented as part of a workflow for Probabilistic Seismic Hazard and Risk Assessment of the seismicity induced by gas production from the Groningen gas field. These are semi-empirical statistical geomechanical models derived from observations of production-induced seismicity, reservoir compaction and structure of the field itself. Initial versions of the seismological model were based on a characterisation of the seismicity in terms of its moment budget. Subsequent versions of the model were formulated in terms of seismic event rates, this change being driven in part by the reduction in variability of the model forecasts in this domain. Our approach makes use of the Epidemic Type After Shock model (ETAS) to characterise spatial and temporal clustering of earthquakes and has been extended to also incorporate the concentration of moment release on pre-existing faults and other reservoir topographic structures.

In the Netherlands, over 190 gas fields of varying size have been exploited, and 15% of these have shown seismicity. The prime cause for seismicity due to gas depletion is stress changes caused by pressure depletion and by differential compaction. The observed onset of induced seismicity due to gas depletion in the Netherlands occurs after a considerable pressure drop in the gas fields. Geomechanical studies show that both the delay in the onset of induced seismicity and the nonlinear increase in seismic moment observed for the induced seismicity in the Groningen field can be explained by a model of pressure depletion, if the faults causing the induced seismicity are not critically stressed at the onset of depletion. Our model shows concave patterns of log moment with time for individual faults. This suggests that the growth of future seismicity could well be more limited than would be inferred from extrapolation of the observed trend between production or compaction and seismicity. The geomechanical models predict that seismic moment increase should slow down significantly immediately after a production decrease, independently of the decay rate of the compaction model. These findings are in agreement with the observed reduced seismicity rates in the central area of the Groningen field immediately after production decrease on 17 January 2014. The geomechanical model findings therefore support scope for mitigating induced seismicity by adjusting rates of production and associated pressure change. These simplified models cannot serve as comprehensive models for predicting induced seismicity in any particular field. To this end, a more detailed field-specific study, taking into account the full complexity of reservoir geometry, depletion history and mechanical properties, is required.

Major efforts are being undertaken to quantify seismic hazard and risk due to production-induced earthquakes in the Groningen gas field as the basis for rational decision-making about mitigation measures. An essential element is a model to estimate surface ground motions expected at any location for each earthquake originating within the gas reservoir. Taking advantage of the excellent geological and geophysical characterisation of the field and a growing database of ground-motion recordings, models have been developed for predicting response spectral accelerations, peak ground velocity and ground-motion durations for a wide range of magnitudes. The models reflect the unique source and travel path characteristics of the Groningen earthquakes, and account for the inevitable uncertainty in extrapolating from the small observed magnitudes to potential larger events. The predictions of ground-motion amplitudes include the effects of nonlinear site response of the relatively soft near-surface deposits throughout the field.

The shallow subsurface of Groningen, the Netherlands, is heterogeneous due to its formation in a Holocene tidal coastal setting on a periglacially and glacially inherited landscape with strong lateral variation in subsurface architecture. Soft sediments with low, small-strain shear wave velocities (VS30 around 200 m s−1) are known to amplify earthquake motions. Knowledge of the architecture and properties of the subsurface and the combined effect on the propagation of earthquake waves is imperative for the prediction of geohazards of ground shaking and liquefaction at the surface. In order to provide information for the seismic hazard and risk analysis, two geological models were constructed. The first is the ‘Geological model for Site response in Groningen’ (GSG model) and is based on the detailed 3D GeoTOP voxel model containing lithostratigraphy and lithoclass attributes. The GeoTOP model was combined with information from boreholes, cone penetration tests, regional digital geological and geohydrological models to cover the full range from the surface down to the base of the North Sea Supergroup (base Paleogene) at ~800 m depth. The GSG model consists of a microzonation based on geology and a stack of soil stratigraphy for each of the 140,000 grid cells (100 m × 100 m) to which properties (VS and parameters relevant for nonlinear soil behaviour) were assigned. The GSG model serves as input to the site response calculations that feed into the Ground Motion Model. The second model is the ‘Geological model for Liquefaction sensitivity in Groningen’ (GLG). Generally, loosely packed sands might be susceptible to liquefaction upon earthquake shaking. In order to delineate zones of loosely packed sand in the first 40 m below the surface, GeoTOP was combined with relative densities inferred from a large cone penetration test database. The marine Naaldwijk and Eem Formations have the highest proportion of loosely packed sand (31% and 38%, respectively) and thus are considered to be the most vulnerable to liquefaction; other units contain 5–17% loosely packed sand. The GLG model serves as one of the inputs for further research on the liquefaction potential in Groningen, such as the development of region-specific magnitude scaling factors (MSF) and depth–stress reduction relationships (rd).

The increase in number and strength of shallow induced seismicity connected to the Groningen gas field since 2003 and the occurrence of a ML 3.6 event in 2012 started the development of a full probabilistic seismic hazard assessment (PSHA) for Groningen, required by the regulator. Densification of the monitoring network resulted in a decrease of the location threshold and magnitude of completeness down to ~ ML=0.5. Combined with a detailed local velocity model, epicentre accuracy could be reduced from 0.5–1km to 0.1–0.3km and a vertical resolution ~0.3km. Time-dependent seismic activity is observed and taken into account into PSHA calculations. Development of the Ground Motion Model for Groningen resulted in a significant reduction of the hazard. Comparison of different implementations of the PSHA, using different source models, based on either a compaction model and production scenarios or on extrapolation of past seismicity, and methods of calculation, shows similar results.

This paper describes the ongoing experimental and analytical activities that are being carried out to develop fatality and consequence models for the estimation of ‘Inside Local Personal Risk’ (ILPR) of buildings within the Groningen field. ILPR is defined as the annual probability of fatality for a hypothetical person who is continuously present without protection inside a building. In order to be able to estimate this risk metric, a robust estimate of the probability of collapse of structural and non-structural elements within a building is needed, as these have been found to be the greatest drivers of fatality risk.

To estimate the collapse potential of buildings in Groningen, structural numerical models of a number of representative case studies have been developed and calibrated through in situ and laboratory testing on materials, connections, structural components and even full-scale buildings. These numerical models are then subjected to increased levels of ground shaking to estimate the probability of collapse, and the associated consequences are estimated from the observed collapse mechanisms.

Earthquakes associated with gas production have been recorded in the northern part of the Netherlands since 1986. The Huizinge earthquake of 16 August 2012, the strongest so far with a magnitude of ML = 3.6, prompted reassessment of the seismicity induced by production from the Groningen gas field. An international research programme was initiated, with the participation of many Dutch and international universities, knowledge institutes and recognised experts.

The prime aim of the programme was to assess the hazard and risk resulting from the induced seismicity. Classic probabilistic seismic hazard and risk assessment (PSHA) was implemented using a Monte Carlo method. The scope of the research programme extended from the cause (production of gas from the underground reservoir) to the effects (risk to people and damage to buildings). Data acquisition through field measurements and laboratory experiments was a substantial element of the research programme. The existing geophone and accelerometer monitoring network was extended, a new network of accelerometers in building foundations was installed, geophones were placed at reservoir level in deep wells, GPS stations were installed and a gravity survey was conducted.

Results of the probabilistic seismic hazard and risk assessment have been published in production plans submitted to the Minister of Economic Affairs, Winningsplan Groningen 2013 and 2016 and several intermediate updates. The studies and data acquisition further constrained the uncertainties and resulted in a reduction of the initially assessed hazard and risk.

Shaking and damage in the province of Groningen, the Netherlands, resulting from production-induced seismicity has caused increased public anxiety. Since 2014, production offtake has been reduced stepwise by over 50% in an attempt to minimise production-induced seismicity. The earthquake catalogue, combined with comprehensive data of the changes in production offtake, shows a clear response of seismic activity following the production measures taken. Associated temporal variations in the proportionality between smaller- and larger-magnitude events (the b-value of the Gutenberg–Richter relation) are observed. Since production measures were imposed, the b-value has tended to increase, thus lowering the probability of a larger-magnitude event. The analysis also shows increases in activity rate and b-value prior to larger-magnitude events. Subsequently, the probability of a larger-magnitude event seems to be decreasing prior to the events occurring. This implies that for short-term earthquake prediction of hydrocarbon-production-induced seismicity, these types of analysis could be misleading. However, regional analysis is necessary to explain the observations in terms of rupture initiation. At present, each event felt still draws the interest of both public and press. As some clustering of events in both time and space is still observed, managing both the seismicity and the public perception provides a continuing challenge.

Further research into seismicity caused by natural gas production from the Groningen field is necessary to improve the assessment of seismic risk and develop means to control and reduce it. Research into subsurface aspects is primarily of relevance to assess the seismic hazard component in the cause-and-effect chain that governs the seismic risk. It requires a wide range of research activities that can be broadly classified as follows:

• • Increasing understanding of the physical mechanisms that govern production-induced seismicity, in particular source mechanisms, compaction behaviour, propagation of energy to the surface, and the effects of fluctuating production.

• • Reducing uncertainty by acquiring additional field data to improve statistical inference, and developing statistical methods and procedures that can cope with the non-stationary nature of the process.

• • Developing tools and techniques to improve risk management, and support operational control and policy measures under uncertainty.

An essential requirement for further research will be the possibility of developing competing theories for many aspects of the modelling chain. This requires an overall hazard and risk assessment methodology that can accommodate multiple models, and an organisational structure that facilitates the comparison of competing approaches while safeguarding their independent development. This will have to be supported by the availability of reliable data via shared databases. Finally, the scientific community should be prepared to make a major effort to translate their research results into popular scientific versions in order to keep stakeholders abreast of progressive insight into the origin, predictability and prevention of induced seismicity.