Developing a model for the prediction of ground motions due to earthquakes in the Groningen gas field

Abstract

Major efforts are being undertaken to quantify seismic hazard and risk due to production-induced earthquakes in the Groningen gas field as the basis for rational decision-making about mitigation measures. An essential element is a model to estimate surface ground motions expected at any location for each earthquake originating within the gas reservoir. Taking advantage of the excellent geological and geophysical characterisation of the field and a growing database of ground-motion recordings, models have been developed for predicting response spectral accelerations, peak ground velocity and ground-motion durations for a wide range of magnitudes. The models reflect the unique source and travel path characteristics of the Groningen earthquakes, and account for the inevitable uncertainty in extrapolating from the small observed magnitudes to potential larger events. The predictions of ground-motion amplitudes include the effects of nonlinear site response of the relatively soft near-surface deposits throughout the field.

Introduction

Induced earthquakes are occurring as a result of reservoir compaction in the Groningen field. These events are quantified by their epicentral locations and local magnitudes, M _L, determined by KNMI assuming a focal depth of 3km, which is the average depth of the Rotliegend sandstone in which the gas reservoir resides. Models have been developed for estimating possible future patterns of seismicity as a function of gas production (e.g. Bourne et al., 2014). To ascertain possible consequences of each future earthquake scenario on buildings and infrastructure at the ground surface, a model is required that predicts parameters characterising the shaking at the ground surface.

The main ground-motion parameters of relevance to the modelling of seismic risk in the Groningen field are 5%-damped elastic response spectral accelerations, SA, for a wide range of oscillator frequencies, from 0.2Hz to 100Hz, the latter being assumed equivalent to peak ground acceleration (PGA). Response spectral acceleration is generally a good indicator of the seismic demand on a structure – for which reason it is the parameter generally used in earthquake-resistant design – but for unreinforced masonry buildings the possibility of collapse can also depend on the duration of shaking (e.g. Bommer et al., 2004). Consequently, some of the building fragility functions for houses in Groningen have been defined as a function of both spectral acceleration and duration (Crowley et al., 2017), creating the need for predictions of the duration, D _S, as well. Additionally, since regulations regarding tolerable levels of vibration from anthropogenic sources are generally framed in terms of peak ground velocity (PGV), models have also been developed for this parameter.

The models predicting these parameters are a function of magnitude, M _L, and distance from the earthquake source. The models also account for influence of the near-surface geology and the amplification of shaking by soft layers of sands, clays and peats that are encountered in the field. Ground-motion models do not predict unique values of SA, PGA and PGV but rather distributions of values; the logarithmic residuals defining the distribution generally follow a normal distribution, characterised by the standard deviation, σ (e.g. Strasser et al., 2009). The value of sigma is as important as the values of the coefficients defining the median (μ) predictions. In any application, the number, ε, of standard deviations above the median must be specified; if this is zero, median values will be predicted, which have a 50% probability of being exceeded for the magnitude–distance scenario in question.

This paper provides an overview of the development of ground-motion models (GMMs) for SA, PGV and D _S in the Groningen field. Following this brief introduction, an overview is given of the available data and resource for this work together with a brief discussion of some of the specific challenges faced in this endeavour. The evolution of the models for SA, PGV and D _S is then summarised in the following sections, and the paper closes with a brief discussion of ongoing work and future refinements of the models.

Opportunities and challenges

Until quite recently, ground-motion modelling has been focused almost exclusively on the estimation of ground motions from large tectonic earthquakes in order to provide inputs to seismic design and loss assessment. Extrapolating such standard GMMs to small-magnitude, shallow-focus induced earthquakes is problematic, requiring the development of bespoke models.

Opportunities in Groningen

The work described in this paper has benefited from excellent databases available in the Groningen field. There is detailed characterisation of the shallow geology in the field (Kruiver et al., 2017a), which, in combination with numerous measurements, has enabled the development of reliable shear-wave velocity profiles that extend from the ground surface to the base of the North Sea supergroup at some 800m depth (Kruiver et al., 2017b). There are also extensive networks of seismic recording instruments in the field (Fig. 1), which have yielded a rich database of accelerograms (Fig. 2); detailed information about the recording networks is provided by Dost et al. (2017). The B-stations refer to surface accelerographs installed and operated by KNMI, which have been operating for many years and have recorded all significant earthquakes in the field. KNMI also operates the G-stations, installed by NAM, which are boreholes with geophones installed at 50, 100, 150 and 200m together with a surface accelerograph. As can be appreciated from Figure 2, the addition of the G-stations has greatly increased the number of recordings captured in each earthquake. The database used to derive the most recent models includes 178 recordings obtained at distances of up to about 25km from 22 earthquakes with M _L≥2.5.

Fig. 1. Strong-motion recording instruments in the Groningen field. The coordinates are given in metres in the Dutch RD system.

Fig. 2. Histogram of accumulation over time of accelerograph recording from the B-stations (red) and G-stations (blue); magnitudes are indicated at the top of the diagram

In addition to these two networks, NAM (Nederlandse Aardolie Maatschappij B.V.) has funded the installation by TNO of some 300 instruments in public buildings and private homes (grey squares in Fig. 1). While these accelerographs have other objectives, the recordings are being analysed for potential inclusion in the database for the GMM development. One issue being addressed is that not all of these instruments are installed at ground level and the degree of influence of the structural response needs to be determined. Another issue is that the instruments are set with a trigger threshold defined by PGV values of 1mms⁻¹; recordings of lower amplitude are excluded, with the result that the recordings may be biased high with respect to those from the KNMI network (Fig. 3).

Fig. 3. Geometric mean horizontal PGV values from the KNMI and TNO accelerographs from the M_L 3.1 2015 Hellum earthquake, showing the censoring effect of the trigger threshold on the household instruments. The green dots are PGV values reported as part of the ‘heartbeat’ monitoring that transmits the highest velocity value recorded in each minute; the time-histories corresponding to these records are not retained.

Another notable resource from which this work has benefited enormously is the body of expertise and experience on which it has been possible to draw. This has included, for example, teams at both Shell and ExxonMobil performing full waveform simulations using finite difference techniques with 3D velocity models for the field. Expert reviews of the model development have also been of great benefit in terms of providing critical feedback and constructive suggestions. In particular, a two-day review workshop was held in October 2015 with the participation of several of the leading names in the fields of ground-motion modelling and site response analysis: Gail Atkinson (Canada), Hilmar Bungum (Norway), Fabrice Cotton (Germany), John Douglas (UK), Jonathan Stewart (USA), Ivan Wong (USA) and Bob Youngs (USA). As noted below, another very constructive workshop was held subsequently on the specific topic of finite fault rupture modelling.

Challenges in Groningen

One of the key challenges in developing a GMM incorporating local site effects for Groningen has been the scale. The study area for risk modelling is the entire onshore portion of the gas field plus a 5km buffer, yielding an area of about 1000km². The site response model is possibly the largest seismic microzonation that has ever been developed at this level of resolution.

At a very early stage of the work it became apparent that a reliable GMM would need to be developed specifically for Groningen rather than imported from any other region. The main reason is the unusual upper crustal profile in the region and specifically the high-velocity Zechstein salt layer situated immediately above the reservoir leading to refraction and reflection of seismic waves (Kraaijpoel & Dost, 2013). This effect was very clearly manifested when comparing the available recordings of Groningen ground motions with the predicted PGA and PGV values from the equations of Dost et al. (2004), which were derived from recordings of induced earthquakes in the Roswinkel gas field: as can be seen in Figure 4, the equation severely over-predicts the recorded peaks. Whereas in Groningen the Zechstein is located above the gas reservoir, in the Roswinkel field the salt layer is below the gas reservoir.

Fig. 4. Residuals of recorded PGA (left) and PGV (right) values from the Groningen field relative to predictions from the Dost et al. (2004) equations derived from recordings of induced earthquakes in the Netherlands. Negative residuals imply over-prediction of the observations.

Another obvious challenge relates to the need to predict ground motions from earthquakes far larger than those from which recordings are available in the Groningen database. As the outcome of an expert workshop held in March 2016, a distribution of maximum earthquake magnitudes was defined that allows for the possibility of events of up to ~7 (e.g. Bommer & van Elk, 2017). Even though the probability assigned to such earthquakes being possible is very low, the GMM must necessarily be calibrated to provide reliable predictions up to this level, an extrapolation of more than three magnitude units beyond the data.

Evolution of spectral acceleration models

For the seismic hazard calculations performed to support the 2013 Winningsplan, a preliminary GMM was developed only for PGA and PGV. This model adopted the European equations of Akkar et al. (2014) for larger magnitudes, assuming a field-wide 30m time-average shear-wave velocity (V _S30) of 200ms⁻¹. For magnitudes lower than ~4, the coefficients of the equations were adjusted so that the equation matched the local data; the derivation is described in Bourne et al. (2015). This conservative model – referred to as the V0 model – derived as a stopgap measure, was subsequently superseded by four phases of development of the model, as required by the regulatory review process. Table 1 summarises the key developmental steps. In each of the developmental stages, the objective was to produce a robust and stable model incorporating one or more advances towards the final goal. In the following sections, the key developments are discussed with respect to the incorporation of site amplification effects and the prediction of motions at the underlying rock horizon.

Table 1. Summary of the evolution of models for the prediction of the amplitudes.

Site response

The V1 model, the first Groningen-specific equation, was obtained through point-source simulations based on inversions of the Fourier amplitude spectra (FAS) of surface recordings (Bommer et al., 2016). The model captured characteristics of the field but in making predictions directly at the ground surface it was only able to model an average field-wide site amplification function and linear site response. The V2 model represented the first important refinement of the site response modelling, by predicting motions at a buried horizon and then convolving the predicted motions with nonlinear site amplification factors derived for the overlying soil layers. In V2, the reference rock horizon was the base of the Upper North Sea formation (NU_B) located at about 350m; in the V3 model, this was moved to the base of the North Sea supergroup (NS_B) at ~800m depth since it corresponds to a much clearer impedance contrast. Using the velocity profiles from NS_B to the surface (Kruiver et al., 2017b) and input motions at this horizon, nonlinear site amplification factors were derived for all spectral accelerations (Rodriguez-Marek et al., 2017). Starting from a geological zonation of the field (Kruiver et al., 2017a), adjustments were made to the boundaries such that each zone represented a relatively consistent behaviour in terms of the amplification factors (Fig. 5).

Fig. 5. The 160 zones for site amplification factors in the V4 model.

In the V4 model, as well as being functions of oscillator frequency and the amplitude of motion at the NS_B horizon, the amplification factors at short periods are also dependent on magnitude and distance (Stafford et al., 2017).

Rock horizon motions

From the V2 model onwards, the first stage of the model construction has been to deconvolve the surface recordings to the reference rock horizon, assuming linear site response because of the low amplitudes of motion (Bommer et al., 2017). Inversions are then performed on the FAS to estimate possible ranges of the Brune stress parameter (Δσ) and Q; kappa is estimated directly from the FAS and the geometric spreading function partly constrained by the finite difference simulations. Optimal values of the source, path and site parameters are then found to minimise the misfit to the SA values at rock. In forward modelling, however, alternative values of Δσ (and magnitude scaling of this parameter) are considered in order to capture the epistemic uncertainty at larger magnitudes. The uppermost branch of the logic-tree for the rock motions is now calibrated to predict motions consistent with those from tectonic ground-motion prediction equations (GMPEs) for the same magnitude–distance combinations.

In the V1 model, only a small number of oscillator periods were considered for convenience. Several periods were added at the V2 stage to cover all possible requirements in terms of structural analyses and fragility functions; the additional periods added subsequently are only at the high-frequency end of the spectrum and were included to facilitate the application of vertical-to-horizontal ratios to obtain vertical response spectra (e.g. Bozorgnia & Campbell, 2016).

Up to V3, the models were based on point-source simulations with the distance being measured relative to the epicentre (since the depth is a constant 3km), and in hazard calculations earthquakes were represented as hypocentres. This is unrealistic for larger earthquakes for which fault rupture dimensions may be several kilometres. A workshop was held in July 2016 to discuss options for introducing extended fault ruptures, with presentations from Christine Goulet of the Southern California Earthquake Center (SCEC) on the benchmarking of fault rupture simulation models. Norm Abrahamson and Bob Youngs presented the experience on using finite fault simulations in GMM-building projects such as NGA-East (http://peer.berkeley.edu/ngaeast/), concluding that the most robust and reliable way to proceed in terms of stable results for the full range of frequencies of interest was to use stochastic simulations. For deriving the V4 model, the program EXSIM (Motazedian & Atkinson, 2005) was selected, in which a finite fault rupture is approximated by a series of point-sources distributed over a plane and initiated at intervals to mimic the rupture propagation.

Figure 6 shows examples of predicted response spectra at the NS_B rock horizon and at the ground surface assuming both linear and nonlinear site amplification. In the V4 logic-tree there are four branches for median predictions of spectral acceleration.

Fig. 6. Predicted median response spectra from the four branches of the V4 model and one of the site amplification zones, for a relatively weak (left) and a stronger (right) earthquake scenario.

Aleatory variability

The variability (σ) associated with the predictions of motion at the NS_B rock horizon is decomposed into between-earthquake (τ) and within-earthquake (φ) variance components. The former is estimated from the residuals of the observed motions and is found to be smaller than values associated with most modern GMPEs for tectonic earthquakes; this is not a surprising result given that the Groningen earthquakes all essentially have the same seismic source. While many modern GMPEs model smaller values of between-earthquake variability at larger magnitudes, it was decided to maintain the value constant in extrapolations to larger magnitude. Although φ can also be estimated from the data, the adopted value is smaller since the site-to-site component of this variability is accounted for in the site amplification model. This non-ergodic, or single-station, value of φ is inferred from other studies that have found this quantity to be stable over different regions (Rodriguez-Marek et al., 2014; Al Atik, 2015).

There is additional variability associated with the site amplification factors, which is based on the variability in the calculated factors over each zone and for the full range of input rock motions. This site response variability within each zone is referred to as φ _S2S.

Implementation of the complete GMM is illustrated in Figure 7. For any earthquake scenario, a fault rupture plane – of dimensions consistent with the magnitude and propagating laterally and downwards from the reservoir – must be defined, from which the rupture distance to any location on the NS_B horizon can be calculated. Spectral accelerations at those locations can be estimated from the predictive equations, sampling the same between-event residual at all locations and randomly sampling from the within-event variability. In the current hazard and risk models, spatial correlation is not explicitly modelled since uniform motion is considered within zones (or sub-areas for the larger zones), hence yielding perfect correlation within the zones and no correlation from zone to zone, approximating more realistic spatial correlation models. The SA value at NS_B is then input to the amplification function and the amplification factor, AF, calculated by also sampling from the site-to-site variability. The final surface motions are obtained multiplying the SA value at NS_B by the AF value.

Fig. 7. Schematic illustration of the calculation of SA at three surface points, in two zones, for an earthquake of magnitude M_a and an event-term of ε_bτ; in this simple example, the within-event variability is sampled without considering spatial correlation.

PGV predictions

The V0 model did include an adaptation of the Akkar et al. (2014) equation for PGV, but subsequently this parameter was not modelled since it was not employed in any of the fragility functions related to estimating partial or complete collapse of buildings (Crowley et al., 2017). Moreover, PGV could always be estimated from the spectral acceleration at a response period of 0.3s since these two parameters are very well correlated in the Groningen data (Fig. 8), which is consistent with other studies that have shown such proportionality between SA and PGV, with the controlling period being a function of the earthquake magnitude (Bommer & Alarcón, 2006).

Fig. 8. Correlation between PGV and SA (0.3s) values from Groningen data.

Nonetheless, the potential use of PGV in estimation of lower levels of structural damage and other applications led to this parameter being modelled. Concurrent with the V3 model, an empirical prediction model was derived from the Groningen data; initially V _S30 was included as an explanatory variable but found to exert very little influence – given the modest variation in this parameter over the field – so the equation was a function of M _L and R _epi only. In the development of the V4 model, PGV was included as an additional predicted parameter and modelled in the same way as the 23 spectral accelerations. Figure 9 shows predicted median surface values of PGV from both the empirical model and the four equations of the V4 model as a function of magnitude and at two rupture distances (hypocentral distance is used for the empirical model). The plot shows, as expected, general consistency between the V4 and empirical equations in the magnitude range of the data, but it also serves to illustrate why simple extrapolation from the data yields inappropriate results, vindicating the decision to use simulations.

Fig. 9. Comparison of median predictions of surface PGV from the four equations of the V4 model and the empirical equation derived by regression on the local data

Duration predictions

The Groningen recordings show very short durations – even for such small-magnitude earthquakes – at short epicentral distances, which rapidly increase with distance due to the refraction and reflection of seismic waves by the overlying high-velocity salt layer (Bommer et al., 2016; Fig. 10). While stochastic simulations provide a powerful tool for extrapolating predictions of ground-motion amplitudes using seismological theory (Boore, 2003), they are not useful for predicting duration since they actually require duration estimates as an input. Consequently, in earlier phases of the model-building, attempts were made to adjust predictive equations for duration, such as those of Kempton & Stewart (2006) and Afshari & Stewart (2016a), to the Groningen data. The results were not very satisfactory because of the strong adjustments necessary to match the Groningen durations at small magnitudes and the strong distance dependence close to the epicentre and the lack of constraints at larger magnitudes and distances.

Fig. 10. Comparison of two accelerograms – obtained at the MID1 and WIN recording stations (Dost et al., 2017) – from the M_L 3.6 Huizinge earthquake in 2012 and their Husid plots showing the accumulation of energy over time; although the WIN record is from an epicentral distance just 4km greater, the duration, based on 5–75% accumulation of the total Arias intensity, is more than seven times longer.

In each stage of development, the duration model from the previous stage was used to define the durations required for the stochastic simulations. At the V4 stage, this same procedure was followed, but now the limitation of the stochastic modelling approach was diminished because the durations are only input to the small-magnitude events – for which the duration model is constrained by the data from the field – representing patches of the fault rupture modelled in EXSIM. The durations from larger earthquakes are then obtained from the resulting time-histories generated at the NS_B horizon as a result of the interactions of the signals from the individual patches. Clearly the predicted durations therefore depend on the assumed rupture velocity and aspect ratio, but the basis for these estimates of the shaking duration are more physical than in previous models and consistent with the geophysical properties of the field. The durations predicted at the NS_B horizon are transformed to the ground surface via terms adapted from the Afshari & Stewart (2016a) model to match the local recordings; the site term is defined as a function of V _S30, for which reason a map of median V _S30 for each zone has been prepared (Bommer et al., 2017). Figure 11 compares predictions from the V4 duration model – in which there are four branches, each twinned with one of the branches for median spectral accelerations – with those from Afshari & Stewart (2016a). For small magnitudes – where all the V4 branches are equivalent – the V4 model captures the shorter durations close to the epicentre and the more rapid increase of duration with distance than modelled by the Afshari & Stewart (2016a) model, which is derived from recordings of tectonic earthquakes. For larger magnitudes, the predictions all show similar trends with distance, but the Groningen model yields durations almost twice as long as those from Afshari & Stewart (2016a). At this stage, it is not clear whether this is a genuine feature of motions from larger earthquakes in the gas field or simply the tendency for EXSIM to overestimate durations from empirical equations, as reported by Afshari & Stewart (2016b).

Fig. 11. Predicted median durations as a function of distance from earthquakes of magnitude 3.5 and 6.5, using the four branches of the V4 model and the equation of Afshari & Stewart (2016a). For M_L 3.5, all of the V4 branches are identical.

In applications to the risk modelling, the durations are predicted as conditional on the spectral accelerations using the correlation coefficients between residuals proposed by Bradley (2011), which have been found to be generally consistent with the trends observed for Groningen motions.

Conclusion

A model for the estimation of ground-motion amplitudes and durations due to induced earthquakes in the Groningen field has been developed through an evolutionary process. The framework of the model is now stable in terms of prediction of motions at a reference rock horizon (the base of the North Sea Supergroup at 800m depth) using extended ruptures to represent earthquake sources, and convolution of these rock motions with nonlinear, frequency-dependent amplification factors to obtain the motions at the ground surface. The model is consistent with available data for the Groningen field but also captures the epistemic uncertainty associated with extrapolations for the larger earthquake scenarios considered in hazard and risk estimation. Work is ongoing to further refine the model, but future changes are now expected to be incremental. The V4 Groningen GMM can be used to estimate losses due to potential future earthquakes in the gas field and to define inputs to earthquake-resistant design of structures.