Heat flow density values and paleoclimate determined from stochastic inversion of four temperature-depth profiles from the Superior Province of the Canadian Shield

Tectonophysics,

345

164 (1989) 345-359

Elsevier Science Publishers

B.V., Amsterdam

- Printed

in The Netherlands

Heat flow density values and paleoclimate determined from stochastic inversion of four temperature-depth profiles from the Superior Province of the Canadian Shield S.B. NIELSEN ’ and A.E. BECK 2 ’ Laboratory for Geophysics, Aarhus University, Finlandsgade ’ Department

of Geophysics,

6-8, DK-8200 Aarhus (Denmark)

University of Western Ontario, London, Ont. N6A 587 {Canada)

(Received

and accepted

January

3,1988)

Abstract Nielsen,

S.B. and Beck, A.E., 1989. Heat flow density

four temperature-depth E.R. Decker

profiles

(Editors),

surface

obtained temperature

history

(1) to provide available history,

temperature

history

how uncertainty

manner

structure

The bounds

about

structure.

is applied

ground

While computationally

and (3) to provide

numerical

during

The pattern of heat flow density at the surface of the earth is one of the principal data for testing evolutionary models of the lithosphere and for probing its present-day thermal structure. For this pattern to be more than just a qualitative constraint, reliable heat flow density observations are necessary. This paper deals with the problem of obtaining reliable heat flow density values from the measurement of borehole thermal properties when the geothermal field is influenced by surface 0 1989 Elsevier Science Publishers

The method

from

B.V.

due to been

models

of the surface

the forward limits,

four

600-m

the Wisconsin

method

fails

(2) to extract

all

of interest.

A

treats the surface

temperature

which have been constrained

unrealistic

to constrain

heat flow from shallow

Introduction

0040-1951/89/$03.50

data

may be biased

limits on the parameters

and prevent

of and

have ~nventionally

and

convenient,

heat flow as unknowns

show how it is possible

temperature

which

with confidence

confidence

is described.

stability

to geothermal

background

from boreholes conduction

data

inversion

L. Rybach

164: 345-359.

corrections

of heat

of all available

and the background

Shield. The examples the average

on the equation

these shortcomings

ensure

and the steady-state

observations

from stochastic

Shield. In: V. cerm&

Tectonophysics,

heat flow density

based

determined

of the Canadian

Structure.

representation

which overcomes

is safe. The technique of the Canadian

Province

This calls for pale~~matic

properties

quantitative

conductivity

with soft bounds.

Province

problem

in an optimum

inverse method

hypothesis

conditions.

and the thermal

a stringent

the thermal

that terrestrial

a forward

information

stochastic priori

temperature

by solving

values and paleoclimate

Heat Flow and Lithosphere

It has long been recognized non-steady

from the Superior

deep

solutions boreholes

simultaneously

geothermal Glaciation

data,

when

a

the a-priori

in the Superior the recent

surface

and furthermore

show

may be taken into account.

temperature variations. In this context, “reliable” (Wang and Beck, 1987) means a heat flow density value accompanied by a statistically valid error bound. In principle, (q)

is obtained

the terrestrial by taking

heat flow density the product

of the

thermal conductivity (A) and the temperature gradient @7’/&) so that 4 = h - H’,&, with the depth z increasing downwards. This procedure yields a local conductive heat flow density value with error bounds derivable from the observational uncertainties in h and iP/az. Under ideal

346

circumstances

this

heat

tegrates lithospheric production decay

thermal

thermal

tectonothermal lithosphere,

processes

value

in-

such as heat

perturbations

mobilization

the

due

of the crust

and heat conducted

through

to

and the the litho-

presence

of, flow,

geneity

for

example,

erosion

in thermal

and

temperature

variations flow density

however,

et

quantitative,

assessment

ground-

answer-an

hetero-

values

and

or surface influence

obscure

their

based

phenomenon: al.,

1988),

on a model

groundwater

flow

Beck, 1977). If there is more than one disturbance-causing phenomenon, it is usual to apply corrections seriatum. An inevitable consequence of this approach is an increase in the variance of the result. This is because slightly different yet equally plausible models of the disturbance produce different corrections. This scatter of corrections, generally around a preferred value, may conveniently be by a variance

of the correction.

of unwanted

effects may

therefore

of the uncertainty knowledge

phenomenon? inverse

problem depth

and

approach

variations

to make a

more

resulting

about

This

of climatic

reliable, from

the

the disturbance-

paper

presents

designed

one

to treat the

perturbations

combined

with

in the thermal

parameters.

The

basic idea follows the major

it possible

Vasseur

difference

et al. (1985)

that a variable

but with

surface

tem-

perature history is included. This leads to differences in the forward procedure outlined in the following section.

of

e.g., topog-

(Mansure and Reiter, 1979) and the influence of past climatic changes (cermak and Jessop, 1971;

described

more

topography,

may significantly

heat flow density

the disturbance-causing (Powell

We pose the question-is

sedimentation, structure,

value

(also when the effects are zero), the variances

causing

relationship with deep-seated lithospheric thermal processes. The conventional approach to these problems is the calculation of corrections to the

raphy

the observed

not have been accumulated.

that the

properties

local heat

observed

stripping

lack of precise

from below.

It has long been recognized, water

density

due to the decay of radioelements, of

sphere

flow

When

The forward problem

Given the thermal properties structure, the background heat flow and the surface temperature history, we obtain the model temperatures at the depth points of observation by a forward calculation. The solution of the forward problem requires the choice of a model and a parameterization. The model should be as simple as possible, yet consistent with the need to represent the information available. In most cases of paleoclimatic corrections, this requirement

is met by a laterally

homo-

earth model with a depth-dependent

ther-

adding the preferred correction to q, the variance of the correction adds to the variance already in q.

geneous

The fact that the preferred correction may be zero does not imply that the variance of the preferred correction is zero. A statement such as “the influence of climatic changes on this heat flow density observation is believed to be negligible” merely says that the preferred correction is zero because the average ground temperature is believed to

background heat flow, and a laterally homogeneous surface temperature history. This model can explain most of the observed conductive temperature-depth profiles. A variety of analytical forward solutions exists in the literature. They may be quite useful, but

have been fairly constant for a sufficiently long time. The possibility that the ground temperature may in fact have deviated from the assumed mean for significant periods of time, however, produces the variance of the zero correction. This leads to the conclusion that the error bounds that frequently accompany published heat flow density values may not be reliable; in the process of

mal

properties

structure,

laterally

homogeneous

always invoke some manner of restriction on the physical model. The case of a piecewise linear thermal conductivity depth dependence was solved by Wang et al. (1986) but requires a periodical surface temperature history. Truly transient surface temperature histories have been treated analytically by Nielsen and Balling (1985), but require the earth to be composed of elementary layers each of constant thermal properties.

347

In this paper we avoid the restrictions imposed by analytical methods and adopt a numerical approach. The depth-dependent thermal properties of the half space is represented by a one-dimensional finite element grid, and the time evolution is performed in discrete time steps. In the following the method is outlined. The one-dimensional time-dependent heat equation with the time-independent dynamic heat capacity pc and heat generated at a rate H per unit volume reads (Carslaw and Jaeger, 1959): 8 (XilT,‘&z),‘~z

- pdT/i%

= -H

(1)

The surface boundary condition is T( z = 0, t) = f(z), where f(t) denotes the surface temperature history. At infinite depth, the requirement is that the heat flow be equal to the time-independent background heat flow: h. aT/& = qbg. The finite element approximation is obtained in the usual way (e.g., Becker et al., 1981) by choosing a set of interpolation functions Ok(z), k = 1,...t N, on the finite depth interval zc(O, Z,,,). Z max is sufficiently large that the surface temperature history has no effect at that depth. By stating the unknown temperature solution to eqn. (1) as Too(z, t) = C,T,(t) +Q,(z) and demanding that the weighted residual integrals (Galerkin’s criterion) izrnax{ a(xar*/a+faz

dDk(z)

- pcazy3z

This algorithm is 0( At’) accurate and unconditionally stable. As the initial condition for the algorithm of eqn. (4) we use the steady-state temperature-depth profile as determined by surface temperature, thermal conducti~ty, heat production and background heat flow. Figure 1 shows a comparison of the numerical results with an analytical solution. The surface temperature history is a 5-K step change at 1000 yrs B.P., and the geological model a half space composed of two layers of constant thermal capao ity but different thermal conductivity and diffusivity. It is clear that linear interpolation on the finite element grid is sufficient to obtain the desired inaccuracy of 1 mK. Since each line element has only two nodes, the system matrix is t&diagonal. This matrix need only be assembled once when

500 YBP

+ 121

= 2 Wm-‘I(-’

Xl

= 8 ’ 10m7m2 5-l

AZ

= 3 Wm-‘K-l

x2

= 12, 10e7m25-l

200

dz, k=l,...,N

(2)

all be zero, a system of ordinary differential equations in the unknown nodal temperatures, T,(t), k=l,..., N, arises. With u = (T,, T,, . . . , TN)‘, A = a matrix of thermal conductivities, B = a matrix of dynamic heat capacities, and g = a vector of heat productions, the inte~ation of eqn. (2) results in: A

A,

.u+B.a2qat=g

TEMPERATURE

(B/At+A/‘2)-u,+,

0 100 200 300 400 600 800 1000

o'c

5.0000 3.2027 1.6820 0.9453 0.4764 0.0859 0.0096 0.0007

Fig. 1. Comparison for a two-layer Carslaw

between half

and Jaeger

for the numerical mesh

extends

temperature for

14

IN

Depth Cm) analytical numerical -----._._______-_____._____._._._____.______________

(3)

A finite difference method is used to propagate u of eqn. (3) in time. We choose the Crank-Nicholson time integration algorithm (e.g., Carey and Oden, 1984) which relates the solution u~+~ at time At(k + 1) to the solution uk at time Atk by

=gk+;+(B,‘At-A/2)-z+

my

1~~12

ma~mum

space.

5.0000 3.2030 1.6819 0.9446 0.4754 0.0851 0.0094 0.0006

analytical The

and numerical

analytical

(1959, pp. 320-322). solution

to depth

z,,

solution

solutions is from

The time step length

is At = 10 yrs. The finite element = 5 km. The distances

between

nodes are AZ = 25 m for 0 5 z I 1 km. AZ = 50 m km, absolute

and

A~=500

m for

2~2~5

error is 1 mK and occurs around

km.

The

z = 400 m.

348

the

time

further time

step

length

is kept

constant.

use of Gauss-Choleski stepping

substitution (4) making

reduces

By the

factorization,

to forward

Some synthetic examples

the In order to illustrate

and backward

sweeps on the right-side vector in eqn. computation very efficient.

to elucidate

inverse method, profiles

Tbe inverse problem A solution and

model

parameters

thermal heat

to the inverse (surface

conductivity

flow) which

problem

reasonable

is a set of

values

temperature

structure

is consistent

and

of the history,

unstable

and non-unique.

The non-lin-

at three

times

background

heat flow. Prior to the sudden

ture was X = 3 Wm-’ K-‘, pc = 2.5. lo6 J me3 K-‘, A = 1 PW rnm3, and q = 40 mW me2. The

pc is a fixed function

ple at 25 m and the deepest

depth

profile

the thermal

of depth.

This non-linearity

changes

smoothly

conductivity-depth

with changes profile.

equidistant choice

in

The lack

of stability means that some model parameters do not depend continuously on the data. Small varia-

homoge-

temperature change, the surface temperature is 2°C. For all cases, the thermal properties struc-

temperature-depth

not very severe, as the temperature-

a

over the entire

and a laterally

earity shows up in the rock thermal parameters, i.e., in our case in the thermal conductivity, since is, however,

change

in the past,

production

data interval,

were gener-

half space with

temperature

different

rate of heat

temperature neous

background

with the observed

data within their statistical accuracy. The inverse climatic perturbation problem suffers from being non-linear,

located constant

isotropic surface

and of the

temperature-depth

Three datasets

ated for a homogeneous

physically

procedure

and limitations

some synthetic

were inverted.

a 5-K step function reliable

the analysis

the possibilities

being

depth

profiles points

dictated

sampled

sample

by typical

the real cases discussed

at 25 sam-

at 600 m, this

field practice

in

later.

Although the stochastic not formally distinguish “parameters”

were

with the shallowest

inverse between

we will use “data”

method “data”

does and

for the observed

tions (e.g., errors) in the observed subsurface temperatures may produce very large changes in the

temperature-depth profile and “parameters” for the generally less well known conductivity-depth

surface temperature history. The lack of uniqueness means that certain characteristics of the data may be explained in more than one way. A curvature in the temperature-depth profile, for example, may arise from surface temperature varia-

profile, ground

tions,

of the measurement

thermal

diogenic

conductivity-depth

heat production,

variations,

ra-

or some combination

of

these. We use a unified approach to process ill-posed and weakly non-linear inverse problems proposed by Tarantola (1987). The

and Valette (1982a, b) and Tarantola technique employs the concept of

using the a-priori knowledge about model parameters to stabilize the numerical procedure and prevent unphysical solutions. In short, the method combines the state of information on the values of the dataset, the state of information on the unknown parameters, and the state of information on the theory relating data and unknowns. Other applications to heat flow studies are given by Vasseur et al. (1985), Nielsen (1986) and Wang and Beck (1987).

surface temperature heat flow.

history

and

back-

Two sources of noise in the data are considered: observational errors and theoretical errors. Observational errors originate in the internal noise sources

that prevent

equipment the precise,

and

from

absolute

other

observa-

tion of a well-defined temperature at a specific depth. For modem temperature logging equipment, this Theoretical This is the effects that

error is approximately k < 10 mK. errors arise from an imperfect theory. case when the data contain physical are not described by the theory. The

thermal conductivity-depth profile, for example, is not a smooth piecewise linear or quadratic function of depth. Taking averages over some meters or tens of meters, the thermal conductivity may be close to constant or may show a smooth trend, but on a more local scale it can be very erratic. These small-scale conductivity fluctuations are not part of the theory and the effect should be considered as noise in the temperature data. The

349

normally distributed noise of u = 18 mK to the generated temperature data. In the following, we describe the assi~ent of a-priori values of parameters and data. The a-priori data information is the actually observed temperatures with associated uncertainties. The temperature data are shown as circles in the temperature-depth plots of Figs. 2A-C. The choice of a-priori model parameters is guided by observations, previous experience and by what is believed to be reasonable. The thermal conductivity was usually measured, and the values may serve as a guide for the choice of an a-priori ~ondu~ti~ty model and its confidence limits. For industrial wells in sedimentary environments, the

amplitude of theoretical errors may vary greatly from one datum to another and from one case to another. For example, in the real cases to be considered later the thermal conductivity-depth profiles exhibit small-scale fluctuations but generally no large-scale trends, and the theoretical errors were estimated to be + 15 mK. At the other extreme, to explain a temperature-depth profile influenced by groundwater flow in terms of a purely conductive mode1 may require theoretical errors in the 1-K range if unphysical solutions are to be avoided. The two noise processes may be considered uncorrelated with normally dist~buted amplitudes and were simulated by adding zero mean pseudo-

A Q = 2225

‘C

mWm_2

,,~~..____.‘~“‘.____--~,

‘\

0

‘,

_----__

__-----

G’

..\_,

_______\

*a---.__

J’

n 6

-5

0

1000

500

1000

500

1500

2000

2500

1500

2000

2500

15 Q

.g

=

C

39 t 5 mWmm2

,/---L__

‘s___-__

:02 IY is -5

_.. ./

D

500 TIME

1000 BEFORE

1500 PRESENT(Y)

2000

2MO CONDUCTIVITY

[Wm.‘K-‘1

Fig. 2. A-C. Inversion results for synthetic examples. The 5-K temperature step is located at 500, 1000 and 1500 yrs B.P. for Figs. 2 A, B and C, respectively. For further explanation, see text.

350

lithology

and porosity

logs

may

and

help

conductivity-depth

may be known in constraining profile

well known

W m-i

K-’

between

which

amplitudes

and was assigned

at the

with a-priori They

conductivity-depth K-‘.

temperature constrained.

confidence are shown

limits

is constant

The amplitudes

simply

weakly

be noted:

of kO.1 in the 1 and

and equal to 3

not

ground

temperature

likely

to

temperature

curve

narrow

separation

that past climatic

changes

indirect

evidence

have occurred

show

and may

indicate how an a-priori climatic model with confidence limits should be chosen. The circles in Fig. 2 show the locations and values of the free amplitudes of the a-priori climatic interpolates linearly between and

is constant

prior

model. The climate the free amplitudes

to 2400 yrs B.P. It is as-

the plot.

reproduce

the

The

the tempera-

i.e., the continuous circles

indicating

conductivity-depth

changed,

profile

as expected

the is not

from its a-priori

limits.

For all cases, the estimated

curacy

other

in

interfor the

climatic

model is a

rather blurred version of the true step temperature change, with a generally decreasing statistical ac-

ice and

and

limits

observations.

confidence

at the measurement site. Tree ring records, pollen analysis, relative abundances of oxygen isotopes in snow

models

through

The thermal significantly

and the

of the three cases may

data satisfactorily,

passes

temperature

disappear

features

The estimated ture-depth

The confidence

profile general

1000

The continuous

the A 1 SD confidence

val of this estimate. following

of the past surface

are generally

to the effective

lines indicate

the

variations are generally a-priori poorly Meteorological records only exist for

the last 100 yrs and

1500 yrs B.P., respectively.

at 500,

broken

1 km,

as circles

located

to be

Only

plots in Fig. 2. Between

10 km, the conductivity

relate

linearly.

and

change

lines show the most likely model estimates

the value 3

0, 0.5 and

it interpolates

K-‘.

W m-i

depths

5-K step temperature

cases,

was assumed

at 0 and 0.5 km were assumed

unknown W m-i

the a-priori

profile. For the synthetic

the conductivity-depth a-priori

from well

oldest

for increasing of the

climatic

well resolved. tively sensitive

age (visible broken

amplitude,

in the increasing

lines

envelope).

however,

The

is relatively

This is because the data are relato this amplitude as it controls the

climate for times > 2000 yrs B.P. For comparison, the amplitude at 1600 yrs B.P. influences the climate only from 1200 to 2000 yrs B.P., and has a much less visible effect in the data and is, there-

sumed that the present-day temperature is relatively well known, with an a-priori value of 7 f

fore, less well resolved. While all three inferred climatic changes are far from step-like, the figure

1” C. Furthermore, it is assumed dence indicates a steady increase

that other eviin the average

illustrates how the more remote climatic events are the less well resolved. This is an inevitable conse-

over the past 2500 yrs. The (at 2400 yrs B.P.) is assigned

quence of the nature of diffusion, the finite amount and limited precision of the data, and the finite

ground temperature oldest free amplitude

the a-priori value of 3 + 4” C; the remaining amplitudes are given values of 6, 5, 5, 4 and 3 & 4°C. For Fig. 2A, two additional amplitudes at 200 and 600 yrs B.P. were introduced to accommodate the relatively recent temperature change at 500 yrs B.P. Initially, the steady-state background heat flow is poorly known but is expected to be in the range 50 + 50 mW m-*. To summarize, the a-priori information indicates that the thermal properties structure known and that the inverse procedure

is well should

primarily modify the surface temperature history and the background heat flow in order to satisfy the temperature data. Figures 2A-C show the inversion results from a

depth

of the borehole.

Short-period

surface

tem-

perature variations, which contain the bulk of information on when the temperature change occurred, are effectively suppressed by the low-pass filtering effect of the earth, while longer periods may suffer aliasing due to the finite length of the borehole. The net result is that shallow temperature-depth profiles relatively quickly “forget” precisely what occurred in terms of surface temperature variations. The shape of the climatic confidence envelopes is a consequence of the linear interpolation of the climate between the free climatic amplitudes and their statistical correlation. Let f(t) =fi . (1 - t) + j2 . t, tc(O, 1) be the interpolated climate between

351

the free amplitudes fi and f2 with standard deviations ui and a, and covariance cov( fi, f2) = a,~,~. The standard deviation of the interpolated climate for tr(0, 1) then reads: SD(f(N = /UT. (1 - t)2 + u; . t2 + 2U,U,~~ (1 - t) * t (5)

The optimum estimate (continuous line in Fig. 2), together with eqn. (5) defines the broken line confidence envelope. In general, $ is negative, indicating that changes in opposite directions of neighbouring climatic amplitudes that maintain the linear connection between the amplitudes within the confidence envelopes will have little influence on the predicted temperature-depth profile. The width of the error envelope emphasizes the high degree of equivalence in the inverse climatic problem. The estimated background heat flow density values with confidence limits are listed in Table 1. Although initially poorly constrained, they are well resolved for all cases. Repeated runs of each case show that the estimated heat flow density typically falls in the interval 40 k 5 mW rne2, depending on the actual realizations of the random data errors. The uncertainty of the estimated heat flow density shows little variation, because it is dominated by the uncertainty of the climatic amplitude at 2400 yrs B.P. Several runs with the incorrectvaluesofA=O~Wm-3andA=2~W m -3 show that the results are not sensitive to the heat production over a data interval of 600 m. Figure 2A also shows how the a-priori climatic amplitudes at 1200, 1600 and 2000 yrs B.P. influence the a-posteriori climatic estimate. The aposteriori estimate converges to the a-priori value because the climate in this time interval is not well

TABLE 1 Synthetic cases Case

Step change

heat flow density

(from Fig. 2)

(yrs B.P.)

(mW m-‘)

A

5cil

42k5

B

1000

41*5

C

1500

39+5

constrained by the temperature data, resulting in an overestimation of the surface temperature. This is compensated in the overall picture by an underestimate of the surface temperature at 800 yrs B.P. Several runs with an additional free conductivity amplitude at the depth of 250 m and a relaxed a-priori conductivity constraint of +0.5 W m-i K-’ show indeed how not everything is permitted to be loosely constrained. A stable solution which satisfies the temperature data is produced, but as a compromise of surface temperature variations and thermal conductivity-depth variations. The reduced near-surface temperature gradient is interpreted partly in terms of a higher thermal conductivity and partly in terms of an increase in surface temperature. It appears that surface temperature variations and thermal conductivitydepth variations to some extent may substitute for each other and that a necessary condition for obtaining more reliable information on past climatic event from borehole temperatures is strong control on the thermal properties structure.

Four real cases from the Canadian Shield

In the summer of 1985, a continuous temperature logging tool (Conaway, 1977; Conaway and Beck, 1977a, b) was used to obtain very accurate temperature-depth profiles in four shallow holes in the Superior Province of the Canadian Shield (Fig. 3). The holes were drilled for heat flow purposes in the winter of 1967-1968 by the Dominion Observatory, Ottawa, as part of the Upper Mantle Project. The obtained temperature logs, thermal conductivities and heat production values were described in Cermak and Jessop (1971) and Jessop and Lewis (1978). The continuous curves in Fig. 4 show the continuous temperature logs. The general features are oscillatory disturbances to a maximum depth of - 25 m due to the yearly temperature cycle, a pronounced curvature between - 25 m and - 250 m (except for Otoskwin River), in two cases resulting in a negative temperature gradient, and a relatively constant gradient from - 250 m to the bottom of the holes. The circles in Fig. 4 show, on the mK scale, the difference between the continu-

352

- 1000 yrs, which

past

observed

temperature

that a prominent

cooling

Ice Age is required Recent

events

Age may produce heat

as the Little

Ice

significant

and

experimentally

A more

such

as the temperature

remote

Glaciation

climatic

change,

- 10,000

nature

because

the heat flow is almost Fig. 3. Location Canadian

of boreholes

Shield. Geographical

W): Hearst--49O 97’

07.9’;

in the Superior

41.4’,

Minchin

coordinates

83”

Lake-50°

River-51’

32.1’;

Province

are (Lat. N, Long.

Cochrane

42.7’, 49.5’, 89’

of the

90”

-49O 28.8’;

06.2’,

upper

however, of the

yrs B.P., is not detectable

trace

its perturbation

constant

with depth

to in the

600 m.

Early depth

Otoskwin

heat

interval

ature gradient

35.9’.

in a 600-m

rise at the end

likely to leave any experimentally of its transient

the

he finds

the data.

such

flow perturbations

hole.

Wisconsin

with

at the time of the Little

to satisfy

climatic

unequivocal

are consistent

logs. In particular,

flow

workers

typically

of a reasonably

selected

constant

a

temper-

for heat flow determination,

in this

way avoiding the influence of recent surface temperature changes. The observed heat flow value ous temperature logs and the discrete logs obtained 16 years earlier. Over the constant gradient interval, the difference logs generally lie in the

would then be corrected for the assumed influence of more remote climatic changes. With the inverse approach it is possible to utilize the total extent of

interval

the temperature log for a heat flow density determination and take into consideration the in-

f25

mK and show a smooth

trend

with

depth. In one case (Minchin Lake), the difference changes from 0 mK at - 250 m to - 50 mK at - 600 m. In the depth

interval

of pronounced

curvature, the differences show a bell-shaped feature with a maximum amplitude of 100-200 mK at a depth of - 100 m. This is not visible for Otoskwin River. Clearly, the two temperature gradient librium. ductivity,

profiles

showing

a

inversion cannot be in conductive equiFurthermore, high values of thermal conwhich have not been observed,

would be

required to reduce the near-surface gradient remaining holes if they were in thermal

fluence

of climatic

the recent climatic

in the equi-

librium. Because, additionally, there is no evidence of water flow, it is natural to interpret the bellshaped feature as the result of the conductive decay of a transient temperature disturbance over the - 16 years that have elapsed since the first logging of the holes. This interpretation is consistent with the work of Cermak (1971). For two holes, one of which (Hearst) is shared by the present study, he interprets the curvature as the result of a general climatic warming over the last - 200 yrs and constructs surface temperature histories for the

changes events,

on any time scale. For the transient

remains

of

which are still visible in the temperature log, the requirement of temperature data satisfaction will cause a modification of the a-priori assigned recent temperature history (if the a-priori assigned standard

deviations

are not too small). The a-priori

assigned remote temperature history, however, will remain unchanged and act as a correction because modifications

are not justified

role

surface

of

the

by the data.

temperature

history

The thus

changes gradually with increasing age from being one of producing a fit to the temperature data to one of acting as a correction. Clearly, the inversion result in terms of heat flow density strongly depends on the choice of the corrective end of the surface temperature history. This has long been acknowledged and great care has been exercised in constructing plausible climatic models for use in heat flow corrections. Cermak and Jessop (19710 arrived at an estimated 4-K temperature rise at the end of the Wisconsin Glaciation by taking the average ground temperature to be close to the melting point of ice under

353

HEARST

COCHRANE 5

0

0 O-

10

15

00

loo-

500 -

600-100

1 100

MINCHIN 5

0

0

I1 0

1

I

I 200

I1

300

OTOSKWIN

LAKE 10 I

-

15

RIVER

0

*

100

500

\

(

600 100

Fig. 4. Continuous between

curves

I 0

I

I 100

show the continuous

the continuous

I

I 200

I

I 300

temperature

logs and the discrete

,,

1oc

logs obtained

logs obtained

0I

I 1,

100

in the summer

,

200 ,

,

300 ,(rnK)

of 1985. Circles

16 years earlier. The difference

indicate

the difference

log refers to the mK scale.

the pressure at the base of the ice cover. Beck (1977) proposed a 12-K temperature rise in particular based on the 1sO/‘6O ratio variations in

With a view to the data analysis technique presented here we let the still unresolved speculations on the average ground temperature at the

ice cores from Camp Century, Greenland (Dansgaard et al., 1971) and the actual observed temperature at the base of the Greenland ice shield. In terms of paleoclimatic corrections to terrestrial heat flow density observations in the shield holes, the difference between the two models amounts to - 10 mW mP2.

base of the ice sheet be reflected in the a-priori information on that parameter. To accommodate both models we take the temperature rise at the end of the glaciation to be 8 K, with a standard deviation of 4 K. This reflects our belief that the most likely temperature rise is 8 K and that there is only one chance in three that the actual temper-

354

TABLE

2

A-priori

and a-posteriori

Time (yrs B.P.)

climate

for the Canadian

Shield holes

A priori

A posteriori

Amplitude

Hearst

Cochrane

(“C)

Minchin

Otoskwin

Lake

River

0

2

*4

5.2 kO.1

5.9 +0.1

4.4

+0.1

50

2

*4

6.2 kO.2

5.4 +0.2

4.5

+0.2

2.9 +O.l 2.6 *0.2

100

2

*4

5.3 kO.4

3.5 +0.4

3.7 kO.3

2.2 +0.3

150

2

+4

3.3 *0.7

1.6 +0.6

2.6 +0.6

2.0 +0.5

200

2

*4

1.8 +0.8

1.3 f0.8

2.3 +0.8

2.0 IfIO.8

250

2

*4

1.4 *1.2

2.3

k1.3

2.7 51.2

2.3 +l.l

300

2

+4

1.8 51.5

3.4 +1.6

3.3 +1.6

2.5 k1.5

350

2

*4

2.2 f1.9

3.9 k1.9

3.6 + 1.9

2.6 *1.8

400

2

+4

2.4 It2.2

3.6 +2.2

3.6 +2.2

2.6 f2.2

500

2

&4

2.3 f2.6

2.5 +2.7

3.3 *2.6

2.4 +2.5

600

2

*4

2.2 f3.1

1.9 f3.2

3.1 f3.1

2.3 f3.0

800

2

*4

2.1 k3.5

1.8 +3.5

3.3 +3.5

2.5 k3.4 2.8 *3.3

1000

2

*4

2.1 f3.4

2.0 *3.4

3.9 k3.4

1500

2

+4

2.1

2.3 k3.5

4.2

2000

2.25 f 4

2.34 f 3.7

2.62 f 3.8

3.97 + 3.7

3.16 f 3.6

2500

3.75 f 1

3.76 + 1.0

3.82 + 1.0

3.97 f 1.0

3.89 f 1.0

8000

3.75 f 1

3.76 rt 1.0

3.78 f 1.0

3.86 + 1.0

3.82 f 1.0

12000

- 4.25 f 4

- 4.21 f 3.9

- 4.05 + 3.9

- 3.32 f 3.9

- 3.82 f 3.7

60000

-4.25

- 4.24 + 4.0

- 4.24 f 4.0

- 4.05 * 4.0

- 4.21 f 4.0

0.00 f 2.0

- 0.03 f 2.0

0.03 * 2.0

- 0.04 f 2.0

70000

* 4

O.OOf2

The climate interpolates

linearly

between

k3.4

the specified

amplitudes.

+3.4

Prior to 70,000 yrs B.P., the surface

3.1 k3.3

temperature

was constant

and

equal to the value at 70,000 yrs B.P.

ature rise lies outside the interval 4-12 K. The a-priori surface temperature amplitudes with confidence limits are listed in Table 2. It is assumed that prior to 70,000 yrs B.P. the climatic variations may be averaged to a constant which we believe is

reduced 340-390 gradient TABLE

3

most likely zero. This a-priori

model

Thermal

properties

braces

and Jessop

the models

of Cermak

and Beck (1977). The a-priori thermal

conductivities

are the values used by Cermak and and Jessop and Lewis (1978) with a dence limit for cases A, B, and D thermal conductivity is allowed to

largely

em-

temperature gradient in the interval m and a further decreasing temperature towards the bottom of the hole. This is

(1971)

(Table

3)

Jessop (1971) + 10% confi(Fig. 2). The vary linearly

with depth between amplitudes at depths of 0, 0.5 and 1 km, the upper two of which are free parameters. The model extends to Z,, = 10 km with a constant conductivity equal to the a-priori assigned conductivity between 1 and 10 km. However, case C (Minchin Lake) appears to have a more complicated conductivity-depth profile. The temperature gradient log of this hole shows a

Hearst

Conductivity

Heat production

(W m-l

(PW m-?

K-l)

A priori

A posteriori

Fixed

3.14*0.31

3.04kO.27

1.8

3.14 f 0.29 Cochrane

2.51 f 0.25

2.42 + 0.22

1.3

2.53 4 0.24 Minchin Lake

Variable

2.89 + 0.27

0.86

3.26 + 0.31 Otoskwin River

2.69 f 0.27

2.67 + 0.26

0.18

2.70 f 0.27 The two a-posteriori conductivity

with

pc = 2.5.106

J me3

conductivity uncertainties K-r

values indicate at depths

for all cases.

the estimated

of 0 and

0.5 km.

355

accommodated

by increasing

tivity between

340-390

the conductivity bottom

m and by further

allowing

linearly

of the hole. With further is linearly

value of 3 W m-l

at 1 km and dynamic

conduc-

to increase

the conductivity surface

the thermal

maintained

heat capacity

be constant

K-‘,

is constant

and Jessop

to

(Beck, for each

(1971) and

3 summarizes

the

B

A

flow estimates Hearst

52

57.6 k6.0 (10.4%)

57.6 f 3.9 (6.6%)

Cochrane

43

45.3 f 5.4 (11.9%)

45.5 f 3.1 (6.8%)

42

50.8 f 5.8 (11.4%)

51.9+ 3.3 (6.3%)

25

29.4+4.5

29.8 + 1.8 (6.2%)

Minchin Lake Otoskwin River

(15.3%)

information.

The results of the inversions Figs. 5A-D. The circles indicate

are displayed in a-priori values of

data and parameters, while the continuous lines are model estimates. Broken lines indicate the k 1 SD confidence limits of the most likely model. From

Present heat flow estimates

heat

to its

= 10 km. The

(mW m-*)

Previous

depths

back

with pc = 2.5. lo6 J rnp3 K-’

properties

Heat flow densities

which is reached

to Z,,

and Lewis (1978). Table

4

the

for all cases is assumed

case, as given by Cermak thermal

increasing

brought

1977) and the heat production Jessop

towards

TABLE

the temperature-depth

plots

it follows

Glaciation. Furthermore, since the present climatic model contains the model of Cermak and Jessop (1971) within its confidence limits, it appears that the present

heat

flow estimates

with their confi-

that the estimated models satisfy the temperature data. The conductivity-depth plots show some, but very little, departure from the a-priori data

dence limits generally embrace the previous estimates. The Minchin Lake case, however, shows a significantly higher heat flow density value. This is

(Table 3). The estimated (Table 2) display

due to our choice of a conductivity-depth profile showing greater conductivity values towards the bottom of the hole in order to accommodate a

surface temperature similar characteristics.

histories Adjust-

ments have occurred over the last 1000-1500 yrs B.P., with greatly reduced standard deviations back to - 500 yrs B.P. The estimates all show a statistically significant low coinciding with the historical climatic event known as the Little Ice Age. Prior to 1500 yrs B.P., the a-posterior-i amplitude estimates are practically identical to the a-priori assigned amplitudes, and the a-priori standard deviations have not been reduced significantly. This means that the effects of surface temperature variations prior to - 1500 yrs B.P. are not essential in explaining the curvature in the present temperature logs and that this part of the climatic model acts more like a correction to the data. Column A in Table 4 shows the resulting heat flow density estimates with absolute confidence limits and a percentage error. Also shown are the accepted heat flow density values from Cermak and Jessop (1971) and Jessop and Lewis (1978). It may be seen how the present estimates are consistently higher than the previous estimates. This is due to the choice of 8 K rather than 4 K for the temperature rise at the end of the Wisconsin

reduced pret

temperature

this hole

vity-depth estimates.

gradient.

in terms

profile Column

Attempts

of a constant

resulted in unrealistic B in Table 4 shows

to interconducticlimatic for com-

parison the result of performing the analysis with zero a-priori uncertainty on the climatic amplitudes of t 2 8000 yrs B.P. in Table 3. While the heat flow estimates are only slightly adjusted, the a-posterior-i standard deviations of the estimates are significantly reduced. This emphasizes the importance of accurate knowledge on the climatic warming at the end of the Wisconsin heat flow density determinations. The results

show how realistic

Glaciation uncertainties

to in

the observed thermal parameters, such as those of Fig. 5, may imply significant error bounds on the heat flow density value obtained. How well heat flow density differences between holes are resolved depends on the lateral variation in the effective ground temperature history. If laterally homogeneous, heat flow density differences are likely to be real rather than caused by geographical variations in exposure.

A: HEARST q =

57.6!:6.0mWmm2

.7 2

TEMPERATURE

2s

3

CONDUCTIVITY

(‘0

3.5

‘

(Wm-‘K-‘)

zoo0

500

TIME

ZOO

BEFORE

10000

PRESENT

(Yl

6: COCHRANE q=45.325.4

mWme2

.6

I 7 i-

0

5

10

TEMPERATURE

Fig. 5. A-D.

15

(“Cl

2

25

CONDUCTIVITY

3

35

‘

(Wm.‘K-l)

Inversion results for the four holes from the Superior Province of the Canadian Shield. A. Hearst. B. Cochrane. Minchin Lake. D. Otoskwin River. For further explanation, see text.

C

357

1OOD

500

1500

2000

TIME

2500

BEFORE

l0000

PRESENT

30000

soooo

70000

(Y)

a .I

C:MINCHIN

.2

LAKE

q= 50.8z5.8mWme2

.3 .L .5-

.7 -0

5

IO

TEMPERATURE

.7' 2

15

(“Cl

2.5

CONDUCTIVITY

.“I' 3 3.5

L

(Wm-‘K-‘1

I’

L________---------------_____

, -----_____*,

500

moo

2000

1500

TIME

2500

BEFORE

moo

PRESENT

30000

0: OTOSKWIN q=29.4~4.SmWm-2

0

5

TEMPERATURE

10

15

(“~3

CONDUCTIVITY

(Wm-‘K-9

50000

(Y)

RIVER

70000

358

Conclusions

The problem of obtaining a reliable heat flow density value from observations of borehole thermal properties and indirect evidence on past climatic changes has been treated within the framework of stochastic inverse theory. The forward problem is solved by using a onedimensional finite element discretization in space and the Crank-Nicholson time integration algorithm. This allows complex thermal parameter-depth profiles, as may be found in sedimentary environments, to be represented. Furthermore, restrictions essential to certain analytical forward methods such as periodicity of the surface temperature history or the representation of the subsurface by elementary layers are avoided. The inverse method allows for free parameters of any kind. The present implementation allows adjustment of surface temperature history, thermal conductivity structure and background heat flow. These are the principal parameters that influence the temperature-depth profile and they are sufficient for explaining most conductive onedimensional temperature-depth profiles. Each parameter of the problem is assigned an a-priori value with a standard deviation that reflects the confidence we have in that value. During the inversion, the parameters are adjusted with a view to the a-priori constraints and the data sensitivity in order to reproduce the temperature data within their statistical accuracy. The resulting aposteriori parameter estimates with standard deviations are the most likely values of the parameters, given the a-priori information on data and parameters and the mathematical structure of the theory. The method is robust in the sense that a superfluous parameter or almost proportionality between different data-parameter relationships does not cause instability. When a parameter cannot be resolved from the data, its a-posteriori value is guided by the a-priori value and possible a-priori covariances or theoretical relationships with better resolved parameters. In particular, this means that climatic events which bestowed no transient remains on a given temperature log and therefore are inherently merged with the steady-state thermal structure we

wish to resolve may still be included in the inverse analysis. The a-priori temperature-time structure of the event, which is estimated from indirect evidence, then acts as a correction to the temperature data while the a-priori uncertainty propagates into uncertainty on, for example, the steady-state background heat flow. This quality of the inverse technique makes it a valuable tool for constructing physically reasonable models that satisfy the data and for investigating quantitatively how uncertainties in model parameters and data affect the principal parameters of interest. The technique has been applied to four temperature profiles obtained with a continuous logging device from the Superior Province of the Canadian Shield. All but one of the temperature logs exhibit a conspicuous curvature in the upper 200 m. Comparison with temperature logs taken 16 years earlier show a characteristic difference over the depth interval of pronounced curvature. The difference is consistent with conductive relaxation of the curved temperature profiles and substantiates the transient nature of the shallow geothermal field. Inversion results show how thermal properties from 600-m deep holes may simultaneously constrain the surface temperature history for the past 500 years and the terrestrial heat flow density. In all cases, a statistically significant low at the time of the Little Ice Age is required to satisfy the data. In two cases (Hearst and Cochrane), the estimated change in average ground temperature over the past 300 years is 5 K. This may in part be explained by changes in surface characteristics due to the felling of trees. Additionally, it was demonstrated how uncertainty concerning the temperature rise at the end of the Wisconsin Glaciation may be included in the analysis and how it may affect the accuracy of heat flow density determinations. Acknowledgements

We gratefully acknowledge the Danish Natural Science Research Council, The World University Service of Canada and the Natural Science and Engineering Research Council of Canada for financial support, K. Wang for help with the

359

equipment and during the field trip, Dr. P.Y. Shen and Dr. N. Balling for helpful discussions and Dr. A.M. Jessop and Dr. M.J. Drury for locating the holes and for supplying thermal properties information.

Jessop, A.M. and Lewis, T., 1978. Heat flow and heat generation in the Superior Tectonophysics,

Province

of the Canadian

Mansure, A.J. and Reiter, M., 1979. A vertical ground water movement correction

for heat flow. J. Geophys.

and applications.

and their effect

on regional

means. Tectonophysics,

and continental

heat-flow

41: 17-39.

Beck, A.E., 1982. Precision logging of temperature and the extraction

of past climate.

gradients

Tectonophysics,

83:

1-11. Becker, E.B., Carey, G.F. and Oden, J.I., 1981. Finite elements. Vol., 1, an Introduction.

Prentice-Hall,

Englewood Cliffs,

N.J., 258 pp. and Oden, J.I.,

Computational

Aspects.

1984. Finite Prentice-Hall,

Elements.

Vol., 3,

Englewood Cliffs,

N.J., 350 pp. Carslaw, H.S. and Jaeger, J.C., 1959. Conduction

of Heat in

Solids. Oxford Univ. Press, 510 pp. term&k,

V., 1971.

climatic

Underground

temperature

Palaeoclimatol.,

temperature

of the past

Palaeoecol.,

and inferred

millenium,

Palaeogr.,

10: 1-19.

Cermak, V. and Jessop, A.M., 1971. Heat flow, heat generation and crustal temperature

in the Kapuskasing

Canadian Shield. Tectonophysics, Conaway, J.G.,

area of the

11: 287-303.

1977. Deconvolution

of temperature gradient

logs. Geophysics, 42: 823-837. Conaway, J.G. and Beck, A.E., 1977a. Continuous

logging of

temperature gradients. In: A.M. Jessop (Editor), Heat Flow and Geodynamics.

Tectonophysics,

41: l-7.

Conaway, J.G. and Beck, A.E., 1977b. Fine-scale between

Ph.D.

temperature log: Method

Thesis,

Univ. Western

Ontario,

London, Ont.

Beck, A.E., 1977. Climatically perturbed temperature gradients

Carey, G.F.

Res., 84

(B7): 3490-3496. Nielsen, S.B., 1986. The continuous

References

Shield,

50: 55-57.

temperature

gradient

correlation

logs and lithology.

Geo-

physics, 42: 140-1410. Dansgaard, W., Johnsen, S.J., Clausen, H.B. and Langway Jr., C.C., 1971. Climate record revealed by the Camp Century ice core. In: K.K.

Turekian (Editor),

The Late Cenozoic

Glacial Ages. Yale Univ. Press, pp. 37-56.

Nielsen, S.B. and Balling, N., 1985. Transient stratified medium. Tectonophysics, Powell, W.G.,

Chapman,

1988. Continental

D.S.,

heat-flow

heat flow in a

121: l-10.

Balling,

N. and Beck, A.E.,

density.

In:

R. Haenel,

L.

Rybach, and L. Stegena (Editors), Handbook of Continental Heat Flow Density Determination,

Riedel, Dordrecht,

pp. 167-122. Shen, P.Y. and Beck, A.E., temperature

1983. Determination

of surface

history from borehole temperature

gradients.

J. Geophys. Res., 88: 7485-7493. Tarantola,

A., 1987. Inverse Problem Theory-Methods

Data Fitting and Parameter Estimation.

for

Elsevier, Amster-

dam, 613 pp. Tarantola,

A. and Valette, B., 1982a. Inverse problems-quest

for information. Tarantola,

J. Geophys., 50: 1599170.

A. and Valette

B., 1982b. Generalized

nonlinear

inverse problems solved using the least squares criterion. Rev. Geophys. Space Phys., 20: 219-232. Vasseur, G., Bernard, Ph., Van de Meulebrouck,

J., Kast, Y.

and Jolivet, J., 1983. Holocene palcotemperatures from geothermal measurements.

deduced

Palaeogr., Palaeoclimatol.,

Palaeoecol., 43: 237-259. Vasseur, G., Lucazeau, F. and Bayer, R., 1985. The problem of heat

flow density

Tectonophysics,

determination

from

inaccurate

data.

121: 25-34.

Wang, K. and Beck, A.E., 1987. Heat flow measurements

in

lacustrine or oceanic sediments without recording bottom temperature

variations.

J.

Geophys.

Res.,

92

(B12):

12837-12845. Wang, K., Shen, P.Y. and Beck, A.E., 1986. On the effects of thermal properties structure and water bottom temperature variation on temperature gradients in lake sediments. Can. J. Earth Sci.. 23: 1257-1264.