Liquefaction Hazard in the Groningen Region of the Netherlands due to Induced Seismicity

The Groningen gas field is located in the northeastern region of the Netherlands and is one of the largest in the world. It has produced over $\mathrm{2,000}\mathrm{billion}{\mathrm{m}}^3$ of natural gas since the start of production in 1963. The first earthquakes linked to the gas production in the Groningen field occurred in December 1991, although earthquakes were linked to production at other gas fields in the region since 1986. To date, the largest induced earthquake linked to production at the Groningen field is the 2012 local magnitude (${\mathrm{M}}_{\mathrm{L}}$) 3.6 Huizinge event. However, the largest recorded peak ground acceleration (PGA: 0.11g) to date is associated with motions recorded during the January 8, 2018, ${\mathrm{M}}_{\mathrm{L}}$ 3.4 Zeerijp earthquake. In response to concerns about the induced earthquakes, the field operator Nederlandse Aardolie Maatschappij (NAM) is leading an effort to quantify the seismic hazard and risk resulting from the gas production operations (van Elk et al. 2017, 2019). Because of the widespread deposits of saturated sands in the region, the risk due to liquefaction triggering was evaluated as part of this effort. Although an almost negligible contributor to earthquake fatalities (e.g., Hakuno 2004; Green and Bommer 2019), liquefaction triggering and associated phenomena are important threats to the built environment and in particular to infrastructure and lifelines (e.g., Bird and Bommer 2004).

The 2017 version of the Netherlands’ National Annex to the Eurocode 8 [NPR 9998 (NPR 2017)] details requirements for assessing the impact of liquefaction and related phenomena and recommends the use of the Idriss and Boulanger (2008) simplified model (IB08) for predicting liquefaction triggering. However, as detailed in Green et al. (2019), the suitability of this model in Groningen is questionable, as is any other existing simplified model. This is because simplified models are semiempirical, with the empirical aspects having been developed primarily for tectonic earthquakes in active shallow-crustal tectonic regimes (e.g., California, Japan, and New Zealand). The seismic hazard and the geologic profiles/soil deposits in Groningen differ significantly from those where existing models were developed. Specifically, and as detailed in subsequent sections of this paper, the suitability of the depth-stress reduction factor (${\mathrm{r}}_{\mathrm{d}}$) and magnitude scaling factor (MSF) relationships inherent to existing models are uncertain for use in Groningen. Accordingly, efforts to assess the liquefaction hazard in Groningen due to induced seismicity first focused on developing a triggering model that is suitable for this task. This actually required an additional step backward to develop a revised liquefaction triggering model for tectonic earthquakes due to potential biases in the ${\mathrm{r}}_{\mathrm{d}}$ and MSF relationships inherent to existing variants of the simplified model (Lasley et al. 2016, 2017; Green et al. 2019). The Groningen-specific ${\mathrm{r}}_{\mathrm{d}}$ and MSF relationships then can be developed following the approaches used to develop the new tectonic earthquake relationships. The consistency in the approaches used to develop the revised ${\mathrm{r}}_{\mathrm{d}}$ and MSF for tectonic earthquakes and the Groningen-specific relationships allows the relationships to be interchanged within the overall liquefaction triggering evaluation procedure (NRC 2016).

To determine whether a liquefaction hazard assessment for the entire Groningen region is warranted, a comprehensive probabilistic liquefaction hazard analysis (PLHA) was performed for an area that encompassed the region of highest shaking hazard (i.e., near the village of Loppersum) and soils with the highest liquefaction susceptibility (i.e., thick Holocene sand deposits that comprise the Naaldwijk formation). Liquefaction damage-potential hazard curves are developed using a Monte Carlo method wherein probability distributions for seismic activity rates (Bourne and Oates 2017), event locations and magnitudes, and resulting ground motions are sampled such that the simulated future seismic hazard is consistent with historical seismic and reservoir-compaction datasets (Bourne et al. 2015). For each event scenario considered, the Green et al. (2019) triggering model is used in conjunction with Groningen-specific ${\mathrm{r}}_{\mathrm{d}}$ and MSF relationships to compute the factor of safety against liquefaction triggering (${\mathrm{FS}}_{\mathrm{liq}}$) as a function of depth for 95 profiles across the study area. For each of the 95 profiles, corresponding Ishihara-inspired Liquefaction Potential Index (${\mathrm{LPI}}_{\mathrm{ish}}$) (Maurer et al. 2015b) hazard curves are computed, where ${\mathrm{LPI}}_{\mathrm{ish}}$ is an index value for the severity of predicted surficial liquefaction manifestations, which has been shown to correlate with liquefaction damage potential (e.g., Iwasaki et al. 1978).

The various components of the Groningen PLHA and their interrelationships are shown in the flowchart in Fig. 1(a), with the portion of the flowchart encompassed by the dashed line being integral to the Monte Carlo simulations of the liquefaction triggering evaluations. As may be surmised from this figure, the study comprehensively and rigorously considered all significant factors that influence the regional liquefaction hazard. Also, as discussed subsequently, a similar process was implemented using IB08 for predicting liquefaction triggering [Fig. 1(b)], in lieu of the Groningen-specific triggering model developed herein, because IB08 is specified in the NPR 9998-2017 (NPR 2017).

Consistent with the requirements of NPR 9998-2017 (NPR 2017), ${\mathrm{LPI}}_{\mathrm{ish}}$ values corresponding to an annual frequency of exceedance (AFE) of $\sim 4\times {10}^{-4}$ (or a 2475-year return period) are of particular interest. However, the approach used in this study, wherein the liquefaction hazard is assessed by determining the ${\mathrm{LPI}}_{\mathrm{ish}}$ values for the specified return period, goes well beyond the requirements of NPR 9998-2017. Specifically, NPR 9998-2017 (as well as all other building codes that the authors are familiar with) allows a pseudo-probabilistic approach to be employed to assess liquefaction hazards [i.e., assessing the liquefaction triggering potential for a ground motion having a 2475-year return period, e.g., Korff et al. (2019)]. Additionally, the integration of the liquefaction hazard study within the broader regional seismic hazard study resulted in one of the most comprehensive liquefaction hazard studies performed globally to date, due to the development and use of a regional stochastic source model, ground motion prediction equation (GMPE), site response model, and geologic model, as well and the region-specific liquefaction triggering model.

Tectonic seismicity in the Netherlands is primarily associated with the Roer Valley graben in the southeast of the country (Fig. 2), where the Roer Valley graben was the source of the largest earthquake known to have occurred in the Netherlands, the moment magnitude, $\mathbf{M}$, 5.8 April 13, 1992 Roermond earthquake (Trifonov et al. 1994). In contrast, induced seismicity resulting from natural gas production mostly occurs in the northern portion of the country, including seismicity resulting from gas production in the Groningen field. This regional gas reservoir exists within the Slochteren-Rotliegend sandstone at a depth of about 3 km; it is overlain by the Zechstein salt layer, an $\sim 1\text{-}\text{km}$ thick layer of chalk, and then the North Sea Formation (Fig. 3). The seismicity resulting from production at the Groningen field from 1996 to 2019 is shown in Fig. 4. Based on projected gas production rates from 2016 to 2021, motions having a 2475-year return period are predicted to have PGAs up to approximately 0.21g (van Elk et al. 2017), with the dominant contributions from $\mathbf{M}$ 4.0–5.5 events.

The velocity model for the entire gas field from the ground surface down to the base of the North Sea Supergroup (designated henceforth as NS_B) is described in detail by Kruiver et al. (2017a, b). The NS_B horizon has an average depth of $\sim 800\mathrm{m}$ across the Groningen field and is the first elevation at which a strong and persistent velocity contrast is encountered. The velocity model from the ground surface to 50 m below the Dutch Ordnance Datum is based on a geostatistical model (i.e., the GeoTOP model) that has a $100\times 100\mathrm{m}$ spatial resolution that assigns a stratigraphic unit and a lithological class to 0.5-m thick voxels (Stafleu et al. 2011). The GeoTOP model correlates each stratigraphic lithological unit to soil parameters (mean and standard deviation). These parameters include small-strain shear wave velocity (${\mathrm{V}}_{\mathrm{S}}$), soil density, coefficients of uniformity, median grain-size diameter, cone penetration test (CPT) tip resistance (${\mathrm{q}}_{\mathrm{c}}$), and undrained shear strength. When observed, the depth dependence of ${\mathrm{V}}_{\mathrm{S}}$ is included in these correlations. Independent measurements of the near-surface ${\mathrm{V}}_{\mathrm{S}}$ profiles at accelerograph recording stations showed the GeoTOP-derived profiles to be accurate (Noorlandt et al. 2018). For depths greater than 50 m, velocities are assigned from the analysis of surface waves collected in field-wide seismic reflection surveys of the Groningen gas reservoir. These measurements extend the ${\mathrm{V}}_{\mathrm{S}}$ profile to a depth of about 120 m and are described in more detail in Kruiver et al. (2017a, b). Below this depth, measurements from sonic logs in the field are used to extend the profiles to the NS_B reference horizon (Kruiver et al. 2017a, b). The uncertainty in ${\mathrm{V}}_{\mathrm{S}}$ for depths greater than about 50 m is ignored because it was determined that these uncertainties have little impact on computed site response. The unit weights for the various stratigraphic lithological units are estimated by correlations with CPT from Lunne et al. (1997). For some of the deeper formations, the unit weights are assumed to be constant, consistent with the logs from two deep boreholes (Kruiver et al. 2017a, b). The NS_B and underlying strata have ${\mathrm{V}}_{\mathrm{S}}$ of $\sim 1.5\mathrm{km}/\mathrm{s}$ and a unit weight of $\sim 21\mathrm{kN}/{\mathrm{m}}^3$. An example of the resulting ${\mathrm{V}}_{\mathrm{S}}$ profiles from the field-wide velocity model is shown in Fig. 5.

As mentioned in the Introduction, the stress-based simplified liquefaction triggering model is central to the liquefaction hazard/risk assessment of the Groningen field. The revised model proposed by Green et al. (2019) is a modification of the Boulanger and Idriss (2014) model, which in turn is an update of the model recommended by NPR 9998-2017 (i.e., Idriss and Boulanger 2008). The modifications proposed by Green et al. (2019) were prompted by potential biases in the ${\mathrm{r}}_{\mathrm{d}}$ and MSF relationships inherent to the Boulanger and Idriss (2014) procedure, with these relationships accounting for the nonrigid response of the soil column during earthquake shaking and the influence of strong ground motion duration on liquefaction triggering. However, given the semiempirical approach used to develop the liquefaction triggering curves (or cyclic resistant ratio, CRR, curves), alternative ${\mathrm{r}}_{\mathrm{d}}$ and MSF relationships cannot simply be used in conjunction with a CRR curve that was developed using other ${\mathrm{r}}_{\mathrm{d}}$ and MSF relationships (NRC 2016). Rather, the new ${\mathrm{r}}_{\mathrm{d}}$ and MSF relationships, along with other needed relationships, must be used to analyze the compiled liquefaction/no liquefaction case histories and a new CRR curve must be regressed. Green et al. (2019) did this using the ${\mathrm{r}}_{\mathrm{d}}$ relationship proposed by Lasley et al. (2016) and an MSF that they developed using the number of equivalent cycles (${\mathrm{n}}_{\mathrm{eq}}$) correlation proposed by Lasley et al. (2017). Both of these relationships are for shallow-crustal events in active tectonic settings, consistent with the liquefaction/no liquefaction case histories analyzed.

Liquefaction hazard curves were calculated for 95 sites across the study area using a Monte Carlo approach. The probability distributions for seismic activity rates (Bourne and Oates 2017), event locations and magnitudes, and resulting ground motions were sampled such that the simulated future seismic hazard is consistent with historical seismic and reservoir-compaction datasets (Bourne et al. 2015) up to a maximum event, which its size is defined by a logic-tree (Bommer and van Elk 2017). The ground motions are predicted using the ground motion model described in Bommer et al. (2017a). The ground motion model was developed following the scheme described in Bommer et al. (2017b) and summarized previously herein is used to generate ground motions for multiple scenarios (e.g., multiple magnitude and distance combinations) at the NS_B horizon. Period-, intensity-, and scenario-dependent amplification factors (median values and standard deviations) for each of the zones in the Groningen region (e.g., Fig. 6) were used to convert the NS_B motions to surface ground motions (Rodriguez-Marek et al. 2017; Stafford et al. 2017).

For each event scenario, the Groningen-specific relationships, ${\mathrm{r}}_{\mathrm{d}\text{-}\mathrm{Gron}}$ and ${\mathrm{MSF}}_{\mathrm{Gron}}$, in conjunction with the revised CRR curve (Green et al. 2019), were used to compute the ${\mathrm{FS}}_{\mathrm{liq}}$ as a function of depth. All other relationships required to compute ${\mathrm{FS}}_{\mathrm{liq}}$ were adopted from Boulanger and Idriss (2014), consistent with the development of the revised CRR curve. The ${\mathrm{LPI}}_{\mathrm{ish}}$ framework (Maurer et al. 2015b) was used to relate the computed ${\mathrm{FS}}_{\mathrm{liq}}$ in strata at depth in a profile to the predicted severity of surficial liquefaction manifestations, which has been shown to correlate to the liquefaction damage potential for level-ground sites (e.g., Iwasaki et al. 1978; Tonkin and Taylor 2013). The ${\mathrm{LPI}}_{\mathrm{ish}}$ framework is a conceptual and mathematical merger of the Ishihara (1985) ${\mathrm{H}}_1\text{-}{\mathrm{H}}_2$ chart and liquefaction potential index (LPI) framework (Iwasaki et al. 1978), resulting in a framework that is relatively easy to implement for profiles having varying stratigraphies (where the major limitation of the Ishihara ${\mathrm{H}}_1\text{-}{\mathrm{H}}_2$ chart is that it can only be applied to profiles having relatively simple stratigraphies, as discussed in the following section) and better accounts for the thickness of the nonliquefiable crust on the severity of surficial liquefaction manifestations than the LPI framework (e.g., Green et al. 2018b). As detailed in Maurer et al. (2015b), ${\mathrm{LPI}}_{\mathrm{ish}}$ index values are functions of the thickness of the nonliquefiable crust, cumulative thickness of liquefied layers, proximity of these layers to the ground surface, and amount by which ${\mathrm{FS}}_{\mathrm{liq}}$ in each layer is less than 1.0.

The optimal ${\mathrm{LPI}}_{\mathrm{ish}}$ thresholds corresponding to different severities of surficial liquefaction manifestations are dependent on the liquefaction triggering procedure used to compute ${\mathrm{FS}}_{\mathrm{liq}}$ and the characteristics of the profile; the same is true for other liquefaction damage severity frameworks (Maurer et al. 2015a). However, without liquefaction case history data to develop Groningen-specific thresholds, the thresholds proposed by Iwasaki et al. (1978) for the LPI framework were conservatively (Maurer et al. 2015a) used in this study with the ${\mathrm{LPI}}_{\mathrm{ish}}$ framework i.e.: ${\mathrm{LPI}}_{\mathrm{ish}}<5$: no to minor surficial liquefaction manifestations; $5\le {\mathrm{LPI}}_{\mathrm{ish}}\le 15$: moderate surficial liquefaction manifestations; ${\mathrm{LPI}}_{\mathrm{ish}}>15$: severe surficial liquefaction manifestations.

The output from the liquefaction study is a set of liquefaction hazard curves for the 95 sites across the liquefaction study area, where the hazard curves show the AFE for varying ${\mathrm{LPI}}_{\mathrm{ish}}$ values for a site. Example curves computed for two sites, one in Zone 602 and one in Zandeweer, are shown in Fig. 11; Zone 602 is representative of most of the zones across the liquefaction study area, except for Zandeweer, which has the highest computed liquefaction hazard of all the zones. Consistent with the requirements of NPR 9998-2017 (NPR 2017), ${\mathrm{LPI}}_{\mathrm{ish}}$ values corresponding to an AFE of $\sim 4\times {10}^{-4}$ (or a 2475-year return period) are of most interest. However, the previous version of the NPR 9998 (NPR 2015) specified an 800-year return period for evaluating the liquefaction potential for CC1b type structures (i.e., 1, 2, or 3-story single family homes, agriculture buildings, greenhouses, and 1- or 2-story industrial buildings), which was used as the basis for a previous liquefaction hazard study of Zandeweer (Kumar 2017). Also, both the 2015 and 2017 versions of the NPR 9998 reference the IB08 liquefaction triggering model. As a result, Fig. 11 also shows the AFE corresponding to an 800-year return period and the ${\mathrm{LPI}}_{\mathrm{ish}}$ hazard curves computed using IB08. The liquefaction hazard curves for all 95 sites across the liquefaction study area, grouped by zones, are presented in Green et al. (2018a).

As may be observed from Fig. 11, use of the Groningen-specific liquefaction model in the PLHA results in a higher computed liquefaction hazard for a 2475-year return period than when IB08 is used; this was the case for all sites evaluated. Additionally, ${\mathrm{LPI}}_{\mathrm{ish}}<5$: no to minor surficial liquefaction manifestations are predicted for an 800-year return period for both sites when both the Groningen-specific and IB08 liquefaction triggering models are used in the PLHA; this again was the case for all sites evaluated. For a 2475-year return period, ${\mathrm{LPI}}_{\mathrm{ish}}<5$: no to minor surficial liquefaction manifestations are predicted for both sites when IB08 is used in the PLHA (this was the case for all sites) and for Zone 602 site when the Groningen-specific model is used in the PLHA. However, ${\mathrm{LPI}}_{\mathrm{ish}}\sim 8$ (i.e., moderate surficial liquefaction manifestations) is predicted for a 2475-year return period for the Zandeweer site when the Groningen-specific model is used in the PLHA; this is the case for 21 of the 27 sites evaluated in Zandeweer.

Fig. 12 shows liquefaction hazard maps for a 2475-year return period computed using the Groningen-specific model in the PLHA for the entire liquefaction study area and for Zandeweer. The only sites in the liquefaction study area that are predicted to have moderate surficial liquefaction manifestations are all located in Zandeweer [Fig. 12(b)]; all other sites in the study area have ${\mathrm{LPI}}_{\mathrm{ish}}<5$.

Although an almost negligible contributor to earthquake fatalities, liquefaction triggering is an important threat to the built environment and in particular to infrastructure and lifelines and is thus being included as part of the seismic hazard and risk study of the Groningen region being directed by NAM. Due to the unique characteristics of both the seismic hazard and the geologic profiles/soil deposits in Groningen, direct application of existing liquefaction triggering models in the study was deemed inappropriate and was shown to be so by comparing ${\mathrm{r}}_{\mathrm{d}}$ and MSF relationships inherent to existing triggering models with Groningen-specific relationships developed as part of this study. A PLHA was performed using a Monte Carlo method wherein probability distributions for seismic activity rates, event locations and magnitudes, and resulting ground motions are sampled such that the simulated future seismic hazard is consistent with historical seismic and reservoir-compaction datasets up to a maximum event, the size of which is defined by a logic-tree distribution. The differences between the Groningen-specific and the NPR 9998-2017 specified IB08 liquefaction triggering models manifest in the computed liquefaction hazard curves (i.e., AFE versus ${\mathrm{LPI}}_{\mathrm{ish}}$) for sites across the study area, with the predicted liquefaction hazard generally being higher when the Groningen-specific liquefaction triggering model is used in the PLHA.

Consistent with the requirements of NPR 9998-2017, ${\mathrm{LPI}}_{\mathrm{ish}}$ values corresponding to an AFE of $\sim 4\times {10}^{-4}$ (or a 2475-year return period) are of particular interest. However, assessing the liquefaction hazard by determining the ${\mathrm{LPI}}_{\mathrm{ish}}$ values for the specified return period goes well beyond the requirements of NPR 9998-2017 (as well as all other building codes that the authors are familiar with), which allows a pseudo-probabilistic approach to be employed (i.e., assessing the liquefaction triggering potential for a ground motion having a 2475-year return period). Additionally, the integration of the liquefaction hazard study within the broader regional seismic hazard study resulted in one of the most comprehensive PLHA performed to date globally due to the development and use of a regional stochastic source model, GMPE, site response model, and geologic model, as well as the region-specific liquefaction model (refer to Fig. 1). The ${\mathrm{LPI}}_{\mathrm{ish}}$ values for the vast majority of the sites across the study area are less than 5, indicating no to minor surficial liquefaction manifestations. The only sites within the study area that had ${\mathrm{LPI}}_{\mathrm{ish}}$ values greater than 5, which is the threshold between no to minor surficial liquefaction manifestations and moderate surficial liquefaction manifestations, were in Zandeweer, with only some of the sites in Zandeweer exceeding this threshold value. No sites across the study are predicted to have severe surficial liquefaction manifestations.

Finally, mention is needed regarding two phenomena that were not considered in the liquefaction study presented in this paper: the influence of sand age on liquefaction triggering resistance (i.e., aging effects) and the influence of thin layer effects on measured CPT tip resistance (i.e., thin layer effects). Investigations into both of these phenomena were being performed in parallel with the work presented in this paper, with the efforts lead by colleagues at Deltares. The Deltares’ studies were not far enough along to be incorporated into the analyses presented herein. However, it should be noted that ignoring both aging and thin layer effects likely results in an overprediction of liquefaction hazard (i.e., the hazard estimates presented herein may be considered conservative), but the extent of the possible overprediction is unknown at this time.

This research was partially funded by Nederlandse Aardolie Maatschappij B.V. (NAM) and the National Science Foundation (NSF) Grants CMMI-1435494, CMMI-1724575, CMMI-1825189, and CMMI-1937984. This support is gratefully acknowledged. This study has also significantly benefited from enlightening discussions with colleagues at Shell, Deltares, Arup, Fugro, Beca, and on the Nederlands Normalisatie Instituut (NEN) Liquefaction Task Force, which is coordinated by Dr. Mandy Korff. However, any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the NSF, NAM, or others that are acknowledged.