Whole mantle convection and the thermal evolution of the earth

Physics of the Earth and Planetary Interiors, 29 (1982) 281—304

Elsevier Scientific Publishing Company, Amsterdam

—

281

Printed in The Netherlands

Whole mantle convection and the thermal evolution of the Earth W.R. Peltier 1 and G.T. Jarvis 2 ‘Department of Physics, University of Toronto, Toronto, Ont. M5S JA 7 (Canada) 2 Department of Geology, University of Toronto, Toronto, Ont. M5S IA 7 (Canada)

(Received and accepted for publication March 22, 1982)

Peltier, W.R. and Jarvis, G.T., 1982. Whole mantle convection and the thermal evolution of the Earth. Phys. Earth Planet. Inter., 29: 28 1—304. Three recent results in the fields of geophysics and geochemistry have combined to rejuvenate the study of the Earth’s thermal history. The first, and perhaps most important, of these is the new understanding of the planetary accretion process which has been realized as a consequence of work by V.S. Safronov and others. This work shows that the terrestrial planets accrete on a time scale of -~ l0~y by the accumulation of planetesimal sized objects and that they probably form sufficiently hot to induce ‘immediate’ differentiation of iron from silicates. The second important recent contribution has been from geochemistry and concerns the net complement of radioactivity of the Earth as a whole. Although the net complement of radioactive elements has usually been inferred on the basis of the assumption of equilibrium with the surface heat flow, a priori analysis based upon elemental condensation temperatures in the primitive solar nebula yields a model in which the net radioactivity is at least a factor of two less than would be required to explain the mean surface heat flow. There are good reasons to believe that this model is valid for the Earth. The third result has been the successful development of a technique for incorporating the influence of convective heat transfer in a thermal evolution model which does not require explicit solution of the dynamical equations. This paper begins with an analysis of the mantle convective circulation which is required to understand plate tectonics and continental drift, an analysis which presents the arguments for and against the crucial notion that the main circulation has a depth scale equal to the mantle thickness. Using the idea of ‘parametenzed’ convection, the radial heat transfer by the mantle convective circulation is included in a thermal history model of the Earth. This model turns out to contain within it an interesting explanation of the fact that no crustal rocks with ages in excess of 3.8 Ga have ever been found, an observation that has heretofore proven enigmatic.

I. Introduction: possible styles of convection in the Earth’s mantle The recent literature on thermal convection in the Earth’s mantle has been characterized by the emergence of two distinctly different models and of two groups of respective model devotees. Although both of these models are certainly oversimplifications of reality they are nevertheless sufficiently different from one another that the question of which is closer to the truth is probably a meaningful one. These two models are illustrated schematically in Fig. 1 and will be referred to as 0031-9201/82/0000—0000/$02.75

the ‘whole-mantle’ model (Fig. la) and the ‘layered’ model (Fig. lb). In the former, the entire mantle is mechanically well-mixed throughout its volume and the mean temperature and density fields vary adiabatically with depth except in thin (of the order of 102 km thick) boundary layers adjacent to the Earth’s surface and to the core. At the Earth’s surface the oceanic lithosphere is perfectly coupled to the motions in the mantle beneath and is coincident with the thermal boundary layer. In this model the continents constitute a combined chemical and thermal boundary layer which may be somewhat thicker than the oldest oceanic litho-

© 1982 Elsevier Scientific Publishing Company

282 TEMPERATURE (K)

~‘°:~j~!L~

3c~o

o~ç

(a)

000

lOOC

~

WHOLE MANTLE 2000

LOWER

E

I—

16

I

CONVECTION

_J

MANTL

18

20

%DENSITY

SOUDUS

I

.~

X2000

CHEMICAL AND

BL

-

CORE

22

4

6

8

10

LOG

3) 10 (VISCOSITY [Pa s~)

DENSITY (Mg/rn Mi

J~_

UPPER MANTLE

~

////////

00020003000

/~/////~ -

000

18

20

/

/

/

~LAYEREDCONVECTION “.-..,

“ “‘i,

.~

1000

%DENSITY

~

i::J:TT 16

TEMPERATURE(K)

-

22

*

OLIDUS

~T:~CMB £

LOG

4

6

8

0

3)

10(VISCOSITY DENSITY (Mg/rn Fig. I. Schematic diagrams[pa of s]) the whole-mantle andCOR layered convection models with associated viscosity, density and temperature profiles. Note the question mark adjacent to the viscosity profile for the layered convection model. We expect the large viscosity drop only subject to the assumptions discussed in the text.

sphere, particularly beneath old continental cratons. This picture differs from that for the layered model in several fundamental ways as illustrated in Fig. lb. In the layered model, mass transfer across the seismic discontinuity at 670 km depth does not occur. Consequently, although the upper and lower mantles are individually well mixed they may have distinct chemical compositions. Because heat can be transferred across the boundary at 670 km depth only by diffusion, the layered model has a thermal boundary layer there through which the temperature must increase by at least 500 K (as shown in Fig. 1) and more likely by 1000 K if the upper mantle is depleted of radioactivity, as recent geochemical data have sometimes been taken to suggest. Although the lower mantle of the layered model is also made adiabatic by its convective circulation this adiabat must everywhere be at a

temperature which is from 500 to 1000 K higher than that of the whole mantle model. In consequence of this, for a fixed known temperature at the core—mantle boundary (CMB) the layered model must have a weaker thermal boundary layer at this interface which implies a smaller heat flux from the core into the mantle. In the upper mantle of the layered model there are also differences between it and the whole mantle alternative. Foremost among these, according to most advocates of this picture, is the notion that there exists a separate scale of convective circulation which is uncorrelated with the motion of the plates and which may be responsible for the radial transport of a significant fraction of the heat flowing across the upper mantle. Which of these two models is closer to the truth will, of course, have rather profound consequences

283

insofar as the thermal history of the Earth is concerned. Most important is the fact that for a given temperature drop across the mantle, and a given complement of internal radioactivity, the layered model will transport heat across the mantie significantly more slowly than will the wholemantle model. This can be understood most simpiy using the empirical relationship between the heat flow effected by convection and the Rayleigh number of the circulation. For a plane layer with constant physical properties and free stress, heated below, boundary conditions the heat flow q is given by: q

—

(R/R~) “3K L~T/L

where K is the thermal conductivity of the layer, L~Tthe temperature drop across it, and L the layer depth. The Rayleigh number R = gcsL~TL3/scv, in which g is the gravitational acceleration, d the coefficient of thermal expansion, and K and v the thermal and viscous diffusivities. R~is the critical Rayleigh number and is inevitably of the order of io~.On the basis of this empirical relation between q and R, which may in fact be derived from boundary layer theory, we may examine the effect upon q due to layering of the convective circulation. It is straightforward to show that if q 1 is the heat flow for a whole mantle model heated from below and q2 is the heat flow for a model which consists of two layers of equal thickness then q1/q2 (1/16)1/3. In fact, the heat flow for a layered model is reduced approximately by the ratio of the thickness of the thinnest layer to the mantle depth. We would therefore expect on this basis that the rate of heat transfer associated with the layered model shown in Fig. lb would be about 5 times smaller than that associated with the whole mantle model illustrated in Fig. Ia. The consequences for the thermal history of which of these models is correct are therefore extremely important. Are there any geophysical data which might allow us to distinguish unambiguously between these two models? In Table I we have listed the major observations which have been invoked by various authors to constrain models of mantle convection. This list is a slightly modified and extended version of that in Peltier (1981b). Of this —

TABLE I Constraints on models of the mantle convective circulation

_____________________________________________ (I) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) (14) (15)

Mean horizontal plate scale 1 4000 km 1 Mean surface plate speed heatuflow ~4 q cm~80 y~mW m2 Mean plate thickness 8 100 km Mean mantle viscosity v—1022—5)< 1022 poise Seismicity ceases below the 670 km seismic discontinuity Heat flow and topography a ~ of the ocean floor Volcanism observed in plate interiors—hot spots poise to Viscosity increases only modestly from 1022 about 3 X 1022 poise through the transition region a Boundary layer nature of the D” region above the core— mantle boundary Deep seismic focal mechanism compressive b Bathymetry deviates from (age) l/2 for ocean floor age~ 7.l0~y Seismic low velocity zone and high attenuation zone exist only beneath young oceanic lithosphere Rb—Sr and Nd—Sm constraints on mantle mixing and crustal growth b Earth sufficient radioactivity 50% ofcontains the observed surface heat flow to account for ~

(16) High pressure chemistry of the olivine binary suggests that both the 120 km and the 670 km seismic discontinuities are phase boundaries and not chemical discontinuities a a b

Observations used in support of the whole-mantle model. Observations used in support of the layered model.

list of 16 entries three (6, 11 and 14) have been invoked by several authors as providing strong support for the layered model; entries 1, 9 and 16 have been taken to imply the validity of the whole-mantle model. In fact, it was only with the discovery of point 9 (Peltier, 1974, 1981 a; Cathles, 1975; Peltier and Andrews, 1976) that the wholemantle model became tenable at all (Peltier, 1972, 1976). As the ideas of mantle convection have developed over the past decade, leading to the impasse illustrated in Fig. 1, discussion has usually been focussed upon the layered model or some simple variant of it. The reason for this had to do with early geophysical understanding of the variation of mantle viscosity as a function of depth. Modern work on this problem essentially began with McConnell’s (1968) analysis of Sauramo’s (1958) shoreline data from Fennoscandia. McConnell’s results included an estimate of the thickness of the

284

lithosphere in the Baltic Shield (— 120 km), and a measurement of the viscosity of the upper mantle which he found to be 1021 Pa s~,in agreement with the original work by Haskell (1937). In order to retrieve information on the viscosity of the mantle beneath the transition region, McConnell was obliged to invoke information beyond that contained in the shoreline data themselves. What he did was to assume (following Munk and MacDonald, 1960) that the Earth possessed a real nonhydrostatic equatorial bulge, i.e. that the polar flattening was in excess of that which would be in equilibrium with the current rotation rate. He further assumed that the excess bulge was produced by glaciation. It then followed from the existence of the bulge that the relaxation time of the second-degree harmonic must be extremely long, On this basis McKenzie (1966, 1967, 1968) and McConnell (1968) were led to construct models in which the viscosity beneath the transition zone was in excess of 1025 Pa s~. This argument for high lower mantle viscosity was completely undermined by Goldreich and Toomre (1969) who pointed out that the non-hydrostatic equatorial bulge did not in fact exist! It had been inferred from a spherical harmonic expansion which was improperly biased towards the second-degree harmonic. In spite of the fact that McKenzie’s original argument in favour of extremely high viscosity in the lower mantle was shown to be incorrect, this has not led to a rejection of the notion itself. In the model of mantle convection advocated by McKenzie Ct al. (1974) and that discussed by Oxburgh and Turcotte (1978) it is assumed that the circulation is confined to the upper mantle alone. In McKenzie and Weiss (1975) it is argued that either the upper mantle convects alone or that the upper and lower mantles are filled with distinct circulations (the layered model). The argument given for the existence of a separate circulation in the upper mantle is again that there exists a large increase of viscosity across the phase transition at 670 km depth, although no observational evidence is given to support this assertion. A speculation appears here that the supposed large increase of viscosity is effected by an increase of the creep activation energy of 2 eV. Observation 6 of ‘-~

Table I is invoked in support of the argument that the upper mantle flow cannot penetrate beneath 670 km depth. In Peltier (1972) an analysis in, terms of linear stability theory was given of the magnitude of the viscosity increase across 670 km depth which would be required to focus a distinct circulation into the upper mantle. Here it was shown that even a contrast of three orders of magnitude would be inadequate to produce such confinement. Peltier (1976) argued on the basis of these results and the new postglacial rebound analyses (Peltier, 1974; Cathles, 1975; Peltier and Andrews, 1976) that: “if the radial variation of mantle viscosity is not ‘extreme’ then the radial mixing length for thermal convection in this system is liable to be on the order of the thickness of the mantle itself”. This argument for the whole-mantle model of convection was later echoed by Davies (1977), O’Connell (1977) and Schubert (1979) and has now gained a considerable following. That the observed viscosity profile of the mantle is known to be rather uniform has not completely won the argument for the whole-mantle model of convection because of the additional fact that a separation into distinct upper and lower mantle circulations could be effected by means other than through a marked variation of viscosity. It has often been argued (e.g. see D.L. Anderson, 1981 for a recent review) that the mantle is chemically stratified. If the seismic discontinuity at 670 km depth were a chemical boundary (as shown in Fig. 1 b) then the stabilizing effect of a change in mean atomic weight across this depth could easily prevent mass transport across the boundary (Richter and Johnson, 1974). This idea has apparently received very strong support from recent studies of mantle differentiation and crustal growth using the Sm—Nd isotopic system (Wasserburg and DePaulo, 1979; O’Nions et al., 1979). These data have been interpreted as requiring the existence of two distinct chemical reservoirs in the mantle and mass balance arguments have been used to argue that the two reservoirs should be associated with the upper and lower mantles. The idea of chemical stratification coupled with the supporting isotopic evidence has therefore served to prolong the life of the layered-mantle model.

285

The layered model is nevertheless vulnerable to criticism because it is supported only by notions which are in themselves extremely contentious. On the question of chemical stratification, for example it has now been demonstrated (Mao et al., 1979; Jeanloz, 1981) through direct static high-pressure experiments on the olivine binary with the diamond anvil apparatus, that the seismic discontinuity at 670 km depth is a phase boundary. The observed change of density across the phase boundary, furthermore, is just that required by the seismic data (Mao et al., 1979). It therefore seems that there is very little room for any substantial density change due to a change of mean atomic weight and therefore that the mantle has no appreciable chemical stratification. This model has long been advocated, for example, by Ringwood (1979). Anderson (1981) has therefore been forced to suggest that the 670 km discontinuity is not only a phase change but also a chemical boundary. This seems rather ad hoc. The only evidence that the phase change may not be sufficient to explain the observed characteristics of the 670 km discontinuity is that the discontinuity must be rather sharp to explain the good seismic reflections which are obtained from it, and somewhat sharper than it has been found to be in the laboratory measurements. As Jeanloz (1981) has pointed out, however, the laboratory data have all been taken at a fixed temperature of 1000°C whereas at 670 km depth the temperature is 2000°C. Since the boundary is expected to be much sharper at 2000°C than it is at 1000°Cthere is as yet no necessity for invoking chemical stratification to explain the good seismic reflections which are observed, In spite of the great weight which was initially attached to the constraint on the degree of mantle mixing provided by the Nd—Sm data, these data have also proven to be rather more ambiguous than either Wasserburg and DePaulo (1979) or O’Nions et al. (1979) were willing to allow. Estimates of the fraction of the mantle which has been differentiated to form the crust vary from 0.25 (DePaulo, 1981), 0.33 (Jacobsen and Wasserburg, 1979), or 0.50 (O’Nions et al., 1979) to 1.0 (Davies, 1981). All of these analyses are based upon the assumption that the initial chemical state of the Earth is reasonably well known and given essen-~

tially by the abundances, listed by O’Nions et al. (1979, Table I). The abundances of the radioactive elements in such models are essentially those derived for the bulk Earth by Smith (1977) on a different basis and are lower by a factor of 2 than the values which would be required to explain the observed mean surface heat flow. The geochemical mass balance arguments which lead to estimates of 0.5 for the fraction of the mantle which has been depleted are based upon the assumption that the volume of the correspondingly enriched crustal layer is known accurately. lithe concentration of LIL’s in the bulk Earth and the corresponding concentrations in the present-day enriched crust are known, then it is very simple to calculate the fraction of the mantle which must have been depleted to form the crust. This is the basis of the calculations which have been performed. The most conspicuous source of error in these calculations of which we are aware is the assumption made concerning the volume of the enriched crustal layer. There are very good reasons, it seems to us, for suspecting that the volume of the crustal layer may be at least double that which has been assumed by the geochemists. These reasons have been reviewed in a different context by Jordan (1981) and may be invoked to argue, using the same mass balance arguments, that the entire volume of the mantle has been depleted to the same extent by the process of crustal differentiation. An additional blow to the original interpretation of the isotopic data has recently been struck by Carlson et al. (1981) who have shown that the continental flood basalts of the Columbia river have variable Sm—Nd isotopic ratios which may be explained by mixing different amounts of enriched continental crust with a depleted parent material, presumably the same depleted parent which is responsible for the production of mid-ocean ridge basalt (MORB). The geochemical constraints upon mantle mixing are therefore rather weak, particularly since the MORB and OIB arrays have considerable overlap (Allegre, 1982). Of greatest importance, however, is the geochemical inference that the net complement of radioactivity in the planet is only about half that which would be required to explain the surface heat flow. We find the arguments for this primary conclusion rather convincing. ‘~-

286

Of the observational constraints upon the mantle convective circulation listed in Table I, that numbered 9 is probably the most unambiguously diagnostic of which of the two models shown in Fig. 1 is nearer the truth. The reason for this has to do with the main difference between these two models which is the presence of the thermal boundary layer at 670 km depth. We may calculate the variation of viscosity to be expected across this internal thermal boundary layer on the basis of a (perhaps non-Newtonian) flow law written in the form = AD(a/~s)” (1) where: p. is the elastic shear modulus; a is the deviatoric stress; E is the strain rate; n is the power law exponent (—3 for non-Newtonian and 1 for Newtonian flow); A is an empirical constant depending upon the absolute temperature T. The diffusion coefficient, D, depends upon the self-diffusion of the rate-limiting species and has the form D(P, T)

=

oxygen diffusion then the zero pressure activation energy E* is related to oxygen ion packing by E*(kcal mol~) (7) = (187 ± 16) (3.8 ±0.8) V (A03) where V~/is the volume per 0oxygen ion at zero pressure and 25°C.It is related to the zero pressure density Po by —

Vo = M/POPNA (8) where M is the molecular weight, NA is Avogadros number, and 4~is the stoichiometric coefficient of oxygen (~4). On the basis of the above analysis and under the assumption that the creep mechanism is indeed rate controlled by the diffusion of oxygen ions, 6E * can be estimated from the seismically observed change in density across the phase boundary (z~p= 0.27 g cm3 in model l066B of Gilbert and Dziewonski (1975) and i~sp= 0.39 g cm3 in the PREM model of Dziewonski and Anderson

D 0

e~/ItT

(2)

where G*(P, T) is the Gibbs free energy for activation G*=E*+PV*_TS*

(3)

in which E*, V”, and S~are activation energy, volume and entropy, respectively. For this flow law the effective viscosity is v = cs/2E = (p.”a’”/2AD 0) e~*~~T (4) Now if T, P and a are continuous across the phase boundary and n is insensitive to structural changes then the ratio of viscosities on the two sides of the phase boundary is (Sammis et al., 1977). 6~’°8~’~ (5) = x e~ Here v 1 is the viscosity above the transition, v2 the viscosity below, and x is a constant near to unity. Since V” is closely related to the atomic volume it must be smaller in the more dense perovskite + magnesiowustite phase beneath 670 km depth so that 8V* <0. An upper bound on the viscosity ratio is then

(1981)). Given 8E* estimated in this way from (7), which has a value of —6 kcal mol’ (Sammis et al., 1977), substitution irito (6) gives (~‘2/’I)whoIemantle~

10

(9)

so that the viscosity is expected to increase by less than a factor of 10 across the phase boundary if the temperature varies slowly and adiabatically 1)through this region. As shown in (—‘ 0.3 K km~ Fig. la, therefore, if the whole mantle is mechanically well-mixed and thermally adiabatic we except a modest increase of viscosity through the transition region. If there is a thermal boundary layer at 670 km above quite eq. depth, prediction however, changes as shown inradically. Fig. lb Using then the 4, a similar equation to (6), and a value of 200 kcal mol’ for G*, we see immediately that with the 500°Cincrease of temperature across the 670 km thermal boundary layer suggested by Jeanloz and Richter (1979), we have (~2/~I)1d mantle ~ l0~ (10)

(6)

The analysis then shows that the viscosity should decrease across the phase transition by about 4

Sammis et al. (1977) showed that if the creep mechanism were rate-limited by intrinsic volume

orders of magnitude. The two models of convection therefore predict strikingly different mantle

v2/v1


287

viscosity profiles and we might hope to determine which model is nearest the truth by using a suitable set of geophysical data to measure the viscosity depth profile accurately. The viscosity depth profiles associated with the whole mantle and layered models are shown in Fig. Ia and b, respectively, An independent method of estimating the viscosity profiles expected in the whole-mantle and layered convection models leads to a similar prediction to that obtained on the basis of the above analysis in terms of the activation parameters. Weertman (1970) and Weertman and Weertman (1975) have given some fundamental justification for an empirical rheological model of the form v(T)v0exp(gT

m

/T)

(11)

in which Tm(P) is the pressure-dependent melting temperature of mantle material, and g is a dimensionless constant which controls the rate at which mantle material weakens as the temperature approaches the melting temperature. For metals g 18, while for silicates g 30. In Fig. 1 the melting temperature profile (solidus) is shown on the plates containing the model geotherms for the whole mantle and layered model and is that given by Kennedy and Higgins (1972). Clearly on the basis of (11) the effect of the thermal boundary layer at 670 km depth is to drastically decrease the Separation in temperature between the geotherm and the solidus and thus to force a sharp decrease of viscosity at this boundary just as was predicted on the basis of the expected depth dependence of the activation parameters. Even with the minimum boundary layer temperature rise of 500 K, the geotherm of the layered model actually crosses the solidus! On the basis, of both of the above arguments, the layered convection model appears to be incompatible with a constant viscosity mantle. In the next section we present the data which require that the mantle viscosity be uniform and will take these as strong circumstantial evidence in favour of the whole mantle model. It should be noted that this argument is rather different, from that employed by Tozer (1972) who suggests that, for any rheological model, the circulation will eventually adjust so as to produce a geotherm such that the viscosity is everywhere equal to 1021 Pa s’. It must be

recognized, however, that our argument relies upon an assumption that the Sammis et a). (1977) empirical analysis of the increase of creep-activation energy across the 670 km phase transition is essentially correct and/or that the Kennedy and Higgins (1972) melting curve is an adequate approximation. If there were a substantial increase of the creep activation energy across 670 km depth (an increase of 45 kcal mol~ (— 2 eV) would be needed to compensate for an increase of temperature of 500°C)and a corresponding sharp increase in the melting temperature in the more dense phase, then we could see a constant viscosity mantle even if an internal boundary layer were present. Such a suggestion strikes one as being somewhat contnved, however, and until laboratory measurements of the critical parameters are forthcoming in support of it the idea is not really credible. Below (section 2) we describe the data which require that the viscosity of the mantle be essentially uniform, and recent results from high Rayleigh number simulations of whole-mantle convection (section 3). An attempt is made to address the question of the extent to which the mantle convective circulation is forced by heating from within. Section 4 concerns an application of these ideas to an analysis of the thermal history of the Earth based upon the idea of parameterized convection developed in Sharpe and Peltier (1978, 1979) and Peltier (1980).

2. Arguments for a mantle with weak viscosity stratification Inference of the variation of mantle viscosity with depth is based upon inversion of the observations of glacial isostatic adjustment in terms of a theoretical model of the adjustment process. Recent detailed reviews of this problem will be found in Peltier (1981, 1982) and references cited therein. We give only a rather cursory discussion here which will concentrate upon the observations and their interpretation rather than upon the theoretical details. The data consist of the following three distinct types of information: (1) relative sea-level data; (2) free-air gravity anomalies; (3) polar motion and acceleration of rotation data.

288

2.1. Relative sea-level data

are located near the centre of postglacial rebound in Hudson’s Bay, Boston is very near the margin of the ancient Laurentide ice sheet, Clinton (CT) is in the peripheral region of submergence, and Recife (Brazil) is in the far field of the main deglaciation event. In Fig. 3 we show a set of comparisons between the observed relative sea level histories at a sequence of six locations which were under the Laurentide ice sheet and the predictions of the theoretical model (Wu and Peltier, l982b). Inspection of these data shows that the model with highest lower mantle viscosity is strongly rejected and that the uniform model satisfies the largest fraction of the information. This picture is reinforced by Fig. 4, which compares relative displacement predictions of a circular disc load model of the ice sheet with three different sets of central Laurentide rsl data for Earth models with 1021 Pa

These observations consist of ‘4C controlled emergence and/or submergence curves. In consequence of the last deglaciation event of the current ice age, which began 18000y ago, the ocean basins were flooded by an amount of meltwater equivalent to a mean global rise of sea level of 80—100 m. Because of the removal of the ice load from continental regions, the unbalanced gravitational force induced uplift of the crust relative to the geoid in such locations. Beyond the ice margin relict beaches are found below sea level, both because of the deglaciation-induced collapse of the peripheral bulge in the solid Earth, and because of the meltwater added to the ocean basins. Figure 2 shows sample relative sea-level curves from four sites at successively greater distances from the ice sheet which covered Canada. The Ottawa Islands

I

200 Boston

-

(a)

-

I

I OtIawo Is.

40

(b)

GX-688

20

-

-

I

oW-99/ ~

/1/9

/24±20

:~

00

20

~

I 2547

0/256/257

0 /474 /475 0/254/265 0 /2761124/250

!

29±4

40

7±2

~

!-490~.-49Ci8

I 6

1

18

2

0 0

-8

~ -4

-6

-2

0

-J (1~

(c)

I

I

I

Clinton

2

I

(d)

8

20±04

8 —a.-

4,

10

—12

I Recife.Braz,I

l_~

‘~
Y 840855

,~

Y~54/055 V 1/75 1/77

4

(FAF&WXE /976)

I I I -6 -4 -2 0 3 Fig. 2. Examples of relative sea level Years data from one site in each of the characteristic regions. The hatched region incorporates the from present ( x i0 experimental error due to the mapping from 14C age to sidereal age (Svess, 1970; Stuiver, 1971). -8

I

-6

I

-4

-2I

8

-8

289

Time ( yx103

Fig. 3. Comparisons of predicted and observed relative sea-level data for six sites under the northern part of the Laurentide ice sheet. The hatched regions represent observations. The dotted, dashed, and solid curves are predictions for models with lower mantle viscosity of 102’ Pa s-‘, 1O22 Pa s-l, and 5X 1O22 Pa SK’, respectively.

s-’ and 1O22 Pa s- ’ lower mantle viscosities. Inspection of curves Ll and L2 suggests that the actual Earth model is intermediate between these two. In order to account for initial disequilibrium effects we require knowledge of the loading prehistory and this may be estimated on the basis of oxygen isotope data taken from deep sea sedimentary cores. These data reveal the saw-tooth history of the load cycle shown in Fig. 5 which is characterized by a timescale of lo5 y between successive interglacials. A more detailed discussion of these ideas will be found in Peltier (1982) and Wu and Peltier (1982a, b).

1

Fig. 4. Relative sea-level predictions for a parabolic load deglaciation history compared with relative sea level data from three sites near the centre of the Laurentide ice sheet. Models Ll and L2 have lower mantle viscosities of 102’ Pa s-i and 10z2 Pa s-‘, respectively. The solid rectangles, circles and triangles show the relative sea-level data points from Ottawa Island, Churchill and Castle Island, respectively. For each of these models the solid curves are calculations made subject to the assumption that isostatic equilibrium prevails initially whereas the dashes curves include the influence of a degree of disequilibrium implied by particular viscosity and model.

KYRBP

TIME

Fig. 5. (a) The variation with time of the oxygen isotope ratio in a typical deep-sea sedimentary core. (b) Sawtooth approximation to the load history based upon isotopic data.

290

2.2. Free-air gravity data Free air gravity maps for the Laurentide and Fennoscandia regions are shown in Fig. 6a and b. The Laurentide map is based upon the data of Innes et al. (1968) as compiled by Walcott (1970) and was constructed by averaging the observed free-air anomalies in 1° X 2° grid elements, each of which typically contained 100 observations over land and —40 observations over water (Hudson Bay). This averaging procedure effectively filters fluctuations with length-scales less than 400 km so that influences due to near surface geological structure and surface topography are removed. The major anomalous structure which remains consists of an elongated northwest trending elliptical trough with peak amplitude between —30 mgal and —40 mgal which correlates perfectly with the topography of the Laurentide ice load, The map for Fennoscandia (Fig. 6b) is from Balling (1980) and is based upon the data of Honkasolo (1964). The data from this region have historically proven to be rather difficult to interpret because the spatial scale of the loaded region is not sufficiently separated from the scales on

.

which near surface geological and surface topographic effects matter. Balling’s reduction of the data is the most useful which is presently available. The map shown in Fig. 6b was obtained by removing from the raw data that part of the anomaly which was linearly correlated with the topography and shows a residual which is very well correlated with the maximum of the original ice sheet topography. On the basis of the fact that the + 5 to + 10 mgal contour coincides with the location of the ancient ice margin, rather than the 0 mgal contour which would obtain if the anomaly were unbiased by larger scale variations of different cause, one may argue that the peak anomaly associated with glacial isostatic disequilibrium should be 15 mgal to —20 mgal. In Fig. 7 we show theoretical predictions of the free-air gravity anomaly which should exist at present over Hudson Bay for two assumed profiles of mantle viscosity. Both calculations have been done on the basis of the assumption that isostatic equilibrium prevails initially. That shown in Fig. 7a is for the model with uniform mantle viscosity of 1021 Pa s~and that in Fig. 7b is for the model with lower mantle viscosity of 1022 Pa s The —

Fig. 6. Free-air gravity anomalies (in mgal) observed over the Laurentide and Fennoscandian depressions.

~.

291

Fig. 7. Predicted free-air gravity maps for the Laurentide region: (a) is the prediction for the uniform viscosity model (10 21 Pa s (b) is for the model with lower mantle viscosity of 1022 Pa s~.

former predicts a peak anomaly of about —30 mgal (which is marginally acceptable) whereas the latter gives a peak anomaly near —60 mgal which

-12C =

is unacceptably high. The effect of initial isostatic disequilibrium on the predicted free air anomaly, however, is quite large and therefore must be

I

-100-

(a)

I);

-60

I

LAURENTIDE

(b) FENNOSCANDIA -

-

g >-

..-.

-80-

~‘

~

.

a.-.--

~

-60

-

,__

~

~ — —

~

•

-30

~77~7~7-20~

20— 0 ~Q22

.

,1I1lpII1

I023

11LOWER MANTLE

lit

11111

(POISE)

I

I024

.400

-

I022

II

1111111

illil

023

LOWER MANTLE

I024

0

(POISE)

Fig. 8. Predicted present-day free-air anomalies using paraboloid disc load approximations to the Laurentide and Fennoscandian ice sheets. The observed peak anomaly is shown as the hatched region. The dashed—dotted, dashed and solid lines are repectively predictions assuming initial equilibrium, including the effect of disequilibrium but with fixed ice sheet radius, and including the effect of disequilibrium with time varying ice sheet radius; the solid curve includes the effect of initial disequilibrium and allows for the continuous variation of ice sheet radius during loading and unloading (see Wu and Peltier, 1982a, b for details). Predicted anomalies are shown as a function of lower mantle viscosity.

292

taken into account when this datum is employed to constrain the viscosity profile. Figure 8 shows the peak free air anomaly predicted for Laurentide and Fennoscandia (Fig. 8a. and b, respectively), on the basis of disc-model deglaciation histories, The present-day anomaly is shown as a function of the viscosity of the lower mantle and three different predictions are included. The glaciation prehistory is constrained by the oxygen isotope data shown in Fig. 5, and the ice sheet is assumed to maintain an equilibrium plastic profile when its radius is allowed to vary. Inspection of Fig. 8 demonstrates several important points. Firstly, it shows that the effect of initial disequilibrium is to reduce the predicted free-air anomaly. The extent of this reduction, of course, increases with the lower mantle viscosity. It is interesting to note that when the variation of ice sheet radius as a function of time is included in the disequilibrium calculation the predicted reduction is considerably reduced. Accepting the solid lines as the most accurate of the free air anomaly predictions, the Laurentide data require 2 X 1021 Pa s~’<~‘ LM < ~ X 10~l Pa s’ whereas the Fennoscandia data suggest 102! Pa s~
therefore agree as to the necessity of a modest increase of viscosity in the lower mantle. Furthermore, this increase is in accord with that expected if the temperature variation across the phase boundary is adiabatic but rules out the presence of substantial thermal boundary layer unless the increase in creep activation energy across the phase boundary is 45 kcal mol~,or larger. 2.3. Polar motion and day-length observations A very useful additional data set exists which can be employed to further check the validity of the mantle viscosity profile inferred from the relative sea level and free air gravity data. The first of these additional observations is shown in Fig. 9

and consists of the ILS—IPMS record of the monthly averaged position of the rotation pole relative to the surface geography during the time interval since about 1900 A.D. The net variation is separated into components x and y relative to the conventional international origin (CIO) which are displayed separately. The inset polar projection shows the orientation of the x and y axes relative to the continents. Although the polar motion is dominated by the 7y beat between the 14 month Chandler wobble and the 12 month annual wobble, the motion which is of interest here is the secular drift upon which these oscillatory motions are superimposed. This drift is at a rate near 1°/106yand in the direction shown by the arrow in the polar projection, i.e. towards Hudson Bay and therefore towards what was the centroid of the ancient Laurentide ice sheet. Although it has often been suggested in the past that the secular drift of the pole was driven by sea-floor spreading (e.g. Dickman, 1979) it has recently been shown that this secular drift is, in fact, forced by the same deglaciation event responsible for the previouslydescribed r.s.l. variations and free-air gravity anomalies (Sabadini and Peltier, 1981; Peltier, 1982). Using the theoretical methods developed in Peltier (1974) to analyze surface loading problems for viscoelastic models, Sabadim and Peltier (1981) showed that the perturbations of planetary inertia forced by the deglaciation event were proportional to the factor (1 +k 2(s)) where k2(s) is the frequency domain form of the second-degree Love number and s is the Laplace transform variable. Solution of the equation for angular momentum conservation of the deformable body is then straightforward and leads to a prediction of the variation of polar wander speed and direction as a function of time for a given deglaciation model and mantle viscosity profile. In Fig. 10 we show the prediction of polar wander speed as a function of time for a disc load deglaciation model of the Laurentide ice sheet and for a homogeneous model of the planetary interior. The observed polar wander speed is shown by the vertical bar at the present time measured relative to a zero taken at the end of the model deglaciation phase. Model predictions are shown for several different strati-

293

PACIFIC O~~EAN

5I~’

~

ATLANTIC OCEAN

47Q11

-

FENNOSC ND/A

-

~

(x) _470h

1900

I

I

I

YEAR

I

I

1975

Fig. 9. Components x andy of the polar motion during the period 1900—1975 based upon ILS—IPMS data. The arrow in the central polar projection shows the direction of the secular part of the polar motion upon which the Chandler and annual wobbles are superimposed.

fied viscoelastic models which differ from one another only in the thickness of their surface lithosphere, as noted in the figure caption. All models have uniform mantle viscosity of 1021 Pa s~ and an elastic structure identical to model 1 066B of

Gilbert and Dziewonski (1975). Extrapolation on the basis of these results (Peltier, 1982) shows that the observed polar wander speed and direction will fit with a model whose lithospheric thickness L is 260 km. Because this observation is as sensitive

294

to lithospheric thickness as it is to the viscosity of the mantle, it cannot be employed to unambiguously constrain the latter. This analysis is contrary to that of Sabadini and Pettier (1981), which was based upon an ‘invalid approximation, as noted in Peltier (1982). There does exist a second observation, however, which does not suffer from this ambiguity of interpretation and which can be used to provide an independent verification of the inference of uniform mantle viscosity which was based upon the relative sea level and free air gravity data. This is the astronomically observed non-tidal acceleration of the Earth’s rotation. Although the length of day is generally increasing due to the torque exerted on the Earth through the agency of the ocean tides, the observed rate of increase is less than would be expected on the basis of the known present-day value of this torque. The discrepancy can be understood as a consequence of the ongoing relaxation of the component of planetary oblateness, which was forced by the accumulation of ice on the polar continents during the last glaciation. In Fig. 11, we compare the predictions of several layered viscoelastic models with the observed

1.5

I I

I

non-tidal acceleration. The characteristics of the models are noted in the figure caption. Detailed analysis of this data in Peltier (1982) and Wu and Peltier (in preparation) shows that it exhibits a preference for the same uniform viscosity mantle as that required by the relative sea level and free air gravity data. On the basis of arguments presented in the introduction, the observation that the viscosity of the mantle is essentially uniform may be taken to imply that the convective circulation is not layered but fills the entire mantle. In the next section we will investigate some properties of the circulation itself.

3. Aspects of the mantle convective circulation: boundary layer structure and the rate of internal heating Convection in the mantle is a high Rayleigh number thermally-forced circulation which may be described to first order in terms of the equations of Newtonian hydrodynamics. These consist of the following partial differential equations which respectively express the pointwise conservation of mass, momentum, and energy, and the equilibrium relation between the thermodynamic variables (12)

1.2

-

E~0.9’ ~

3

P(-~7+u.vu)

-

I

I

~c~(~+u.vT)

-

-

0.3

-

~

—vP+pg+v’r

(13)

(kvT)+H

(14)

=~.

P_Pd[l~t(TTd)] (15) The symbols in (12)—(lS) have their conventional meanings. Introduction of the high Rayleigh number scaling first discussed in Peltier (1972) reduces the field equations to the respective forms

~

0.6

=

-

_____________________________ 00

0

I

2

4

6

I

8

(16) 10

TIME(KYRS.) Fig. 10. Observed speed of polar wander and predictions for three different Earth models which differ from one another only in lithospheric thickness. Models 1, 6 and 7, have lithospheric thicknesses of 120, 195 and 245 km, respectively.

Rdul

2

—

PrPdt~P0~+V

dT dt p

=

—

[1

dP p dt —

ô( T

—

—

~ +

pR ‘- V T + h,

Td)]

(17)

U 2

—

—

e,j e,1

18) (19)

295

where the Prandtl number Pr is Pr = v/ac. The Rayleigh number R is: R = ga ~Td3/acv, if heated below; and R = gad5(H/pc~)/ac2v if heated within, where d is the vertical depth of the layer, L~Tthe temperature drop across it, and H the rate of internal heating. The dissipation number T is T =

I

I

I

I

-

-

3-

-

2

-

o -~

U)

gad/c p

>-

Since the thermal diffusivity ac 102 cm2 s’ for mantle material and the kinematic viscosity 1021 cm2 s1 the Prandtl number Pr 1023, and therefore the inertial force on the left-hand side of (17) is completely negligible for any reasonable value of the Rayleigh number. High Prandtl number flows are characterized by the absence of any pronounced shear in the velocity field even though ~

~

4( F + Hd )/pocp ,‘2v Rand (gad the fractional heating from within

(20)

Hd/F+ Hd (21) in which F is the heat flux through the lower boundary of the convecting region so that in steady state F + Hd is the total heat flux through the upper boundary. In Fig. 12 we show solutions for the two-dimensional temperature field in aspect ratio 1 boxes which are heated entirely from below at Rayleigh

ii. =

-

~

C 0

the temperature field may have a well developed boundary layer character. Jarvis and Peltier (1982) have recently described a sequence of solutions of the infinite Prandtl number, zero dissipation number version of eqs. 16— 19 which are expected to provide a reasonably accurate first order description of the mantle convective circulation. Their solutions are the highest Rayleigh number numerical analyses yet performed and in particular are the only calculations to achieve those values of the Rayleigh number which we expect to be characteristic of flow in the Earth. Their analyses, furthermore, are unique in that the properties of the circulation have been explored as a function of the fraction of the surface heat flow which is generated within the fluid (i.e. simulating the influence of radioactive heating). These simulations are characterized by two nondimensional parameters which are the Rayleigh number

C L

-

2

4

6

8

10

TIME (KYRS) Fig. II. Observed non-tidal acceleration andone predictions of four different Earth models, which differ from another only in the value of the viscosity beneath 670 km depth. These values are 1021,3.1021, 1022 and 1023 Pa ~‘-I for models 1,2, 3 and 4, respectively. Results are for forcing by the Laurentide ice sheet only. Including Fennoscandia and Antarctica lowers the implied lower mantle viscosity to the value mentioned in the text.

numbers which are l0~and iO~times the critical value R~which marks the onset of convection. On each frame the dashed line is a single reference isotherm and the solid lines are vertical and horizontal cross-sections through the temperature field at the corresponding vertical and horizontal levels. Inspection of these plates clearly shows that at high Rayleigh number the flow is dominated by sharp thermal boundary layers adjacent to the horizontal boundaries and by thin thermal plumes which bound the cell laterally. Awaylayers, from the vertical plumes and cores horizontal boundary the cell core is everywhere close to isothermal (adiabatic). An important question which we wish to address here is whether there are geophysical data which might be used to test the model of the interior temperature field suggested by the high Rayleigh number simulations of convection. There are in fact several such observations which we will proceed to discuss below. The first of these is from Anderson (1981) and consists of upper mantle temperature profiles in-

296

—

~

~

~

______

______

______

~

~

T(x)

_____ ____

!

Ti

______I

________ il II II

__________

/

1~z)——.

/

ii II

/ I I

_____ I

I

~~7”_’~’~\ I’ / _____ ____

I

I

______

________________________ ____________________________ _____________________________I

I

_____

I

~

—

______

______

I

_____

II. I

/

1

T(x)

II

____ 1

_____ ____

_____

I/

_____I

I’

I0~R~’ Fig. 12. Internal temperature fields and vertical and horizontal cross sections for heated from below convection at: (a) R = l0~R~ and (b) R = l0~R~.Note the ‘overshoots’ of the vertical temperature profiles above and below the isothermal core temperature.

ferred from the seismic structure of models 1 066AB of Gilbert and Dziewonski (1975). The profiles shown in Fig. 13 were computed on the basis of different assumptions concerning composition and the degree of partial melting. It will be noted that these profiles are characterized by an ‘overshoot’ of the temperature above that which obtains at depths of 300—400 km. Comparison with Fig. 12 shows that this overshoot is entirely expected on the basis of convection theory if we assume as suggested in Peltier (1980) that the oceanic lithosphere is the cold thermal boundary layer of the mantle convective circulation. Furthermore if mid-mantle temperatures are 1400—1500°C then the magnitude of the overshoot is also predicted to be 100°C to 200°Cas implied by the seismic data. Further support for the notion that the oceanic lithosphere should be interpreted as the surface boundary layer of the convective circulation has recently come from the study of magnetotelluric data taken on the sea floor. A recent analysis of the small set of such data which currently exists, by Oldenburg (1981), has concerned the inversion of the induction data to give electrical conductivity depth profiles as a function of age. These data were in turn employed to infer temperature vs. depth profiles on the basis of a particular model of partial melting. In Fig. 14 we show a reproduction

160C

-

1500 (a)

~ ~ = ~

(b) 1450 _______________

1400

-

___________

a E

~

1350

1300

(a) Eclogite. dry (increasing garnet) (b) Garnet lherzolite dry (increasing garnet) (c) b with slight partial

1250

1200

meft in

LVZ

-

(d) bwith

100

large partial melt in LVZ

200

300

400

DEPTH 1km)

Fig. 13. Geotherms based upon four different models of mantle composition and the elastic structure from modelare: l066A of Gilbert and Dziewonski (1975). The compositions (a) dry eclogite; (b) dry lherzolite; (c) lherzolite with a small amount of partial melting; (d) lherzolite with a large amount of partial melting. After Anderson (1981).

297 RESISTIVITY Lim 01

0

50 100

02

3

I0

I

lO

IO~

0

1~2,l?3

( \ \

,....‘

(.~

\

200 250

% MELT

300

o0 2

4

6 8 0

~

~

hot plume rising in the mantle beneath the ridge crest. Models of plate tectonics which are is a

I

/J

/ / /

50 ~

02

~

8

10

50

k 100 50 200

based the the assumption thatthe there is noand correlation upon between motion of plates the large scale circulation of the mantle are incapable of explaining the characteristic overshoot of the sub lithospheric geotherm and the observation of the overshoot may be taken to imply that these models are invalid. The surface of the Earth is fully involved in the mantle general circulation. The fact that at least the oceanic lithosphere is an intrinsic part of the mantle convective circulation implies that its properties may be diagnostic of the way the circulation is forced, that is whether the motion is maintained predominantly by heating from below or predominantly by heating from within. At a Rayleigh number R = 6 l0~,Fig. 15 X

Fig. 14. Possible models of resistivity (upper) and inferred proportion of melt (lower) beneath three ages of ocean floor in the Pacific: (a) 72 Ma, (b) 30 Ma, (c) 1 Ma. After Oldenburg (1981).

(after Jarvis and Pettier, 1982) shows the bathymetry vs. distance (age) curves and the heat flow vs. age curves for several values of the fractional internal heating ~s.For heating entirely from within (~s= 1) the bathymetry shows no variation with age above the region of rising motion (x = 0) since in this case the flow contains no hot lower boundary layer and corresponding rising plume.

of the three resistivity and partial melt vs. depth profiles inferred by Oldenburg for sea-floor ages of 72, 30 and 1 Ma. These data imply temperature vs. depth profiles which are also characterized by a sharp overshoot beneath the surface lithosphere of the temperature above the temperature at greater depth. They furthermore imply that the magnitude of the overshoot decreases with age and that the depth to the temperature maximum increases with age. Comparison with Fig. 12 shows that even this detailed structure is predicted by the numerical convection model. The numerical model of thermal convection predicts the existence of a temperature overshoot beneath the thermal boundary layer because of the presence of the hot plume coincident with the rising limb of the heated from below convection cell. If the oceanic lithosphere is associated with the thermal boundary layer of the mantlewide convective circulation we then expect a temperature overshoot beneath the lithosphere only if there

Only the cold top boundary layer and corresponding sinking plume exist and find expression in the surface bathymetry. As .s 0 and heating from below becomes progressively more important the bathymetry profile gradually assumes a ~ë dependence above the rising stream and across most of the width of the cell (%/ã~ë reference curves are shown as dashed lines adjacent to each of the individual curves). Figure 16 shows the same data for a fixed value of ~tt = 0.2 as a function of the Rayleigh number. For the Earth the Rayleigh number is near tO7. As the Rayleigh number increases for fixed p. the initial a~ë~ decay of topography becomes more extreme and for R l0~is able to deliver the 3 km of relief between the oceanic ridge and deep basins. Also evident on this Figure is the fact that for sufficiently great age the bathymetry is predicted to flatten away from the ~‘ä~ëreference curve. In the limit p. 0 this flattening entirely disappears. Since the sea floor is observed to flatten away from the %/ä~curve for ocean floor age greater than about 70 myr (Parsons and Sclater, 1977) this suggests that this

250

(a)

(b)

Cc)

—~

—

—p

298

1 (a)

,RR=6.9x105

0.5

0

,

(b)

/

1.0

RR = 6.011 IO5

I

o-

,

05 X’--

0

XI-

IO

Fig. 15. Horizontal variations of: (a) bathymetry; (b) surface heat flow; and (c) scaled bathymetry across mantle convection cells at Rayleigh number, R =6X IO5 and a range of c as indicated on each plate. Curves are offset at the origin for clarity. Solid curves are model predictions and dashed curves are best fit ,/L@ curves constrained at the points indicated by the arrows.

1982) shows our attempt to do this. On the domain diagram defined by the parameters R and ~1 we have plotted contours of constant percentile flattening- of the bathymetry profile. Since the

observation might be employed to infer the value of n and therefore the fraction of the surface heat flok which is due to internal heating. Figure 17 (from Jarvis and Peltier 1980, 1981,

(b)

E g

t ._ ‘i

RR =1.25x IO’

1,5_

3

a

l.O-

I

2.00

0.5 x’ -

Fig. 16. Same as in Fig. 15 but for fixed p =0.2

OoI and variable

I.5

xRayleigh

I 0

1 Q5 x’-

number

R.

I I.0

299

—‘

I

I

I

I

I

2.\\ io

T ~

I0~.~\

~\ \

\

8

~ 64.-

‘

Flow

~

\

~ lo~.\

\ \

5s

64 4

\\\

\ \~~\ \\\ o6. I

-

6 5xl0.~

Flow

~ .2

6

~

.4 .6 ~U. ~

8

~

.8 1.0

.1

-

~

\ \\ ~

\/.2

~

2

-

4. The cooling of the Earth and crustal differentiation

I

~~en~ent

~

8

‘o ~o\ 2

\ ~

2\

Dependent

,~

I

~

Whole Mantle

‘

~o~°\

\\ \

\

1~

I

i .2

~‘°~

,

~

.4

.6 .8 1.0

/.L “k

Fig. 17. Domain diagrams showing contours of constant percentage flattening as a function of R and ~sbased on: (a) an exact calculation of bathymetry while plate; (b) an approximation discussed in Jarvis and Peltier (1982). The broken diagonal roughly divides domains of steady and unsteady single cell flow. The horizontal dashed line at R = 2.5 X l0~indicates our estimate of the value of R required in whole mantle convection to fit the observed mean surface heat flow,

mantle Rayleigh number must be near lO~to explain the observed mean surface heat flow and since the percentile flattening of the sea floor is 20% (Parsons and Sclater, 1977) we need p. 0.2 to deliver the observed flattening. This suggests that the mantle convective circulation may be forced to a very substantial extent by heating from below. It should be kept clearly in mind that this inference is based upon two major assumptions: (I) that the convective circulation is whole mantle in vertical scale; and (2) that the oceanic lithosphere is the thermal boundary layer of this circulation. The implication that the mantle circulation does not contain sufficient radioactivity to sustain itself is in good accord with the a priori geochemical argument mentioned in the Introduction (e.g. O’Nions et al., 1979) that the Earth contains far too little radioactivity to explain the observed surface heat flow. The heat flow deficit must obviously be made up by secular cooling of the planet in bulk and this is the subject which we turn to in the next section.

A considerable revitalization of the field of thermal history modelling has taken place in the past five years due to the development of a scheme for the evolution effect of thermal convection ‘parameterizing’ upon the thermal which does not require explicit solution of the full set of hydrodynamic equations. This was first introduced by Sharpe and Pettier (1978, 1979) and is based upon the empirical relation which obtainsnumber betweenwhen the Nusselt number and the Rayleigh the convective circulation is vigorous The Nusselt number Nu of a convective circulation is simply a measure of the heat transport nondimensionatized by a measure of the conductive transport. For heated below convection

(

(d \ d(T> ~j~7”/ k dz

Nu

\

zd

122

/

where ~T> is the depth-dependent horizontal average of the temperature. In Fig. 18 we show a plot of Nu, of characteristic horizontal and vertical velocities ii and W of the circulation, and of two measures of boundary layer thickness 61, 6~all as functions of the Rayleigh number (after Jarvis and Pettier, 1982). For R ~ 50 R~each of these variables follows a power law relation of the form 6 = fi, (~ ) d( R B/R )SI (23a) ti=

/J 2(L~)(ac/d)(RB/RC)2 ~ ‘R

(23b) ~23 ~ c (23d)

—= ~ ~ W

~

~

~

fl/

Cl

Nu = $4(~)(RB/RC)~~~

q=$4(L~)(acL~T/d)(R~/RC) (23e) where $ = (0.3, 0.118, 0.093, 2.30), S = 0.278, 0.645, 0.684, 0.313), and ~ is the aspect ratio of the circulation. Pettier (1980, 1981b) has shown how these relations may be employed to construct an argument for the whole mantle convection model under the assumption that the flow has the usually preferred aspect ratio of i~ = 1. Here we will discuss the way in which empirical relations (23), and similar relations to them for flows which are partly heated from within, may be employed to analyse the thermal history. (—

300 I

I

I

I

I

I

I

1

60’

I

I~

lO~

RB/RC Fig. 18. Power law relations for the Nusselt number Nu, the boundary layer thickness 8’, and characteristic horizontal (0) and vertical (W) velocities in high Rayleigh number heated from below convection with free—free boundary conditions.

Although some qualitative notions of how (23) might be so employed have appeared in Tozer (1967, 1972, 1974) and have been further explored by McKenzie and Weiss (1975), none of these articles made any attempt to construct a rational thermal history. Tozer’s argument has always been essentially that the circulation will adjust itself to deliver a geotherm such that the interior maintains a uniform value of the viscosity near 1021 Pa s~1. Although this idea of the fundamental role played by the negative feedback due to the dependence of viscosity upon temperature is useful, it does not answer the question crucial to the thermal history which concerns the time scale on which the mean mantle temperature changes. How long does it take for an initially hot system, for example, to cool into the state with i = 1021 Pa s’? Sharpe

characteristic time for the mantle convective circulation is much shorter than the timescale on which significant changes in mean temperature occur, the effect of convection upon the evolutionary history of the mean temperature may be ‘parameterized’ in a way which does not require explicit solution of the hydrodynamic equations. They derive a model of mantle ‘climate’ in which the details of the convective circulation are ‘weather’ which may be safely ignored insofar as the long timescale variations of the mean temperature are concerned. This model is simply based upon separate energy balance equations for the mantle and core which have the respective forms (e.g. Peltier, 1980, l98lb) 8 çrp 2 2 2 PmcPmr ~1r r
and Pettier (1978, 1979) argued that this question could be answered rather simply because the thermal history problem is governed by two quite widely separated timescales. To the extent that the

+ (r,~ ‘~ 2q~ + r~3Q~ ~-

p0c~‘~-

f

~r2<1~>dr

=

‘~

j-

—

)

Qm

(24) (25)

301

in which
Q~

Q~

(Q

m

= 0 or all heat sources contained in a thin ‘crustal’

layer at the surface), then the steady state heat flow across the CMB would be

j

T 1 (rC + r~) (r~ r~) ~ 4r2 2 1I C S Nu

(T q~—k

— —

)

(26)

where km is the effective thermal conductivity of ~the mantle and T~and T 5 are the temperatures at the core mantle boundary and the Earth’s surface, respectively. The Nusselt number is given by the empirical expression (23d) in which the Rayleigh number is defined as 3(~T— ~ (27) ~ R = ga(r~ r~) in which ~ is the adiabatic temperature drop across the mantle and L~Tis the actual temperature drop. In Sharpe and Peltier (1978, 1979) the effective mantle viscosity ~m was represented as —

i”

to chemical differentiation of the mantle can also be accommodated in the model through appropriate adjustment of the power law exponent. In order to run the above thermal history model we require an initial geotherm for in thethe time immediately subsequent to advances the Earth’s formation. Rather important recent understanding of accretion processes (e.g. Wetherill, 1976; Safronov, 1981) have fortunately established that the accretion time scale is -~108 y (a small fraction of the Earth’s age) and that due to the size of the accreting bodies the planet forms sufficiently hot to lead to the gravitational differentiation of iron from silicates and thus to large further increases of temperature. Thisassumed picture isinin with the initial thermal states theaccord thermal histories of Sharpe and Peltier (1978, 1979) in which the core consists of liquid iron and the mantle has been cooled to the solidus by freezing outwards from the CMB, due to the efficiency of convective cooling in the partially molten state

exp( gT

(28) using the previously discussed empirical expression of Weertman (1970), in which Tm is a midmantle temperature and 7, is the melting temperature at the same depth. The introduction of the temperature dependence of viscosity in this way incorporates the important negative feedback which is crucial to the thermal history. If Qm = 0 then (24) may be rearranged to give the surface heat flow as 0/Tm)

=

so that the surface heat flow includes a contribution from the secular cooling of the mantle on the long time scale T. Aside from the derivation of this parameterized convection model, the most important result in Sharpe and Peltier (1978, 1979) was the demonstration that secular cooling could play a significant role in driving the mantle convective circulation. An additional consequence of this secular cooling is a slight modification of the power law exponent in (23d) such as to impede the cooling rate. Circumstances in which Qm ~ 0 and in which Qm is a function of i- due both to decay of the strength of the radioactivity with time and

+ PmCpm~ f’r2
(29)

and the fact that the mantle solidus is everywhere superadiabatic. In Fig. 19 we show the evolution with time of the mean mantle viscosity, the mantle Rayleigh number, and the mean mantle heat flow for three different thermal history models with parameterized convection which have been run for 4.5 Ga. The curves marked =0 are for a model which has no internal radioactivity. Inspection shows that this model cools in the first I Ga into a very cold quasi-steady state in which the mean mantle viscosity is -~1023 Pa s~ (1024 Poise) Q

which is obtained 2 orders from of magnitude thandata. the viscosity postglacialhigher rebound

302

total complement of radioactivity equal to half the and which present day heat slowly flow (—40 mWthe m2) is differentiated to form crust over the

010

I

loll

I0~

~

~ I0~

~ 10~

“

/ — — ~/‘~‘~

51022

~ ~ I0

/

I~IOe

1020

age of the Earth. The presence of radioactivity in the mantle also reduces the cooling rate since the

I IO’’

——

/

~I0I5

heat flux across CMB is lower for a given temperature drop the across the mantle. The main idea to be discussed here concerns the

02 lllllllll~

______________

TiME (Ga)—.-

~

4.8

24(,

T~E Ga—.- 4.0 48

_____________

50 AREA FRACTION

________________________

210

~o

-

0~

00

180 i 150

20

~ Z

0M~0

— —— —

60 30 .5

10

1.5

20

25

TIME I Gal

35

~

45

~

—

Fig. 19. Variations of mantle Rayleigh number, mean mantle viscosity, and surface heat flow as a function of time for three different thermal histories with parameterized convection. The arrows at about 3.8 Ga indicate the time at which continental crustal matenal is first stabilized at the Earth’s surface based upon the age of the oldest known surface rocks,

Furthermore, the present day surface heat flow is only about 15 mW m2 compared to the observed mean surface head flow of about 80 mW m2. In order to inhibit this too rapid cooling we must add internal heat sources. The curves marked ~ 0 are for a model in which the mantle contains no radioactivity but the core contains 0.2 wt. % potassium. Although this model must be considered unlikely since current geochemical evidence suggests that the core should be devoid of radioactivity, it nevertheless predicts a present day state which is compatible with the observations. The non-monotonic nature of the evolution is due to the fact that the core heat source concentration is so high that the heat transferred by mantle convection cannot initially keep us with the generation. The preferred history is that which follows the dashed curves which are for a model which has a

Q~

meaning of the relaxation time of 1 Ga which these models predict for the cooling history. The models without core heat-sources (aside from the latent heat of freezing of the inner core) both take . . 1 Ga to cool into a high viscosity quasi-steady state. The arrows on each plate mark the time in the thermal evolution at which continental crust appears first to have been irreversibly differentiated from the mantle. This is determined by the age of the Isua sequence of about 3.8 Ga (e.g. Moorebath, 1975). The area fraction vs. time curve in the third plate shows a typical crustal growth curve from Taylor and McLennan (1981). The thermal history calculations discussed here very strongly suggest a close connection between the ability of the mantle reservoir to irreversibly differentiate crustal material and the vigour of the mantle convective circulation. Not until the system .

.

has cooled into the high viscosity quasi-steady state is it possible to stabilize continental crust. In the first 1 Ga the secular cooling of the interior is so rapid that the circulation must be highly transient in space and time. Only if rising and sinking regions have strongly time dependent locations will it be possible to effect the required uniform secular cooling of the interior of the mantle. In this phase mixing is extremely efficient and crustal material is simply reingested into the main flow. It is interesting to note that the whole mantle convection model predicts this relaxation time correctly. On the basis of the discussion in the Introduction it is quite clear that a layered convection model would predict a much longer characteristic relaxation time for the cooling history.

5. Conclusions With the possible exception of the isotopic abundance data mentioned previously, most geo-

303

physical observations are best explained by the whole mantle model of the convective circulation. Furthermore, since evidence for the necessity of a chemical discontinuity coincident with the phase transition at 670 km depth is unconvincing (Jeanloz, 1981) and since it has already been established that the phase change itself cannot inhibit penetration by convection (Richter, 1973; Christensen, 1982) there appears to be no physical mechanism which could prevent mechanical mixing across this boundary. If convection in the lower mantle were decoupled from1 convection in the upper mantle then there would be a thermal boundary layer at 670 km depth which should be evident in the viscosity vs. depth profile; but there is no evidence in this parameter for any boundary layer structure at this depth. If the convective circulation were layered then the efficiency of convective heat transfer would be much reduced from that for the whole-mantle model and it is not clear that one could then fit the thermal history constraints. Again, on the basis of the best observations available, it would appear that the whole-mantle model of the convective circulation is strongly favoured.

References Allegre, C.J., 1982. Chemical geodynamics. Tectonophysics, 81: 109—132 Anderson, DL., 1981. Hotspots, basalts and the evolution of the mantle. Science, 213: 82—89. Anderson, 0., 1981. Temperature profiles in the earth. In: Ri. O’Connell and W.S. Fyfe (Editors), Evolution of the Earth. Geodynamics Series Vol. 5. Am. Geophys. Union, Washington, DC, pp. 19—27. Balling, N., 1980. The land uplift in Fennoscandia, gravity field anomalies and isostasy. In: N-A. MUmer (Editor), Earth Rheology, Isostasy and Eustasy. John Wiley, New York, NY. Carlson, R.W., Lugmair, G.W. and MacDougall, J.D., 1981. Crustal influence in the generation of continental flood basalts. Nature (London), 289: 160—162. Cathles, L.M., 1975. The Viscosity of the Earth Mantle. Princeton University Press, Princeton, NJ, 386 pp. Christensen, U., 1982. Phase boundaries in finite amplitude mantle convection. Geophys. JR. Astron. Soc.; 68: 487—497 Davies, G.F., 1977. Whole mantle convection and plate tectonics. Geophys. JR. Astron. Soc., 49: 459—486. Davies, G.F., 1981. Earth’s neodynium budget and structure and evolution of the mantle. Nature (London), 290: 208— 213.

DePaulo, D.J., 1981. Nd isotopic studies: some new perspectives on Earth structure and evolution. EOS, Trans. Am. Geophys. Union 62: 137—140. Dickman, S.R., 1979. and true polar wandering. Geophys. J.R. Continental Astron. Soc.,drift 57: 41—50. Dziewonski, AM. and Anderson, D.L., 1981. Preliminary reference earth model. Phys. Earth Planet. Inter., 25: 297—356. Gilbert, F. and Dziewonski, A, 1975. An application of normal mode theory to the retrieval of structural parameters and source mechanisms spectra. Philos. Trans. R. Soc. London, Ser. A,from 278:seismic 187—269. Goldreich, P. and Toomre, A., 1969. Some remarks on polar wandering. J. Geophys. Res., 74: 2555—2567. Haskell, NA., 1937. The viscosity of the asthenosphere. Am. J. Sci., 33: 22—28. Honkasolo, T., 1964. theupheavel use of gravity measurements for investigation of theOn land in Fennoscandia. Fennia, 89: 21—23. Innes, M.J.S., Goodacre, AK., Weston, A. and Webber, JR., 1968. Gravity and isostasy in the Hudson Bay region. In: CS. Beals (Editor), Science, History and Hudson Bay, Part V. Jacobsen, SB. and Wasserburg, G.J., 1979. The mean age of mantle and crustal reservoirs. J. Geophys. Res., 84: 7411— 7427.

Jarvis, G.T. and Peltier, W.R., 1980. Oceanic bathymetry profiles flattened by radiogenic heating in a convecting mantle. Nature (London), 285:Peltier, 649—651. Jarvis, G.T., and W.R. 1981. Effects of lithospheric rigidity on ocean floor bathymetry and heat flow. Geophys. Res. Lett., 8: 857—860. Jarvis, G.T. and Peltier, W.R., 1982. Mantle convection as a boundary layer phenomenon. Geophys. JR. Astron. Soc., 68: 389—427. Jeanloz, R. and F.M. Richter, 1979. Convection, composition and the thermal state of the lower mantle. J. Geophys. Res., 84: 5497—5504. Jeanloz, R., 1981. Phase transitions and mantle discontinuities. Rev. Geophys. Space Phys. (in press). Jordon, T.H., 1981. Continents as a chemical boundary layer. Philos. Trans. R. Soc. London 5cr. A, 301: 359—373. Kennedy, G.C., and G.H. Higgins, 1972. Melting temperatures in the Earth’s mantle. Tectonophysics, 13: 221. Mao, H.K., Bell, P.M. and Yagi, T., 1979. Iron magnesium fractionation model for the Earth. Carnegie Inst. Wash. Yearb. 78: 621—625. McConnell, R.K., 1968. Viscosity of the mantle from relaxation time spectra of isostatic adjustment. J. Geophys. Res. 73: 7089—7105. McKenzie, D.P., 1966. The viscosity of the lower mantle. J. Geophys. Res., 71: 3995—4010. McKenzie, D.P., 1967. The viscosity of the mantle. Geophys. JR. Astron. Soc., 14: 297—305. McKenzie, D.P., 1968. The geophysical importance of high temperature creep. In: R.A. Phinney (Editor), The History of the Earth’s Crust. Princeton University Press, Princeton. NJ.

304 McKenzie, D.P. and Weiss, N.0., 1975. Speculations on the thermal and tectonic history of the earth. Geophys. JR. Astron. Soc., 42: 131—174. McKenzie, D.P., Roberts, J.M. and Weiss, N.O., 1974. Convection in the Earth’s mantle: towards a numerical simulation. J. Fluid Mech. 62: 465—538. Moorebath, S., 1975, The geological significance of early Precambrian rocks. Proc. Geol. Assoc. London, 86: 259—279. Munk, W.H. and MacDonald, G.J.F., 1960. The Rotation of the Earth. Cambridge University Press, Cambridge. O’Connell, R.J., 1977. On the scale of mantle convection, Tectonophysics, 38: 119—136. Oldenburg, D., 1981. Conductivity structure of oceanic upper mantle beneath the Pacific ocean. Geophys. J.R. Astron. Soc., 65: 359—394. O’Nions, R.K., Evenson, N.M. and Hamilton, P.J., 1979. Geochemical modelling of mantle differentiation and crustal growth. J. Geophys. Res., 84: 6091—6101. Oxburgh, ER. and Turcotte, D.L., 1978. Mechanisms of continental drift. Rep. Prog. Phys., 41: 1249—13 12. Parsons, B. and Sclater, J.G., 1977. An analysis of the variation of ocean floor bathymetry and heat flow with age. J. Geophys. Res., 82: 803—827. Peltier, W.R., 1972. Penetrative convection in the planetary mantle. Geophys. Fluid Dyn., 5: 47—88. Peltier, W.R., 1974. The impulse response of a Maxwell Earth, Rev. Geophys. Space Phys., 12: 649—669. Peltier, W.R., 1976. Glacial isostatic adjustment. II. The inverse problem. Geophys. JR. Astron. Soc., 46; 669—706. Peltier, W.R., 1980. Mantle convection and viscosity. In: A. Dziewonski and E. Boschi (Editors), Physics of The Earth’s Interior. North Holland, Amsterdam, pp. 362—431. Pettier, W.R. l981a. Ice age geodynamscs. Annu. Rev. Earth Planet. Sci., 9: 199—225. Peltier, W.R., 1981b. Surface plates and thermal plumes: separate scales of the mantle convective circulation. In: R. O’Connell and W. Fyfe (Editors), Evolution of the Earth, Am Geophys. Union, Washington, DC, pp. 134—143. Peltier, W.R., 1982. The dynamics of the ice age earth. Adv. Geophys. 24 (in press). Peltier, W.R. and Andrews, J.T., 1976. Glacial isostatic adjustment. I. The forward problem. Geophysics JR. Astron. Soc., 46: 605—646. Peltier, W.R., Farrell, W.E. and Clark, J.A., 1978. Glacial isostasy and relative sea level: a global finite element model. Tectonophysics, 50: 81—110. ~. Richter, F.M., 1973. Finit~aniplitudeconvection through a ,phase boundary. Geophys. JR. Astron. Soc., 35: 365—376. Richter, F.M. and Johnson, C., 1974. Stability of a chemically layered mantle. J. Geophys. Res., 79: 1635—1637. Ringwood, A.E., 1979. Origin of the Earth and Moon. Springer Verlag, New York, NY 295 pp. Sabadini, R. and Peltier, W.R., 1981. Pleistocene deglaciation and the earth’s rotation: implications for mantle viscosity. Geophys. JR. Astron. Soc., 66: 552—578. Safronov, V., 1981. Initial state of the Earth and its early

evolution. In: R.J. O’Connell and W.S. Fyfe (Editors), Evolution of the Earth. Geodynamics Series, Vol. 5. Am. Geophys. Union, Washington, DC, pp. 249—255. Sammis, C.G., Smith, J.C., Schubert, G., and Yuen, D.A., 1977. Viscosity depth profile of the Earth’s mantle: effects of polymorphic phase transition. J. Geophys. Res., 82: 3747— 3761. Sauramo, M., 1955. Land uplift with hinge lines in Fennoscandia. Ann. Acad. Sci. Fenn. Ser. AlIl, 44: 1—25. Schubert, 0., 1979. Subsolidus convection in the mantles of terrestrial planets. Annu. Rev. Earth Planet. Sci., 7: 289—342. Sharpe, H.N. and Peltier, W.R., 1978. Parameterized mantle convection and the Earth’s thermal history. Geophys. Res. Lett., 5: 737—749. Sharpe, H.N. and Peltier, W.R., 1979. A thermal history model for the Earth with parameterized convection. Geophys. JR. Astron. Soc., 59: 171—205. Smith, J.V., 1977. Possible controls on the bulk composition of the Earth: implications for the origin of the Earth and Moon. Proc. Lunar Sci. Conf., 8: 333—369. Stuiver, M., 1971.in Evidence the variation atmospheric 4C content the late for Quaternary. In: of K.K. Turekian ‘ (Editor) The Late Cenozoic Glacial Ages. Yale University Press, New Haven, CT. Suess, HE., 1970. Bnstlecone pine calibration of the radiocarbon timescale 5200 B.C. to the present. In: Radiocarbon Variations and Absolute Chronology. Proc. 21st Nobel Symp., Wiley, New York, NY 652 pp. Taylor, S.R. and McLennan, S.M., 1981. The composition and evolution of the continental crust: rare earth element cvidence from sedimentary rocks. Philos. Trans. R. Soc. London Ser. A, 301. Tozer, D.C., 1967. Towards a theory of thermal convection in the mantle. In: T.F. Gaskell (Editor) The Earth’s Mantle Academic Press, New York, NY, pp. 325—353. Tozer, D.C., 1972. The present thermal state of the terrestrial planets. Phys. Earth Planet. Inter., 6: 182—197. Tozer, D.C., 1974. The internal evolution of planetary-sized objects. Moon, 9: 167-182. Walcott, R.I., 1970. Isostatic response to loading of the crust in Canada. Can. J. Earth Sci., 7: 716—727. Wasserburg, G.J. and DePaulo, D.J., 1979. Models of Earth structure inferred from neodymium and strontium isotopic abundances. Proc. NatI. Acad. Sci. U.S.A., 76: 3594—3598. Weertman, J., 1970. The creep strength of the Earth’s mantle. Rev. Geophys. Space Phys., 8: 145—168. Weertman, J. and Weertman, J.R., 1975. High temperature creep of rock and mantle viscosity. Annu. Rev. Earth Planet. Sci., 3: 293—315. Wetherill, G., 1976. The role of large bodies in the formation of the Earth and Moon. Proc. Lunar Sci. Conf., 7: 3245—3257. Wu, P. and Peltier, W.R., 1982a. Viscous gravitational relaxation. Geophys. JR. Astron. Soc., (in press). Wu, P. and Pettier, W.R., 1982b. Glacial isostatic adjustment and the free air gravity anomaly as a constraint on deep mantle viscosity. Geophys. J.R. Astron. Soc., (in press).