Construction of earthquake location uncertainty maps for
the Netherlands
KNMI, R&D Seismology and Acoustics
August 24, 2023

KNMI

Samenvatting
Een groot deel van Nederland is bedekt met een dikke laag geconsolideerde sedimenten. In
het onderliggend vastere deel van de bovenkorst vinden verschillende antropogene processen
plaats. Olie en gas worden geproduceerd, zout wordt gewonnen met oplossingsmijnbouw,
vloeistoffen en gassen worden opgeslagen, en geothermische energie wordt gewonnen. Al deze
activiteiten kunnen het spanningsveld rondom bestaande breuken veranderen en daarmee
mogelijk leiden tot seismiciteit.
Gedurende de afgelopen decennia is een uitgebreid sensoren netwerk opgebouwd om
seismiciteit te monitoren. Dit netwerk heeft een hoge sensordichtheid op plekken waar veel
gasvelden zijn geproduceerd. Op andere plekken zijn er weinig sensoren. In dit rapport
is de capaciteit van het netwerk onderzocht met de configuratie van september 2021. Dit
is gedaan om te achterhalen in hoeverre het netwerk voldoet om huidige en toekomstige
ondergrondse activiteiten seismisch te monitoren. Ook worden de resultaten getoond voor
de september 2022 staat van het netwerk. Voor dit onderzoek zijn nieuwe gereedschappen
ontwikkeld om vast te kunnen stellen wat de locatie limieten zijn en de verwachte locatie
onzekerheden.
De belangrijkste parameter om de capaciteit van een seismisch monitoringsnetwerk te
kwantificeren is magnitude van completie (MoC). Deze parameter beschrijft de spatiële variatie van de minimale magnitude waarvan bijna elke beving gelokaliseerd kan worden. De
locatie van een beving kan vastgesteld worden als de beving op tenminste 3 sensoren is
gedetecteerd. Wanneer een beving gelokaliseerd kan worden, is het van belang om te weten
met welke nauwkeurigheid dit gedaan kan worden. Kennis van de diepte-onzekerheid is belangrijk, bijvoorbeeld, om tektonische bevingen te kunnen onderscheiden van geı̈nduceerde
bevingen. Kennis van de epicentrale onzekerheid (onzekerheid in het horizontale vlak) is
belangrijk om een beving te kunnen relateren aan een bepaalde breuk en mijnbouwoperatie. Met kennis van de verwachte locatie-onzekerheid van het huidige netwerk kan bepaald
worden of extra sensoren nodig zijn.
In dit rapport zijn kaarten ontwikkeld die de verwachte locatie-onzekerheid tonen voor
bevingen in de bovenkorst in Nederland. Een belangrijk onderliggend model is de onzekerheidsverdeling van seismische reistijden. Dit model is afgeleid van vele metingen van
bevingen in Nederland. Met een Bayesiaanse analyse zijn onzekerheden in reistijden gepropageerd naar onzekerheden in locatie. Er is vanuit gegaan dat locatie gedaan wordt
met verschillen in aankomsttijden van de primaire (P) golf bij verschillende sensoren, en
met reistijdverschillen tussen de P en schuif (S) golf aankomsten. Deze data attributen
bevatten complementaire informatie. Een gezamenlijke inversie van beiden geeft de kleinste
locatie-onzekerheden. Met deze gezamenlijke inversie zijn de kaarten berekend voor lokale
magnitudes van 0,5 tot 4,0 met stappen van 0,5.
De kaarten zijn berekend over een grid van 1x1 km. Voor elke potentiële beving is de
locatie kansendichtheidsfunctie bepaald, als functie van diepte en in het epicentrale vlak.
Vervolgens is aangenomen dat deze kansendistributies goed kunnen worden benaderd als
zijnde (multivariate) normaalverdelingen, zodat ze beschreven kunnen worden met standaard
deviaties. Deze standaard deviaties worden getoond op de kaarten. Daarmee kan dan weer
de grootte van het gebied, of dieptebereik, berekend worden dat bijvoorbeeld 68% of 95% van
de kansen bevat. De kaarten zijn doorgerekend voor een beving diepte van 3 km. Daarnaast
is een gevoeligheidsanalyse gedaan om te evalueren hoe de locatie onzekerheid verandert als
functie van verschillende parameters. Bij een grotere diepte dan 3 km neemt de epicentrale
onzekerheid toe en zal ook de diepte-onzekerheid iets toenemen.
Signaal-ruisverhoudingen spelen een belangrijke rol in de bepaling welke sensoren een
mogelijke beving zouden detecteren. Voor het bepalen van het verwachte signaalniveau
is een P-golf grondbewegings-vergelijking afgeleid. Deze vergelijking beschrijft de P-golf
amplitudes als functie van brondiepte, epicentrale afstand en magnitude. Het ruisniveau per
ontvanger is verkregen met de lokale meting van de seismische achtergrondsruis in dezelfde

1

KNMI

frequentieband die gebruikt wordt voor detectie (5 tot 40 Hz). Als ruisniveau is de 90ste
percentiel van de ontvanger-specifieke ruismeting gebruikt. Dit is het ruisniveau dat 10%
van de tijd wordt overschreden.
De MoC en locatie-onzekerheids kaarten tonen een grote variatie van waarden over Nederland. De belangrijkste onderliggende oorzaak is de sterke variatie in sensordichtheid. Een
andere belangrijke reden is de sensorplaatsing, te weten nabij het aardoppervlakte of in
200 m diepe boorputten. Gemiddeld genomen zijn de ruisniveaus op 200 m diepte een factor
26 keer lager dan aan het aardoppervlak, in de detectie frequentieband. Anderzijds zijn
ook de signaalamplitudes een factor 4 lager op 200 m diepte. Beide effecten samen leveren
een winst in signaalruisverhouding van een factor 6.5 met plaatsing op diepte. Daarnaast
is de lokale setting van belang. Een regionaal netwerk van sensoren in een druk stedelijk
deel van Zuid-Holland presteert minder goed dan een vergelijkbare ontvangerverdeling in
een rustig ruraal deel van Twente. De sensoren die het beste presteren wat betreft detectie
van seismiciteit zijn sensoren geplaatst op hard gesteente in Zuid-Limburg.
Het rapport bevat magnitude van completie (MoC) kaarten voor 2021 en 2022. De
onderliggende waardes zijn te raadplegen via KNMI (2023a). De komende jaren zullen er
updates blijven verschijnen op het KNMI Data Platform. De kaarten zijn gemaakt met
gemeten stations-specifieke ruiscondities en zijn daarmee een verbeterslag ten opzichte van
eerdere kaarten die gemaakt zijn met een aanname van gemiddelde ruiscondities over het
seismische netwerk. Uit de nieuwe 2021 en 2022 kaarten volgt dat de MoC lager is of gelijk
is aan 2,0 voor heel Nederland. In gebieden met hoge sensordichtheden (Groningen, delen
van Friesland, Drenthe, Twente, Zuid-Limburg en de regio rondom Alkmaar) is de MoC
lager dan 0,5. Binnen het relatief dichte netwerk in Zuid-Holland wordt 0,5 niet bereikt
vanwege hoge ruiscondities. In het midden van het land zijn er niet veel seismische sensoren
aangesloten op het KNMI netwerk en bereikt de MoC waardes tussen de 1,5 en 2,0. Voor de
berekening van de MoC is de 90ste percentiel van de ruiscondities gebruikt. Dientengevolge
zou in ongeveer 10% van de gevallen de MoC (net) niet gehaald worden.
Voor hogere magnitudes is de locatie onzekerheid lager, omdat dan een groter netwerk
van sensoren gebruikt kan worden om de locatie van de beving vast te stellen. Binnen
de dichtbemonsterde regionale netwerken in Nederland is de epicentrale locatie-onzekerheid
veelal kleiner dan een paar honderd meter. Deze onzekerheid groeit tot een paar kilometer,
aan de randen van gebieden waar locatie nog net mogelijk is.
De onzekerheid in dieptebepaling gedraagt zich anders dan de epicentrale onzekerheid.
Over het algemeen kan de diepte met minder zekerheid bepaald worden dan het epicentrum.
Voor het vaststellen van de diepte is het van belang om sensoren te hebben nabij het epicentrum, zodat het hyperbolische deel van de reistijdencurve bemonsterd wordt. Voor het
vaststellen van het epicentrum is het belangrijker om een goede azimutale distributie van
ontvangers te hebben. Daarnaast is voor het vaststellen van diepte van belang om informatie
te hebben over de lokale snelheidsopbouw van de ondergrond. Voor het vaststellen van het
epicentrum kan veelal een gemiddeld of een aardbevings-specifiek reistijdenmodel gebruikt
worden.
De berekende locatie onzekerheid wordt minder correct weergegeven voor bron-ontvanger
configuraties waarbij het azimutale gat groter is dan 250o . Het azimutale gat is de grootste
hoek tussen twee ontvangers, zoals gezien vanuit het epicentrum. Wanneer het azimutale
gat groter is dan 250o , is de locatie kansenverdeling niet meer goed te beschrijven als een
(multivariate) normaalverdeling. De grootste onzekerheden treden op wanneer ontvangers
bijna op een lijn staan, hetgeen overeenkomt met een azimutaal gat van bijna 360o . In dat
geval kan de locatie ambigu zijn en is het belangrijk om de polarisatie van P-golven mee
te nemen om deze ambiguı̈teit te verhelpen. In de appendix zijn kaarten opgenomen met
daarop de grootte van het azimutale gat.
In de appendix staan de verkregen locatie-onzekerheidskaarten voor het Nederlands seismisch netwerk van 2021 en 2022. Ook is de onderliggende database publiek beschikbaar met
daarin de waardes voor elke gridpunt (KNMI , 2023b). Wanneer precieze waarden op een

2

KNMI

bepaalde plek van belang zijn, kan deze database geraadpleegd worden. De komende jaren
zullen er updates blijven verschijnen op het KNMI Data Platform.
De ontwikkelde hulpmiddelen zijn ook in te zetten om de capaciteiten van scenario
netwerken door te rekenen. Voor een mogelijk toekomstige sensorendistributie kan uitgerekend worden wat de verwachte impact is op MoC en locatie-onzekerheid.

3

KNMI

Summary
A large part of the Netherlands is covered by a thick blanket of unconsolidated sediments.
In the underlying more consolidated part of the upper crust, various anthropogenic activities
are taking place. Oil and gas are produced, salt is extracted through solution mining, liquids
and gasses are stored and geothermal heat is exploited. All these activities can alter the
stress field around existing faults and therewith potentially induce seismicity.
Over the years, an extensive seismic network has been installed to monitor seismic activity. This network has a dense distribution of sensors in regions with past subsurface activity
and is more sparse in other regions. In this report, the network capabilities are assessed
for the September 2021 network configuration, to find out whether this network suffices to
monitor current and future subsurface operations. Also, results are shown for the September 2022 network configuration. New tools are developed to assess location thresholds and
expected location uncertainties.
The main seismic network-assessment parameter is magnitude of completeness. This
parameter describes the spatial variation of the minimum magnitude for which almost every
earthquake can be located. The location of an earthquake can be determined if the signal is
detected on at least three sensors. As soon as an event can be located, it is important to know
with which precision. Knowing the depth uncertainty is important, e.g., to discriminate
induced events from tectonic events. Knowing the epicentral uncertainty is important, e.g.,
to assign the event to one of multiple nearby mining operations. By knowing the expected
location uncertainties pertaining to the current network, decisions can be made on the need
of additional sensors.
In this report, maps are developed that show the expected location uncertainty over
the Netherlands for upper-crustal seismicity. An important prerequisite is a travel-time
uncertainty model, which is derived from a large collection of records from induced events.
Uncertainty in travel time is mapped to uncertainty in location using a Bayesian framework.
Location is considered using differential P-wave arrival times over the network (P-delays),
and also by using the travel-time difference between first P and S waves at each receiver
(P-S delays). It is shown that both data attributes yield complimentary information on the
location problem and that a joint inversion with both P- and P-S delays yields smallest location uncertainty. For this specific combination of data attributes, the location-uncertainty
maps are computed for local magnitudes of 0.5 until 4.0, with steps of 0.5.
The maps have been computed over a grid of 1x1 km. For each scenario event, the location
probability density function has been determined, both as function of depth as as function of
the horizontal coordinates. The assumption is made that these probability density functions
can be approximated well as being (multivariate) normal distributions, such that they can
be described with standard deviations. These standard deviations are shown on the maps.
From the standard deviations the size of, for example, the 68% and 95% confidence zones
can be computed. The maps have been computed for a source depth of 3 km. Besides, a
sensitivity analysis has been performed to evaluate how the location uncertainty varies as
function of different parameters. When the source is placed at depths larger than 3 km, the
epicentral uncertainty increases, and also the depth uncertainty somewhat increases.
Signal-to-noise levels play an important role in the determination at which sensors a
potential earthquake is likely detected. The signal is modeled through deriving P-wave
ground-motion prediction equations that describe the expected P-wave amplitudes as function of magnitude and distance from the source. The noise level is based on receiver-specific
recordings of background seismic noise in the frequency band that is used for detecting induced events (5 to 40 Hz). As the location-specific noise level, the 90th percentile of the
recorded noise distributions is used, meaning a noise level that is exceeded 10% of the time.
The magnitude-of-completeness and location-uncertainty maps show large variations over
the Netherlands. The main underlying factor is the large variation in sensor density. Another
important factor is the position of the sensor, i.e., on the Earth’s surface or 200 m below the

4

KNMI

surface in a borehole. On average, in the detection frequency band, noise levels are a factor
of 26 lower at 200 m depth than at the Earth’s surface. On the other hand, signal levels are a
factor of 4 lower at 200 m depth. Combining the two effects, this results in a signal-to-noise
gain of 6.5 at the 200 m depth sensors. Moreover, the local setting is important. A regional
network in South-Holland performs worse than a similar network in Twente because of much
higher seismic noise conditions near Rotterdam and the Hague than in the north of Twente.
The sensors with the best detection capabilities are located on hard-rock conditions in the
southernmost part of the Netherlands.
The report contains magnitude-of-completeness (MoC) maps for 2021 and 2022. The underlying values can be obtained at KNMI (2023a). The coming years, updates will continue
to appear at the KNMI Data Platform. Previous maps were made with assumed average
noise conditions over the seismic network. The new maps have been made with recorded
station-specific noise conditions, which yields a more accurate representation of the MoC
that can regionally be achieved. On the new 2021 and 2022 maps, MoC is smaller or equal to
2.0 in the Netherlands. In areas with dense sensor distributions (Groningen, parts of Friesland, Drenthe, Twente, south of Limburg and region around Alkmaar) the MoC is lower
than 0.5. The network in South-Holland does not reach MoC=0.5 due to high local noise
conditions. In the middle of the country, there are not many seismic sensors and the MoC
has values between 1.5 and 2.0. For the computation of the MoC, the 90th percentile of the
noise distribution has been used. As a consequence, in about 10% of the cases the stated
MoC can (just) not be reached.
Location uncertainty reduces with event magnitude as a better spatial distribution of
receivers detects the event for higher magnitudes. Within the dense regional networks, the
epicentral uncertainty is within a few hundred meters. The uncertainty grows to a few
kilometers at the edges of the regions where location is just possible.
Depth uncertainty behaves quite differently from epicentral uncertainty. In general, the
depth is less well constrained than the epicenter. For constraining the depth, it is important
to have receivers nearby the epicenter, so that the hyperbolic part of the travel-time curve is
sampled well. For constraining the epicenter it is more important to have a good azimuthal
distribution of receivers. For accurate depth estimation, local velocity information needs to
be available. For estimating epicenter, one can mostly rely on background or event-derived
travel-time models.
The computed location uncertainty is less well represented for source-receiver configurations for which the azimuthal gap is larger than 250o . The azimuthal gap is the largest
angle between two receivers as seen from the epicenter. When the azimuthal gap is larger
than 250o , the location probability density function cannot be approximated well as being
a multivariate normal distribution. The largest uncertainties occur when the receivers are
nearly on one line, which corresponds to an azimuthal gap of nearly 360o . In this case, location can be ambiguous and the polarisation of P-waves is needed as an additional constraint
to remove this ambiguity.
In the appendix, the location-uncertainty maps can be found for the Netherlands seismic
network state of 2021 and 2022. Also the underlying database, with the values at each grid
point, is publicly available (KNMI , 2023b). This database is to be consulted if the precise
values at a certain location are of interest. The coming years, updates will continue to
appear at the KNMI Data Platform.
The developed tools also provide a means to study design options for future extensions
of the network. A scenario layout can be checked for its capability in terms of the minimum
magnitude of locatable earthquakes and their associated uncertainty in location.

5

Contents
1 Introduction

8

2 Prerequisites
2.1 List of stations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 P-wave ground-motion prediction equations . . . . . . . . . . . . . . . . . . .
2.2.1 Database . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 Ground motion model . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.3 Model calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.4 Quality checking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Station-specific noise distributions . . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Power-spectrum density . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.2 Root-mean-square velocity . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.3 Noise reduction with depth . . . . . . . . . . . . . . . . . . . . . . . .
2.4 Picking threshold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5 Hard-rock correction factors . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6 Travel-time models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6.1 Background polynomial model . . . . . . . . . . . . . . . . . . . . . .
2.6.2 Event-specific polynomial model . . . . . . . . . . . . . . . . . . . . .
2.6.3 Hyperbolic travel times . . . . . . . . . . . . . . . . . . . . . . . . . .

11
12
15
15
16
17
20
25
25
27
31
35
39
41
41
44
45

3 Location-uncertainty method
3.1 Bayesian framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.1 P delay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.2 P-S delay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.3 Joint inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Uncertainty quantification . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Comparison of data attributes . . . . . . . . . . . . . . . . . . . . . . . . . . .

47
47
47
49
50
50
53

4 Sensitivity analysis
4.1 Travel-time uncertainty . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Source depth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Source-receiver distance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4 Number of receivers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5 Azimuthal gap . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

56
56
59
60
63
64

5 Location-uncertainty maps
5.1 Settings and quality checks . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.1 Travel-time uncertainty . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.2 Grid-search settings . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.3 Discontinuity in uncertainty . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Epicentral uncertainty . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3 Depth uncertainty . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

67
68
68
69
72
73
77

6

CONTENTS

5.4

KNMI

Magnitude of completeness . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4.1 Sensitivity tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4.2 2021 Map . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4.3 2022 Map . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

78
79
79
83

6 Discussion and Conclusions
6.1 Detection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2 Location PDF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.3 Knowledge of travel-time uncertainty . . . . . . . . . . . . . . . . . . . . . . .
6.4 Magnitude of completeness . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.5 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.6 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

85
85
86
89
90
90
91

A 2021 Station list

96

B 2021 NL-wide maps with location-uncertainty parameters
101
B.1 Largest-horizontal-location-uncertainty (σ1 ) maps . . . . . . . . . . . . . . . . 101
B.2 Vertical-location-uncertainty (σZ ) maps . . . . . . . . . . . . . . . . . . . . . 111
B.3 Azimuthal-gap (ψ) maps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
C 2022 Station list

127

D 2022 NL-wide maps with location-uncertainty parameters
132
D.1 Largest-horizontal-location-uncertainty (σ1 ) maps . . . . . . . . . . . . . . . . 132
D.2 Vertical-location-uncertainty (σZ ) maps . . . . . . . . . . . . . . . . . . . . . 142
D.3 Azimuthal-gap (ψ) maps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150

7

1 Introduction
In the Netherlands there is a growing number of mining operations. Comprehensive knowledge has been built up of the upper crust. E.g., detailed velocity models (Van Dalfsen
et al., 2006) and fault maps (Ten Veen et al., 2022) are openly available. This allows precise
planning of subsurface operations, ranging from hydrocarbon production to salt extraction.
These subsurface operations perturb the existing stress field and therewith have a potential
to induce seismic activity on pre-existing faults (e.g., Buijze et al., 2023). This seismicity is
to be monitored, especially if it reaches levels that can be felt or inflict damage. Monitoring
is the first step to find, and eventually mitigate, the underlying cause. This cause could be
natural —there is tectonic stress build up and release in the Netherlands– or could have an
anthropogenic origin. It is important that the location accuracy is sufficient to pinpoint the
natural process or mining activity responsible for the seismic event. The Royal Netherlands
Meteorological institute (further: KNMI) was asked by the The Dutch State Supervision of
Mines (further: SodM) to assess the location accuracy of the current seismic network.
Some gasfields are amenable to seismicity (e.g. Roswinkel and Groningen) while others
have not shown any seismic event over their lifetime (e.g. gasfields in Twente). The occurrence or absence of seismicity over many different gasfields has been used to identify reservoir
parameters that correlate with seismicity (Van Eijs et al., 2006; van Thienen-Visser et al.,
2012). Despite the knowledge accrued, however, forecasting the occurrence of seismicity
remains a challenge. For that reason, mining operations are to be seismically monitored to
some degree, even if the chance of seismicity is thought to be small.
The monitoring capabilities largely depend on the seismic network and noise conditions.
The Netherlands has an extensive network of borehole geophones, surface accelerometers
and broadband seismometers (Dost et al., 2017). Prior to 1990, most sensors were located
in the south of the Netherlands, where tectonic activity occurs related to the Lower Rhine
Graben (e.g., Hinzen et al., 2021). During the end of the eighties, the first induced seismic
events were detected in the north of the Netherlands, which were attributed to gas extraction
(Haak , 1993). Since then, there has been a continuous expansion of the network to monitor
induced seismicity.
Most sensors that are deployed in the Netherlands Seismic and Acoustic Network (KNMI ,
1993) have instrument noise that is lower than the seismic background noise, at least in the
frequency band that is relevant for detecting induced seismicity. As such, a bottleneck in
detection of events is the seismic noise, which varies largely from site to site. For example,
seismic noise levels are much higher near a highway than in the middle of a forest. To
properly assess the network capabilities, this location-specific background noise needs to be
taken into account.
For the Netherlands, a relatively new potential source of induced seismicity is related to
geothermal energy production (Buijze et al., 2019). Up to date, a magnitude 0.0 event has
been detected near an operation in South-Holland (Naranjo et al., 2022) and a magnitude

8

KNMI

N

Z

Figure 1.1: In this report, for each point in the Netherlands (left-hand figure) the epicentral
probability-density function (PDF, middle figure) and depth PDF is derived (right-hand
figure). The epicenral PDF is described with three parameters, σ1 , σ2 and θ. The depth
PDF is parameterized with σZ . The coloured areas of the epicentral and depth PDF denote
the 95% confidence areas.
1.7 event near an operation in Limburg (Vörös and Baisch, 2022). Exploration and exploitation also takes place in areas of the Netherlands where the seismic network is sparse
(Mijnlieff , 2020), because of an absence of gas production or tectonic activity. SodM is involved in the permitting process as advisor to the Ministry of Economic Affairs and Climate
Policy. One of the questions that needs to be answered in light of safe operations is whether
the nationwide existing network suffices in terms of monitoring, or that new stations are
needed to adequately monitor a geothermal operation. SodM has asked the seismic network operator KNMI to assess the current network capabilities. To monitor seismicity and
to connect detected seismicity to an underlying cause, two network capabilities are relevant:
Epicentral uncertainty. This is a description of location uncertainty in the horizontal plane (map view). With this uncertainty it can be assessed whether events from nearby
operations can be well separated. Or, it can be assessed whether the uncertainty is low
enough to pinpoint the specific fault on which the rupture must have taken place.
Depth uncertainty. This parameter describes the location uncertainty in the vertical
direction. Knowing the likely depth uncertainty is important to find out whether the mining operation (in the upper crust) caused the event, or whether a tectonic release of stresses
in the deeper crust is more likely. Another option is that the mining activity on a shallow
depth triggered an event on a somewhat larger depth. Moreover, in case multiple mining
operations are stacked, it can be assessed whether the depth-location accuracy would be
sufficient to distinguish events from these different operations.
The location accuracy that can be achieved is a function of event magnitude. For higher
magnitudes, more stations detect the event above the noise floor, which enables a better constrained source solution. Fig. 1.1 illustrates the type of maps that are developed within this
report. For each point in the Netherlands and for a range of magnitudes (M=0.5, 1.0, 1.5,
2.0, 2.5, 3.0, 3.5 and 4.0), the location uncertainty is modelled and mapped. This amounts
to generating maps for the Netherlands showing the epicentral uncertainty (which is parameterized with σ1 , σ2 and θ) and maps showing depth uncertainty (which is parameterized
9

KNMI

Input:
- Station list
- P-wave GMPE
- Noise per station
- Traveltime uncertainty
- Picking threshold
- Hardrock correction factor

Engine:

Output:

For each ML
For each source gridpoint
- Determine pick list
- Compute epicentral uncertainty
- Compute depth uncertainty

- Maps depicting epicentral
location uncertainty
- Maps depicting depth
location uncertainty

Figure 1.2: A flowchart summarizing the different components that are developed within
this project. The required ’input’ or prerequisites are listed on the left. The ’Engine’ is
an implementation of a Baysian method to compute location probability density functions.
The input and engine are developed to yield the desired ’Output’: epicentral and depth
uncertainty maps for the Netherlands.
with σz ). The exact meaning of these uncertainty parameters is discussed in Chapter 3.
Fig. 1.2 shows a flow-chart for producing the location uncertainty maps. The outline of
the report follows this flow-chart. In Chapter 2 (Prerequisites) all the different inputs for
the analyses are derived and discussed. E.g., a P-wave ground-motion prediction equation
is derived in order to model the amplitude at a station at a given distance from a scenario
event. For that station, the noise conditions need to be known to determine whether it is
likely that this station detects the scenario event. In Chapter 3 the methods are developed
to propagate errors in timing of P- and S-waves (due to unknowns in the velocity structure, picking errors or modelling errors) to location uncertainty. The implementation of this
’method’ forms the ’engine’ for computation of the location-uncertainty maps. In Chapter
4 it is assessed what to be the influence on location uncertainty of different parameters, like
depth of the source, azimuthal gap and number of receivers. In Chapter 5 the produced
maps are shown and discussed, followed by a final Discussion and Conclusions chapter and
an appendix that contains the NL-wide location-uncertainty maps.
As a by-product, in this report the magnitude of completeness (MoC) map is updated.
In order to find out that an earthquake has happened, it first needs to be detected. Detection can be done already with a single sensor. However, a single detection is not sufficient
to constrain both depth and epicenter of the event, or to be sure that the recorded ground
motion was caused by a fault rupture in the upper crust. In order to determine a location
with some confidence, detection at a minimum of three receivers is required. The MoC
at a certain point is the minimum magnitude for which these three detections are likely
reached. A MoC map shows how this minimum magnitude varies over an area. For the
Netherlands, a first MoC map was derived in Van Eck et al. (2004) an later updated in Dost
et al. (2012), Dost et al. (2017) and Kruiver et al. (2021) to incorporate network extensions.
These maps were based on average noise conditions at 200 m deep geophones and surface
receivers that existed before 2004. Besides modeling the MoC spatially, one can also, on
hindsight, evaluate which MoC was achieved over a certain region. Such estimates, using
the KNMI induced event catalogue, can be found in van Thienen-Visser et al. (2016), Dost
et al. (2012) and Muntendam-Bos and Grobbe (2022). In this report, a thorough update
is made of the anticipated MoC by incorporating recent developments of the network, the
newly derived P-wave ground-motion prediction equations and site-specific noise curves.

10

2 Prerequisites
The construction of a location-uncertainty map requires many elements. These prerequisites
are assembled and developed in this chapter:
• List of stations that are included in the analysis. This must be a list of stations that
can actually be processed within the operational chain at KNMI. Then still, different
choices are possible. Are Belgian and German stations included to cover the border
regions? Do we include temporary networks? Do we only include stations that yield
real-time data? Etc.
• P-wave ground-motion prediction equations or GMPE are relations that model
the expected maximum amplitude for incoming primary (P) waves, as function of distance and magnitude. Detection of induced seismicity is done based on P-waves using
relatively high frequencies (5 to 40 Hz). P-wave ground-motion prediction equations
(GMPEs) are derived, specifically for this frequency band and (induced) seismicity in
the upper crust.
• Station-specific noise distributions are needed to assess whether incoming earthquakerelated waves would likely stand out from the noise. For the detection frequency band
and using at least half a year of data, the seismic noise distribution is obtained for
each station from the list.
• Picking threshold is the average minimum factor needed between signal and noise
for an automatic detection to be made. This factor is empirically derived using a
collection of recorded seismograms.
• Hard-rock correction factors are factors used to adapt the P-wave GMPEs to
hard-rock sites. They are needed because the P-wave GMPEs are calibrated purely
with recordings on soft sites (unconsolidated sediments). A few stations in the Netherlands record on hard rock, which results in lower recorded amplitudes for the same
earthquake and distance. Based on assumed average sediment and hard-rock elastic
properties, it is derived how much lower these amplitudes would be.
• Travel-time models describe the arrival times of P- and S-waves as function of
distance, and the variability therein for induced seismicity in the Netherlands. This
variability, leading to uncertainty of the actual traveltimes, is derived from a large
database of recordings. In the following chapters, this uncertainty in time is mapped
uncertainty in location.

11

2.1. LIST OF STATIONS

KNMI

Figure 2.1: Stations in the Netherlands and direct surroundings, which data can be obtained
through https://www.orfeus-eu.org/fdsnws/ or http://eida.bgr.de/fdsnws/. The overview is
from September 2021. Green triangles denote broadband stations, orange triangles are shortperiod stations and purple squares are accelerometers. Below many of the accelerometers,
there is a co-located string of geophones for which the deepest node is at 200 m depth.

2.1

List of stations

In the Netherlands and surroundings, different seismic networks are active. We use FDSN
web services (https://www.fdsn.org/webservices/) to collect all publicly available stations
for the following networks
• BE: Belgian Seismic Network
• BQ: Bensberg Earthquake Network
• GE: GEOFON Seismic Network
• GR: German Regional Seismic Network
• NL: Netherlands Seismic and Acoustic Network
• NR: NARS - The Network of Autonomously Recording Seismographs
• RN: RuhrNet – Seismic Network of the Ruhr-University Bochum
The locations of these stations are depicted in Fig. 2.1. Within the NL network, not all
stations are operated by KNMI. Data from stations operated by the industry (Vermilion,
TAQA, Nouryon and NEDMAG) is flowing to the KNMI in real time and further distributed
within the NL network. Most of the (former) NAM stations have been transferred to KNMI.
In the first iteration, all stations with a temporary nature were removed. From the overview,
this only concerns the NARS network stations, which is a research network that moves every
12

2.1. LIST OF STATIONS

KNMI

few years to a new location. An S-wave crustal model was obtained with a former NARS
deployment in the Netherlands (Yudistira et al., 2017).
In the second iteration, it was decided to keep as close as possible to the station list that
is used at KNMI for detection of earthquakes. All these stations provide real-time data.
Moreover, for the KNMI detection list already a selection has been made of the stations
with the highest signal-to-noise ratios (SNRs). For instance, only the deepest of the five
sensors in a borehole setup is used during detection and is therefore included in the station
list. If sensor issues are known for this specific geophone, it will be the second deepest, etc.
Stations from the BQ network are not included in the detection list, because the data is not
available in real-time. It takes seven days before data becomes available to KNMI.
Next, the list is restricted to stations within 100 km from the Dutch onshore border. This
then includes a set of stations in Germany and Belgium that are sometimes used for location
of events in the Netherlands. Stations further away may still detect the larger events in the
Netherlands. However, these events are then already detected by a good azimuthal coverage
of stations within the Netherlands and close to border, which obviates the need for more
distant stations. Finally, sensors are removed where still multiple sensors exist at the same
(borehole) location. This is the case for many of the pre-2015 borehole stations.
In the third iteration, the list is made specific for the state of the network in 2021. This is
done by
1. Removing stations that have been closed by September 1, 2021. This concerns stations
NL.GUR1 and NL.BING.
2. Removing stations that were installed prior to 2021 but have had less then 50% data
availability over the time period September 1, 2020 until August 31, 2021. This
concerns stations NL.G024, NL.STR, NL.SUH4, NL.T014, NL.WDB4, NL.WMH4,
NL.WYN4 and GR.AHRW.
3. Removing stations that were (re-) installed in 2021 but have delivered less than two
months of data by September 1, 2021. This concerns station NL.G314.
The final list yields 200 stations with unique locations. The sensor locations are shown
in Fig. 2.2 and are listed in Appendix A. The sensor locations and specifications of the
September 1, 2022 network are listed in Appendix C.

13

2.1. LIST OF STATIONS

KNMI

Figure 2.2: Stations that are used for earthquake detection by the KNMI, within 100 km of
the Dutch onshore border (red line). The overview is from September 2021. Only stations
are included that provided sufficient real-time data over the previous year (counting from
September 2021). This same list of stations is used for computation of earthquake location
uncertainty maps. Green triangles denote broadband stations, orange triangles are shortperiod stations placed at depth (> 40 m) and purple squares are accelerometers or geophones
placed at the Earth’s surface.

14

2.2. P-WAVE GROUND-MOTION PREDICTION EQUATIONS

2.2

KNMI

P-wave ground-motion prediction equations

To determine whether earthquake waves are likely detected, a model is needed that describes the maximum event-related amplitudes as function of (hypocentral) distance and
magnitude. A model that does exactly this is called a ground-motion prediction equation
(GMPE). Detection is performed with P-waves and hence P-wave amplitudes are important
when studying detection capabilities. In the literature, however, the emphasis is on S-wave
GMPEs also sometimes called attenuation relations. S-wave attenuation models find their
use, e.g., in local magnitude determination (e.g., Kanamori et al., 1993) and S-wave GMPEs
are used especially for hazard studies (e.g., Douglas et al., 2013). P-waves typically play a
minor role in hazard studies and, as such, P-wave GMPEs usually receive little attention.
Similarly, most magnitude estimation is done with S-waves, which obviates the need of a
P-wave attenuation model.
In the following, we derive a P-wave GMPE for seismicity in the upper crust of the Netherlands, specifically for a detection frequency band of 5 to 40 Hz. In the Netherlands, detection
is done either with sensors installed at, or close to, the Earth’s surface, or with the lowermost
sensor in a vertical array of geophones. These lowermost geophones are typically installed
at 200 m depth. Therefore, two empirical relationships are derived, one for motions at the
Earth’s surface (Surface Level or SL model) and one for motions at 200 m depth (Depth
Level or DL model). Calibration is done with more than 20,000 recordings from the Netherlands Seismic and Acoustic Network (KNMI , 1993).
As a measure, the peak-ground velocity (PGV) is used. In the remainder of the report, PGV
refers to the vertical-component PGV related to incoming P waves. Since P-waves arrive
at (unconsolidated-sediment) recording sites with small angles of incidence, the verticalcomponent PGV is a good approximation of the P-wave related maximum PGV that would
be obtained from the three-dimensional ground motion.

2.2.1

Database

The ground motions database is expressed in vertical-component peak-ground velocities
(PGVs) in the detection frequency band of 5 to 40 Hz. The maximum absolute amplitude is
selected within a small time window around the expected P-wave arrival. The starting time
of this time window is based on the average velocity structure. For distances smaller than
30 km, the time window is quite short (1.0 s) to make sure that first arriving P-waves are
selected instead of, for instance, an S to P conversion. At distances larger than 30 km, the
time window is increased to a duration of 2.5 s to account for heterogeneity in the velocity
structure that can amass to sizeable timing differences between expected (mean model) and
actual P-wave arrival times. Signal time windows are only used if their SNR is equal to, or
higher than, 12 dB. This large threshold value would reject some valid earthquake ground
motions. However, because of the large size of the database, this is not a big concern.
Despite a high SNR threshold, some noise-related motions do remain in the automated
database. The SNR threshold works well in the presence of semi-continuous noise. The
SNR does not suffice when strong transient noise signals are present, which are stronger
than the earthquake P-wave arrivals in the signal time window. This especially occurs at
accelerometers, many of which are located close to local roads with infrequent car traffic. A
noise-source pollution of the signal time window may result in an unexpectedly strong PGV
value. These can be identified by computing the misfit with the P-wave GMPE. Values that
have an absolute misfit |ln(P GV ) − ln(model)| > 2, are removed, where ln is the natural
logarithm. Note that this step can only be done in a second iteration, i.e., after fitting the
first version of the GMPE. The rejection step results in the removal of 192 values (2.43%)
from the SL database, with 7709 measurements remaining. It results in the removal of 26
15

2.2. P-WAVE GROUND-MOTION PREDICTION EQUATIONS

KNMI

Figure 2.3: Induced events that have resulted in usable PGV records at a) the Earth’s surface
and b) at around 200 m depth. The size of the circles is indicative for the local magnitude,
which varies from 0.4 to 3.6. The background map is obtained from www.openstreetmap.org.
values (0.16%) from the DL database, with 15753 measurements remaining, illustrating the
lower noise levels at the deep geophones. Fig. 2.3 shows the source locations of the final PGV
database, which contains 683 and 1382 events for the SL and DL database, respectively.

2.2.2

Ground motion model

Bommer et al. (2019) developed a GMPE specifically for the Groningen gas field (further:
BMR2019). Their geometrical spreading term is partly based on numerical simulations that
use the detailed 3D velocity model that is available in the Groningen region (Romijn, 2017).
The functional form of our GMPE is largely adopted from Bommer et al. (2019). The main
difference is that we include the event depth as an additional term, similar to Atkinson et al.
(2016). Our model is parameterized as
ln(Ymod ) = c1 + c2 M + g(R∗ ),

(2.1)

where Ymod represents the modelled PGV expressed in mm/s, M the local magnitude and
g(R∗ ) a function that describes the geometrical spreading. c1 and c2 are fitting parameters.
The effective point-source distance R∗ is defined as
q
R∗ = R2 + D2 + [ee1 M +e2 ]2 ,
(2.2)
where R is the epicentral distance in km and D is the event depth in km. Eq. 2.2 contains
the magnitude-dependent distance saturation
term ee1 M +e2 with fitting parameters e1 and
√
2
e2 ; for small hypocentral distances ( R + D2 ) the exponential term is significant which
results in a non-linear scaling of Ymod with magnitude.
For the S-wave as developed in Bommer et al. (2019) the attenuation is divided in three distance ranges with different attenuation. P-waves, however, behave differently than S-waves
(as in BMR2019). Moreover, the BMR2019 model has been calibrated for the Groningen
region only, whereas we include all induced seismic events in the Netherlands and thus sample over more varying velocity structures. For our database we find that two distance ranges
suffice. That is, the misfit is not degrading when we go from three to two distance ranges.
With this simplification into two distance ranges, the attenuation terms read
g(R∗ ) = c4 ln(R∗ )
∗

∗

g(R ) = c4 ln(d) + c4a ln(R /d)
16

R∗ ≤ d
R∗ > d,

(2.3)

2.2. P-WAVE GROUND-MOTION PREDICTION EQUATIONS

KNMI

where c4 and c4a are fitting parameters. d is the effective point-source distance at which the
transition occurs from strong attenuation near the source to milder attenuation at larger
offset. In the following section the model is calibrated with the (P-wave) PGV database.
This results in an estimate of all fitting parameters that appear in Eqs. 2.1-2.3: c1 , c2 ,
c4 , c4a , e1 , e2 and d. All fitted parameters are listed in Table 2.1. More details on the
calibration procedure are provided in the next section.
Coefficient
c1
c2
c4
c4a
e1
e2
d [km]

SL model
-0.20
1.96
-3.44
-1.62
0.45
-0.80
8.0

DL model
-1.60
1.96
-3.44
-1.62
0.45
-0.80
8.0

Table 2.1: Parameters of the P-wave GMPEs for PGV values at the Earth’s surface (SL
model) and at approximately 200 m depth (DL model). Both models are for unconsolidatedsediments sites.
A model is derived for recordings at the free surface (SL model) and a model for motions
at approximately 200 m depth (DL model). Both models are the same, except for the static
offset term c1 . This term results in a 4.06 times larger PGV for the free-surface model
than for the deep model. A factor of two can be explained by the free-surface effect; at the
Earth’s surface the up-going and down-going waves are recorded at the same time, which
results in an amplitude doubling. The other factor of approximately two is caused by nearsurface amplification, which is mainly caused by a reduction of seismic impedance towards
the Earth’s surface. Fig. 2.4 shows an example borehole recording where both effects can
be seen on the first up-going P-wave arrival. van Ginkel et al. (2022b) describe the behavior
of P-waves in the unconsolidated near-surface setting in Groningen.

2.2.3

Model calibration

In Ruigrok and Dost (2020) an S-wave GMPE was derived for induced seismicity in the
Netherlands, next to the Groningen-only model (BMR2019). The Netherlands-wide S-wave
GMPE model for PGV has nearly the same functional form as the P-wave model introduced
in the previous section. Also, the geological setting is the same for both models. For these
reasons, the parameters found for the S-wave GMPE are used as a starting point for the
P-wave model calibration.
Fig. 2.5(a)-(b) shows the recorded PGV values within 2 km distance and an approximate
source depth of 3 km. Maximum P-wave amplitudes are expected within 2 km epicentral
distance. The maximum amplitude may occur at the epicentrum, it can also occur at a
larger distance, depending on the source mechanism. The model is fitted such that it well
follows the mean of the recorded values. The magnitude range of the measured PGVs is
sufficiently broad to estimate magnitude-dependent saturation (ee1 M +e2 in Eq. 2.2). It is
clear from Fig. 2.5(a)-(b) that only curved lines fit the data points adequately. The strong
saturation is probably related to the relatively high frequencies ([5 40] Hz). At epicentral
distances of 10 km the magnitude-dependent scaling is approximately linear (Fig. 2.5(c)(d)). Fig. 2.5(c) has no PGV datapoints below 1.0 µm/s. Due to high noise levels on the
surface instrument, lower ground motions have been rejected through the SNR threshold.
In the calibration process, care is taken that the models perform well for different magnitude

17

2.2. P-WAVE GROUND-MOTION PREDICTION EQUATIONS

KNMI

Figure 2.4: Seismic recording of the Westeremden M=1.7 event at November 8, 2011 at
borehole G22, which is situated at 2.98 km epicentral distance from the source. The red,
blue and green wiggles show the vertical, radial and transverse components, respectively.
To all seismograms, band-pass filtering has been applied between 5 and 40 Hz.
and distance bins. In the parameter-fitting exercise, the DL database has more weight, because it contains most PGV values. Moreover, the amplitudes recorded at depth show less
variability because the site response plays a relatively minor role at these depths. Fig. 2.6
shows the distance-dependent misfit between the PGV databases and the fitted models for a
source depth of 3 km. The DL model shows very little bias until 40 km distance. The mean
misfit from the SL model is slightly undulating as function of distance. The SL database
has a lack of measurements for epicentral distances larger than 32 km, which explains the
larger variability from that point onward. The DL database still has records beyond 120 km.
The amplitude decay at distances beyond 40 km has been solely based on the DL database.
Until 40 km the model performs very well, between 40 and 70 km the model on average
overpredicts PGV whereas between 70 and 100 km there is on average an underprediction.
Fig. 2.7 shows the estimated mean models. The model values at M=4.0 are an extrapolation
since the largest event in the database has M=3.6 (Huizinge event at August 16, 2012). The
SL model is shown together with the average 90th percentile of the vRM S noise as recorded
over surface accelerometers (Section 2.3). This gives, for each magnitude, an indication of
the distance range at which a detection is likely. The DL model is shown together with the
vRM S noise recorded at vertical component geophones at 200 m depth.
Variabilities in the model are split up into within-event variability φ andpbetween-event
variability τ . This yields the total standard deviation of the model σ = φ2 + τ 2 . The
variabilities are expressed in ln and are listed in Table 2.2.
τ expresses, among others, that there is a variability in the isotropic component of the
sources. τ is thus a property of the sources, although also the network configuration has
some influence on how well it can be estimated. For example, the network configuration
determines how well the radiation pattern is averaged out. Both the DL and SL networks
yield an independent estimate of τ . The DL dataset contains a larger amount of sources
(Fig. 2.3). Using recordings at 200 and 0 m depth, values of 0.2936 and 0.2872 are obtained,
18

2.2. P-WAVE GROUND-MOTION PREDICTION EQUATIONS

KNMI

Figure 2.5: Measured PGV values (red crosses) and fitted mean models (black lines) as
function of local magnitude. Panels (a) and (b) are for recordings within 2 km epicentral
distance (and a model for a source at 3 km depth and 1 km epicentral distance), panels (c)
and (d) are for recordings between 8 and 12 km epicentral distance (and a model for a source
at 3 km depth and 10 km epicentral distance). The left panels (a,c) show surface recordings
and models, the right panels (b,d) show recordings at approximately 200 m depth.

19

2.2. P-WAVE GROUND-MOTION PREDICTION EQUATIONS

KNMI

Figure 2.6: Misfit (ln(Yobs ) − ln(Ymod ) distributions between surface level (SL) and depth
level (DL) recordings and their respective PGV database and the surface level (SL) and
depth level (DL) GMPEs. The mean misfit (circles) and standard deviation (the bars are
from -1 to +1 standard deviation) are displayed over distance bins of 2 km for the left panels.
The right panel shows distance bins of 10 km. Only for the DL database there is a sufficient
number of recorded PGV values to show misfits until 120 km epicentral distance.
respectively. These values are averaged to find the τ as listed in Table 2.2.
φ expresses the variation of amplitudes for a single earthquake and distance. These variations
are largely caused by the radiation pattern of the source and variability in site amplification.
Prior to obtaining φ, first the total variability σ is determined by 1) computing the residual
between each data point and the mean model and 2) by taking the standard deviation of
this residual distribution. The within-event variability φ is then obtained by de-averaging
τ from the total variability σ. This total variability is higher for the 0 m depth level and
hence φ is larger at the Earth’s surface. This is to be expected since the variability in the
site effect is an important constituent in the SL model, which is hardly present in the DL
model.

Coefficient
τ
φ
σ

SL model
0.2904
0.5307
0.6050

DL model
0.2904
0.4650
0.5482

Table 2.2: Total σ, intra-event φ and inter-event τ standard deviations in the PGV empirical
model.

2.2.4

Quality checking

As a quality check of the model, we plot the measured and modelled PGV levels for several
selected events. An example is shown on Fig. 2.8. It shows a recent event induced at
20

2.2. P-WAVE GROUND-MOTION PREDICTION EQUATIONS

KNMI

Figure 2.7: The PGV GMPEs for (a) surface level (SL) motions and (b) depth level (DL)
motions. The models have been calibrated with seismicity between local magnitudes of 0.4
and 3.6. Dashed lines show, on (a), the average 90th percentile surface (accelerometer) vRM S
noise level and, on (b), the average 90th percentile 200 m depth (geophone) noise level.

21

2.2. P-WAVE GROUND-MOTION PREDICTION EQUATIONS

KNMI

the Groningen gas field. Nearly all measured PGVs lie within the 98th percent confidence
area of the SL model (which is the area in Fig. 2.8(a) from the 1st to the 99th percentile),
whereas all measured values at 200 m depth lie within the 98th percent confidence area of
the DL model (Fig. 2.8b).
The derived GMPEs have been calibrated with events inside and outside Groningen (Fig.
2.3). Most of the measurements, however, are from Groningen events, which could result in
a bias to Groningen-type events and subsurface conditions. If there is a profound difference
in induced events inside and outside Groningen, this would result in the models being less
suitable outside of Groningen. In Fig. 2.9 we check the model performance for events
outside of Groningen. Non-Groningen data points are shown in black; on the background,
Groningen-event data points are shown in red. The GMPEs are displayed for M=1.5. All
data points have been mapped to M=1.5. This is done by
1. computing Ymod (Eqs. 2.1–2.3) for the actual distance and magnitude pertaining to
each data point
2. computing the Ymod in case it were an M=1.5 event, and
3. subtracting the difference in the values computed in Step 1&2 from the measured
PGV.
Fig. 2.9 shows that both the Groningen and non-Groningen data points appear well described by the derived GMPEs. Groningen and non-Groningen data points have a similar
(distance-dependent) mean and variance.

22

2.2. P-WAVE GROUND-MOTION PREDICTION EQUATIONS

KNMI

Figure 2.8: Recorded vertical-component PGV values (red crosses) and the GMPE percentiles (dashed lines) for (a) sensors at the Earth’s surface and (b) sensors at depth. Data
is shown for the Zeerijp M=2.51 event at October 4, 2021.

23

2.2. P-WAVE GROUND-MOTION PREDICTION EQUATIONS

KNMI

Figure 2.9: Recorded vertical-component PGV values mapped to M=1.5 using the derived
P-wave GMPEs (Eqs. 2.1–2.3 and Table 2.1) for (a) recordings at the Earth’s surface and
(b) recordings at 200 m depth. At the Earth’s surface there are 7592 recordings of Groningen
events (red crosses) and 117 recordings of non-Groningen events (black crosses). At 200 m
depth there are 15078 recordings of Groningen events (red crosses) and 675 of non-Groningen
events (black crosses). The mean PGV models (GMPEs) are shown in grey.

24

2.3. STATION-SPECIFIC NOISE DISTRIBUTIONS

2.3

KNMI

Station-specific noise distributions

There are many different types of sources that produce seismic waves: cars, trains, explosions, wind turbines, generators, trees moving back and forth in a wind field, earthquakes,
etc. From these different sources, a part is classified as transient sources, meaning that they
have a short duration. Sources with a long duration, on the other hand, are called noise
sources. E.g., a car passing by a seismic sensor can be seen as a transient source near this
sensor. If the sensor is standing close to a highway, there will be a continuous succession of
transient sources that merge into a constant source. As such, the highway is a noise source.
In the current work, we are only interested in earthquakes. In this context, all other sources
that produce seismic waves are characterized as noise.
In this study, we are only interested in vertical-component noise levels, because only the
vertical component is used for detection. The noise levels are expressed in root-mean-square
particle velocities vRM S . We are especially interested in the 90th percentile levels (P90)
levels. These are the levels are exceeded 10% of the time and which are used later on for
generating the maps. In the following we derive, in a few steps, the P90 vertical-component
vRM S levels of all stations on the detection list.

2.3.1

Power-spectrum density

Seismic sensors in the NL network record continuously. The recorded wavefield is quantified
by computing the power-spectrum density (PSD) (e.g., McNamara, 2004). At KNMI, PSDs
are routinely computed, for example to quality check the instrumentation (Koymans et al.,
2021). The approach described in the aforementioned references is used to compute and
efficiently store the PSDs. These PSDs form the basis of our noise analysis.
The database described by Koymans et al. (2021) has since been updated to a more
precise version. The former database was computed using a 1-octave smoothing window,
sampling the decibel dimension with steps of 1 dB and sampling the frequency dimension
with steps of 1/8 octave. The new database is computed using a 1/8-octave smoothing
window, sampling the decibel dimension with steps of 1 dB and sampling the frequency
dimension with steps of 1/24 octave. This updated database is used in this study.
First we retrieve continuous vertical-component time-amplitude recordings, remove the
instrument response and transform the data to particle acceleration: a(t)[m/s2 ]. Second, we
use the Fast Fourier Transform (FFT) to compute the one-sided (positive frequency only)
amplitude spectrum |A(f )|. Subsequently, the PSD is obtained with
P SD(f ) = 2|A(f )|

∆t2
,
T

(2.4)

where ∆t is the sample duration of the seismic recording and T is the duration of the seismic (noise) record. The factor of two is needed when only positive frequencies are kept
with |A(f )|. To compute a stable PSD, it is averaged over multiple partly overlapping time
segments. Additionally, smoothing over the frequency dimension is applied. This results in
one estimate of P SD(f ) for each hour. The process is repeated for all hours in the day,
using a 50% overlap. Hence, 48 P SD(f ) functions are found (using also the first half hour
of the next day). These 48 PSDs give the distribution of frequency-dependent noise levels
over that day. By up-scaling this approach, also the distribution of PSDs can be found for
a week, a month and a year of data. From these distributions, the probability of occurrence
of a certain power level can be computed for each frequency. This yields the probability
density function (PDF) of the PSD.

Fig. 2.10 show examples of PDFPSDs for one year of data. The sites chosen show broadly
the noise level variation over the Netherlands.
25

2.3. STATION-SPECIFIC NOISE DISTRIBUTIONS

KNMI

Figure 2.10: PDFPSDs (probability-density functions of power-spectrum density) for the
vertical component at stations (a) VPR (surface accelerometer at Rotterdam harbour), (b)
G204 (borehole geophone at 200 m depth in Groningen) and (c) TERZ.01 (broadband sensor
at 250 m depth in South Limburg). Data from 2020-09-01 until 2021-08-31 is shown. From
the seismic noise distribution the 10th percentile and 90th percentile are depicted with black
lines, the white dotted line denotes the mean (P50). The grey lines denote a higher and
lower noise reference model as derived in Peterson (1993).
26

2.3. STATION-SPECIFIC NOISE DISTRIBUTIONS

KNMI

The upper plot is from an accelerometer that is located at the surface at an industrial
site in Pernis, within the Rotterdam harbour area. At this site, ground vibration levels are
very high due to nearby roads (A4 within a few hundred meters), water ways (Nieuwe Maas)
and large oil refineries. It can be seen that for frequencies larger than 0.4 Hz, the average
noise condition (white-dotted line) is higher than the upper grey curve, which depicts the
higher noise reference model from Peterson (1993). This shows that noise conditions in the
Netherlands can be very severe in comparison with a worldwide distribution of sites with
(broadband seismic) instrumentation.
The middle plot is from a sensor in Groningen, at 200 m depth, recording much lower
vibration levels. A large part of the anthropogenic noise (> 1 Hz) does not reach the 200 m
depth level.
The lower plot is from a broadband station in the south of Limburg, at 250 m depth
and situated within hard rock. At this site, very low background ground-motion levels are
recorded, which makes it quite ideal to detect (low-magnitude) seismicity.
The lower black lines in Fig. 2.10 denote the 5th percentile (P05) of the PDF; the upper black lines denote the 95th percentile (P95). At frequencies where both lines have a
large separation, there is much variation in noise level over the year. Fig. 2.10 also serves
to show the main differences between the sensor capabilities. Both the accelerometer and
the geophone (upper two panels) reach their recording limits around 0.2 Hz. Below this
frequency, the recording is dominated by instrumental noise. This instrumental noise is
nearly constant over the year. This results in the P05 and P95 lines coming close together
and the PDF showing a large probability (red) for a single power level. The broadband
station (lower panel) still records well the actual ground-motion below 0.2 Hz, resulting in
a continuing wide band of power levels below 0.2 Hz. The geophone and the broadband
station show a more confined distribution at the upper end of the frequency range. This is
due to the 50 Hz grid distortion that is picked up at many instruments.
The PSD is computed over a broad frequency band. For detection, a more restricted frequency band is used in which a high SNR may be expected for local and regional seismicity.
For borehole geophones and (surface) accelerometers, 4th-order Butterworth filters are applied to select a band between 5 and 40 Hz. In the current operational chain, the recording
from many of the broadband stations is band-pass filtered between 0.7 and 5.0 Hz prior
to running a detection algorithm. These broadband stations are primarily located in the
south (Fig. 2.2) with the aim to detect seismicity of a tectonic origin, both inside and
outside the Netherlands. When more geothermal operations are planned in the south, a
parallel production chain will be activated to gear the broadband stations also towards detecting induced seismicity, while keeping the ’tectonic’ chain running at the same time. In
the following, we use the ’induced’ detection band-pass filter from 5.0 to 40 Hz for all sensors.

2.3.2

Root-mean-square velocity

The operational PSDs are expressed in squared acceleration per Hz and plotted in dB:
10log10 [(m/s2 )2 /Hz]. For detection, the time-amplitude recording is used which is expressed
in particle velocity (m/s). In Section 2.2 models were derived that estimate a PGV level,
which is also expressed in particle velocity. Thus, in order to find whether an incoming Pwave likely stands out from the noise, the same unit needs to be used, i.e. particle velocity.
First, the PSD is expressed again in a linear scale: (m/s2 )2 /Hz. Subsequently, a division by
ω 2 is needed to go from acceleration squared to particle velocity squared, then an integration
over the relevant part of the frequency domain is done, and a square-root to go back from

27

2.3. STATION-SPECIFIC NOISE DISTRIBUTIONS

KNMI

Figure 2.11: Spectrogram showing one week of power spectrum densities for the vertical
component of geophone T061. Data is shown for Monday the 14th of January until Sunday
the 20th of January, 2019. UTC time used.
a power to an amplitude value:
sZ

f2

vrms =
f1

P SD(f )
df ,
ω2

(2.5)

where vrms is the root-mean-square particle velocity. A similar formula can be found in
Bormann (1998). The main difference is that we do not include a factor of 2 in Eq. 2.5
since we already implemented this factor in Eq. 2.4 to compute the PSD from a one-sided
FFT.
The minimum noise that is to be expected is the instrumental noise. The newer geophones
in the network (>2014) are of the type SM6H. In the detection frequency band, their instrumental noise level is 0.008 µm/s. The newer accelerometers that are used in the network
are of type Etna2. Their instrumental noise (including digitization etc.) is 0.048 µm/s in
the 5 to 40 Hz band. These values are found by substituting the instrumental noise models
of these instruments expressed in P SD(f ) (as published in Koymans et al. (2021)) in Eq.
2.5.
In the 5 to 40 Hz band, most of the noise has an anthropogenic origin. A large part can
be attributed to traffic, which has a clear diurnal pattern as shown in Fig. 2.11, with quiet
hours during the night and higher noise conditions during the day. For some stations, the
weekend days look markedly different from week days. For that reason, at least one week
of data is required to capture the variability in on-site vibration levels. One week of data
does not suffice, however, for stations close to wind turbines or large trees. Wind conditions
show a seasonal variability. By sampling both a more quiet and a more windy season, much
of the resulting variability in vibration levels is captured. For that reason, we use at least
half a year of data. For most stations, we use (almost) one year of data for quantifying the
seismic noise distribution (PDFPSD), from which we compute e.g. a mean vrms .
There are currently seven broadband stations not having a channel sampling the seismic
record at 100 or 200 Hz. These stations in the NL network (ARCN, HGN, HRKB, MAME,
OPLO, ROLD and WTSB) do have a channel that samples at 40 Hz, yielding a Nyquist
frequency at 20 Hz. Due to the anti-alias filter and its effect to lower frequencies (due to
28

2.3. STATION-SPECIFIC NOISE DISTRIBUTIONS

KNMI

Figure 2.12: PDFPSD for the year 2020 for the vertical component at station TERZ (lower
sensor). The PDFPSD is expressed in terms of particle velocity instead of particle acceleration (Fig. 2.10c). From the seismic noise distribution the 10th percentile and 90th percentile
are shown with black lines, the white dotted line denotes the mean (P50). The grey lines
denote a higher and lower noise reference model as derived in Peterson (1993).
smoothing when computing PSDs) the frequency band for vrms computation is restricted
to 13 Hz. Stations G81B, TERZ.01 and DBN are used to test the influence of the restricted
P 90
. For the year
frequency range on the 90th percentile of the root-mean-square velocity vrms
2020, noise-level reductions of 3.7, 10.6 and 3.4% are found, respectively.
In Fig. 2.12 a PDF is plotted in terms of particle velocity (instead of particle acceleration). It shows that most of the power resides at frequencies lower than 13 Hz. This
P 90
is reduced only with a modest percentage when the frequency range
explains why the vrms
13 to 40 Hz is left out from the 5 to 40 Hz band. Although the frequency limitation does
reduce the noise level, we do not to correct for this. The reason is twofold. On the one
hand, it is unknown what the proper correction for the few different broadband sites should
be. On the other hand, the effect is limited.
For all stations the vrms is determined using data of September 1, 2020 until August 31,
2021. Stations that have less than 50% data availability during this time period, are left
out from the station list (see Iteration 3 in Section 2.1). These stations would have a less
representative vrms . Yet, more importantly, these stations could not be used for detection
and location for more than 50% of the time.
For stations installed in 2021, the station-rejection criterion is less strict. The new
stations are still included when they have produced at least two months of data. New
stations that are hence included for the 2021 version are BLIJE, DON, NLDW, ZDL, ZH04.
Note that the vrms of especially these stations might change somewhat in 2022 when more
data becomes available.
For all stations from the list (Section 2.1) we compute and store the P05, P10, P50, P90 and
P95 levels of vrms . Fig. 2.13 shows histograms of the resulting values. The accelerometer and
geophone distributions are clearly separated. On average, noise levels on accelerometers are
much higher. This louder noise condition can be ascribed to placement of the accelerometers
at the Earth’s surface, whereas all the geophones in the detection list are deeply buried,

29

2.3. STATION-SPECIFIC NOISE DISTRIBUTIONS

KNMI

Figure 2.13: Vertical-component vRM S distribution for stations from the detection list (Fig.
2.2), with the exception of stations J01 and J02. In (a) the 10th percentile of the rootmean-square particle velocity is shown, in (b) the 50th percentile (mean) and in (c) the 90th
percentile.

30

2.3. STATION-SPECIFIC NOISE DISTRIBUTIONS

KNMI

most of them at 200 m depth. The surface short-period sensors J01 and J02 have been
left out from the distribution. The broadband sensors show a mixed picture. The most
P 90
quietly situated is TERZ.01, which is buried deep and in hard-rock (vrms
= 0.004 µm/s,
Fig. 2.10a). The highest seismic noise levels, at a broadband station, are observed at DBN,
which is at the Earth’s surface, on top of soft sediments and close to several busy roads
P 90
(vrms
= 1.812 µm/s). The overall highest noise levels are observed at a station in Pernis
P 90
= 19.55 µm/s, Fig. 2.10c) which is situated within a few hundred meters of an oil
(vrms
production plant, a refinery, a logistical hub and a highway.
Fig. 2.13(a) shows the P10 vrms distributions. The noise levels on the geophones do not
show a normal distribution, but a skewed distribution with a large peak around 0.01 µm/s.
This can be ascribed to the instrumental-noise limit of geophones of 0.008 µm/s. This implies
that for some of the instruments, for a few percent of the time, not the environmental noise
is the limiting factor for measuring very small earthquakes, but the instrumental noise.
Fig. 2.13(c) shows the P90 vrms distributions. Both the accelerometer and geophone
noise levels over the network show approximately a log-normal distribution. For the 61 accelerometers, the logarithmic mean, plus-minus one standard deviation is 0.423±0.350 log10 (µm/s),
which corresponds to 2.646 µm/s for the mean. The logarithmic mean over the 119 geophones, plus-minus one standard deviation is −1.063±0.308 log10 (µm/s), which corresponds
to 0.086 µm/s for the mean. Thus, on average there is a factor of 30.6 difference (29.7 dB)
in seismic noise level between the surface accelerometer and 200 m depth geophone installations, for the 5 to 40 Hz frequency band.
Fig. 2.14 shows the spatial distribution of vrms . Because of multiple orders of magnitude in
which noise levels vary, the amplitude scale (colorbar) is chosen as the 10-based logarithm
of the noise levels. Again, mainly a difference can be seen between deep geophone sensors
(denoted with triangles) and surface accelerometer installations (squares). A part of the 23
broadband stations are in a dedicated vaults and in a hard-rock setting, resulting in low
noise conditions. Broadband stations in Utrecht, the east of Noord-Brabant and Limburg
(DBN, OPLO, ARCN, HRKB) are on top of soft sediments and have higher seismic noise
levels. Also station ROLD has quite high noise levels, which can likely be attributed to the
urban setting of this station. The Groningen geophone network (tight triangle distribution
at the northeastern edge of Fig. 2.14) has a dense spatial distribution. This makes it possible to notice regional differences in noise levels. The region around Groningen city and the
Eemshaven is significantly louder than other parts of the network. In the Groningen city
area, roads are the main source of seismic noise. In the Eemshaven area, wind turbines are
the main cause of near constant seismic vibrations.

2.3.3

Noise reduction with depth

The accelerometers and geophones in the detection list are on different locations. Most of
the geophone stations have been chosen outside urban areas, where there is more space to
flush a temporary borehole. Some of the accelerometer sites have been chosen on purpose
within urban areas, to have a measurement of peak-ground acceleration (PGA) and PGV
that buildings in the area experience. Therefore, the accelerometers are likely on average
on more noisy locations than the geophones. As a consequence, the noise reduction values
found in the previous section are likely not caused only by differences in depth. In the
following we isolate the depth dependence on the noise by analyzing only data from the
G-network stations. At each such station there is an accelerometer at the Earth’s surface
and a string with geophones with nodes at 50, 100, 150 and 200 m depth (Dost et al., 2017).
We select one year of vertical-component PDFPSD data for all 69 borehole stations of
the G-network. Numbering in the G-network goes from G01 to G70, but G15 has never
been constructed. G010 denotes the 0 m depth level at station G01 and G011-G014 denote
31

2.3. STATION-SPECIFIC NOISE DISTRIBUTIONS

KNMI

Figure 2.14: Spatial distribution of vRM S . The 90th percentile of the vertical-component
is shown, for stations from the detection list (Fig. 2.2). Accelerometers and surface shortperiod instruments are shown with squares, deep geophones with triangles and broadbands
with circles. Only the names of the broadband stations are indicated. Names of other
stations can be found at http://rdsa.knmi.nl/network/NL/. Thin black lines denote
province borders within the Netherlands. The coordinates are in the local rectangular grid
(Rijksdriehoekstelsel) in kilometers.
32

2.3. STATION-SPECIFIC NOISE DISTRIBUTIONS

KNMI

P 90
the 50 to 200 m depth levels. As in the previous section, vrms
are determined in the 5 to
40 Hz frequency band. Stations are included in the analysis when at least 50% of the data
is available from September 1, 2020 until August 31, 2021. Four vertical-component sensors
have issues (G050, G061, G320 and G330) or insufficient data over 2021 (G190 and G31).
Therefore, data of all depth levels at G05, G06, G19, G31, G32 and G33 are excluded in the
P 90
analysis, leaving 63 stations with vrms
values established at all five depth levels.

Figure 2.15: Vertical-component vRM S levels for 63 stations of the G-network. The distribution is shown of the 90th percentile level, computed over one calendar year of continuous
recordings, for recordings at 0 m (green) and 200 m (yellow) depth.
P 90
Fig. 2.15 shows vrms
distributions for the 63 stations included in the analysis, for two depth
levels. The noise levels at the Earth’s surface and at 200 m depth are approximately lognormal distributions. The logarithmic means at 0 and 200 m are 2.293 and 0.088 µm/s,
respectively. This yields a noise-reduction gain of 25.9 (corresponding to 28.3 dB). As anticipated, this value is smaller than the factor of 30.6 noise reduction we found in the previous
section between accelerometers at the Earth’s surface and geophones at 200 m depth in the
detection list. Table 2.3 lists the average noise values at all different depth levels, as well as
the noise-level reduction gains when installing the sensor at depth. Note that these values
are only representative for unconsolidated-sediments sites and the detection frequency band.

33

2.3. STATION-SPECIFIC NOISE DISTRIBUTIONS

Depth [m]
0
50
100
150
200

µ
0.360
-0.696
-0.863
-0.966
-1.053

σ
0.324
0.299
0.305
0.292
0.292

10µ [µm/s]
2.293
0.201
0.137
0.108
0.088

KNMI

Gain [-]
1.0
11.4
16.7
21.2
25.9

Gain [dB]
0.0
21.1
24.5
26.5
28.3

Table 2.3: Average vertical-component 90th percentile vrms levels over different depth levels of G-network stations, in the detection frequency band. Also the standard deviation
per depth level σ is listed. µ and σ are the log-normal mean and standard deviation
[log10 (µm/s)]. The fourth column shows the log-normal mean expressed in µm/s. The
last two columns show the noise-reduction gain by placing the sensor at depth with respect
to placement at the Earth’s surface, expressed as a dimensionless factor and in decibels.

34

2.4. PICKING THRESHOLD

2.4

KNMI

Picking threshold

At the KNMI, continuous data is flowing in from sensors of the NL network and nearby
countries. An automatic detection algorithm is running on data streams from sensors on
the detection list (Section 2.1). First it is identified whether there is a sudden rise in
amplitudes. If this is found, the precise onset of this rise is found with another algorithm.
The timing of this onset is called a pick. An automatic event-location module is only started
up when picks are made at multiple sensors within a certain time window. In this section, we
determine the picking signal-to-noise ratio (SNR) threshold. This is the average minimum
factor needed between signal and noise for an automatic detection to be made. We use the
following definition of SNR:
SN R = 20log10

max(|s(t)|)
,
max(|n(t)|)

(2.6)

where max(|s(t)|) is the maximum absolute amplitude of the signal (within a certain time
window). In the current application that would be synonymous to the vertical-component
PGV (Section 2.2.2). max(|n(t)|) is the maximum absolute amplitude of the background
noise. In Section 2.3.2, we expressed the measured noise as root-mean-square values vrms .
With the assumption that the noise is sinusoidal in nature so that vrms = 0.707max(|n(t)|),
the minimum amplitude needed to make a pick can be expressed as
P GV =

10Φ/20
vrms ,
0.707

(2.7)

where Φ is the automatic picking SNR threshold in decibel.
An automatic detection is made with a trigger algorithm. The STA/LTA algorithm is used
for the detection of induced earthquakes in the Netherlands. STA/LTA is an abbreviation
of short-time-average over long-time-average (e.g., Havskov and Alguacil , 2016). With the
SeisComP (GEMPA, 2022) implementation at KNMI, the STA is the average absolute amplitude over the last 0.5 s of waveform data flowing in, and the LTA is the average absolute
amplitude over the last 10 s. An automatic detection is made when the STA/LTA reaches a
level of 3.
The STA/LTA algorithm detects a rise of amplitudes with respect to the background-noise
amplitude levels. The threshold of 3 cannot straightforwardly be translated to a SNR because the corresponding SNR depends on the mix of waveforms (signal amplitudes) that
is measured within the short-time window. We determine Φ empirically by using a large
collection of seismograms with different waveform distributions.
Fig. 2.16 shows the different steps in the derivation of the picking threshold. 77 seismograms
of the Westerwijtwerd M=3.4 event are used, recorded at 200 m depth. Due to the different
distances, azimuths, near-surface conditions, etc., the seismograms show a wide variation
in recorded waveforms. All seismograms are bandpass filtered in the detection band and
inspected. Seismograms with noise levels visible without zooming in, are disregarded. This
leaves 75 seismograms in the analysis. 15 levels of uniform noise are added to each remaining
seismogram (Fig. 2.16a), corresponding to SNR levels varying from 1 to 15 dB. For each
noisy seismogram, subsequently, the STA/LTA is computed and the first SNR level is determined for which the STA/LTA reaches the threshold of 3 (Fig. 2.16b). Repeating this
for the 75 seismograms recorded at different stations, a distribution of picking thresholds is
obtained (Fig. 2.16c). This distribution is by approximation a Gaussian distribution with
mean 9.33 and standard deviation 2.33. This corresponds to the PGV being, on average,
2.92 times higher than the (peak) noise and 4.14 times higher than the vrms (Eq. 2.7).

35

2.4. PICKING THRESHOLD

KNMI

Figure 2.16: Workflow illustration for deriving an automatic picking signal-to-noise threshold. (a) G054 recording of the Westerijtwerd M=3.4 event (May 22, 2019) with different
levels of uniform noise in the detection frequency band, yielding signal-to-noise ratios (SNRs)
ranging from 1 to 15 dB. (b) Seismogram with lowest SNR level (10 dB) for which an automatic detection (upper panel) is obtained, meaning that the STA/LTA (lower panel) reaches
a level of 3 (for details, see main text). (c) Histogram of minimum detection SNR levels at
75 different recordings of the Westerwijtwerd event.

36

2.4. PICKING THRESHOLD

KNMI

Fig. 2.17 shows, for each magnitude, the distance until which automatic picks are likely
feasible. As noise conditions vary from station to station, this distance range is station
dependent. Moreover, for each station also the noise conditions vary over time. Per station,
5 lines are shown to visualise the impact of different noise conditions.
P 05
With low noise conditions (vrms
resulting in the the dark green lines) the epicentral
distance at which a pick can be made is considerably larger than for high-noise conditions
P 95
resulting in the yellow lines). Making a pick at a distance larger than indicated by the
(vrms
P05 line is very unlikely, because noise conditions lower than P05 are very unlikely. On the
other hand, having a non-pick at distances smaller then indicated by the P95 line is very
unlikely, as noise conditions worse than P95 are very unlikely.
The most important line for our analysis is the P90 model as we are using this model to
decide (in the next chapter) whether a pick can, or cannot be made. To find out whether
the model predictions make sense, we compare them with actual picks. Fig.2.17 shows,
additionally, distances and magnitudes combinations where detections could be made (black
crosses) and where not (grey circles), using the STA/LTA algorithm. When checking the
picks that could actually be made (black crosses) there is indeed no pick at distances larger
than forecasted with the P05 lines. Most picks are present below the P50 lines. Below the
distances indicated by the P90 lines, almost all picks can be made (i.e. many crosses and
little circles). This shows that the P90 model yields little false positives. That is, it is only
rarely the case that this model predicts that a pick can be made, while in reality it cannot
be made.

37

2.4. PICKING THRESHOLD

KNMI

Figure 2.17: Distance-versus-magnitude picking models, together with actual picks (black
crosses) and non-picks (grey circles). The models are shown for five different noise levels,
represented by the percentile levels of the noise distribution: P05, P10, P50, P90 and P95.
At distances larger than indicated by the models (colored lines), likely no pick can be
made, given the specific noise condition. At smaller distances, a pick is likely feasible. The
model predictions are compared with distance-magnitude combinations at which a pick could
actually be made, with a STA/LTA algorithm, on recorded Z-component waveforms (black
crosses) or not (grey circles). The comparison is shown for four recording sites at different
distances from the centre of gravity of induced events (which lies in Groningen). Models
and picks at 200 m depth geophones are used as they have less false-positive picks due to
transient noise. Picks and non-picks are shown for waveform data from the years 2016, 2017
and 2018. At sensor T044, the P05 and P10 lines are on top of each other.

38

2.5. HARD-ROCK CORRECTION FACTORS

2.5

KNMI

Hard-rock correction factors

In Section 2.2 two P-wave GMPEs were developed, one for motions at the surface (surfacelevel or SL model), and one for motions at about 200 m depth (depth-level or DL model).
Both models have been calibrated with measurements in the north of the Netherlands. All
the used sensors are situated in, or on top of, unconsolidated sediments. Hence, these models
are not valid for hard-rock sites. Factors need to be found that express how much weaker
ground motions are on typical hard-rock sites in the Netherlands than on soft-rock sites (in
the 5-40 Hz frequency band). Here we determine two hard-rock correction factors, one for
the deep- and one for the surface-level motions. Fig. 2.18 illustrates the rationale behind
these factors.
The average P-wave velocity at 200 m depth in Groningen is vsd =1900 m/s (Hofman et al.,
2017) and density is ρsd =1700 kg/m3 (Kruiver et al., 2017). The subscript sd is used to
denote it are values for the typical unconsolidated sediment. For non-Groningen northNetherlands buried sensors, similar Neogene sediments are present and hence similar elastic
conditions exist at 200 m depth (van Ginkel et al., 2022a). In the Netherlands, hard-rock
conditions are only found within 200 m depth in the Achterhoek (Fig. 2.14, area around station WTSB) and in the south of Limburg (area around station VKB and further south from
this station). For hard-rock conditions at depth (sensors buried more than 40 m deep) we
assume a sandstone with a P-wave velocity vhr =3700 m/s and a density ρhr =2200 kg/m3
(e.g., Romijn, 2017). Assuming further that the impedance difference is reached over a
gradient, the hard-rock correction factor Ψ can be found with
r
ρhr vhr
.
(2.8)
Ψ=
ρsd vsd
Filling in the aforementioned values results in ΨDL = 1.6, where the subscript DL shows
that this is for deep-level recordings. This hard-rock correction factor is applied to all deeply
buried Belgium stations in the detection list (BOST, OPTB, RCHB and UCC) and TERZ.01
in the Netherlands. Hence, the modeled motion (Ymod in Eq. 2.1) is divided by ΨDL to
obtain an equivalent PGV at a deep hard-rock site.
In Section 2.2.2 it was derived that the SL model has amplitudes 4.06 times larger than
the DL model. At hard-rock sites, the amplitude difference between DL and SL motions is
different. The free-surface effect remains present, yielding a factor of two. The amplification
caused by impedance reduction towards the surface and resonances of near-surface layers,
is mostly less at hard-rock sites (without a severely weathered layer near the surface) than
at soft-rock sites. Therefore, a total factor of 2.5 is assumed between the hard-rock depthlevel motion and the hard-rock surface-level motion. As a consequence, there is a factor
of 2.6 difference in amplitudes between the SL sedimentary sites and the hard-rock sites:
ΨSL = 2.6. This factor is implemented for stations at, or close to, the surface in Belgium
(BEBN and MEM), the German broadband stations in the detection list (IBBN and BUG)
and a few NL broadband stations (WTSB, VKB, MAME and HGN).
For station NL.HGN we assess the impact of applying the hard-rock correction factor. Fig.
2.19 shows distance-magnitude combinations for which automatic picks could, or could not,
be made. The station is a few hundred kilometers south from induced seismicity in the
north of the Netherlands. In Fig. 2.19(a) the distance-magnitude picking models (Section
2.4) are plotted without applying ΨSL = 2.6. Beyond the P90 curve there are 16 picks and
12 non-picks. This means that this model is much too optimistic in forecasting that a pick
can be made. Fig. 2.19(b) shows the models after applying the hard-rock correction factor.
Beyond the P90 curve, there are now 10 picks that could be made, and two that could not
be made. Hence, this model yields a more realistic forecast of the likelihood that a detection
can be made.
39

2.5. HARD-ROCK CORRECTION FACTORS

KNMI

Figure 2.18: Sensors (green triangles) at hard-rock sites experience smaller ground motions
than at unconsolidated sedimentary sites, for the same earthquake (blue star). This is due
to larger seismic impedance (darker grey colour) at the hard-rock sites together with smaller
near-surface amplification. The figure shows factors used to translate the derived models
for sedimentary sites (Section 2.2), to hard-rock sites (at horizontal arrows) and the mean
amplification factors in both settings (at vertical arrows).

Figure 2.19: Distance-versus-magnitude picking models, together with actual picks (black
crosses) and non-picks (grey circles) for (a) station HGN and (b) station HGN after applying
the hard-rock correction factor. Induced events from the years 2016—2021 are used with
M ≥ 2.0. An STA/LTA algorithm is used for detection. The P05 and P10 models are
visually identical for station HGN.

40

2.6. TRAVEL-TIME MODELS

2.6

KNMI

Travel-time models

In this section we discuss travel-time models that well fit the observations from induced
seismicity in the Netherlands. Like for the P-wave GMPEs in Section 2.2, both mean models and their variability (or uncertainty) can be derived. The mean models are useful for
different applications. E.g, for defining time windows in which the direct P- or S-wave arrival is expected. In this report, however, we are especially interested in the uncertainty of
the travel-time model. It is the uncertainty in travel time that is mapped to uncertainty
in source location. Hence, the travel-time uncertainty is one of the main ingredients for
deriving the earthquake uncertainty maps,
For locating earthquakes, P- and/or S-wave velocity models are needed. In the operational practise at KNMI, 1D velocity profiles are used and —for the Groningen region—
hypocenters are fixed at 3 km depth. After picking arrival times of P- and S-waves (tobs )
from an induced event, the actual location algorithm is started up. The algorithm uses
pre-computed travel times (tmod ). E.g, raytracing is used to derive travel-time models from
the velocity models. Subsequently, it is assessed which source location yields the best fit
between the observed and modeled travel times. Even for the best-fitting source-location,
there is a misfit between the modeled travel times and the picks. The remaining residuals
(tobs -tmod ) are a measure of the uncertainty in travel times. The travel-time model cannot
perfectly explain the observations. Possible causes are related to:
• Lateral heterogeneity. In the assumed model there are no lateral variations in
velocity structure. In reality there are, of course.
• Picking errors. The picked arrival times are not necessarily the true arrival times
of first P- or S-waves. In the presence of noise, it could be difficult to pick the precise
onset of the wave. Furthermore, it could be that a secondary arrival was picked instead
of the assumed first arrival.
In this section we derive empirical travel-time models, and their uncertainty. We take
advantage of the rich database of induced events in the Netherlands. The advantage of
these empirical relations is that they are intrinsically based on arrival times that can be
well picked. Since the models are directly expressed in travel times, instead of velocity, no
forward modelling is needed to derive travel times. Hence, numerical modeling error1 is
left out of the equation. We estimate travel-time uncertainties for two situations; for the
situation that there are not sufficient local picks and one needs to resort to a background
travel-time model (Section 2.6.1) and for the situation that there is sufficient data to extract
an event-specific travel-time model (Section 2.6.2). In the last section, a travel-time model
for the near range is discussed.

2.6.1

Background polynomial model

For three years of data (2016-2018) we select all induced events with magnitudes of 1.5 and
higher. This includes events from Groningen, Friesland, Drenthe and North Holland. We
use the operationally hand-picked P-wave arrival times (3048 picks). This gives an overview
of actually observed arrival times as function of epicentral distance (Fig. 2.20). For all
these events, the approximate source depth is 3 km. Subsequently, these arrival times are
least-squares fitted with a polynomial of degree 2:
TP = 0.725 + 0.204R − 2.55 · 10−4 R2 ,

(2.9)

1 Numerical techniques to forward model arrival times have their limitations. E.g., raytracing cannot
forward model all possible phases. With a more advanced technique, phases might be forward modeled that
in reality are too weak to be picked.

41

2.6. TRAVEL-TIME MODELS

KNMI

Figure 2.20: Ppick distributions for induced events in 2016-2018. The upper panel shows
the Ppicks (blue circles) as function of epicentral distance and time, together with a linear
travel-time model (straight red line) and a degree-2 polynomial least-squares fit (curved
yellow line). Until about 40 km the travel times can be well explained with a linear model.
The lower panel shows the residual distribution with respect to the degree-2 polynomial
model.

42

2.6. TRAVEL-TIME MODELS

KNMI

Figure 2.21: Spick distributions for induced events in 2016–2018. The upper panel shows
the Spicks (blue circles) as function of epicentral distance and time, together with a linear
travel-time model (straight red line) and a degree-2 polynomial least-squares fit (curved
yellow line). Until about 20 km distance the travel times can well be explained with a
linear model. The lower panel shows the residual distribution with respect to the degree-2
polynomial model.
where R denotes epicentral distance. The lower panel in Fig. 2.20 shows the residual
distribution, which has a standard deviation of 0.137 s and a mean that is practically zero
(6.23 · 10−11 s).
For the same collection of events, we add 2066 hand-picked S-wave travel times. Also the
S-wave picks are fitted with a degree-2 polynomial function (Fig. 2.21):
TS = 1.855 + 0.344R − 3.89 · 10−4 R2 ,

(2.10)

The lower panel in Fig. 2.21 shows the residual distribution, which has a standard deviation
of 0.248 s.
In the preceding we estimated average travel-time uncertainty for a large epicentral distance range, from 0 to 160 km. In practise, however, travel-time uncertainty is a function
of distance. At small distances, the travel-time residuals remain small as the travel times
are small. With larger distances, a few percent velocity deviation results in larger residuals.
For this reason, smaller location uncertainty can be reached when a local network is used
instead of a regional or national network.
For the quantification of the distant-dependent travel-time uncertainty we split up the distance range in three bins: local [0 20], regional [20 60] and national [60 160] km. Subsequently, using only picks that fall in the respective bins and using the derived P- and S-wave
travel-time models (Eqs. 2.9&2.10), we compute the residuals and their standard deviations.
The results are listed in Table 2.4, together with the previous result for the entire distance
range. As expected, the uncertainty increases with distance.
43

2.6. TRAVEL-TIME MODELS

KNMI

Figure 2.22: Picks of (a) P-wave arrival times and (b) S-wave arrival times for the Froombosch M=2.4 event on February 25, 2016. The picks are fitted with a degree-2 polynomial
function, from which the intercept time T0 , the apparent velocity at small offset v1 and the
standard deviation of the residuals σt are shown in the panels.

Travel-time uncertainty in Table 2.4 is determined with picks until 160 km distance. Picks
are sometimes possible for even larger distances. For example, the larger events in the
Groningen gas field can be detected at stations in the south of Limburg. These picks are
not used, however, for location. There is a sufficient number of picks available at smaller
distances with lower travel-time uncertainty.
Naming
Local
Regional
National
All

Distance range [km]
0 – 20
20 – 60
60 – 160
0 – 160

σP [s]
0.115
0.162
0.295
0.137

σS [s]
0.186
0.322
0.568
0.248

Table 2.4: P-wave travel-time uncertainty σP and S-wave travel-time uncertainty σS for four
distance bins.

2.6.2

Event-specific polynomial model

The uncertainty of the travel-time models (Eqs.2.9&2.10) is largely caused by lateral heterogeneity in the subsurface velocity structure. This uncertainty is larger when using an
average travel-time model for different regions in the Netherlands, like in the previous section. If a sufficient number of local or regional picks is available, then an event-specific
travel-time model can be derived. A good sampling of the regional velocity structure, and
its uncertainty, is achieved when the source is well surrounded by stations. As a measure
of ’well-surroundedness’ the largest azimuthal gap ψ is used. Only when ψ < 120o (Bondár
et al., 2004) and at least 12 picks are available at distances smaller than 40 km, an eventspecific model is derived.
Fig. 2.22 shows examples of event-tailored travel-time models, both a P-wave and S-wave
model. For the years 2016–2018, 53 induced M ≥ 1.5 events are available of which 50 satisfy
the pick-number and ψ < 120o condition. For each of these events, the standard deviation
of the residuals is determined. The average P uncertainty (average over the standard deviations) is 0.0899 s. The S uncertainty is again almost twice as large, with a value of 0.170 s.

44

2.6. TRAVEL-TIME MODELS

2.6.3

KNMI

Hyperbolic travel times

The above travel-time models were derived for sources at approximately 3 km depth. For
the entire distance range, it was assumed that the travel-time models can be described by
degree-2 polynomials. These polynomial models would be fair to describe refracted and
diving waves. Polynomial models are less suitable, however, to describe up-going waves
near the source. The travel-times for these waves have a hyperbolic move-out and can be
approximately described with
p
(2.11)
T = (D/v0 )2 + (R/v0 )2 ,
where D is the depth of the event, R the epicentral distance and v0 an average overburden
velocity near the source. The hyperbolic move-out describes a strong curvature of the traveltime function in the first few kilometers epicentral distance. In the following we show both
with actual travel-time picks and modelling that a hyperbolic behavior is a fair description
at near offset.
Fig. 2.23 shows picks for the near-offset range (same picks as in Figs. 2.20&2.21). The
picks show a somewhat curved line from zero to approximately 3 km epicentral distance. At
further distances ([3 20] km) the travel times are well approximated with the nearly linear
trajectory of a degree-2 polynomial travel-time model.
If the hyperbolic model (Eq. 2.11) is used to fit the picked travel times, then for the [0
20] km distance range the obtained uncertainty is similar to the polynomial fit (Table 2.4).
Only for the first few kilometers distance, the hyperbolic model yields a better description.
The variability in observed travel times in the first few kilometers (the hyperbolic part)
epicentral distance is still similar to the variability at somewhat larger distances (e.g., around
8 km distance in the linear part, Fig. 2.23).
The picks in Fig. 2.23 are compared with modelling results generated with a 1D velocity
profile at Huizinge, a small village in Groningen. This profile is extracted from a detailed
velocity model of the region (NAM 3D velocity model; Romijn (2017)). Subsequently, for
sources at varying depths and receivers at 200 m depth, the wavefield is forward modelled
with a visco-elastic finite-difference scheme (Robertsson et al., 1994). The direct P- and
S-waves are then picked with a kurtosis-based picker (Langet et al., 2014) and stored in a
travel-time table. Fig. 2.23 shows the resulting travel-time curves for sources between 1
and 5 km depth. These travel-time curves are shifted with 0.07 seconds so that the 3 km
depth curve fits quite well with the mean of the travel-time picks. Thus, the used Huizinge
region profile results in a somewhat biased travel-time model for the collection of picks from
varying regions.
For depths other than 3 km, we do not have a good empirical model. There may be enough
data points to estimate the polynomial part of the travel-time model, as in the previous
sections. However, it is unlikely that the near-offset range is well covered with seismic
stations. Generally, one needs to resort to a local 1D profile to estimate the travel times at
short offsets, and specifically to find the influence of source depth on these travel times.

45

2.6. TRAVEL-TIME MODELS

KNMI

Figure 2.23: (a) Ppick and (b) Spick distribution for induced events in 2016–2018, together
with travel-time curves. The curves have been generated for a realistic 1D velocity profile
in the Groningen region (Huizinge), for source depth varying between 1 and 5 km. By
approximation, the curves are hyperbola.

46

3 Location-uncertainty method
3.1

Bayesian framework

Different data attributes can be used to locate a source: arrival times of different phases (P,
S, Pg, PmP, etc.), polarization of P waves and relative amplitudes over the network for different phases. In this report we restrict ourselves to timing data attributes, and specifically
to the timing of the first arriving P and S waves. It is also this timing that is used in an
operational setting, and for which the uncertainty has been quantified. In Section 2.6 we
derived mean travel-time models and quantified the variability thereof. Since the residual
distribution was approximately a normal distribution, this variability could be quantified
with a standard deviation, one for P-wave timing σp and one for S-wave timing σs . These
uncertainties are used in the inversion and are placed in a data covariance matrix Cd .
The location of an event is found by taking a misfit between observed and forward modeled
travel times, or travel-time differences. When doing this in a grid search over 3D space, a
most likely location is found where the misfit reaches a minimum value. Moreover, a cloud
of points is found where the misfit is still reasonably low. Because the timing is uncertain,
locations around the most likely location cannot be discarded; they remain likely given the
chance that, e.g., the actual velocity structure is a bit different from the model that was
assumed. In a Bayesian framework, the misfit together with a data-error model are used to
compute the probability of each point being the actual source. After knowing these probabilities and normalizing the sum of all probabilities to one, source-confidence volumes can
be extracted. The assumption of normally distributed data errors allows using the following formulation from Tarantola (2005) to map errors in timing to the ’a posterori’ location
probability density function (PDF, σM (x)):
σM (x) = k exp(

−1 calc
calc
[d (x) − dobs ]T C−1
(x) − dobs ],
d [d
2

(3.1)

where k is a normalization constant, x = (x1 , x2 , x3 ) are the Cartesian coordinates spanning
up the model space and T denotes the transpose. dobs is a vector that contains the data
attributes for the different receivers, or receiver combinations. dcalc (x) is the forward modeled data attribute for a source at position x in the model. Cd is the data covariance matrix
and C−1
d is the inverse thereof. σM (x) is computed numerically by doing a grid search over
the three space variables.
In the following sections we work out Eq. 3.1 for two different data attributes, and a
combination of these two. Similar approaches of computing location error ellipsoids are
worked out in Lomax et al. (2009), Freudenreich et al. (2012) and Kinnaert et al. (2016).

3.1.1

P delay

From the travel-time measured at one receiver, an equidistant curve can be drawn. This is
a surface in 3D with possible source locations. With three receivers, three such equidistant
47

3.1. BAYESIAN FRAMEWORK

KNMI

curves can be drawn and their intersection (a single point in case there is no noise) is the
hypocenter location in 3D space (x1 , x2 , x3 ). What is measured, though, are no (P-wave)
travel times, but arrival times. The arrival time tp is the addition of the earthquake origin
time t0 and the travel-time ttp . Thus there are four unknowns, (to , x1 , x2 , x3 ) for which a
4D PDF could be constructed. However, location and time are linked variables. Given a
certain velocity model, there is one t0 that corresponds to each x. Thus, there is no need to
compute a 4D PDF. t0 can be obtained in a post-processing step.
In this report, we constrain ourselves to the location problem; we do not invert for t0 . Taking
arrival-time differences is equivalent to taking travel-time differences and it is an effective
approach to take t0 out of the equation, e.g., ∆tp21 = tp2 − tp1 = (t0 + ttp2 ) − (t0 + ttp1 ) =
ttp2 − ttp1 . Zhou (1994) used these travel-time differences over a network of stations to find
a hypocenter location and coined it the equal differential time (EDT) method. The same
data attribute was used by Spetzler and Dost (2017) for location of events in the Groningen
region. To put the EDT method in a Bayesian framework in which multiple data attributes
can be combined, in this report we refer to it as location using the P-delay data attribute.
In the following we illustrate the terms in Eq. 3.1 with the P-delay data attribute. Using a minimum of three receivers to find a PDF in x, the data attribute vector is composed
as


tp2 − tp1
d=
,
(3.2)
tp3 − tp1

Containing only the travel-time differences of the P-wave arriving at the second and first
receiver (top row) and third and first receiver (bottom row).
For location, we use the travel-time models as derived in Section 2.6. It is assumed that
the arrival times can be modeled using a Gaussian probability density distribution with a
mean travel time and standard deviation σp . For three P-wave detections we have three
measurements including noise, yielding the following noise model:
 
tp1
N = tp2 
(3.3)
tp3
and corresponding noise covariance matrix
 2
σp
Cn =  0
0

0
σp2
0


0
0 .
σp2

The matrix A that relates the data-attribute vector to the noise vector (d=AN) is


−1 1 0
A=
.
−1 0 1

Using that Cd = ACn AT , the data covariance matrix reads


2 1
.
Cd = σp2
1 2

48

(3.4)

(3.5)

(3.6)

3.1. BAYESIAN FRAMEWORK

KNMI

In the data-attribute vector (Eq. 3.2) only n − 1 linearly independent receiver combinations
are taken. Instead, one could choose to use all n(n − 1)/2 unique receiver combinations.
This results in larger matrices, and thus longer computation times. Moreover, in this case
the inverse of Cd does not exist; different rows of d can be written as a linear combination
of other rows. And hence one would need to resort to a pseudo inverse. Given the used
noise model, using n(n − 1)/2 unique station combinations yields the same hypocenter PDF
as with the n − 1 linearly independent station combinations. So, there is no advantage in
extending the data-attribute vector with linearly dependent receivers combinations and thus
we use the n − 1 implementation as outlined above.

3.1.2

P-S delay

An equidistant curve can also be constructed from the arrival-time difference of the first Pand S-wave arrival at a single receiver. With P-S delay measurements at three stations, a
location can be found in 3D which would be —if there was no noise— the intersection of
the three equidistant curves. We use the minimum amount of three receivers to illustrate
the terms in Eq. 3.1 when the P-S delay data attribute is used. The data-attribute vector
looks as follows


ts1 − tp1
d = ts2 − tp2  ,
(3.7)
ts3 − tp3
where ts and tp are S-wave and P-wave arrival times and the subscript number refers to the
receiver number.
We use the travel-time models as derived in Section 2.6. It is assumed that the arrival
times can be modeled using Gaussian probability density functions with mean arrival times
and standard deviation thereof, σp for P-waves and σs for S-waves. For three P-wave detections and three S-wave detections we have six measurements including noise, yielding the
following noise model:
 
ts1
ts2 
 
ts3 

N=
(3.8)
tp1 
 
tp2 
tp3
and corresponding noise covariance matrix
 2
σs 0
0
 0 σs2 0

0
0 σs2
Cn = 
0
0
0

0
0
0
0
0
0

0
0
0
σp2
0
0

0
0
0
0
σp2
0


0
0

0
.
0

0
σp2

(3.9)

The matrix A that relates the data-attribute vector to the noise vector (d=AN) is


1 0 0 −1 0
0
A = 0 1 0 0 −1 0  .
(3.10)
0 0 1 0
0 −1
49

3.2. UNCERTAINTY QUANTIFICATION

KNMI

Using that Cd = ACn AT , the data covariance matrix reads

 2
σs + σp2
0
0
σs2 + σp2
0 .
Cd =  0
0
0
σs2 + σp2

3.1.3

(3.11)

Joint inversion

When errors may be assumed to be normally distributed, the use of more measurements
reduces the location uncertainty. Likewise, also the usage of multiple data attributes per
station, or station-combination, may be advantageous. Many different combinations are
possible. Below we work out the combination of the P-S delay and P-delay data attribute.
These two data attributes yield location PDFs with very different shapes. That is, both
data attributes have complimentary sensitivity to the location problem. The joint inversion
thus results in a significantly smaller PDF, as will be shown in Section 3.3.
We use again the minimum amount of three receivers to illustrate the terms in Eq. 3.1.
When using both P-S delays and P delays, we have the following data-attribute vector:


ts1 − tp1
ts2 − tp2 



d=
(3.12)
ts3 − tp3  .
tp2 − tp1 
tp3 − tp1

In the above vector, the same measurements are used as in the previous section. Hence, the
same noise model (Eq. 3.8) and noise covariance matrix (Eq. 3.9) remain valid.
The matrix A that relates the data-attribute vector to the noise vector (d=AN) is


1 0 0 −1 0
0
0 1 0 0 −1 0 


0 −1
(3.13)
A=
.
0 0 1 0
0 0 0 −1 1
0
0 0 0 −1 0
1
Using that Cd = ACn AT , the data covariance matrix reads
 2
0
0
σp2
σs + σp2
2
2
 0
σs + σp
0
−σp2

2
2
0
σs + σp
0
Cd = 
 02
 σp
−σp2
0
σp2 + σp2
σp2
0
−σp2
σp2

3.2


σp2
0 

−σp2 
.
σp2 
σp2 + σp2

(3.14)

Uncertainty quantification

For finding a source location, Eq. 3.1 is implemented. The resulting solution is not just
the most likely source location, but the entire location PDF (σM (x)) which gives the probability of the source as function of a range of locations in 3D space. This PDF is —by
approximation— a multivariate normal distribution which can be described with its mean
µ=[µx µy µz ] and standard deviation σ=[σ10 σ20 σ30 ], which are the square root of the variances in the three principal directions. The mean of a normal distribution corresponds to the
location in space with the highest probability density (the mode). This most likely location
50

3.2. UNCERTAINTY QUANTIFICATION

KNMI

is communicated as the source location, and the uncertainty is communicated through the
standard deviations of the location PDF.
Fig. 3.1(a) shows an example of a 95% confidence ellipsoid. There is a probability of
95% that the source location lies within this volume. Generally, the principal directions of
the uncertainty ellipsoid (x01 , x02 , x03 ) do not correspond to the directions of the grid (X, Y, Z).
Hence, nine parameters are needed to describe the ellipsoid: three to position the center
point of the ellipsoid in 3D space, three rotation angles to map from (X, Y, Z) to (x01 , x02 , x03 )
and three sigmas to denote the size of the ellipsoid. The ellipsoid can be constructed by
1. Taking the mode of the PDF as the center point of the ellipsoid.
2. Outlining the directions of the principal axes.
3. Drawing the ellipsoid with the correct size of the principal axes.
This size is found by multiplying 2 times 2.7955 to σ10 , σ20 and σ30 . That is, for a tri-variate
normal distribution, 95% of the probability lies within 2.7955 standard deviations from the
mean (Wang et al., 2015). Conversely, e.g. σ10 can be obtained by dividing the distance
span by the 95% confidence in the x01 direction by 5.5910.
With the above parameterisation, three parameters are used to denote the most likely location and 6 to describe the uncertainty. In seismology, however, we are especially interested
in epicentral uncertainty (uncertainty in the horizontal plane) and depth uncertainty (uncertainty in depth direction). Instead of parameterising the uncertainty ellipsoid, we quantify
the epicentre uncertainty and depth uncertainty separately. This yields a simplified description of the uncertainty with only four parameters: σZ (Fig. 3.1b), σ1 , σ2 and θ (Fig. 3.1c).
The depth uncertainty is found by —at the epicentre— computing the location probability
in the vertical direction. This gives a 1D PDF as depicted in Fig. 3.1(b), which is a univariate normal distribution. Subsequently, the depth range is computed that contains 95%
of the location probability, Z 95 . In 1D, Z 95 = 2 · 1.96σZ (e.g., Wang et al., 2015), wherein
σZ is the standard deviation in the vertical direction. The above formula is often simplified
by stating the 95% of the probability is reached within two sigmas from the mean. We use
the more exact
Z 95
Z 95
=
.
(3.15)
σZ =
2 · 1.96
3.92
The epicentral uncertainty is found by computing the 95% confidence region in a X-Y
(or East-North) plane through the assumed earthquake depth. This gives a 2D PDF as
depicted in Fig. 3.1(c), which is a bi-variate normal distribution. Subsequently, the ellipse
is computed that contains 95% of the location probability:
x21
x2
+ 2 = κ2 ,
σ1
σ2

(3.16)

which ellipse has principal directions x1 and x2 . θ is the clockwise rotation angle of x1 with
respect to north (Fig. 3.1c). When the ellipse describes the bounds of the 95% confidence
region computed in 2D, k = 2.4477 (Wang et al., 2015). From this ellipse, the major axis
and minor axis are X195 and X295 , respectively. Using the above formula, it follows that the
uncertainty in the two principal directions in the X-Y plane can be found with
σ1,2 =

95
95
X1,2
X1,2
=
.
2 · 2.4477
4.8954

51

(3.17)

3.2. UNCERTAINTY QUANTIFICATION

KNMI

Figure 3.1: Location-uncertainty assessment in different dimensions. (a) 95% confidence
ellipsoid (grey shaded) in 3D shown together with the acquisition geometry of receivers
(green triangles) on the Earth’s surface (yellow). (b) Event-depth PDF, of which the 95%
confidence region is shaded blue. The arrow depicts the size of σz (Eq. 3.15). The PDF
has been computed at the epicenter. (c) Event-epicenter PDF of which the 95% confidence
region is bounded by the black ellipse. The PDF has been computed at the event depth.
σ1 and σ2 are the largest and smallest uncertainties in the epicentral plane and the angle θ
describes the directivity of σ1 with respect to north. σ2 is orthogonal to σ1 .

52

3.3. COMPARISON OF DATA ATTRIBUTES

KNMI

The following steps are taken to quantify the epicentral uncertainty in terms of σ1 , σ2
and θ:
• Compute a PDF in the 2D plane through the (assumed) depth of the source.
• Normalize the total probability to 100%.
• Sort the probabilities from high to low.
• Find the grid point, and corresponding probability, at which 95% of the probability is
reached (P95).
• Interpolate a contour through 2D space where the P95 value is reached.
• Fit the contour with an ellipse.
• Determine the major (X195 ) and minor axes (X295 ) of the ellipse and its rotation angle
with north θ.
• Use Eq. 3.17 to compute σ1 and σ2 .
In Wang et al. (2015) the ratio between the size of the confidence region and σ are called
magnification ratio’s. In 2D, the magnification ratio corresponds to κ in Eq. 3.16. Note
that these ratios are not fixed, but depend on the space dimensionality in which the PDF
is computed. The higher the dimension, the larger the confidence zone and thus the larger
the magnification ratio. This can be seen, e.g., when comparing Fig 3.1(a) and (b). In 3D
the depth extend of the 95% confidence region extends to the free surface, whereas in 1D it
ends about 1400 m below the surface.
At Fig. 3.1 it can be seen that the PDF is a (multivariate) normal distribution only by
approximation. For instance, note the (small) asymmetry of the PDF in the depth extend.
The confidence zone is larger at the shallow end of the most likely source location than at
the deep end. Due to the presence of the free surface, the mean does not correspond to
the mode. The mode still depicts the actual source depth, whereas the mean would give an
underestimate of the source depth. Therefore, the mode is used to communicate the source
location.

3.3

Comparison of data attributes

In this section we illustrate the location uncertainty that is obtained using the different data
attributes as introduced in Section 3.1. For this illustration, we use a local network with
six receivers (Fig. 3.2 left column). The source is placed at 3 km depth below the network.
There is a good azimuthal coverage. We use the low travel-time uncertainties of σp = 0.0893
and σs = 0.170 s as estimated in Section 2.6.2. Fig. 3.2 shows the resulting PDFs.
The top row in Fig. 3.2 shows the PDFs (in 3D, epicenter and depth) when using P delays.
This data attribute gives a quite confined epicentral 95% confidence zone. However, it yields
a large uncertainty in depth. This is because shifting the source location in X-Y results in
(significant) time shifts of the incoming P wave at the receivers. Changing the depth of the
source, however, leads to much smaller differences in the P-wave front that is sampled with
the local network.
The middle row in Fig. 3.2 shows the PDFs (in 3D, epicenter and depth) when using
P-S delays. This data attribute gives a less confined epicentral 95% region than when using
P delays. However, it yields a much better confined depth uncertainty. The travel-time
differences between incoming P- and S-waves contain especially information of the distance
53

3.3. COMPARISON OF DATA ATTRIBUTES

KNMI

Figure 3.2: 95% location confidence regions for a setup with six receivers (green triangles) on
the Earth’s surface (yellow surface) in 3D (left column), in the epicentral plane (middle column) and in the depth direction (right column). The top row shows results when using the Pdelay data attribute, yielding (σ1 , σ2 , σZ )=(489, 350, 1837) m and θ = 33.7o . The middle row
shows results when using the P-S delay data attribute, yielding (σ1 , σ2 , σZ )=(746, 622, 561) m
and θ = 7.1o . The lower row shows results when doing a joint inversion with both data attributes, yielding (σ1 , σ2 , σZ )=(205, 153, 434) m and θ = 31.4o .

54

3.3. COMPARISON OF DATA ATTRIBUTES

KNMI

of the source with respect to the network. When a high depth resolution is desired, it is
elementary to include this data attribute in an inversion.
The bottom row in Fig. 3.2 shows the PDFs (in 3D, epicenter and depth) when using
a joint inversion with both P- and P-S delays. The joint inversion leads to smaller PDFs
than with any of the two data attributes separately. The near-orthogonal sensitivity of both
data attributes (depth versus epicentral) are combined to obtain the best overall hypocentre
estimate. Its estimate has the smallest values for σ1 , σ2 and σZ .

55

4 Sensitivity analysis
In this chapter we assess how location uncertainty varies as a function of different variables.
The variables we consider are uncertainty of the underlying measurements (σp and σs ), systematic changes in the network layout (azimuthal gap, number of receivers, array aperture
and distance to the nearest station) and source depth. Understanding sensitivity of location
uncertainty to the different variables paves the way of finding effective means of reducing
the location uncertainty. Also, the sensitivity analyses help in understanding the details
of the Netherlands-wide location-uncertainty maps that follow in the next chapter and in
Appendices B and D.
In the previous chapter, specifically in Fig. 3.2, it was shown that a joint inversion with P-S
and P delays is most effective in constraining the source location, both in the depth extend
(σZ ) and in the epicentral plane (σ1 and σ2 ). All numerical analyses in this chapter are only
performed for this joint inversion with multiple data attributes.
The local velocity structure, and the resulting travel-time function, has an influence on
how the travel-time errors are propagated to event-location uncertainty.
On the one hand, errors are relative to an average propagation velocity. When the
propagation velocity is 6 km/s instead of 3 km/s, a 0.1 s error leads to twice the location
error. Hence, it is important to take realistic background velocities. In this section, the error
propagation is done through a homogeneous background medium with P-wave and S-wave
velocities of vP =4.9 and vS =2.9 km/s. These are typical apparent velocities as found for
upper-crustal events in the Netherlands (Section 2.6.1 and Jagt et al. (2017)).
On the other hand, complex details in the velocity structure lead to perturbations of the
location-uncertainty PDF. This last aspect is ignored in this report. The perturbations in
the location PDF can only be computed well if detailed local velocity models are available.
Moreover, these perturbations will be different in different regions, which make them hard
to include in a NL-wide framework. When possible complexities of the location PDF are
of interest, (additional) local or regional studies are needed. The modelling outcomes in
the coming sections serve to show the general trend of how location uncertainty behaves as
function of different variables.

4.1

Travel-time uncertainty

In Section 2.6 we estimated the uncertainty of the measurements that are used for location:
P-wave and S-wave travel times. In specific cases, the travel-time uncertainty may in fact be
lower. For instance, when travel times are forward modeled through a detailed local 3D velocity model. In other cases, the travel-time uncertainty may be higher. For instance, when
the source is on a fault with a large offset, with very different velocity structure on either
side of the fault. Here we assess how the location uncertainty changes when the travel-time
uncertainty is varied.

56

4.1. TRAVEL-TIME UNCERTAINTY

KNMI

Figure 4.1: Configuration of stations (orange triangles) that likely detects a M=0.5 source
placed at 3 km depth (black dot). Within the red box, a grid search is done to find the
epicentral PDF. Background map is from www.openstreetmap.org.
We use a source with M=0.5 at 3 km depth within the Twente seismic network (Fig. 4.1)
and forward model travel times through a homogeneous background medium. We use the
tools developed in Chapter 2 to identify the stations where this scenario earthquake is likely
detected: we use the P-wave GMPEs to model the amplitude at each station of the detection
list and compare the modeled amplitude with the P90 local noise conditions and account
for the picking threshold and possibly a hard-rock correction factor. For this M=0.5 event
this yields a pick list that contains the six nearest stations of the Twente array. From these
stations we use the P-wave and S-wave picks.
We first assess the situation in which σs is kept fixed at 0.248 s and σp is varied. The
resulting location uncertainties (σ1 , σ2 and σZ ) are shown in Fig. 4.2. For small σp , the
model uncertainties depend nearly linearly on σp . For larger σp , the epicentral uncertainties are reaching a plateau. The plateau starts at σp = σs = 0.0248 s. Beyond that point,
the epicentral location is predominantly constrained by S-wave travel-time information, and
since the uncertainty thereof is fixed, the epicentral uncertainty nearly stops increasing.
In Fig. 4.2 it is also shown how the epicentral (lower panels) and depth (higher panels)
PDFs look like for σp = 0.12 and σp = 0.36 s. The epicentral PDFs do not notably change
shape when σp increases. The depth PDF shows different behavior: for small σp there is
limited interaction of the PDF with the free surface. For larger σp , the PDF is bounded on
the upper end by the free surface, which limits its growth. Thus, the σZ function (yellow
line) has a higher derivative (dσZ /dσp ) for σp < 0.18 than for σp > 0.18.
We have derived in Section 2.6 that σs is about twice as large as σp for induced seismicity in
the Netherlands. We assess how the location uncertainty changes as function of simultane-

57

4.1. TRAVEL-TIME UNCERTAINTY

KNMI

Figure 4.2: Location uncertainty as a function of P-wave travel-time uncertainty σp . The
model uncertainty in the epicentral plane is expressed by σ1 and σ2 . The uncertainty in
the depth direction is given by σZ . The location uncertainty is computed using the station
configuration in Fig 4.1 and a source at 3 km depth. The S-wave travel-time uncertainty is
fixed at 0.248 s. The insets show the shape of the epicentral and depth PDFs for σp = 0.12 s
(left-hand panels) and for σp = 0.36 (right-hand panels).

58

4.2. SOURCE DEPTH

KNMI

Figure 4.3: Location uncertainty (σ1 , σ2 and σZ ) as function of travel-time uncertainty (σp
and σs ). σ1 and σ2 are varied at the same time, with σp annotated on the horizontal axis
and σs = 2σp . The location uncertainty is computed using the station configuration in Fig
4.1 and a source at 3 km depth.
ous changes in both σp and σs , when σs = 2σp . Fig. 4.3 shows the resulting curves. In this
case, the epicentral uncertainty (orange and blue lines) varies approximately linearly with
travel-time uncertainty. The depth uncertainty (yellow line) is significantly larger than the
epicentral uncertainty. σZ shows a similar character as in Fig. 4.2, with a large dσZ /dσp
derivative at small σp and a smaller derivative as soon as the PDF becomes constrained by
the free surface.

4.2

Source depth

In this section we numerically compute the influence of source nucleation depth on location
uncertainty. We use the same six stations as in the previous section (Fig. 4.1), with the
same epicentral source location and vary the source depth. The smallest (epicentral) sourcereceiver distance within the configuration is 1.82 km and the largest is 10.48 km. We use low
travel-time uncertainties of σp = 0.0893 and σs = 0.170 s as estimated in Section 2.6.2. Fig.
4.4 shows the variation of σ1 , σ2 and σZ as function of source depth.
Fig. 4.4 shows that epicentral uncertainty (σ1 , σ2 ) increases nearly linearly with source
depth. As we have seen in Section 3.3, the epicentral uncertainty is predominantly constrained by P-delays over the network. The deeper the source, the smaller the P-delay times
(that is, the more the wavefront approximates a plane wave over the network) and hence
the larger the epicentral uncertainty.
The depth uncertainty σZ shows a different character (yellow line in Fig. 4.4). As function
of depth, it peaks at shallow depths and becomes stable at larger depths. When the source
is shallow, the PDF is bounded by the free surface. This results in a reduction of the PDF

59

4.3. SOURCE-RECEIVER DISTANCE

KNMI

Figure 4.4: Location uncertainty (σ1 , σ2 and σZ ) as function of source depth for the sourcereceiver configuration as shown on Fig. 4.1.
and its 95% confidence zone when the source is shallower than 1200 m. It results in an
extended 95% confidence zone when the source depth is between 1200 and 3600 m. At larger
depths, there is no interaction of the PDF with the free surface and σZ remains stable at
approximately 590 m.
Generally, velocities increase with depth, which results in an increase of the location error with depth (or likewise a decrease of the error at shallow depth). This additional effect
is not reflected in Fig. 4.4) because propagation velocities were kept stable.

4.3

Source-receiver distance

In this section we numerically illustrate the influence of source-receiver distance on location
uncertainty. We use a setup of five receivers that surround a source at 3 km depth. We use
the travel-time uncertainties of σp = 0.0893 and σs = 0.170 s as estimated in Section 2.6.2.
We consider two different configurations.
Fig. 4.5(a) shows the setup for the first experiment. All five receivers are placed with a
source-receiver distance of 0.4 km and the location uncertainty is computed. Successively,
the source-receiver distances are increased to 20 km, with steps of 0.4 km, and for each configuration the location uncertainty is computed. Fig. 4.6 shows the resulting variation of
location uncertainty as function of array aperture.
When the aperture becomes infinitely small, it would be impossible to estimate a location based on travel-time differences. In this case, the PDF would be a half-sphere with
many epicenter-depth combinations with equal likelihood. Hence, the numerical illustration
starts for source-receiver distances of 0.4 km.
Fig. 4.6 shows that σ1 and σ2 decrease from 0.4 until 1.2 km distance and remain stable

60

4.3. SOURCE-RECEIVER DISTANCE

KNMI

Figure 4.5: Configurations for a numerical experiments in which (a) the aperture of a local
network is varied and (b) only the source-receiver distance for a single receiver is varied.
For each experiment, the same source location (red dot) is used, whereas positions of the
five receivers (triangles) are varied in experiment (a) and only the position of one receiver
is varied in experiment (b).
for larger distances. The uncertainties are relatively high when the local network has an
aperture that is considerably smaller than the source depth. From 1.2 km distance, the
epicentral location precision is not affected by an increase in aperture. The epicentral uncertainty is primarily dictated by the azimuthal coverage of receivers. Since this coverage
remains identical, also the uncertainty remains identical.
In practise, the epicentral uncertainty would grow with larger aperture as the travel-time
uncertainty does grow with epicentral distance (Section 2.6.1). This effect is not included
in this section. Here we merely model the geometric influence of the source-receiver configuration.
The depth uncertainty σZ behaves differently. For small array aperture it approximates
600 m. For larger array aperture it increases first with a large derivative, and for distances
beyond ∼3.5 km with a smaller derivative. As with the previous experiments, the trajectory
with smaller growth in σZ with aperture occurs when the PDF is bounded on the upper
side by the free surface.
Fig. 4.5(b) shows the setup for the second experiment. Four receivers are placed with a
source-receiver distance of 10 km. A fifth receiver is placed at a distance of only 0.4 km and
the location precision is computed. Successively, the source-receiver distances is increased
to 20 km, with steps of 0.4 km, whereas the other four receivers remain at 10 km distance.
Fig. 4.7 shows the resulting variation of location uncertainty.
Fig. 4.7 shows that σ1 and σ2 are approximately stable as function of source-receiver distance. As with the previous experiment, the azimuthal coverage remains identical and hence
there is no variation in epicentral uncertainty. Since four of the receivers are at 10 km distance, the case with a network aperture smaller than the source depth, does not occur.
The depth uncertainty σZ in Fig. 4.7 behaves similarly as σZ in Fig. 4.6 for the first
10 km source-receiver distance. This means that placing only a single receiver close to the
source zone is nearly as effective as placing all receivers near the source, in reducing the
depth uncertainty.
The depth uncertainty on Fig. 4.7 reaches a maximum when all receivers have a distance
of 10 km from the source. The depth uncertainty decreases again when the single receiver
61

4.3. SOURCE-RECEIVER DISTANCE

KNMI

Figure 4.6: Location uncertainty (σ1 , σ2 and σZ ) as function of network aperture. The
corresponding source-receiver configuration is shown in Fig. 4.5(a) for a source at 3 km
depth.

Figure 4.7: Location uncertainty (σ1 , σ2 and σZ ) as function of source-receiver distance
variation of one (out of five) receivers. The other four receivers are fixed at a source-receiver
distance of 10 km. The corresponding source-receiver configuration is shown in Fig. 4.5(b)
for a source at 3 km depth.

62

4.4. NUMBER OF RECEIVERS

KNMI

Figure 4.8: Location uncertainty (σ1 , σ2 and σZ ) as function of number of receivers surrounding a 3 km deep source at 5 km epicentral distance.
is placed at larger distances. This shows that it is favourable to have receivers at varying distances from the source, to better sample the P-S-delay-versus-distance function and
therewith to have a better estimate of source depth.
With a point-symmetric configuration, the 95% epicentral confidence zone is a circle, meaning that the ellipse minor axis is equal to the major axis, or σ1 =σ2 . In order to show
non-overlapping lines for σ1 and σ2 we have slightly perturbed the configurations in Fig.
4.5, resulting in σ1 and σ2 functions that can just be distinguished from each other in Figs.
4.6&4.7.

4.4

Number of receivers

In this section we numerically illustrate the influence of the amount of receivers on location
uncertainty. We use the low travel-time uncertainties of σp = 0.0893 and σs = 0.170 s as
estimated in Section 2.6.2. We use a setup with a source at the origin, at 3 km depth,
which is surrounded with receivers at a distance of 5 km. The receivers are placed at equal
azimuths. The minimum amount of receivers is three, in which case the azimuthal difference
(with respect to the source) between any two receivers is 120o . For the next case with four
receivers, the azimuthal difference is 90o , etc. The location uncertainty is continued until
the case with 40 surrounding receivers. The results are shown in Fig. 4.8.
Fig. 4.8 shows that both the epicentral (σ1 , σ2 ) and depth uncertainty σZ keep reducing with
an increasing amount of receivers. However, the returns are diminishing as the reduction in
uncertainty becomes smaller at larger receivers numbers. In this example, there is a limited
reduction in epicentral uncertainty beyond ∼10 receivers and a limited reduction in depth
uncertainty beyond ∼20 receivers.
For this experiment σ1 = σ2 , because point-symmetric configurations of receivers are chosen.

63

4.5. AZIMUTHAL GAP

KNMI

In reality, receiver configurations would be more chaotic, resulting in larger azimuthal gaps
(next section) and, as a consequence, σ1 > σ2 . The uncertainties in Fig. 4.8 are a lower
bound of what can be achieved.

4.5

Azimuthal gap

In this section we numerically illustrate the influence of azimuthal gap on location uncertainty. The azimuthal gap is defined as the largest angle between two stations as seen from
the epicenter (Scudero et al., 2021). Thus, it is the largest azimuthal difference between
any two stations. We use the travel-time uncertainties of σp = 0.0893 and σs = 0.170 s as
estimated in Section 2.6.2. We use a setup with a source at the origin, at 3 km depth and
receivers at 5 km distance.
Fig. 4.9 shows the receiver configurations and modelling results for three different azimuthal
gaps. Fig. 4.10 shows the location uncertainties as function of all considered azimuthal gaps.
The epicentral uncertainties (σ1 , σ2 ) increase with azimuthal gap. The depth uncertainty
σZ remains stable as function of azimuthal gap. σZ is primarily a function of the nearest
source-receiver distance (Fig. 4.7). Because this distance remains stable, also the uncertainty remains stable.
Fig. 4.9 shows in the left column the 3D location uncertainty 95% confidence regions.
For the upper two configurations (azimuthal gaps of 140 and 220o , respectively) the 3D
PDF can still be approximated as multivariate normal distributions. That is, the 95% confidence regions are still by approximation ellipsoids which size can be approximated with
σ1 , σ2 and σZ . For the lower configuration (azimuthal gap of 300o ) the 95% confidence
region does not look like an ellipsoid. It is a twisted volume in 3D which cannot be easily
parameterized. With these complicated PDFs, location may be ambiguous. By inspection
of the PDFs for all azimuthal gaps, we find that complicated PDFs start to occur from an
azimuthal gap of about 250o . This corresponds to the point Fig. 4.10 where σ1 seizes to
grow further. We mark this as the azimuthal gap beyond which the assumptions made with
the PDF parameterization are increasingly unvalid, and location may be problematic.

64

4.5. AZIMUTHAL GAP

KNMI

Figure 4.9: 95% location confidence regions for a setup with 5 receivers (green triangles)
on the Earth’s surface (yellow surface) in 3D (left column), in the epicentral plane (middle
column) and in the depth direction (right column). The top row shows results when the
receivers are placed such that there is an azimuthal gap ψ of 140o . The middle row shows
results for an azimuthal gap of 220o . The lower row shows results for an azimuthal gap of
300o .

65

4.5. AZIMUTHAL GAP

KNMI

Figure 4.10: Location uncertainty (σ1 , σ2 and σZ ) as function of azimuthal gap, for 5
receivers at 5 km epicentral distance from a 3 km deep source. The grey zone shows the area
where the parameterization of the location PDF with σ1 , σ2 and σZ is increasingly unvalid.
In this zone, the found uncertainties are likely underestimates.

66

5 Location-uncertainty maps
In this chapter, the location uncertainty maps are presented. These maps have been computed for the Netherlands and border regions, using the station distribution as shown in
Fig. 2.2.
The computation is done using a grid of (scenario) events with a grid-cell spacing of 1 km.
A local rectangular grid (Rijksdriehoekstelsel, RD) is used, which has as origin a location
close to Paris. This ensures that the whole of the Netherlands is in the first quadrant. The
computation grid spans RDx=[-20 300] and RDy=[270 670] km. (RDx,RDy)=(-20, 270) is
the grid center of the grid cell furthest to the southwest. In wgs84, this point corresponds
to (lat, lon)=(50.3940, 2.9262). (RDx,RDy)=(300, 670) is the grid center of the grid cell
furthest to the northeast. In wgs84, this point corresponds to (lat, lon)=(53.9951, 7.5980).
In this chapter, uncertainty maps are shown for parts of the Netherlands. In the appendix,
one-page maps are included that show location uncertainty for the whole of the Netherlands.
Even when these maps are printed on a single page, not all details are clearly visible. Therefore, also the underlying data with the computed values for each grid cell are distributed
via KNMI (2023b).
The azimuthal gap ψ is often used as a proxy for epicentral uncertainty. In this research we
explicitly compute epicentral uncertainty, so there is no need to show a proxy. Nevertheless,
we also derive ψ maps in this chapter. The reason to include ψ maps is that we found
that the used parametrization of location uncertainty ceases to hold for ψ ≥ 250o (Section
4.5). The maps in Appendices B.3&D.3 serve to indicate where this may occur. Moreover,
ψ could be used to assess the possibility to implement a well-constrained moment-tensor
inversion, which is possible when multiple quadrants of the radiation pattern are sampled.
At least two different quadrants are sampled when ψ ≤ 120o .
This chapter starts with a description of the detailed settings that have been chosen in
the implementation. The province of North Holland has been used for quality checking. In
the following two sections, the epicentral- and depth-uncertainty maps are presented. In the
last section, the updated magnitude-of-completeness map is published.
Examples in this report have been made with a station distribution that was in operation
in September 2021 (Table A.1). By September 2022 there have been several changes (Table
C.1). These changes are described in the magnitude-of-completeness section, which contains
both the 2021 and 2022 map. The 2021 and 2022 versions of MoC maps can be found at
KNMI (2023a) whereas the location-uncertainty maps are published at KNMI (2023b) and
can be found in the appendix. Yearly updates of the maps are foreseen. Overall, the station
coverage has been improved from 2021 to 2022, resulting in reductions in uncertainty.

67

5.1. SETTINGS AND QUALITY CHECKS

5.1

Settings and quality checks

5.1.1

Travel-time uncertainty

KNMI

When a distance-independent uncertainty can be assumed, the more picks are included, the
better (Fig. 4.8). There are diminishing returns, however. Beyond 20 picks there is hardly
an improvement. We restrict the amount of picks that is used to the 40 nearest detections
to limit computation time.
In Section 2.6 we found the travel-time uncertainty for different cases. When the azimuthal
coverage is good (ψ < 120o ) and there are many picks nearby (at least 12 at distances
smaller than 40 km) an event-specific travel-time model can be derived from these picks.
This yields the smallest uncertainties as estimated in Section 2.6.2.
For all other cases, a background travel-time model is used, which is an average over different regions in the Netherlands. This leads to larger uncertainties as estimated in Section
2.6.1. We found that the travel-time uncertainty is distance dependent. The uncertainty
was estimated for three different distance bins (local, regional and national) and listed in
Table 2.4. The effective uncertainty is implemented as a weighted average of picks available
in different distance bins. For P-waves we use
σp = (0.115n1 + 0.162n2 + 0.295n3 )/nT ,

(5.1)

where σp is the effective P-wave travel-time uncertainty for the collection of picks, n1 , n2
and n3 are the number of picks that fall within the first, second and third distance bin,
respectively, and nT is the total number of picks. Similarly, for effective S-wave travel-time
uncertainty we use
σs = (0.186n1 + 0.322n2 + 0.568n3 )/nT .
(5.2)

Eqs. 5.1&5.2 entail that the effective travel-time uncertainty increases when more distant
picks are used in the location assessment. In these cases, a more constrained source location
might be achieved when leaving out more distant picks. Whether it is favourable to include
more distant picks depends on the station coverage that is already achieved with local or
regional picks.
We test one approach of removing more distant picks. We assume that a good coverage
is achieved when there are at least 9 picks available in smaller distance bins. Thus, when
at least 9 picks are available in the local distance bin, picks are removed from the regional
and national bin. Likewise, if at least 9 picks are available in the local and regional bin
combined, picks are removed from the national bin. If say 20 picks are available in the local
distance bin, they are all used, because the larger amount of picks (with the same intrinsic
uncertainty) only helps to constrain the source location.
Fig. 5.1 shows a case in which the pick-restriction approach is favourable. A source is
chosen in the province North Holland ([RDx, RDy, depth]=[110, 520, 3] km) with M=2.0.
Picks can be made at 12 stations, 8 of which are local, 2 regional and 2 national (Fig. 5.1a).
The last two picks are at distant low-noise stations: T064 (Twente) and TERZ.01 (South
Limburg). Using all the available picks leads to the epicentral uncertainty ellipse as shown
in Fig. 5.1(b) in green. After removing the two most distant picks, a smaller epicentral
location uncertainty is obtained (purple ellipse in Fig. 5.1b).
Also for both cases, the depth uncertainty is computed. When all picks are used, σZ =472 m.
When the two most distant picks are omitted, σZ =446 m.
68

5.1. SETTINGS AND QUALITY CHECKS

KNMI

Figure 5.1: (a) Stations (orange triangles) where likely picks can be made for a M=2.0 source
in North Holland (black dot) at 3 km depth. (b) Size of the 95% confidence area when all
stations indicated in (a) are used (green ellipse) or when the two most-distant stations are
omitted (purple ellipse). This figure is for illustration purposes only. It was made with an
P 90
values.
obsolete version of the vrms

The above pick restriction is tested for all locations in the Netherlands. Fig. 5.2 shows
σ1 results for a part of North Holland, when (a) the pick restriction is implemented and (b)
when the pick restriction is not implemented. It can be seen that small gains are indeed
achieved at some place. That is, Fig. 5.2(a) contains areas with colder colours than Fig.
5.2(b). On the other hand, Fig. 5.2(a) also contains areas where the epicentral location
deteriorates when the pick restriction is implemented, leading to warmer colours than in
Fig. 5.2b). Overall, the effect of the picking-restriction approach (i.e., difference between
Fig. 5.2(a) and Fig. 5.2(b)) is small.
For implementing an ideal pick restriction, for each source location a grid search would need
to be done on many possible omissions of more distant picks. This additional grid search
could lead to a hundred-fold computation increase, while bringing only modest gains (Fig.
5.1). For that reason, a more refined pick-restriction approach is omitted and simply all 40
closest picks (when available) are included for location.

5.1.2

Grid-search settings

The epicentral search is done within a square that is defined around the event location. For
magnitudes 0.5 until 1.5, the square is defined with parameter a = 2(11 − 9/27p), where p is
the amount of P-wave detections for the event. The square has sides with size a in km and
the actual event location is always in the middle of the square. One side is at least 2 km.
Hence, in the limiting case when location is done with only 3 detections, the size of search
grid is 20x20 km. With 30 or more detection, the area is reduced to 2x2 km. Irrespective
of the size of the area, the zone is sampled with 10000 equally spaced points to obtain the
epicentral PDF.
For magnitudes 2.0 until 4.0, the square is defined with parameter a = 2(31 − p). For
69

5.1. SETTINGS AND QUALITY CHECKS

KNMI

Figure 5.2: σ1 location uncertainty in North Holland and surroundings when (a) removing
picks from large-distance bins if nine local or regional picks are available and when (b) leaving
out the most distant picks if more than 40 picks are available. The epicentral uncertainty
has been computed for M=2.0. Thick black lines denote the coastline. Dark red depicts
uncertainties of 2800 m and higher. This figure is for illustration purposes only. It was made
P 90
with an obsolete version of the vrms
values.
70

5.1. SETTINGS AND QUALITY CHECKS

KNMI

Figure 5.3: Grid search within (a) 60 km rectangular search grid (as depicted with the red
square) and (b) 30 km search grid, yielding σ1 =8552 and σ1 =1516 m, respectively. Note that
in (a) the 95% epicentral confidence zone includes a small area in the southeastern edge of
the grid-search area, which is responsible for the large extend of the fitted σ1 value Examples
are for M=2.0 with [rdx, rdy]=[112, 558] km. This figure is for illustration purposes only. It
P 90
values.
was made with an obsolete version of the vrms
these magnitudes, more extreme cases occur with detections at only a few distant stations.
This gives the need to extend the grid-search area, to cover potentially very large PDFs.
For these magnitudes, when there are only 3 detections, the size of the grid-search area is
56x56 km. When there are 30 or more detections, the area is reduced to 2x2 km.
The above settings have been tested for many different sensor configurations, to make sure
that the grid search contains the area where nearly all location chances lie. In the testing,
one case was found where these settings may not be desired. This is a case with a nearlinear array of receivers as shown in Fig. 5.3(a). In this example there is an ambiguity. A
large part of the probability lies around the actual source location, but some part of the
95% confidence zone lies in a separated blob that is a blurred mirror image (with the line
extending from the receivers acting as the mirror) from the actual event location. Through
these separated areas, an ellipse is fitted, which results in an ellipse with a large major axis
P 90
and hence a large apparent σ1 value. When using an outdated vrms
list, this occured for
M=2.0 in a zone at the southern side of Texel (Fig. 5.2). In these cases that the confidence
area is split up in two zones that are well separated, it would be possible to find the actual
source area by using the polarization of the incoming P waves. In order to show only the
parameters of the main blob, for this case the grid-search area could be reduced to a square
of 30x30 km yielding the results as shown in Fig. 5.3(b).
To keep a fixed procedure, we have not implemented the above fix. This means that for
some areas with very large modeled σ1 values, the uncertainty could be reduced in reality
by including a polarization data attribute.
The discretization settings for computing the depth PDF are as follows. A search in depth
is done from 0 to 20 km with steps of 50 m. The actual source is placed at 3 km depth. As
such, the PDF is well sampled for any receiver-distribution scenario. As only 401 points are
evaluated in the computation, no further optimizations are needed to save computational
time.

71

5.1. SETTINGS AND QUALITY CHECKS

KNMI

Figure 5.4: Large difference in uncertainty of nearby source locations (a) [rdx, rdy]=[155,610]
and (b) [rdx,rdy]=[155,618] because of large difference in receiver (orange triangles) azimuthal coverage. σ1 is (a) 1132 and (b) 11314 m, respectively. Examples are for a M=2.0
at 3 km depth. This figure is for illustration purposes only. It was made with an obsolete
P 90
version of the vrms
values.

5.1.3

Discontinuity in uncertainty

When generating test maps, cases were found where the uncertainty varies largely from one
to the other grid point. These discontinuities were seen as possible implementation errors.
However, after inspection, it turned out that these discontinuities would actually occur. One
of the cases is shown in this section.
In Fig. 5.4 an example is shown of σ1 increasing about ten fold by moving the scenario
source location just a few kilometers. In Fig. 5.4(a) the uncertainty is relatively small with
5 receivers picking up the signal from this event. When the event is moved a few kilometers
to the north (Fig. 5.4b) only 3 receivers, in a more or less linear configuration, pick up
the event. Like in Fig. 5.3(a), the probabilities lie in two mirrored confidence zones, but
in this case these two zones have merged, which makes it harder to disentangle them. This
merging of two confidence zones occurs when the source location lies in the extension of a
linear configuration of receivers and the actual source location is also close to the line that
can be fitted through the receivers. In other words, these extended confidence zones occur
if the azimuthal gap is very large.

72

5.2. EPICENTRAL UNCERTAINTY

5.2

KNMI

Epicentral uncertainty

In this section, epicentral uncertainty is computed for hypocentres at 3 km depth. For this
depth, the epicentral probability density function (PDF) is computed. This PDF is approximately a bi-variate normal distribution, of which a confidence zone typically looks like an
ellipse. This PDF is parameterized with σ1 , σ2 and θ. σ1 and σ2 describe the uncertainty in
the major-axis and minor-axis directions of the ellipse, respectively (Fig 3.1c). θ describes
the orientation of the ellipse. It is defined as the smallest clockwise angle between north
and the major axis.
In this section we show the computed values for specific areas in the Netherlands for an
event depth of 3 km. The effect of event depth on epicenter uncertainty has been studied in
Section 4.2. σ1 -maps for the Netherlands can be found in Appendices B.1&D.1 .
Fig. 5.5 introduces the plotting style that has been chosen. Fig. 5.5(a) shows the grid
cells for an area around North Holland. The colour of the grid cell denotes the modeled
σ1 value for a scenario M=1.5 event in the centre of the grid cell. A white grid cell means
that for the center point there is an insufficient number of detections (< 3) and hence no
location, nor location uncertainty, can be computed.
In Fig. 5.5(b) the visualisation is simplified by removing the grid lines and discretizing σ1
with 200 m levels. Hence, the lowermost uncertainty level, denoted by dark blue, contains all
values from 0 to 200 m. Higher uncertainties are shown with warmer colours. Dark orange
denotes uncertainties between 2600 and 2800 m. Red denotes values of 2800 m and higher.
In practise, there are little places where σ1 > 3000 is reached.
In Fig. 5.5(c) the visualisation is made more eye-pleasing by interpolating the borders between different uncertainty levels. This is the plotting style that is chosen in this report.
The exact place where borders are drawn depends on the interpolation scheme that is used.
This post-processing leads to small perturbations. Thus, this map gives a good overview of
the uncertainties, but for the actual value at a position, one needs to resort to the value of
the underlying grid cell.
Fig. 5.6 gives an overview of the different epicentral uncertainty parameters. Again the
same part of North Holland is shown as in Fig. 5.5 for scenario events with M=1.5. Fig.
5.6(a)&(b) show the σ1 and σ2 distribution, respectively. For the entire area where a location can be made, σ2 remains small to modest (< 1200 m). Outside of the area covered by
the local network, however, there are areas where σ1 becomes large (≈ 3000 m). Fig. 5.6(c)
shows the aspect ratio σ1 /σ2 . In the same areas where σ1 is large, also the aspect ratio is
significant. This indicates that the uncertainty ellipse is elongated. Also Fig. 5.6(d) shows
the size of the aspect ratio, by the size of the stripes. At the same time, the orientation θ of
the uncertainty ellipse is shown through the orientation of the stripe. In an area surrounded
by geophones, the aspect ratio is very small, resulting in points instead of stripes. At some
places, the orientation of the ellipse, and the aspect ratio, vary smoothly as function of
distance from a station. At places where there is a discontinuity in detectability (addition
or removal of stations that can be used for location) there is also a discontinuity in aspect
ratio and/or ellipse orientation. Altogether, the interplay of irregularly placed stations with
different detection capabilities (noise levels) results in a complicated spatial pattern of epicentral uncertainty.
From the three parameters describing epicentral uncertainty, one would generally be interested in the largest uncertainty σ1 . Overview maps of σ1 for the complete zone of investigation (the Netherlands and border areas) can be found in Appendices B.1&D.1. They have
been computed for M=0.5 to M=4.0 with steps of 0.5. For each magnitude and each grid
point, the values of θ, σ1 , σ2 and ψ are available at KNMI (2023b).

73

5.2. EPICENTRAL UNCERTAINTY

KNMI

Figure 5.5: Comparison of plotting styles for showing the σ1 distribution over North Holland
for scenario events with M=1.5; (a) With color gradient and showing the underlying grid, (b)
with one colour for each 200 m bin and (c) with interpolated contours between the different
σ1 zones. The coast is depicted with a thick black line, province boundaries with a thin
black line, surface sensors are shown with squares and deep geophones with inverted orange
triangles. Coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

74

5.2. EPICENTRAL UNCERTAINTY

KNMI

Figure 5.6: Epicentral uncertainty plots for North Holland, for scenario events with M=1.5.
(a) Uncertainty in direction with largest uncertainty σ1 , (b) uncertainty in direction with
smallest uncertainty σ2 , (c) aspect ratio between the two orthogonal uncertainty directions
in the horizontal plane: σ1 /σ2 . (d) a quiver plot showing the aspect ratio (through size of
the stripes) and orientation of the uncertainty ellipses (through orientation of the stripes),
for a 1x1 km grid of scenario events. Coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

75

5.2. EPICENTRAL UNCERTAINTY

KNMI

When inspecting the σ1 figures for the different magnitudes (Fig. B.1–B.8) one can see
that, first of all, the area where a location can be made grows with magnitude. Also,
generally, the uncertainty reduces with magnitude. This is to be expected, since higher
magnitudes lead to higher amplitudes, more detections and thus a better receiver coverage.
Within the onshore part of the Netherlands, the location uncertainty keeps reducing until
M=3.0. For larger magnitudes, the location accuracy remains practically stable.
Higher magnitudes lead to inclusion of picks from further distances, which results in an
increase in the effective travel-time uncertainty. At some places, the increased station coverage is undone by the increased travel-time uncertainty, leading to a (small) increase in
location uncertainty with magnitude. These places can, e.g. be noted by comparing Fig.
B.7 and Fig. B.8. In Section 5.1.1 we discussed this effect. It would be possible to undo the
(small) increase in σ1 by finding an optimal combination of picks to include for the location.
Because the gain is small, and the additional computational expense large, this fix has not
been implemented.

76

5.3. DEPTH UNCERTAINTY

5.3

KNMI

Depth uncertainty

In this section it is shown how the expected depth uncertainty varies over the Netherlands.
The depth PDF is computed at the epicentral location for sources at 3 km depth. The last
panel in Fig. 1.1 shows an example. The depth PDF is parameterized with σz which is approximately one-fourth of the depth interval that contains 95% of the location probabilities.
σz -maps for the Netherlands can be found in Appendices B.2&D.2. The underlying values can be obtained at KNMI (2023b). When inspecting the σz figures for the different
magnitudes (Figs. B.9–B.16) one can see that, first of all, the area where a location can
be made grows with magnitude. Also, generally, the uncertainty reduces with magnitude.
However, the uncertainty reduction with magnitude is more limited than for σ1 , because
σz is especially dependent on the sampling of the travel-time curve at near offsets, which
sampling does not further improve for higher magnitudes. The σz maps for M=3.5 and
M=4.0 are nearly identical.
The general pattern of σz can be studied on, e.g., Fig. B.12. The further an event location is from the nearest station, the larger the depth uncertainty. The smallest location
uncertainties (< 500 m) are reached in areas with the densest station distributions. This
occurs for the dense Groningen network with receiver spacings of around 4 km, but also
for dense accelerometer distributions in Friesland and Drenthe. For the area covered by
the networks in Twente, South-Holland and North-Holland there is a larger depth uncertainty (500 < σz ≤ 1000) than in Groningen. In the former areas there would, overall, be
less travel-time picks available to constrain the location than in Groningen and surroundings.
Fig. 5.7 shows M=1.5 σz maps for two areas in the Netherlands. σz is shown in more
detail than in the appendix; in Fig. 5.7 the σz values are discretized with steps of 200 m, in
Fig. B.11 with steps of 500 m. In Groningen (Fig. 5.7a), the quite uniform receiver layout
results in a quite uniform σz < 400 m over most of the network. In the southeast of the area,
there is an additional infill network to monitor gas storage in salt caverns: the Zuidwending
network. In the area covered with this very dense network, σz < 200 m is achieved.
In Fig. 5.7(b) the station distribution is more patchy. Near two geophones in Twente
σz < 600 m is achieved. The smallest uncertainties (σz < 400 m) are obtained here near the
accelerometers in South Drenthe. Because these accelerometers are still relatively close to
the G-network, one obtains more travel-time picks for event locations in this area, than for
locations in Twente.

77

5.4. MAGNITUDE OF COMPLETENESS

KNMI

Figure 5.7: M=1.5 depth-uncertainty σZ distribution for two areas in the Netherlands: (a)
Groningen and (b) an area covering the south of Drenthe and Twente. Accelerometers are
shown with purple squares, deep geophones with orange triangles. Uncertainties of 1400 m
and higher are shown in red.

5.4

Magnitude of completeness

In this section we use the tools developed in Chapter 2 to update the magnitude-of-completeness
(MoC) map for the Netherlands. This new map supersedes the map published in Kruiver
et al. (2021). The MoC at a certain point is the minimum magnitude for which an earthquake can be located in severe noise conditions. This corresponds to a minimum of three
detections that are likely to occur in the presence of P90 noise conditions. The map is
computed for a typical induced-seismicity focal depth of 3 km. The map is made using a
1x1 km grid. A contoured version of this map is published in this report. The computed
values for each grid point can be found at KNMI (2023a).
Computation of a MoC (for one point) entails the following steps:
• Compute the expected amplitude at all receivers in the detection list, using the P-wave
GMPEs (Section 2.2).
• Repeat this computation for a range of magnitudes.
• Implement the picking threshold (Section 2.4).
• For hard-rock sites, implement a hard-rock correction factor (Section 2.5).
• Compare the resulting amplitudes with the station P90 noise levels (Section 2.3) to
find the minimum detection magnitudes at each receiver.
• Order this detection-magnitude-per-receiver list from low to high. The third-lowest
magnitude is the MoC for the evaluated source location.
• MoC values lower than 0.4 are replaced with 0.4. Clipping at MoC=0.4 is done to
make sure that the modeled outcome falls within the range for which the underlying
P-wave GMPE has been calibrated.
The term magnitude of completeness implicates that the catalogue is complete from the
stated MoC. With other words, if the MoC for a certain region is 0.5, it implicates that all
events with M=0.5 and higher are detectable and locatable. However, one needs to take into

78

5.4. MAGNITUDE OF COMPLETENESS

KNMI

account that also the modeled MoC is a probabilistic value. It is computed using the 90th
percentile of the noise distribution, meaning that it yields MoC values that can very often,
but not always, be reached in practise. If the P90 MoC value is 0.5, then 90 percent of the
time a M=0.5 event would be locatable. However, there remains a 10 percent chance (e.g.,
during stormy days) that noise levels are so high that a M=0.5 event would be missed.

5.4.1

Sensitivity tests

Fig. 5.8 shows the impact of noise-percentile choice on the estimation of MoC. Using the
P10 from the noise, a very optimistic magnitude distribution is obtained as shown on Fig.
5.8(a). When noise is very low (during quiet summer nights) a magnitude of 0.8 (cyan
contour) can be detected and located over a large part of North Holland. For average noise
conditions (P50, Fig. 5.8b) the area where a magnitude of 0.8 can be located is significantly
smaller. When the 90th percentile of the noise is used (P90, Fig. 5.8c) an 0.8 event can only
be located in an area around Alkmaar. This last map shows the magnitude distribution
we refer to as MoC. For the zone indicated with MoC=0.8 on this map, a nearly complete
catalogue could be reached for M≥0.8. So, it is the P90 level we further use when making
MoC maps.
Another point of attention is the choice of three detections. The number of three makes
sense from a location perspective, since it is the minimum number for which a 3D location of
the source can be found. It is not the amount of receivers, however, that is currently used in
an operational sense at KNMI. At KNMI the software package SeisComP (GEMPA, 2022)
is used, which only starts up a location algorithm after detection at six stations. Fig. 5.9
shows, for the Groningen region, the MoC that is obtained when requiring (a) detection at
six receivers and (b) at three receivers. Clearly, lower MoC levels can be reached when only
three detections are needed. Especially for monitoring the gas-storage plants it was deemed
important to reach these lower MoC levels also in an operation sense. This was achieved by
placing also the second-deepest geophones in the detection list, for the geophone strings at
the Norg and Grijpskerk networks. For the Groningen gas field, Fig. 5.9(b) shows the MoC
that could be reached by post-processing the automated picks.

5.4.2

2021 Map

Fig. 5.10 shows the MoC map for the Netherlands, with a contouring of 0.5 magnitude levels. At regions where dense distributions of borehole geophones exist, MoC≤0.5 is reached
(dark-blue areas). An exception is for the network in South Holland, primarily due to high
noise levels (Fig. 2.14). In the south of Limburg, MoC≤0.5 is reached with broadband stations, with some of them having very low noise conditions. Around the regional networks,
there are zones in which 0.5<MoC≤1.0 is reached (blue), and at further distance from the
dense station distributions, a zone in which 1.0<MoC≤1.5 is reached (light blue). In areas
far from station clusters, like most of Utrecht and Gelderland, 1.5<MoC≤2.0 (yellow). In
the onshore part of the Netherlands MoC≤2.0 with the exception of a tiny part of Zeeland.
Most of the offshore region on the map has 2.0<MoC≤2.5 (orange).
Note that the map does show parts of Germany and Belgium, but would not be correct
for these surrounding countries. The seismological authorities across the border use additional stations of their own networks. Also take into consideration that we did not include
British stations (nor Danish and Norwegian stations) in this study. Hence, further offshore
the MoC could be lower than modeled in this study.
Within the dense borehole network in Groningen, a MoC of 0.4 is achieved (Figs. 5.9

79

5.4. MAGNITUDE OF COMPLETENESS

KNMI

Figure 5.8: Magnitude-of-completeness (MoC) estimation at the northwestern part of the
Netherlands, using (a) the 10th percentile, (b) the 50th percentile and (c) the 90th percentile
of the noise distribution (Section 2.3). The actual MoC is obtained using the P90 noise
distribution. Contours are shown for MoC values of 0.4 (blue), 0.8 (cyan), 1.2 (light green),
1.6 (yellow) and 2.0 (orange). Surface sensors are shown with purple squares, deep geophones
with orange triangles. Background map is from www.openstreetmap.org.

80

5.4. MAGNITUDE OF COMPLETENESS

KNMI

Figure 5.9: Magnitude of completeness (MoC) for the Groningen area using a minimum of
(a) six detections and (b) three detections. Contours are shown for MoC values of 0.4 (
blue), 0.8 (cyan) and 1.2 (light green). Background map is from www.openstreetmap.org.
and 5.10). When MoC is computed with historical events from the KNMI induced event
catalogue, a higher MoC would be obtained. There are different reasons for this:
• The modeled MoC would always be slightly lower since it is based on P90 noise levels.
This means that events with magnitudes close to MoC would just be missed in severe
noise conditions.
• For determination of a local magnitude, the records are first convolved with a WoodAnderson response function. This convolution acts as a high-pass filter for frequencies
larger than 1.0 (Chung and Bernreuter , 1981). In the Groningen area there is severe
low-frequency noise (0.4 to 2.0 Hz) due to the nearness of the Wadden Sea and North
Sea. As a result, a part of the low-magnitude events has sufficient SNR to be detected
and located using the 5-to-40 Hz frequency band, but cannot be assigned a local magnitude due to insufficient SNR at frequencies between 0.9 and 5 Hz. Currently, an
analogue local-magnitude relation is derived that is especially suited for the smallmagnitude events. For the time being, the KNMI induced data catalogue misses some
of the locatable events for approximately M≤0.6.
• The modeled MoC is for the network layout of September 2021, whereas measured
MoC is based on events in the past in which the network was less well developed.
From the above list, it is expected that the second item is the most important; for the
Groningen region, over the years a long list was built up with events that were too small
to assign a local magnitude. For other regions, there are no such lists. It is hence expected
that for other regions, the modeled MoC would come quite close to an MoC as derived
from the published KNMI catalogue. For example, an M0.3 event in Twente —if it were to
happen— can always be assigned a magnitude since there is no strong noise source between
0.9 and 5 Hz. In South-Holland there is a similarly strong low-frequency noise source as
in Groningen (i.e., the active coast). However, in South-Holland the MoC≥0.8. An M0.8
would yield sufficiently strong low-frequency amplitudes to assign a local magnitude.

81

5.4. MAGNITUDE OF COMPLETENESS

KNMI

Figure 5.10: Magnitude of completeness over the Netherlands, computed using P90 noise
levels in 2021, detections at three stations, and the station list of September 2021 (Fig. 2.2).

82

5.4. MAGNITUDE OF COMPLETENESS

5.4.3

KNMI

2022 Map

Also for the network state of September 2022 the MoC is computed. For the 2022 map, data
from September 1 2021 until August 31 is used for computing the vrms values. The generally
small changes in noise conditions result in small changes in the MoC. Larger changes in MoC
are caused by chances in the effective network layout. The following changes exist from the
September 2021 to September 2022 detection list:
• Addition of 9 stations that had insufficient data over 2021, but sufficient data over
2022 (G024, G314, STR, SUH4, T014, WDB4, WMH4, WYN4, AHRW).
• Addition of 2 new borehole geophones (ZH074 and VDM24).
• Addition of 8 new accelerometers (DOK, KMP, OOTH, SNB, EVK, PMR, BRTL and
OLK).
• Removal of 1 accelerometer that was closed in 2022 (BHAR).
• Removal of geophone ZWE2 because of bad data quality.
• Replacement of borehole geophone CHA5 with surface broadband station NE116.
The locations of new stations can be found at, e.g., http://rdsa.knmi.nl/network/NL/.
Fig. 5.11 shows the MoC for 2022. Almost everywhere the MoC is either stable or improves with respect to the 2021 version (Fig. 5.10). Only in north of Friesland, the MoC
somewhat deteriorates due to the removal of geophone ZWE2 from the detection list. It can
be seen that both the MoC=1.5 area (cyan) and MoC=1.0 area (blue) become somewhat
smaller here.

83

5.4. MAGNITUDE OF COMPLETENESS

KNMI

Figure 5.11: Magnitude of completeness over the Netherlands, computed using P90 noise
levels in 2022, detections at three stations, and the station list of September 2022.

84

6 Discussion and Conclusions
In this report, location-uncertainty maps were derived for the Netherlands, using the network
configuration as existed in September 2021 (Fig. 2.2 and Table A.1) and September 2022
(Table C.1). The resulting maps can be found in Appendices B&D and the underlying values
can be obtained at KNMI (2023b). Yearly future updates are anticipated.
As a side-product, magnitude-of-completeness maps were updated. The 2021 and 2022
versions of these maps can be found in this report (Figs. 5.10&5.11) and the values for each
grid point can be obtained at KNMI (2023a).
In this chapter, the produced maps, and the underlying results from the previous chapters, are summarized and discussed.

6.1

Detection

For modelling whether a detection can likely be made, P-wave ground-motion prediction
equations (GMPEs) were derived (Section 2.2). These equations were specifically derived
for seismicity in the upper crust, for a frequency band between 5 and 40 Hz, which is used
for detecting and locating induced seismicity. They are empirical models that are constrained with a database of more than 20,000 PGV recordings from induced seismicity in
the Netherlands. The models have an empirical base between local magnitudes of 0.4 and
3.6, epicentral distances between 0 and 120 km and earthquake depths between 2 and 3 km
(or hypocentral distances between about 2 and 120 km). A model was derived for motions
at the free surface (SL model) and at around 200 m depth (DL model). Both models are
the same, except for the static-offset fitting term (c1 in Eq. 2.1). This term is 4.06 times
larger for the SL model than for the DL model (Table 2.1).
The derived GMPEs are specifically for sensors at soft-sediment settings. For callibration,
only sensors in or at unconsolidated sediments were used. Most sensors at depth are in a
P-wave velocity setting between 1.6 and 2.2 km/s (Hofman et al., 2017). Near the Earth’s
surface, P-wave velocities can be significantly lower (van Ginkel et al., 2022b). For sensors
installed at bedrock conditions, a correction factor was derived (Section 2.5) to account for
the higher seismic impedance and hence lower recorded amplitudes.
S-wave GMPEs were not used in this project. We assumed that, when a P-wave can be
detected and picked, also an S-wave can be picked and used for locating the event. As
nearly always S-waves are stronger than P-waves, this is mostly a valid assumption. Still, it
can be argued that this assumption is both optimistic and pessimistic. We think that the
net effect remains that it is a good assumption to make.
When a station happens to be in or close to a nodal point of the S-wave radiation pattern,
the S-wave amplitude would be smaller than the P-wave amplitude. So, in a few percent

85

6.2. LOCATION PDF

KNMI

of the cases it could be that we assume that an S-wave can also be picked, whereas this is
not possible. On the other hand, it can also happen that the station is in the nodal plane
for the direct P-wave. This situation is implicitly included in the derivation of the P-wave
GMPE. The implication is that, according to the model, both a P-wave and S-wave pick
cannot be made. In reality, it might still be possible to pick S and include it for the location.
An important choice we made when modelling detection, is the noise percentile. Throughout
the report, we assumed the 90th percentile, meaning a noise level that is only exceeded 10%
of the time. Therewith, the generated maps show location accuracy that can be achieved
in quite severe noise conditions. In low or average noise conditions, there would be more
detections, leading to a better location.
The somewhat conservative noise-condition choice is partly offset by another factor that
have not been included. It has been assumed that all stations in the detection list provide
data. In practise, about 10% of the stations has issues that result in no real-time data being
available at KNMI, and hence no (automatic) detection and pick being made.
In Section 2.3.3 we found that the average noise reduction with depth is a factor of 26,
when sensors are placed at 200 m instead of at the Earth’s surface. In Section 2.2.2 we
found that P-wave amplitudes at the Earth’s surface are, on average, a factor of 4 higher
than at 200 m depth. Hence, effectively a factor of 6.5 gain in detection is reached by placing
sensors at 200 m depth. This better detectability also resulted in a much larger empirical
base of the derived GMPE for motions at 200 m depth: 15,078 vs 7,592 recordings that
passed the (automated) QC criteria (Fig. 2.9). Note that these numbers are only valid for
the investigated picking frequency band (from 5 to 40 Hz).

6.2

Location PDF

We used a Bayesian framework (Section 3.1) to map uncertainties in travel time (Section 2.6)
to uncertainties in location. We considered location using differential P-wave arrival times
over a network (P-delays). Also, we considered location using the travel-time difference
between first P and S waves at each receiver (P-S delays). We found that a joint inversion
with both P- and P-S delays yields smallest location uncertainty (Section 3.3). Both data
attributes have complimentary information on the location problem. With sensors at or
close to the Earth’s surface, P-S delays are more sensitive to source depth, whereas P-delays
are more sensitive to source epicenter.
The location probability density function (PDF) can in many cases be approximated well
as a multivariate normal distribution (Section 3.2). We expressed the uncertainty in the
different directions as the standard deviations of the PDF in these directions. Most of the
probability lies within one standard deviation from the most likely source location (i.e., the
mode of the PDF). Hence, this standard deviation gives a good measure of the uncertainty
to take into account. However, bear in mind that the PDF describes a wide distribution of
locations and that also less likely locations remain possible.
The 3D PDF may be very complicated and many parameters are needed to describe its
shape. We have simplified its description to 4 parameters only (Fig. 3.1):
• σ1 standard deviation describing the largest uncertainty in the horizontal plane, at the
depth of the source.
• σ2 standard deviation describing the smallest uncertainty in the horizontal plane, at
the depth of the source.

86

6.2. LOCATION PDF

KNMI

• θ (smallest clockwise) angle of the σ1 direction with respect to north.
• σZ standard deviation describing the depth uncertainty, at the epicenter of the source.
σ1 and σ2 are computed at the most likely depth of the source, σZ is computed at the
most likely epicenter. Note that for an actual event, at least one iteration is needed to
obtain these parameters. In the KNMI operational flow, first a solution is found with the
Hypocenter method (Lienert et al., 1986) using an average northern Netherlands or southern
Netherlands (Lower Rhine Graben) velocity model. This first solution yields an approximate
depth and epicenter. Only in a second iteration, an optimized travel-time model is used and
location uncertainties are computed. If the second solution is quite different from the first,
another iteration is needed.
Depth uncertainty behaves differently from epicentral uncertainty. In general, the depth
is less well constrained than the epicenter (e.g., Fig. 4.7). For constraining source depth
very well, receivers would be needed both above and below the source, which situation does
not occur in the Netherlands seismic network. Because depth uncertainty is larger, also the
dynamic range for displaying σZ (Appendices B.2&D.2) has been chosen larger than for σ1
(Appendices B.1&D.1). For constraining the depth, it is most important to have receivers
nearby the epicenter, so that the curved part of the travel-time function is sampled well
(Section 4.3). For constraining epicenter it is more important to have a good azimuthal
distribution of receivers (Section 4.5).
In seismology, different institutes use different ways to communicate uncertainty. With
the assumption that a multivariate normal distribution is described, the most-common uncertainty measures can be translated to the standard deviations used in this report:
• The European-Mediterranean Seismological Centre (EMSC) uses the 95% confidence
zone when communicating epicentral uncertainty. For earthquakes in their catalogue
they list the semi-major and semi-minor axis of the 95% confidence ellipse. This
corresponds to 2.45 times the σ1 and σ2 values used in this report, respectively (Eq.
3.17).
• The EMSC also communicates depth uncertainty. However, it is not clear how this
parameter is defined. We assume that their depth-uncertainty parameter is analogous
to their epicentral-uncertainty parameters; half of the depth-location 95% confidence
interval. In that case, σZ needs to be multiplied by a factor of 1.96 to obtain their
uncertainty parameter.
• In the IASPEI Seismic Format (ISF; http://www.isc.ac.uk/standards/isf) also epicentral and depth uncertainty are treated separately. An epicentral ellipse is used that
describes the 90% confidence area. In ISF, the semi-major and semi-minor axis of that
ellipse are reported. In Wang et al. (2015) it is derived that for a bi-variate normal
distribution 90% of the probabilities lie within 2.1460 times the standard deviation
from the mean. Hence, a factor of 2.15 needs to be applied to translate the σ1 and σ2
values in this report to the ISF format.
• In the ISF format, the depth error is communicated as the extend of the 90% confidence
interval. Following Wang et al. (2015), for a (1D) normal distribution, 90% of the
probabilities lie within 1.6449 times the standard deviation from the mean. Hence, a
factor of 3.29 needs to be applied to translate σZ from this report to the ISF format.
We found that the simplified parameterization of the location PDF (in σ1 , σ2 , σZ and
θ) starts to break down when the azimutal gap is larger than 250o (Section 4.5). For this
reason we have also included maps displaying the azimuthal gaps for the various magnitudes
(Appendices B.3&D.3). For ψ ≥ 250o location uncertainty is likely underestimated with the
87

6.2. LOCATION PDF

KNMI

values shown on the maps in Appendices B.1&D.1 and B.2&D.2. When ψ ≥ 250o , it would
be important to find an additional constraint on the location, e.g., the polarisation of the
direct P-wave.
Generally, ψ (and σ1 , σ2 , σZ ) become smaller with increasing magnitude. The larger the
magnitude, the more stations detect the event, the better the azimuthal coverage and hence
the smaller the azimuthal gap. This is also reflected on the maps in Appendices B.3&D.3.
Note, however, that we built in a pick restriction. Only the 40 nearest picks are included
in the location. For that reason, at some places ψ (slightly) increases with magnitude. At
such places, ψ does remain quite low; with already 40 picks, chances are very slim that a
large ψ occurs.
The maps published in this report show contoured location-uncertainty levels. The levels at which these contours are drawn are somewhat subjective. Moreover, for drawing
contours, an interpolation is done between the computed values on a grid of 1x1 km. This
interpolation can result in small deviations from the actual computed values. Hence, if local
values are needed with high precision, the values at the underlying grid are to be retrieved
from KNMI (2023b).
The location-uncertainty maps were specifically derived for upper-crustal seismicity. For
calibration, events were used that were caused by depleting gas reservoirs. These events
have been found to be no pure slip events, but also to contain a negative isotropic component (Kühn et al., 2020; Dost et al., 2020). As a consequence, their P-wave content is
somewhat larger than for tectonic events. Given the large variability in ground motions
however (Table 2.2) also tectonic and other induced events are probably still quite well described with the derived GMPEs and hence the uncertainty maps can also be used for other
upper-crustal events.
We have computed the uncertainty maps for a typical induced-seismicity case with a source
depth at 3 km. In Section 4.2 we assessed what the effect would be if the source were to
occur at a different depth. For larger depths, the epicentral uncertainty would be larger,
whereas the depth uncertainty would be similar. For that modelling, however, we did not
include that we have less knowledge of the velocities at larger depths, and thus a larger
travel-time uncertainty and location uncertainty would result (Section 4.1). Combining the
insights from Sections 4.1 and 4.2 both the epicentral and depth uncertainty are expected to
grow for depths larger than 3 km. Likewise, the epicentral uncertainty would be somewhat
smaller when the source is at, e.g., 2 km depth instead of 3 km. The depth uncertainty
somewhat undulates at depths smaller than 3 km due to interaction of the PDF with the
free surface (Fig. 4.4).
In this report we made the case that P-S delays need to be included in a joint inversion
in order to improve the depth estimation of the source (Fig. 3.2). In the operational workflow for Groningen seismicity, however, location is done using only P-delays. To maintain an
induced-earthquake catalogue with a consistent workflow from the nineteen eighties onwards,
also for the Groningen region an average North-Netherlands (NN) 1D velocity model is used
(Kraaijpoel and Dost, 2013). The NN S-wave model is not a good average for the Groningen
region. The North-Sea Group (top 800 m) S-wave velocity in the NN model is 1.16 km/s
whereas average velocities in Groningen are 0.53 km/s (van Ginkel et al., 2020). Using the PS data attribute together with the NN model would result in placing sources further from the
network than they actually are. Thus, in the first iteration, only P-delays are used and the
depth is fixed to 3 km. This is the location solution that is published in the online catalogue
(https://www.knmi.nl/kennis-en-datacentrum/dataset/aardbevingscatalogus). In following
iterations, P-S delays are included together with better local velocity or travel-time models,

88

6.3. KNOWLEDGE OF TRAVEL-TIME UNCERTAINTY

KNMI

like local 1D profiles from the detailed 3D velocity model (Romijn, 2017). The operator
of the Groningen gas field (NAM) uses the latter 3D model to locate events in Groningen (Willacy et al., 2019). The use of a 3D model has the potential to further reduce the
travel-time residuals and therewith the location uncertainty.

6.3

Knowledge of travel-time uncertainty

The maps generated in this research are conditioned on the travel-time uncertainty model
as estimated in Section 2.6. The estimated travel-time uncertainty models are based on
deviations from a degree-2 polynomial model. These polynomial models act well to describe
mean travel times from a few kilometers epicentral distance onwards. Travel-time picks at
these distances are used especially to constrain the epicenter. Hence, the modeled epicentral
location uncertainty is likely reasonable. And also the depth uncertainty that is achieved
with picks from a few kilometers epicentral distance, is likely reasonable.
For more accurate depth location one needs picks at small epicentral distances. At these distances, the wavefield travels directly upwards from the source, leading to a travel-time curve
that is hyperbolic (Section 2.6.3). This hyperbolic part of the travel-time model cannot
easily be estimated or checked from the picks. For accurately modelling the arrival times at
the nearby stations, a good 1D velocity profile of the source region is needed. We assumed
that also in the hyperbolic part of the travel-time curve, similar uncertainty exists as with
polynomial part. In reality, there will be regional variation in the quality of the available
velocity information, or travel-time curves. The produced maps (near station locations)
show likely realistic σZ values for regions for which already quite good models exist for Pand S-wave velocities in the top few kilometers of the crust1 , like most parts of north of the
Netherlands. The model would have less variability than assumed when locating a source
in, e.g., the Zuidwending salt dome, for which the velocity is well known and local stations
exist. Hence, in practise a smaller σZ could be achieved than indicated, e.g., on Fig. 5.7(a).
On the other hand, in regions of, e.g., Limburg where there are still large unknowns in the
velocity structure of the upper crust, the σZ error (near receivers) would be larger than
indicated on the maps in Appendices B.2&D.2.
For the polynomial part of the travel-time model, one can well check whether it is a good
average description of the travel-time picks. For the hyperbolic part of the model, one often
has too little picks to do this check. This increases the chance that the used local 1D model
is a biased model and that hence the actual depth of the hypocenter is an outlier to the
modeled depth PDF.
For constraining the source depth, the P-S delay plays an important role (Section 3.2).
Hence, the depth uncertainty can be reduced by reducing both the P-wave and S-wave
travel-time uncertainties. S-wave travel-times, however, show a larger variability and hence
largest gains can be achieved by better modelling S-wave travel times. Generally, local Pwave velocities in the overburden are quite well known from seismic exploration campaigns.
Local S-wave velocity is much less well known. For reducing depth uncertainty, it is elementary to improve the local S-wave models, e.g., with tomographic methods or S-wave
logging.
1 Note that most of the travel time, and hence also most of the travel-time error, occurs where the wave
propagates slowest, which is in the upper crust.

89

6.4. MAGNITUDE OF COMPLETENESS

6.4

KNMI

Magnitude of completeness

The computed magnitude-of-completeness (MoC) map (Fig. 5.10) is probabilistic. In the
strict definition, MoC is the minimum magnitude for which the catalogue is complete. There
is no such thing as complete, however. Exceptional circumstances may exist (e.g., a large
storm) that causes seismic events with magnitudes above the MoC to be undetected. We
computed MoC as the lower magnitude that can very likely be detected and located. We
used the 90th percentile of the noise distribution to model whether at least 3 detections can
be made. When computing the MoC in hindsight, from a catalogue of events in a region, it
is thus likely that one finds a slightly higher MoC than indicated on Figs. 5.10 and 5.11.
The MoC is computed for a quite high noise level (P90). This means that when the noise
levels are average (P50) or low (P05), events can be detected and located below the MoC
threshold (Fig. 5.8).
The MoC-map was computed for a source depth of 3 km. If sources are shallower, the
amplitudes in the near field are higher, which would result in a somewhat smaller MoC.
On the other hand, for deeper seismicity, the MoC would be somewhat larger than shown
on Fig. 5.10. For regions with a small number of seismometers, the MoC would be almost
independent of source depth since amplitudes at larger distances are nearly source-depth
independent.
The previous versions of the MoC map that were published for the Netherlands (Kruiver
et al., 2021; KNMI , 2021) were based on different models and assumptions. Whether a
detection could be made, was based on an S-wave attenuation relation, in combination with
average noise levels that existed over the old seismic network (< 2004) in the Groningen region. The new map is based on newly derived P-wave GMPEs together with station-specific
noise distributions. Also it has been determined how much higher, on average, the signal
should be than the noise for an automatic pick to materialize (Section 2.4) and hard-rock
correction factors have been implemented where relevant (Section 2.5).
The latest MoC map based on old assumptions (KNMI , 2021) and the one derived in
this report (Fig. 5.10) are similar within the dense networks in Groningen, Twente and
North-Holland. In most other places, the new map shows slightly higher values. There
are two exceptions. First, the area encompassed by the network in South-Holland shows
significantly larger MoC values on the new map, due to very high noise conditions that were
found. Second, the southeast of the Netherlands shows smaller MoC values on the new MoC
map, due to very low noise levels that were found for stations in south Limburg and across
the border.

6.5

Future work

In Section 2.5 average hard-rock correction factors were estimated for sensors at, or near,
the Earth’s surface and at > 40 m depth. In this study, the same correction factors were
used for various sites. In reality, there is likely quite some variability between the hard-rock
sites. In a future study, site-specific correction factors could be established. E.g., station
WTSB is likely on much weaker rock than HGN and hence a smaller correction factor would
need to be used here.
At most sites with a vertical array of sensors, only the sensor with the best signal-to-noise
ratio (SNR) is used for detection. In the current operations, no advantage is taken of noisereduction capabilities of these vertical arrays. E.g, the arrays can be used to increase the
SNR of incoming body waves by forming ’enhanced seismograms’ (Ruigrok et al., 2016). A
topic of future research would be how much gain can be achieved by this processing and how
90

6.6. ACKNOWLEDGEMENTS

KNMI

the MoC and location accuracy maps would improve if the enhanced seismogram were to be
implemented. Similarly, also the gain of horizontal phased arrays, like within the LOFAR
setup (Drijkoningen et al., 2006) can be studied.
For seven of the broadband stations in the NL network, the sampling frequency has not
been optimized for detecting induced seismicity. As a consequence, we could not have a
complete picture of the noise distribution in the detection frequency band (Section 2.3.1).
When the sampling rate is increased at these stations, the noise levels need to be re-assessed.
It could be assessed whether changing the detection frequency band is favorable for improving the detection capabilities of the network.
In this study, the travel-time uncertainty was found to be distance dependent (Section 2.6.1).
This dependence was taken into account by computing an average travel-time uncertainty
over the collection of available picks (Section 5.1.1). Instead, the uncertainty pertaining
to each individual pick could be used within the covariance matrix (e.g., Eq. 3.14). This
likely solves the issue that sometimes the location uncertainty grows if more distant picks
are included (Fig. 5.1).
The noise conditions in the Netherlands change over time. The seismic noise levels are
expected to increase, with a population still growing and more anthropogenic activities near
the Earth’s surface. Especially the size of onshore wind turbines is expected to grow the
coming years. This makes it necessary to re-assess the station-specific noise levels in the
future. This will be taken into account with the yearly updates, as currently already exist
for the years 2021 and 2022.

6.6

Acknowledgements

The authors are grateful for the financial support provided by the Dutch State Supervision
of Mines through the project “Kaart voor Nederland van de onzekerheid op de lokalisatie van
aardbevingen door het landelijke KNMI netwerk”. This project was initiated together with
Denise de Vos and Niels Grobbe. They remained involved during the execution. They gave
valuable feedback during various project meetings and on the draft of this report. Läslo
Evers gave managerial support. Jesper Spetzler delivered the finite-difference generated
travel times as used for Fig. 2.23. Also, we would like to thank Hugues Djikpesse. He
supervised ENR during an internship back in 2007 and was instrumental in developing the
error-propagation approach as outlined in Section 3.1. Finally, we would like to thank
different members of the RDSA agile team who built the PDFPSD database and shaped
the detection list: Mathijs Koymans, Gert-Jan van der Hazel, Jordi Domingo-Ballesta, João
Paulo Pereira Zanetti and Luca Trani.

91

Bibliography
Atkinson, G. M., E. Yenier, N. Sharma, and V. Convertito (2016), Constraints on the neardistance saturation of ground-motion amplitudes for small-to-moderate induced earthquakes, Bulletin of the Seismological Society of America, 106 (5), 2104–2111.
Bommer, J. J., P. J. Stafford, and M. Ntinalexis (2019), Updated emperical GMPEs for
PGV from Groningen earthquakes — March 2019, NAM Study Report.
Bondár, I., S. C. Myers, E. R. Engdahl, and E. A. Bergman (2004), Epicentre accuracy
based on seismic network criteria, Geophysical Journal International, 156 (3), 483–496.
Bormann, P. (1998), Conversion and comparability of data presentations on seismic background noise, Journal of Seismology, 2 (1), 37–45.
Buijze, L., H. Veldkamp, and B. Wassing (2023), Comparison of hydrocarbon and geothermal
energy production in the Netherlands: reservoir characteristics, pressure and temperature
changes, and implications for fault reactivation, Netherlands Journal of Geosciences, 102,
e7.
Buijze, L., et al. (2019), Review of induced seismicity in geothermal systems worldwide
and implications for geothermal systems in the Netherlands, Netherlands Journal of Geosciences, 98, e13.
Chung, D., and D. Bernreuter (1981), Regional relationships among earthquake magnitude
scales, Reviews of Geophysics, 19 (4), 649–663.
Dost, B., F. Goutbeek, T. van Eck, and D. Kraaijpoel (2012), Monitoring induced seismicity
in the north of the Netherlands: status report 2010, KNMI Scientific Report, pp. 2012–03.
Dost, B., E. Ruigrok, and J. Spetzler (2017), Development of probabilistic seismic hazard assessment for the Groningen gas field, Netherlands Journal of Geosciences, 96 (5),
s235–s245, doi:10.1017/njg.2017.20.
Dost, B., A. van Stiphout, D. Kühn, M. Kortekaas, E. Ruigrok, and S. Heimann (2020),
Probabilistic moment tensor inversion for hydrocarbon-induced seismicity in the Groningen gas field, the Netherlands, part 2: Application, Bulletin of the Seismological Society
of America, 110 (5), 2112–2123.
Douglas, J., B. Edwards, V. Convertito, N. Sharma, A. Tramelli, D. Kraaijpoel, B. M.
Cabrera, N. Maercklin, and C. Troise (2013), Predicting ground motion from induced
earthquakes in geothermal areas, Bulletin of the Seismological Society of America, 103 (3),
1875–1897.
Drijkoningen, G., J. Brouwer, J. Kooijman, G. Steenbergen, B. Dost, and A. Huijgen (2006),
4C seismic monitoring on land in the Netherlands – Results for acquisition design, 68th
EAGE Conference and Exhibition, extended abstracts, P-291.

92

BIBLIOGRAPHY

KNMI

Freudenreich, Y., S. J. Oates, and W. Berlang (2012), Microseismic feasibility studies–
assessing the probability of success of monitoring projects, Geophysical Prospecting, 60 (6),
1043–1053.
GEMPA (2022), SeisComp, https://www.seiscomp.de/doc/, accessed: 2022-07-06.
Haak, H. (Ed.) (1993), Summary of the final report on a multidisciplinary study of the
relation between Gas production and earthquakes in the northern part of the Netherlands,
p. 16, Royal Netherlands Meteorological Institute (de Bilt), ISBN 90-369-2052-3.
Havskov, J., and G. Alguacil (2016), Instrumentation in Earthquake Seismology, Springer
International Publishing, doi:10.1007/978-3-319-21314-9.
Hinzen, K.-G., S. K. Reamer, and C. Fleischer (2021), Seismicity in the Northern Rhine
Area (1995–2018), Journal of Seismology, 25 (2), 351–367.
Hofman, L., E. Ruigrok, B. Dost, and H. Paulssen (2017), A shallow velocity model for the
Groningen area in the Netherlands, Journal of Geophysical Research: Solid Earth, 122,
8035–8050, doi:10.1002/2017JB014419.
Jagt, L., E. Ruigrok, and H. Paulssen (2017), Relocation of clustered earthquakes
in the Groningen gas field, Netherlands Journal of Geosciences, 96 (5), s163–s173,
doi:10.1017/njg.2017.12.
Kanamori, H., J. Mori, E. Hauksson, T. H. Heaton, L. K. Hutton, and L. M. Jones (1993),
Determination of earthquake energy release and ML using TERRAscope, Bulletin of the
Seismological Society of America, 83 (2), 330–346.
Kinnaert, X., E. Gaucher, U. Achauer, and T. Kohl (2016), Modelling earthquake location
errors at a reservoir scale: a case study in the Upper Rhine Graben, Geophysical Journal
International, 206 (2), 861–879.
KNMI (1993), Netherlands Seismic and Acoustic Network, Royal Netherlands Meteorological Institute (KNMI), Other/Seismic Network, doi:10.21944/e970fd34-23b9-3411-b366e4f72877d2c5.
KNMI
(2021),
Netherlands
earthquake
magnitude
completeness
map,
https://dataplatform.knmi.nl/dataset/
netherlands-earthquake-magnitude-completeness-map-1, KNMI Data Platform,
Accessed: 2023-03-01.
KNMI
(2023a),
Magnitude
of
Completeness
Maps
for
the
Netherlands,
https://dataplatform.knmi.nl/dataset/
netherlands-earthquake-magnitude-completeness-1-0, KNMI Data Platform,
Accessed: 2023-07-24.
KNMI
(2023b),
Earthquake
Location
Uncertainty
Maps
for
the
Netherlands,
https://dataplatform.knmi.nl/dataset/
netherlands-earthquake-location-uncertainty-1-0,
KNMI Data Platform,
Accessed: 2023-07-24.
Koymans, M., J. D. Ballesta, E. Ruigrok, R. Sleeman, L. Trani, and L. Evers (2021), Performance assessment of geophysical instrumentation through the automated analysis of
power spectral density estimates, Earth and Space Science, e2021EA001675.
Kraaijpoel, D., and B. Dost (2013), Implications of salt-related propagation and mode conversion effects on the analysis of induced seismicity, Journal of Seismology, 17 (1), 95–107.

93

BIBLIOGRAPHY

KNMI

Kruiver, P., B. Dost, E. Ruigrok, J. Spetzler, G.-J. van den Hazel, T. Campmans, and L. Evers (2021), Ontwerp basismonitoringsnetwerk voor geı̈nduceerde seismiciteit ten gevolge
van mijnbouwactiviteiten, KNMI Technical Report, TR-391.
Kruiver, P. P., et al. (2017), An integrated shear-wave velocity model for the Groningen
gas field, the Netherlands, Bulletin of Earthquake Engineering, doi: 10.1007/s10518-0170105-y.
Kühn, D., S. Heimann, M. P. Isken, E. Ruigrok, and B. Dost (2020), Probabilistic moment tensor inversion for hydrocarbon-induced seismicity in the Groningen gas field, the
Netherlands, part 1: Testing, Bulletin of the Seismological Society of America, 110 (5),
2095–2111.
Langet, N., A. Maggi, A. Michelini, and F. Brenguier (2014), Continuous kurtosis-based migration for seismic event detection and location, with application to Piton de la Fournaise
Volcano, La Reunion, Bulletin of the Seismological Society of America, 104 (1), 229–246.
Lienert, B. R., E. Berg, and L. N. Frazer (1986), HYPOCENTER: An earthquake location method using centered, scaled, and adaptively damped least squares, Bulletin of the
Seismological Society of America, 76 (3), 771–783.
Lomax, A., A. Michelini, A. Curtis, and R. Meyers (2009), Earthquake location, direct,
global-search methods, Encyclopedia of complexity and systems science, 5, 2449–2473.
McNamara, D. E. (2004), Ambient noise levels in the continental united states, Bulletin of
the Seismological Society of America, 94 (4), 1517–1527, doi:10.1785/012003001.
Mijnlieff, H. F. (2020), Introduction to the geothermal play and reservoir geology of the
Netherlands, Netherlands Journal of Geosciences, 99, e2.
Muntendam-Bos, A., and N. Grobbe (2022), Data-driven spatiotemporal assessment of the
event-size distribution of the Groningen extraction-induced seismicity catalogue, Scientific
Reports, 12 (1), 1–11.
Naranjo, D., D. Draganov, K. Polychronopoulou, M. de Bas, and C. Weemstra (2022),
Seismic monitoring of Nature’s Heat geothermal reservoir in Kwintsheul (Netherlands).
Peterson, J. R. (1993), Observations and modeling of seismic background noise, doi:10.3133/
ofr93322.
Robertsson, J. O., J. O. Blanch, and W. W. Symes (1994), Viscoelastic finite-difference
modeling, Geophysics, 59 (9), 1444–1456.
Romijn, R. (2017), Groningen velocity model 2017, Tech. rep., NAM (Nederlands Aardolie
Maatschappij).
Ruigrok, E., and B. Dost (2020), Advice on the computation of peak-ground-velocity confidence regions for events in gas fields other than the Groningen gas field, KNMI Technical
Report, TR-386.
Ruigrok, E., P. Paulssen, and B. Dost (2016), Enhanced seismograms from borehole arrays,
European Seismological Commission 35th General Assembly, eSC:2016-497.
Scudero, S., C. Marcocci, and A. D’Alessandro (2021), Insights on the italian seismic network
from location uncertainties, Journal of Seismology, 25 (4), 1061–1076.
Spetzler, J., and B. Dost (2017), Hypocentre estimation of induced earthquakes in Groningen, Geophysical Journal International, 209 (1), 453–465.

94

BIBLIOGRAPHY

KNMI

Tarantola, A. (2005), Inverse Problem Theory and Methods for Model Parameter Estimation,
SIAM, Philadelphia.
Ten Veen, J., M. Den Dulk, S. Van Gessel, R. Van Ede, H. Doornenbal, J. Hettelaar, and
E. Hintersberger (2022), A novel classification system and database for faults and tectonic
features in the netherlands, in 83rd EAGE Annual Conference & Exhibition, vol. 2022,
pp. 1–5, European Association of Geoscientists & Engineers.
Van Dalfsen, W., J. Doornenbal, S. Dortland, and J. Gunnink (2006), A comprehensive
seismic velocity model for the netherlands based on lithostratigraphic layers, Netherlands
Journal of Geosciences, 85 (4), 277.
Van Eck, T., F. Goutbeek, H. Haak, and B. Dost (2004), Seismic hazard due to small shallow
induced earthquakes, KNMI Scientific Report, WR-2004-01.
Van Eijs, R., F. Mulders, M. Nepveu, C. Kenter, and B. Scheffers (2006), Correlation between
hydrocarbon reservoir properties and induced seismicity in the netherlands, Engineering
Geology, 84 (3), 99–111.
van Ginkel, J., E. Ruigrok, and R. Herber (2020), Using horizontal-to-vertical spectral ratios
to construct shear-wave velocity profiles, Solid Earth, 11 (6), 2015–2030.
van Ginkel, J., E. Ruigrok, J. Stafleu, and R. Herber (2022a), Development of a seismic siteresponse zonation map for the Netherlands, Natural Hazards and Earth System Sciences,
22 (1), 41–63.
van Ginkel, J., E. Ruigrok, R. Wentinck, and R. Herber (2022b), Amplification behaviour
of compressional waves in unconsolidated sediments, Frontiers in Earth Science, p. 475.
van Thienen-Visser, K., M. Nepveu, and J. Hettelaar (2012), Deterministische hazard analyse voor geı̈nduceerde seismiciteit in nederland, TNO-rapport, R10198.
van Thienen-Visser, K., D. Sijacic, J.-D. van Wees, D. Kraaijpoel, and J. Roholl
(2016), Groningen field 2013 to present, gas production and induced seismicity, Tech.
rep., TNO 2016 R10425, https://www.nlog.nl/sites/default/files/final_tno%
20report_2016_r10425.pdf.
Vörös, R., and S. Baisch (2022), Induced seismicity and seismic risk management–a showcase from the Californië geothermal field (the Netherlands), Netherlands Journal of Geosciences, 101, e15.
Wang, B., W. Shi, and Z. Miao (2015), Confidence analysis of standard deviational ellipse
and its extension into higher dimensional euclidean space, PloS one, 10 (3), e0118,537.
Willacy, C., E. van Dedem, S. Minisini, J. Li, J.-W. Blokland, I. Das, and A. Droujinine
(2019), Full-waveform event location and moment tensor inversion for induced seismicity,
Geophysics, 84 (2), KS39–KS57.
Yudistira, T., H. Paulssen, and J. Trampert (2017), The crustal structure beneath the
Netherlands derived from ambient seismic noise, Tectonophysics, 721, 361–371.
Zhou, H.-w. (1994), Rapid three-dimensional hypocentral determination using a master
station method, Journal of Geophysical Research: Solid Earth, 99 (B8), 15,439–15,455.

95

A 2021 Station list
Table A.1: The year 2021 list of 200 stations that is used in this study. For details on the
selection procedure, see Section 2.1. The locations and sensor types are shown on Fig. 2.2.
In the last column it is indicated whether the sensor is situated in hardrock conditions.
Station name
NL.ALK1..HGZ
NL.ALK2..HGZ
NL.ARCN..BHZ
NL.ASS1..HGZ
NL.ASS2..HGZ
NL.BAPP..HGZ
NL.BER1..HGZ
NL.BER2..HGZ
NL.BFB2..HGZ
NL.BGAR..HGZ
NL.BHAR..HGZ
NL.BHKS..HGZ
NL.BKMZ..HGZ
NL.BLIJE..HGZ
NL.BLOP..HGZ
NL.BMD2..HGZ
NL.BOWW..HGZ
NL.BSTD.00.HGZ
NL.BWIR..HGZ
NL.BWSE..HGZ
NL.COE2..HGZ
NL.COE3..HGZ
NL.COE4..HGZ
NL.COE..HGZ
NL.DBLS..HGZ
NL.DBN.01.BHZ
NL.DON..HGZ
NL.DR014..HHZ
NL.DR024..HHZ
NL.DR034..HHZ
NL.DRA..HGZ
NL.EETW..HGZ
NL.ELE..HGZ
NL.ENM4.01.HHZ
NL.ENV2.01.HHZ
NL.FDKD..HGZ

Latitude [deg]
52.6338
52.6510
51.5013
52.9778
53.0143
53.3148
52.6614
52.6478
53.1875
53.3679
53.2292
53.2919
53.2901
53.3603
53.3339
53.3533
53.3405
53.3122
53.3226
53.3444
52.6328
52.6704
52.6598
52.6606
52.8391
52.1017
53.1005
52.7987
52.7034
52.8976
53.1202
53.3348
52.9566
53.4064
52.8944
52.8483

Longitude [deg]
4.7444
4.7234
6.1942
6.5574
6.5706
6.8354
4.7037
4.7150
6.7655
6.7136
6.7090
6.7850
6.3189
5.8492
6.7466
6.6472
6.8125
6.6921
6.7861
6.7099
6.7057
6.6525
6.7451
6.7836
6.0497
5.1767
6.5429
6.1785
6.5081
6.3084
6.1033
6.1054
6.5776
6.4817
6.6337
6.1788

96

Depth [m]
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
200
200
200
0
0
0
200
100
0

Sensor
accelerometer
accelerometer
broadband
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
broadband
accelerometer
geophone
geophone
geophone
accelerometer
accelerometer
accelerometer
geophone
geophone
accelerometer

Hardrock
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

KNMI

NL.FR014..HHZ
NL.FR024..HHZ
NL.FR034..HHZ
NL.FR044..HHZ
NL.FR054..HHZ
NL.FR064..HHZ
NL.FSW5.01.HHZ
NL.G014..HHZ
NL.G034..HHZ
NL.G044..HHZ
NL.G054..HHZ
NL.G064..HHZ
NL.G074..HHZ
NL.G084..HHZ
NL.G094..HHZ
NL.G104..HHZ
NL.G114..HHZ
NL.G124..HHZ
NL.G134..HHZ
NL.G144..HHZ
NL.G164..HHZ
NL.G174..HHZ
NL.G184..HHZ
NL.G194..HHZ
NL.G204..HHZ
NL.G214..HHZ
NL.G224..HHZ
NL.G234..HHZ
NL.G244..HHZ
NL.G254..HHZ
NL.G264..HHZ
NL.G274..HHZ
NL.G284..HHZ
NL.G294..HHZ
NL.G304..HHZ
NL.G324..HHZ
NL.G334..HHZ
NL.G344..HHZ
NL.G354..HHZ
NL.G364..HHZ
NL.G374..HHZ
NL.G384..HHZ
NL.G394..HHZ
NL.G404..HHZ
NL.G414..HHZ
NL.G424..HHZ
NL.G434..HHZ
NL.G444..HHZ
NL.G454..HHZ
NL.G464..HHZ
NL.G474..HHZ
NL.G484..HHZ
NL.G494..HHZ

53.1424
52.8559
52.9815
52.8761
53.0705
53.0007
53.2135
53.4425
53.4159
53.4147
53.4184
53.4111
53.3946
53.3954
53.3878
53.3785
53.3669
53.3699
53.3489
53.3586
53.3385
53.3362
53.3306
53.3378
53.3353
53.3134
53.3103
53.3095
53.2900
53.3081
53.2821
53.2832
53.2722
53.2814
53.2792
53.2498
53.2499
53.2525
53.2549
53.2651
53.2740
53.2005
53.2253
53.2250
53.2234
53.2217
53.2319
53.2042
53.1995
53.1978
53.1964
53.1973
53.1765

5.8868
5.9304
5.7126
6.1187
5.8685
6.0237
7.1195
6.7312
6.5976
6.6747
6.7713
6.8709
6.5458
6.6439
6.7245
6.8044
6.8891
6.5714
6.6686
6.7708
6.5314
6.6363
6.7217
6.8170
6.8952
6.5937
6.6776
6.7675
6.8562
6.9433
6.5458
6.6281
6.7292
6.7870
6.8951
6.5792
6.6708
6.7651
6.8621
6.9245
7.0350
6.5476
6.7241
6.8099
6.8919
6.9886
7.0791
6.6840
6.7645
6.8498
6.9413
7.0290
6.7306

97

200
200
200
200
200
200
300
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200

geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone

0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

KNMI

NL.G504..HHZ
NL.G514..HHZ
NL.G524..HHZ
NL.G534..HHZ
NL.G544..HHZ
NL.G554..HHZ
NL.G564..HHZ
NL.G574..HHZ
NL.G584..HHZ
NL.G594..HHZ
NL.G604..HHZ
NL.G614..HHZ
NL.G624..HHZ
NL.G634..HHZ
NL.G644..HHZ
NL.G654..HHZ
NL.G664..HHZ
NL.G674..HHZ
NL.G684..HHZ
NL.G694..HHZ
NL.G704..HHZ
NL.G710..HGZ
NL.G720..HGZ
NL.G730..HGZ
NL.G740..HGZ
NL.G750..HGZ
NL.G760..HGZ
NL.G770..HGZ
NL.G780..HGZ
NL.G790..HGZ
NL.G800..HGZ
NL.G81B..HHZ
NL.G82B..HHZ
NL.G83B..HHZ
NL.G84B..HHZ
NL.GK014..HHZ
NL.GK024..HHZ
NL.GK034..HHZ
NL.GK044..HHZ
NL.GRP..HGZ
NL.HGN.02.BHZ
NL.HLB.11.EPZ
NL.HLG.11.EPZ
NL.HLH.11.EPZ
NL.HLK.11.EPZ
NL.HOEV..HGZ
NL.HRDB..HGZ
NL.HRKB..BHZ
NL.HWF4.01.HHZ
NL.J01..HHZ
NL.J02.01.HHZ
NL.KLWL..HGZ
NL.LUTT..HGZ

53.1745
53.1732
53.1699
53.1674
53.1243
53.1433
53.1484
53.1351
53.1103
53.1160
53.2979
53.3819
53.3870
53.3314
53.2207
53.1692
53.1794
53.3213
53.2351
53.2128
53.2407
53.1483
53.2462
53.2719
53.3043
53.3078
53.3507
53.3403
53.3648
53.3836
53.4154
53.3103
53.2253
53.1484
53.2207
53.2906
53.3010
53.2827
53.2462
53.1788
50.7640
53.1465
53.1491
53.1404
53.1450
52.9025
52.5770
51.1879
53.0710
52.7650
52.7255
51.7755
52.6133

6.8093
6.8848
6.9837
7.0801
6.7580
6.8461
6.9309
7.0222
6.9071
6.9877
7.0752
6.6813
6.8419
6.9109
7.0380
6.9543
6.5852
6.8228
6.6347
6.6438
6.5379
6.6507
6.4792
6.4301
6.3894
6.4875
6.3901
6.3304
6.4894
6.3319
6.4876
6.6776
6.7241
6.9309
7.0380
6.2471
6.2927
6.3056
6.3261
5.9779
5.9317
7.0036
7.0071
7.0065
7.0096
6.1055
6.6259
6.1678
6.3512
5.0545
5.1948
4.4406
6.5715

98

200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
0
0
0
0
0
0
0
0
0
0
100
100
100
100
200
200
200
200
0
4
60
58
60
60
0
0
0
200
0
0
0
0

geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
accelerometer
broadband
geophone
geophone
geophone
geophone
accelerometer
accelerometer
broadband
geophone
geophone
geophone
accelerometer
accelerometer

0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0

KNMI

NL.MAME..BHZ
NL.MAR..HGZ
NL.N014..HHZ
NL.N024..HHZ
NL.N034..HHZ
NL.NGA..HGZ
NL.NH014..HHZ
NL.NIW4.01.HHZ
NL.NLDW..HGZ
NL.ODBK..HGZ
NL.ODHP..HGZ
NL.OPLO.01.BHZ
NL.OSTM..HGZ
NL.OTL2.01.HHZ
NL.PET..HGZ
NL.PPB3.01.HHZ
NL.ROLD..BHZ
NL.SL014..HHZ
NL.SPY4.01.HHZ
NL.SRHV..HGZ
NL.T024..HHZ
NL.T034..HHZ
NL.T044..HHZ
NL.T054..HHZ
NL.T064..HHZ
NL.T084..HHZ
NL.TERZ.01.BHZ
NL.TWZL..HGZ
NL.UTRP..HGZ
NL.VBG4.01.HHZ
NL.VHP..HGZ
NL.VKB..BHZ
NL.VLD..HGZ
NL.VLW3.01.HHZ
NL.VNDM4..HHZ
NL.VNWD..HGZ
NL.VPR..HGZ
NL.VRS..HGZ
NL.WSVN..HGZ
NL.WTSB.01.BHZ
NL.ZDL..HGZ
NL.ZH014..HHZ
NL.ZH024..HHZ
NL.ZH034..HHZ
NL.ZH044..HHZ
NL.ZH064..HHZ
NL.ZLV4.01.HHZ
NL.ZWE2.00.HHZ
GE.IBBN..HHZ
GR.BUG..HHZ
GR.RAST..EHZ
BE.BOST..HHZ
BE.MEM..HHZ

50.8000
53.1540
53.1185
53.0720
53.1110
53.1364
52.5657
53.3548
52.0050
52.9405
52.8852
51.5888
53.1875
52.6289
52.7706
52.6526
50.8694
52.8551
53.4098
53.1833
52.4634
52.4297
52.4375
52.3806
52.3432
52.3893
50.7568
53.2325
53.0969
52.5440
53.0957
50.8669
52.8668
52.9682
53.1151
53.2368
51.8915
53.0360
52.8204
51.9663
53.0791
51.9137
52.0744
51.9783
51.8417
51.9781
53.0921
53.1885
52.3063
51.4455
52.9096
51.2382
50.6087

5.9727
6.2631
6.4733
6.4530
6.4045
6.0257
4.9965
6.0430
4.1957
6.1339
6.0249
5.8121
6.0857
4.8227
4.6595
4.6700
6.0847
4.8683
6.7838
6.1849
6.9787
6.7918
6.8905
6.9057
6.9097
6.9965
5.9061
6.0927
6.1842
6.6693
5.9553
5.7847
6.2157
7.0972
6.8390
5.9984
4.3715
6.6046
6.1983
6.7989
6.7026
4.2881
4.4884
4.4001
4.3626
4.2204
6.7533
5.6045
7.7592
7.2643
7.7164
2.9387
6.0096

99

0
0
200
200
200
0
200
120
0
0
0
0
0
100
0
165
0
200
120
0
200
200
200
200
200
200
250
0
0
192
0
0
0
150
200
0
0
0
0
0
0
200
200
200
200
200
200
60
1
0
0
304
0

broadband
accelerometer
geophone
geophone
geophone
accelerometer
geophone
geophone
accelerometer
accelerometer
accelerometer
broadband
accelerometer
geophone
accelerometer
geophone
broadband
geophone
geophone
accelerometer
geophone
geophone
geophone
geophone
geophone
geophone
broadband
accelerometer
accelerometer
geophone
accelerometer
broadband
accelerometer
geophone
geophone
accelerometer
accelerometer
accelerometer
accelerometer
broadband
accelerometer
geophone
geophone
geophone
geophone
geophone
geophone
geophone
broadband
broadband
geophone
broadband
broadband

1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
1
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
1
1
0
1
1

KNMI

BE.BEBN..HHZ
BE.OPTB..HHZ
BE.RCHB..HHZ
BE.UCC..HHZ
NL.CHA5.01.HHZ

50.7970
51.1115
50.1552
50.7972
51.5043

5.6778
5.6360
5.2268
4.3604
4.9212

100

0
379
43
141
141

broadband
broadband
broadband
broadband
geophone

1
1
1
1
0

B 2021 NL-wide maps with locationuncertainty parameters
B.1

Largest-horizontal-location-uncertainty (σ1 ) maps

101

B.1. LARGEST-HORIZONTAL-LOCATION-UNCERTAINTY (σ1 ) MAPS

KNMI

Figure B.1: Spatial distribution of σ1 in 2021 for M=0.5 (color map). White means that
events with M=0.5 can likely not be detected. Surface sensors are shown with squares,
deep geophones with orange triangles and broadbands with green triangles. Thin black lines
denote province borders within the Netherlands. Thick black lines are country borders. The
coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in kilometers. Dark red
depicts uncertainties of 2800 m and higher.

102

B.1. LARGEST-HORIZONTAL-LOCATION-UNCERTAINTY (σ1 ) MAPS

KNMI

Figure B.2: Spatial distribution of σ1 in 2021 for M=1.0 (color map). White means that
events with M=1.0 can likely not be detected. Surface sensors are shown with squares,
deep geophones with orange triangles and broadbands with green triangles. Thin black lines
denote province borders within the Netherlands. Thick black lines are country borders. The
coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

103

B.1. LARGEST-HORIZONTAL-LOCATION-UNCERTAINTY (σ1 ) MAPS

KNMI

Figure B.3: Spatial distribution of σ1 in 2021 for M=1.5 (color map). White means that
events with M=1.5 can likely not be detected. Surface sensors are shown with squares,
deep geophones with orange triangles and broadbands with green triangles. Thin black lines
denote province borders within the Netherlands. Thick black lines are country borders. The
coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

104

B.1. LARGEST-HORIZONTAL-LOCATION-UNCERTAINTY (σ1 ) MAPS

KNMI

Figure B.4: Spatial distribution of σ1 in 2021 for M=2.0 (color map). White means that
events with M=2.0 can likely not be detected. Surface sensors are shown with squares,
deep geophones with orange triangles and broadbands with green triangles. Thin black lines
denote province borders within the Netherlands. Thick black lines are country borders. The
coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

105

B.1. LARGEST-HORIZONTAL-LOCATION-UNCERTAINTY (σ1 ) MAPS

KNMI

Figure B.5: Spatial distribution of σ1 in 2021 for M=2.5 (color map). White means that
events with M=2.5 can likely not be detected. Surface sensors are shown with squares,
deep geophones with orange triangles and broadbands with green triangles. Thin black lines
denote province borders within the Netherlands. Thick black lines are country borders. The
coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

106

B.1. LARGEST-HORIZONTAL-LOCATION-UNCERTAINTY (σ1 ) MAPS

KNMI

Figure B.6: Spatial distribution of σ1 in 2021 for M=3.0 (color map). Surface sensors
are shown with squares, deep geophones with orange triangles and broadbands with green
triangles. Thin black lines denote province borders within the Netherlands. Thick black lines
are country borders. The coordinates are in the local rectangular grid (Rijksdriehoekstelsel)
in kilometers.

107

B.1. LARGEST-HORIZONTAL-LOCATION-UNCERTAINTY (σ1 ) MAPS

KNMI

Figure B.7: Spatial distribution of σ1 in 2021 for M=3.5 (color map). Surface sensors
are shown with squares, deep geophones with orange triangles and broadbands with green
triangles. Thin black lines denote province borders within the Netherlands. Thick black lines
are country borders. The coordinates are in the local rectangular grid (Rijksdriehoekstelsel)
in kilometers.

108

B.1. LARGEST-HORIZONTAL-LOCATION-UNCERTAINTY (σ1 ) MAPS

KNMI

Figure B.8: Spatial distribution of σ1 in 2021 for M=4.0 (color map). Surface sensors
are shown with squares, deep geophones with orange triangles and broadbands with green
triangles. Thin black lines denote province borders within the Netherlands. Thick black lines
are country borders. The coordinates are in the local rectangular grid (Rijksdriehoekstelsel)
in kilometers.

109

B.1. LARGEST-HORIZONTAL-LOCATION-UNCERTAINTY (σ1 ) MAPS

110

KNMI

B.2. VERTICAL-LOCATION-UNCERTAINTY (σZ ) MAPS

B.2

KNMI

Vertical-location-uncertainty (σZ ) maps

Figure B.9: Spatial distribution of σZ in 2021 for M=0.5 (color map). Surface sensors
are shown with squares, deep geophones with orange triangles and broadbands with green
triangles. White means that events with M=0.5 can likely not be detected. Thin black lines
denote province borders within the Netherlands. Thick black lines are country borders. The
coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in kilometers. Dark red
depicts uncertainties of 3500 m and higher.
111

B.2. VERTICAL-LOCATION-UNCERTAINTY (σZ ) MAPS

KNMI

Figure B.10: Spatial distribution of σZ in 2021 for M=1.0 (color map). Surface sensors
are shown with squares, deep geophones with orange triangles and broadbands with green
triangles. White means that events with M=1.0 can likely not be detected. Uncertainties
of 3500 m and higher are shown in red. Thin black lines denote province borders within
the Netherlands. Thick black lines are country borders. The coordinates are in the local
rectangular grid (Rijksdriehoekstelsel) in kilometers.

112

B.2. VERTICAL-LOCATION-UNCERTAINTY (σZ ) MAPS

KNMI

Figure B.11: Spatial distribution of σZ in 2021 for M=1.5 (color map). Surface sensors
are shown with squares, deep geophones with orange triangles and broadbands with green
triangles. White means that events with M=1.5 can likely not be detected. Thin black lines
denote province borders within the Netherlands. Thick black lines are country borders. The
coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

113

B.2. VERTICAL-LOCATION-UNCERTAINTY (σZ ) MAPS

KNMI

Figure B.12: Spatial distribution of σZ in 2021 for M=2.0 (color map). Surface sensors
are shown with squares, deep geophones with orange triangles and broadbands with green
triangles. White means that events with M=2.0 can likely not be detected. Thin black lines
denote province borders within the Netherlands. Thick black lines are country borders. The
coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

114

B.2. VERTICAL-LOCATION-UNCERTAINTY (σZ ) MAPS

KNMI

Figure B.13: Spatial distribution of σZ in 2021 for M=2.5 (color map). Surface sensors
are shown with squares, deep geophones with orange triangles and broadbands with green
triangles. White means that events with M=2.5 can likely not be detected. Thin black lines
denote province borders within the Netherlands. Thick black lines are country borders. The
coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

115

B.2. VERTICAL-LOCATION-UNCERTAINTY (σZ ) MAPS

KNMI

Figure B.14: Spatial distribution of σZ in 2021 for M=3.0 (color map). Surface sensors
are shown with squares, deep geophones with orange triangles and broadbands with green
triangles. Thin black lines denote province borders within the Netherlands. Thick black lines
are country borders. The coordinates are in the local rectangular grid (Rijksdriehoekstelsel)
in kilometers.

116

B.2. VERTICAL-LOCATION-UNCERTAINTY (σZ ) MAPS

KNMI

Figure B.15: Spatial distribution of σZ in 2021 for M=3.5 (color map). Surface sensors
are shown with squares, deep geophones with orange triangles and broadbands with green
triangles. Thin black lines denote province borders within the Netherlands. Thick black lines
are country borders. The coordinates are in the local rectangular grid (Rijksdriehoekstelsel)
in kilometers.

117

B.2. VERTICAL-LOCATION-UNCERTAINTY (σZ ) MAPS

KNMI

Figure B.16: Spatial distribution of σZ in 2021 for M=4.0 (color map). Surface sensors
are shown with squares, deep geophones with orange triangles and broadbands with green
triangles. Thin black lines denote province borders within the Netherlands. Thick black lines
are country borders. The coordinates are in the local rectangular grid (Rijksdriehoekstelsel)
in kilometers.

118

B.3. AZIMUTHAL-GAP (ψ) MAPS

B.3

KNMI

Azimuthal-gap (ψ) maps

Figure B.17: Spatial distribution of azimuthal gap ψ in 2021 for M=0.5 (color map). White
means that events with M=0.5 can likely not be detected. Thin black lines denote province
borders within the Netherlands. Thick black lines are country borders. The coordinates are
in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

119

B.3. AZIMUTHAL-GAP (ψ) MAPS

KNMI

Figure B.18: Spatial distribution of azimuthal gap ψ in 2021 for M=1.0 (color map). White
means that events with M=1.0 can likely not be detected. Thin black lines denote province
borders within the Netherlands. Thick black lines are country borders. The coordinates are
in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

120

B.3. AZIMUTHAL-GAP (ψ) MAPS

KNMI

Figure B.19: Spatial distribution of azimuthal gap ψ in 2021 for M=1.5 (color map). White
means that events with M=1.5 can likely not be detected. Thin black lines denote province
borders within the Netherlands. Thick black lines are country borders. The coordinates are
in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

121

B.3. AZIMUTHAL-GAP (ψ) MAPS

KNMI

Figure B.20: Spatial distribution of azimuthal gap ψ in 2021 for M=2.0 (color map). White
means that events with M=2.0 can likely not be detected. Thin black lines denote province
borders within the Netherlands. Thick black lines are country borders. The coordinates are
in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

122

B.3. AZIMUTHAL-GAP (ψ) MAPS

KNMI

Figure B.21: Spatial distribution of azimuthal gap ψ in 2021 for M=2.5 (color map). White
means that events with M=2.5 can likely not be detected. Thin black lines denote province
borders within the Netherlands. Thick black lines are country borders. The coordinates are
in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

123

B.3. AZIMUTHAL-GAP (ψ) MAPS

KNMI

Figure B.22: Spatial distribution of azimuthal gap ψ in 2021 for M=3.0 (color map).
Thin black lines denote province borders within the Netherlands. Thick black lines are
country borders. The coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in
kilometers.

124

B.3. AZIMUTHAL-GAP (ψ) MAPS

KNMI

Figure B.23: Spatial distribution of azimuthal gap ψ in 2021 for M=3.5 (color map).
Thin black lines denote province borders within the Netherlands. Thick black lines are
country borders. The coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in
kilometers.

125

B.3. AZIMUTHAL-GAP (ψ) MAPS

KNMI

Figure B.24: Spatial distribution of azimuthal gap ψ in 2021 for M=4.0 (color map).
Thin black lines denote province borders within the Netherlands. Thick black lines are
country borders. The coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in
kilometers.

126

C 2022 Station list
Table C.1: The year 2022 list of 217 stations. For details on the selection procedure, see
Sections 2.1 and 5.4.3. In the last column it is indicated whether the sensor is situated in
hardrock conditions.
Station name
NL.ALK1..HGZ
NL.ALK2..HGZ
NL.ARCN..BHZ
NL.ASS1..HGZ
NL.ASS2..HGZ
NL.BAPP..HGZ
NL.BER1..HGZ
NL.BER2..HGZ
NL.BFB2..HGZ
NL.BGAR..HGZ
NL.BHKS..HGZ
NL.BKMZ..HGZ
NL.BLIJE..HGZ
NL.BLOP..HGZ
NL.BMD2..HGZ
NL.BOWW..HGZ
NL.BRTL..HGZ
NL.BSTD.00.HGZ
NL.BWIR..HGZ
NL.BWSE..HGZ
NL.COE2..HGZ
NL.COE3..HGZ
NL.COE4..HGZ
NL.COE..HGZ
NL.DBLS..HGZ
NL.DBN.01.BHZ
NL.DOK..HGZ
NL.DON..HGZ
NL.DR014..HHZ
NL.DR024..HHZ
NL.DR034..HHZ
NL.DRA..HGZ
NL.EETW..HGZ
NL.ELE..HGZ
NL.ENM4.01.HHZ
NL.ENV2.01.HHZ

Latitude [deg]
52.6338
52.6510
51.5013
52.9778
53.0143
53.3148
52.6614
52.6478
53.1875
53.3679
53.2919
53.2901
53.3603
53.3339
53.3533
53.3405
53.2380
53.3122
53.3226
53.3444
52.6328
52.6704
52.6598
52.6606
52.8391
52.1017
53.3273
53.1005
52.7987
52.7034
52.8976
53.1202
53.3348
52.9566
53.4064
52.8944

Longitude [deg]
4.7444
4.7234
6.1942
6.5574
6.5706
6.8354
4.7037
4.7150
6.7655
6.7136
6.7850
6.3189
5.8492
6.7466
6.6472
6.8125
6.3932
6.6921
6.7861
6.7099
6.7057
6.6525
6.7451
6.7836
6.0497
5.1767
5.9853
6.5429
6.1785
6.5081
6.3084
6.1033
6.1054
6.5776
6.4817
6.6337

127

Depth [m]
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
200
200
200
0
0
0
200
100

Sensor
accelerometer
accelerometer
broadband
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
broadband
accelerometer
accelerometer
geophone
geophone
geophone
accelerometer
accelerometer
accelerometer
geophone
geophone

Hardrock
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

KNMI

NL.EVK..HGZ
NL.FDKD..HGZ
NL.FR014..HHZ
NL.FR024..HHZ
NL.FR034..HHZ
NL.FR044..HHZ
NL.FR054..HHZ
NL.FR064..HHZ
NL.FSW5.01.HHZ
NL.G014..HHZ
NL.G024..HHZ
NL.G034..HHZ
NL.G044..HHZ
NL.G054..HHZ
NL.G064..HHZ
NL.G074..HHZ
NL.G084..HHZ
NL.G094..HHZ
NL.G104..HHZ
NL.G114..HHZ
NL.G124..HHZ
NL.G134..HHZ
NL.G144..HHZ
NL.G164..HHZ
NL.G174..HHZ
NL.G184..HHZ
NL.G194..HHZ
NL.G204..HHZ
NL.G214..HHZ
NL.G224..HHZ
NL.G234..HHZ
NL.G244..HHZ
NL.G254..HHZ
NL.G264..HHZ
NL.G274..HHZ
NL.G284..HHZ
NL.G294..HHZ
NL.G304..HHZ
NL.G314..HHZ
NL.G324..HHZ
NL.G334..HHZ
NL.G344..HHZ
NL.G354..HHZ
NL.G364..HHZ
NL.G374..HHZ
NL.G384..HHZ
NL.G394..HHZ
NL.G404..HHZ
NL.G414..HHZ
NL.G424..HHZ
NL.G434..HHZ
NL.G444..HHZ
NL.G454..HHZ

53.0797
52.8483
53.1424
52.8559
52.9815
52.8761
53.0705
53.0007
53.2135
53.4425
53.4490
53.4159
53.4147
53.4184
53.4111
53.3946
53.3954
53.3878
53.3785
53.3669
53.3699
53.3489
53.3586
53.3385
53.3362
53.3306
53.3378
53.3353
53.3134
53.3103
53.3095
53.2900
53.3081
53.2821
53.2832
53.2722
53.2814
53.2792
53.2962
53.2498
53.2499
53.2525
53.2549
53.2651
53.2740
53.2005
53.2253
53.2250
53.2234
53.2217
53.2319
53.2042
53.1995

6.8116
6.1788
5.8868
5.9304
5.7126
6.1187
5.8685
6.0237
7.1195
6.7312
6.8058
6.5976
6.6747
6.7713
6.8709
6.5458
6.6439
6.7245
6.8044
6.8891
6.5714
6.6686
6.7708
6.5314
6.6363
6.7217
6.8170
6.8952
6.5937
6.6776
6.7675
6.8562
6.9433
6.5458
6.6281
6.7292
6.7870
6.8951
7.0218
6.5792
6.6708
6.7651
6.8621
6.9245
7.0350
6.5476
6.7241
6.8099
6.8919
6.9886
7.0791
6.6840
6.7645

128

0
0
200
200
200
200
200
200
300
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200

accelerometer
accelerometer
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone

0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

KNMI

NL.G464..HHZ
NL.G474..HHZ
NL.G484..HHZ
NL.G494..HHZ
NL.G504..HHZ
NL.G514..HHZ
NL.G524..HHZ
NL.G534..HHZ
NL.G544..HHZ
NL.G554..HHZ
NL.G564..HHZ
NL.G574..HHZ
NL.G584..HHZ
NL.G594..HHZ
NL.G604..HHZ
NL.G614..HHZ
NL.G624..HHZ
NL.G634..HHZ
NL.G644..HHZ
NL.G654..HHZ
NL.G664..HHZ
NL.G674..HHZ
NL.G684..HHZ
NL.G694..HHZ
NL.G704..HHZ
NL.G710..HGZ
NL.G720..HGZ
NL.G730..HGZ
NL.G740..HGZ
NL.G750..HGZ
NL.G760..HGZ
NL.G770..HGZ
NL.G780..HGZ
NL.G790..HGZ
NL.G800..HGZ
NL.G81B..HHZ
NL.G82B..HHZ
NL.G83B..HHZ
NL.G84B..HHZ
NL.GK014..HHZ
NL.GK024..HHZ
NL.GK034..HHZ
NL.GK044..HHZ
NL.GRP..HGZ
NL.HGN.02.BHZ
NL.HLB.11.EPZ
NL.HLG.11.EPZ
NL.HLH.11.EPZ
NL.HLK.11.EPZ
NL.HOEV..HGZ
NL.HRDB..HGZ
NL.HRKB..BHZ
NL.HWF4.01.HHZ

53.1978
53.1964
53.1973
53.1765
53.1745
53.1732
53.1699
53.1674
53.1243
53.1433
53.1484
53.1351
53.1103
53.1160
53.2979
53.3819
53.3870
53.3314
53.2207
53.1692
53.1794
53.3213
53.2351
53.2128
53.2407
53.1483
53.2462
53.2719
53.3043
53.3078
53.3507
53.3403
53.3648
53.3836
53.4154
53.3103
53.2253
53.1484
53.2207
53.2906
53.3010
53.2827
53.2462
53.1788
50.7640
53.1465
53.1491
53.1404
53.1450
52.9025
52.5770
51.1879
53.0710

6.8498
6.9413
7.0290
6.7306
6.8093
6.8848
6.9837
7.0801
6.7580
6.8461
6.9309
7.0222
6.9071
6.9877
7.0752
6.6813
6.8419
6.9109
7.0380
6.9543
6.5852
6.8228
6.6347
6.6438
6.5379
6.6507
6.4792
6.4301
6.3894
6.4875
6.3901
6.3304
6.4894
6.3319
6.4876
6.6776
6.7241
6.9309
7.0380
6.2471
6.2927
6.3056
6.3261
5.9779
5.9317
7.0036
7.0071
7.0065
7.0096
6.1055
6.6259
6.1678
6.3512

129

200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
200
0
0
0
0
0
0
0
0
0
0
100
100
100
100
200
200
200
200
0
4
60
58
60
60
0
0
0
200

geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
accelerometer
broadband
geophone
geophone
geophone
geophone
accelerometer
accelerometer
broadband
geophone

0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0

KNMI

NL.J01..HHZ
NL.J02.01.HHZ
NL.KLWL..HGZ
NL.KMP..HGZ
NL.LUTT..HGZ
NL.MAME..BHZ
NL.MAR..HGZ
NL.N014..HHZ
NL.N024..HHZ
NL.N034..HHZ
NL.NGA..HGZ
NL.NH014..HHZ
NL.NIW4.01.HHZ
NL.NLDW..HGZ
NL.ODBK..HGZ
NL.ODHP..HGZ
NL.OLK..HGZ
NL.OOTH..HGZ
NL.OPLO.01.BHZ
NL.OSTM..HGZ
NL.OTL2.01.HHZ
NL.PET..HGZ
NL.PMR..HGZ
NL.PPB3.01.HHZ
NL.ROLD..BHZ
NL.SL014..HHZ
NL.SNB..HGZ
NL.SPY4.01.HHZ
NL.SRHV..HGZ
NL.STR..HGZ
NL.SUH4.01.HHZ
NL.T014..HHZ
NL.T024..HHZ
NL.T034..HHZ
NL.T044..HHZ
NL.T054..HHZ
NL.T064..HHZ
NL.T084..HHZ
NL.TERZ.01.BHZ
NL.TWZL..HGZ
NL.UTRP..HGZ
NL.VBG4.01.HHZ
NL.VDM24..HHZ
NL.VHP..HGZ
NL.VKB..BHZ
NL.VLD..HGZ
NL.VLW3.01.HHZ
NL.VNDM4..HHZ
NL.VNWD..HGZ
NL.VPR..HGZ
NL.VRS..HGZ
NL.WDB4.01.HHZ
NL.WMH4.01.HHZ

52.7650
52.7255
51.7755
53.3206
52.6133
50.8000
53.1540
53.1185
53.0720
53.1110
53.1364
52.5657
53.3548
52.0050
52.9405
52.8852
52.6293
52.7295
51.5888
53.1875
52.6289
52.7706
52.5116
52.6526
50.8694
52.8551
52.6842
53.4098
53.1833
51.7437
53.2113
52.4643
52.4634
52.4297
52.4375
52.3806
52.3432
52.3893
50.7568
53.2325
53.0969
52.5440
53.0844
53.0957
50.8669
52.8668
52.9682
53.1151
53.2368
51.8915
53.0360
53.2083
52.7096

5.0545
5.1948
4.4406
6.1920
6.5715
5.9727
6.2631
6.4733
6.4530
6.4045
6.0257
4.9965
6.0430
4.1957
6.1339
6.0249
4.8495
6.7534
5.8121
6.0857
4.8227
4.6595
5.0106
4.6700
6.0847
4.8683
6.8719
6.7838
6.1849
4.5573
6.2110
6.7401
6.9787
6.7918
6.8905
6.9057
6.9097
6.9965
5.9061
6.0927
6.1842
6.6693
6.8336
5.9553
5.7847
6.2157
7.0972
6.8390
5.9984
4.3715
6.6046
6.7355
4.7498

130

0
0
0
0
0
0
0
200
200
200
0
200
120
0
0
0
0
0
0
0
100
0
0
165
0
200
0
120
0
0
120
200
200
200
200
200
200
200
250
0
0
192
200
0
0
0
150
200
0
0
0
197
200

geophone
geophone
accelerometer
accelerometer
accelerometer
broadband
accelerometer
geophone
geophone
geophone
accelerometer
geophone
geophone
accelerometer
accelerometer
accelerometer
accelerometer
accelerometer
broadband
accelerometer
geophone
accelerometer
accelerometer
geophone
broadband
geophone
accelerometer
geophone
accelerometer
accelerometer
geophone
geophone
geophone
geophone
geophone
geophone
geophone
geophone
broadband
accelerometer
accelerometer
geophone
geophone
accelerometer
broadband
accelerometer
geophone
geophone
accelerometer
accelerometer
accelerometer
geophone
geophone

0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
1
0
0
0
0
0
0
0
0

KNMI

NL.WSVN..HGZ
NL.WTSB.01.BHZ
NL.WYN4.00.HHZ
NL.ZDL..HGZ
NL.ZH014..HHZ
NL.ZH024..HHZ
NL.ZH034..HHZ
NL.ZH044..HHZ
NL.ZH064..HHZ
NL.ZH074..HHZ
NL.ZLV4.01.HHZ
GE.IBBN..HHZ
GR.AHRW..HHZ
GR.BUG..HHZ
GR.RAST..EHZ
NR.NE116..BHZ
BE.BOST..HHZ
BE.MEM..HHZ
BE.BEBN..HHZ
BE.OPTB..HHZ
BE.RCHB..HHZ
BE.UCC..HHZ

52.8204
51.9663
53.1995
53.0791
51.9137
52.0744
51.9783
51.8417
51.9781
52.1037
53.0921
52.3063
50.5410
51.4455
52.9096
51.5042
51.2382
50.6087
50.7970
51.1115
50.1552
50.7972

6.1983
6.7989
5.4578
6.7026
4.2881
4.4884
4.4001
4.3626
4.2204
4.2946
6.7533
7.7592
7.0760
7.2643
7.7164
4.9209
2.9387
6.0096
5.6778
5.6360
5.2268
4.3604

131

0
0
120
0
200
200
200
200
200
200
200
1
0
0
0
0
304
0
0
379
43
141

accelerometer
broadband
geophone
accelerometer
geophone
geophone
geophone
geophone
geophone
geophone
geophone
broadband
broadband
broadband
geophone
broadband
broadband
broadband
broadband
broadband
broadband
broadband

0
1
0
0
0
0
0
0
0
0
0
1
1
1
0
0
1
1
1
1
1
1

D 2022 NL-wide maps with locationuncertainty parameters
D.1

Largest-horizontal-location-uncertainty (σ1 ) maps

132

D.1. LARGEST-HORIZONTAL-LOCATION-UNCERTAINTY (σ1 ) MAPS

KNMI

Figure D.1: Spatial distribution of σ1 in 2022 for M=0.5 (color map). White means that
events with M=0.5 can likely not be detected. Surface sensors are shown with squares,
deep geophones with orange triangles and broadbands with green triangles. Thin black lines
denote province borders within the Netherlands. Thick black lines are country borders. The
coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in kilometers. Dark red
depicts uncertainties of 2800 m and higher.

133

D.1. LARGEST-HORIZONTAL-LOCATION-UNCERTAINTY (σ1 ) MAPS

KNMI

Figure D.2: Spatial distribution of σ1 in 2022 for M=1.0 (color map). White means that
events with M=1.0 can likely not be detected. Surface sensors are shown with squares,
deep geophones with orange triangles and broadbands with green triangles. Thin black lines
denote province borders within the Netherlands. Thick black lines are country borders. The
coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

134

D.1. LARGEST-HORIZONTAL-LOCATION-UNCERTAINTY (σ1 ) MAPS

KNMI

Figure D.3: Spatial distribution of σ1 in 2022 for M=1.5 (color map). White means that
events with M=1.5 can likely not be detected. Surface sensors are shown with squares,
deep geophones with orange triangles and broadbands with green triangles. Thin black lines
denote province borders within the Netherlands. Thick black lines are country borders. The
coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

135

D.1. LARGEST-HORIZONTAL-LOCATION-UNCERTAINTY (σ1 ) MAPS

KNMI

Figure D.4: Spatial distribution of σ1 in 2022 for M=2.0 (color map). White means that
events with M=2.0 can likely not be detected. Surface sensors are shown with squares,
deep geophones with orange triangles and broadbands with green triangles. Thin black lines
denote province borders within the Netherlands. Thick black lines are country borders. The
coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

136

D.1. LARGEST-HORIZONTAL-LOCATION-UNCERTAINTY (σ1 ) MAPS

KNMI

Figure D.5: Spatial distribution of σ1 in 2022 for M=2.5 (color map). White means that
events with M=2.5 can likely not be detected. Surface sensors are shown with squares,
deep geophones with orange triangles and broadbands with green triangles. Thin black lines
denote province borders within the Netherlands. Thick black lines are country borders. The
coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

137

D.1. LARGEST-HORIZONTAL-LOCATION-UNCERTAINTY (σ1 ) MAPS

KNMI

Figure D.6: Spatial distribution of σ1 in 2022 for M=3.0 (color map). Surface sensors
are shown with squares, deep geophones with orange triangles and broadbands with green
triangles. Thin black lines denote province borders within the Netherlands. Thick black lines
are country borders. The coordinates are in the local rectangular grid (Rijksdriehoekstelsel)
in kilometers.

138

D.1. LARGEST-HORIZONTAL-LOCATION-UNCERTAINTY (σ1 ) MAPS

KNMI

Figure D.7: Spatial distribution of σ1 in 2022 for M=3.5 (color map). Surface sensors
are shown with squares, deep geophones with orange triangles and broadbands with green
triangles. Thin black lines denote province borders within the Netherlands. Thick black lines
are country borders. The coordinates are in the local rectangular grid (Rijksdriehoekstelsel)
in kilometers.

139

D.1. LARGEST-HORIZONTAL-LOCATION-UNCERTAINTY (σ1 ) MAPS

KNMI

Figure D.8: Spatial distribution of σ1 in 2022 for M=4.0 (color map). Surface sensors
are shown with squares, deep geophones with orange triangles and broadbands with green
triangles. Thin black lines denote province borders within the Netherlands. Thick black lines
are country borders. The coordinates are in the local rectangular grid (Rijksdriehoekstelsel)
in kilometers.

140

D.1. LARGEST-HORIZONTAL-LOCATION-UNCERTAINTY (σ1 ) MAPS

141

KNMI

D.2. VERTICAL-LOCATION-UNCERTAINTY (σZ ) MAPS

D.2

KNMI

Vertical-location-uncertainty (σZ ) maps

Figure D.9: Spatial distribution of σZ in 2022 for M=0.5 (color map). Surface sensors
are shown with squares, deep geophones with orange triangles and broadbands with green
triangles. White means that events with M=0.5 can likely not be detected. Thin black lines
denote province borders within the Netherlands. Thick black lines are country borders. The
coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in kilometers. Dark red
depicts uncertainties of 3500 m and higher.
142

D.2. VERTICAL-LOCATION-UNCERTAINTY (σZ ) MAPS

KNMI

Figure D.10: Spatial distribution of σZ in 2022 for M=1.0 (color map). Surface sensors
are shown with squares, deep geophones with orange triangles and broadbands with green
triangles. White means that events with M=1.0 can likely not be detected. Uncertainties
of 3500 m and higher are shown in red. Thin black lines denote province borders within
the Netherlands. Thick black lines are country borders. The coordinates are in the local
rectangular grid (Rijksdriehoekstelsel) in kilometers.

143

D.2. VERTICAL-LOCATION-UNCERTAINTY (σZ ) MAPS

KNMI

Figure D.11: Spatial distribution of σZ in 2022 for M=1.5 (color map). Surface sensors
are shown with squares, deep geophones with orange triangles and broadbands with green
triangles. White means that events with M=1.5 can likely not be detected. Thin black lines
denote province borders within the Netherlands. Thick black lines are country borders. The
coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

144

D.2. VERTICAL-LOCATION-UNCERTAINTY (σZ ) MAPS

KNMI

Figure D.12: Spatial distribution of σZ in 2022 for M=2.0 (color map). Surface sensors
are shown with squares, deep geophones with orange triangles and broadbands with green
triangles. White means that events with M=2.0 can likely not be detected. Thin black lines
denote province borders within the Netherlands. Thick black lines are country borders. The
coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

145

D.2. VERTICAL-LOCATION-UNCERTAINTY (σZ ) MAPS

KNMI

Figure D.13: Spatial distribution of σZ in 2022 for M=2.5 (color map). Surface sensors
are shown with squares, deep geophones with orange triangles and broadbands with green
triangles. White means that events with M=2.5 can likely not be detected. Thin black lines
denote province borders within the Netherlands. Thick black lines are country borders. The
coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

146

D.2. VERTICAL-LOCATION-UNCERTAINTY (σZ ) MAPS

KNMI

Figure D.14: Spatial distribution of σZ in 2022 for M=3.0 (color map). Surface sensors
are shown with squares, deep geophones with orange triangles and broadbands with green
triangles. Thin black lines denote province borders within the Netherlands. Thick black lines
are country borders. The coordinates are in the local rectangular grid (Rijksdriehoekstelsel)
in kilometers.

147

D.2. VERTICAL-LOCATION-UNCERTAINTY (σZ ) MAPS

KNMI

Figure D.15: Spatial distribution of σZ in 2022 for M=3.5 (color map). Surface sensors
are shown with squares, deep geophones with orange triangles and broadbands with green
triangles. Thin black lines denote province borders within the Netherlands. Thick black lines
are country borders. The coordinates are in the local rectangular grid (Rijksdriehoekstelsel)
in kilometers.

148

D.2. VERTICAL-LOCATION-UNCERTAINTY (σZ ) MAPS

KNMI

Figure D.16: Spatial distribution of σZ in 2022 for M=4.0 (color map). Surface sensors
are shown with squares, deep geophones with orange triangles and broadbands with green
triangles. Thin black lines denote province borders within the Netherlands. Thick black lines
are country borders. The coordinates are in the local rectangular grid (Rijksdriehoekstelsel)
in kilometers.

149

D.3. AZIMUTHAL-GAP (ψ) MAPS

D.3

KNMI

Azimuthal-gap (ψ) maps

Figure D.17: Spatial distribution of azimuthal gap ψ in 2022 for M=0.5 (color map). White
means that events with M=0.5 can likely not be detected. Thin black lines denote province
borders within the Netherlands. Thick black lines are country borders. The coordinates are
in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

150

D.3. AZIMUTHAL-GAP (ψ) MAPS

KNMI

Figure D.18: Spatial distribution of azimuthal gap ψ in 2022 for M=1.0 (color map). White
means that events with M=1.0 can likely not be detected. Thin black lines denote province
borders within the Netherlands. Thick black lines are country borders. The coordinates are
in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

151

D.3. AZIMUTHAL-GAP (ψ) MAPS

KNMI

Figure D.19: Spatial distribution of azimuthal gap ψ in 2022 for M=1.5 (color map). White
means that events with M=1.5 can likely not be detected. Thin black lines denote province
borders within the Netherlands. Thick black lines are country borders. The coordinates are
in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

152

D.3. AZIMUTHAL-GAP (ψ) MAPS

KNMI

Figure D.20: Spatial distribution of azimuthal gap ψ in 2022 for M=2.0 (color map). White
means that events with M=2.0 can likely not be detected. Thin black lines denote province
borders within the Netherlands. Thick black lines are country borders. The coordinates are
in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

153

D.3. AZIMUTHAL-GAP (ψ) MAPS

KNMI

Figure D.21: Spatial distribution of azimuthal gap ψ in 2022 for M=2.5 (color map). White
means that events with M=2.5 can likely not be detected. Thin black lines denote province
borders within the Netherlands. Thick black lines are country borders. The coordinates are
in the local rectangular grid (Rijksdriehoekstelsel) in kilometers.

154

D.3. AZIMUTHAL-GAP (ψ) MAPS

KNMI

Figure D.22: Spatial distribution of azimuthal gap ψ in 2022 for M=3.0 (color map).
Thin black lines denote province borders within the Netherlands. Thick black lines are
country borders. The coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in
kilometers.

155

D.3. AZIMUTHAL-GAP (ψ) MAPS

KNMI

Figure D.23: Spatial distribution of azimuthal gap ψ in 2022 for M=3.5 (color map).
Thin black lines denote province borders within the Netherlands. Thick black lines are
country borders. The coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in
kilometers.

156

D.3. AZIMUTHAL-GAP (ψ) MAPS

KNMI

Figure D.24: Spatial distribution of azimuthal gap ψ in 2022 for M=4.0 (color map).
Thin black lines denote province borders within the Netherlands. Thick black lines are
country borders. The coordinates are in the local rectangular grid (Rijksdriehoekstelsel) in
kilometers.

157