Ground-motion prediction models for induced earthquakes in the Groningen gas field, the Netherlands

Small-magnitude earthquakes induced by gas production in the Groningen field in the Netherlands have prompted the development of seismic risk models that serve both to estimate the impact of these events and to explore the efficacy of different risk mitigation strategies. A core element of the risk modelling is ground-motion prediction models (GMPM) derived from an extensive database of recordings obtained from a dense network of accelerographs installed in the field. For the verification of damage claims, an empirical GMPM for peak ground velocity (PGV) has been developed, which predicts horizontal PGV as a function of local magnitude, ML; hypocentral distance, Rhyp; and the time-averaged shear-wave velocity over the upper 30 m, VS30. For modelling the risk due to potential induced and triggered earthquakes of larger magnitude, a GMPM for response spectral accelerations has been developed from regressions on the outputs from finite-rupture simulations of motions at a deeply buried rock horizon. The GMPM for rock motions is coupled with a zonation map defining frequency-dependent non-linear amplification factors to obtain estimates of surface motions in the region of thick deposits of soft soils. The GMPM for spectral accelerations is formulated within a logic-tree framework to capture the epistemic uncertainty associated with extrapolation from recordings of events of ML ≤ 3.6 to much larger magnitudes.


Introduction
Conventional gas production in several fields in the north and northeast of the Netherlands has resulted in induced earthquakes (van Eijs et al. 2006;van Eck et al. 2006). Natural seismicity in the Netherlands is generally confined to the south of the country (Fig. 1), a notable example being the 1992 Roermond earthquake of local magnitude M L 5.8 and moment magnitude M 5.3 (Camelbeeck et al. 1994;Haak et al. 1994;Dufumier et al. 1997). The Groningen gas field, which is by far the largest in the Netherlands, has produced the greatest number of induced earthquakes and the largest induced event in the Netherlands, namely the Huizinge earthquake of August 2012, assigned magnitudes of M L 3.6 and M 3.5.
The Groningen gas reservoir is situated in the Rotliegend-Slochteren sandstone, a formation of up to ~ 300 m in thickness, located at a depth of about 3 km. The reservoir is offset in numerous locations by faults that have not been active in geologically recent times. Pressure depletion as a result of gas production has led to compaction of the reservoir that leads to stress increases across the offsets, which has caused the reactivation of many of the faults (e.g. Buijze et al. 2017;Bourne et al. 2018). The first induced earthquake to be detected (M L 2.4) occurred in 1991, about 28 years after production began, when the maximum reservoir compaction had reached about 20 cm (Bourne et al. 2014).
In response to the induced earthquakes, the field operator NAM (Nederlandse Aardolie Maatschappij) initiated the development of risk analysis model to estimate the impact of potential larger induced earthquakes and even larger triggered events (van . The express purpose of the a risk model was to inform decision-making regarding risk mitigation options, through the capability to model the consequences of reduced production levels (and hence lower additional reservoir compaction) on the seismicity rates and the effects of building strengthening on reducing fragility (van . The vital link between these two elements-the seismicity model and the exposure and fragility models-is a model for the prediction of ground-motion levels due to each earthquake scenario. An obvious candidate for the ground-motion prediction, at least for peak ground acceleration (PGA) and velocity (PGV), were the equations derived by Dost et al. (2004) using recordings from induced earthquakes in the Roswinkel gas field supplemented by a few recordings from a swarm of tectonic earthquakes in the most southern region of the Netherlands near the city of Voerendaal (there were no ground-motion recordings in the Roermond earthquake within 50 km of the epicentre). However, the recorded motions from the Groningen field indicate systematic differences in the amplitudes of ground shaking compared to those from Roswinkel, with the Groningen motions being consistently lower, as shown in Fig. 2. The earthquakes in Roswinkel occurred at depths of about 2 km, whereas the Groningen events occur at the reservoir depth of 3 km. However, this difference in the ground motions observed from induced earthquakes in the two fields is believed to be mainly due to differences in the geology overlying the gas reservoirs. In the Groningen field, the gas-bearing sandstone, where the earthquakes originate, is situated below a thick layer of Zechstein salt (including thin layers of very highvelocity anhydrites), which has been shown to exert a pronounced influence on the propagation of seismic waves (Kraaijpoel and Dost 2013). In the Roswinkel the empirical model additionally account for this factor. • For risk analyses in terms of structural damage and injury, ground-motion prediction models for response spectral accelerations are required for magnitude 2.5 to above 7, to account for triggered seismicity. Finite fault simulations, calibrated to the local recordings, are used to predict motions at a buried rock horizon, and these are combined with a field zonation for frequency-dependent, non-linear amplification factors to obtain surface motions. The model is developed in a logic-tree framework to capture the uncertainty associated with the large extrapolation in magnitude. gas field, however, the reservoir is situated above the high-velocity Zechstein formation. The surface soils in the Roswinkel field are generally stiffer than in Groningen, with Holocene soils (clays and peats) being largely absent; hence, greater site amplification would be expected in Groningen than in Roswinkel. These observations reinforced the decision that it was necessary to develop a Groningen-specific ground-motion prediction model (GMPM). In addition to requiring a model that captures the influence of the upper crustal structure in the Groningen field, it was also considered highly desirable to develop a GMPM to reflect the specific site response characteristics of the thick soil layers, including soft clay and peat deposits, that overlie most of the field. Consequently, ground-motion models have been constructed that are exclusively applicable to the Groningen gas field and surrounding areas and cannot be transported to other regions. However, the framework in which the model has been developed and some of the insights obtained through this iterative process of model development may be of value to others faced with predicting ground motions from induced and triggered earthquakes.

An evolutionary development of ground-motion models
The GMPMs presented in this paper are the final products of an evolutionary process that commenced in 2013 and has passed through several iterations. The multiple stages of development were driven by three parallel processes, the first being the regulatory requirement for regular updates of the hazard and risk estimates for the Groningen field. This resulted in fully documented versions of the GMPM being generated at relatively short intervals, which brought both advantages and disadvantages. On the positive side, the incremental development provided opportunities for a rigorous review of each version of the GMPM and testing of models through the implementation in the hazard and risk calculations. On the negative side, the sprints to produce each version of the model, combined with the computational burden of the ground-motion simulations and site response analyses (see Section 5) meant that there was little opportunity to iterate on individual elements of the model as each complete version of the GMPM was generated. Moreover, two versions of the model were affected by errors that arose due to the tight schedule: one version was influenced by a preliminary relationship between local and moment magnitude, which was subsequently revised, and another by errors in processing the large number of recordings from an earthquake that occurred just as the model-building process was starting. The second driver for the generation of multiple versions of the GMPM was the growth of the database, both in terms of the number of ground-motion recordings available and the geotechnical characterisation of the recording stations. The groundmotion database expanded very considerably during the course of the model development, due to ongoing seismicity and the installation of new recording networks (see Section 3). Campaigns of in situ field measurements also led to most of the recording stations eventually having measured near-surface shearwave velocity, V S , profiles.  PGA (upper) and PGV (lower) values recorded at epicentral distances of less than 2 km from induced earthquakes in the Roswinkel and Groningen fields, plotted against magnitude. The points are plotted at values fractionally above (triangles) and below (circles) the exact magnitude in order to separate the symbols for clarity; the data shown are limited to the very short distances at which ground motions were recorded in Roswinkel The third driver of the model development was related to a progressive improvement of the modelling process to better capture the physical characteristics of the earthquakes, the wave propagation, and the dynamic response of the soft soil deposits. These refinements were prompted both by the insights of the authors and the suggestions and challenges of the review panel established to oversee the work (see Acknowledgements). For obvious reasons, the opportunities to refine the model were often closely linked to the availability of data.
Another factor in the development of the GMPMs, affecting the left-hand side of the equations rather than the functional form and parameters on the right, has been adapting to the evolution of the fragility functions defined for the building typologies encountered in the Groningen region (Crowley et al. 2017a, b;Crowley et al. 2019) and various potential applications of the model. At one time, this meant that the model included prediction of the duration of motion, but this is no longer required.
The GMPM development has also had to be compatible with the seismicity model that it links to the exposure and fragility models for the calculation of risk. The seismological model is based on local magnitude, M L , as determined by the Dutch seismological service, KNMI. However, it has been demonstrated that for magnitudes larger than 2.5 (the lower limit for which the GMPM is applicable), M L in Groningen can be considered equivalent, on average, to moment magnitude M (Dost et al. 2018. The very first GMPM (version 0 or V0) was developed very rapidly and only predicted PGA and PGV. The model was essentially based on the equations of Akkar et al. (2014) developed for tectonic seismicity in the Mediterranean and Middle East regions, adjusted at smaller magnitudes (below 3.8 for PGV and below 4.2 for PGA) to match the small database of Groningen ground motions available when this work commenced (Bourne et al. 2015). At that time, estimates of V S30 (the time-average V S over the uppermost 30 m) were not available for any of the recording stations; the model was applied assuming a fieldaverage V S30 of 200 m/s, which has subsequently been confirmed as an accurate estimate.
The first complete GMPM (V1) was developed by inverting Fourier amplitude spectra of recorded surface motions to estimate the source ( Δ ), path (Q), and average site parameters (amplification and site damping, 0 ). These parameters were then used in point-source stochastic simulations using SMSIM (Boore 2005) with different values of the stress parameter Δ in order to capture the epistemic uncertainty in extrapolating to larger magnitudes. While this model did reflect Groningen conditions, it had two important weaknesses in terms of the site response: the amplification corresponded to a network average (the network being located in the northern part of the field, where V S30 values are generally lower), and only linear site response was modelled, even up to large magnitudes. These shortcomings were addressed in the V2 GMPM, which established the basic framework used for all subsequent models. Fundamental to being able to implement this framework were the development of a field-wide V S model (Kruiver et al. 2017a, b) and in situ measurements of V S profiles at the recording stations (Noorlandt et al. 2018). In this GMPM, the surface motions were deconvolved to a buried rock horizon (at ~ 350 m depth), assuming linear site response because of the relatively low levels of motion. The same procedure as used for the V1 model was then applied to these rock motions to provide the inputs to stochastic simulations of motions at the rock horizon for a wide range of magnitudes. Frequency-dependent, non-linear amplification factors were then determined and assigned to zones, originally defined on the basis of near-surface geology, to transform the rock motions to the surface.
All subsequent versions of the model represented modifications to the V2 framework. In the V3 model, the reference rock horizon was moved to the base of the North Sea Supergroup at ~ 800 m depth, which is a more persistent and more pronounced impedance contrast . At this stage, the models were based on point-source simulations, despite the fact that the hazard and risk calculations considered a maximum magnitude, Mmax, of 6.5-determined on the basis of all the strain associated with full production of the gas reservoir being released seismically in a single event (Bourne et al. 2014). Epicentral distance, R epi , was used as the distance metric because KNMI assigned a common focal depth of 3 km to all Groningen earthquakes. To compensate for the point-source approximation, a magnitude-and distance-dependent increase was imposed on the sigma value representing the aleatory variability associated with the predicted motions (Bommer et al. 2016). A workshop dedicated to the estimation of Mmax for Groningen held in March 2016 resulted in a distribution of values peaking at 4.5, but with low-weighted upper tail exceeding magnitude 7 . This prompted a move from V4 onwards to finite rupture simulations, using the software EXSIM (Motazedian and Atkinson 2005) as modified by Boore (2009). In conjunction with this change, the distance metric was changed from R epi to rupture distance, R rup , although for the small-magnitude events in the database, R Rup was considered equivalent to hypocentral distance, R hyp . At this point, the previously mentioned adjustment to the sigma values was dropped. Another change that was implemented in the V4 model was the inclusion of magnitude-and distance-dependence in the short-period linear site amplification factors, AFs (Stafford et al. 2017).
The framework of the V4 model has remained unchanged, with each subsequent version representing a refinement based on improved modelling, for example of the soil damping model , and the continued growth of the groundmotion database. Up to the V5 model, the derivation used surface recordings from the B-network and geophone recordings at 200 m depth from the G-network (see Section 3), because measured V S profiles were available at the former but not the latter, and the intention was to avoid the additional uncertainty in the transfer functions and AFs used to deconvolve the G-station recordings. However, doubts had arisen regarding the consistency of the surface and borehole recordings, possibly related to the operational issues with the borehole geophones, prompting the decision to use only surface accelerograms in the V6 model; this also meant that the deconvolution and forward modelling were along common paths. For the V7 model, measured V S profiles had become available for nearly all the surface accelerographs in the G-network, removing entirely the original motivation for using borehole records. Additionally, the V7 model made use of more accurate depth determinations by Spetzler and Dost (2017), which reflected the convex shape of the reservoir over the field, in the calculation of R hyp .
The V7 model presented herein represents the culmination of this evolutionary process. The V7 development benefited from an intervention by the international expert review panel, who counselled against running through yet another sprint to produce another incremental modification of the V6 model. This resulted in the V6 model being used for another round of hazard and risk calculations, but this created the space for individual elements of the V7 model to be developed, reviewed, and refined prior to initiating the full model production process. Consequently, this is a model whose development has benefited from multiple iterations and enhancements and likely represents an asymptote in terms of potential improvements although we highlight some issues that would be worthy of additional investigation.
In this paper, we present the current groundmotion models for Groningen but do not make comparisons with the predictions from earlier versions of these models. Such comparisons are not straightforward in some cases because of changes of independent variables such as the distance metrics, but we believe that they are also of little value since each new model supersedes the previous development stages. Moreover, each model is judged on its own merit with respect to the data and the modelling decisions, not on whether it predicts higher or lower amplitudes of motion than the previous version. We begin by briefly summarising the ground-motion database and then we present an empirical model for PGV, which provides some useful insights into the nature of the recorded ground motions in the field, before summarising the current GMPM for seismic risk analysis in Groningen.

Groningen ground-motion database
The characteristics of the various seismic recording networks in and around the Groningen field are described by Dost et al. (2017) and Ntinalexis et al. (2019). The networks that provide the data used in the derivation of the GMPMs include the B-network, which comprises accelerographs installed by KNMI in houses and farm buildings mainly in the northern part of the field. Following the 2012 Huizinge earthquake, this network was expanded from 10 to 18 stations, and the instruments were upgraded to Kinemetrics Episensors; for reasons that will be made apparent in Section 4, the networks prior to and after the upgrade are referred to as B_old and B_new, respectively. The G-network, which is also operated by KNMI, was installed by NAM between 2014 and 2017, comprising 79 surface accelerographs (also Episensors), 69 of which are co-located with a borehole instrumented with geophones at 50, 100, 150, and 200 m depths. The magnitude-distance distribution of the recordings, classified by recording network and by V S30 , is shown in Fig. 3. As can be appreciated, the stations are generally rather soft, with most records corresponding to V S30 values in the range from 160 to 260 m/s.
The records have all been processed uniformly following the procedure described by Edwards and Ntinalexis (2021). The database is described in greater detail by Ntinalexis et al. (2022), which also provides a link to a repository from which the entire collection of raw and processed accelerograms can be freely downloaded.

Empirical model for peak ground velocity
Several early versions of the GMPM included predictions for PGV, as well as response spectra accelerations and durations, for potential applications to risk analysis of some infrastructure lifelines, but this was discontinued since we were not aware of any such applications in practice. However, the need arose for a PGV model to be used in the evaluation of damage claims, the cumulative value of which has reached unprecedented levels Bommer 2022). Since this model was to be applied only within the magnitude range of the data, an empirical groundmotion prediction equation (GMPE) was considered sufficient. For the derivation of this GMPE, a lower magnitude of M L 1.8 was selected to cover the full range of felt events in the field. The recorded PGV values, grouped by magnitude ranges and plotted against hypocentral distance, are shown in Fig. 4. The largest PGV value in the database is 3.5 cm/s, which was recorded during the 2012 Huizinge earthquake.
The data in Fig. 4 display a distinct feature of the Groningen ground motions in terms of geometric spreading: rapid decay over short distances, a flattening between about 7 and 12 km, and more gradual decay at greater distances. This is consistent with the results of finite difference waveform modelling using a detailed 3D velocity model of the field , which has been used to constrain the distances at which these transitions occur.
Modern GMPEs usually predict the geometric mean value of the two horizontal components of each accelerogram or other component definitions that result in approximately the same values at most periods (e.g. Boore 2010). For the PGV GMPE for Groningen, however, coefficients were also derived for the larger horizontal components and maximum rotated. This component definition was chosen because it is consistent with the definition of the V top parameter defined in Dutch guidelines for tolerable vibration thresholds; the model is presented only for the larger component.
The GMPE for the larger component obtained from regression on the empirical database is as follows: with PGV in cm/s, M being local magnitude M L , and V S30 is the harmonic average shear-wave velocity over the top 30 m (in m/s). The normalising V S30 value of 200 m/s was chosen on the basis of the approximate field average, but using another value would exert no influence on the predictions (it only changes the constant term). The distance term, R (in km), is an effective distance that accounts for magnitude-dependent near-source saturation effects, is a function of R hyp , and is defined as follows: with the magnitude-dependent saturation term set equal to as follows: (1) In Eq. (3), the multiplier on magnitude is fixed to the constant value of 1.1513 so that the saturation distance scales in proportion to the expected rupture area, whereas the first coefficient is solved for in the regression. The distance dependence is defined by the following: The rapid decay of amplitude at short distances is interpreted to be the result of the refraction of the seismic waves by the high-velocity Zechstein salt layer (Kraaijpoel and Dost 2013).
The total variability in the predictions, σ, is decomposed into between-event (τ), site-to-site ( S2S ), and within-event ( ss ) standard deviations, all of which are related by the following expression (e.g. Al Atik et al. 2010): Fig. 4 As-recorded horizontal PGV values from the Groningen ground-motion database, grouped by magnitude, plotted against hypocentral distance The estimates of these variance components as well as the coefficients of the GMPEs are also influenced by considering the uncertainty in magnitude determinations for each earthquake. The resulting values are summarised in Table 1.
An interesting observation that can be made is that the influence of V S30 is rather modest: decreasing the V S30 value from 260 m/s to 160 m/s results in only a 17% increase in the predicted median value. This is consistent with findings in other projects-including NGA-Sub (Parker and Stewart 2021) and NGA-East (Parker et al. 2019)-which find a flattening of the relative amplification as V S30 values become very low (< ~ 300 m/s). As a result of damage claims being submitted from a much larger area than that defined by the boundaries of the gas field, to facilitate the implementation of the PGV GMPE, a map of V S30 , defined at the level of 4-digit postcodes, has been generated ( Fig. 5) following essentially the same procedures applied by Kruiver et al. (2017b) to construct the uppermost part of the field-wide V S profiles.
Analysis of the residuals revealed that the models provide a good fit to the data, with no discernible trends with respect to magnitude, distance, or V S30 (Fig. 6). However, the residuals do show that the motions recorded at the B_new stations are, on average, lower than those from the G-network and the B_old stations. An apparent tendency for shortperiod motions-PGV in Groningen has been found to correlate well with the spectral acceleration at 0.2-0.3 s (Bommer et al. 2017b)-to be lower at the B-stations than the G-stations has been attributed by some researchers to soil-structure interaction (SSI) effects, which is an obvious candidate given that the B-network stations are installed in buildings and the G-network stations are free-field accelerographs housed in small instrument shelters. The same researchers recommended that all the B-station recordings should therefore be removed from the database, which would mean the loss of many of the strongest motions and the majority of the near-source recordings of the larger events in the database (Fig. 3). Cavalieri et al. (2021), however, using advanced dynamic SSI analyses, demonstrated robustly that SSI effects do not explain any differences between B-network (not including those with basements) and G-network recordings in the short-period range. The same study showed that the most likely explanation is the presence of a layer of improved soil below the buildings housing the stations of the B-network, which are commonly created by pre-loading the soil with a layer of compacted sand prior to construction. No systematic correction for this effect can be made, however, since no information is available regarding the presence or absence of such soil layers at all the B-network stations. At the same time, it is noted that there are clear SSI effects at three of the B-network stations where the buildings include deep basements (BUHZ, BWIN, and BZN1), but corrections for the embedment effect at these stations are applied to all the recordings following the procedures recommended in NIST (2012). Ongoing investigations are underway in an attempt to identify and quantify the exact cause of the apparent differences between the B-station and G-station recordings, taking account of recent insights into the influence of instrument installation characteristics on the nature of recorded high-frequency motions (Hollender et al. 2020).
The interesting question that then arises is why this effect should manifest only in the upgraded B-stations (B_new) and not in B_old recordings. The recordings from B_old instruments, which operated on a triggering basis, were generally selected by KNMI for processing and distribution, and this manual process would naturally have favoured the stronger recordings, and lower amplitude motions will have been excluded. Analyses have indicated that the attenuating effect at the B-network stations being a function of amplitude and manifesting most consistently for lower levels of motion, although it should be clearly noted that no physical explanation for this observation has been put forward. The GeoSig accelerographs of the B_old stations displayed much higher noise levels than the new Episensor accelerographs ).  Consequently, an average effect of high-frequency suppression at the B-network stations may be genuine, but it could have been largely removed from the B_old recordings as a result of data selection procedures. In view of the apparent influence of the recording network, an alternative form of Eq. (1) was also derived, including an additional term with a dummy variable F NB that takes a value of 0 for B_new recordings and 1 otherwise: The near-source saturation term of Eq. (3) and the distance dependence defined in Eq. (4) become as follows: The inclusion of this term to distinguish the contributing accelerograph network leads to a small reduction in the total variability, as shown in Table 1. Figure 7 compares the median prediction from Eqs.
(1) and (6), which shows the rather pronounced reduction in amplitudes at the B_new stations.
This raises the question of which equation should be used in applications. Since the objective is to estimate the level of shaking experienced by typical buildings in the region, the most appropriate model would be Eq. (6) with F NB set to 0. However, there are two considerations that deter us from this recommendation, the first being that it is not known at what proportion of buildings in the Groningen region do the effects of attenuated high-frequency motions manifest. The current hypothesis is that the effect is due to improved soil layers below the foundations, but there are no data available to determine how pervasive this feature is among the entire building stock. The second consideration, which is not scientific but rather related to optics (and taking account of the regulatory environment and the nature of the scrutiny to which the hazard and risk modelling is subjected), is that this would mean adopting the model yielding the lowest predictions, which would likely be treated with suspicion. Our recommendation, therefore, is to use the model derived from all the available recordings, combining the free-field G-stations and the B-stations within buildings, namely Eq. (1). The implicit assumption in this choice is that collectively, the recordings obtained from the stations of both networks approximate the distribution of buildings with and without the high-frequency attenuation effect.
Additional details on the derivation of the PGV models are presented in Bommer et al. (2021), where the coefficients for the prediction of other horizontal component definitions can also be found.
A final point that can be made here is that the variability components presented in Table 1 allow the event term to be calculated for any earthquake using the recordings that become available on the KNMI portal soon after each earthquake. The event term can then be subtracted from the constant term in the equation, and predicted motions can be estimated sampling only from the within-event variability.

Models for the prediction of response spectral accelerations
The GMPM developed for seismic risk analysis predicts 5%-damped ordinates of pseudo-spectral acceleration, Sa(T), at 10 oscillator periods, T (0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.85, and 1.0 s); the fragility functions for buildings in the Groningen exposure model are defined in terms of the geometric mean value of Sa over these periods. The model is developed to predict the geometric mean horizontal component of motion, for application to hazard assessments, but the appropriate definition for the risk analyses is the arbitrary horizontal component of motion, which requires an adjustment to the standard deviation to account for the component-to-component variability (Baker and Cornell 2006). The GMPM needs to be applicable over the full range of magnitudes considered in the risk calculations, with a lower limit set at 2.5 and upper limit at 7.25, also accounting for the uncertainty associated with this large extrapolation beyond the upper limit of the data. As a result of considering much stronger earthquake scenarios than those that have occurred in the field, the GMPM was also required to consider the non-linear site response of the soft soil deposits that cover most of the region of interest. As with the PGV predictions, our choice was to use all the available ground-motion data, regardless of the network from which it was derived, provided there was a measured V S profile for the recording station. This is a less critical choice than in the derivation of the empirical PGV model since the model is not obtained by direct regression on the recordings, but the data are used in the estimation of the variance components.

Model overview
The overall scheme for the derivation of the V7 GMM is illustrated schematically in Fig. 8. The essential elements of the GMM are equations to predict Sa at the NS_B reference rock horizon and non-linear amplification factors (AF) that transform the NS_B motions to surface accelerations. The derivation of the model for the prediction of Sa at the NS_B horizon and of the AFs is explained in Sections 5.2 and 5.3, respectively, including the quantification of both aleatory variability and epistemic uncertainty in both parts of the model. The basic framework is presented below, but the reader interested in greater detail regarding the model derivation can refer to Bommer et al. (2022). Fig. 7 Predicted median PGV values against R hyp for three magnitudes from the models with and without the term to distinguish the contributing accelerograph network Using the V S profiles down to the NS_B horizon at each recording station, together with the damping profiles , transfer functions (TFs) are developed to deconvolve the Fourier amplitude spectra (FAS) to the reference rock horizon; here, as in all the site response modelling, we make the standard assumption of vertically propagating waves. The high-frequency damping parameter is estimated using the method of Anderson and Hough (1984) and then the FAS are inverted to estimate the geometric spreading coefficients (using the hinge distances constrained by finite difference simulations of the wave propagation), the average quality factor Q (using a layered Q model), and the stress parameter Δσ. Due to the inevitable trade-offs among these parameters, there are many combinations that are found to be consistent with the FAS. The next step is to calculate scenario-dependent AFs using the site profiles and the available recordings from each station. These AFs are then used to deconvolve the Sa to the NS_B horizon. With the attenuation functions fixed (R −x and Q), point-source simulations (using SMSIM) are used to generate response spectra at the NS_B horizon for the M-R pairs of the recordings. The average bias between the predictions and the deconvolved response spectra, averaged over the 10 target oscillator periods, is calculated to determine the optimal combination of Δσ and to match the recorded motions (black cross in left-hand panel of Fig. 9). To capture the epistemic uncertainty, two alternative Δσ-combinations that also lead to low bias (blue crosses in left-hand panel of Fig. 9) are also selected. The central model is defined by Δσ = 22 bars and equal to 0.002 s; the upper and lower values of the stress parameter are 33 and 15 bars.
On the basis of conclusions from the expert panel appointed to develop the Mmax distribution for Groningen (Bommer and van Elk 2017), any earthquakes of magnitude greater than 5 would be assumed to be Fig. 8 Schematic illustration of the process to derive the V7 GMM. TF is transfer function, AF is amplification factor, and R −x represents geometric spreading with R being distance. AH84 Anderson and Hough (1984) to estimate . The triangles represent surface accelerograph stations of both the B-and G-networks. R EFF is a magnitude-dependent effective distance (Boore 2009) triggered events that rupture largely outside of the reservoir (which would imply rupture propagation downwards into the Carboniferous since faults would not rupture the overlying salt deposits). In order to create a branch of the logic-tree to accommodate such tectonic earthquakes (even though large ruptures initiating at 3 km depth and propagating downwards would be unusual), values of the stress parameter were determined to match, over the 10 target periods, Sa predicted from tectonic GMPMs. Using the models of Abrahamson et al. (2014), Akkar et al. (2014), Bindi et al. (2014), Boore et al. (2014), Campbell and Bozorgnia (2014), and Chiou and Youngs (2014), median Sa values were predicted at the NS_B using V S30 = 1400 m/s for magnitudes 5, 5.5 and, 6 for logarithmically spaced distances up to 20 km. Simulations were then performed for the same M-R combinations using the central value of to explore the Δσ value that best matched the spectra from the tectonic earthquakes. For efficiency, these simulations were performed with SMSIM using the magnitude-dependent effective distance to account for finite-source effects, which has been shown to be reasonable approximation (Boore 2009). The optimal value of Δσ, giving the smallest bias in the predicted spectra at NS_B, was found to be 50 bars (see right-hand panel of Fig. 9).
To model the epistemic uncertainty in the predictions, a logic tree is established with three branches for the value of the stress parameter in the magnitude range of the data (i.e. 2.5 to 3.6), with the dominant weight (0.6) on the central branch and weights of 0.2 on the upper and lower alternative branches. The range of epistemic uncertainty for the predictions at magnitudes of 5 and greater must be large by virtue of the extrapolation beyond the data and the fact that nothing is known about the fault ruptures that would be associated with such events. The end members were therefore defined by maintaining the lowest stress parameter branch constant across the magnitude range, at the lower end, and converging to the branch calibrated to tectonic GMPMs at the upper end. To adequately cover the space between these two limits, the central branch bifurcates for magnitudes above 3.6 and converges to branches that are equally spaced in terms of log(Δσ) for magnitudes of 5 and greater (Fig. 10, upper). In the large magnitude range, the lowest branch is considered feasible but somewhat unlikely, so its weight is reduced to 0.1, and the distribution of weights becomes skewed towards the higher branches; consequently, the branch weights are not constant but vary with magnitude (Fig. 10, lower).
The four values of Δσ, each coupled with a value of 0 at the NS_B, are used to model the sub-events in finite fault simulations using EXSIM, which are performed for events covering the full magnitude range in increments of 0.25 units. For each earthquake, a fault rupture is defined with a dip consistent with the observed geometry of faults in the field and dimensions calculated from Wells and Coppersmith (1994). The hypocentre is located randomly along the strike Fig. 9 Left: bias of predicted Sa at the NS_B with respect to the recorded Sa deconvolved to this horizon, as a function of kappa and the stress parameter, the black cross indicating the best-fit solution and the blue crosses alternative values for the upper and lower branches. Right: similar plot for the bias with respect to predicted Sa from several tectonic GMPMs as a function of the stress parameter with kappa fixed to the central value of the fault along the top of the rupture (i.e. within the reservoir) and the fault propagates downwards (Fig. 11). Waveforms are generated at receiver locations directly above the top of the rupture and at logarithmically spaced distances out to 60 km at various azimuths, with additional recording stations situated at the hinge distances that define the different segments of geometric spreading. Response spectral accelerations are calculated from the waveforms to create a suite of simulated ground-motions for all magnitude and distance combinations.
One important feature of the simulations to be noted is that although the target is the NS_B at 800 m depth, the path is extended to the surface with infinite Q (i.e. no damping) in the layers above the reference rock horizon. The motivation for this choice is the rather unusual situation where a significant proportion of the travel path is modelled through site response analyses, which assume 1D vertical wave propagation and therefore do not account for geometrical spreading in the uppermost portion of the travel path. This was considered to be preferable to the alternative of only propagating the motions to the NS_B horizon, which would ignore the geometric spreading over about one-quarter of the total travel path to the surface (for epicentral motions) and also result in higher motions at the NS_B, leading to more nonlinear response in the overlying layers.
To facilitate the implementation of the groundmotion predictions, parametric equations are derived through regression on the spectral accelerations generated from the simulations. The equations are a function only of magnitude and distance, but the functional form is somewhat showing the values of stress parameter (upper) and the associated branch weights (lower) as these vary with magnitude complex in order to capture the non-linear magnitude scaling over the range of magnitudes for which the model is applicable and the segmented geometric spreading. In contrast to the PGV models, the geometrical spreading is also magnitude dependent. The overall form of the equation for the prediction of Sa at the NS_B horizon is as follows: in which the source-related function is as follows: The hinge magnitude separating the two quadratic magnitude scaling functions, M m , is 4.75. The pathrelated function is as follows: Consistent with Eq. (8), the hinge distances separating different rates of spreading are defined as R h,0 = 3 km, R h,1 = 7 km, R h,2 = 12 km, and R h,3 = 25 km, while the rates r 0 , r 1 , r 2 , r 3 are defined by the expressions in Eqs. (12) to (15): These equations require the determination of up to 21 coefficients at the shortest periods, with slightly fewer being required at longer periods. However, we also have a very large dataset of EXSIM simulations (12) Fig. 11 Schematic illustration of the forward simulations using EXSIM. M 0 is the seismic moment and T S the rise time, which is taken as the reciprocal of the corner frequency; β is the shear-wave velocity and the rupture propagates at 0.8β; the right-hand shows the V S and Q profiles down to 4 km, including the top of the Carboniferous; and W and L are the width and length, respectively, of the fault rupture and empirical observations that provide an excellent constraint on the initial estimates of these coefficients. Examples of the predicted median Sa values from the four models that occupy the logic-tree branches are plotted in Fig. 12. The plots clearly show the strong magnitude dependence of the geometric spreading and the increase of epistemic uncertainty with increasing magnitude.
The aleatory variability is defined in the same way as for the PGV model (Eq. (5)), but in this case, the site-to-site variability is accounted for in the AFs determined from site response analyses for the overlying layers (see Section 5.3); hence, the variability at the NS_B horizon is given by the single-station sigma (Atkinson 2006) defined by the following: The model for 2 SS is based upon work originally conducted by Rodriguez-Marek et al. (2013) that was subsequently updated by Al Atik (2015). The approach assumes that the variance 2 SS is χ 2 -distributed and to select quantiles of this χ 2distribution as discrete logic-tree nodes, two equally weighted branches are adopted corresponding to the 16th and 84th percentiles of the distribution. The within-event variability is independent of magnitude and distance but varies with the oscillator period.
The between-event variability is determined from the residuals of the Sa of the ground-motion records deconvolved to the NS_B with respect to the median models, accounting for the uncertainty in the magnitude estimate for each earthquake. The results are found to display a small dependence on period, but this is sufficiently small to be ignored; hence, a constant value was adopted for the entire period range (Fig. 13). The final model is based on three branches that correspond to the 5th, 50th, and 95th percentiles, with weights assigned based on standard discrete approximations. The logic tree for the variance components therefore has two nodes with three (τ) and two branches ( SS ), respectively; values of single-station intra-event have been found to be relatively stable globally, whereas there is little or constraint-and hence greater uncertainty-on the inter-event variability in the model at larger magnitudes. These six branch combinations are each combined with each of the four branches for the median predictions, leading to a total of 24 (16) ss = √ 2 + 2 ss separate predictions for each realisation of the NS_B motions.
For risk calculations, as noted previously, it is necessary to transform the predicted spectral accelerations from the geometric mean component to the arbitrary component, for which the component-tocomponent variability σ c2c must be added to the variance. A magnitude-and distance-dependent model for σ c2c has been derived from the highly polarised Groningen data ; the model is calibrated to converge to the values found by Campbell and Bozorgnia (2007) from tectonic earthquake recordings for distances of 25 km and greater and for magnitude of 5.6 and larger. To correctly calculate the average spectral acceleration required as input to the fragility functions, the period-to-period correlation structure of Baker and Jayaram (2008), which has been found to be consistent with the Groningen data (computing inter-period correlations for residuals at the NS_B horizon, which is equivalent to looking at the correlation structure of the surface motions since deterministic site response effects are used to map from the surface to NS_B), is imposed.

Site response model
The V S and lithological model for the entire study region-which is defined by the boundary of the gas field plus a 5 km buffer onshore-is developed in the form of voxel stacks of area 100 × 100 m that extend from the ground surface to the NS_B horizon. The construction of the profiles is described in Kruiver et al. (2017a, b), with some minor modifications and refinements that have been made subsequently. The modulus damping and reduction (MRD) curves of Darendeli (2001) and Menq (2003) are selected for clays and sands, respectively, whereas the MRD curves derived by Zwanenburg et al. (2020) from tests of local soil samples are applied to the Holland peat; for older and more consolidated deposits at depths greater than ~ 350 m (the Lower North Sea formation), a linear response is assumed. While this is perhaps an unusually large depth over which to assume a nonlinear response, the soft nature of the Groningen soils means that large strains are possible for some cases, even at considerable depths. The site response analyses are performed using the STRATA program by Kottke and Rathje (2008) using the random vibration theory (RVT; Rathje and Ozbey 2006) option, which requires input in the form of FAS at the NS_B, which are obtained from the simulations. Although refinements have been made to the model and to the procedures, the process for obtaining the AFs is essentially as described by Rodriguez-Marek et al. (2017).
The amplification factors are grouped into 160 zones over the field (Fig. 14). The zonation was originally defined on the basis of near-surface geology ) and subsequently modified to achieve a high degree of consistency of the AFs within each zone . These modifications to the zonation boundaries have been relatively minor, showing the geological zonation was a surprisingly good surrogate for site amplification, and in the transition from the V6 to V7 GMPM, it has remained unchanged.  Figure 15 shows the AFs obtained from the site response analyses at nine oscillator periods for Zone 1801, which has an average V S30 very close to the average over the entire field. The plots show that the response is highly non-linear, and for short oscillator periods, the linear response is also strongly dependent on magnitude (and, to a lesser extent, on distance).
The computed AFs are used to derive an intensity and scenario (magnitude and distance) dependent model. The selected model is given by  where Sa NS_B,g is the spectral acceleration at the NS_B horizon and is given in units of g (the acceleration of gravity), f 2 and f 3 are model parameters, is a standard normal random variable, lnAF is a parameter that represents the standard deviation of the data with respect to the median prediction of the model. The parameter f * 1 is magnitude-and distance-dependent and is given by the following:  Fig. 16. The standard deviation σ lnAF is allowed to vary with the input spectral acceleration (i.e. a heteroskedastic model) and is fitted using a trilinear function given by the following: where low , high , x l , and x h are model parameters. This represents the variability across the voxel stacks within a zone and across the range of input motions, but the full site-to-site variability also needs to account for the additional uncertainty due to uncertainty in the MRD curves, which were not varied in the site response analyses, and for the modelling uncertainty. For the former, the model of Bahrampouri et al. (2019) was implemented. The modelling uncertainty was estimated from comparisons between the empirical amplification factors determined from the borehole and surface recordings at the G-stations and the calculated AFs for the same 200 m interval. The resulting terms for the modelling error are slightly lower than those found by Stewart and Afshari (2021), for example, but this may be due to the site profiles in Groningen being more amenable to 1D site response modelling than some of the Californian sites analysed in that study. The total site-to-site variability is then given by the following: where Sa rock,low and Sa rock,high are model parameters and where lnAF,MRD is the additional uncertainty due to MRD for low and high intensity, lnAF,high and lnAF,low are the parameters of the model in Eq. (20), and model is the model error. Since this represents both the uncertainty in the V S profile at each voxel and the uncertainty related to the AF at any specific location relative to the zone average, both of which could be reduced by in situ measurements, the siteto-site variability in Eqs. (21) and (22) is modelled as epistemic uncertainty. This is added as a node in the logic tree (see Section 6) with three branches defined by samples of the distribution and associated weights that approximate the normal distribution (Keefer and Bodily 1983), following tests to confirm that the values conform to such a distribution.
For buildings located on dwelling mounds (locally known as wierden), an additional penalty factor is imposed on the site amplification factors (Kruiver et al. 2022a). For risk analyses using any spatially aggregated measure of damage or loss, models for the spatial correlation of the ground motions have also been developed Kruiver et al. 2022b).

Model performance
Total residuals have been calculated for the recorded motions with respect to the V7 GMM predictions at the ground surface . The fit to the data is generally very good, especially when considering that the model is developed by deconvolving motions to the NS_B using station-specific AFs but uses AFs assigned to broad zones in the forward predictions. The exception to this is the spectral acceleration at 0.01 s (which is equivalent to PGA), where the model significantly underpredicts the observed motions a consequence of calibrating the model to match the average spectral acceleration. At longer periods (0.7 s and greater), the model tends to slightly overpredict the recorded motions. These trends are of some concern, especially if the model is used to generate hazard maps of PGA, but the specific purpose of the GMM is for application in the calculation of seismic risk in the Groningen field. The parameter used to define the fragility functions is the average (geometric mean) spectral acceleration over the period range from 0.01 to 1.0 s. Figure 17 shows the residuals of the average spectral acceleration relative to the V7 GMM, from which it can be appreciated that the model performance is very good indeed. At the same time, it is important to emphasise that the primary purpose of the GMM is not to reproduce the observed ground motions but rather to provide reliable predictions of spectral accelerations over the range of earthquake magnitude that control seismic risk in the Groningen field. The range of earthquake magnitudes contributing to the risk estimates is effectively mutually exclusive with the range of magnitudes contributing to the database. The most important feature of the model, therefore, is not that 1 3 Vol:. (1234567890) Fig. 15 Calculated AFs from site response analyses performed for the voxels within site response zone 1801, grouped by the magnitude of the scenario used to generate the input motions it provides a perfect fit to the data-although such an outcome may be desirable-but rather that it estimates ground motions from potential larger earthquakes that reflect the source, path, and site characteristics of the field, while also capturing the epistemic uncertainty in these estimates.

Discussion and conclusions
Through an iterative process conducted over several years, with multiple iterations and rigorous review, ground-motion models have been developed for application to induced earthquakes in the Groningen gas field. These models include simple equations to predict PGV in the magnitude range of the events that have occurred in the field (up to M L 3.6) and a more elaborate GMPM for application to risk analyses that consider both larger induced events and even larger triggered earthquakes. To accommodate the epistemic uncertainty that is inevitably associated with the large extrapolations over more than three-and-half units of magnitude, the GMPM is formulated in a logic-tree framework. The complete logic tree for the model implementation is shown in Fig. 18. The ground-motion models are very well calibrated to the characteristics of the Groningen field by virtue of extensive and intensive data collection activities that have generated a large ground-motion database and detailed characterisation of the V S profile above the gas reservoir, with a particular focus on the near-surface layers of soft soils. While the general approaches that have been used could be adopted or adapted to other cases of induced seismicity, without such concerted efforts to gather similar datasets, there would inevitably be much larger epistemic uncertainty that would need to be accounted for in any ground-motion models. We believe that the models can be used with confidence for their intended purposes although their performance should be checked against new recordings generated by any future earthquakes in the field. There is probably little that can be done to improve the empirical GMPEs for PGV (which is based on an excellent distribution of data in the magnitudedistance to which it is applicable), although it would clearly be helpful to investigate the exact cause of the systematic differences between the motions recorded by the two different accelerograph networks. There may, however, be scope for improving or refining the GMPM for spectral accelerations. In particularly, the underestimation of PGA needs to be addressed for applications that use this parameter, and the minor overprediction at longer periods could also be usefully investigated. The alternative procedure could also be explored for the development of the model, for example basing the NS_B reference rock model entirely on Fourier amplitude spectra. To facilitate exploration of any of these issues by interested parties, the entire database, including all the V S , Fig. 17 Total residuals of average spectral acceleration (between 0.01 and 1.0 s) with respect to predictions of ground motions at the surface plotted against distance Fig. 18 Logic tree for the Groningen GMPM for response spectral ordinates damping and lithological profiles and the associated transfer functions and amplification factors have been made freely available ) so that others may now undertake the work of improving this model. Some readers may question the relevance of this work in light of the political decision by the Dutch government to close the Groningen gas field (e.g. Bommer 2022). There are several reasons, however, why seismic hazard and risk assessments, and hence the ground-motion models, for the field, remain relevant. Firstly, if production ceases, it is likely that the seismicity will continue for some time until the system reaches equilibrium. Secondly, the closure of the field has been delayed because the nitrogen-mixing facility required to render imported gas similar to the low-calorific gas from the field (to which all gas-burning applications in the Netherlands are calibrated) is behind schedule. And thirdly, in view of the primary source of imported gas being Russia, recent geo-political developments may result in the decision to close the field being reversed.