Holocene paleotemperatures deduced from geothermal measurements

Palaeogeography, Palaeoclimatology, Palaeoecology, 43 (1983): 2 3 7 - - 2 5 9

237

Elsevier Science Publishers B.V., A m s t e r d a m -- Printed in The Netherlands

HOLOCENE PALEOTEMPERATURESDEDUCED FROM GEOTHERMAL MEASUREMENTS

G. V A S S E U R l, Ph. B E R N A R D 1, J. V A N DE M E U L E B R O U C K l, Y. K A S T 1 and J. J O L I V E T 2

1Centre G~ologique et G~ophysique, Universit~ des Sciences et Techniques du Languedoc, Place E. Bataillon, 34060 Montpellier Cedex (France) :Institut de Physique du Globe, 4, Place Jussieu, 75230 Paris Cedex 05 (France) (Received and accepted April 28, 1983)

ABSTRACT Vasseur, G., Bernard, Ph., Van de Meulebrouck, J., Kast, Y. and Jolivet, J., 1983. H o l o c e n e paleotemperatures deduced f r o m geothermal measurements. Palaeogeogr., Palaeoclimatol., Palaeoecol., 43 : 237--259. Detailed measurements o f t e m p e r a t u r e and rock thermal properties were carried o u t in a 1000 m deep borehole drilled in massive granite. Using a suitable reduction, these measurements reveal a c o h e r e n t variation of apparent heat flow with depth, which is attributed to past variations of the ground t e m p e r a t u r e during the last 10,000 years. T w o inversion techniques of the data are described in order to reconstruct these paleotemperatures. A classical inversion m e t h o d allows a discussion of the i n f o r m a t i o n c o n t e n t o f the data, but leads to unreasonable t e m p e r a t u r e variations. A n o t h e r m e t h o d using b o u n d s on p a l e o t e m p e r a t u r e variations provides m o r e acceptable t e m p e r a t u r e models. The main features o f these models are the presence of a cold event (identified as the Little Ice Age) and of warm episodes around 1500 BC and at the beginning of the present century. The amplitude of the inferred variations would be at least 1.09°C. INTRODUCTION

Due to the development of heat flow studies, considerable attention has been paid in recent years to temperature measurements in boreholes. From a physical point of view, the ground generally behaves as a semi-infinite conducting medium. Its surface temperature is imposed b y external phenomena and stands for an upper boundary condition; if thermal equilibrium was to occur, the temperature variation with depth would only reflect internal heat process: internal heat flow, heat production, etc. However, any variation with time of the surface temperature diffuses slowly downward while attenuating and therefore modifies the internal temperature field. For example, seasonal surface temperature variations of large amplitude can penetrate in the upper levels of the ground: they are characterized b y an exponential decrease of amplitude with depth with a scale depth of several meters, and below a few tens of meters their effect becomes negligible. However, 0031-0182/83/$03.00

© 1983 Elsevier Science Publishers B.V.

238

temperature variations with a longer period can penetrate deeper; it can be shown (Carslaw and Jaeger, 1959) that variations of the surface temperature with a period of some tens to several thousands of years remain recorded in the ground at a depth of a few tens to a few hundreds of meters. The importance of this effect on temperature measurements in boreholes was soon recognized (Birch, 1948; Cermak, 1971); in particular the possible influence of the last glacial epoch during the Pleistocene. Indeed, using numerous indirect data, it is possible to state that the world climate has not remained unchanged during the past 100,000 years (see for example Lamb, 1977). During the Pleistocene it is possible to distinguish several glacial and interglacial epochs with a time scale of several tens of thousands of years. The last glacial episode (Wfirm in the alpine.domain), lasting roughly from 70,000 to 10,000 yr. B.P., is characterized by a mean annual temperature lower than the present one by some 10°C. During the past glacial epoch (Holocene), variations of this mean annual temperature also occurred; however the amplitude of these variations is much smaller (some tenths of °C to a few °C) and their actual geographical and temporal extents are not well documented. The fact that the corresponding variations of the mean surface temperature are recorded in the subsurface suggests a possible use of geothermal data to solve an inverse problem, i.e. to reconstruct past surface temperatures from measurements of the present temperature as a function of depth. It is the purpose of the present paper to apply such an inversion procedure to good quality data collected in a 1000 m borehole drilled (with core) in a granitic area. Using these measurements a systematic depth variation of the temperature gradient and of the apparent heat flow is evidenced; this variation can quite certainly be attributed to the paleoclimatic effect and can therefore be used as basic data for solving the inverse problem. Several authors have already attempted this kind of approach (Beck and Judge, 1969; Cermak, 1971): for example, Cermak (1971) used measurements in several 600 m deep boreholes of the Canadian shield to propose a paleotemperature model for the past millennium. In our case, however, we make use of recent progress in the inversion techniques which are now c o m m o n l y used in geophysics (Jordan, 1973; Sabatier, 1977a). The data that we obtained are processed by two complementary techniques which enable us to point out the information contained in the data, and to c o m p u t e realistic temperature models. DATA DESCRIPTION

Temperature data An exploration borehole has been drilled in the granitic massif of Auriat near Limoges, reaching a depth o f 1000 m (Fig.l). Various geological studies (Duthou, 1977; Chenevoy, 1958) have shown that this granitic batholith,

239

i,

.

f]/l

Fig.1. Geographical location of the Auriat borehole. Hatched areas are b a s e m e n t outcrops.

20 km wide at outcrop, is of Late Hercynian (330 m.y.) age and characterized b y a homogeneous porphyroidic monzogranite. The borehole, entirely cored, was drilled from March to June 1980, and various measurements and tests have been performed. The thermal logging was executed in December 1980; at that time the borehole had been at rest for several weeks. The temperature was measured continuously during descent, using a low thermal capacity probe with an absolute accuracy of 0.01°C and a sensitivity of 0.001 ° C. The resulting temperature profile as a function of depth is plotted in Fig.2. Below 20 m (depth of the water level in the borehole) a very regular increase of temperature is observed; however, when substracting the main linear trend, slight fluctuations are obtained. As will be shown later, these fluctuations are n o t necessarily a signature of paleoclimatic variations. The field of temperature in the ground depends on many factors and in particular on the thermal conductivity, K. Inhomogeneities of conductivity result in a modification of the temperature gradient. The classical way to remove this effect is to evaluate the vertical heat flow which is defined as: q =g~O/~z

(1)

240 10

20

3o

~ (-c)

0

,,

.

200

400

6OO

a)

\

I

b)

800

z(m)

- 0,5

0

0,5

°C

Fig.2. a. T e m p e r a t u r e profile as a f u n c t i o n of d e p t h , b. T e m p e r a t u r e d e v i a t i o n s t o a linear

law (0 = 10--0.029z).

Thermal conductivity measurements and heat flow evaluation Cores were sampled every 10 m for thermal conductivity measurements. Two laboratory techniques were used: the standard divided-bar equipment and a transient conductivity probe (Kast and Gallice, 1980). The measured conductivities K are plotted in Fig.3b as a function of depth. In Fig.3a, b conductivity measurements can be compared with temperature gradients obtained at the same level. The gradient is calculated as the slope of the straight line fitting in the least square sense the temperatures observed at depth intervals of 20 m around each conductivity measurement. Fig.3c displays the heat flow estimates obtained as the product of the temperature gradient and the conductivity.

241

a) i

2O

i

i

25

3O

r

W m 1 "C _1

b)

°Cm-1 r

i

i

i

i

i

i

i

~ i

]

I

i

[

:

,

i

3

!

~

i

c) t

4

i

t

80

100

m W m -2 I

120

,/

100

j~

=--r

200

-

300

400 J

I

500

5--

600

?

700

800

/

k

I

(

,t

900 2

/

1000

Zm

Z

Z

Fig.3. Variation as a f u n c t i o n o f depth o f : a, the temperature gradient o n 2 0 - m intervals, a e / a z ; b, the thermal c o n d u c t i v i t y , K; c, the heat f l o w , q = K a e / a z . H o r i z o n t a l bars indicate the standard deviations.

242 Other thermal parameters

Two other parameters are needed for this study: the thermal capacity and the heat production. Classical calorimetry experiments were used to obtain the specific heat of samples; a value of 2 • 106 J m -3 o C_1 was obtained. Volume heat production is mainly due to the decay of natural isotopes such as 23SU, 238U, 232Th and 4°K. The resulting heat production is given by Birch et al. (1968): A~wm-3 = 1.33" 1 0 - 6 d [ 0 . 7 ( U ) + 0 . 2 0 ( T h ) + 0.27(K)]

(2)

with (U) and (Th) in p.p.m., and (K) in percent, d being the density in kg m -3. The concentration of these elements was measured on ten samples yielding the following mean value: A = 5.8 + 0.7 pW m -3 Data reduction

In Fig.3 we observe small-amplitude variations of the temperature gradient (1 to 2 . 1 0 -3 °C m -~ amplitude or 3--6% of the mean value) and larger fluctuations of the conductivity (10--20% of the mean value). There is no correlation between gradient and conductivity fluctuations. In fact such conductivity fluctuations are mainly related to the small-scale heterogeneity of the measured samples; t h e y are n o t representative of the mean rock conductivity at the decametric scale. In order to obtain more significant quantities, it is possible to calculate a mean value of the various quantities (temperature gradient, conductivity, heat flow) in some depth interval; we choose an interval of 70 m which results in a mean over seven or eight individual quantities. In each interval, we obtain a mean value of each quantity, ~ = Z x l / n , and also its variance, s 2 = Z (xi -- ~)2/(n -- 1), and the variance of the mean, s2/n (Hamilton, 1964). These twelve mean quantities are plotted as a function of depth in Fig.4. They present a smooth variation which seems significant on account of the estimated variance. The temperature gradient increases whereas the mean conductivity decreases from 3.2--3.35 W m -~ °C -1 in the upper level to a smaller value, 3.0 W m -~ °C -1, below 700 m. The heat flow reflects these two variations: it increases regularly to 300 m, then decreases by 6 mW m -2 between 380 and 775 m and then again increases slowly. R ole o f the various disturbing p h e n o m e n a

As stated before, the temperature gradient is influenced by conductivity contrasts. In the absence o f any other p h e n o m e n o n than steady geothermal heat flow and in a horizontally stratified medium, the product K aO /az remains constant whereas the conductivity and the temperature gradient

243 • 025

.Cm-1

.03

I

T

,

J

2.5

i

,

a)

I

3.5

3 '

'

=

'

I

200

200

40O

400

600

600

800

800

\

80

90

I

I

I

100

Wm

-'1 " C -1 1

,::m') 110

'

L ....

b/

z(m) 70

4 . . . .

....

l

'

m Wm

5

-2

mWm-

l

200

200 I.

c~

400

400

600

600

all;, •800

800

z(m) z(m)

Fig.4. a. S m o o t h e d t e m p e r a t u r e g r a d i e n t profile over 70-m intervals, b. S m o o t h e d conductivity profile over 70-m intervals, c. S m o o t h e d h e a t flow profile over 70-m intervals, d. Heat flow c o r r e c t i o n s for t o p o g r a p h y (topoi) , h e a t p r o d u c t i o n (Wi) a n d last glaciation (Gli).

244

may vary with depth. The question arises to know whether the vertical variation of the mean conductivity (Fig.4b) can be assimilated to a horizontal stratification. No information concerning lateral variations is available. However, it is geologically reasonable to assume that the mean thermal properties of the medium have horizontal gradients much smaller than their vertical ones. Using this assumption, it is necessary to discuss in detail the various causes which could affect the vertical variations of the observed heat flow: (a) There is no hint for the existence of natural water movements. In fact, the rocks were observed to have a very low permeability (3 • 10 -'1 to 10 -'; m s-l).

(b) The effect o f erosion is certainly very small. The studied area belongs to a geological unit -- Massif Central -- which was affected by some uplift during the Villafranchien (Derruau, 1971). In the eastern part of the Massif Central, the thickness o f eroded layers is small (less than 100 m) and results in a negligible effect. (c) Around the borehole, the topography presents some undulations; the corresponding topographic effect can be estimated by the m e t h o d of Bullard (1940). The resulting correction on the heat flow (Fig.4d) is positive and decreases with depth from 4.5 mW m -2 at 100 m to 1.4 mW m -2 at 900 m. (d) The local heat production with rate A modifies the gradient. Assuming A to be constant (5.8 p W m-3), it results in a linear decrease with slope --A. In order to account for this effect -- i.e., to reduce the observed heat flow to its surface value -- the observed values must receive a positive correction of about 0.5 mW m -2 at 100 m to 5.3 mW m -2 at 900 m depth (Fig.4d). (e) The two previous effects result in a small variation which can be removed by the corrections described above. The total correction is almost constant with depth and the corrected heat flow profile still presents a systematic variation with depth. This variation is ascribed to paleoclimatic effects. THEORY OF T H E R M A L DISTURBANCE PROPAGATION

The effect o f a temperature disturbance generated at the surface of a homogeneous semi-infinite half medium is classically obtained as the solution of the one-dimensional heat equation. Let K and K be the conductivity and the diffusivity of the half medium at a constant initial temperature, ~ (z, To) = 0. If, starting from this initial time r0, an external disturbance ¢(r) [= 6 (0, r)] is imposed at the surface o f the medium, then, for r > r0, the resulting temperature inside the medium is given by (Carslaw and Jaeger, 1959): T

0(z, r) = ( z / 2 v ~ ) ~ ¢(X) exp [--z2/4K ( r -

~,)] [ r -

~]-3n dX

(3)

0

In this solution, no heat flow from internal origin and no distributed heat production are considered. In order to account for an initial heat flow qo

245

constant with time, it is necessary to add to eq. 3 the term qoz/K. Similarly, a supplementary term --Aoz2/2K accounts for the effect of a constant heat production (in this modified solution, the initial state is assumed to be that of equilibrium). Now let us count positively past time as scaled from the present epoch (i.e. 1980 A.D., which is the year of measurement) and shift the initial time to +oo. Then, performing suitable variable changes on eq. 3, we obtain the following relation between the present
(4)

o

By derivation of this quantity with respect to z, we obtain the temperature gradient: ~0 /3z = ( 1 / 2 ~ - ~ ) ; 0(t) exp (--z:/4K t)(1 -- z2/2Kt) t-3/:dt + (qo -- A o z ) / K o

(5) and the corresponding heat flow: q(z) = K ~O/~z = qo - - A o z + ( K / 2 x / ' ~ ) ; ~(t) exp (--z2/4Kt) o

X (1 --

z2/2Kt) t-3/2dt

(6)

Eqs. 4, 5 and 6 are valid only for a medium with homogeneous conductivity and diffusivity. In the case of a horizontally stratified medium (K and K functions of z), the temperature and gradient profiles are modified b y conductivity contrasts. This effect disappears on the steady-state heat-flow profile; however, the integral term of eq. 6 which accounts for surface temperature variations is no longer valid. In fact, using numerical simulations for a layered medium with conductivity variations of + 10%, it can be shown that conductivity variations have only secondary-order effects on this term. INVERSION WITHOUT BOUNDS

Analysis o f the inverse problem The inverse problem to be solved can be stated as follows: using eq. 6 find o u t ¢(t) -- i.e. the variation of paleotemperature -- from observations of q(z). This problem belongs to a class o f classical linear inverse problems. Assume that q0 is known; each observation of the apparent heat flow qi observed at level zi -- eventually corrected for other effects -- is related to the u n k n o w n function ~ (t) b y a linear relation:

246

7i = q ( z i ) - q o =

] Gi(t)¢(t)dt

fori =1,...,12

(7)

0

Gi(t) is the function defined by eq. 6 for z = zi. The general method of solving this type of problem, where the data are in finite number and the unknown is a function, was described by Backus and Gilbert (1967). However, due to the fact that in our case f~" Gt(t)dt = 0, the original method proposed by these authors cannot be applied here. Therefore we applied a slightly different inversion technique based on the same principles; this technique, developed by Courtillot et al. (1978; see also Sabatier, 1977a), allows an interesting discussion of the informational content of the data. In this method, a "projected" solution ¢(t) is investigated in the form of a linear combination of functions Gi(t): 12

¢(t) = Y+ aiVi(t)

(8)

i=l

From eqs. 7 and 8 the coefficients a~ can be obtained by solving the linear system:

7, = ~ ~i ; Vi(t) Gj(t)dt j=l

f o r i = 1 , . . . , 12

(9)

0

It is possible to show (Sabatier, 1977a) that this projected solution is also the solution with least L~ norm [/~" ~2(t)dt minimum]. In fact, the direct numerical solution of the linear system (eq. 9) is generally unstable; thus it seems better to use the formalism described by Courtillot et al. (1978) and summarized in Appendix I. Following these authors, the data space ~ 12 and the functional model space are both orthonormalized. Eq. 8 is then replaced by: 12

¢(t) = ~. yiVi(t)

(10)

t=1

where V+(t) is the ith eigenfunction of the operator GG* given in Appendix I and corresponding to a positive eigenvalue X~ (with hl/> ~,~ I> . . . ~ ~2 > 0). The computation of Yi and Vi(t) is straightforward and the advantage of this procedure is the following: It is possible to show that, when the data are affected by independent errors with the sarhe variance, a s, the variance of Yl is given by a2/Xp and becomes greater for smallest eigenvalues, X~. This observation leads to retaining truncated solutions of the form:

247 Q
~Q(t) =

Z

yi Vi(t)

(11)

i=l

which result in solutions with smaller variance but where some misfit with the data is accepted (the higher Q is, the smaller this misfit is).

Application to actual data The previous technique was applied to the data of Fig.4. A mean conductivity of 3.2 W m -1 ° C -1 and a mean diffusivity of 1.6 • 106 m 2 s -~ are chosen (K = K / p c ) . The ten first eigenfunctions, Vi(t), resulting from the orthonormalizing process on functions Gi (t) are shown on Fig.5. For each eigenfunction, the corresponding eigenvector o f the data space ~ 2 is also shown; it can be interpreted as the result on the heat flow of a surface temperature disturbance ¢(t} similar to Vi (t). The three first eigenfunctions -- corresponding to higher eigenvalues -obviously express the influence of recent events; in the data space they involve a negligible contribution of the deepest data. Higher-order eigenfunctions are characterized, on a logarithmic plot, by a greater number of oscillations and also a higher influence of older events and deepest data. The last two eigenvalues are obtained as numerically null. Using the graphs of Fig.5, it is possible to discuss the t y p e of information which can be expected from the data. Since all eigenfunctions vanish for t < 5 yr. and t > 8000 yr., no information about these time intervals can be obtained; in particular the influence of the last glacial epoch appears to be negligible. Moreover, the resolving power of the m e t h o d can be evaluated from the "half-period" of highest-order eigenfunctions: for temperature events taking place around 30, 300 and 3000 yr. B.P., the resolving power is o f the order of 2 0 , 2 0 0 and 2000 yr., respectively. The next step is the computation of the coefficients Yi of the projected solution. As discussed above, the data used in this computation are corrected for various effects and are written as: ~/i = qi + topoi + Wi + Gli -- qo

(12)

where qi is the observed heat flow at level zi; topoi is the topographic correction at the same level; Wi is the correction due to heat production (A0zi); Gli is the correction arising from thermal events located outside the studied epoch (e.g. last glaciation); and q0 is the actual heat flow at ground surface. The two corrections, topoi and Wi, are given above (Fig.4d). For the correction Gli, only the effect of the last glacial epoch is considered; since the area has not been ice-covered during the glaciation, this effect is approximated b y a --10°C step f u n c t i o n (temperature difference with respect to the present temperature) between 65,000 yr. and 10,500 yr. (Lamb, 1977; Woillard, 1979). As shown on Fig.4d, this correction Gli is positive and

248 10 I

102

10 3

I

10 4

0

500

IY

I

I

mWm

j • 2 1 = 1444

1000 -Im

-2

I

1 0 -1

T

0

'~ ;4 ,-

:I -,....

~ .~--I:I,

I

Z

-1 E l 1

~ 2 2 = 98

1 0

0

'i

I . . .- i -. . .

-J-

-1 -1 1

j~23.15

0 d~.24 - 3 . 3 2

I.i

e

-1

1 0 -1

-1 I j~25=,

1 0 -1

/-7~

845

/-i~_/J~26 - , 188

,i i

-I

~. . . .

11 o

.

L -,•. . . .

e,

•

d~27= . 0293

d~28 - . 0 0 3 1 8

L

~

,

,

. . ,-i.,..

-1

'V o

1 0 -1

I I•,•,

. i l< L

i,

OL-I-.,

1 0 -1

,

-1

11

1 0 -1

le

: ";" ;

-1L

I .1

•:

"

'P

'~ .1° I

°

;

t

d~29 - - 00025 1 0 -1

1[

~'"

IL. . . . . .

d~210- • 00001

1

'r 0

0 -1

- I - i,,= i III

_IL . . . . .

Model

space

" '"

•

e

"~'

I i~ J i II I

i.

I

, ;,

Data space

Fig. 5. Left part: plot of the model space eigenfunctions Vi(t ) (arbitrary units) corresponding to decreasing eigenfunctions k~. Right part: plot of the corresponding data space eigenvectors u i as a function of depth.

249

decreases slightly from 12 to 8 mW m -2 in the depth range 100 to 900 m. The proper choice of q0 is rather delicate. From other geothermal measurements in the western part of the Massif Central, the regional heat flow value would a m o u n t to 100 to 120 mW m -2 (Lucazeau et al., 1981). However, its actual value is not known and several assumptions have to be made. For example, assuming a value of 106 mW m -2 for q, the corresponding solution is shown in Fig.6. In this figure, various truncated solutions, ~Q, using more and more eigenfunctions are plotted with their confidence domain: Q

660(t) =

V~

Y~ G~ ~ ( t ) 1

and the resulting 7i values are compared with the actual data. For Q = 1, 2, 3 these truncated solutions have small confidence domains; but t h e y cannot explain the observed 7 i values. The fit of the data is improved for Q > 4. For Q > 6, this improvement also results in solutions with very large variance and large-amplitude oscillations. The effect of our assumption for the q0 value is illustrated in Fig.7 where the truncated solution for Q = 6 is obtained for three different values of q0. In spite of some difference in amplitude, the three different values of solutions clearly show c o m m o n features: three maxima around 30 yr., 150 yr. and 1500 yr. before origin. This stability leads us to think that these features may be significant. However, these solutions are characterized by a fundamental ambiguity: it is possible to add any function which is orthogonal to the subspace of functions, Gi (t). On the other hand, the amplitude of temperature oscillations obtained by this technique are larger than +2°C, whereas during the Holocene actual climatic disturbance would not likely exceed 1 or 2 ° C. For this reason we now develop another inversion technique which can take into account more realistic constraints. INVERSION

USING BOUNDS

Principles of the method We now discretize the considered time range into M intervals and assume that ~(t) can be approached by a stepwise constant function: ¢(t) = q~ for t; < t < t~. + 1. Then the inverse problem is to obtain the u n k n o w n values of @j such that: M

7 i = ~, j=l

with:

giJCY

(13)

250 m W m -2

"C

5-

0

0

,/"

-2

4-+~ ~: + "4~I' ~r d-: ',(+ I Z

- f

~2 2 0

5-

~" ~r~--..~.~.:..i~._

++44-

0

-2

V , ,~P~,

:

,4-~4- 4 -i-''l i+ I

-5

(/L) 3

/r-.%+,

0 -2

0

,~ ~ I ~ . ,

-s/¢"7t .....

,4- ~

'

~*'"

0

d-:' ' ]

jp

0 -2

.

,

^

P.~..--'M..

,~..~

.>-~~

o

10 0 .10

I

I

I

I

1

10

102

103

Model

IY

Q

104

}

0

I -

500

pm

1000

Data

Fig.6. L e f t part: t r u n c a t e d s o l u t i o n s 0Q(t) = ~ YiVi(t) as a f u n c t i o n o f t i m e in years •=1 b e f o r e 1980 (log scale), a n d c o n f i d e n c e d o m a i n s (1 s.d.) for Q = 1 t o 8 (q0 = 106 mW m-2). R i g h t p a r t : c o m p a r i s o n o f t h e data (crosses) t o t h e c o m p u t e d value (line) as a function of depth.

251 mWm-2

Q o = 103 m W m - 2

"C

it

i "

4 ~'<~''' "

-~

I

.~>

'-~

; 7

........ _

Q o = 106 m W m

-2

! 2 *

0

I

li

-2

Oo=

109 m W m - 2

2 O-

0 -2

--k

- 5-

q 10

I" I

I 102

t 103

-H 104 Y

F 0

I SO0

--I m 1000

Fig.7. Left part: truncated solutions and confidence domains assuming q0 = 103, 106, 109 mW m -~. Right part: data and computed values. tj+l

gu =

~

G~(t)dt : K [ ( 1 / ~ )

exp (--z~/4Kt i + ,)

ti -- ( 1 / x / - ~ s ) exp (--z~/4K tj)]

(14)

Using this formulation and relevant numerical algorithms, it is possible to impose further constraints on the solution in order to obtain more realistic results: for example, upper and lower bounds can be settled for each u n k n o w n ~s. Besides, the surface heat flow can be taken as an additional unknown. The effect of the glaciation is simply handled by an unknown corresponding to the relevant time interval (65,000 yr. to 10,500 yr. before origin). Moreover, the errors in the data can be taken into account. Mathematically, this problem consists of an investigation of possible solutions of the following inequality system relating the unknowns (4, . . . . q~M, q0):

252

qi+topo~+Aoz~--Sqi<.qo+ ~ g ~ i . ~ j < q ~ + t o p o i + A o z ~ + S q i

(15)

J

q i being the error on data q~ evaluated above. The unknowns are b o u n d e d by: ¢i

rain

< ¢i < ~i max

(16)

where quantities labelled b y min or max indicate a priori bounds which have to be chosen. The following set o f bounds is proposed: (1) for temperature disturbances during the Holocene, bounds of +2°C seem adequate since Holocene temperature variations are expected to be quite small (Lamb, 1977); and (2) for the glaciation, bounds of --12°C and --6°C seem to cover the extreme range of accepted values (Emiliani, 1955; Washburn, 1980). In general the above problem (eqs. 15 and 16) has no unique solution b u t a set of solutions ("feasible solutions") (Dantzig, 1963). Some of them can be obtained b y using the algorithms of linear programming (Dantzig, 1963; Cuer and Bayer, 1980).

Results of the linear programming analysis A particularly interesting solution is the one which presents the smallest amplitude of variation during the Holocene. This "ideal solution" (Parker, 1975; Sabatier, 1977b) is defined as the one which gives the minimum of

Sup

t i e H ol o o m e

{¢i[

and produces a lower b o u n d for the amplitude of temperature variation during this period. For our data this lower b o u n d is 1.09°C; therefore, in order to account for the data, at least some values of I¢i[ must reach and even exceed 1.09 ° C. The corresponding solution is plotted in Fig.8a. It is negative (--1.09°C) for 100 < t < 900 yr. B.P. (1080--1880 A.D.) and positive for 20 < t < 80 yr. B.P. and 1000 < t < 10500 yr. B.P. (between 1900 and 1960 A.D. and between the postglacial warming and 980 A.D.). Since the temperature variation of this last solution has a quite realistic amplitude, it is interesting to obtain a slightly different related solution using the following procedure: "ideal solutions" are sought for all the Holocene period except for some time intervals when it is not possible to decrease the I¢il values further (mathematically the corresponding free unknowns are the ones which correspond to the highest values of "marginal costs"; Cuer and Bayer, 1980). In this way we obtain new solutions with an upper b o u n d (Sup I¢il) lower than the previous one, except for those time periods where ¢i is assigned as free (but b o u n d e d b y +-2°C). Proceeding in this w a y we obtain solutions where more and more ¢t's are assigned as free; for the other ~i's, the obtained upper b o u n d significantly decreases and finally reaches a very small value. In Fig.8b we obtain a solution where ¢i reaches + 2 ° C

253 mWm

--~L - - - J I

o

'~L

Ii . . . . . . .

~'

free

b

2~

~---

1°°I

]

i .

o

-2

1°°1

"C

.

.

.

.

.

.

.

.

.

.

r___l

.

L__.

I Ii

_2~

1°°t

free C


'*,

i L__

-2

_J

- >

I

TT'

90 1 80

i t

'

free t 1

10

--

t 10 2

t

i

10 3

10 4

Y

I

t

0

500

f

rfl

1000

Fig.8. Left part: "ideal solutions" obtained by minimization of Sup I¢tL for the whole Holocene (a); for the Holocene except from 2200 to 5000 yr. ago Ib); for the Holocene except from 2200 to 5000 yr. ago and from 200 to 450 yr. ago (c). Right part: data and computed values.

between 1800 and 5000 yr. B.P., the other ~i being bounded by +0.6°C. In Fig.8c we interpret the data with t w o time steps: --2°C between 200 and 450 yr. B.P. (1530 to 1780 A.D.) and +2°C between 1800 and 5000 yr. B.P. (3020 B.C. to 180 A.D.), the other ~ values being very near to 0. These two time periods correspond to the unknown quantities which are more significantly constrained by the data. In these solutions the data are taken into account b y inequality constraints. A larger weight can be conceded to the data b y computing the solution which gives the best fit for the LI or the L= norms: in the case of L1 we minimize the sum of absolute differences between the data and the computed values and for the L~ norm, their greatest difference. These t w o solutions are shown in Fig.9a, b. The temperature variation of Fig.9a ( L ) has three positive periods, which agrees quite well with the maximum of the projected solution (Fig.6), whereas Fig.9b (L1) has only two; in b o t h cases the amplitude of variation is more realistic {bounds of +1.3°C are imposed). As a conclusion of this investigation, the second technique does not lead to a single paleotemperature variation b u t it allows us to point out the main features of the temperature variations which are comparable with the data. These features are the lower b o u n d on the temperature variations and the presence of several significant positive and negative events.

254

1oo m

"C a

il

L

J

.

h

,

,

W

m

°2

90

•

83

2 :-~'- .....

t_.J

0

i

~"i-i

' I L.,',

!

10090

"~ .............

'

~

T

!'

,

80

-2 Y

I

1

10

102

103

104

I 0

~ 500

m 1000

Fig.9. a. Solution giving the best fit to the data for the Loonorm. b. S o l u t i o n giving the best fit to the data for the L~ norm. Bounds o f 1.3°C are imposed. CONCLUSION

From our observations in this 1000 m borehole drilled in granite, it is quite obvious that the measured heat flow is significantly modified by past disturbances of the ground temperature. Conversely, the observed modification may be used as a distorted record of past ground temperatures during the Holocene period. However, in order to carry out a quantitative interpretation in terms of paleotemperature, it is necessary to discuss in detail the role of various other disturbing phenomena and also to develop a relevant mathematical analysis. From a methodological point of view, the linear inverse problem of paleotemperature reconstructions is largely underdetermined. The classical inverse theory provides projected solutions which can explain the observed data but leads to an unrealistic amplitude of paleotemperature variations during the Holocene period (several °C). Therefore, we developed a second inversion method in order to obtain more realistic solutions; in this method, based on linear programming techniques, a priori bounds can be imposed on the unknown temperature. Using these two methods, the following results are obtained: (a) The data can be explained by disturbance of the ground temperature during the Holocene reaching an amplitude of at least +1.09°C. This amplitude is reasonable when compared to various estimates for the same period (La Marche, 1974; Markgraf, 1974; Matthews, 1976; MSrner, 1977; etc.). (b) For all the solutions obtained, there exist two time periods, around 1600 A.D. and 1500 B.C., when the ground temperature is systematically and respectively below and above the present temperature. The first event is well known: it corresponds to the so-called "Little Ice Age" which affected Western Europe during the late medieval epoch (Lamb,

255

1965). From our solutions its amplitude is difficult to specify: it could be as low as --0.6°C (Fig.8b) but then the event would have lasted 900 yr. Other models with shorter duration of this cold episode (250 yr.) can be obtained with the assumption that lower temperatures (--2°C) are admitted (Fig.8c). The second event is a warm episode which ended around 900 yr. ago. Due to the resolving power of our data it is impossible to specify accurate variations during this early Holocene period. However, the time period which is more significantly positive runs from about 3000 B.C. to the birth of Christ. This period covers the Subboreal and Subatlantic periods which are not characterized by significant warming (Lamb, 1977). For the more recent period we obtain a greater scatter in the solutions. Most solutions indicate the presence of a relative m a x i m u m around 20--60 yr. ago (Figs.7 and 9b); this event seems to match with a well-documented global warming reaching its m a x i m u m around 1940 A.D. (Mitchell, 1961). The same solutions indicate the presence of another maximum around 1850 A.D., which could correspond to the end of the Little Ice Age. However, with other solutions (e.g. Fig.9a) this fluctuation disappears. The general features described above give interesting results. However, for a better comparison of our solutions with other indirect data about paleoclimates, the following points must be taken into account: (1) The relationship between the mean air temperature and the mean ground temperature is not an obvious one; in particular it could depend upon rain and snow precipitation. (2) Whereas some temperature events are well recognized all over the earth or at least on a large geographical scale, there could also exist a large variability for some climatic events with a short duration. In our case, it seems difficult to explain the observed heat flow variation as a function of depth by other phenomena than paleotemperature variations similar to the ones presented. However, the paleotemperature reconstruction is based on a single borehole observation and more data are needed to ascertain the results. Anyhow, it is clear that detailed geothermal measurements in favourable sites and relevant mathematical analysis can provide valuable direct information on the ground temperature. APPENDIX I

The scalar product of two quadratically integrable functions f and g on the inverval (0, + ~ ) is defined by the integral: oo

(f, g) = ~ f(t)g(t)dt

(A1)

0

Let us introduce the vectorial operator G which associates with function

re(t) the vector p (which belongs to ~N) defined by:

256

Pt = j G l ( t ) m ( t ) d t = ( G i , m ) 0

fori=l,...N

(A2)

and G*, adjoint of the operator G, which associates with any vector q = (q 1, q2, • •., qN) (which belongs to . ~ ) the function n ( t ) defined by: N

n ( t ) = G*q = Y. q i G i ( t )

(A3)

i=1

The GG* operator is a classical matrix operator in ~ g characterised by a symmetrical matrix. The elements are given by: g~j = f V ~ ( t ) V j ( t ) d t = [K2K/4e] [384(Z~ Z2~~/Z2Z'2~,41, ~ i, 0

- 4 S ( Z ~ i + Z~J)-2]

(A4)

If the G~(t)'s are linearly independent functions, the g~j matrix is a definite positive one. Let k~ and ui (i = 1, . . . , N) be the eigenvalues and the eigenvectors (belonging to .~N) of the gij matrix (with ~ / > ~ / > . . . / > k~v > 0). Let us define the Vi (t) functions by: G*ui = E u~Gi(t) = k i V i ( t ) i

(A5)

it follows that: G Vi = 1/kiGG* ui = X~/ki ui = k~ui

(A6)

thus: (V~, Vi) = (1/k/k1) • u~uJ(Gk, Gl) = 5ij

(A7)

k,!

The Vi functions (i = 1, . . . , N) generate an orthonormal base of the subspace which is defined by the Gi functions. One can demonstrate that Vi functions are also eigenfunctions of the G*G operator (with corresponding k~ eigenvalues). This problem can be projected in data space (on ui vectors belonging to ~N ) and in functions space (on Vi(t ) functions); then it allows us to separate the effect of each eigenfunction. From eq. A6 it follows that to an assumed model Vi(t) corresponds a set of data which are the components of ktu i. Then one seeks the solutions of thiainverse problem o f the form: N

~(t) = Y

y~Yi(t)

i=l

It follows that:

(AS)

257 N

(A9)

7i = (Gs, ¢ ) = ~. Y i ( G j , V i ) = ~. YiXiu~ i

i=1

Or if we consider the matricial formulation: 7 = UAy

(A10)

where U is an N × N matrix whose columns are the eigenvectors ui, A the diagonal matrix of eigenvalues Xi. Since U is an orthogonal matrix (that is to say satisfying U T u = U U T = IN) the solution is obviously: y = A-1 uT'/

(All)

From eqs. A8 and A l l one can calculate the solution. However, it is often interesting to truncate the solution to the form: Q ¢*(t) = ~

y~ V i ( t )

withQ
(A12)

i=1

Indeed, suppose that the data errors di (i = 1, . .., N) are independent random variables of variance 02, one can demonstrate easily that the y vector is associated with a covariance matrix: By = 02A -I

(A13)

Thus the variance of each yl is independent and increases with i like 1/X~ i

= o lX ).

The truncated solution does not satisfy the data set: '/* = G ¢ *

(A14)

We have: N

71" - - 7j = ( G j , ¢* - - 49) =

~. i=Q+I

N

Y i ( G s , Vi) =

~

XiYiu~

(A15)

i=Q+l

The quadratic deviation between data and fitted model is given by: N

117"--711 = Y~ X2.2 i~Vi

(AI6)

Q+I

which decreases when Q increases. This property joined with the estimation of the truncated solution variance corresponds to the well-known "trade-off" between the solution variance and the misfit to the data. ACKNOWLEDGEMENTS The borehole near Auriat was executed by the French Commissariat l'Energie Atomique as a contribution to a research program on the feasibility

258 o f nuclear waste burial in d e e p geologic media. This w o r k is also part o f a research a n d d e v e l o p m e n t p r o g r a m o f t h e C o m m i s s i o n o f E u r o p e a n Communities. We are grateful to t h e I n s t i t u t de P r o t e c t i o n et de Suret~ Nucl~aire p o u r les D~chets R a d i o a c t i f s f o r t h e p e r m i s s i o n f o r access t o t h e b o r e h o l e and for its help during t h e m e a s u r e m e n t s . We t h a n k MM. Sousselier, Barbreau, Derlich and Peruss for their valuable c o l l a b o r a t i o n . REFERENCES Backus, G. E. and Gilbert, F., 1967. Numerical applications of a formalism for geophysical inverse problems. Geophys. J. R. Astron. Soc., 13: 247--276. Beck, A. E. and Judge, A. S., 1969. Analysis of heat flow data: detailed observation in a single borehole. Geophys. J. R. Astron. Soc., 18: 145--158. Birch, F., 1948. The effects of Pleistocene climatic variations upon geothermal gradients. Am. J. Sci., 246: 729--760. Birch, F., Roy, R. F. and Decker, E. R., 1968. Heat flow and thermal history in New England and New York. In: E. An Zen (Editor), Studies of Appalachian Geology. Interscience, New York, N.Y., pp. 437--451. Bullard, E. C., 1940. The disturbance on the temperature gradient in the earth crust by inequalities of height. Mon. Not. Astron. Soc., Geophys. Suppl., 4: 300--362. Carslaw, H. S. and Jaeger, J. G., 1959. Conduction of Heat in Solids. Oxford University Press, London, 2nd ed. Cermak, V., 1971. Underground temperature and inferred climatic temperature of the past millenium. Palaeogeogr., Palaeoclimatol., Palaeoecol., 10: 1--19. Chenevoy, M., 1958. Contribution ~ l'4tude des schistes cristallins de la partie Nord Ouest du Massif Central franqais. M6m. Explicatif Serv. Carte Geol. Fr. Courtillot, V., Ducruix, J. and Le Mouel, J. L., 1978. Inverse methods applied to continuation problems in geophysics. Lecture Not. Phys., 85: 48--82. Cuer, M. and Bayer, R., 1980. Fortran routines for linear inverse problems. Geophysics, 45: 1706--1719. Dantzig, G. B., 1963. Linear Programming and Extensions. Princeton University Press, Princeton, NJ. Derruau, M., 1971. Sur la morphologie du Massif Central franqais. Symp. J. Jung, Plein Air Serv. Ed., Clermont-Ferrand, 33--44. Duthou, J. L., 1977. Chronologie Rb-Sr et g6ochimie des granitoi'des d'un segment de la chafne varisque, relation avec le m6tamorphisme: le Nord Limousin (Massif Central franqais). Thesis, Clermont-Ferrand, Ann. Fac. Sci. Clermont, 6 3 : 2 9 4 pp. Emiliani, C., 1955. Pleistocene temperatures. J. Geol., 63: 538--573. Hamilton, M., 1964. Statistics in Physical Science. Ronal Press Co., New York, N.Y. Jordan, J. H., 1973. Estimation of the Radial Variation of Seismic Velocities and Density in the Earth. Thesis, California Institute of Technology, Pasadena, CA, 199 pp. Kast, Y. and Gallice, J., 1980. Mise au point et construction d'un conductivim~tre rapide. In: Abstr. 26th Int. Geol. Congr., Paris, p. 1122. La Marche, V. C., 1974. Paleoclimatic inferences from long tree ring records. Science, 183: 1043--1048. Lamb, H. H., 1965. The early medieval warm epoch and its sequel. Palaeogeogr., Palaeoclimatol., Palaeoecol., 1: 13--67. Lamb, H. H., 1977. Climate: Present, Past and Future, 2. Climatic History and the Future. Methuen and Co., London, 833 pp. Lucazeau, F., Vasseur, G., Kast, Y. and Jolivet, J., 1981. Donn~es de flux de chaleur darts le Massif Central francais. Ann. G$ophys., 37: 481--491. Markgraf, V., 1974. Paleoclimatic evidence derived from timberline fluctuations. Coll. Internes CNRS, 219: 67--76.

259 Matthews, J.A., 1976. Little Ice Age paleotemperatures from high altitude tree growth in south Norway. Nature, 264: 243--245. Mitchell, J. M., 1961. Recent secular change of global temperatures. Ann. N.Y. Acad. Sci., 95(I): 235--250. MSrner, N. A., 1977. Climatic framework of the end of the Pleistocene and the Holocene: paleoclimatic variation during the last 35,000 years. Geogr. Phys. Quat., 31: 23--35. Parker, R. L., 1975. The theory of ideal bodies for gravity interpretation. Geophys. J. R. Astron. Soc., 42: 315--334. Sabatier, P. C., 1977a. On geophysical inverse problems and constraints. J. Geophys. Res., 43: 115--137. Sabatier, P. C., 1977b. Positivity constraints in linear inverse problem. Geophys. J. R. Astron. Soc., 48: 415--422. Washburn, A. L., 1980. Permafrost features as evidence of climatic change. Earth-Sci. Rev., 15: 327--402. Woillard, G., 1979. The last interglacial~lacial cycle at Grande Pile in Northeastern France. Bull. Soc. Belge G4ol., 88: 51--69.