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, Reference Dost, Haak, Wong, Batjes and de Jager2007). 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., Reference Dost, Goutbeek, Van Eck and Kraaijpoel2012). 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 Reference Hofman, Ruigrok, Dost and Paulssenpress).

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., Reference Bourne, Oates, Bommer, Dost, van Elk and Doornhof2015). 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., Reference Lienert, Berg and Frazer1986), the EDT method by Lomax (Reference Lomax2006) was used (Spetzler & Dost, Reference Spetzler and Dost2017) 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 (Reference Pickering2015) 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 (Reference Kraaijpoel and Dost2013) 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, Reference Spetzler and Dost2017), 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, Reference Waldhauser and Ellsworth2000). The relocation was tested for a subset of the seismicity (e.g. Jagt et al., Reference Jagt, Ruigrok and Paulssen2017). 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., Reference Bourne, Oates, van Elk and Doornhof2014). 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. (Reference Van Eck, Goutbeek, Haak and Dost2006) 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 (Reference Dost and Kraaijpoel2013) showed a change around 2003, also reflected in the FM relations for the period before and after the change. Nepveu et al. (Reference Nepveu, van Thienen-Visser and Sijacic2016) 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 (Reference Tinti and Mulargia1987). Marzocchi & Sandri (Reference Marzocchi and Sandri2003) 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., Reference Volant, Grasso, Chatelain and Frogneux1992).

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. (Reference Paleja, Bierman and Jones2016) concluded that activity rate in this area was reduced with respect to the year preceding the change in production. Nepveu et al. (Reference Nepveu, van Thienen-Visser and Sijacic2016) 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., Reference Van Gent, Back, Urai, Kukla and Reicherter2009) 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. (Reference Bourne, Oates, van Elk and Doornhof2014) 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., Reference Dost, van Eck and Haak2004). 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., Reference Dost, Edwards and Bommer2016). At this moment a reassessment of these results is ongoing.

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 (Reference Paleja and Bierman2016) report on an initial analysis of temporal variations in M _c across the entire Groningen field. The maximum curvature method (Wiemer & Wyss, Reference Wiemer and Wyss2000) 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 (Reference Paleja and Bierman2016) 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 (Reference Kraaijpoel and Dost2013) 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. (Reference Van Eck, Goutbeek, Haak and Dost2006). For the north of the Netherlands, Dost et al. (Reference Dost, van Eck and Haak2004) 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., Reference Bommer, Dost, Edwards, Stafford, van Elk, Doornhof and Ntinalexis2016). An overview of the development of GMMs for Groningen is given in Bommer et al. (Reference Bommer, Dost, Edwards, Kruiver, Ntinalexis, Rodriguez-Marek, Stafford and Van Elk2017). 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. (Reference Dost, Caccavale, van Eck and Kraaijpoel2013), 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., Reference Akkar, Sandikkaya and Bommer2014) 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. (Reference Bourne, Oates, Bommer, Dost, van Elk and Doornhof2015) 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, Reference Dost and Spetzler2015) and v2 (Spetzler & Dost, Reference Spetzler and Dost2016). 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 (Reference Bazzurro and Cornell2004). 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 (Reference Spetzler and Dost2016) 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., Reference Bourne, Oates, Bommer, Dost, van Elk and Doornhof2015).

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, Reference Spetzler and Dost2016). 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, Reference Spetzler and Dost2016).

The next update is expected mid-2017 and will show results of the implementation of GMM-v4 (Bommer et al., Reference Bommer, Dost, Edwards, Kruiver, Ntinalexis, Rodriguez-Marek, Stafford and Van Elk2017). 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.