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

The operator of the Groningen gas field is leading an effort to quantify the seismic hazard and risk of the region due to induced earthquakes, including overseeing one of the most comprehensive liquefaction hazard studies performed globally to date. Due to the unique characteristics of the seismic hazard and the geologic deposits in Groningen, efforts first focused on developing relationships for a Groningenspecific liquefaction triggering model. The liquefaction hazard was then assessed using a Monte Carlo method, wherein a range of credible event scenarios were considered in computing liquefaction damage-potential hazard curves. This effort entailed the use of a regional stochastic seismic source model, ground motion prediction equation, site response model, and geologic model that were developed as part of the broader regional seismic hazard assessment. No to minor surficial liquefaction manifestations are predicted for most sites across the study area for a 2475-year return period. The only sites where moderate surficial liquefaction manifestations are predicted are in the town of Zandeweer, with only some of the sites in the town being predicted to experience this severity of liquefaction for this return period. DOI: 10.1061/(ASCE)GT.1943-5606.0002286. This work is made available under the terms of the Creative Commons Attribution 4.0 International license, https://creativecommons.org/licenses/by/4.0/.


Introduction
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 2,000 billion 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 (M 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, M 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 . 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;, 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 , 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 (r d ) and magnitude scaling factor (MSF) relationships inherent to existing models are uncertain for 1 Professor, Dept. of Civil and Environmental Engineering, Virginia Tech, Blacksburg, VA 24061 (corresponding author). ORCID: https:// orcid.org/0000-0002-5648-2331. Email: rugreen@vt.edu 2 Senior Research Investigator, Dept. of Civil and Environmental Engineering, Imperial College London, London SW7 2AZ, UK. ORCID: https://orcid.org/0000-0002-9709-5223 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 r d and MSF relationships inherent to existing variants of the simplified model (Lasley et al. 2016(Lasley et al. , 2017. The Groningen-specific r 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 r 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 , 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  triggering model is used in conjunction with Groningen-specific r d and MSF relationships to compute the factor of safety against liquefaction triggering (FS 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 (LPI ish ) (Maurer et al. 2015b) hazard curves are computed, where LPI 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 Groningenspecific triggering model developed herein, because IB08 is specified in the NPR 9998-20179998- (NPR 2017. Consistent with the requirements of NPR 9998-2017 (NPR 2017), LPI ish values corresponding to an annual frequency of exceedance (AFE) of ∼4 × 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 LPI 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.

Background
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, 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 ∼1-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 , with the dominant contributions from 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 ∼800 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 × 100 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 (V S ), soil density, coefficients of uniformity, median grain-size diameter, cone penetration test (CPT) tip resistance (q c ), and undrained shear strength. When observed, the depth dependence of V S is included in these correlations. Independent measurements of the near-surface V 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 V 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 V 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 V S of ∼1.5 km=s and a unit weight of ∼21 kN=m 3 . An example of the resulting V S profiles from the field-wide velocity model is shown in Fig. 5.   (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 r d and MSF relationships cannot simply be used in conjunction with a CRR curve that was developed using other r d and MSF relationships (NRC 2016). Rather, the new r 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.  did this using the r d relationship proposed by Lasley et al. (2016) and an MSF that they developed using the number of equivalent cycles (n 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.
Groningen-Specific r d , n eq , and MSF Relationships

Geologic Profiles
As mentioned previously, Groningen-specific r d , n eq , and MSF needed to be developed following the approaches used by Lasley et al. (2016Lasley et al. ( , 2017 and  to develop the revised r d , n eq , and MSF for tectonic earthquakes. This allows the Groningenspecific relationships to be used in conjunction with the CRR curve developed by . Accordingly, a series of site response analyses needed to be performed using region-specific ground motions and geologic profiles. The GeoTOP model was used to develop profiles for use in the site response analyses.
The NS_B was chosen as the reference rock horizon and was treated as the top of the elastic half-space in the site response analyses performed in this study . Although another velocity contrast is encountered at a depth of ∼400 m at the Brussels Sands formation, a velocity reversal occurs at a depth of ∼500 m that is inconsistent with the properties of an elastic half-space and thus prevents the top of the Brussels Sands formation from being used as the reference horizon. Moreover, the Brussels Sands formation is not consistently mapped across the entire field; in contrast, the NS_B is well defined over the entire study area.
The liquefaction study area shown in Fig. 6 crosses over several zones used to develop site amplification factors for the region (Rodriguez-Marek et al. 2017); for the most part, the site amplification zones coincide with the geological zonation presented in Kruiver et al. (2017a). From NW to SE, the liquefaction study area crosses over zones 821, 801, 603, 604, 1001, 602, 1032, and 2001. Given that r d , n eq , and MSF relationships are functions of the response characteristics of a geologic profile, separate relationships were developed for each of the site amplification zones within the liquefaction study area. Additionally, previous studies raised concerns about the liquefaction hazard in the downtown region of Zandeweer (Kumar 2017), which lies within the liquefaction study area. As a result, this region was treated as a separate zone.
In selecting sites used to develop the Groningen-zone-specific r d , n eq , and MSF relationships, preference was given to sites that were characterized by CPT; the locations of the CPT soundings within the study area that were available at the start of this study are shown in Fig. 6(b). For each CPT, the corresponding profile realization (i.e., GeoTOP voxel stack) in which the sounding was performed was identified. In several cases, two or more soundings were performed within the same 100 × 100 m GeoTOP voxel stack. If the number of unique GeoTOP voxel stacks in which the CPT soundings were performed in a zone was less than ten, then additional GeoTOP voxel stacks were chosen to ensure that at least ten unique GeoTOP voxel stacks per zone were used in this study. The locations of these additional GeoTOP voxel stacks were selected to ensure a relatively uniform spatial distribution across a zone. Also, if the number of unique GeoTOP voxel stacks in which the CPT soundings were performed in a zone was at least ten but the spatial distribution of these stacks was not relatively uniform across the zone, three additional GeoTOP voxel stacks were chosen for that zone. The locations of the additional GeoTOP voxel stack profiles used in the site response analyses are shown in Fig. 6(b) and listed as "Extra GeoTOP locations" in the figure legend. The numbers of profiles per zone are listed in Table 1, with a total of 110 profiles across the liquefaction study area used to develop the zone-specific r d , n eq , and MSF relationships.

Ground Motions
The ground motions at the NS_B reference horizon used in the site response analyses to develop the Groningen-zone-specific r d , n eq , and MSF relationships were generated using finite-fault, stochastic simulation using the EXSIM code (Motazedian and Aktinson 2005;Boore 2009). Each of the distributed subfaults in this technique is assumed to be a point source (effectively a small magnitude earthquake) and can be characterized using the seismological parameters observed in events recorded in the Groningen gas field. This process is detailed by Bommer et al. (2017b) and Edwards et al. (2019) and briefly summarized in the following section.
The first stage in calibrating the source model is to transform the surface recordings to the NS_B reference horizon, central to which is the velocity model from the NS_B horizon to the ground surface discussed previously. The Fourier amplitude spectra (FAS) from the surface and from 200-m boreholes were transformed to the NS_B horizon using a one-dimensional transfer function, as implemented in the site response software STRATA (Kottke and Rathje 2008). Because the motions recorded to date are very weak (i.e., M L 2.5 to 3.6 and distances from 0 to 20 km), the near-surface layers were assumed to have responded essentially linearly to the excitations; therefore, the deconvolution was made assuming linear site response. Linear amplification factors were also calculated for the response spectra using the random vibration theory (RVT) procedure implemented in STRATA. The factors were calculated for the same V S , damping, and unit weight profiles, with input motions at the NS_B obtained from simulations, using an earlier version of the model. The amplification factors for the response spectra were found to be scenario-dependent (Stafford et al. 2017).
The FAS at the NS_B horizon were then inverted, following the approach of Edwards et al. (2019), for source, path, and site  parameters. Three parameters were then fitted to each FAS: the event-specific source corner-frequency (f 0 ) of a parametric source spectrum (e.g., Brune 1970;Boatwright 1978), record-specific signal moment (long-period spectral displacement plateau), and highfrequency attenuation term, kappa (κ). As a result of the relatively small dataset of recorded motions available for these inversions, some elements were constrained independently. First, κ was estimated from individual FAS plotted on log-linear axes following Anderson and Hough (1984). Second, the form of the geometric spreading was obtained from full waveform finite difference simulations performed using a 3D velocity model for the field. The inversions were then used to estimate the stress parameter (Δσ), decay rates in each segment, quality factor (Q) value, and amplification factor at the NS_B (which was found to be close to 1). For the inversions, the use of both the Brune (1970) and Boatwright (1978) spectra was explored, and because neither performed consistently better, the former was used.
Based on the geometric spreading model and the results of the inversions of the FAS, many combinations of Δσ, Q, and site kappa (κ 0 ) were explored to identify the combination of parameters that best fit the response spectral ordinates at the NS_B horizon. The duration model used to link the FAS and response spectral ordinates was an empirically-based Groningen-specific model (Bommer et al. 2016). The model that best fit the motions at the NS_B horizon in this case had a Δσ of 60 bar, frequency independent Q of 200, and κ 0 of 0.015 s, which were then used in forward simulations.
The software EXSIM (Boore 2009, based on Motazedian and Aktinson 2005) was used in conjunction with the Groningenspecific stochastic point source model parameters to generate motions at the NS_B for magnitudes ranging from M 3.5 to 7.0, with ΔM ¼ 0.5, and horizontal distances to the fault center (R) ranging from 0.1 to 60 km, with Δ logðRÞ ¼ 0.2, resulting in a total of 120 magnitude-distance combinations. Each M-R scenario was recorded at eight azimuths radially around the center of the strike of the finite fault, leading to 960 motions. Fault dimensions were based on Wells and Coppersmith (1994) for M > 5, otherwise, Brune (1970), with each event sampling one realization of the distribution (with a standard deviation of 0.15 log-units) of possible fault dimensions. Hypocenters were randomly located along strike, but always at z ¼ 3 km (the reservoir depth), with propagation downwards (at 0.8β with 50% subfault pulsing) until the seismogenic depth of 13 km. For larger events, fault ruptures then grow horizontally, bounded by the reservoir and the seismogenic depth. Subfault durations were sampled from the Bommer et al. (2016) duration model, accounting for between-and within-event variability. The maximum magnitude of M 7.0 was determined by an expert panel . In extrapolating to magnitudes so much larger than the largest recorded event (M 3.6), the inevitable epistemic uncertainty in the generated motions for larger magnitude events needed to be considered. This was performed by both introducing magnitude-dependence of Δσ into the model and also creating alternative (higher and lower) models, as had been done previously (Bommer et al. 2017b). The weights assigned to these branches are 0.1 for the lower branch since it is unclear whether low stress drop values would persist at larger magnitudes and 0.3 each to the two central and upper models. In total, 3840 motions were generated and used in the site response analyses to develop the Groningen-zone-specific r d , n eq , and MSF relationships.

Regression Analysis
In total, 422,400 site response analyses were performed to develop Groningen-zone-specific r d , n eq , and MSF relationships (110 profiles × 3,840 motions = 422,400 analyses). The site response analyses were performed using the equivalent linear code, ShakeVT2 (Thum et al. 2019), which outputs r d and n eq values for the liquefiable layers, directly, for each analysis as a function of depth in the profile. The depth range considered for the regression analyses was restricted to 20 m and shallower. The reason for this was twofold. First, the scaling of r d was very stable over this depth range, and including more data for deeper depths only acted to negatively influence the regression fit over the 0-20 m depth range. Second, and more importantly, liquefaction is a near-surface phenomenon and none of the compiled liquefaction case history databases (e.g., Cetin 2000; Moss et al. 2003;Kayen et al. 2013;Boulanger and Idriss 2014) include cases deeper than 20 m. The regression analyses were performed on well over 1,000,000 data points.
Groningen-Zone-Specific r d Relationships For ease of model development, ease of implementation, and consistency with previous work, a relatively simple regression approach was used to develop the Groningen-zone-specific r d relationships. Traditional nonlinear least-squares regression could not be used because of the complicated heteroscedastic standard deviation (i.e., the standard deviation is a function of the predictor variables: depth, z, and moment magnitude, M). Therefore, maximum likelihood estimation was used directly and error estimates for each parameter were obtained from the Hessian matrix of the likelihood function (i.e., negative of the second derivative of the logarithm of the likelihood function).
Several functional forms were considered during the model development, including the functional form adopted in the worldwide model developed by Lasley et al. (2016). The final functional form was based upon a sigmoid shape with the main variable being logarithmic depth. However, the location and scale parameters of the sigmoid are functions of M, V S , and a max (i.e., PGA at the surface of the soil profile). There is an apparent break in scaling for relatively large values of a max ; consequently, this effect was modeled. The dependence on the time-weighted average V S of the upper 12 m of the profile (V S12 ) that was included in the worldwide r d relationship developed by Cetin (2000) was also observed in the Groningen data and was thus included in the relationships developed herein. The same functional form of the regression equation was used for all zones, and the regression coefficients were found to be statistically significant in all cases. The final functional form was selected based on likelihood ratio tests performed on the alternative models considered.
The final functional form for r d-Gron is where z has units of m; β i = regression coefficients (Table 2); and A rd = Eq. (1b) and (1c) and represents the asymptotic level for r d-Gron at depth (that is, r d-Gron ¼ 1-A rd as z → ∞).
The model also has a heteroscedastic standard deviation that is defined by where z is in units of m; and β i = regression coefficients ( Table 2). Note that in Eq. (1d) the depth is defined as the maximum of the actual depth, z, and 5 m to prevent the standard deviation from becoming too small at very shallow depths. The residual plots against the predictor variables used in Eq.
(1) are shown in Fig. 7 for Zone 602, as an example. The very large number of residuals for each analysis case poses challenges because small deviations are statistically significant and can lead to over-fitting. Therefore, in fitting the model, we focused on the main trends that were observed. The binned residuals are shown in Fig. 7 to better illustrate the center and the variance of the data. For example, in Fig. 7(a), there is a larger deviation of the data from the center near logða max Þ ¼ −1 from the residuals, but the binned residuals show that there is not significantly more variance here, but rather, it is just that there is such a large sample that more points are seen in the tails of the distribution. Additionally, it is important to note that the residual plots do not show any clear trends against the variables shown; thus, the model appears to capture the locationspecific scaling of r d well. Fig. 8 shows a comparison of the Groningen-specific r d relationship for Zone 602, as an example, and the worldwide r d  relationships proposed by Idriss (1999) (adopted by both Idriss and Boulanger 2008;Boulanger and Idriss 2014) and Lasley et al. (2016). As may be observed from this figure, for small magnitude events (e.g., M 3.5), the Idriss (1999) and Lasley et al. (2016) relationships significantly overpredict and underpredict, respectively, r d for the Groningen liquefaction study area. However, the lack of accuracy of these relationships for M 3.5 is not altogether surprising because their intended use was for M ∼ 5.0 and greater. For moderate-sized events (e.g., M 5.25), there is still a general trend for the Idriss (1999) and Lasley et al. (2016) relationships to respectively overpredict and underpredict r d , albeit not significantly so, with the Lasley et al. (2016) relationship providing more accurate predictions for the Groningen region than the Idriss (1999) relationship. For large magnitude events (e.g., M 7.0), the three r d relationships predict similar values, with the general trend for the Idriss (1999) and Lasley et al. (2016) relationships to respectively overpredict and underpredict r d still persisting but at a much less significant level.

Groningen-Zone-Specific n eq and MSF Relationships
The functional form for the Groningen-zone-specific n eq relationship is a slight modification of the worldwide model developed by Lasley et al. (2017), with the modifications made to reflect a clear break in scaling at large values of a max and to include the observed dependence on V S12 . The functional form for the Groningen-zonespecific n eq relationship is Because r d-Lea16 and r d-Gron are functions of V S12 , a representative value for the liquefaction study area of 140 m=s was assumed for the plots. Additionally, because r d-Gron is also a function of a max , curves for a max ¼ 0.1, 0.2, and 0.3g are shown in the plots.
where a max is in units of g; V S12 is in units of m=s; α i = regression coefficients (Table 3); and σ lnðneqMÞ for each zone is listed in the last column of Table 3. This model was developed under the tested assumption of log-normality of the equivalent number of cycles. Fig. 9 shows the residuals against the predictor variables used in Eq. (2) for Zone 602, as an example. The same comments made in regards to the residual plots for the r d-Gron model (Fig. 7) also apply to the residual plots for the n eqM-Gron model shown in Fig. 9. In interpreting the residual plots shown in both Figs. 7 and 9, it needs to be realized that any of the small biases that may be observed are smaller than the standard deviations of the binned residuals.
Although the worldwide model for n eq developed by Lasley et al. (2017) included a depth dependence of the standard deviation, there is no compelling evidence of this in the Groningen datasets [ Fig. 9(d)]. Accordingly, a homoscedastic standard deviation (i.e., independent of predictor variables) was adopted for the Groningen-zone-specific model.
Consistent with the revised MSF relationship given in , the Groningen-zone-specific MSF relationship is given by the following expression: where the denominator in Eq. (3a) is computed using Eq. (2). Also, note that the numerator in Eq. (3a) differs from that used by  for the revised worldwide MSF (i.e., 7.25 versus 14). The reason for this is that the revised worldwide MSF relationship was derived using both horizontal components of motions recorded at a site during tectonic events and Eq. (3a) was derived using single components of motions generated using a seismic source model, as discussed. Accordingly, 7.25 represents the reference number of equivalent cycles for one horizontal component of motion for an M 7.5 event in active shallow tectonic seismic zones and is computed using the relationship given in Lasley et al. (2017), Approach 1, WUS: M ¼ 7.5 and a max ¼ 0.35g. The cap of 2.04 on MSF Gron corresponds to a motion composed of one, low amplitude pulse in one of the horizontal components of motion. Fig. 10 shows a comparison of the Groningen-specific MSF relationship for Zone 602, as an example, and the worldwide relationships proposed by Idriss and Boulanger (2008) (MSF IB08 ) and  (MSF WUS ). As shown in Fig. 10, the MSF IB08 relationship predicts significantly larger values for magnitudes less than M ∼ 7.0, with the overprediction being more significant as magnitude decreases. However, the MSF WUS and MSF Gron do not differ significantly across magnitude, with the MSF Gron showing a slightly less dependence on a max .
Correlations between Groningen-Zone-Specific r d and lnn eq Relationships The results show that r d-Gron is correlated across depths and that lnðn eqM-Gron Þ and r d-Gron are negatively correlated at a given depth. The correlation coefficient of r d-Gron at depths z i and z j is given by the following expression where ε rd-Gron = number of standard deviations from the median r d-Gron value; z i;j are in m; and α rd-Gron is dimensionless (Table 4). Also, the correlation coefficients for lnðn eqM-Gron Þ and r d-Gron (i.e., ρ lnðneqM-GronÞ;rd-Gron ) are listed in Table 4. In the next section, the Groningen-specific relationships, r d-Gron and MSF Gron , in conjunction with the revised CRR curve , are used to assess the liquefaction hazard across the Groningen liquefaction study area shown in Fig. 6.

Groningen Liquefaction Hazard
Liquefaction hazard curves were calculated for 95 sites across the study area using a Monte Carlo approach. The probability distributions for seismic activity rates , 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 . 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 Stafford et al. 2017). For each event scenario, the Groningen-specific relationships, r d-Gron and MSF Gron , in conjunction with the revised CRR curve , were used to compute the FS liq as a function of depth. All other relationships required to compute FS liq were adopted from Boulanger and Idriss (2014), consistent with the development of the revised CRR curve. The LPI ish framework (Maurer et al. 2015b) was used to relate the computed FS 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 LPI ish framework is a conceptual and mathematical merger of the Ishihara (1985) H 1 -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 H 1 -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), LPI 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 FS liq in each layer is less than 1.0.
The optimal LPI ish thresholds corresponding to different severities of surficial liquefaction manifestations are dependent on the liquefaction triggering procedure used to compute FS 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 LPI ish framework i.e.: LPI ish < 5: no to minor surficial liquefaction manifestations; 5 ≤ LPI ish ≤ 15: moderate surficial liquefaction manifestations; LPI 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 LPI 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), LPI ish values corresponding to an AFE of ∼4 × 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 LPI 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, LPI 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, LPI 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, LPI ish ∼ 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 LPI ish < 5.

Conclusions
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 r 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 Groningenspecific and the NPR 9998-2017 specified IB08 liquefaction triggering models manifest in the computed liquefaction hazard curves (i.e., AFE versus LPI 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, LPI ish values corresponding to an AFE of ∼4 × 10 −4 (or a 2475-year return period) are of particular interest. However, assessing the liquefaction hazard by determining the LPI 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 LPI 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 LPI 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.