Introduction
Since 1986, earthquakes have been recorded near producing gas fields in the provinces of Groningen, Drenthe and Noord-Holland and in northern Germany. In the early 1990s, a multidisciplinary study of these events concluded that the observed earthquakes were of non-tectonic origin and induced by reservoir depletion (i.e. gas production). Following up on this conclusion, a shallow borehole seismometer network was installed in the Groningen area. The network was designed to detect earthquakes, pinpoint their locations and quantify their magnitudes, and has been operational since 1995. Additional accelerometers were installed in areas with highest earthquake density. This seismic monitoring network showed a gradual increase in seismic activity, particularly after 2003.
A renewed focus on the issue of seismicity induced by gas production in Groningen started in 2012 following the earthquake near Huizinge (16 August 2012) with magnitude M L = 3.6. This earthquake was felt as more intense, with a longer duration, and caused significantly more building damage than previous earthquakes. People living in the Groningen field area have been confronted with increasing intensity of the effects of induced earthquakes (Tomale, Reference Tomale2015), and the Huizinge event has left an especially deep imprint. A programme to study seismicity induced by gas production from the Groningen field and the associated effects on buildings and people living above the field was initiated (NAM et al., Reference van Elk and Doornhof2012, Reference van Elk and Doornhof2014, 2016 Reference van Elk and Doornhofa,Reference van Elk and Doornhofb). This required studies into domains where NAM did not have sufficient expertise. Therefore, the plan mobilised support from Dutch and international universities, knowledge institutes and laboratories and sought the assistance and advice of external recognised experts.
Research into induced seismicity in Groningen
The seismicity in Groningen differs from typical tectonic earthquakes in three aspects: (1) it is induced by the withdrawal of gas, (2) the shallow subsurface in the Groningen area consists of thick layers of very soft deposits and (3) the buildings in the Groningen area have been designed and constructed without consideration for the horizontal loads typically experienced during an earthquake. Acquisition of data specific to the Groningen area is therefore essential, and data gathering in the field and laboratory experiments are an important part of the research programme. The existing geophone and accelerometer network has been extended to cover the full field, and geophones have been placed at the reservoir level in deep wells. The programme was first described in the ‘Study and Data Acquisition Plan’ of November 2012 (NAM et al., Reference van Elk and Doornhof2012), which has since been updated several times (NAM et al., 2014, Reference Graziotti, Penna and Magenes2016 Reference van Elk and Doornhofa,Reference van Elk and Doornhofb).
The main objectives of the Study and Data Acquisition Plan are:
-
1. To understand the impact of the earthquake hazard on buildings and other structures and the subsequent impact on safety of the community;
-
2. To perform a fully integrated probabilistic hazard and risk assessment for the Groningen region, with all uncertainties fully and consistently recognised and quantified;
-
3. To identify, evaluate and develop mitigation options to reduce safety risk:
-
a. Production measures, i.e. changes in the production from the field
-
b. An optimised Structural Safety Upgrading programme
-
c. Identify buildings and/or building elements that pose a safety risk
-
d. Establish optimal structural upgrading methodologies;
-
-
4. Measures for industry and infrastructure.
Other important objectives are to:
-
5. Discuss the merits of alternative scientific views, and initiate additional studies and/or data acquisition to promote consensus amongst the knowledge institutes;
-
6. Monitor compaction, subsidence and seismicity;
-
7. Continuously improve our understanding of the physical mechanisms leading to induced seismicity and the resulting hazard;
-
8. Reduce the uncertainty in the hazard and risk assessment.
To achieve these objectives the research programme covered the full causal chain (Fig. 1) from the production of gas to the effects of the earthquakes at the surface, damage to buildings and risk for people. Updates on the impact of the data acquisition and studies of the assessment of hazard and risk were prepared regularly during the execution of the research programme. This also allowed a review of the programme, and identification of additional study requirements. In practice, almost every 6 months an updated hazard and risk assessment was published (NAM, 2013, 2014, Reference Bommer, Crowley and Pinho2015 a,b, 2016). The aim of the programme was to achieve a well-calibrated hazard and risk assessment for the update of the Winningsplan in July 2016. The intermediary assessments aimed to be conservative in the sense that future updates were more likely to give a reduced assessment of hazard and risk rather than an upwards adjustment (for the same level of gas production).
Fig. 1. The causal chain from gas production to safety of people was depicted as a bridge in the ‘Study and Data Acquisition Plan’ (NAM et al., 2014).
The model used to assess hazard and risk resulting from induced seismicity in Groningen, is based on the probabilistic seismic hazard and risk assessment (PSHA) methodology introduced by Cornell (Reference Cornell1968), which has since become the internationally recognised standard method for seismic risk assessments. This probabilistic risk assessment method fits well with the methods used to determine other risks, such as those resulting from floods, plane crashes around airports and industrial accidents. The Monte Carlo method is used to implement the model.
The PSHA methodology has been further developed since its inception in the 1960s (Calvi et al., Reference Calvi, Pinho, Magenes, Bommer, Restrepo-Vèlez and Crowley2006; Scawthorn, Reference Scawthorn, Bostrom, French and Gottlieb2007, Reference Scawthorn2008; Bommer et al., Reference Bommer, Crowley and Pinho2015), and forms the basis of seismic risk assessments developed in the USA under the auspices of regulators like the National Research Council (NRC) and US Geological Survey (USGS). Also, the HAZUS risk assessment method and tools of the Federal Emergency Management Agency (FEMA) are based on a similar probabilistic method for seismic risk assessment. The method allows the assessed risk to be compared to safety norms set by society. For induced seismicity in the Netherlands, the norms were set by the ‘Committee Meijdam’ (Commissie-Meijdam, 2015).
Although the methodology used for the assessment of hazard and risk (PSHA) and the implementation (Monte Carlo) are based on standardised accepted and established processes, the hazard and risk assessment for induced seismicity in Groningen also contains innovations and adaptations to the local situation. Central to the probabilistic assessment of hazard and risk is the identification, quantification and consistent treatment of all uncertainties throughout the calculation chain of the Monte Carlo procedure. The data acquisition activities are critical for the assessment of the uncertainties and reduction of the epistemic uncertainty.
Detailed models were built of the Groningen gas reservoir (de Jager & Visser, Reference De Jager and Visser2017; Visser & Solano Viota, Reference Visser and Solano Viota2017) and these were calibrated with available and newly acquired pressure data, observations of water ingress into the gas reservoir and subsidence measurements (Van Eijs et al., Reference Van Eijs, van der Wal and Doornhof2017; Van Oeveren et al., Reference Van Oeveren, Valvatne, Geurtsen and van Elk2017;). Based on geomechanical models of fault movement and the historical earthquake catalogue, seismological models were built (Bommer & Van Elk, Reference Bommer and van Elk2017; Dempsey & Suckale, Reference Dempsey and Suckale2017; Zöller & Holschneider, Reference Zöller and Holschneider2017) for induced seismicity in response to gas production scenarios. An expert panel using the SSHAC (Senior Seismic Hazard Analysis Committee) method estimated the maximum magnitude of an induced earthquake in Groningen (Harris & Bourne, Reference Harris and Bourne2017; Zöller & Holschneider, Reference Zöller and Holschneider2017).
Additionally, geomechanical studies (Burnett et al., Reference Burnett, Gans, Gist, Terrell, Reilly, Tomic, Anderson, Tan, Dedontney and Pais2016; Lele et al., Reference Lele, Hsu, Garzon, DeDontney, Searles, Gist, Sanz, Biediger and Dale2016; Buijze et al., Reference Buijze, van den Bogert, Wassing, Orlic and ten Veen2017;) and experiments on core material acquired in the Zeerijp-3A well (Hunfeld et al., 2017) have been carried out to gain a better understanding of the mechanism causing the earthquake ruptures.
Based on detailed models of the shallow subsurface (Kruiver et al., Reference Kruiver, van Dedem, Romijn, de Lange, Korff, Stafleu, Gunnink, Rodriguez-Marek, Bommer, van Elk and Doornhof2017 Reference Kruiver, van Dedem, Romijn, de Lange, Korff, Stafleu, Gunnink, Rodriguez-Marek, Bommer, van Elk and Doornhofa,Reference Kruiver, Wiersma, Kloosterman, de Lange, Korff, Stafleu, Busscher, Harting, Gunnink, Green, van Elk and Doornhofb; Rodriguez-Marek et al., Reference Rodriguez-Marek, Kruiver, Meijers, Bommer, Dost, van Elk and Doornhof2017; R.P. Noorlandt et al., unpublished work), ground motion prediction models have been developed for application to the Groningen earthquakes (Bommer et al., Reference Bommer, Dost, Edwards, Stafford, van Elk, Doornhof and Ntinalexis2016, Reference Bommer and van Elk2017 Reference Bommer, Dost, Edwards, Kruiver, Ntinalexis, Rodriguez-Marek, Stafford and van Elka,Reference Bommer, Stafford, Edwards, Dost, van Dedem, Rodriguez-Marek, Kruiver, van Elk, Doornhof and Ntinalexisb). Besides the hazard posed by ground movement, the potential for liquefaction is also studied (Lasley et al., Reference Lasley, Green and Rodriguez-Marek2016, Reference Lasley, Green and Rodriguez-Marek2017).
To be able to assess the response of buildings to the induced earthquakes in Groningen, the study programme also included a civil engineering component. The properties of building material and building elements like single leaf walls and cavity walls have been tested in laboratory experiments (Graziotti et al., Reference Graziotti, Penna and Magenes2016a,Reference Graziotti, Tomassetti, Penna and Magenesb; Malome et al., Reference Malome, Pinho and Pennein press). Several buildings have been tested at the shake tables of EUcentre (Pavia, Italy) (Graziotti et al., 2017) and LNEC (Lisbon, Portugal). Results of these experiments were used to calibrate computer models for different buildings and prepare fragility functions for different building typologies encountered in the Groningen area (Graziotti et al., Reference Graziotti, Tomassetti, Kallioras, Penna and Magenes2017).
To ensure the research is carried out as objectively as possible, it is subjected to rigorous assurance. This process consists of seven layers, the most important of which are the assurance by the Scientific Advisory Committee (SAC) of the Minister of Economic Affairs and independent international experts (NAM et al., 2016 Reference van Elk and Doornhofa,Reference van Elk and Doornhofb). Also transparency of the studies is important to allow external experts to form an opinion of the quality and objectivity of the studies. Study reports are published on a publicly accessible dedicated website (www.nam.nl/feiten-en-cijfers/onderzoeksrapporten), and seismic data are shared with interested academics. Currently some 70,000 downloads (as at 31 July 2017) have been made from the website, while data have been shared with 10 universities from all over the world. After publication on the website, the research is also published in peer-reviewed journals.
Probabilistic assessment of the hazard resulting from induced seismicity
Hazard metric
Because peak ground acceleration (PGA) is a widely used metric for ground-shaking intensity, it was chosen as the most appropriate hazard metric for this seismic hazard assessment. When extending the assessment to encompass risk (i.e. the response of buildings to ground shaking), spectral acceleration (SA) will be used; this takes into account the response period of the building being considered. Figure 2 shows the measured acceleration near the epicentre during the Huizinge earthquake of 16 August 2012. In addition to the PGA values the duration of the events has also been incorporated in the ground motion prediction methodology, as this is also important for the seismic risk.
Fig. 2. Accelerogram of the earthquake near Huizinge recorded on 16 August 2012 by the accelerometer located near Westeremden (near the epicentre).
Peak ground acceleration and hazard maps
For the probabilistic description of the ground accelerations (PGA, or generalised to peak spectral acceleration, PSA), a hazard map is used. On this map for each location the acceleration is plotted that could occur, with a prescribed annualised probability of exceedance (exceedance level), during a prescribed analysis period. Hazard levels are shown using a gradual colour scale.
The hazard maps shown in the section were constructed according to the following procedure. Each location in the analysis area during the analysis period is subjected to ground motion accelerations resulting from induced earthquakes. At some locations, e.g. near Loppersum, the chance of exceeding a given peak ground acceleration threshold is higher than at the periphery of the field. Equally, at any one location, the chance of exceeding some value of peak ground acceleration decreases with increasing peak ground acceleration. An example of a set of hazard curves is shown for a number of locations in Figure 3. Each declining line indicates the hazard curve for a single location in the field.
Fig. 3. An example of a set of hazard curves showing average annual exceedance rate for peak ground acceleration at different locations in the field. Each line corresponds to a location in the field. The bold line indicates the maximum PGA anywhere within the field for a given exceedance level (bounding envelope). The red line indicates that for an exceedance level of 0.2%a−1 the highest PGA in the field is 0.21g.
To prepare a hazard map, an exceedance level needs to be chosen. This is not a purely technical choice. Inspired by Eurocode 8Footnote 1 , part of the current technical standards for structural design in Europe, it has become common practice to prepare hazard maps for an exceedance level of 0.2%a−1. This exceedance level is equivalent to a 475-year return period for stationary seismicity. The same exceedance level is also used by KNMI for their hazard maps. Hazard maps can be made for different production scenarios. This return period is only used for representation of the hazard. The actual assessment of hazard and risk is not affected by this choice.
The hazard assessment critically depends on the modelling of reservoir pressure in response to gas withdrawal from the reservoir, the resulting reservoir compaction, assessed seismicity and the method to predict the motion at surface. During the study programme, progress was made in each of the areas of study. Together with the impact of lower annual gas production levels, this is reflected in the hazard maps that have been published while the studies progressed (Fig. 4).
Fig. 4. Hazard maps published by NAM from 2013 to 2016, illustrating the development of the hazard maps while the research programme progressed.
The hazard map included in Winningsplan 2013 (NAM, 2013) used a method based on tectonic earthquakes in southern Europe to predict the ground motion. This map also assumed a soft soil to be present throughout the Groningen area. The hazard map published in November 2015 included the impact of the detailed mapping of the soil geology. Instead of smooth concentric hazard contours, the new contours now reflect variation in the soil properties. Table 1 gives an overview of the most important parameters for the hazard maps. Furthermore, with version 2 also assessment of the response spectra and seismic event duration was added. This allowed for more accurate estimation of the building response to a seismic event.
Table 1. Most relevant background data for the hazard maps published by NAM from 2013 to 2016.
The impact of lateral heterogeneity in the composition of the shallow subsurface can clearly be seen in the hazard map. The hazard map for the scenario of an annual field production of 27bcma−1, and the current offtake distribution is shown in Figure 5.
Fig. 5. PGA hazard maps. Period: 2016–2021; production: 27 bcma−1; compaction: inversion; activity rate model: Version V2, 3.5≤M≤6.5; metric: 0.2%a−1 chance of exceedance (10% chance in 50 years). Mean hazard from logic tree. The maximum PGA in this map is 0.21g.
Disaggregation of seismic hazard
Disaggregation of the Monte Carlo results for the hazard assessment provides insight into which earthquakes have most impact on the hazard assessment at a given location. Examples of a disaggregation are shown in Figure 6 (for the Loppersum area) and Figure 7 (for the city of Groningen).
Fig. 6. (A) Occurrence rates for peak spectral acceleration at 0.01s as a function of magnitude, distance and GMPE epsilon, ε, for a single surface location directly above the region of maximum reservoir compaction. Grey denotes no occurrence in any of the simulations. (B) The fractional contribution to the ground motion with a 0.2% annual probability of exceedance from January 2016 to January 2021.
Fig. 7. As previous figure, except for a surface location in the centre of Groningen city.
The disaggregation of the hazard for Loppersum shows that the largest contribution to the hazard is from earthquakes within the Loppersum area (small distance of less than 5km away, with a magnitude ranging from 4 to 5). In contrast, the largest contribution to the hazard in the Groningen city is from earthquakes with an epicentre approximately 10km away from the city (towards the Loppersum area). To cause significant ground acceleration in the city of Groningen, these earthquakes located further away require a larger magnitude or GMPE epsilon to cause similar ground motions.
Sites with poor response (ε≥0) contribute most to the hazard for the Loppersum area and the city of Groningen, with the largest contribution for ε=1.
Probabilistic assessment of the risk resulting from induced seismicity
Risk metrics
The results from the probabilistic hazard and risk analysis (PHRA) are summarised via risk metrics which are related to the annualised probability of fatality for an individual person or for groups of people, taken as an average across the forecast period of the PHRA. These risk metrics – ‘Inside Local Personal Risk’, ‘Object-bound Individual Risk’, ‘Number of People at Risk’, ‘Community Risk, and ‘Social Risk’ – are defined below. ‘Inside Local Personal Risk’ and ‘Object-bound Individual Risk’ are individual risk metrics (related to probability of fatality for an individual), whereas the remaining three metrics are measures of aggregate risk (related to probability of fatality for multiple people or for groups of people).
When measuring risk, it is important to select a risk metric that is appropriate given the purpose of the risk measurement. In many cases there is more than one option available. An advisory committee, Commissie Meijdam, was established in early 2015 to advise on risk policy related to Groningen earthquakes, including the selection of risk metrics. In December 2015 the Commissie-Meijdam shared its third and final advice with the Minister of Economic Affairs. The selection of risk metrics for this PHRA reflects the final advice published by Commissie-Meijdam (2015).
Risk metric
Both individual (related to probability of fatality for an individual) and aggregated (related to probability of fatality for multiple people or for groups of people) metrics have been used to describe the risk resulting from induced seismicity in Groningen. Individual risk metrics include ‘Inside Local Personal Risk’ and ‘Object-bound Individual Risk’. Aggregated risk metrics include ‘Community Risk, ‘Group Risk’ and ‘Social Risk’. In this paper we will focus on Inside Local Personal Risk (ILPR); the annual probability of fatality for a fictional person, who is continuously present without protection, inside a building.
Risk assessment for Inside Local Personal Risk (ILPR)
With knowledge of the presence of people in the buildings, the number of people exceeding an ILPR can be estimated. The solid black line in Figure 8 shows the number of people exposed to a certain level of local personal risk. During this 5-year period, there are no buildings where the occupants are exposed to a mean local personal risk larger than 10−4a−1. Occupants of some 100 buildings are exposed to a mean local personal risk exceeding 10−5a−1 in the period 2016–2021. Over the period 2021–2026 this increases with some 100 additional buildings. The uncertainty in the risk assessment indicates this could be as many as 1,000 buildings. As risk in this context is often plotted as a logarithmic quantity, the mean log local personal risk is also shown. The shaded grey areas indicate the norm set by the Commissie-Meijdam (2015).
Fig. 8. Number of buildings and people exceeding a given ILPR shown on (top) a linear scale and (below) a log scale for the 27bcm production scenario and the 2016–2021 assessment period. The grey areas indicate the norm advised by the Commissie-Meijdam.
The spatial distribution of buildings within given ranges of ILPR is shown in Figure 9 for an optimized production scenario of 27 billion m3 (bcm) a−1. Each of the approximately 160,000 occupied buildings within the exposure area is represented by a single dot. These are plotted in order of increasing risk so that the largest risks plot on top. Grey dots denote risks smaller than 10−6a−1.
Fig. 9. ILPR for every individual building within four equal risk bands from 10−7 to 103a−1 for the 5-year assessment period 2016–2021 under the 27bcm production scenario.
When comparing these numbers with the norms advised by the Commissie-Meijdam, the second map is relevant. It shows the location of the buildings that do not meet this norm in the 27bcma−1 production scenario. As this is the result of a probabilistic assessment, it must be validated through inspections of buildings. These estimates do not include risk from non-structural elements, which has been assessed through a separate methodology and is described in a separate risk metric for falling objects.
Disaggregation of risk (ILPR)
A disaggregation of contributions to the base-case ILPR was performed for magnitude, distance from the epicentre, the ground motion variability measure ε, and spectral acceleration causing building collapse. As an example, Figure 10 shows the results for the residential apartment buildings of unreinforced masonry with silica–calcium load-bearing walls (type B) in the Loppersum area (typology RESA-URM-B).
Fig. 10. The fractional contribution to ILPR for the RESA-URM-B building typology at two locations: Loppersum (top row) and Groningen city centre (bottom row). This result was obtained for the 2016–2021 assessment period under the 33bcma−1 production scenario and the base-case scenario of the risk logic tree. Fluctuations between neighbouring points are due to finite sampling effects of the Monte Carlo procedure; nonetheless the underlying trends are clear.
As for hazard, earthquakes in the Loppersum area (i.e. at epicentral distances less than 5km) contribute most to the risk for this area. For Groningen city, earthquakes at an epicentral distance of 10km (i.e. in the Loppersum area) are the most important contribution to the risk.
Structural upgrading as a mitigating option
The estimate of the number of buildings exceeding an ILPR of 10−5 does not directly translate into an estimate of the scope of the building-strengthening effort. Reasons for a larger scope of the structural upgrading effort are:
-
■ Efficiency of identifying buildings with ILPR >10−5 has not yet been proven. This is a probabilistic assessment and does not directly indicate every individual building that needs to be included in the structural upgrading plan.
-
■ Remaining uncertainty in hazard and risk assessment. Future updates of the hazard and risk assessment could result in a different mean value of the risk, for instance based on the results of shake table tests of a new typology of building, or where the small number of buildings subject to ‘special circumstances’ (e.g. buildings located on ‘Wierden’) is taken into account.
-
■ Differences between the hazard and risk assessment and NEN-NPR building code. Ultimately the structural upgrading scope will be based on the NEN-NPR building code.
Conclusions
As part of the research programme started by NAM in 2013, assessments of hazard and risk resulting from the exposure to induced earthquakes have been prepared. While more data were collected in the Groningen area as the studies progressed, the models could be refined and the uncertainty in these assessments could be reduced. These assessments have been publicly shared and supported both Winningsplan 2013 and Winningsplan 2016.
Acknowledgements
Many people have contributed to this work, either directly or through very constructive and helpful discussions. In particular, the comments from the Scientific Advisory Committee led by Lucia van Geuns have been instrumental in improving this work. We also wish to thank the anonymous reviewers for helpful comments on the original version of this paper, and guest editor Karin van Thienen-Visser for assistance with processing this contribution.
Open access
