Development of seismicity and probabilistic hazard assessment for the Groningen gas field

Abstract

The increase in number and strength of shallow induced seismicity connected to the Groningen gas field since 2003 and the occurrence of a M_L 3.6 event in 2012 started the development of a full probabilistic seismic hazard assessment (PSHA) for Groningen, required by the regulator. Densification of the monitoring network resulted in a decrease of the location threshold and magnitude of completeness down to ~ M_L=0.5. Combined with a detailed local velocity model, epicentre accuracy could be reduced from 0.5–1km to 0.1–0.3km and a vertical resolution ~0.3km. Time-dependent seismic activity is observed and taken into account into PSHA calculations. Development of the Ground Motion Model for Groningen resulted in a significant reduction of the hazard. Comparison of different implementations of the PSHA, using different source models, based on either a compaction model and production scenarios or on extrapolation of past seismicity, and methods of calculation, shows similar results.

Introduction

Since 1991, induced earthquakes have occurred in the province of Groningen that are linked to gas production from the Groningen gas field. A monitoring system for induced seismicity in the north of the Netherlands became operational in 1995 (Dost & Haak, 2007). This network was designed to detect and locate earthquakes of magnitude 1.5 and larger. The reason for this threshold is the fact that people in the region reported felt tremors starting around M 1.8. The Groningen gas field remained at a low activity rate until 2003. In this time period, larger events (M 3–3.5) occurred in smaller onshore gas fields (e.g. Roswinkel and Bergermeer).

Since 2003 seismic activity has increased in Groningen, coinciding with an increase in production of the field. Due to the small number of events recorded in the Groningen field, initially seismicity related to all gas fields in the north of the Netherlands was combined to gain insight into the characteristics of seismicity. Van Eck et al. (2006) showed differences in frequency–magnitude relation between Groningen and the other gas fields and published a first probabilistic seismic hazard assessment (PSHA) for the region. Output of the PSHA is given in the form of a map of predicted ground motions for a specified probability of exceedance, which is selected at 10% probability in 50 years or a return period of 475 years.

Further development of the PSHA for Groningen proceeded along two lines. The first approach consists of a seismic source model based on compaction of the reservoir (Bourne et al., 2014) combined with a Monte Carlo approach to the hazard calculations (Bourne et al., 2015). The source model allows comparison of the effects of different production scenarios on seismic activity. We will refer to this model as the NAM model. The second approach, followed by the KNMI (Royal Netherlands Meteorological Institute), uses a seismic source model based on recorded seismicity and performs an integration over seismic zones. The KNMI seismicity model is stationary, but uses a limited time frame of 5 years to calibrate (Dost & Spetzler, 2015). This time frame was selected as a compromise between (1) a short duration to minimise the effect of variations in activity rate and (2) the availability of a dataset large enough to calculate statistical parameters for each seismic zone. In the NAM model, annual activity rates are calculated for different production scenarios (e.g. NAM, 2016a). Observed and simulated annual event counts show a good correspondence within the 95% confidence interval, and for the production scenarios most relevant for the next 5 years (21 and 27bcm (billion cubic metres)) a nearly constant level of the median annual event rate is predicted for the next 5-year period. Averages over a 3- to 5-year period vary smoothly. Both methods share the same Ground Motion Model (GMM; e.g. Bommer et al., 2017), and result in hazard maps expressed in ground acceleration. In this paper we give an overview of the development of the monitoring network, seismicity, earthquake location and seismic hazard analysis.

Development of the monitoring network

The high-noise conditions in the north of the Netherlands prevent the detection of small events in the region. Therefore a borehole string was designed to improve the signal-to-noise ratio (Dost & Haak, 2007). For noisy sites the improvement is up to 20–30dB in the 1–20Hz band, while for quiet sites it may be limited to 10dB. Figure 1A shows the development of the network. The first pilot borehole, station FSW, was constructed in 1991, near the village of Finsterwolde and just east of the Groningen field. This station consists of five levels of 4.5Hz geophones with a 75m spacing (Fig. 2). Lessons learned from this station were used to build a regional network in 1995. It was designed to cover a 40×80km region with an average station spacing of 20km (Fig. 1A). For the seven new borehole stations, the vertical geophone spacing was reduced to 50m and the surface sensor was omitted (Fig. 2). In 2010 three borehole stations were added to increase the aperture of the network. These stations (NIW, SUH, SPY) have a 30m vertical spacing between sensors.

Fig. 1. Overview of the development of the permanent seismic network around the Groningen gas field: (A) borehole network and (B) accelerometer B-network. Station names are indicated, except for the G-network (G01–G70). The fill colour of triangles and squares denotes the year of installation of the geophones and accelerometers, respectively. The blue line depicts the border of the Groningen field. Coordinates are shown in the Dutch National Triangulation Grid (Rijksdriehoekstelsel).

Fig. 2. Overview of five borehole station configurations in the Groningen area. Green triangles denote depth levels of 4.5Hz geophones, green squares denote accelerometers. Station names are depicted on top of the borehole configurations, years of installation on the bottom.

In 2012 the largest event in the region was recorded, M _L=3.6 near Huizinge, causing damage to structures in the province of Groningen. After this event it was decided to extend the monitoring network with the aim of improving earthquake location accuracy and lowering the detection and location threshold. A total of 70 stations were added in 2014–2016 to cover the Groningen field at an average station spacing of 4–5km with the aim of improving hypocentre accuracy and source mechanism determination. Detailed velocity models are available for the first 4km of the crust, allowing the application of new location methods. Multiple levels of sensors are used to discriminate between P and S arrivals and also for redundancy in case of an instrument failure, since the boreholes are non-cased. These new stations are similar in their set-up to the 1995 ones, with a surface accelerometer added at the ground surface (Fig. 2). Together they form the G-network. In 2016 another three stations were installed near the village of Norg to monitor a gas-storage facility southwest of the Groningen field. These stations are coined the N-network and are also used for constraining event locations in the Groningen field. As of the end of 2016, a total of 337 geophones are in place in the Groningen region.

In addition to the borehole network, a surface accelerometer network was installed in the region (Fig. 1B). The locations of the accelerometers were planned where felt events were reported. Data are used in the evaluation of damage and for the development of the GMM. Initially, in 1997, five accelerometer stations were installed in the central-north part of the field near Middelstum. The accelerometer network gradually developed over time. In 2017, 17 accelerometers are in operation apart from the borehole network. Together with the borehole accelerometers from the G- and N-networks, a total of 90 accelerometers are in place.

Though all the sensors are in place, the G-network is not completely operational at the beginning of 2017. Station G15 could not be connected to the electricity grid. For most other stations, real-time communication is only possible through wireless communication over a 3G network, since a fixed connection through a digital subscriber line (DSL) could not be realised in a short time frame. However, some parts of Groningen have a poor 3G signal, resulting in varying data availability of some of the stations. Figure 3 shows the temporal development of data completeness. Installation began in the end of 2014. Over most of 2015 there is an increasing trend as more and more stations are built and connected to 3G. Not until the beginning of 2017 did the completeness rise further when a part of the stations was connected to DSL. This rising trend is expected to continue into 2017 when the remaining stations are connected.

Fig. 3. Temporal development of data completeness for the G-network (Fig. 1A). 100% completeness corresponds to continuous data being available from all 69×5 sensors.

Detection and location thresholds were calculated based on measured average noise levels in the borehole sensors at 200m depth and the attenuation relation derived for the region used to calculate local magnitude (Dost et al., 2012). For the location threshold, a minimum of three stations contributing phase readings are required. Figure 4 shows the location threshold for the original and the extended network. For Groningen, the location threshold decreased from M _L~1.0–1.5 to M _L~0.5. Until 2010 all stations were operating in a triggered mode and only events were selected and transmitted to the data centre. Since 2010 instruments have been connected to the KNMI data centre in real time, and continuous data transmission is realised. Data can be accessed at http://rdsa.knmi.nl/dataportal. In the full SEED volumes metadata is available and an overview of stations is available at www.knmi.nl/nederland-nu/seismologie/stations.

Fig. 4. Location threshold for the national seismic network in the Netherlands before (left) and after (right) the network upgrade.

The orientation of the horizontal geophone components is unknown by the time the sensors are lowered down the hole. For most pre-2014 geophones, the orientations were determined later using checkshots at short distances. Orientation of the new network is being determined using checkshots and through cross-correlation between the geophones and the co-located surface accelerometer (Hofman et al., in press).

All stations mentioned thus far are part of the national network operated by the KNMI (network code NL). Additionally, two broadband stations exist in the Groningen area, which are part of the NARS network run by Utrecht University. Station NR.NE117 is co-located with NL.FSW; the other station (NR.NE018) is co-located with NL.VLW (Fig. 1A). In 2017 the NL network will be further expanded in the Groningen area by installing four broadband sensors at 100m depth.

The Dutch technology institute TNO constructed another seismic network in the Groningen area. This accelerometer network became operational in 2015, specifically to study the seismic response of different buildings in Groningen. It encompasses a total of 300 accelerometers, 280 of which reside in private dwellings and the remaining 20 in public buildings.

In addition to the networks mentioned above, various geophone strings have been installed at reservoir depth (NAM, 2016c). This deep-geophone network is run by the operator of the Groningen field (NAM). Two strings have been placed in former production wells near the villages Zeerijp and Stedum (stations SDM-1 and ZRP-1). In 2016, these two geophone strings were replaced with equipment lowered into dedicated monitoring wells (stations ZRP-2 and ZRP-3). These new strings have 15 geophones each, with sensor spacing varying between 25 and 65m. In 2016 another borehole geophone string was installed below the village Harkstede (station HRS-2A) to monitor seismic activity near the city of Groningen.

Earthquake location

The accuracy of the locations of earthquakes strongly depends on the geometry of the station network and knowledge of the velocity structure in the subsurface. The layout of the first borehole network required the use of an average 1D velocity model for the north of the Netherlands. The resulting accuracy of event locations was 0.5–1km in the horizontal plane. Resolution in depth was even more limited (e.g. Bourne et al., 2015). In practice, event depth was fixed in the location procedure at 3km for Groningen, the average depth of the reservoir. Location accuracy could be improved by (1) decrease in average station distance and (2) introduction of a local velocity model. The hypocentre locations are published on www.knmi.nl/nederland-nu/seismologie/aardbevingen.

The upgrade of the network in 2014–2016 and the availability of a detailed 3D velocity model for Groningen, provided by NAM, allowed for the implementation of other source location methods. As an alternative to the standard HYPOCENTER method (Lienert et al., 1986), the EDT method by Lomax (2006) was used (Spetzler & Dost, 2017) and the resolution in epicentre location improved to 0.1–0.3km. Also uncertainty in depth has been reduced to about 0.3km, making it possible, for the first time with the shallow borehole network, to distinguish between events within and outside the reservoir. For relocated events recorded since 2014 it was found that 90% of the events take place within the reservoir or the top of the underburden, with the possibility of a few weak events (M _L<1.0) whose locations correlate with the location of anhydrite layers within the Zechstein layer above the reservoir. In a first analysis of micro-earthquakes, recorded in two deep downhole tools (ZRP-1 and SDM-1), Pickering (2015) also concluded that seismicity is confined to the reservoir. The hypocentre of micro-seismic induced events in the central north (Loppersum) area, located using data from the deep-geophone stations SDM and ZRP, is published on the NAM platform webpage (www.namplatform.nl/feiten-en-cijfers.html).

The subsurface structure of the north of the Netherlands is characterised by a thick deposit of Zechstein evaporates. The Groningen field is overlain by the Zechstein layer, a high-velocity layer, varying in thickness between a few hundred metres and about 1.5km. Kraaijpoel & Dost (2013) investigated the influence of the subsurface in Groningen on recorded waveforms and combined recordings from the sparse borehole network with a local accelerometer network. Two effects were observed that are specific for Groningen: (1) defocusing of seismic energy and (2) relatively strong conversions from P- to S-energy at the bottom of the Zechstein layer, around 2800m depth, leading to S-wave precursors and the presence of many multiples between direct P and direct S. The latter effect results in a complex wave train. These observations, together with the large uncertainties in the S-wave velocity model, led us to limit the hypocentre inversion procedure to an inversion of the P-waves only.

The 3D velocity model provides detailed information for the upper 3–4km of the crust. The velocity structure of the deeper part of the Carboniferous layer below the reservoir is much less known. The upper few hundred metres of the Carboniferous is characterised by a linear velocity gradient (Spetzler & Dost, 2017), but it is yet unknown to what depth this gradient continues. The first-arriving seismic waves sample the deeper part of this layer at epicentral distances larger than 10km. A more detailed investigation of the deeper velocity structure is expected to further improve earthquake locations in the Groningen gas field.

Locations of seismicity prior to the densification of the network could be updated using post-2014 events as reference events and applying a double-difference algorithm (Waldhauser & Ellsworth, 2000). The relocation was tested for a subset of the seismicity (e.g. Jagt et al., 2017). In order to update event locations, events with high similarity in waveforms, representing rupture from the same fault, or nearby parallel fault segments, are needed. For the period analysed, 2010–2014, it turned out that only a fraction of the events were clustered and could be relocated successfully. By the end of 2016, the cluster analysis has been expanded to cover the entire catalogue (1991–2017). Many more clusters are found, with large time lags between repeating events.

Development of seismicity

Seismicity in the Groningen field is thought to be caused by an increase in total strain due to compaction of the producing gas layer (Bourne et al., 2014). This means that production changes will influence observed seismicity and the earthquake generation process is non-stationary. Both temporal and spatial variations in seismicity have been observed.

Temporal variations

Based on a small dataset of 179 events recorded in the period 1991–2004 and associated with the Groningen field, Van Eck et al. (2006) show the difference between the frequency–magnitude (FM) relations for the Groningen field and all other gas fields in the north of the Netherlands. The slopes of the curves (the b-values) are similar, while the activity rate (a-value) for Groningen, expressed as the cumulative annual frequency, varies in time. At the end of 2016 a total of 1035 events have been recorded in Groningen, 279 events of M _L≥1.5 (Fig. 5A). From an assessment of the cumulative seismic energy release, Dost & Kraaijpoel (2013) showed a change around 2003, also reflected in the FM relations for the period before and after the change. Nepveu et al. (2016) investigated a possible change in activity rate using a Bayesian change-point model and concluded that a change occurred in January 2003. This change could not be related to any physical changes in the reservoir.

Fig. 5. (A) Annual seismic activity in Groningen as a function of time. (B) Annual cumulative frequency as a function of magnitude for the Groningen gas field for 1996–2003 and 2003–2017.

Figure 5B shows an update of the FM relation for Groningen, comparing the periods before and after January 2003. The b-values for the two periods were calculated using the maximum likelihood method of Tinti & Mulargia (1987). Marzocchi & Sandri (2003) showed that this method is simple and also has no significant biases for small datasets. A small difference in b-values is observed, from 1.06±0.07 for the first period to 0.93±0.03 for the second period. The number of earthquakes in these periods is 99 and 887 respectively. Activity rate increased from 2.16 to 2.60 and is expressed as the logarithm of the cumulative annual number of events. The range of measured b-values is similar to results at other gas fields, e.g. Lacq where b-values vary between 0.6 and 1.1. (Volant et al., 1992).

Since 2010 a total of 80–100 events per year are recorded. This allows a more detailed look at changes in the FM relation calculated over non-overlapping intervals of ~2 years. Time periods were chosen to contain around 200 events. Results are shown in Table 1. For the period 2010–13, the b-value is slightly higher than for the other period and similar to the period before 2003, but for all other intervals remains around 0.9. Activity rate, however, shows a rapid increase since 2003 and a drop in the last time period. Changes in production were made in early 2014, imposed by the government.

Table 1. Measured FM parameters and cumulative production for the listed time periods.

In a statistical assessment of the effects of production shut-in in the Loppersum area in 2014, Paleja et al. (2016) concluded that activity rate in this area was reduced with respect to the year preceding the change in production. Nepveu et al. (2016) came to a similar conclusion and proposed a delay of 2–8 months between the change in production and the change in activity rate.

Spatial variations

The fault structure related to the hydrocarbon reservoir is well known (e.g. Van Gent et al., 2009) and, although location accuracy is limited and therefore events occurring before mid- 2014 cannot be uniquely connected to individual faults, seismicity seems to follow the regions of high fault density (Fig. 6A). Most seismic energy has been released in the Loppersum area in the northwest of the field (Fig. 6B), where there are primarily northwest–southeast-striking faults. Also the southwestern part of the field, with similar fault orientations, has been quite active. There has been little energy release in the northeastern part of the field, where fault density is smaller than in other parts of the field. Almost no activity has been found in the southeastern part of the field, where fault orientations are primarily north–south and east–west

Fig. 6. Seismicity development over the Groningen field. (a) Earthquake epicentres of detected seismicity from 1991 (dark green dots) to the end of 2016 (yellow dots), 0.5<M_L<3.6. In the background the outline of the Groningen field (blue line) and faults at reservoir level (grey lines) are shown. (B) Spatial distribution of the total seismic moment release from 1991 until 2016. (C) Development of the centre of gravity of the seismic moment from 2010 until 2016 (dots). The radius of the circle for every year is scaled with the root of the seismic moment for that year. Lines connect the centre of gravity from year to year. (D) Seismic moment release in 2016 in relation to gas production. The red circles are scaled with the total production for different production units. Red arrows point at the production units that have been largely shut in since the beginning of 2014.

In order to gain insight into the spatial variation of seismic moment, the centre of gravity of the annual seismic moment was calculated and is plotted in Figure 6C. After a strong concentration in the Loppersum area, there is a decreasing and southward trend of the seismic moment from 2013 until 2016. Figure 6D shows the seismic moment release in 2016 combined with information on production. Seismic moment release has largely dropped in the Loppersum area since the end of 2013. However, the area remains somewhat active, as evidenced by the seismic moment distribution in 2016.

Bourne et al. (2014) showed a possible dependence of the b-value on compaction, leading to a lower b-value at higher compaction. The b-values for Groningen vary between 0.8 and 1.4. Since 2014 seismicity shifted from a region with high compaction in the central-north (Loppersum) region to regions of lower compaction south of Loppersum, an increase in b-value was expected. This has not been observed in our present results (Table 1).

Magnitude

Shortly after the first borehole network was installed, empirical relations were derived to determine the attenuation in the region, necessary to obtain a local magnitude (Dost et al., 2004). These relations are based on recorded amplitudes at 200m depth on the boreholes and calibrated with respect to the national network for events larger than M _L=3. At magnitudes important to PSHA (M>2.5) it was assumed that local magnitudes and moment magnitudes are equal.

Later, moment magnitudes for Groningen events were estimated from the spectra of surface accelerometer recordings. Since most recordings are at short distances (<10km), effects of geometrical spreading and damping play an important role. First results on the scaling between M _w and M _L indicated a relation M _w=M _L–0.2 (Dost et al., 2016). At this moment a reassessment of these results is ongoing.

Maximum magnitude

An estimate of the maximum magnitude for induced events in the Netherlands was based on (1) modelling of the frequency–magnitude relation, assuming a double truncated Gutenberg–Richter model, and (2) a physical model based on estimates of available fault surface (e.g. Dost et al., 2012). Bourne et al. (2014) show that there is no statistical evidence for an upper bound of the FM curve and argued for a maximum magnitude, 6.5, based on the seismic release of all strain accumulated over the life cycle of the field at the end of production in a single event. Dost et al. (2013) discussed the selection of the M _max for Groningen, concluded that M _max could not be determined based on statistics alone and adopted M _max=5. This value came from a literature study on the maximum observed magnitudes of induced events at gas fields worldwide. The most recent statistical analysis for Groningen was published by Zöller & Holschneider (2016), who proposed a M _max=4.4. Both statistics and modelling did not provide a unique value for the M _max for Groningen. Therefore, a workshop was held on the issue of the M _max for Groningen lead by a panel of international experts. In their report the panel proposed a distribution of M _max values, peaked at M _max= 4.5 and based on expert judgement (Bommer & Van Elk, 2017).

Magnitude of completeness

The FM curves in Figure 5 show a deviation from the linear part of the curve towards the lower magnitudes. The lowest magnitude at which all events in the field could be observed is called the magnitude of completeness, M _c, and is equal to the location threshold (Fig. 4). Assessment of M _c from the FM curve provides a check on the validity of the network design and a boundary condition for statistical analysis. Until 2014 M _c=1.5 was adopted for the north of the Netherlands.

Paleja & Bierman (2016) report on an initial analysis of temporal variations in M _c across the entire Groningen field. The maximum curvature method (Wiemer & Wyss, 2000) was used and four time periods selected. They conclude that, for the period 2003–2009, M _c=1.2 and, for the period September 2014–September 2016, M _c =0.6. In between these periods M _c=0.9. Since the FM curves are not smooth, the authors did not calculate synthetic data, but determined the maximum number of recorded events. We did implement the full method and calculated the goodness-of-fit parameter R (Table 2) and came to comparable conclusions, although the second interval showed a higher value of M _c. In the last period, M _c dropped to lower values, but the procedure is not very stable. Paleja & Bierman (2016) showed uncertainty by application of bootstrapping and concluded that the results over the last period are different from the earlier epochs.

Table 2. M_c for four non-overlapping periods.

Earthquake mechanism

Kraaijpoel & Dost (2013) published source parameters for four events derived from P-wave first motion directions and amplitude ratios The mechanisms show normal faulting at steep angles (60–70°) and a strike of around 290°. In addition, there was some indication of a small rake component (−105 to −120°). However, azimuthal distribution of the boreholes with respect to earthquakes in the Groningen field was limited, which also applied to the number of accelerometers deployed at the time (2005–2010). With the installation of the new network these limitations are resolved and it is expected that the accuracy of the source mechanism solutions will be improved. The source mechanisms inferred from the data are in line with a reactivation of existing faults in the reservoir.

Probabilistic seismic hazard assessment

First simple PSHA calculations for Groningen were presented by Van Eck et al. (2006). For the north of the Netherlands, Dost et al. (2004) derived a GMM based on measured acceleration, mainly from the Roswinkel gas field. After surface acceleration data became available in Groningen, it became clear that the model provided an overestimation and a new application-specific GMM should be developed (Bommer et al., 2016). An overview of the development of GMMs for Groningen is given in Bommer et al. (2017). Other important parameters for the PSHA include the seismic activity rate, b-value and the maximum magnitude.

A first update of the PSHA for Groningen, presented in Dost et al. (2013), was based on the v0 version of the GMM for the region and a new zonation estimated from recorded seismicity, known faults and information on compaction. An existing GMM for the European–Mediterranean region (Akkar et al., 2014) was adapted and adjusted for magnitudes below M=4 in order to fit the model to Groningen conditions. Results for a 475-year return period showed high values for maximum peak ground acceleration (PGA) (0.42g; 1g corresponds to 9.81ms⁻²) and peak ground velocity (PGV) (16cms⁻¹), mainly due to large uncertainties in the model. Bourne et al. (2015) present results for the NAM model for the same region. For a 10-year period (2013–23) they found a 2% chance of exceedance, a maximum PGA of 0.57g and a maximum PGV of 22cms⁻¹. The spatial patterns of both models are comparable.

Next versions of the PSHA for Groningen were published after the development of v1 (Dost & Spetzler, 2015) and v2 (Spetzler & Dost, 2016). The v1 model was fully based on accelerometer recordings from only induced earthquakes in Groningen. Site amplification effects were accounted for through a simple field-wide factor that remained linear under all loading levels. Implementation of v1 led to a reduction of the maximum PGA in the hazard map to 0.36g. Since the introduction of a laterally varying site effect in v2, the PSHA method is combined with the convolution approach described in Bazzuro & Cornell (2004). The NAM model also computes the surface hazard including the amplification factor at the specific site but by means of a Monte Carlo approach.

Spetzler & Dost (2016) compared the two PSHA approaches, not only for PGA, equal to spectral acceleration at T=0.01s, but also for spectral acceleration at other periods at two locations (Groningen and Loppersum). For the implementation of the v2 model, the NAM model uses a distribution of M _max values (M 5, 5.75 and 6.5) with equal weights, while the KNMI model uses M _max=5. Deaggregation shows that for a return period of 475 years the main contribution to the hazard comes from events of M 4–5 (Bourne et al., 2015).

Despite the differences in methodology, the hazard maps for the V2 GMM calculated by the KNMI and NAM are remarkably similar. A comparison of the two hazard maps is shown in Figure 7. The two maps show the same patterns inherent to the zone-specific amplification factors and the same maximum PGA,0.22g.

Fig. 7. Probabilistic seismic hazard map for Groningen for the period T=0.01s. The return period is 475 years according to Eurocode 8. The black solid line indicates the contours of the Groningen gas field. Left: KNMI hazard map, max PGA = 0.22g (Spetzler & Dost, 2016). Right: NAM hazard map for 33bcma⁻¹ production, max PGA=0.21g (NAM, 2016b).

Comparison between KNMI and NAM results for the spectral accelerations is shown in Figure 8. Both models provide similar results for T<1s. At longer periods the two models start to differ, and the cause of this difference is being investigated. A possible reason for this difference is the influence of the choice of M _max.

Fig. 8. Spectral acceleration for locations in Groningen city (left) and Loppersum (right). KNMI results in red, NAM results in green (Spetzler & Dost, 2016).

The next update is expected mid-2017 and will show results of the implementation of GMM-v4 (Bommer et al., 2017). Apart from the spectral acceleration at 23 periods, this model also provides equations for PGV, which was lacking in previous versions (v1–3). Due to the non-stationarity of the seismic activity in Groningen, zonation and hazard parameters must routinely be updated.

Conclusions

Long-term monitoring of induced seismicity in the north of the Netherlands and the recent expansion of the network in Groningen forms the basis for understanding the processes responsible for the generation of earthquakes. An accurate estimate of the seismic hazard is essential for the risk assessment and the subsequent hazard mitigation planning. Since induced seismicity is a non-stationary process, understanding temporal and spatial variations in seismicity patterns and their relation to changes in production is essential. However, due to the relatively small number of events recorded, it takes time to detect statistically significant changes. The new borehole network in Groningen decreases the magnitude of completion (M _c), strongly increasing the database of events that could be used. Seismological studies are focused on the further reduction of uncertainty in the hazard assessment. The use of the borehole strings enables the calculation of interval shear velocities in the upper 200m (Hofman et al., in press) and the same for the damping (Q), both essential parameters in the calculation of the site amplification factors as part of the next version of the GMM.