2015 Scienti􀅮ic Paper

Phase 1 update May 2015:
trend changes in ground
subsidence in Groningen

The views expressed in this paper are those of the author and do not necessarily re􀅭lect the
policies of Statistics Netherlands.

F.P. Pijpers, D.J. van der Laan

Contents
1

Introduction

3

2

The 􀅮iltered time series for each station

3

Drift trends

4

Alternative model

5

Re-processing of GPS data by TU Delft

4

7
8
9

5.1 ﬁtting to the TU Delft data without additional ﬁltering

6

10

Conclusions 13

Statistics Netherlands | Scienti􀅭ic paper 2015

2

Nederlands
Deze rapportage is een verslag van onderzoek dat is uitgevoerd in het kader van fase 1 van
een onderzoeksproject door het CBS in opdracht van Staatstoezicht op de Mijnen (SodM).
Dit onderzoek is ten behoeve van een statistische onderbouwing van het meet- en
regelprotocol voor gasexploitatie in de provincie Groningen. Het onderwerp van dit rapport
is een heranalyse van trends in de bodemdaling in de provincie Groningen, gerapporteerd in
fase 0, in December 2014. Voor de voorliggende analyse is de tijdreeks voor de GPS
gegevens aangevuld tot aan 3 Januari 2015. Daarnaast zijn alle GPS data opnieuw
gecalibreerd, in een onafhankelijk onderzoek uitgevoerd door de TU Delft.
Zoals ook bleek uit de eerdere analyses is er een statistisch signiﬁcante trendbreuk in de
bodemdaling ongeveer 2 maanden nadat de productie sterk was gereduceerd. De
trendwijziging kan zich geleidelijk gemanifesteerd hebben over een periode van enkele
weken, en er is daarom een onzekerheid van ruwweg een week of twee over de centrale
datum van deze overgang.
English
This is a report on the research that has been carried within phase 1 of a project being
carried out by Statistics Netherlands and commissioned by State Supervision of Mines
(SodM). This research is part of the underpinning of the statistical methods employed to
support the protocol for measurement and regulation of the production of natural gas in
the province of Groningen. The subject of this report is to re-analyse the trends in the
ground subsidence in and around that region, ﬁrst reported in phase 0, in december 2014.
For this re-analysis the time series of the GPS data has been updated with more recent
measurements up to Jan. 3 2015. The entire GPS data set has also be re-calibrated,
independently from the previous calibrations, by researchers at the TU Delft.
In accordance with the earlier report, it is found that there has been a signiﬁcant change in
the ground subsidence. The changeover in the trend of ground subsidence can have
become manifest gradually over a period of several weeks, which implies that the central
date of the transition is also uncertain by roughly a week or two.

1 Introduction
In this update for the most part the same methods are applied to the new time series of GPS
data. For this reason the reader is referred to the earlier report for an overview of the methods
and the methodological background. The time series being analysed in this report contains the
entire period that was analysed in the previous report, but is extended forward to cover the
epoch up to January 3 2015. Also, the most up-to-date information is being used regarding
reference stations. In practice this means that the entire set of GPS data is re-calibrated for all
of the stations, and for the entire period.

Statistics Netherlands | Scienti􀅭ic paper 2015

3

Figure 1 The relative locations with names of each of the GPS stations from which
data are available. The red diamond symbols are locations within the area of the
production 􀅮ield `Middle', where production has been reduced. The stations indicated
with blue squares are locations outside of the main 􀅮ield, collectively designated as
`Edge', although at the southwestern station Eemskanaal there are wells in
production. The two stations indicated with green triangles are also production
locations in a separate part of the 􀅮ield `South' where production has not been
reduced. Well locations where production has been reduced are shown as black
crosses, with in addition three wells very close to this region (AMR, TJM, SDB). The
indicated scales are in km
5945

Usquert

5940

't Zand
ZND

5935
Zeerijp
5930

LRM

Stedum

Delfzijl
POS
ten Post

5925

PAU
5920

Overschild

AMR

OVS

Middle

TJM
SDB

Edge

Tjuchem

South

Eemskanaal

Wells
5915
Zuiderveen

Froombosch
5910

5905

Veendam

5900
435

440

445

450

455

460

465

2 The 􀅮iltered time series for each station
The ﬁrst processing step of the time series of the GPS stations is to apply the ﬁltering procedure
with the ﬁlter factors as derived in the previous report. As is done in the previous report, the
average height of each station is subtracted. The results are shown in ﬁgure 2
The quantity of interest is the ‘sagging’ of the stations in group ‘Middle’ and in group ‘South’
with respect to those in group ‘Edge’. It is this diﬀerential displacement over a spatial scale of
the order of, or smaller than, the ﬁeld that will lead to the build-up of stresses. Therefore the
displacements, averaged for each group, are subtracted in the sense ‘Middle’ minus ‘Edge’ and
‘South’ minus ‘Edge’. The result is shown in ﬁgure 3.
To assess the extent to which correlations might still be present in the diﬀerential displacement,
two additional diﬀerential time series are shown in ﬁgure 4.
One measure concerns the diﬀerential displacement between stations that are all within group
‘Middle’, in the sense (Zeerijp+ten Post) − (Stedum+Overschild) which is shown as an orange
Statistics Netherlands | Scienti􀅭ic paper 2015

4

0.005

sted

−0.005

Date

tenp
Date

0.005

tjuc

−0.005

Date

zand
Date

−0.005

0.005

zeer
Date

M

A

M

J

J

A

S

O

N

D

J

F

M

A

M

J

J

A

S

O

N

D

2014

J

2015

Date
Date

0.008

usqu

veen

−0.004

0.002

Date

dzyl

0.008

Date

eems
Date

0.002

Relative displacement (m)
Relative
displacement
(m)
Relative displacement
Relative displacement (m)
Relative
displacement
(m) (m)

over

−0.004

Relative displacementRelative
(m)
displacementRelative
(m) Relative
displacement
Relative
(m)
displacement
(m)
displacementRelative
(m)
displacement (m)
displacement
(m) Relative

Figure 2 The time series of the GPS height, after 􀅮iltering out intermediate time scale
variations. For presentational purposes for each station the average of its time series
over the period from day 50 to day 214 is subtracted. The vertical scale is in m. Top
Panel: GPS stations near the centre of the production 􀅮ield, bottom panel: GPS
locations in group edge.

M

A

M

J

J

A

S

O

N

D

J

F

2014

M

A

M

J

J

A

S

O

N

D

J

2015

Date
Date

line. The other concerns the diﬀerential displacement between stations that are all within group
‘Edge’, in the sense (Usquert+Veendam) − (Delfzijl+Eemskanaal) shown as a purple line. In the
absence of remaining biases these diﬀerential measurements should not show any trends.
Statistics Netherlands | Scienti􀅭ic paper 2015

5

0.000

0.002

0.004

Middle − Edge
South − Edge

−0.002
−0.004

Relative displacement (m)

0.006

Figure 3 The time series of the 􀅮iltered GPS height for the differences between the
group averages `Middle' − `Edge' (in red), and `South' − `Edge' (in yellow).

M

A

M

J

J

A

S

O

N

D

J

F

M

A

M

J

J

A

S

O

N

2014

D

J

2015

Date

While it is clear that for most of 2014 no trend is present in these sets of diﬀerential

0.000

0.002

0.004

(udqu+veen)−(dzyl+eems)
(zeer+tenp)−(sted+over)

−0.002
−0.004

Relative displacement (m)

0.006

Figure 4 The differential time series of the 􀅮iltered GPS height within groups. Purple
line: differential measurements between stations all within group `Edge'. Orange line:
differential measurements between stations all within group `Middle'.

M

A

M

J

J

A

S

O

N

D

J

F

M

2014

A

M

J

J

A

S

O

N

D

J

2015

Date

displacements, in October there is a sudden change in particular in the group ’Edge’, but also to
some extent in the group ’Middle’. At this stage it is unclear whether this is caused by further
calibration issues that require resolving, or whether the altered spatial pattern of production
could have inﬂuenced these measures.

Statistics Netherlands | Scienti􀅭ic paper 2015

6

3 Drift trends
Since the reduction of the gas production at some locations near stations in the group ‘Middle’
has been rapid, it is of interest to explore whether a break can be found in the linear component
of the downward drift for the time series for ‘Middle’ minus ‘Edge’. To this end one can ﬁt not
only a single straight line, using standard least-squares ﬁtting, but also introduce a break-point
with a diﬀerent straight line ﬁt before and after that point. Given that the latter ﬁt has more
degrees of freedom than the former, this should be accounted for in assessing whether the
two-section ﬁt with a break is genuinely better than a single straight line. The Akaike
Information Criterion (AIC) is an equivalent measure to the reduced 𝜒 􀍯 measure used in the
previous reports.

−184000
−186000

AIC

−182000

−180000

Figure 5 AIC score as a function of the position of a break in the linear model 􀅮it of
the difference between the group averages of `Middle' and `Edge'. The minimum for
`jump' is at 2014-02-06 and that of `kink' at 2014-04-02.

Jump
Kink
M

A

M

J

J

A

S

O

N

D

J

F

M

A

M

J

J

A

S

O

N

2014
Date

Statistics Netherlands | Scienti􀅭ic paper 2015

7

0.002

0.004

0.006

Middle−Edge
Trend
Trend w. kink
Trend w. jump

0.000

Relative displacement (m)

0.008

Figure 6 Filtered times series for the group averages of `Middle' and `Edge' with the
􀅮itted trends.

M

A

M

J

J

A

S

O

N

D

J

F

M

A

M

J

J

A

S

O

2014

N

D

J

2015

Date

4 Alternative model
One possible issue with the previous ﬁt of the break in the trend, is that the number of
measurement contributing to each of the average lines changes with time. Before Februari
2014 there is only one measurement station contributing to each of the lines. The previous
model does not take this into account: one would expect the model to give more weight to
diﬀerences at later dates as more measurement stations are contributing to the averages.
One possibility to avoid this problem is to not calculate the diﬀerence, but to model all time
series combined. The diﬀerence in slopes between the two groups can be modelled explicitly
and it is possible to test if this diﬀerence is signiﬁcant.
A piecewise linear continuous ﬁt can be written as
(1)

𝑦􀐕 = 𝛽􀍭 + 𝛽􀍮 𝑡􀐕 + 𝛽􀍯 𝑠(𝑡 − 𝑡􀍭 )𝑡􀐕 + 𝑒􀐕 ,

with 𝑠(𝑡) the step function (𝑠(𝑡) = 1 for 𝑡 ≥ 0 and 0 otherwise). We can then ﬁt each of the time
series with a piecewise linear continuous ﬁt:
𝛽􀍭 + 𝛽􀍮 𝑡􀐕 + 𝛽􀍯 𝑠(𝑡 − 𝑡􀍭 )𝑡􀐕 + 𝑒􀐕
𝑦􀐕 = {
(𝛽􀍭 + 𝛽􀍰 ) + (𝛽􀍮 + 𝛽􀍱 )𝑡􀐕 + (𝛽􀍯 + 𝛽􀍲 )𝑠(𝑡 − 𝑡􀍭 )𝑡􀐕 + 𝑒􀐕

if ‘Edge’,
if ‘Middle’.

(2)

We can then test if the coeﬃcient 𝛽􀍲 is zero. In ﬁgure 7 the coeﬃcient signiﬁcantly diﬀers from
zero. Also in ‘Edge’ the discontinuity is signiﬁcant (𝑏𝑒𝑡𝑎􀍱 ≠ 0). Therefore, both ‘Edge’ and
‘Middle‘ decrease in time and for both this decrease decreases from approximately 2014-03-14.
However, for ‘Middle’ this decrease is much stronger than for ‘Edge’, suggesting that this
diﬀerence is related to the gas production.

Statistics Netherlands | Scienti􀅭ic paper 2015

8

Figure 7 Filtered time series of each of the measurement stations in `Edge' and
`Middle' with the 􀅮itted model. The break is located at 2014-03-14.

0.005
−0.005

0.000

Height (m)

0.010

Edge
Middle

M

A

M

J

J

A

S

O

N

D

J

F

2014

M

A

M

J

J

A

S

O

N

D

J

2015

Date

5 Re-processing of GPS data by TU Delft
The GPS data used has been processed from the original satellite measurements using one
standard algorithm and associated method. In order to assess some of the assumptions
regarding sources of error and their eﬀect, a diﬀerent algorithm has been applied, developed at
the TU Delft. Their dataset of reduced GPS signal per GPS ground location has been made
available, so that the same diﬀerence technique can be applied, as has been done for the
previous reduced GPS datasets. The data are shown in ﬁgure 8, which is the equivalent of ﬁgure
2.
Fig. 9 shows the ‘sagging’ of the stations in group ‘Middle’ and in group ‘South’ with respect to
those in group ‘Edge’ , ie. the equivalent of what is shown in ﬁgure 3.
The diﬀerential displacement between stations that are all within group ‘Middle’, in the sense
(Zeerijp+ten Post) − (Stedum+Overschild) is shown as an orange line in ﬁgure 10. The other
measure concerns the diﬀerential displacement between stations that are all within group
‘Edge’, in the sense (Usquert+Veendam) − (Delfzijl+Eemskanaal) shown as a purple line. In the
absence of remaining biases these diﬀerential measurements should not show any trends.
Comparing ﬁgure 10 with ﬁgure 4 shows that for the diﬀerential time series between the
stations that are all in the group ’Edge’ the behaviour after October 2014 is very diﬀerent. This
would appear to indicate that in that epoch there are indeed some calibration issues in the
earlier dataset which the reprocessing by the TU Delft have resolved. The same epoch does still
show some anomalies for the time series involving stations within the group ’Middle’.

Statistics Netherlands | Scienti􀅭ic paper 2015

9

over

10

sted

0

5

Date

tenp
Date

10

tjuc

0

5

Date

zand
Date

10

zeer

0

5

Date

M

A

M

J

J

A

S

O

N

D

J

F

M

A

M

J

J

A

S

O

N

2014

D

J

F

2015

Date
Date

0 2 4 6 8

usqu

veen

−4

Date

dzyl

0 2 4 6 8

Date

eems
Date

−4

Relative displacement (mm)
Relative
displacement
(mm)
Relative displacement
Relative displacement (mm)
Relative
displacement
(mm)(mm)

Relative displacementRelative
(mm) displacementRelative
(mm)
displacement
Relative
(mm) displacement
Relative
(mm) displacementRelative
(mm) displacement (mm)
Relative
displacement
(mm)

Figure 8 The time series of the GPS height as determined by the TU Delft method,
and after 􀅮iltering out intermediate time scale variations. Note that the vertical scale
is in mm. Top Panel: GPS stations near the centre of the production 􀅮ield, bottom
panel: GPS locations in group edge.

M

A

M

J

J

A

S

O

N

D

J

F

M

A

2014

M

J

J

A

S

O

N

D

J

F

2015

Date
Date

5.1 􀅮itting to the TU Delft data without additional 􀅮iltering
While the ﬁtting procedures used are robust to the precise properties of the errors, it is
worthwhile to explore to what extent the point-to-point correlation introduced by the ﬁltering
Statistics Netherlands | Scienti􀅭ic paper 2015

10

6

Figure 9 The time series of the TU Delft processed GPS height data, and after
􀅮iltering, for the differences between the group averages `Middle' − `Edge' (in red),
and `South' − `Edge' (in yellow). Note that the vertical scale is in mm.

2
0
−4

−2

Relative displacement (mm)

4

Middle − Edge
South − Edge

M

A

M

J

J

A

S

O

N

D

J

F

M

A

M

J

J

A

S

O

N

2014

D

J

F

2015
Date

6

Figure 10 The differential time series of the TU Delft processed GPS height data, and
after 􀅮iltering, within groups. Purple line: differential measurements between
stations all within group `Edge'. Orange line: differential measurements between
stations all within group `Middle'.

2
0
−4

−2

Relative displacement (mm)

4

(udqu+veen)−(dzyl+eems)
(zeer+tenp)−(sted+over)

M

A

M

J

J

A

S

O

N

D

J

F

M

2014

A

M

J

J

A

S

O

N

D

J

F

2015
Date

process aﬀects the conclusions drawn from the trend ﬁtting. At this stage it is worth
remembering that the choice of ﬁts does not have an underlying physical model for the
subsidence. The intention is merely to assess whether one can reject the hypothesis that there
is no change in average subsidence rate over the entire period of observation.
Figure 11 is the equivalent of ﬁgure 5 but now the ﬁtting is done on GPS data as processed by

Statistics Netherlands | Scienti􀅭ic paper 2015

11

the TU Delft without any additional ﬁltering by CBS. While there are diﬀerences in the AIC for

39000
37000

AIC

41000

43000

Figure 11 AIC score as a function of the position of a break in the linear model 􀅮it of
the difference between the group averages of `Middle' and `Edge'. The minimum for
`jump' (red curve) is at 2014-04-05. For the `kink' (yellow curve) there are several
choices which are equally valid, which can be seen from the three broad minima. The
􀅮inal minimum on 2014-03-01 coincides within the uncertainties with what is
obtained for the 􀅮iltered data (􀅮ig. 5)

Jump
Kink
M

A

M

J

J

A

S

O

N

D

J

F

M

A

M

J

J

A

S

O

N

2014
Date

the ﬁts as determined for the ﬁltered and unﬁltered data, as one would expect, in both cases
one concludes that ﬁts for which the average subsidence rate changes over the period of
observation are superior and statistically signiﬁcantly so. The data appear to point to a trend
change in late March early April 2014, with an uncertainty of a few weeks either way, possibly
because the change is not quite as abrupt as the ﬁts are.

Statistics Netherlands | Scienti􀅭ic paper 2015

12

Figure 12 The 􀅮itted models, with either a jump or a kink, to the data as processed by
the TU Delft without any additional 􀅮iltering for the best-􀅮itting position of that break
in the linear model 􀅮it of the difference between the group averages of `Middle' and
`Edge'.

4
2
0
−2

Relative displacement (mm)

6

Middle−Edge
Trend
Trend w. kink
Trend w. jump

M

A

M

J

J

A

S

O

N

D

J

F

M

2014

A

M

J

J

A

S

O

N

D

J

F

2015
Date

6 Conclusions
From the GPS data it can be concluded that there is continued subsidence of the ground in the
area of the wells where production was reduced in the month of January of 2014. However the
rate of subsidence is lower after a date. The date of this break has been redetermined with the
more extended and re-calibrated set of GPS data. In addition a completely independent data
reduction of the GPS data, developed at the TU Delft, leads to the same conclusions in our
analysis. The location of the break is around mid-March or very early in April, ie. about 9 weeks
after the reduction in production, but there is a fair margin of uncertainty (some 2 weeks)
around the exact value of the time gap.
While there is clear statistical evidence that a break has occurred in the subsidence rates, and
the reduction in subsidence speeds is measurable with a high degree of precision, research is
now in progress to perform a more systematic analysis of the correlation between the
production time series and the GPS height time series.

Statistics Netherlands | Scienti􀅭ic paper 2015

13

References
Nederlandse Aardolie Maatschappij BV, 2013, A technical addendum to the winningsplan
Groningen 2013 Subsidence, Induced Earthquakes and Seismic Hazard Analysis in the
Groningen Field.
Bracewell, R.N., 1965, The Fourier transform and its applications, 5ﬀ, 24ﬀ, McGraw-Hill
(paperback edition: 1999).
Pijpers, F.P., 2014, Phase 0 report 1 : signiﬁcance of trend changes in ground subsidence in
Groningen, Statistics Netherlands

Statistics Netherlands | Scienti􀅭ic paper 2015

14