KEM-24: Effect of pressure maintenance by
fluid injection on seismic hazard, Hazard
analysis
Report | KEM-24 Report C
172147-REP01-FNV_SHA_KEM24 03 | Fugro NL Land B.V.

Fugro NL Land B.V.

Document Control
Document Information
Project Title

KEM-24 “Effect of pressure maintenance by fluid injection on seismic hazard”

Document Title

KEM-24: Effect of pressure maintenance by fluid injection on seismic hazard, Hazard analysis

Fugro Project No.

172147

Fugro Document No.

172147-REP01-FNV_SHA_KEM24

Issue Number

03

Issue Status

Final Report

Fugro Legal Entity

Fugro France SAS

Issuing Office Address

115 Avenue de la Capelado, 34160 CASTRIES - France

Client Information
Client

Client: Ministry of Economic Affairs and Climate Policy

Client Address

Procurement Office, P.O. Box 93144, 2509 AC THE HAGUE

Revision History
Issue

Date

Status

Comments on Content

Prepared
By

Checked
By

Approved
By

01

15 November
2021

For Review

Awaiting client comments

MLA, JMA,
CDU, RSE

BWH

RSE

02

20 May 2022

For Review

Awaiting client comments

MLA, JMA,
CDU, RSE

BWH

RSE

03

30 Sep 2022

Final

MLA, JMA,
CDU, RSE

BWH

RSE

Project Team
Initials

Name

Role

CDU

Cedric DUVAIL

Senior Geologist, Project Manager

MLA

Manel LABIDI

Geologist

JMA

Judith MARINIERE

Seismologist

RSE

Ramon SECANELL

Senior Seismologist

BWH

Beau WHITNEY

Senior Geologist

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic hazard,
Hazard analysis

Fugro NL Land B.V.

Executive Summary
Fugro conducted an induced seismic hazard assessment of the Groningen gas field (Figure 1.1)
taking into account the intermediate results of Dynafrax UG report for KEM-24 on the effect of
pressure maintenance by fluid injection on seismic risk.
The goal was to establish, using synthetic seismic catalogs generated by Dynafrax using different
injection gas hypothesis, which final fluid injection scenarios should be considered for the
purposes of reservoir pressure maintenance and lowering the seismic risk profile.
Two site-specific probabilistic seismic hazard analyses (PSHA) were performed, using a
methodology consistent with current international practice for industrial facilities and in
agreement with EC8 recommendations. One PSHA considering only seismotectonic sources and
one including induced seismic sources.
The hazard model used includes:
◼

Two (2) seismotectonic models;

◼

Five (5) Ground-Motion Prediction Equations (GMPEs);

The soil conditions considered correspond to standard rock with a VS30 = 200 m/s, which
correspond to generic soil conditions in the Netherlands, where the standard rock conditions
(Vs30=800 m/s) are not usual and soft soil conditions are often found.
Note that the calculations in this report are based on an existing probabilistic model (catalogue
and seismotectonic models) used for an Offshore Windfarm project, in Netherlands, report
P904711/SHA 01 (FUGRO, 2020). The regional model used corresponds to that probabilistic model
and then, we included the induced seismicity produced by the different injection scenarios
considered.
We considered the 4 different scenarios provided by Dynafrax (arbitrary hypothetical scenarios
selected), associated with different injection rates. The four scenarios considered (see Figure 6.1)
were:
◼

Scenario 1: 60 years of depletion from 1960 to 2020, 10 years of injection (up to
2030) and 20 years of shut-in (up to 2050)

◼

Scenario 2: 60 years of depletion from 1960 to 2020, 1 years of injection (up to
2022) and 29 years of shut-in (up to 2050)

◼

Scenario 3: 60 years of depletion from 1960 to 2020, 30 years without activity
(called reference period), 10 years of injection (up to 2060) and 20 years of shut-in
(up to 2080)

◼

Scenario 3: 60 years of depletion from 1960 to 2020, 30 years without activity
(called reference period), 5 years of injection (up to 2055) and 25 years of shut-in
(up to 2080)

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

For each one of these scenarios and for two periods of time (injection time only or injection with
a shut-in period), a seismic occurrence distribution (Gutenberg-Richter law) was determined and
the seismic hazard, taking into account the new induced zone (corresponding to Groningen
reservoir) in the seismotectonic model and its seismic occurrence distribution, was estimated.
Conclusion on induced seismicity
Independently of the scenario used, the seismic occurrence distribution of the induced zone
(Groningen reservoir) showed higher seismic activity levels than the natural seismotectonic zone
if we consider a period equivalent to the injection time only. This means that the seismic activity
rate (ind(M)/km2 value) of the induced zone is always higher than the seismic activity rate
(nat(M)/km2) of the natural seismotectonic zone where Groningen is situated. However, we note
that the natural seismic activity of the Groningen region is very low. So that even a low induced
seismicity contribution can become significant compared to the background seismic activity. We
should also note that if we consider a longer period of 30 years (such as injection time + shut-in),
scenarios 1, 2 and 4 will have similar seismic activity levels than the natural seismotectonic zone.
As a result of that, considering a period equivalent to the injection time only, we found that the
introduction of the induced seismicity in the seismic hazard model, independently of the scenario,
increases in a non-negligible way the seismic hazard of the Groningen region. However, we note
that:
1) Despite the non-negligible increase in the seismic hazard, it remains very low, less than
0.04 g for the worst case scenario at 475 years return period. However, we have to note
that the consideration of a 475 years return period assumes the hypothesis of stationarity
along the time, which would be a very conservative hypothesis. This value is significantly
lower than the PGA predicted in past TNO report, where the PGA was fixed at 0.117g (see
chapter 2.5). Very likely the reason is the low seismic activity predicted by the modelling
used here, compared with the really observed seismicity used in TNO report (see chapter
2.5)
2) The injection induced change in seismic hazard should only occur for the 1, 5 or 10 years
of injection simulation. After these 1, 5 or 10 years (depending on the scenario chosen) the
seismic hazard returns to its natural baseline, without any induced seismicity consideration
(with mean PGA around 0.015 g).
3) The induced seismicity only produces small earthquakes (the Mmax<4.0 based on the
simulations performed) and the associated seismic accelerations from these events would
be lower than 0.1 g, except in extreme cases (i.e., considering the median acceleration
predicted by the GMPE plus 3 standard deviations).
4) When the injection of fluid is stopped, the simulations indicate that the induced seismicity
diminishes to close to zero.
172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

5) Based on the seismic hazard associated to each scenario, the best scenario would be to
follow the procedure of Scenario 1 (considering an injection period of 10 years) that leads
to an increase of only 30 % at 475 years return period and can even be reduced to 11 %
considering a period of 30 years (injection time + shut-in). The increase of seismic activity
rates would be lower with longer periods of time. In any case, the introduction of injection
increases the seismic hazard for the 4 scenarios tested.
Therefore, we can conclude that, from the seismic hazard point of view the “best action” for
induced seismicity is to allow the system to re-equilibrate pressures in a natural way. The
injection of fluids, considering several different scenarios tested, leads to an increase in the
seismicity rate and, therefore, to an increase of the seismic hazard. This increase is temporal being
controlled by the injection duration and affects predominantly small ground accelerations.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

Contents
1.

Introduction

8

1.1

Scope and purpose

10

1.2

Organisation of the report

12

2.

Review of previous seismic hazard assessment studies

13

2.1

Seismic zoning map conforming to Eurocode 8 for the Netherlands (Crook, 1996)

13

2.2

Global Seismic Hazard Assessment Program

15

2.3

The 2013 European Seismic Hazard Model: Key Components and Results (SHARE Project)16

2.4

Global Earthquake Model

17

2.5

TNO Model Chain

18

3.

Methodology

21

3.1

Methods for determining ground motions

21

3.2

Probabilistic approach

21

4.

Seismological database

24

4.1

Earthquake catalogue with seismotectonic origin

24

4.1.1
4.1.2
4.1.3

24
25
25

Data compilation
Final compiled declustered catalogue
Completeness analysis

4.2

Induced earthquake catalogue

27

5.

Seismic sources

31

5.1

Seismotectonic Model

31

5.2

Induced Model

34

6.

Seismic activity model

35

6.1

Presentation

35

6.2

Gutenberg-Richter definition

36

6.2.1
6.2.2
6.2.3

36
36
39

Seismotectonic Model 1
Seismotectonic Model 2
Induced Model

7.

Ground Motion Prediction Equations

43

7.1

Selection of GMPEs for seismotectonic sources

43

7.2

Selection of GMPEs for induced seismicity sources

43

8.

Modelling Induced Seismicity

46

8.1

Global Probabilistic Logic Tree

46

8.2

Common hypotheses

48

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

8.3

Probabilistic Seismic Hazard Results

48

8.3.1
8.3.2
8.3.3
8.3.4

49
54
54
55

Seismic Hazard Curve at PGA for Groningen Gas Field
Sensitivity analysis for Mmax
Sensitivity analysis for GMPEs
Sensitivity analysis for Mmin

9.

Conclusion

57

10.

References

59

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

1.

Introduction
Fugro conducted an induced seismic hazard and risk assessment for the Groningen gas
field (Figure 1.1) following the induced seismic scenarios defined by Dynafrax UG for the
KEM-24 project. The scenarios tested were defined to test the effect of pressure
maintenance by fluid injection on seismic risk.
The goal was to establish, using numerical/synthetic seismicity catalogues generated by
Dynafrax, the effect on seismic hazard of several fluid injection scenarios. To achieve this,
a PSHA was performed that considered natural and induced seismic sources. For the
PSHA, one calculation point was used for simplicity and comparison purposes. It
corresponds to the central point of the Groningen reservoir, close to the city of
Loppersum.
Coordinates of the calculation point used for the PSHA are specified in Table 1.1.
Table 1.1: Coordinates of the calculation point for the PSHA.
Coordinates (WGS 84)
Latitude (°)

Longitude (°)

53.328367°N

6.762381°E

The Groningen gas field (GGF) is located in an area of very low natural seismic activity
according to GSHAP seismic hazard maps (Giardini et al., 1999), with Peak Ground
Acceleration (PGA) values below 0.2 m/s2, for rock conditions at the 475-year return
period. The area of interest is far from moderate natural seismic sources of the
Netherlands which are situated in the south of country in the Roer Valley Graben (Limburg
area).

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

Figure 1.1: Location of the Groningen gas field.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

1.1

Scope and purpose
The objective of this study was the evaluation of induced seismicity using modelling of
different injection scenarios to determine its influence on the seismic hazard of the region.
To perform this evaluation, we performed 2 different calculations:
1) A seismic hazard assessment taking into account only natural seismotectonic sources
2) A seismic hazard assessment taking into account the seismotectonic sources and the
induced seismicity.
In this study, we were mainly interested in the comparison of these seismic hazard
calculations. Typically, for regions where the seismic activity is very high (i.e., in Indonesia,
Chile or other tectonically active areas), induced seismicity from human activities is
negligible for the local seismic hazard because the relatively low magnitudes produced
by induced sources are not significant compared to larger and more frequent natural
seismic activity.
However, in zones with less natural seismicity, the consideration of induced seismic
sources (i.e., case of Groningen) can produce a non-negligible increase to the seismic
hazard.
For this study, firstly, a probabilistic seismic hazard assessment (PSHA) at the Groningen
gas field that considers only natural seismicity was conducted. The probabilistic model
used was defined in a PSHA performed by Fugro for an Offshore Windfarm project, in
2020.
Uncertainties related to the delineation of seismic sources and ground motion models
are considered according to state-of-the-art methods and propagated in the probabilistic
approach through logic trees.
Then, a second PSHA was the performed considering the previous probabilistic model
(considering only natural seismotectonic seismicity) with the addition of a new induced
seismic source associated with injection activities in the Groningen field.
The induced seismicity was assumed to occur at 3 km depth, equivalent to average
Groningen reservoir depth. The VS30 value adopted was 200 m/s and one calculation point
was considered at the center of the Groningen reservoir area, close to the city of
Loppersum.
Then, we compared the ground motion values derived from each assessment. The
comparison is done using seismic hazard curves (annual exceedance rates of different
acceleration levels). Due to the time dependency (the injection will only happen for 10
years or less with shut-in periods from 20 to 29 years), we do not define response spectra
associated with longer return periods.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

We have to note that the seismic hazard assessment was performed using the current
industrial practice for seismic hazard assessment of critical facilities. It means that the
seismic activity rates are supposed to be stationary. However, we know that the induced
seismicity is not stationary. To take into account the non-stationarity, some temporal
slices were considered. In each slice, the induced seismicity is supposed to be stationary.
Therefore, the seismic hazard assessment here presented should be considered as a first
approximation of the reality. In Groningen, TNO developed more sophisticated
probabilistic models using nonstationary models. However, at the time of the realization
of this study the models were not public available and we decided the use of standard
models for seismic hazard assessment.
We have to note also that the central part of the KEM-24 project was the development of
fluid injection models and not a very detailed seismic hazard assessment (which was
already performed in Groningen). The main objective was the comparison between the
seismic hazard with and without consideration of the induced seismicity. To do this
comparison, we used one of the output data of the fluid injection modelling: the induced
seismic catalogues. Then, these seismic catalogues associated to different fluid injection
scenarios were used to perform the seismic hazard assessment considering induced
seismicity. Even if the seismic hazard assessment can be considered only as a first
approximation to the reality, we consider that the seismic hazard assessment performed
allows to identify if the fluid injection scenarios used could lead to significant increments
(or not) of the seismic hazard. Therefore, the seismic hazard assessment results, even if
they were defined using stationary models of seismicity, allow the identification of seismic
hazard increments and facilitate further decisions by authorities.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

1.2

Organisation of the report
The purpose of this report is to present the results of the probabilistic seismic hazard
analyses and a comparison of their results. It is organized as follows:
◼

Chapter 2 concerns the review and evaluation of previous seismic hazard studies
conducted by others in the region. It provides background information on the ground
motion estimates of Groningen area and provides a means to compare the PSHA
results obtained in this study. As the objective is the comparison between the seismic
hazard with and without induced seismicity, the review of seismic hazard assessments
considering only natural seismicity demonstrates the range of values (or uncertainty)
achieved by others. This provides a frame of reference for the impact of including
induced seismicity in the assessment;

◼

Chapter 3 describes the methodology for determining the seismic ground motions
in both seismic hazard assessment performed (with and without induced sources);

◼

Chapter 4 presents the declustered earthquake catalogue used;

◼

Chapter 5 describes the alternative seismotectonic models used;

◼

Chapter 6 presents the methodology to derive:

◼

•

the seismic activity parameters for all the zones of each seismotectonic model;

•

the seismic activity parameters of the induced zone;

Chapter 7 presents the ground motion prediction equations (GMPEs) used in this
study;

◼

Chapter 8 presents the description of the logic tree, global parameters for PSHA
calculation and the computed ground motions for the natural and induced seismicity
models;

◼

Chapter 9 provides the conclusions of this study;

◼

Chapter 10 includes the references used for the study.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

2.

Review of previous seismic hazard assessment
studies
PSHA results, considering only natural seismicity, at the Groningen reservoir are quite
similar for all of the publications analysed. The region where the gas field is situated
corresponds to a low natural seismicity region, a so called ”Stable Continental Region”
(SCR) as opposed to “Active Crustal Region“.
In terms of PGA, the results previously published generally range from 0.01 g (0.1
m/s2) to 0.02 g (0.2 m/s2), for rock conditions at a return period of 475 years and
considering only natural seismicity.

2.1

Seismic zoning map conforming to Eurocode 8 for the Netherlands
(Crook, 1996)
References:
CROOK, Th. de, (1996), A seismic zoning map conforming to Eurocode 8, and practical
earthquake parameter relations for the Netherlands, Geologie en Mijnbouw, 75, pp 1118.
BROUWER J.W.R. (2010). The meaning of Eurocode 8 and Induced Seismicity For Earthquake
Engineering in the Netherlands. Fifth International Conference on Recent Advances in
Geotechnical Earthquake Engineering and Soil Dynamics, May 24-29, 2010, San Diego,
California.
De Crook (1996) performed the latest seismic hazard study for the Netherlands,
considering natural seismicity, based on the earthquake catalogue up to 1993. Figure 2.1
and Figure 2.2 show the maps of seismic hazard zones in the Netherlands based on this
study. These maps are used as the current hazard zonation map in the Netherlands. An
update of this study using a revised and extended catalogue of earthquakes and new
ground motion prediction relations is being prepared at the KNMI (Koninklijk Nederlands
Meteorologisch Instituut). Ground motions are used instead of intensities, which is more
convenient for engineering purposes and can be compared with actual measurements
(Brouwer, 2010).
The seismic zonation map of the Netherlands is based on a seismic hazard study with a
10 % probability of exceedance in 50 years (return period of 475 years).
Four seismic zones are defined, where for stiff soil type, the PGAs are respectively:
◼

Seismic zone A: PGA = 0.10 m/s² (0.01 g);

◼

Seismic zone B: PGA = 0.22 m/s² (0.022 g);

◼

Seismic zone C: PGA = 0.50 m/s² (0.05 g);

◼

Seismic zone D: PGA = 1.00 m/s² (0.1 g).

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

On this seismic hazard map, the Groningen Gas Field is located in a region associated to
an estimated PGA value lower than 0.1 m/s2 (zone A) along the northern coast of the
Netherlands, at the return period of 475 years (Figure 2.2). This seismic level
corresponds to a very low seismic hazard.

Figure 2.1: Currently available seismic hazard
zonation in the Netherlands in terms of intensity
(Crook, 1996).

Figure 2.2: Seismic zonation map of the Netherlands
based on a seismic hazard study with a 10% of exceedance
in 50 years (return period 475 years).

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

2.2

Global Seismic Hazard Assessment Program
References:
GRÜNTHAL, G., BOSSE, C., SELLAMI, S., MAYER-ROSA, D. & GIARDINI, D. (1999). Compilation
of the GSHAP regional seismic hazard for Europe, Africa and the Middle East.
www.seismo.ethz.ch/static/GSHAP/index.html
A seismic hazard map was published by GSHAP (Global Seismic Hazard Assessment
Program) for the region (Grünthal et al., 1999). The objective of the GSHAP program was
to propose a global seismic hazard map at the world scale, based on homogenized and
coordinated probabilistic approach.
Figure 2.3 is an extract of the global GSHAP hazard map focusing on the project region.
This map shows the PGA expected with 10% probability of exceedance in 50 years.
The area of interest of the Groningen Gas Field corresponds on the GSHAP map to a
region where estimated PGA values are below 0.2 m/s2 along the northern coast of
the Netherlands, at the return period of 475 years (Figure 2.3). This seismic level
corresponds to a very low seismic hazard, considering only natural seismicity.

Figure 2.3: Seismic hazard map for a return period of 475 years according to GSHAP (Grünthal et al., 1999)
(data and plotting tool from http://gmo.gfz-Potsdam.de/) and location of the Groningen Gas Field.
172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

2.3

The 2013 European Seismic Hazard Model: Key Components and Results
(SHARE Project)
References:
GIARDINI, D., J. WOESSNER, L. DANCIU (2014) Mapping Europe’s Seismic Hazard. EOS,
95(29): p. 261-262.
WOESSNER, J., DANCIU L., D. GIARDINI and the SHARE consortium (2015), The 2013
European Seismic Hazard Model: key components and results, Bull. Earthq. Eng.,
doi:10.1007/s10518-015-9795-1.
The 2013 European Seismic Hazard Model (ESHM13) results from a community-based
probabilistic seismic hazard assessment supported by the EU-FP7 project “Seismic Hazard
Harmonization in Europe” (SHARE, 2009–2013). The ESHM13 is a consistent seismic
hazard model for Europe and Turkey which overcomes the limitation of national borders
and includes a thorough quantification of the uncertainties. It corresponds to the first
completed regional effort contributing to the “Global Earthquake Model” initiative.
Seismic hazard calculations were computed for a standard rock condition (i.e., VS30= 800
m/s, Pagani et al., 2014). Several maps at different return periods (i.e., 475, 975, 2 475 and
4 975 years) were produced.
In this publication, and at the return period of 475 years (10 % exceedance in 50
years), the PGA is lower than 0.01 g (0.1 m/s2) along the northern coast of the
Netherlands (Figure 2.4).

Figure 2.4: Extract from SHARE map showing the PGA at 475-year return period (10% chance of
exceedance in 50 years) and location of the Groningen Gas Field.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

2.4

Global Earthquake Model
References:
M. PAGANI, J. GARCIA-PELAEZ, R. GEE, K. JOHNSON, V. POGGI, R. STYRON, G. WEATHERILL,
M. SIMIONATO, D. VIGANÒ, L. DANCIU, D. MONELLI (2018). Global Earthquake Model
(GEM) Seismic Hazard Map (version 2018.1 - December 2018), DOI: 10.13117/GEMGLOBAL-SEISMIC-HAZARD-MAP-2018.1
The Global Earthquake Model (GEM) is a non-profit foundation that works to develop and
disseminate global earthquake hazard information to inform the public and policy
makers. The most recent global seismic hazard map (version 2018.1) was created by
collating and aggregating a global database of regional and national probabilistic seismic
hazard models. These models are then used to calculate hazard values for the globe for
reference rock conditions (VS30 of 760-800 m/s).
Similar to the GSHAP model, the GEM hazard map presents an overview of hazard
variability at large scale.
In this publication, the PGA at the return period of 475 years is lower than 0.01 g
(0.1 m/s2) along the northern coast of the Netherlands, which is consistent with a
region of low seismicity (Figure 2.5).

Figure 2.5: Extract from GEM 2018 showing the 10% chance of exceedance in 50 years (PGA) and location of
the Groningen Gas Field.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

2.5

TNO Model Chain
References:
J.J BOMMER, B. EDWARDS, P.P. KRUIVER, A. RODRIGUEZ-MAREK, P.J. STAFFORD, M.
NTINALEXIS, E. RUIGROK, J. SPETZLER (2019). V6 Ground-Motion Model (GMM) for
Induced Seismicity in the Groningen Field with Assurance Letter. NAM report.
TNO (2020a). Probabilistic Seismic Hazard and Risk Analysis in the TNO Model Chain
Groningen, TNO 2020 R11052.
TNO (2020b). TNO Model Chain Groningen: Update and quick scan comparison of 2020
HRA model, TNO 2020 R11659.
TNO were commissioned by the Ministry of Economic Affairs and Climate of the
Netherlands (EZK) to develop a probabilistic seismic hazard and risk assessment (PSHRA)
procedure (i.e. a model chain) for application to target sites within the public domain
surrounding and within the Groningen gas field. The TNO Model Chain was developed to
replicate the hazard and risk assessment (HRA) models of the independent Groningen
operator NAM using an independent software implementation. The seismic hazard
component of the TNO Model Chain comprises of the seismic source model (SSM) and
the ground motion model (GMM).
The SSM forecasts the temporal-spatial distribution of the gas depletion induced
earthquakes and their magnitudes, conditional on the specified production scenario. The
use of the SSM implicitly assumes that the Groningen seismicity is treated as nonstationary, reflecting the non-Poisson temporal distribution of these earthquakes due to
resulting from gas depletion. Induced earthquakes in the Groningen gas field are
assumed to be caused by differential compaction on existing faults.
The SSM uses static data such as reservoir geometry, reservoir compressibility and fault
data to model the relationship between gas depletion, reservoir compaction and
subsurface faulting. Historic induced earthquake and field pore pressure changes are
required to train this compaction model to forecast future induced earthquakes. For a
given gas production scenario, the compaction model computes vertical strains based on
the associated pore pressure changes. These vertical strains are combined with the fault
data to compute spatial and temporal distributions of Coulomb stress changes. Main
shock seismicity distributions in time and space are then computed from the spatial and
temporal distributions of Coulomb stress changes. The SSM implements an EpidemicType Aftershock Sequence (ETAS) model to compute aftershock seismicity distribution
corresponding to the main shock seismicity distribution. The main shock and aftershock
seismicity distributions are combined to provide a total seismicity distribution. A b-value
is then assigned to these seismic events to define the Gutenberg-Richter magnitude
distribution. The magnitude distribution is bounded by a maximum probable magnitude.
The total seismicity distribution is now defined in time, space and magnitude. A rupture
model then translates the generated hypocenter locations into rupture plane distributions

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

with associated magnitudes. The final output of the SSM is therefore a statistical
distribution of total seismicity of a specific magnitude, at a specific distance and within a
specific year, for each point at the surface given a specific production scenario.
The TNO Model Chain GMM has recently been updated to implement GMM V6 of
Bommer et al. (2019). Given a specific magnitude and rupture distance, the Bommer et al.
(2019) GMM V6 provides statistical distributions of horizontal ground motions for 23
spectral periods, PGV and ground motion duration. Ground motions were computed for
bedrock corresponding to the base of the North Sea Group at a depth of approximately
800 m. A site-response zonation model defined within Bommer et al. (2018) is used to
translate the bedrock motions to the surface.
Within the seismic hazard component of the TNO Model Chain the SSM and GMM are
convolved to produce the seismic hazard estimates for a given production scenario in the
form of annual probabilities of exceedance for given intensities of spectral acceleration
for each grid point at the surface. These results can be presented as hazard curves, or as
hazard maps for a specific year following the onset of the production scenario. Hazard
curves and hazard maps are computed for specific years rather than for specific return
periods because non-stationary earthquakes are considered in the TNO Model Chain,
rather than stationary earthquakes as within conventional PSHA.
The TNO Model Chain reports provide a hazard map for PGA (treated as equivalent to T
= 0.01 s spectral accelerations) for the gas production year 2020/2021 assuming time zero
at the beginning of 2020 (Figure 2.6). Then, using the induced seismicity predicted during
the 2020/2021 period and assuming a return period of 475 years (implicitly, it means that
a stationarity seismicity is considered), the maximum PGA computed for this production
scenario (0.117 g) is stronger than our most conservative scenario (Scenario 2 with an
injection period of 1 year) giving 0.04 g. We have to note that the TNO model forecasts
7.86 events of magnitude 1.5 or above for one year. This annual activity rate is higher than
the annual activity rates predicted in any of the scenarios considered in this study.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

Figure 2.6: Extract from TNO (2020b) showing the PGA values in the Groningen field and the surrounding
regions observed for the 2020/2021 gas production year.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

3.

Methodology

3.1

Methods for determining ground motions
Two methods are commonly applied in seismic hazard assessment studies: deterministic
and probabilistic methods. These two approaches provide complementary results. The
deterministic approach usually focuses on a single earthquake scenario, ‘worst-case
scenario’, while the probabilistic approach combines all possible scenarios taking into
account their probabilities of occurrence. The two approaches can lead to different results
depending on the seismotectonic environment of the site of interest. International
guidelines and standards nowadays recommend the probabilistic approach since it allows
for a better treatment of the uncertainties. However, the deterministic approach is still
used to complement the probabilistic results and is still the reference method in a couple
of national regulations.
For this study, a probabilistic seismic hazard assessment has been conducted.

3.2

Probabilistic approach
The probabilistic method follows the traditional calculation process developed by Cornell
(1968) (Figure 3.1). The basic assumptions are that the seismicity in a region can be
modeled by a set of independent sources and a Poisson occurrence of events with time.
The estimated annual rate at which the ground motion, A, will exceed a particular value,
a, is computed by (Cornell 1968):
𝑁𝑠𝑜𝑢𝑟𝑐𝑒

𝜆[𝐴 > 𝑎] =

∑ 𝑁(𝑀𝑚𝑖𝑛 ) ∬ 𝑃[𝐴 > 𝑎|𝑚, 𝑟] 𝑓𝑀 (𝑚) 𝑓𝑅 (𝑟) 𝑑𝑚 𝑑𝑟
𝑖=1

where Nsource is the total number of sources, N(Mmin) is the annual rate of earthquakes with
magnitude greater than or equal to Mmin, 𝑃[𝐴 > 𝑎|𝑚, 𝑟] is the probability of the ground
motion, A, exceeding the threshold value a, given the earthquake magnitude m and
distance r from the source, and 𝑓𝑀 (𝑚) and 𝑓𝑅 (𝑟) are probability density functions
describing magnitude and distance. Integration over all the magnitudes and distances is
performed. The computation of this integral is carried out numerically. Under the
Poissonian assumption, the probability of exceedance in a specified exposure period t
(typically corresponding to the useful life of a project) 𝑃[𝐴 > 𝑎, 𝑡] is related to the annual
rate of exceedance 𝜆[𝐴 > 𝑎] by:
𝑃[𝐴 > 𝑎, 𝑡] = 1 − 𝑒 −[𝜆[𝐴>𝑎] 𝑡]
In this study the seismic sources are modeled as area sources in which the probability of
future earthquakes occurrence is homogeneous where no specific seismotectonic
structures can be identified.
172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

Then, the seismic activity was modeled using the well-known Gutenberg-Richter model.
The calculation process is based on a logic tree in order to propagate epistemic and
aleatory uncertainties.
The PSHA was performed using the following steps:
◼

Compilation and processing of the earthquake catalogue (including declustering and
definition of completeness periods for magnitude ranges);

◼

Identification and characterization of seismic sources of the seismotectonic models;

◼

Definition of the ground motion prediction equations (GMPEs, sometimes called
‘attenuation laws’);

◼

The construction of a logic tree based on the consideration of a set of alternative
seismotectonic models and GMPEs;

◼

Hazard calculation results in a set of hazard curves (one for each branch of the logictree). A hazard curve gives the probabilities of exceedance at the site for given
acceleration levels;

The software used to derive the seismicity parameters is the Fugro proprietary code
GEOSIS (2018) developed by Fugro and, for seismic hazard calculations, a version of RCRISIS specifically developed by Mario Ordaz et al. (2017) for Fugro. R-CRISIS is the
newest version of CRISIS, a worldwide and well-known solution for performing PSHAs
accounting for both aleatory and epistemic uncertainties using logic-tree computations.
R-CRISIS has been verified and validated during the last PEER project “Probabilistic
Seismic Hazard Analysis Code Veriﬁcation” (Hale et al. 2018). The resulting hazard model
will be the probabilistic combination of branches representing different assumptions and
the weights assigned to them.
The PSHA is performed assuming soft soil conditions in terms of VS30 defined as
representative of the Groningen Gas Field (VS30 = 200m/s).

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

Figure 3.1: Probabilistic seismic hazard method flowchart.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

4.

Seismological database
The earthquake catalogue (seismotectonic origin) used in this study correspond to that
defined for an Offshore Windfarm project (Fugro, 2020). Here below we only present only
a summary of the characteristics of the seismic catalogue developed in that project.
For a detailed explanation on the seismological database, please refer to the previous
work done for the Offshore Windfarm project, in the report P904711/SHA 01 (FUGRO,
2020).
The induced catalogue used is the synthetic seismicity catalogue provided by Dynafrax
during this project considering 4 different injection scenarios. In PSHA calculations
considering induced seismicity, the induced catalogue is “superposed” to the natural
seismic catalogue.

4.1

Earthquake catalogue with seismotectonic origin

4.1.1

Data compilation
The earthquake catalogue covering the historical and instrumental periods was compiled
for the study area from different sources and includes records spanning almost 640 years
(from 1382 to 2019). Considering the tectonic setting and the seismicity around the study
area, the compiled catalogue is made in a radius of 300 km around the site covering an
area between 0° to 15°E of longitude and 48°N to 58°N of latitude.
The following databases were considered to develop the project catalogue:
◼

The International Seismological Centre (ISC) which is a global catalogue compiling
information from various global and local networks. 19 195 events with reported
magnitudes

were

compiled

for

the

period

from

1900

to

2019.

(http://www.isc.ac.uk/iscbulletin/search/bulletin/).
◼

The Royal Netherlands Meteorological Institute (KNMI) catalogue provides
information on date, hour, coordinates, depth, local magnitude and event type of
each event. 1,408 events were compiled for the period from 1911 to 2019.
(https://data.knmi.nl/datasets/aardbevingen_catalogus/1).

◼

The British Geological Survey (BGS) catalogue which provides information on recent
and historical earthquakes. 11,595 events were compiled for the period from 1382 to
2019. (http://www.earthquakes.bgs.ac.uk/earthquakes/dataSearch.html).

◼

The Royal Observatory of Belgium (ROB) catalogue covers the period from 1382 to
2018 and includes 7 530 events. (http://seismologie.be/fr/seismologie/sismicite-enbelgique/banque-de-donnees-en-ligne).

The final compiled catalogue for this study, based on all catalogues mentioned above,
contains 33 688 earthquakes, before duplicate removal, covering the period from 1382 to
end of December 2019.
172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

4.1.2

Final compiled declustered catalogue
After the removal of duplicate events, induced seismicity, the magnitude homogenization,
we obtained a final processed catalogue with the following characteristics:
◼

23,368 events

◼

Period 1382 – 2019

◼

0.1 ≤ Mw ≤ 6.09

◼

0 ≤ hypocentral depth ≤ 87.7 km

For the declustering, the Gardner & Knopoff (1974) algorithm was used, together with the
space and time windows of Burkhard & Grünthal (2009),
The declustered catalogue is mapped on Figure 4.1.
4.1.3

Completeness analysis
In order to compute activity rates, one needs to estimate completeness periods which
correspond to dates after which all events of a given magnitude range are detected and
included in the seismicity catalogue.
The estimated completeness years for the different magnitude bins are given in Table 4.1,
including best-estimate, lower and upper bounds to consider uncertainty on these
estimations.
Table 4.1: Completeness period (best-estimates, lower and upper bounds) by magnitude bins for the
project catalogue.
Mw

Best estimate

Lower Bound

Upper Bound

[1.5 – 2.0)

1990

1980

2000

[2.0 – 2.5)

1980

1965

1995

[2.5 – 3.0)

1975

1960

1990

[3.0 – 3.5)

1970

1950

1985

[3.5 – 4.0)

1950

1920

1980

[4.0 – 4.5)

1900

1850

1950

[4.5 – 5.0)

1850

1800

1900

≥ 5.0

1800

1700

1900

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

Figure 4.1: Seismicity of the region of interest (period 1352 to 2019) – Declustered catalogue.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

4.2

Induced earthquake catalogue
The seismicity catalogues considered are presented in the Table 4.2, Table 4.3, Table 4.4
and Table 4.5 and Figure 4.2, Figure 4.3, Figure 4.4 and Figure 4.5 and they were provided
by Dynafrax (see Figure 6.1). The different injection scenarios (see Figure 6.1) are also
shown in Figure 4.2, Figure 4.3, Figure 4.4 and Figure 4.5.
Table 4.2: Seismic catalogue associated to Scenario 1 (source: Dynafrax)

ID

Longitude

Latitude

Magnitude

Occurrence date

1
2
3
4

6.6309
6.7568
6.8761
6.8259

53.2955
53.2354
53.2655
53.2838

1.93
2.22
2.44
2.39

2021.03.18
2021.07.07
2021.07.16
2024.04.28

Table 4.3: Seismic catalogue associated to Scenario 2 (source: Dynafrax)

ID

Longitude

Latitude

Magnitude

Occurrence date

1
2
3
4
5
6
7
8
9

6.8761
6.7568
6.6309
6.8664
6.8259
6.8397
6.8828
6.8528
6.8615

53.2655
53.2354
53.2955
53.1973
53.2838
53.2734
53.2934
53.2531
53.2792

2.45
2.21
1.93
1.45
2.39
1.48
2.01
1.96
2.02

2021.01.24
2021.02.04
2021.03.19
2021.06.09
2021.06.24
2021.10.02
2021.10.18
2021.12.21
2021.12.28

Table 4.4: Seismic catalogue associated to Scenario 3 (source: Dynafrax)

ID

Longitude

Latitude

Magnitude

Occurrence date

1
2
3

6.8761
6.8259
6.5647

53.2655
53.2838
53.3013

2.43
2.38
2.20

2052.02.02
2053.02.26
2063.01.24

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

Table 4.5: Seismic catalogue associated to Scenario 4 (source: Dynafrax)

ID

Longitude

Latitude

Magnitude

Occurrence date

1
2
3
4
5
6
7
8
9
10
11
12
13

6.8761
6.8259
6.8664
6.8828
6.8397
6.8254
6.7070
6.8335
6.5647
6.8728
6.9440
6.8836
6.8490

53.2655
53.2838
53.1973
53.2934
53.2734
53.2106
53.2089
53.2273
53.3013
53.2199
53.2070
53.2887
53.2578

2.41
2.38
1.43
1.99
1.41
2.18
2.25
1.47
2.20
1.99
2.24
2.15
2.06

2051.03.19
2051.07.31
2052.01.11
2052.07.11
2053.07.26
2054.02.22
2054.07.03
2055.01.18
2055.02.05
2055.08.28
2055.11.04
2055.12.26
2056.01.10

Figure 4.2: Seismic catalogue associated to Scenario 1 (source: Dynafrax).

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

Figure 4.3: Seismic catalogue associated to Scenario 2 (source: Dynafrax).

Figure 4.4: Seismic catalogue associated to Scenario 3 (source: Dynafrax).
172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

Figure 4.5: Seismic catalogue associated to Scenario 4 (source: Dynafrax).

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

5.

Seismic sources
The seismotectonic models used in this study correspond to that defined for an Offshore
Windfarm project. Below we present only a summary of the characteristics of the
seismotectonic models developed in that project.
The development of seismotectonic models and the seismic source parameters are
explained in detail in the FUGRO report P904711/SHA 01 (FUGRO, 2020).

5.1

Seismotectonic Model
Two seismotectonic models defined in a 300 km radius around the Groningen site were
used:
•

Model 1 (M1) is composed by 16 areal seismic sources (Figure 5.1). This model is
mainly derived from the French seismotectonic zonation developed for
probabilistic seismic hazard assessment for the French Metropolitan territory
published by Le Dortz et al. (2019). The Groningen site is situated in WCB zone,
which dominates the seismic hazard at the site.

•

Model 2 (M2) is composed by 17 areal seismic sources (Figure 5.2). The model is
mainly based on the seismic model published by Verbeek et al. (2009) for Belgium.
Groningen site is situated in zone 20, which dominates the seismic hazard at the
site.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

Figure 5.1: Seismotectonic model developed by Fugro (2021), modified from Le Dortz et al., 2019 for the
Groningen site – Model 1.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

Figure 5.2: Seismotectonic model 2 used of the Groningen site (After Verbeek et al., 2009, modified Fugro,
2021).

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

5.2

Induced Model
In the Groningen Gas Field, we define the model as an area of 1216 km2 (thick black line
on Figure 5.3) representative of the oval shape of the reservoir (grey area on Figure 5.3),
from the report of Dynafrax UG. The depth considered for the occurrence of earthquakes
was 3 km.

Figure 5.3: Area chosen (thick black line) to represent the Groningen Gas Field (grey area) used for the
seismic source model.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

6.

Seismic activity model

6.1

Presentation
For the natural seismotectonic sources, the computation of seismic activity rates is
based on the Gutenberg-Richter (GR) model, assuming stationary seismic rates (a
declustered catalogue was used for PSHA calculations). The TNO Model Chain was not
publicly available at the time of redaction of this report. A classical approach using
stationary models is used which is common in international practice.
The method used to determine the coefficients a and b of the model is based on the EPRI
(2012) method which is an updated and improved version of the Weichert (1980) method.
In particular, it allows:
◼

The consideration of non-uniform magnitude bins;

◼

The use of regionally varying completeness periods;

◼

The introduction of a prior on the slope of the GR model (b-value).

The regional completeness is accounted for directly in the calculation. Each zone of the
seismotectonic models can overlap with multiple zones of completeness, in which case
the observed activity rates are computed counting the number of events in each subzone and taking into account the corresponding completeness periods. It is useful for
seismotectonic models with significant regional variation of completeness periods.
However, for this study, the regional variation of completeness was not found to be
significant. Hence the single set of completeness periods is sufficient for all the zones of
the seismotectonic models considered.
The fit is performed if at least two magnitude bins are populated (implying at least two
events with different magnitudes). Otherwise, the b-value is fixed at the prior value and
only the a-value is estimated.
The prior b-value is set to 1.0 in this study.
For the induced sources, the hypothesis is to use a stationary model that will only be
valid for a specific period of time (injection time only or injection time + shut-in) to
simulate the non-stationarity of the induced model. The results will be expressed in terms
of annual exceedance rate for a specific time window.
The process used to define the GR law of the induced source is presented in Chapter 6.2.3.
In PSHA calculations using R-CRISIS (2019), the assignment of the seismic parameters for
each areal source is based on parameters λ and β which essentially corresponds to avalue and b-value respectively following the minimum magnitude considered for PSHA
calculations. Here, we used Mmin=2.0 for the induced seismicity model and Mmin=3.5 for

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

the seismotectonic model to take into account the low magnitudes associated with
induced seismicity.

6.2

Gutenberg-Richter definition

6.2.1

Seismotectonic Model 1
The classical Gutenberg-Richter (GR) fit parameters (a and b-value of the typical GR laws,
annual activity rate of the Poisson distribution considered, λ, and β value of the
exponential distribution assumed for magnitudes, and their associated errors) for all the
zones associated to seismotectonic Model 1 are shown in Table 6.1.

6.2.2

Seismotectonic Model 2
The GR fit parameters (a, b, λ, β and associated errors) for all the zones associated to
seismotectonic Model 2 (Verbeeck, 2009; modified Fugro, 2020) are shown in Table 6.2.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

Table 6.1: Gutenberg-Richter parameters for Model 1 (modified FUGRO, 2020)

Zone

Mmin
(fit)

a

a

(normalized to
106 km2)

(scaled to
area of
zone)

σa

b

σb

λ
(Mmin=3.5)

σλ
(Mmin=3.5)

β

σβ

BBR

2.0

3.6975

2.2219

0.0468

0.9642

0.0763

0.0704

0.0076

2.2201

0.1757

BOP

J94*

3.4582

1.7210

0.0300

0.9750

0.0470

0.0203

0.0048

2.2450

0.1082

EBP

3.5

3.8148

2.6586

0.1773

1.0946

0.1767

0.0672

0.0275

2.5203

0.4069

FRN

1.5

4.3774

2.6669

0.0170

0.9298

0.0288

0.2586

0.0101

2.1410

0.0662

GBP

2.0

4.2801

2.3412

0.0371

0.9248

0.0590

0.1272

0.0109

2.1294

0.1359

GEM

1.5

2.9631

1.6229

0.0450

0.8090

0.0604

0.0619

0.0064

1.8628

0.1391

HGR

2.0

2.8720

1.5272

0.1086

0.9812

0.1398

0.0124

0.0031

2.2593

0.3218

LCG

3.0

3.8728

2.5698

0.1448

1.1322

0.1650

0.0405

0.0135

2.6069

0.3798

LSH

2.0

3.8831

2.2527

0.0678

1.1341

0.1205

0.0192

0.0030

2.6113

0.2774

MCE

2.0

3.3646

2.0472

0.0523

0.9266

0.0799

0.0637

0.0077

2.1335

0.1839

MEL

2.0

4.5743

2.2144

0.0448

0.9423

0.0711

0.0825

0.0085

2.1696

0.1638

MSA

2.5

3.8210

2.5921

0.0608

1.0345

0.0977

0.0936

0.0131

2.3821

0.2249

NOP

2.0

3.4608

1.9233

0.0678

0.9752

0.1037

0.0324

0.0051

2.2455

0.2389

RHA

2.0

3.4638

1.7973

0.0670

0.9104

0.0957

0.0408

0.0063

2.0962

0.2204

SEA

2.0

3.5277

2.3049

0.0450

0.9876

0.0757

0.0705

0.0073

2.2740

0.1744

WCB (Zone of GGF)

3.5

3.5805

2.3567

0.1942

1.0336

0.1778

0.0548

0.0245

2.3799

0.4094

ZFM

2.0

4.2720

2.0348

0.0547

0.9395

0.0840

0.0558

0.0070

2.1632

0.1935

ZML

2.0

3.8205

2.3235

0.0486

1.0279

0.0850

0.0532

0.0059

2.3668

0.1956

* values from Johnston et al. (1994)

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic hazard, Hazard analysis

Fugro NL Land B.V.

Table 6.2: Gutenberg-Richter parameters for Model 2 (Verbeeck, 2009; modified Fugro, 2020)
Zone

Mmin
(fit)

a (normalized to
106 km2)

a (scaled to
area of zone)

σa

b

σb

λ
(Mmin=3.5)

σλ
(Mmin=3.5)

β

σβ

Zone 1

2.0

3.6000

2.3209

0.0393

0.9390

0.0633

0.1082

0.0098

2.1622

0.1457

Zone 2

2.5

3.5447

1.9329

0.1121

0.9858

0.1398

0.0304

0.0078

2.2700

0.3220

Zone 3

3.5

5.2466

2.6100

0.1448

1.0331

0.1645

0.0987

0.0329

2.3787

0.3788

Zone 4

1.5

4.0573

2.2339

0.0325

1.0098

0.0596

0.0501

0.0037

2.3251

0.1372

Zone 5

3.0

4.4333

2.1579

0.1535

1.0179

0.1617

0.0394

0.0139

2.3438

0.3723

Zone 6

2.5

4.7589

1.8999

0.1121

0.9734

0.1375

0.0311

0.0080

2.2412

0.3167

Zone 7

2.0

4.4081

2.2301

0.0396

0.8995

0.0605

0.1207

0.0110

2.0712

0.1394

Zone 8

2.5

3.9399

1.2908

0.2171

0.9601

0.1749

0.0085

0.0043

2.2106

0.4027

Zone 9

2.0

3.7263

1.3160

0.1254

0.9401

0.1454

0.0106

0.0031

2.1647

0.3348

Zone 10

2.0

3.2631

1.6333

0.0906

0.9567

0.1236

0.0193

0.0040

2.2028

0.2846

Zone 11

1.5

4.6191

1.8102

0.0279

0.9535

0.0479

0.0297

0.0033

2.1954

0.1102

Zone 12

1.5

4.6191

2.0958

0.0279

0.9535

0.0479

0.0574

0.0045

2.1954

0.1102

Zone 13

1.5

3.5499

1.2510

0.0887

0.9410

0.1221

0.0091

0.0019

2.1668

0.2812

Zone 14

2.5

4.3578

2.2338

0.0806

0.9914

0.1152

0.0581

0.0108

2.2827

0.2652

Zone 15

3.5

3.6900

2.4904

0.1942

1.0702

0.1790

0.0555

0.0248

2.4642

0.4121

Zone 16

2.5

3.7817

2.5356

0.0620

1.0196

0.0978

0.0927

0.0132

2.3477

0.2251

Zone 17

2.5

3.7312

2.6013

0.0486

0.9636

0.0756

0.1693

0.0189

2.2188

0.1741

Zone 18

2.0

4.8563

3.1952

0.0253

1.1764

0.0571

0.1196

0.0070

2.7087

0.1314

Zone 19

2.5

3.7842

1.7088

0.1448

0.9850

0.1560

0.0183

0.0061

2.2680

0.3592

1.8

3.6224

2.5094

0.0287

0.9861

0.0509

0.1143

0.0076

2.2706

0.1172

Zone 21

2.5

2.4941

1.1417

0.3071

1.0178

0.1880

0.0038

0.0027

2.3435

0.4329

Zone 22

2.5

3.2211

1.6385

0.1942

1.0556

0.1755

0.0088

0.0039

2.4305

0.4040

Zone 23

2.5

3.3843

2.2117

0.0806

0.9830

0.1136

0.0590

0.0110

2.2634

0.2616

Zone 20
(Zone of GGF)

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic hazard, Hazard analysis

Fugro NL Land B.V.

6.2.3

Induced Model
Dynafrax modelled the possible earthquakes that might happen depending on four
different injection scenarios (Figure 6.1).

Figure 6.1: Timeline showing the depletion, reference (when applicable indicate that no injection was
performed), injection and shut-in periods for injection scenarios 1-4 (from Figure 7-4 of the Dynafrax report
(October 2021)).

Figure 6.1 shows the different scenarios modelled in this study. This figure is coming from
the report of Dynafrax UG. They explain that the four injection scenarios, denoted cases
1-4, were simulated, each beginning after the depletion period had run for its entire 60
year duration. Injection was simulated from nine wells. For cases 3 and 4, injection was
preceded by a 30-year reference period, where pressure was allowed to naturally recover.
Depending on the case, the rates of injection will be faster (scenario 2 or 4) or slower
(scenario 1 or 3).
As explained previously, we will calculate from the simulated earthquakes a GutenbergRichter model which will be valid only for a considered period of time.
Due to the scarcity of data, a classical fit of the GR law is not possible therefore we used
an alternative approach:
•

For the b-value, we assume a b-value equal to 1 (global average value observed
in seismotectonic and induced sources around the world).

•

For the a-value calculation, we initiate our model at a Mfit of 1.0, 1.5 and/or 2.0.
And we use the mean a-value of the three Mfit.

We obtain 8 possible recurrence models as shown on Figure 6.2 depending if we choose
a period representing the injection time only (1, 5 or 10 years depending on the scenarios)

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

or a period representing the injection time plus the shut-in period (30 years for all
scenarios).

Figure 6.2: Gutenberg-Richter models for the four scenarios (1-2-3-4) with a mean a-value of all the different
Mfit (1.0, 1.5 or 2.0 depending on the scenarios). Top: rates calculated over the injection time only (1, 5 or 10
years depending on the scenarios). Bottom: rates calculated over the injection time + shut-in period (30
years for all scenarios).

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

The mean a-value for each scenario and for each period of injection is:
◼

Scenario 1: avalue=0.989 (T=10 years) and avalue=0.512 (T=30 years)

◼

Scenario 2: avalue=1.934 (T=1 year) and avalue=0.457 (T=30 years)

◼

Scenario 3: avalue=1.477 (T=10 years) and avalue=1.0 (T=30 years)

◼

Scenario 4: avalue=1.361 (T=5 years) and avalue=0.583 (T=30 years)

We compare our scenarios with the recurrence models from each host zone (WCB and
Zone 20) for the seismotectonic models SSC1 and SSC2, respectively. We normalize our
a-value to 1x106 km2 for comparison.
We observe that we have higher normalized rates for the recurrence models based on
induced seismicity using the injection period of time only (Figure 6.3). This means, the
number of earthquakes of M=2.0 per km2 is higher (between 2 and 20 times) in the
induced source than in the seismotectonic zone where Groningen is situated (Figure 6.3).
The annual activity rate of earthquakes with magnitude equal or greater than 2.0 per
km2, λ2.0 , is:
◼

SSC1 (WCB): 3.26069 x 10-5

◼

SSC2 (Zone 20): 4.46874 x 10-5

◼

Scenario 1: 8.01801 x 10-5

◼

Scenario 2: 7.06426 x 10-4

◼

Scenario 3: 2.46642 x 10-4

◼

Scenario 4: 1.88828 x 10-4

On Figure 6.4, where the period is increased (injection time + shut-in period), we notice
that scenarios 1, 2 and 4 have similar normalized rates for the recurrence models based
on induced seismicity than those based on the seismotectonic models SSC1 and SSC2.
Indeed, we increase the period from 1, 5 or 10 years up to 30 years, significantly reducing
the rates.
The value of λ2.0 per km2 is:
◼

SSC1 (WCB): 3.26069 x 10-5

◼

SSC2 (Zone 20): 4.46874 x 10-5

◼

Scenario 1: 2.67342 x 10-5

◼

Scenario 2: 2.35541 x 10-5

◼

Scenario 3: 8.22368 x 10-5

◼

Scenario 4: 3.14823 x 10-5

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

Figure 6.3: Comparison between the mean Gutenberg-Richter models of the four scenarios (1-2-3-4) with
the GR models of the seismotectonic models SSC1 and SSC2 for a period of 1, 5 or 10 years (injection time
only) depending on the scenarios.

Figure 6.4: Comparison between the mean Gutenberg-Richter models of the four scenarios (1-2-3-4) with
the GR models of the seismotectonic models SSC1 and SSC2 for a period of 30 years (injection time + shutin).

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

7.

Ground Motion Prediction Equations

7.1

Selection of GMPEs for seismotectonic sources
Ground Motion Prediction Equations (GMPEs or attenuation relationships) appropriate
for the seismotectonic context are selected with the objective to capture the epistemic
uncertainties and to fulfill the criteria of inclusiveness and exhaustiveness of the
alternatives.
The area of interest being located in a region of low seismic activity seen in the SHARE
project as a transition zone between active crustal regions of southern Europe and stable
continental regions of northern Europe, GMPEs from these two seismotectonic contexts
have been used. The set of selected GMPEs is given below with a description of the
models:
◼

Akkar et al. (2014). This model corresponds to one of the most recent GMPEs
developed for the Euro-Mediterranean regions and it represents the evolution of the
Akkar & Bommer (2010) GMPE selected in both SHARE and GEM projects.;

◼

Bindi et al. (2017): This model was specifically developed for an application of PSHA
in low to moderate seismicity areas, such as northwestern Europe.;

◼

Campbell & Bozorgnia (2014): This model is also one of the five NGA models (New
Generation of Attenuation relationship) of second generation. It can be considered
as the natural evolution of Campbell and Bozorgnia (2008).;

◼

Cauzzi et al. (2015): This is an update of Cauzzi & Faccioli (2008) model resulting
from the exploitation of a worldwide database recorded in crustal active regions and
mainly composed of Japanese data;

◼

Yenier and Atkinson (2015): This model is one of the most recent GMPE adapted
to stable continental seismic zones. It has been developed in two distinct steps: first
a generic model has been developed based on stochastic data, with parameters
calibrated on records of Californian events.

The use of the 5 selected GMPEs allows a full propagation of the epistemic uncertainties
associated with the GMPEs.
Table 7.1 summarizes the main characteristics of the selected GMPEs and their validity
domains.

7.2

Selection of GMPEs for induced seismicity sources
The attenuation of the ground motion originated by induced sources is typically different
from the attenuation observed in seismotectonic sources. Moreover, in Groningen, sitespecific GMPEs were fitted using local seismic data recorded.
Then, we tested 4 GMPEs developed in the context of induced seismicity:

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

◼

Atkinson (2015): induced GMPE from the NGA-West2 database for earthquakes
between 3.0<MW<6.0 and at distances below 40 km

◼

Bindi et al. (2017): developed for stable areas and low to moderate seismicity areas

◼

Bommer et al. (2016): updated from Dost et al. 2004 and developed specifically for
the Groningen Gas Field.

◼

Dost et al. (2004): induced GMPE developed in Netherlands (Groningen area) for
earthquakes between 1.0<ML<5.0 with the Bommer et al. 2013 adaptation

Their comparison, through a Trellis Plot is represented on Figure 7.1.

Figure 7.1: Trellis Plots at PGA for the induced GMPEs. Top: Acceleration versus epicentral distance for
MW=2.0 (left) and MW=4.0 (right). Bottom: Acceleration versus magnitude for Repi=1 km (left) and Repi=10
km (right). Black line: Atkinson 2015, blue line: Bindi et al. 2017, green line: Bommer et al. 2016 and red line:
Dost et al. 2004 (Bommer et al. 2013 adaptation)

Finally, we performed the seismic hazard calculations using the Bommer et al. 2016 GMPE,
which corresponds to the most recent and detailed GMPE specifically developed for
induced seismicity in Groningen.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

Table 7.1: Characteristics of the selected GMPEs for the site-specific PSHA to the GGF.
Magnitude

Distance

Period

Site

Region of

range

range

range

classification

application

Function of VS30

Active shallow
crustal region

Function of VS30

Stable region
Active shallow
crustal region

GMPE

Region database

Type

Akkar et al. (2014)

Euromediterranean area,
Middle East

Empirical

Mw=4.0-7.6

Bindi et al. (2017)

Worldwide

Empirical

Mw=3.0-8.0

Campbell &
Bozorgnia (2014)

Worldwide

Empirical

Cauzzi et al. (2015)

Worldwide – mostly
Japan

Yenier & Atkinson
(2015)

US (stochastic data)
(CEUS)

Rhypo, Rrup, Repi*

PGA;

= 1-200 km.

0.01-4 s

Rhypo, Rjb

PGA;

= 0-300 km.

0.01-4 s

Mw=3.3-8.5

Rrup = 0-300 km.

PGV, PGA;
0.01-10 s

Function of VS30

Empirical

Mw=4.5-7.9

Rrup =0-150 km

PGA;

EC8 site class or

0.01-10 s

Function of VS30

Active shallow
crustal region

Stochastic

Mw=3.0-8.0

Rrup =1-600 km

Function of VS30

Stable region

PGA;
0.01-10 s

*Rhypo = Hypocentral distance, Rrup = closest distance to the rupture, Repi = Epicentral distance

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic hazard, Hazard analysis

Fugro NL Land B.V.

8.

Modelling Induced Seismicity

8.1

Global Probabilistic Logic Tree
The approach using logic trees allows using alternative models. A weighting factor is
assigned for each alternative model. It is interpreted as the relative confidence of that
specific model being the correct one. The logic tree consists of a series of nodes,
representing points at which models are specified, and branches that represent the
different models specified at each node. The sum of the probabilities of all branches
connected to a given node is 1.
The weights assigned to the main branches of the logic tree associated with
seismotectonic models and GMPEs are as follows (Figure 8.1):
◼

For the seismotectonic models: an equal weight of 1/2 is assigned to each
seismotectonic model;

◼

For the GMPEs:
•

those that can be applied more specifically for the European context are
weighted to represent 50% of the logic tree. An equal weight of 25 % is assigned
to each of these GMPEs (Akkar et al. 2014, Bindi et al. 2017);

•

a 40% weight is associated for the two models (Campbell and Bozorgnia 2014,
Cauzzi et al. 2015) developed in Active Crustal Shallow regions (ASCRs) which
are calibrated on a large database of worldwide records. An equal weight of 20%
is assigned to each one.

•

a 10% weight is associated to the remaining model developed for Central Eastern
North America (Yenier and Atkinson, 2015) expressing less confidence in its
applicability in the region of interest.

This weighting scheme was adopted from an Offshore Windfarm PSHA, performed by
Fugro in 2020.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

Figure 8.1: Logic tree developed for the probabilistic seismic hazard assessment of seismotectonic model
at the Groningen site.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

8.2

Common hypotheses
The common hypothesis taken into account, for all branches of the logic tree are:
◼

A minimum magnitude Mmin = 3.5 (MW) for the calculation of the seismic hazard for
the seismotectonic models. For regions with low seismic activity (or when low return
periods are analyzed), this is a usual value and it is retained here in order to take into
account all magnitude contributions.

◼

A minimum magnitude of Mmin = 2.0 (MW) for the calculation of the seismic hazard
for the induced model. This value was selected to take into account the induced
earthquakes considered able to produce damages. The lower magnitudes are not
considered able to produce damages in the structures and were not considered in
the PSHA calculations. Moreover, magnitudes higher than 2.0 are in the domain of
validity of GMPEs used for induced sources. Finally, we have to note that this value
(2.0) is significantly lower than the 3.5 considered in seismotectonic sources because
the induced earthquakes typically occur close to the surface (hypocenters with small
depth) and, therefore, lower magnitudes of induced earthquakes can produce more
seismic ground motion at surface than same magnitudes of natural earthquakes. In
any case, the choice of Mmin = 2.0 (MW) is always subjective and based in the
previously cited considerations.

◼

Integration up to 3 standard deviations of the GMPEs (-3σ to 3σ). This is a commonly
accepted truncation value in PSHA practice for industrial projects.

◼

The integration distance used is 300 km.

◼

A shear wave velocity over the top 30 meters VS30 equal to 200 m/s.

◼

For the induced model, we specify the Mmax=4.5. The uncertainty in this value is large.
The maximum observed magnitude in Groningen is 3.6 and the maximum magnitude
coming from the different scenarios never exceeds 3.0. Therefore, the adoption of
Mmax=4.5 can be considered as a conservative approach and it correspond to other
Mmax considered in similar projects producing induced seismicity (i.e geothermal
projects).

We compute the hazard at PGA only.
The hazard is computed using the software R-CRISIS (Ordaz et al., 2018). R-CRISIS has
been tested and validated in the PEER project Hale et al. (2018), where different seismic
hazard codes were analysed.

8.3

Probabilistic Seismic Hazard Results
The results of the probabilistic seismic hazard analysis consist of seismic hazard curves
for the PGA, at the center of the Groningen Gas Field. The statistical analysis of the hazard
curves from all the logic tree branches leads to the definition of the mean hazard curves.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

8.3.1

Seismic Hazard Curve at PGA for Groningen Gas Field
Figure 8.2 (top and bottom) shows the mean seismic hazard curves of the seismotectonic
model at PGA for the Groningen Reservoir (black line). We compare this result with the
impact of adding the induced scenarios with a period of injection time only (top of Figure
8.2) or with a period of injection time + shut-in (bottom of Figure 8.2) for the gas field
(green, blue, red and magenta lines).
We observe significant differences between the seismotectonic models and the induced
models for the reduced period of injection (T=1, 5 or 10 years). For example, at 475-year,
for the smallest scenario (# 1), the acceleration is increased 1.3 times and for the strongest
scenario (# 2), the acceleration is increased 2.6 times, while remaining at a very low value
of 0.04 g at PGA if we assume a return period of 475 years for the worst case scenario
(then, assuming stationarity of the seismicity).
However, for the case of a longer period of injection (T=30 years, Figure 8.2 – bottom), at
475 years return period, for the smallest scenario (# 2), the acceleration is increased 1.1
times and for the strongest scenario (# 3), the acceleration is increased 1.3 times only.
Therefore, the introduction of an induced zone source (even with very few and weak
earthquakes) will increase the seismic hazard.
However, we emphasis 3 important considerations:
1. The comparison between the seismotectonic model and scenarios of induced
seismicity is valid only for the next 1, 5 or 10 years (depending on the period of
injection of each scenario) or for the second case, 30 years (period of injection +
shut-in for all scenarios). After those periods, the seismicity is supposed to return
to its natural period. The induced seismicity is time dependent, without injection,
there is no seismicity.
2. Even considering the induced seismicity, the acceleration levels higher than 0.1 g
have annual probabilities of exceedance lower than 0.0002. This means that the
seismic hazard is significantly increased but the probability to have higher
accelerations remains low. The increase in seismic hazard is basically associated
with the very low natural seismicity of Groningen. Then, even a little induced
seismicity produces a non-negligible increase of seismic hazard.
3. The hazard curves from the different scenarios in Figure 8.2 are only valid for the
next 1, 5, 10 or 30 years. In the context of induced seismicity, those hazard curves
are not stationary. The accelerations at 475 years return period (and others) are
only informative and valid only for the seismotectonic seismic hazard.
We should also note that compared to the previous offshore windfarm project performed
by Fugro, we obtain a higher value of acceleration at PGA at 475 years for the mean model

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

because we reduced the minimal magnitude considered in the PSHA, Mmin, at 3.5 instead
of 4.5 in order to take into account the induced seismicity. It means that in our
calculations, we are taking into account that magnitudes from 3.5 to 4.5 can produce
accelerations. In the offshore wind farm project, we considered a Mmin=4.5, the most
typical value used in industrial projects, because we consider than magnitudes lower than
4.5 cannot produce damages in new structures.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

Figure 8.2: Groningen Reservoir site – Top: Mean seismic hazard curves calculated with a period related to
the injection time only (T=1, 5, 10 years). Bottom: Mean seismic hazard curves calculated with a period
related to the injection time + shut-in period (T=30 years) – Horizontal component, VS30 = 200 m/s, 5 %
damping. Induced GMPE: Bommer et al. 2016. Note: the return periods are shown only for information
and they can only be considered for the seismotectonic model.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

Other representation of the results by scenarios directly comparing the difference in
period of injection (10 vs 30 years for scenario 1, for example).
The previous Figure 8.2 can also be represented through Figure 8.3, comparing the
periods of injection scenario by scenario.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

Figure 8.3: Mean seismic hazard curves for the Groningen Reservoir site. From Top to Bottom: Scenario 1,
Scenario 2, Scenario 3, Scenario 4. Comparison with the period of injection only (10, 1, 10 and 5 years
respectively) and the period of injection + shut-in (30 years) – Horizontal component, VS30 = 200 m/s, 5 %
damping. Induced GMPE: Bommer et al. 2016. Note: the return periods are shown only for information
and they can only be considered for the seismotectonic model.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

8.3.2

Sensitivity analysis for Mmax
In Figure 8.4 we compare for the scenario 1, the hazard curves using different Mmax from
4.5 to 5.0. We don’t not a major difference between the two different Mmax. Small
differences will be visible for lower annual probability of exceedance.

Figure 8.4: Hazard curves for the Groningen Reservoir at PGA for VS30=200 m/s. Comparison of different
Mmax used for the induced scenario 1. In green solid line: Mmax=4.5, dashed line: Mmax=5.0.

8.3.3

Sensitivity analysis for GMPEs
In Figure 8.5 we compare for the scenario 1, the hazard curves using different induced
GMPEs.
We note for induced GMPEs that Dost et al. 2004 gives, with Bindi et al. 2017, the
strongest results whereas Bommer et al. 2016 predict smaller accelerations. This is
coherent with the Trellis Plot in Figure 7.1.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

Figure 8.5: Hazard curves for the Groningen Reservoir at PGA for VS30=200 m/s. Comparison of different
induced GMPEs with the seismotectonic model. Dotted line: Atkinson 2015, dashed line: Bindi et al. 2017,
solid line: Bommer et al. 2016, dashdot line: Dost et al. 2004 (Bommer et al. 2013 adaptation).

8.3.4

Sensitivity analysis for Mmin
The Figure 8.6 presents the comparison between the use of Mmin=2.0 and Mmin=3.5 for
both the seismotectonic models and the induced seismicity models.
The black arrow represents the difference in Mmin between the seismotectonic models
only. We observe that including earthquakes between 2.0 and 3.5 has an impact on the
acceleration only for small probability of exceedance.
The blue arrow also represents the difference in Mmin between the seismotectonic models
but with the induced model. As the induced model doesn’t dominate the hazard, we can
find almost the same differences between the seismotectonic models with a Mmin=2.0 or
3.5.
The green arrow represents the difference in Mmin for the induced model between 2.0 and
3.5, while the seismotectonic models have a Mmin=3.5. We see that adding earthquakes
of magnitude 2.0<MW<3.5 has an impact until a return period of 475 years. After 475
years, we don’t see any differences.
To sum up this figure shows that the use of a weak Mmin increases the annual probability
of exceedance only for the weak accelerations (PGA< 0.06g). For PGA>0.06, the seismic
hazard curves defined using both Mmin are almost equivalent (the small differences

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

observed are produced mainly for numerical differences occurring during the integration
process of the hazard integral).

Figure 8.6: Comparison of different Mmin used for the scenario 1 using 10 years of injection time. Groningen
Reservoir site - Mean seismic hazard curves at PGA - Horizontal component, VS30 = 200 m/s, 5 % damping.
Induced GMPE: Bommer et al. 2016.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

9.

Conclusion
To explore the effect of the induced seismicity produced by the induced seismicity
predicted by the injection modelling used (using 4 scenarios), 2 site-specific probabilistic
seismic hazard analysis (PSHA) with and without consideration of induced seismicity were
performed for the Groningen Gas Field in Netherlands
The seismotectonic hazard model includes:
◼

◼

2 seismotectonic models:
•

A model specifically developed for this project (FUGRO, 2020);

•

A model adapted from the Verbeeck et al. (2009) model developed for Belgium;

5 Ground-Motion Prediction Equations (GMPEs):
•

Akkar et al. (2014);

•

Bindi et al. (2017);

•

Campbell & Bozorgnia (2014);

•

Cauzzi et al. (2015);

•

Yenier & Atkinson (2015).

The induced model considered includes:
◼

A zone source covering the Groningen Gas field (1216 km2), where we consider that
the induced seismicity is homogeneous

◼

The GMPE developed by Bommer et al. 2016 specifically developed for Groningen
site.

The comparison of seismic hazard assessment performed with and without induced
seismicity allows us to conclude that:
◼

The induced seismicity predicted by the different scenarios developed by Dynafrax
predict only few and low earthquakes with magnitudes lower than 3.0.

◼

The increase in seismic hazard produced by the induced seismic source can be
significant depending on the scenario chosen, even if the induced earthquakes are
few. The reason is because the natural seismic hazard of the Groningen area is very
low. Then, the addition of a new induced seismic source, even if this zone is not very
active and can only produce low magnitudes, produces a non-negligible increase to
the seismic hazard.

◼

The increase of seismic hazard (annual probability of exceedance of an acceleration
threshold) is basically concentrated in the years of injection. When the injection stops,
the seismicity returns to the natural (seismotectonic) situation. Considering long
periods of time and diluting the increase of seismic hazard observed during the years
of injection into a long period of time, the increase of seismic hazard is not significant.

◼

The increase in seismic hazard is not stationary and should be considered only for
the next 1, 5, 10 or 30 years considered in the Dynafrax injection simulation. After

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

these years, the modelling shows that the induced seismicity tends to cease and
should no longer be considered.
◼

The increase in seismic hazard affects mainly low accelerations (< 0.06 g).

From the seismic hazard point of view, some recommendations could be done:
◼

The best action would be to not inject any kind of fluid, independently of the rates
and duration of the injection. The models used show that if the injection is zero, the
induced seismicity tends also to zero.

◼

However, the seismicity shouldn’t be the only criteria used to take a decision. For
example, maybe the injection of fluids could allow a significant reduction of expected
settlements and the benefits of the reduction of settlement could compensate the
increment of seismic hazard.

◼

If finally, the injection is done, the preferred option would be the use of long periods
of injection with small injection rates instead of short periods of injection using higher
injection rate (the scenario 1 leads to annual exceedance rates lower than scenario 2,
Figure 8.2 top).

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

10.

References
AKKAR S. and BOMMER J. J. (2010). Empirical Equations for the Prediction of PGA, PGV,
and Spectral Accelerations in Europe, the Mediterranean Region, and the Middle
East, Seismological Research Letters, 81(2), 195-206.
AKKAR S., SANDIKKAYA M. A., and BOMMER J. J. (2014). Empirical Ground-Motion
Models for Point- and Extended-Source Crustal Earthquake Scenarios in Europe and
the Middle East, Bulletin of Earthquake Engineering, 12(1): 359 – 387
ARCADIS GERMANY GMBH & GEO-ENGINEERING.ORG GMBH (2018). Geological Desk
Study Hollandse Kust (west) Wind Farm Zone. Project ID RVO.nl: WOZ2180087, 55
p.
ATKINSON, G. M. (2008). Ground-motion prediction equations for eastern North America
from a referenced empirical approach: Implications for epistemic uncertainty,
Bulletin of the Seismological Society of America, 98:3, 1304–1318
ATKINSON, G. A. (2015). Ground-Motion Prediction Equation for Small-to-Moderate
Events at Short Hypocentral Distances, with Application to Induced-Seismicity
Hazards. Bulletin of the Seismological Society of America. 105(2).
ATKINSON, G. A. and BOORE D. M. (2006). Earthquake Ground-Motion Prediction
Equations for Eastern North America, Bulletin of the Seismological Society of
America, Volume 96, No. 6, pages 2181-2205.
ATKINSON, G. M. and BOORE, D. M. (2011). Modifications to existing ground-motion
prediction equations in light of new data, Bulletin of the Seismological Society of
America, 101:3, 1121-1135
BASILI R., et al., (2013). The European Database of Seismogenic Faults (EDSF) compiled in
the framework of the Project SHARE. http://diss.rm.ingv.it/share-edsf/
BINDI D., MASSA M., LUZI L., AMERI G., PACOR F., PUGLIA R. and P. AUGLIERA (2014),
Pan-European ground motion prediction equations for the average horizontal
component of PGA, PGV and 5 %-damped PSA at spectral periods of up to 3.0 s
using the RESORCE dataset, Bulletin of Earthquake Engineering, 12(1), 391 – 340
BINDI D., COTTON F., KOTHA S.R., BOSSE C., STROMEYER D. and GRUENTHAL G. (2017)
Application-driven ground motion prediction equation for seismic hazard
assessments in non-cratonic moderate-seismicity areas", J. Seismology, 21(5), 1201
– 1218
BOMMER J. J., DOUGLAS J., SCHERBAUM F., COTTON F., BUNGUM H. and FAH D. (2010).
On the selection of ground-motion prediction equations for seismic hazard analysis.
Seismological Research Letters, Vol. 81, n° 5, p. 783-793BOSWORTH, W. (1989).
Basin and range style tectonics in East Africa.
BOMMER J. J., DOST B., EDWARDS B., STAFFORD P. J., VAN ELK J., DOORNHOF D.,
NTINALEXIS M. (2016). Developing an Application‐Specific Ground‐Motion Model
for Induced Seismicity. Bulletin of the Seismological Society of America; 106 (1):
158–173

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

BOMMER J. J., B. EDWARDS, P.P. KRUIVER, A. RODRIGUEZ-MAREK, P.J. STAFFORD, M.
NTINALEXIS, E. RUIGROK, J. SPETZLER (2019). V6 Ground-Motion Model (GMM) for
Induced Seismicity in the Groningen Field with Assurance Letter. NAM report.
BROUWER J.W.R. (2010). The meaning of Eurocode 8 and Induced Seismicity For
Earthquake Engineering in the Netherlands. Fifth International Conference on
Recent Advances in Geotechnical Earthquake Engineering and Soil Dynamics, May
24-29, San Diego, California.
BURKHARD M. & GRÜNTHAL G. (2009). “Seismic source zone characterization for the
seismic hazard assessment project PEGASOS by the Expert Group 2 (EG1b)”. Swiss
Journal of Geosciences, Volume 102, Issue 1, p. 149–188.
CAMPBELL K.W. (2003) Prediction of Strong Ground Motion Using the Hybrid Empirical
Method and Its Use in the Development of Ground Motion (Attenuation) Relations
in Eastern North America" (Bulletin of the Seismological Society of America, Volume
93, Number 3, pages 1012-1033).
CAMPBELL K. W. and BOZORGNIA Y. (2008). NGA Ground Motion Model for the
Geometric Mean Horizontal Component of PGA, PGV, PGD and 5 % Damped Linear
Elastic Response Spectra for Periods Ranging from 0.01 to 10s" (Earthquake Spectra,
Volume 24, Number 1, pages 139 - 171)
CAMPBELL K. W. and BOZORGNIA Y. (2014). NGA-West2 Ground Motion Model for the
Average Horizontal Components of PGA, PGV, and 5 % Damped Linear Acceleration
Response Spectra" (Earthquake Spectra, Volume 30, Number 3, pages 1087 - 1115).
CAUZZI C. and FACCIOLI E. (2008). Broadband (0.05 to 20 s) prediction of displacement
response spectra based on worldwide digital records. Journal of Seismology, Vol.
12, n° 4, p. 453-475.
CAUZZI C., FACCIOLI E., VANINI M. et BIANCHINI A. (2015). Updated predictive equations
for broadband (0.01 to 10 s) horizontal response spectra and peak ground motions,
based on a global dataset of digital acceleration records. Bulletin of Earthquake
Engineering, vol. 13, Issue 6, p. 1587-1612.
CHIOU B. S.-J. and Youngs R. R. (2008). An NGA Model for the Average Horizontal
Component of Peak Ground Motion and Response Spectra" (Earthquake Spectra,
Volume 24, No. 1, pages 173-215)
CORNELL, C. A. (1968). Engineering seismic risk analysis. Bulletin of the seismological
society of America, 58(5), 1583-1606.
COTTON F., SCHERBAUM F., BOMMER J.J., BUNGUM H. (2006). Criteria for selecting and
adjusting ground motion models for specific target regions: application to central
Europe and rock sites. J. Seismol., 10, p. 137–156.
CROOK, Th. de, (1996), A seismic zoning map conforming to Eurocode 8, and practical
earthquake parameter relations for the Netherlands, Geologie en Mijnbouw, 75, pp
11-18.
DELAVAUD, E., F. SCHERBAUM, N. KÜHN and T. ALLEN (2012). Testing the Global
Applicability of Ground- Motion Prediction Equations for Active Shallow Crustal
Regions. Bulletin of the Seismological Society of America. Vol. 102, n°2, p. 707-721
172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

DERRAS B., BARD P. Y., COTTON F. (2014). Toward fully data driven ground-motion
prediction models for Europe", Bulletin of Earthquake Engineering 12, 495-516
DOST, B., van Eck, T. and HAAK, H. (2004) Scaling of peak ground acceleration and peak
ground velocity recorded in the Netherlands. Bollettino di Geofisica Teorica ed
Applicata. 45(3), 153 – 168
DOUGLAS, J. (2010). Consistency of ground-motion predictions from the past four
decades, Bulletin of Earthquake Engineering, 8, 1515-1526.
EC8 (2003). Eurocode 8: Design of structures for earthquake resistance - Part 1: General
rules, seismic actions and rules for buildings. Ref. No. prEN 1998-1:2003 E.
FUGRO REPORT P904711/01 “Geotechnical Report Investigation Data Seafloor In Situ
Test Locations Hollandse Kust (west) Wind Farm Zone Dutch Sector, North Sea”
GARDNER J.K. & KNOPOFF L. (1974). Is the sequence of earthquakes in Southern
California, with aftershocks removed, Poissonian? Bulletin of the Seismological
Society of America, 64, p. 1363–1367.
GIARDINI, D., GRUNTHAL, G., SHEDLOCK, K. M & ZHANG, P. (1999). The GSHAP global
seismic hazard map. Ann. Geofisica 42: 1225-1230;
GIARDINI, D., J. WOESSNER, L. DANCIU (2014). Mapping Europe’s Seismic Hazard. EOS,
95(29): p. 261-262.
GRÜNTHAL, G., BOSSE, C., SELLAMI, S., MAYER-ROSA, D. & GIARDINI, D. (1999).
Compilation of the GSHAP regional seismic hazard for Europe, Africa and the Middle
East. www.seismo.ethz.ch/static/GSHAP/index.html.
GRÜNTHAL, G., WAHLSTRÖM, R., STROMEYER, D. (2013): The SHARE European
Earthquake Catalogue (SHEEC) for the time period 1900–2006 and its comparison
to the European-Mediterranean Earthquake Catalogue (EMEC). - Journal of
Seismology, 17, 4, 1339-1344.
GRÜNTHAL G., STROMEYER D., BOSSE C. (2017). The Data set of the earthquake model
for the probabilistic seismic hazard assessment of Germany, version 2016. Report
on supplementary material for the respective publication. Scientific technical report
STR – Data; 17/05, PostDam
GRÜNTHAL, G., STROMEYER, D., BOSSE, C., COTTON, F. and BINDI, D. (2018). “The
probabilistic seismic hazard assessment of Germany — version 2016, considering
the range of epistemic uncertainties and aleatory variability “. Bulletin of Earthquake
Engineering
HALE C., ABRAHAMSON N. and BOZORGNIA Y. (2018). Probabilistic Seismic Hazard
Analysis Code Verification. PEER Report 2018/03
JOHNSTON, A.C., COPPERSMITH, K.J., KANTER, L.R. and CORNELL, C.A. (1994). “The
Earthquakes of Stable Continental Regions”. EPRI Technical Report TR-102261, 5
Volumes, Electric Power Research Institute, Palo Alto, California.
KNMI seismic catalogue. http://rdsa.knmi.nl/dataportal/.

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

LE DORTZ K., COMBES P., CARBON D. (2019). An alternative seismotectonic zonation for
probabilistic and deterministic seismic hazard assessment for Metropolitan France.
Submitted to the Bull. Soc. Geol. France.
ORDAZ, M., MARTINELLI, F., AGUILAR, A., ARBOLEDA, J., MELETTI, C., & D'AMICO, V.
(2017). R-CRISIS. Program and platform for computing seismic hazard.
ORDAZ, M., MARTINELLI, F., AGUILAR, A., ARBOLEDA, J., MELETTI, C., & D'AMICO, V.
(2018). R-CRISIS. Program and platform for computing seismic hazard.
ORDAZ, M. and SALGADO-GALVEZ M. A. (2019). R-CRISIS. Validation and Verification
Document. ERN Technical Report. Mexico City, Mexico.
PAGANI M.; MONELLI D.; WEATHERILL G.; DANCIU L.; CROWLEY H.; SILVA V.; HENSHAW
P.; BUTLER L.; NASTASI M.; PANZERI L.; SIMIONATO M.; VIGANO D. (2014).
OpenQuake Engine: An Open Hazard (and Risk) Software for the Global Earthquake
Model. Seismological Research Letters 85 (3): p. 692-702.
PAGANI M., J. GARCIA-PELAEZ, R. GEE, K. JOHNSON, V. POGGI, R. STYRON, G.
WEATHERILL, M. SIMIONATO, D. VIGANÒ, L. DANCIU, D. MONELLI (2018). Global
Earthquake Model (GEM) Seismic Hazard Map (version 2018.1 - December 2018).
STUCCHI M., et al. (2013). The share European earthquake catalogue (SHEEC) 1000–1899.
Journal of seismology, volume 17, issue 2, p. 523–544.
TAVAKOLI B. and PEZESHK S. (2005) Empirical-Stochastic Ground-Motion Prediction for
Eastern North America" (Bull. Seism. Soc. Am., Volume 95, No. 6, pages 283-2296)
TNO (2020a). Probabilistic Seismic Hazard and Risk Analysis in the TNO Model Chain
Groningen, TNO 2020 R11052.
TNO (2020b). TNO Model Chain Groningen: Update and quick scan comparison of 2020
HRA model, TNO 2020 R11659.
TORO G. R., ABRAHAMSON N. A., SCHNEIDER J. F. (1997) Model of Strong Ground
Motions from Earthquakes in Central and Eastern North America: Best Estimates
and Uncertainties. (Seismological Research Letters, Volume 68, Number 1)
TORO, G.R. (2002). Modification of the Toro et al. (1997) attenuation equations for large
magnitudes and short distances, Technical Report, Risk Engineering.
VERBEECK K., VANNESTE K., CAMELBEECK T. (2009). “Seismotectonic zones for
probabilistic seismic-hazard assessment in Belgium”. ONDRAF/NIRAS report,
Belgium agency for radioactive waste and enriched fissile materials. NIROND-TR
Report 2008-31E. 61 pp.
WEICHERT, D.H. (1980). Estimation of the earthquake recurrence parameters for unequal
observational periods for different magnitudes, Bulletin of the Seismological
Society of America, 70, p. 1337-134
WOESSNER, J., DANCIU, L., GIARDINI, D. and the SHARE consortium. (2015). “The 2013
European Seismic Hazard Model: key components and results”. Bulletin of
Earthquake Engineering.
YENIER E. and ATKINSON G. M. (2015). Regionally Adjustable Generic Ground-Motion
Prediction Equation to Central and Eastern North America" BSSA, vol 105
172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis

Fugro NL Land B.V.

ZHAO, J.X., ZHANG, J., ASANO, A., OHNO, Y., OOUCHI, T., TAKAHASHI, T., OGAWA, H.,
IRIKURA, K., THIO, H.K., SOMERVILLE, P.G., FUKUSHIMA, Y., FUKUSHIMA, Y. (2006)
Attenuation relations of strong ground motion in Japan using site classification
based on predominant period, Bulletin of the Seismological Society of America,
96:3, 898–913

172147-REP01-FNV_SHA_KEM24 03 | KEM-24: Effect of pressure maintenance by fluid injection on seismic
hazard, Hazard analysis