Introduction
The Groningen field, with an estimated total volume of more than 2800bcm (billion cubic metres) of gas, was discovered in 1959 and is still the largest known gas field in western Europe (De Waal et al., Reference de Waal, Muntendam-Bos and Roest2015). Induced seismicity in the Groningen field was first registered in December 1991. Until 2003 the seismic activity rate was low, reasonably constant and located at the centre of the field. Between 2003 and 2013 seismic activity increased exponentially (Bourne et al., Reference Bourne, Oates, Bommer, Dost, van Elk and Doornhof2015; Nepveu et al., Reference Nepveu, van Thienen-Visser and Sijacic2016; Pijpers & de Waal, Reference Pijpers2016b,Reference Pijpersc). While the rate increased, so did the magnitudes of the largest events. The first magnitude 3 events were recorded in 2003, followed by the first magnitude 3.5 event in 2006. The largest event to date (M l=3.6) was recorded on 12 August 2012 near the town of Huizinge (Dost & Kraaijpoel, Reference Dost and Kraaijpoel2013). The event caused significant non-structural damage throughout the region. In an area devoid of any known natural seismicity, the shaking and damage led to anxiety and a great deal of public turmoil.
A State Supervision of Mines (Staatstoezicht op de Mijnen (SodM)) investigation (Muntendam-Bos & De Waal, Reference Muntendam-Bos and de Waal2013) following the Huizinge earthquake showed an increased risk of larger-magnitude events occurring due to gas extraction. Based on the results of this study, the regulator advised reducing gas production as much and as fast as realistically possible. Starting in January 2014, the Dutch Minister of Economic Affairs (MEA) gradually limited production offtake in the Groningen gas field in five steps, from a production of 54bcma−1 in 2013 to 24bcma−1 in its latest decision on 1 October 2016. At a later stage, studies initiated by the regulator (Pijpers, Reference Pijpers2016a) suggested that production-rate changes in themselves play an important role in triggering seismicity. Hence, in the current production philosophy the large summer–winter production rate fluctuations are avoided as much as possible.
Seismic risk cannot be measured, but only quantified through seismic risk assessment. This assessment is associated with a sequence of models all containing large uncertainties and many assumptions (NAM, 2016; SodM, 2016). As a consequence the impact of the measures imposed on the assessed seismic risk appears to be minor (NAM, 2016). Of the sequence of models, only the seismological model could be affected by changes in production rates. In this paper, we demonstrate that the impact of the measures on seismic activity and the probability of larger-magnitude events is significant. We make a systematic assessment of the seismic response to all five measures imposed. In addition, we assess the impact of the measures on the relation between larger- and smaller-magnitude events and draw conclusions for the probability of larger-magnitude events based on this assessment.
Risk posed by induced seismicity
In the Netherlands, seismicity is induced directly by the extraction of natural gas from permeable sand reservoirs (Van Eck et al., Reference van Eck, Goudbeek, Haak and Dost2006; Van Eijs et al., Reference van Eijs, Mulders, Nepveu, Kenter and Scheffers2006; Van Wees et al., Reference van Wees, Buijze, van Thienen-Visser and Fokker2014; Muntendam-Bos et al., Reference Muntendam-Bos, Roest and de Waal2015). Reported local magnitudes (M l) range from below 0 to 3.6. Although seismicity is observed in numerous smaller gas fields, especially in the north of the country, the large Groningen gas field in the northeast of the Netherlands is the most seismically active field, with more than 1000 registered events to date (Fig. 1).
(A) The Groningen gas field on a map of the region and its seismicity as reported by the Dutch Meteorological Institute (KNMI). The colour coding of the seismicity indicates the magnitude of the events. The ellipse indicates the central graben (CG) structure in the field, where the initial production reduction measures were imposed. (B) The development of the Groningen seismicity with time. The solid line denotes the growth of the cumulative number of events of magnitude 1.4 and higher. The colour bars denote the annual number of earthquakes in different magnitude classes.
For production-induced seismicity the primary physical mechanism is the decreasing pore pressure within the gas field due to gas withdrawal. The decrease in pressures causes the reservoir to compact, which induces subsidence at the earth's surface. The decrease in pressures also increases the stresses on pre-existing faults, which in general would tend to stabilise faults. However, the Netherlands being a normal faulting regime, the increase of the effective vertical stresses is faster than the increase of the effective horizontal stresses. Hence, stress conditions on the faults may shift into the shear-failure regime (Roest & Kuilman, Reference Roest, Kuilman and Balkema1994; Zoback, Reference Zoback2007). The effect is significantly enhanced by offsets of the reservoir layer over the fault zone (Mulders, Reference Mulders2003). Studies indicate that the largest impact on shear stresses on a fault results from differential compaction (Roest & Kuilman, Reference Roest, Kuilman and Balkema1994; Mulders, Reference Mulders2003; van den Bogert, Reference van den Bogert2015). As at a certain time these stresses may act on a larger fault area, seismicity could commence at a comparatively high-magnitude level, as is observed at some fields in the Netherlands (Van Eck et al., Reference van Eck, Goudbeek, Haak and Dost2006).
A reduction in production rate may influence the rate of compaction and could lead to a reduced shear stress loading on the faults. Additionally, a reduction in production rate influences the rate at which stress increases on the existing faults. If seismicity is stress-rate dependent, this will result in lower seismic activity rates.
The earthquake catalogue
The available earthquake dataset for the Groningen gas field contains in total 987 events recorded between 1 January 1991 and 1 January 2017. Detailed information on the Groningen network and dataset can be found in Dost et al. (Reference Dost, Ruigrok and Spetzler2017). The network of seismometers in the Groningen region was designed to be complete for earthquakes above magnitudes of 1.5. The network was only fully operational from 1995 onwards. To ensure a reasonably uniform quality of the catalogue, it is preferable to exclude all data prior to 1 January 1995. A recent improvement of the seismic network has had an impact primarily on detection of earthquakes with magnitudes below 1 during and after 2015. Paleja and Bierman (Reference Paleja and Bierman2016) report on an initial analysis of temporal variations in completeness magnitude (M c) for the entire Groningen field. A reanalysis of the M c by Dost et al. (Reference Dost, Ruigrok and Spetzler2017) confirms this temporal variation. They conclude that for the two initial time periods April 2003–February 2009 and August 2009–August 2012 M c=1.2, for the time period August 2012–August 2014 M c=0.8 and for the period September 2014–September 2016 M c=0.5. However, they remark that for the last period the procedure was not very stable.
Induced activity rate
In January 2014, production in five wells surrounding the most active, central region of the Groningen gas field was reduced by 80% (Fig. 1). In order to minimise the reduction of the full field production, the MEA allowed an increase in production in the other areas of the field. In response, especially production in the southwestern clusters of the field was increased. Subsequently, the seismic activity associated with the main NW–SE-orientated central graben structure in the centre of the field reduced significantly (Fig. 2). Support for the significance of the observed changes is provided by independent statistical analysis of the changes in seismicity patterns (Nepveu et al., Reference Nepveu, van Thienen-Visser and Sijacic2016; Pijpers, Reference Pijpers2016b). At the same time, an increase in seismic activity was observed in a second graben structure in the southwest of the Groningen field. The analysis of Pijpers (Reference Pijpers2016b) shows the tentative significance of this increase, which due to the limited data is not conclusive. The increase is clearly correlated to the increase in production in the southwestern clusters.
Event density of the seismicity in the Groningen gas field, plotted on a map of the region, in three successive time periods: (A) March 2013–March 2014, (B) March 2014–March 2015 and (C) March 2015–March 2016. The colour scale runs from 0 (blue) to 0.5 events km−2 (red).
Along with the increase in seismic activity, local public anxiety increased. Public concern focused on the apparent shift of the seismic activity from the centre of the field to the southwest. In January 2015, authorisation of the increase in production in the southwestern clusters was revoked. In each of the three subsequent production reduction steps, further limitations on production offtake were evenly distributed over all producing clusters except the five central clusters, which were kept at near zero production.
Throughout the field history, production from the Groningen field has been subject to seasonal swing. Production rates have been a factor of 2–7 higher in winter than in summer. A medium-strength statistical correlation between the seasonal fluctuation and seismic activity was identified by Nepveu et al. (Reference Nepveu, van Thienen-Visser and Sijacic2016). By analysing the local gas pressure history at the locations of the induced earthquakes, Pijpers (Reference Pijpers2016a) investigated whether a characteristic set of circumstances provokes the occurrence of the seismic events. The analysis shows that the timing of seismic events correlates with increases in local depletion rates.
Hence, as an additional operational measure, the seasonal fluctuation has been minimised effectively since March 2015. The analysis of Pijpers (Reference Pijpers2016c) shows that the seismic activity since March 2015 is consistent with a reversal of the exponentially increasing trend, reducing the tremor rates to the 12-year average rate between 1 September 2004 and 1 September 2016. In the centre of the field, the tremor rate has decreased to half the 12-year average rate, while depletion rates have declined to approximately half the 2013 rate. The reduction in tremor rates in the field centre is stronger than would be expected purely on the basis of the reduction in gas depletion rates.
Gutenberg–Richter
In addition to the above assessment of the activity rate, we have analysed the earthquake catalogue for possible changes in the frequency–magnitude, or Gutenberg–Richter (GR), relation. This relation addresses both changes in activity rate and the ratio between smaller- and larger-magnitude events.
The frequency–magnitude relation (Gutenberg & Richter, Reference Gutenberg and Richter1944, Reference Gutenberg and Richter1954) is given by:
$$\begin{equation*}
{\log _{10}}N\left( M \right) = a - bM
\end{equation*}$$
where N(M) is the annual number of earthquakes which occur in a given area having a magnitude ≥M, the constant a is a measure of the level of seismicity, while the constant b describes how the number of earthquakes in the given area and time period varies for different magnitudes (it is the negative of the slope of the GR relationship).
In order to avoid ambiguities due to non-consistent datasets or problems with the magnitude of completeness, the earthquake catalogue has first been divided into bins, each containing the same number of N(M≥M c)-events. Although Dost et al. (Reference Dost, Ruigrok and Spetzler2017) derived a time dependence in M c due to the extension of the network, for consistency and ease of comparison we used M c=1.2 for the full catalogue between 1995 and 2017. The total number of events N (M≥M c) in each bin has been constrained to the total number of 56 (M≥M c)-events that occurred in the years 2015 and 2016 combined. By constraining the bin size to the number of events in 2015 and 2016 and determining subsequent bins backwards in time (hence start counting in 2017), we ensure that possible changes in activity rate and b-value due to the production reductions in 2014 are captured. Using a maximum likelihood estimator (Mignan & Woessner, Reference Mignan and Woessner2012), the b-value for each bin-dataset has been determined. Figure 3A shows the annualised frequency magnitude distributions (FMDs) for six of the eight time periods (bins). The intermediate bins for the time periods September 2013–January 2015 and April 2009–April 2011 have been omitted for clarity of the figure. However, the FMDs for these periods look very similar to the ones presented. The derived maximum likelihood b-values of all eight time periods including 1σ uncertainty estimate (Marzocchi & Sandri, Reference Marzocchi and Sandri2003) are given in Figure 3B, and the corresponding annualised cumulative number of events of M≥M c in Figure 3C.
(A) FMDs of six out of the eight time period bins (colour-coded). (B) Maximum likelihood estimate of the b-value of the GR relation, including the 1σ uncertainty estimate for each bin-dataset. The dotted lines indicate the bin-boundaries in time. (C) The annualised cumulative number of events of M≥1.2. The black dashed lines indicate the bin-boundaries in time.
A clear variation in the b-value over the different time periods is observed. Periods of higher b-value tend to be followed by periods of lower b-value. Since the production measures were imposed in January 2014, a subsequent continuous increase in b-value is obtained for the periods September 2013–January 2015 and January 2015–January 2017. The variation in the b-value is significant on a 1σ uncertainty level. However, the variations in b-value are not statistically significant on a 2σ or 3σ uncertainty level. As a consequence of the small bin-datasets (56 events), the uncertainties are rather large, which results in a non-conclusive analysis.
To avoid the dependence on the choice of end time/event, we follow Gulia et al.’s (Reference Gulia, Tormann, Wiemer, Hermann and Seif2016) analysis for deriving a continuous b-value time series. We use a window length of both 56 events and 100 events, which is moved through the catalogue event-by-event, thus exploring the full range of variability in the data. The window length of 56 events is consistent with the bin size used. Figure 4A and B show the derived maximum likelihood b-values (Mignan & Woessner, Reference Mignan and Woessner2012) and their 1σ standard deviation (Marzocchi & Sandri, Reference Marzocchi and Sandri2003). To avoid confusion about causative relations, the values are plotted at the end of each time window. We find quite some variation of the b-value for the Groningen field seismic history. Remarkably, the field-wide b-value seems to increase significantly prior to the two largest seismic events in the field: M=3.5 in August 2006 and M=3.6 in August 2012 (grey dashed lines in Fig. 4). This signal is robust with respect to different window lengths.
(A) Temporal b-value evolution throughout the Groningen earthquake catalogue with a 56-event window length. (B) Temporal b-value evolution throughout the Groningen earthquake catalogue with a 100-event window length. (C) Temporal evolution of the annualised cumulative number of events M≥1.2 with a 56-event window length. (D) Temporal evolution of the annualised cumulative number of events M≥1.2 with a 100-event window length. Grey dashed lines: the two M≥3.5 events that occurred in the Groningen field.
The interpretation of the observed increase in b-value in terms of rupture initiation is not straightforward. The increase in b-value is opposite to the distinct low pre-main shock b-value found for natural seismicity (e.g. Tormann et al., Reference Tormann, Wiemer and Hardebeck2012; Schurr et al., Reference Schurr2014; Gulia et al., Reference Gulia, Tormann, Wiemer, Hermann and Seif2016) and observations from laboratory experiments (Meredith et al., Reference Meredith, Main and Jones1990; Main, Reference Main1996; Goebel et al., Reference Goebel, Becker, Schorlemmer, Stanchits, Sammis, Rybacki and Dresen2012). While it is feasible that the b-value is increased due to a build-up of energy on a fault due to asperities on the fault, an alternative model suggests that the b-value increases because the activity rate in an inherent higher b-value area increases. In addition, because in this study we have focused on an analysis of the full field, the observation could be due to an increase in seismicity in a different part of the gas field due to e.g. a local increase in stressing rate.
Hence, in order to properly explain the observed changes in b-value it is crucial to extend our current analysis to a more regional analysis of the changes in b-value. Bourne et al. (Reference Bourne, Oates, Bommer, Dost, van Elk and Doornhof2015) obtained a possible decrease in b-value with increasing compaction. Since 2014, seismicity in the centre of the field, where most compaction has occurred, has decreased and seismicity has mainly been focused in regions of lower compaction. Hence in the Bourne et al. (Reference Bourne, Oates, Bommer, Dost, van Elk and Doornhof2015) model, an increase in b-value for a full field analysis is to be expected. However, an increase in the b-value in the higher-compaction areas would not be feasible. A reanalysis including the observations after production measures were implemented may provide more insight on the development of the b-value in the Groningen gas field in relation to the compaction in the field. An alternative explanation of the variation in b-value could come from the much lower stressing rates which are currently being imposed on the faults especially in the centre of the field. A regional analysis of the changes in b-value and a possible relation with stressing rate is the subject of ongoing research.
The variation in annual activity rate (Figs 3C, 4C, D) corresponds to the observations and analysis in the previous section, with a strongly increasing activity rate up to mid-2014 peaking in 2012 prior to the Huizinge M=3.6 earthquake, and a mainly decreasing rate since. Also prior to the M=3.5 earthquake in August 2006 an increase in activity rate is observed, especially for the shorter window length of 56 events (Fig. 4C).
Implications for seismic hazard and risk
The occurrence probability of an event with magnitude larger than M in a specified time interval is given by:
$$\begin{equation*}
P\left( {0\left| M \right|t} \right) = 1 - \exp ( - N\left( {M\left( t \right)} \right) = 1 - \exp \left( { - {{10}^a}^{ - bM(t)}} \right)
\end{equation*}$$
The occurrence probability of events quickly decreases with increasing magnitudes. Figure 5 shows the annual probability for events with magnitudes larger than a given one to occur in each time period bin of the Groningen gas field, where the annual activity rates and b-values of the previous section have been adopted. Taking into consideration the changes in both the b-value and the activity rate, the annual probability of a magnitude M4+ event was largest (~0.26) in the period April 2012–September 2013: the period in which the M=3.6 Huizinge event occurred. Since production measures have been imposed, the annual probability of an M4+ event has decreased to ~0.04. The bin containing the M=3.5 event of August 2006 (July 2006–April 2009) shows an M4+ annual probability of ~0.07.
The annual occurrence probability of an event with magnitude larger than a given one to occur in each bin-dataset of the Groningen gas field utilising the annual event rate and the maximum likelihood estimates of the b-values derived in the time period analysis.
The time series analysis of the occurrence probability of a M4+ event is given in Figure 6. The derived probabilities are strongly dependent on the changes in activity rate and b-value, as well as on the window length used in the time series. For instance, the analysis using the 56 events window length shows a clear increased probability for an M4+ event by the end of 2013, while the analysis based on 100 events shows a less pronounced increase. However, these high probabilities of an M4+ event in this region of very damage-susceptible buildings does support the necessity of taking measures to reduce the seismic risk. As would be expected, the increases in b-values prior to the two largest events result in low probabilities for a M4+ event just prior to their occurrence and a strong increase in probability just after the events. This implies that for short-term earthquake prediction of hydrocarbon-production-induced seismicity these types of analysis could be misleading. Hence, much care needs to be taken and much further study is needed.
(A) Temporal evolution of the annual probability of an M4+ throughout the Groningen catalogue estimated from the temporal evolution of the b-value and the temporal evolution of the annualised cumulative number of events M≥1.2 both with a 56-event window length. (B) Same as (A), but utilising the temporal evolutions with a 100-event window length. Grey dashed lines: the two M≥3.5 events that occurred in the Groningen field.
Conclusions
Our analysis shows that the changes made both in the production pattern and in the production rates are followed by changes in seismic activity. Along with the decrease in seismic activity rate over the last 2 years a tendency for a decrease in larger-magnitude events is observed. Temporal changes in the slope of the GR relation (b-value) are clearly observed on a 1σ uncertainty level. We note that our observations are based on limited data and at present no model exists to validate our observations.
Strikingly, the changes indicate increases in activity rate and b-values prior to the larger events. Hence, while based on these observations the probability of a larger-magnitude event seems to be decreasing prior to a larger event, evidence shows it may actually be more likely. This implies that for short-term earthquake prediction of hydrocarbon-production-induced seismicity these types of analysis could be misleading. However, further investigation of regional changes of the b-value and activity rate is necessary to exclude the fact that our observations are driven by increases in seismicity in a different part of the gas field.
Along with the decrease in seismic activity, the public commotion related to the seismic risk has also declined. Currently, public displeasure is focused mainly on the process of damage handling and compensation. However, 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 (Fig. 2C), managing both the seismicity and the public perception provides a continuing challenge.
Acknowledgements
We thank two anonymous reviewers for constructive comments that improved the content of this paper. Seismicity data used in this report were taken from the database of the Dutch Meteorological Institute (KNMI) which is accessible online through www.knmi.nl. Production information on the Groningen gas field is available online through www.nam.nl. The maps in Figure 2 were made by TNO at the request of SodM.
Open access
