This special issue of the Netherlands Journal of Geosciences / Geologie en Mijnbouw contains research papers offering important new insights into aspects of the occurrence of gas-production-induced seismicity in Groningen, the Netherlands. The issue is a ‘work in progress’ and provides a snapshot of our knowledge at the end of 2016, which is when most of the contributions were written.

After more than half a century of production and with some 350 wells, the Groningen gas field must be one of the best-studied gas fields in the world. Initially, it was considered to be relatively simple and behaving like one big tank. Now that it is entering a phase of declining production it has become clear that many subtleties are not fully understood yet. Prediction and management of subsidence and induced earth tremors require a detailed understanding of the field geology. In addition, an optimum gas recovery is only possible if details of, for example, reservoir quality distribution and faulting, that did not appear relevant before, are well understood.

The large Groningen field comprises a structurally high block during much of its history, probably already from Devonian times onwards. The desert sandstones of the Rotliegend reservoir exhibit a strong south-to-north proximal–distal relationship. Whilst diagenesis has in many fields led to deterioration of reservoir properties, this effect is small in the Groningen field. The field is dipping to the north, and bounded by a series of normal faults in the west, south and east. Almost all faults are normal extensional faults, but locally inverse reactivation has led to small pop-up structures. Reactivation of older faults must have resulted in oblique movements along most faults.

The challenges for future development of the Groningen field are immense. Managing the risks associated with induced seismicity and recovery of the remaining gas will continue to require an increasingly detailed knowledge and understanding of its geology.

Prediction of gas-production-induced subsidence and seismicity is much more difficult and uncertain than generally recognised in the past. It is now widely accepted that uncertainties in predicted subsidence and seismicity are large prior to and during the initial stages of production. At later stages, predictions remain highly uncertain for periods more than three to five years into the future. This requires a different regulatory framework to ensure that associated risks remain within accepted boundaries. Previously, single-scenario operator predictions were checked against field measurements. When subsidence or seismicity started to deviate beyond claimed uncertainties, the operator was asked to provide prediction updates. The practice was long considered acceptable, as structural damage to buildings and infrastructure or personal risk to people was not expected. This all changed following the 2012 Huizinge seismic event, necessitating better identification, assessment and ranking of risks, the use of scenarios, probabilistic forecasting and a much intensified field monitoring and control loop. It requires that the regulator becomes actively involved in assuring the integrated control loop of risk identification, predictions, monitoring, updating, mitigation measures and the closing of knowledge gaps, to ensure that subsidence (rate) and induced seismicity remain within acceptable limits. And it requires that this increased involvement of the regulator is supported in the mining law and by appropriate conditions in the Production Plan assent.

Depletion of the Groningen gas field has induced earthquakes, although the north of the Netherlands is a tectonically inactive region. Increased seismic activity raised public concern which led the government to initiate a number of studies with the aim of understanding the cause(s) of the earthquakes. If the relationship between production and seismicity were understood then production could be optimized in such a way that the risk of induced seismicity would be minimal. The main question remains how production is correlated with induced seismicity. The Minister of Economic Affairs of the Netherlands decided to reduce production starting from 17 January 2014, specifically in the centre of the gas field as it has the highest rates of seismicity, the largest-magnitude events and the highest compaction values of the field.

A reduction in production could possibly lead to a reduced rate of compaction. Additionally a reduction of production rate could lead to a reduced stress rate increase on the existing faults and consequently fewer seismic events per year. One might envisage a ‘bonus effect’ in the events reduction in the sense that the total number of events will be less, with the same total production smeared out over a longer period. This is as yet unclear.

In this paper we apply different statistical methods to look for evidence supporting or disproving a decrease in the number of seismic events due to production reduction.

The assessment of the seismic hazard and risk associated with the extraction of gas from the Groningen field involves a chain of modelling efforts. The first step is a description of the 3D distribution of reservoir properties in the reservoir – the static reservoir model – and is the subject of this paper. Consecutive steps in the chain of models are described elsewhere in this volume. The construction of a static reservoir model is not strictly a scientific endeavour, but many of the applied modelling techniques are underpinned by extensive scientific research. This paper aims to give a general introduction to the approach followed by NAM to build static models for the Groningen field. More detailed accounts of the applied modelling techniques, the assessment of associated uncertainties or the usage of multiple modelling scenarios are beyond the scope of the current paper, but are referenced in the text.

This paper presents the method applied to history-match the Groningen field dynamic model to conventional data (pressure data and water influx data) and to subsidence data, which is a novelty. Modelled subsidence is matched to subsidence data based on a simplified geomechanical model, which was built into the dynamic simulator.

A two-tier method was used to first achieve a match on a field-wide scale using field-average history-match quantifiers, which was then further improved at a regional/well level using regional history-match quantifiers. The history match was assisted by a space-filling experimental design. The resulting model has a field-average match to pressure data of ±2.17bar with a measurement uncertainty of ±0.4bar, to water influx data of ±2m with a measurement uncertainty of ±0.5m, and to subsidence data of ±4cm with a measurement uncertainty of ±1cm.

The output from this model is used as input for compaction, subsidence and production forecasts feeding into the hazard and risk assessment completed by NAM for the Groningen Winningsplan 2016.

This paper describes a research programme recently initiated at Utrecht University that aims to contribute new, fundamental physical understanding and quantitative descriptions of rock and fault behaviour needed to advance understanding of reservoir compaction and fault behaviour in the context of induced seismicity and subsidence in the Groningen gas field. The NAM-funded programme involves experimental rock and fault mechanics work, microscale observational studies to determine the processes that control reservoir rock deformation and fault slip, modelling and experimental work aimed at establishing upscaling rules between laboratory and field scales, and geomechanical modelling of fault rupture and earthquake generation at the reservoir scale. Here, we focus on describing the programme and its intended contribution to understanding the response of the Groningen field to gas production. The key knowledge gaps that drive the programme are discussed and the approaches employed to address them are highlighted. Some of the first results emerging from the work in progress are also reported briefly and are providing important new insights.

The Groningen field is the largest onshore gas field in Europe. The gas-bearing section comprises aeolian and fluvial Rotliegend sandstones of Permian age and fluvial sandstones of Carboniferous age. Continuous production since 1963 has led to induced seismicity starting in the early 1990s.

Faults at reservoir level play a major role in the seismicity in the Groningen field. Fault slip is expected when shear traction is sufficient to overcome frictional resistance on the fault surface. Clear insights into which faults and fault segments are most susceptible to seismicity could be used to optimise production and minimise the seismic risk. To gain these insights, a detailed and realistic fault model is required as input to both statistical analyses on seismicity and deterministic geomechanical modelling of seismogenic behaviour along fault planes. Geometrical seismic attributes and, subsequently, fault planes were extracted from a reprocessed and depth-imaged 3D seismic volume. This resulted in a detailed visualisation of the faults at reservoir level, with extension into the deeper strata below the reservoir in many cases. They represent fault planes with realistic dimensions and shapes. The fault map based on seismic attributes suggests the presence of faults that have not been included in studies on Groningen seismicity before. The improved fault definition correlates with recent earthquake hypocentres. We conclude that a detailed fault model of the Groningen field can be created using 3D seismic attributes and that detailed 3D fault planes can be extracted from these attributes. The results can be used as input to statistical and geomechanical analyses on seismicity.

The presence of salt in dilatant normal faults may have a strong influence on fault mechanics in the Groningen field and on the related induced seismicity. At present, little is known of the structure of these fault zones. This study starts with the geological evolution of the Groningen area, where, during tectonic faulting, rock salt may have migrated downwards into dilatant faults. These fault zones therefore may contain inclusions of rock salt. Because of its rate-dependent mechanical properties, the presence of salt in a fault may introduce a loading-rate dependency into fault movement and affect the distribution of magnitudes of seismic events. We present a first-look study showing how these processes can be investigated using a combination of analogue and numerical modelling. Full scaling of the models and quantification of implications for induced seismicity in Groningen require further, more detailed studies: an understanding of fault zone structure in the Groningen field is required for improved predictions of induced seismicity. The analogue experiments are based on a simplified stratigraphy of the Groningen area, where it is generally thought that most of the Rotliegend faulting has taken place in the Jurassic, after deposition of the Zechstein. This suggests that, at the time of faulting, the sulphates were already transformed into brittle anhydrite. If these layers were sufficiently brittle to fault in a dilatant fashion, rock salt was able to flow downwards into the dilatant fractures. To test this hypothesis, we use sandbox experiments where we combine cohesive powder as analogue for brittle anhydrites and carbonates with viscous salt analogues to explore the developing fault geometry and the resulting distribution of salt in the faults. Using the observations from analogue models as input, numerical models investigate the stick-slip behaviour of fault zones containing ductile material qualitatively with the discrete element method (DEM). Results show that the DEM approach is suitable for modelling the seismicity of faults containing salt. The stick-slip motion of the fault becomes dependent on shear loading rate with a modification of the frequency–magnitude distribution of the generated seismic events.

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.