Zgryźliwość kojarzy mi się z radością, która źle skończyła.
Chapter 17
Ground Surface Temperature Histories
Reconstructed from Boreholes in Poland:
Implications for Spatial Variability
Darius Mottaghy, Jacek Majorowicz, and Volker Rath
17.1 Introduction
In recent years, geothermal logs have been used to reconstruct ground surface
temperature histories (GSTH) (Nielsen and Beck
1989
; Dahl-Jensen et al.
1998
;
Pollack et al.
1998
; Huang et al.
2000
; Šafanda and Rajver 2001; Kukkonen and
Joeleht
2003)
. In spite of the fact that the temporal resolution of the reconstruction
is low and the coupling of surface air temperature to ground temperature is still not
well understood, such reconstructions represent the most direct record of past tem-
peratures and are potentially valuable in analyzing past climate conditions. Many
authors performed reconstructions on various timescales, from a few 100 to 100,000
years. The short-time reconstructions serve to investigate natural and anthropogenic
climate change, whereas the longer reconstructions aim to understand temperature
variations back to the last Glacial and its transition to our current climate.
In this study we present results from ground surface temperature histories recon-
structed from two deep boreholes in Poland. Our main interest lies in the magnitude of
the Pleistocene-Holocene Warming (PHW). It is well known that there are large differ-
ences in the PHW depending on the location, this is studied for example by Demezhko
et al.
(2007)
for Northern Eurasia. However, there is still a lack of data for a densely
sampled spatial distribution of reliable reconstructions of the PHW magnitude.
One data set origins from the Udryn borehole, located in northeastern Poland. It is
a low heat flow region due to the presence of a norite-anorthosite massif, with values
around 40 mW m
−2
. The second data set originates from Czeszewo which is situated
further west (Fig.
17.1
). Here, heat flow is considerably higher (>80 mW m
−2
).
D. Mottaghy (
)
Geophysica GmbH, Lütticherstr, 32, D-52064, Aachen, Germany
e-mail: d.mottaghy@geophysica.de
J. Majorowicz
Northern Geothermal, Edmonton, Canada
V. Rath
PalMA Research Group, Dpt. of Astrophysics and Atmospheric Sciences, Faculty of Physical
Sciences, Universidad Complutense de Madrid, Madrid, Spain
R. Przybylak et al. (eds.),
The Polish Climate in the European Context:
An Historical Overview
, DOI 10.1007/978-90-481-3167-9_17,
© Springer Science + Business Media B.V. 2010
375
376
D. Mottaghy et al.
Fig. 17.1
Location of the boreholes Czeszewo and Udryn on a heat flow map (Majorowicz et al.
2007)
of eastern Europe (mW m
−2
), uncorrected for paleoclimatic effects
In both cases there is a continuous temperature, porosity, and thermal conductivity
log available. These are favorable conditions for a reliable reconstruction of GSTH
by inversion.
17.2 Inverse Method
We invert borehole temperatures for GSTH using the Tikhonov method. In the fol-
lowing we shortly describe the equations and relations involved in the forward
model, the inverse technique, as well as the regularization parameters, which are
necessary for the type of inversion applied here.
17 Ground Surface Temperature Histories Reconstructed from Boreholes in Poland
377
17.2.1 Forward Model
The solution of any inverse problem requires that the corresponding forward model
is known and adapted to the physical processes involved. As the use of analytical
solutions does not allow for appropriately complex situations, we use a straightfor-
ward numerical solution of the heat equation by finite difference methods. This
gives us considerable freedom of introducing more realistic physics of the systems,
as the non-linear dependencies of the system parameters on temperature and pres-
sure, and, in particular, the effects of latent heat.
The one-dimensional, purely conductive heat equation in a porous medium can
be written as:
∂
∂
+=
l
T
hc
() ,
∂
T
(1)
∂ ∂
e
z
e
t
where
l
is thermal conductivity (W m
−1
K
−1
),
r
is density (kg m
3
),
c
is specific heat
capacity (J kg
−1
K
−1
), and
h
is volumetric heat production (W m
−3
). The subscript
e
marks effective parameters of the porous medium, and can be interpreted as proper-
ties of a two-phase mixture between solid rock and fluid-filled pore space. This
mixture is characterized by porosity
F
. Usually, the geometric mean is chosen for
thermal conductivity, that is,
l
e
= l
w
F
l
m
(1
−
F)
, with the indices
w
and
m
denoting the
fluid and matrix contribution. Volumetric quantities like (
rc
) are averaged arith-
metically, taking
(rc)
e
= F r
w
c
w
+ (1
−
F) r
m
c
m
.
For the paleoclimate application we have in mind, Equation (1) usually is solved
with the appropriate boundary conditions, namely a fixed but time-dependent tem-
perature at the top (
z
=
z
0
),
T Tt z z
0
= =
()
at
,
(2)
and fixed heat flow at the bottom,
∂
∂
T
=− =
q
b
at
zz
.
(3)
z
λ
b
Equation (1) is understood to allow all coefficients, boundaries, and sources, to be
nonlinearly dependent on temperature. Most rocks matrix properties exhibit moder-
ate dependencies on temperature, while pressure can safely be neglected for the
depths under consideration here (Chapman
1986)
. As we are aiming at deep bore-
holes constraining the history of ground surface temperature for some 10,000 years,
we have to extend the numerical model to depths of several 1,000 m, and tempera-
tures of up to 200°C accordingly. This requires taking the temperature dependen-
cies of the thermophysical properties into account. In many practical situations, the
upper layers of the subsurface are of sedimentary origin and show considerable
porosities up top 40%. In this case the physical properties of the fluid have a dis-
tinct influence, as they change with temperature, but unlike the matrix, they are also
pressure-sensitive. We implemented the relations given in Table
17.1
, defining the
temperature (and pressure in the case of fluid density) dependence of the physical
r
∂
z
378
D. Mottaghy et al.
Table 17.1
Dependencies luid and ice properties on temperature and pressure
Property
Reference
λ
f
T
Phillips et al.
(1981)
c
f
T
derived from enthalpy given by Zylkovskij et al.
(1994)
, salinity neglected
r
f
(
T,P
) from Zylkovskij et al.
(1994)
, salinity neglected
λ
i
(
T
),
c
i
(
T
),
r
i
(
T
)
Lide
(2000)
parameters within the forward module. For the domain characteristic of the
boreholes investigated in this study (<5,000 m depth), a comparison of the fluid
properties with the ones from Wagner and Pruß
(2002)
assures that deviations are
in an admissible range (<2%). Certainly, these errors are much smaller than the
ones induced by neglecting of salinity which can play a major role in some
permafrost environments (Ippisch
2001)
.
Equation (1) is solved by a standard one-dimensional finite difference scheme,
allowing for irregular grids, variable coefficients, time-dependent boundary condi-
tions, and source terms. The resulting non-linearity is handled by a simple fixed-
point iteration, which proved to be adequate to the problem. Time-stepping is done
by a general two-level scheme, including the well-known Forward/Backward Euler
and Crank-Nicholson schemes (see, e.g.,Wood
1990)
as special cases.
The simplified form of the heat equation used in the forward problem of palaeo-
climate inversion is easily modified to incorporate the effects of the presence of the
ice phase, which turned out to be of first-order importance for the Udryn study (see
Section 17.3). Details of theory, implementation, and the validation of the approach
can be found in Mottaghy and Rath
(2006)
. The approach chosen is of course a
simplification of a process which is far more complicated than described here
(Ippisch
2001)
. For the application to palaeoclimate inversion we have restricted
ourselves to the simplified approach, as mainly processes with large time constants
are important, and thus the processes in the upper few meters of the subsurface do
not play a significant role for the reconstruction of palaeotemperatures at scales
larger then a few years.
17.2.2 Inverse Technique
Given observed data, that is, recent borehole temperature measurements as a function
of depth,
d
i
= T(z
i
, t
0
)
, the GSTH,
T(0, t)
, can be estimated by a regularized least-
squares procedure. To achieve this, an objective function is set up to be
minimized:
Wdgp Wpp Wpp
(
( ))
2
τ
0
(
)
2
τ
1
(
)
2
(4)
d
2
0
p
a
1
p
a
2
2
Here, d − g(p) ≡ r is the residual vector between the data
d
and the solution of the
forward problem g(p) at these points for a given parameter vector p. The weighted
Θ= − + − + −
17 Ground Surface Temperature Histories Reconstructed from Boreholes in Poland
379
norm of this residual represents the data fit. Data weighting is introduced by W
d
,
which is usually used to standardize the residuals, that is, its diagonal is set to the
inverse square root of the data covariance. The second term in Equation (4) is
defined by the application of a linear operator W
d
on the deviations of the model
parameters p from their preferred values p
a
.
To solve the inverse problem we try to find the minimum of functional (4) by
formally differentiating equating the result to zero. This leads to the well-known
normal equation,
(
WJ WJ W W W W p
)
T
− −− −
τ τ δ
( )
0
T
0
( )
1
T
1
=
(5)
d
d
0
p
p
1
p
p
(
WJr WWpp WWpp
)
T
τ
( )
0
T
0
(
)
τ
( )
1
T
1
(
),
d
0
p
p
a
1
p
p
a
which can be iterated as a modified Gauss-Newton iteration with p
k
+1
= p
k
+
d
p
k
,
where
d
p is determined for each iteration by Equation (5). The derivative matrix
J
(
J
acobian) results from the differentiation of the residuals d – g(p). It is defined as:
J
∂
∂
=,
g
i
ij
p
j
where
p
j
are the model parameters. Differentiation is done by an adaptive method
using divided differences.
To solve the linear system (5), an equivalent formulation can be found, which is
very flexible and allows using a variant of the well-known conjugate gradient
method for its solution:
WJ W d gp
T
(
−
( ))
(6)
d
d
1
0
1
0
τ δ τ
2
W p Wpp
=
− −
2
(
)
0
p
0
p
a
τ τ
1
W Wpp
1
− −
1
1
(
)
2
2
1
p
1
p
a
This rectangular system of linear equations is efficiently solved in each itera-
tion by conjugate-gradient type methods like CGLS or LSQR (Björck
1996
;
Hansen
1997)
. Methods using only the gradient of an objective function like the
one defined in Equation (4), that is, the Jacobian only enters through the matrix-
vector product –J
T
r, may be a good alternative choice for the minimization tech-
nique. The most prominent of these are the variants of Nonlinear Conjugate
Gradients (NLCG), or Quasi-Newton (QN) algorithm (see, e.g., Nocedal and
Wright
1999)
. However, the above-mentioned methods make it difficult to
optimize regularization parameters.
In the case of GSTH inversion, we parameterized the surface temperatures as a
series of step functions for
p
. Number and temporal spacing of steps are set a priori,
leaving the temperature values for each period as inversion parameters. These
parameters are associated to the time steps indirectly. The time discretization of the
forward problem thus can be chosen following numerical requirements, indepen-
dently from the inverse grid employed. Due to the diffusive character of the
+ +
zanotowane.pl doc.pisz.pl pdf.pisz.pl hannaeva.xlx.pl
Ground Surface Temperature Histories
Reconstructed from Boreholes in Poland:
Implications for Spatial Variability
Darius Mottaghy, Jacek Majorowicz, and Volker Rath
17.1 Introduction
In recent years, geothermal logs have been used to reconstruct ground surface
temperature histories (GSTH) (Nielsen and Beck
1989
; Dahl-Jensen et al.
1998
;
Pollack et al.
1998
; Huang et al.
2000
; Šafanda and Rajver 2001; Kukkonen and
Joeleht
2003)
. In spite of the fact that the temporal resolution of the reconstruction
is low and the coupling of surface air temperature to ground temperature is still not
well understood, such reconstructions represent the most direct record of past tem-
peratures and are potentially valuable in analyzing past climate conditions. Many
authors performed reconstructions on various timescales, from a few 100 to 100,000
years. The short-time reconstructions serve to investigate natural and anthropogenic
climate change, whereas the longer reconstructions aim to understand temperature
variations back to the last Glacial and its transition to our current climate.
In this study we present results from ground surface temperature histories recon-
structed from two deep boreholes in Poland. Our main interest lies in the magnitude of
the Pleistocene-Holocene Warming (PHW). It is well known that there are large differ-
ences in the PHW depending on the location, this is studied for example by Demezhko
et al.
(2007)
for Northern Eurasia. However, there is still a lack of data for a densely
sampled spatial distribution of reliable reconstructions of the PHW magnitude.
One data set origins from the Udryn borehole, located in northeastern Poland. It is
a low heat flow region due to the presence of a norite-anorthosite massif, with values
around 40 mW m
−2
. The second data set originates from Czeszewo which is situated
further west (Fig.
17.1
). Here, heat flow is considerably higher (>80 mW m
−2
).
D. Mottaghy (
)
Geophysica GmbH, Lütticherstr, 32, D-52064, Aachen, Germany
e-mail: d.mottaghy@geophysica.de
J. Majorowicz
Northern Geothermal, Edmonton, Canada
V. Rath
PalMA Research Group, Dpt. of Astrophysics and Atmospheric Sciences, Faculty of Physical
Sciences, Universidad Complutense de Madrid, Madrid, Spain
R. Przybylak et al. (eds.),
The Polish Climate in the European Context:
An Historical Overview
, DOI 10.1007/978-90-481-3167-9_17,
© Springer Science + Business Media B.V. 2010
375
376
D. Mottaghy et al.
Fig. 17.1
Location of the boreholes Czeszewo and Udryn on a heat flow map (Majorowicz et al.
2007)
of eastern Europe (mW m
−2
), uncorrected for paleoclimatic effects
In both cases there is a continuous temperature, porosity, and thermal conductivity
log available. These are favorable conditions for a reliable reconstruction of GSTH
by inversion.
17.2 Inverse Method
We invert borehole temperatures for GSTH using the Tikhonov method. In the fol-
lowing we shortly describe the equations and relations involved in the forward
model, the inverse technique, as well as the regularization parameters, which are
necessary for the type of inversion applied here.
17 Ground Surface Temperature Histories Reconstructed from Boreholes in Poland
377
17.2.1 Forward Model
The solution of any inverse problem requires that the corresponding forward model
is known and adapted to the physical processes involved. As the use of analytical
solutions does not allow for appropriately complex situations, we use a straightfor-
ward numerical solution of the heat equation by finite difference methods. This
gives us considerable freedom of introducing more realistic physics of the systems,
as the non-linear dependencies of the system parameters on temperature and pres-
sure, and, in particular, the effects of latent heat.
The one-dimensional, purely conductive heat equation in a porous medium can
be written as:
∂
∂
+=
l
T
hc
() ,
∂
T
(1)
∂ ∂
e
z
e
t
where
l
is thermal conductivity (W m
−1
K
−1
),
r
is density (kg m
3
),
c
is specific heat
capacity (J kg
−1
K
−1
), and
h
is volumetric heat production (W m
−3
). The subscript
e
marks effective parameters of the porous medium, and can be interpreted as proper-
ties of a two-phase mixture between solid rock and fluid-filled pore space. This
mixture is characterized by porosity
F
. Usually, the geometric mean is chosen for
thermal conductivity, that is,
l
e
= l
w
F
l
m
(1
−
F)
, with the indices
w
and
m
denoting the
fluid and matrix contribution. Volumetric quantities like (
rc
) are averaged arith-
metically, taking
(rc)
e
= F r
w
c
w
+ (1
−
F) r
m
c
m
.
For the paleoclimate application we have in mind, Equation (1) usually is solved
with the appropriate boundary conditions, namely a fixed but time-dependent tem-
perature at the top (
z
=
z
0
),
T Tt z z
0
= =
()
at
,
(2)
and fixed heat flow at the bottom,
∂
∂
T
=− =
q
b
at
zz
.
(3)
z
λ
b
Equation (1) is understood to allow all coefficients, boundaries, and sources, to be
nonlinearly dependent on temperature. Most rocks matrix properties exhibit moder-
ate dependencies on temperature, while pressure can safely be neglected for the
depths under consideration here (Chapman
1986)
. As we are aiming at deep bore-
holes constraining the history of ground surface temperature for some 10,000 years,
we have to extend the numerical model to depths of several 1,000 m, and tempera-
tures of up to 200°C accordingly. This requires taking the temperature dependen-
cies of the thermophysical properties into account. In many practical situations, the
upper layers of the subsurface are of sedimentary origin and show considerable
porosities up top 40%. In this case the physical properties of the fluid have a dis-
tinct influence, as they change with temperature, but unlike the matrix, they are also
pressure-sensitive. We implemented the relations given in Table
17.1
, defining the
temperature (and pressure in the case of fluid density) dependence of the physical
r
∂
z
378
D. Mottaghy et al.
Table 17.1
Dependencies luid and ice properties on temperature and pressure
Property
Reference
λ
f
T
Phillips et al.
(1981)
c
f
T
derived from enthalpy given by Zylkovskij et al.
(1994)
, salinity neglected
r
f
(
T,P
) from Zylkovskij et al.
(1994)
, salinity neglected
λ
i
(
T
),
c
i
(
T
),
r
i
(
T
)
Lide
(2000)
parameters within the forward module. For the domain characteristic of the
boreholes investigated in this study (<5,000 m depth), a comparison of the fluid
properties with the ones from Wagner and Pruß
(2002)
assures that deviations are
in an admissible range (<2%). Certainly, these errors are much smaller than the
ones induced by neglecting of salinity which can play a major role in some
permafrost environments (Ippisch
2001)
.
Equation (1) is solved by a standard one-dimensional finite difference scheme,
allowing for irregular grids, variable coefficients, time-dependent boundary condi-
tions, and source terms. The resulting non-linearity is handled by a simple fixed-
point iteration, which proved to be adequate to the problem. Time-stepping is done
by a general two-level scheme, including the well-known Forward/Backward Euler
and Crank-Nicholson schemes (see, e.g.,Wood
1990)
as special cases.
The simplified form of the heat equation used in the forward problem of palaeo-
climate inversion is easily modified to incorporate the effects of the presence of the
ice phase, which turned out to be of first-order importance for the Udryn study (see
Section 17.3). Details of theory, implementation, and the validation of the approach
can be found in Mottaghy and Rath
(2006)
. The approach chosen is of course a
simplification of a process which is far more complicated than described here
(Ippisch
2001)
. For the application to palaeoclimate inversion we have restricted
ourselves to the simplified approach, as mainly processes with large time constants
are important, and thus the processes in the upper few meters of the subsurface do
not play a significant role for the reconstruction of palaeotemperatures at scales
larger then a few years.
17.2.2 Inverse Technique
Given observed data, that is, recent borehole temperature measurements as a function
of depth,
d
i
= T(z
i
, t
0
)
, the GSTH,
T(0, t)
, can be estimated by a regularized least-
squares procedure. To achieve this, an objective function is set up to be
minimized:
Wdgp Wpp Wpp
(
( ))
2
τ
0
(
)
2
τ
1
(
)
2
(4)
d
2
0
p
a
1
p
a
2
2
Here, d − g(p) ≡ r is the residual vector between the data
d
and the solution of the
forward problem g(p) at these points for a given parameter vector p. The weighted
Θ= − + − + −
17 Ground Surface Temperature Histories Reconstructed from Boreholes in Poland
379
norm of this residual represents the data fit. Data weighting is introduced by W
d
,
which is usually used to standardize the residuals, that is, its diagonal is set to the
inverse square root of the data covariance. The second term in Equation (4) is
defined by the application of a linear operator W
d
on the deviations of the model
parameters p from their preferred values p
a
.
To solve the inverse problem we try to find the minimum of functional (4) by
formally differentiating equating the result to zero. This leads to the well-known
normal equation,
(
WJ WJ W W W W p
)
T
− −− −
τ τ δ
( )
0
T
0
( )
1
T
1
=
(5)
d
d
0
p
p
1
p
p
(
WJr WWpp WWpp
)
T
τ
( )
0
T
0
(
)
τ
( )
1
T
1
(
),
d
0
p
p
a
1
p
p
a
which can be iterated as a modified Gauss-Newton iteration with p
k
+1
= p
k
+
d
p
k
,
where
d
p is determined for each iteration by Equation (5). The derivative matrix
J
(
J
acobian) results from the differentiation of the residuals d – g(p). It is defined as:
J
∂
∂
=,
g
i
ij
p
j
where
p
j
are the model parameters. Differentiation is done by an adaptive method
using divided differences.
To solve the linear system (5), an equivalent formulation can be found, which is
very flexible and allows using a variant of the well-known conjugate gradient
method for its solution:
WJ W d gp
T
(
−
( ))
(6)
d
d
1
0
1
0
τ δ τ
2
W p Wpp
=
− −
2
(
)
0
p
0
p
a
τ τ
1
W Wpp
1
− −
1
1
(
)
2
2
1
p
1
p
a
This rectangular system of linear equations is efficiently solved in each itera-
tion by conjugate-gradient type methods like CGLS or LSQR (Björck
1996
;
Hansen
1997)
. Methods using only the gradient of an objective function like the
one defined in Equation (4), that is, the Jacobian only enters through the matrix-
vector product –J
T
r, may be a good alternative choice for the minimization tech-
nique. The most prominent of these are the variants of Nonlinear Conjugate
Gradients (NLCG), or Quasi-Newton (QN) algorithm (see, e.g., Nocedal and
Wright
1999)
. However, the above-mentioned methods make it difficult to
optimize regularization parameters.
In the case of GSTH inversion, we parameterized the surface temperatures as a
series of step functions for
p
. Number and temporal spacing of steps are set a priori,
leaving the temperature values for each period as inversion parameters. These
parameters are associated to the time steps indirectly. The time discretization of the
forward problem thus can be chosen following numerical requirements, indepen-
dently from the inverse grid employed. Due to the diffusive character of the
+ +