Subsolidus convective cooling histories of terrestrial planets

IcAaus 38, 192-211 (1979)

Subsolidus Convective Cooling Histories of Terrestrial Planets G. S C H U B E R T

Department of Earth and Space, Sciences, University of California, Los Angeles, California 90024 P. C A S S E N ANn R. E. Y O U N G

NASA Ames Research Center, Moffett Field, California 9~035 Received August 17, 1978; revised December 19, 1978 The subsolidus convective cooling histories of terrestrial planets evolving from hot initial states are investigated with a simple analytic model which simulates the average heat transport in a vigorously convecting mantle devoid of internal heat sources. The temperature dependence of the effective viscosity of mantle rocks is the single most important factor controlling thermal history. It is responsible for the growth of the rigid lithosphere, a rheological and thermal boundary layer, and serves the function of a thermostat, regulating the rate of cooling by the negative feedback between viscosity and temperature. Except for a relatively short period of time when mantle temperature decreases rapidly during the early stages of cooling, a planet cools mainly by thickening its lithosphere; the underlying mantle temperature decreases relatively slowly. On one-plate planets, the growth of a rigid lithosphere involves an imbalance between the surface heat flux and the heat flow from the mantle; the former is always larger than the latter. Primordial heat can contribute substantially, e.g., as much as about a fourth or a third, to the present surface heat flux of a planet. For both these reasons, the radiogenic heat source content of a planet is likely to be overestimated by inferences from surface heat flow observations. INTRODUCTION The fact that the terrestrial planets appear to be differentiated bodies implies that they have been extensively heated. With the possible exception of Mars, there is evidence t h a t this heating and differentiation occurred early in their histories (Solomon and Chaiken, 1976), perhaps during formation or by the release of gravitational potential energy a c c o m p a n y ing core formation (Birch, 1965; Tozer, 1965). I t is of interest therefore to study the thermal evolutions of terrestrial planets cooling from hot initial states. At some point in the cooling of a terrestrial planet, its mantle would solidify, if it were ever molten or partially molten to begin with,

and subsequent cooling of the planet would be controlled mainly b y subsolidus convection (Runcorn, 1962 ; Tozer, 1967 ; T u r c o t t e and Oxburgh, 1969; Schubert et al., 1969). The primary objective of this paper is to quantitatively investigate this process. Specifically, we address the following questions: If the mantle of a terrestrial planet were hot to begin with, how hot would it be after 4.5 billion years of subsolidus convective cooling without the benefit of internal radiogenic heating or heat input from the core? Is subsolidus convection so efficient at cooling a planetary mantle that primordial heat makes a negligible contribution to the present thermal state? Are the terrestrial planets

192 0019-1035/79/050192-20502.00/O Copyright O 1979 b y Academic Press, Inc. All rights of reproduction in any form reserved.

TERRESTRIAL PLANET COOLING HISTORY in thermal steady state so that measurements of surface heat flux can be used to infer concentrations of internal heat sources ? A closely related question, not specifically addressed here, concerns the possibility that convection could be so efficient as to freeze the core of a terrestrial planet (Schubert and Young, 1976). The only way to shed some light on these issues is by modeling the subsolidus convective cooling process. However, rigorous numerical modeling of highly nonlinear mantle convection is probably not feasible at present. The formidable nature of such an effort argues for the development of simplified approaches to the problem, a point forcefully argued by Tozer (1972a). The very large uncertainties in the properties of the mantles of the terrestrial planets reinforces the requirement for a simplified treatment of convection; the calculations cannot be so time consuming on present day computers as to preclude systematic investigations of parameter variations. With these ideas in mind, we have developed a simple analytic model for the subsolidus convective cooling of a planet. We assume that the average heat flux from a vigorously convecting mantle can be estimated from the temperature of the mantle using a formula which has strong though not rigorous support from theoretical scaling arguments, numerical calculations, and laboratory experiments. The formula connecting heat flux through a convecting layer with temperature can be validated in a number of circumstances where experimental data or detailed numerical results are available. Therefore it probably represents the essential physics of the average behavior of a convecting region, at least qualitatively, even in situations where it has not been tested. The model incorporates a mantle viscosity which is proportional to the exponential of the inverse absolute temperature of the mantle. The temperature dependence of the mantle viscosity is the single most

193

important factor controlling the thermal evolution of the planet (Tozer, 1967, 1972a, 1974). It acts as a thermostat to regulate the mantle temperature. Initially, when the planet is hot, mantle viscosity is low, and extremely vigorous convection rapidly cools the planet. Later in its history, when the planet is relatively cool, its viscosity is higher and more modest convection cools the planet at a reduced rate. Finally, the model also includes a lithosphere which thickens as the planet cools. This essential characteristic of a planet cooling by subsolidus convection is another aspect of the temperature dependence of the rheology. The lithosphere is the near surface thermal boundary layer of a planet (Turcotte and Oxburgh, 1967). Because of the strong temperature dependence of the planet's viscosity, it is also the near surface rheological boundary layer. The main variation of temperature and viscosity in the planet at a given time in its evolution is between the cold, rigid lithosphere and the interior. Except for an initial period of rapid decrease in mantle temperature from a hot initial state, a planet cools mainly by thickening its lithosphere, while the temperature of the underlying convective mantle falls relatively slowly. The growth of a rigid lithosphere leads to a lag between the decaying surface heat flux and the decreasing heat flux into the base of the lithosphere from the underlying cooling mantle; the heat escaping from a planet's surface exceeds the heat escaping from its mantle, a result also found by Cassen et al. (1979) in their recent study of lunar thermal history. The time-dependent heat conduction problem in the lithosphere is simplified by assuming that the temperature in the lithosphere immediately adjusts to the instantaneous steady state profile as the lithosphere thickens. This is a good approximation for the cases considered in this paper, as will be shown later.

194

SCHUBERT, CASSEN, AND YOUNG

We begin with a detailed description of the model followed by justifications of the simplifications employed to simulate convection in the mantle and conduction in the lithosphere. The model is then used to calculate particular cooling histories for Earth, Moon, MaTs, and Mercury. We have not carried out similar calculations for Venus since it so closely resembles Earth in all the relevant parameters except for surface temperature. From these specific thermal evolutions we arrive at a number of broad conclusions regarding the important consequences of the subsolidus convective cooling of the terrestrial planets. DEVELOPMENT OF THE MODEL The essential features of the model are sketched in Fig. 1. A rigid lithosphere of thickness l(t) overlies a convecting mantle of thickness D -- l(t). The total thickness of the lithosphere and convecting mantle is constant with time t. However, the lithosphere itself thickens with time as the planet cools and the convecting mantle thins accordingly. The surface of the planet is kept at the constant temperature T, which is taken to be 0°C in our subsequent calculations for Earth, Moon, and Mars and 100°C for Mercury. The bottom of the lithosphere is defined as the 800°C isotherm. At temperatures less than 800°C

we assume that the viscosity of rocks is effectively infinite on geologic time scales. The temperature distribution in the lithosphere is assumed to be linear throughout the thermal history of the planet. We justify this assumption in a later section. Thus the surface heat flux q.(t) is equal to the conductive heat flux q~(t) transported upward from the base of the lithosphere and both quantities are given by qs(t) =q,(t) = [-k(800°C - T,(°C))/l(t)-],

where k is the thermal conductivity. The surface heat flux decreases with time as the lithosphere thickens. The major temperature variation in a planet is between the lithosphere and the interior. In a small planet, the average temperature in its convecting mantle is relatively uniform, because the adiabatic rise of temperature with depth is small. Thus, for such a planet-(the Moon, for example), it is reasonable to characterize the mantle with a single temperature To(t). In a large planet (Earth and Venus, for example), the adiabatic temperature increase across the mantle can be substantial. Nevertheless, we also characterize the mantles of the large terrestrial planets with a single temperature T¢(t). For the large planets, Tc is no longer interpretable as the isothermal temperature of the mantle since the spherically averaged

qs(O)

•

T

t Convecting Mantle

~-T:800~1:

p, c, k, K, a Tc (o)

(1)

/ qc(t)

1g

A

Tc(t), v=~e Tc

/

1

FIG. 1. Sketch of model showinggeometryand properties of lithosphere and underlyingconvective mantle. The symbols are explainedin the text.

TERRESTRIAL PLANET COOLING HISTORY

195

profile in the lithosphere and a single temperature in the mantle. Such a simplified representation ignores the structure of the thin transition region below the base of the lithosphere in which the temperature increases with depth from 800°C to T¢. Although we do not consider the details = ~ exp(A/T¢), (2) of the convection, we do use a semiempirical relation relating heat flow from the mantle where ~ is a constant and A is an activation q¢(t) to T¢. This relation does not account energy for subsolidus creep. The viscosity for effects of adiabatic compression and of a convecting mantle should be uniform viscous dissipation on convection. Peltier even if a planet is large. That v may indeed (1972), Turcotte et al. (1974), Hewitt be constant throughout the mantles of et al. (1975), and Jarvis and McKenzie even the large planets is indicated by the (1979) have discussed how viscous dissipainference of the Earth's mantle viscosity tion and adiabatic compression influence from glacial rebound data (Cathles, 1975; convection. In a large planet like the Earth, Peltier, 1976) and the theoretical considera- these effects are important, but not tions of Sammis et al. (1977) and O'Connell dominant, since the dissipation number, (1977). Thus it is reasonable to associate a dimensionless parameter measuring the a single viscosity with a convecting mantle importance of the effects, is only of order and, in turn, a single representative unity. temperature with the uniform mantle Other properties of the mantle are viscosity through a rheological law. For density p, specific heat c, thermal conduca small planet, T¢ is the actual mantle tivity k, thermal diffusivity ~ = k / p c , and temperature and (2) is the relevant thermal expansivity a; they are all assumed rheological law since the effect of pressure to be constants. There are no radiogenic on viscosity is negligible. Thus the interpre- heat sources in either the mantle or the tation of v and To is unambiguous. For the lithosphere, and no heat enters the mantle large planet, however, ~ depends on both through its lower boundary. Our interest is pressure and temperature, and while a in how efficiently subsolidus convection single value of v unambiguously character- cools a mantle from an initial high temperaizes the mantle, the single associated ture T~(0) when the mantle is deprived temperature is somewhat artificial. Indeed, of any other source of heat and in how a single temperature cannot really char- rapidly this cooling thickens a lithosphere acterize the mantle viscosity because of from its initial value l (0). the important influence of pressure and A simple energy balance for the convectthe large variation of both thermodynamic ing mantle relates the time rate of change variables through the mantle. Thus for of its internal energy to the heat flowing the large planets one must attribute more physical significance to ~ than to To. For into the lithosphere this reason, our conclusions stress the p c ( D - 1 ) d T c / d t = --q¢. (3) surface heat flow, mantle viscosity, and variation of Tc with time, rather than the In writing (3) we have neglected the complications of spherical geometry, convalue of T~ itself. Thus, independent of the size of the sistent with the overall philosophy of our planet, our model represents the thermal approach. In any case, such geometrical structure in terms of a linear temperature effects would only become important for

mantle temperature must vary at least as strongly as the adiabat. Rather, Tc is that temperature which will yield a mantle kinematic viscosity ~ of the correct magnitude according to a rheological law of the form

196

SCHUBERT, CASSEN, AND YOUNG

very thick lithospheres. An energy balance at the base of the lithosphere provides an equation for the rate of lithosphere thickening. If q~ exceeds ql there is a net flow of heat to the base of the lithosphere. This must be balanced by heating lithosphere material at its base from 800°C to T¢ thereby thinning the lithosphere. If ql exceeds q: there is a net flow of heat away from the base of the lithosphere. This must be balanced by cooling mantle material just below the base of the lithosphere from T¢ to 800°C thereby thickening the lithosphere. As the vigor of convection decreases during a cooling history, there is, in fact, a decrease in the flow of heat q¢ into the base of the lithosphere. The lithosphere attempts to supplement this reduced heat flux by feeding on the internal energy of the mantle. Accordingly, the lithosphere thickens by "freezing" relatively hot mantle material to its base thereby releasing the internal energy of mantle material at a rate which just balances the difference between q~ and q~. This process is analogous to the growth of a solid layer at the top of a cooling liquid. The analogy is not complete, however, since no latent heat is released upon formation of the lithosphere. There is no phase change involved in lithosphere growth. The lithosphere is distinguished from underlying mantle material only b y its contrasting theological behavior. The equation for lithosphere thickening is pc(T¢ -- TOdl/dt

= q~ - - q¢,

(4)

where T1 is the temperature at the base of the lithosphere (T1 = 800°C). There are two essential unknowns in our model of the cooling history of a planet, T ~ ( t ) and l ( t ) . However, the integration of (3) and (4) requires that we know ql and q~ in rerms of T¢ and l. In (1), we have already given the equation for calculating the heat flux through the lithosphere. The

heat flux from the convecting mantle is given by q~ = B : (T¢ - - T~) ] / ( D × l[ag(Tc ~'

-

- - l) T~)(D -- l)3]/

cxp(A/T¢)Racr]l

~,

an expression which is shown below to be equivalent to a well-known result for convective heat transport. All quantities in (5) have been defined except for g, the acceleration of gravity, Ra¢r, the critical value of the Rayleigh number for the onset of convection, and B a constant. The quantity in brackets in (5) is R a / Ra¢r, where the Rayleigh number Ra is Ra = c~g(Tc - - T ~ ) ( D -- l)3/K~.

(6)

The Rayleigh number measures the vigor of convection. For Ra ~ RaCr convection is relatively weak, whereas for Ra >> P~¢, convection is strong. It is common practice in studies of convection to define a dimensionless heat flux, the Nusselt number Nu, as Nu = q ¢ ( D - 1 ) / E k ( T c - - T,)].

(7)

The Nusselt number is the ratio of the heat flux carried by convection and conduction to the heat flux that would be carried by conduction alone if the system were subjected to the same temperature difference across its boundaries. At the onset of convection Nu --- 1. Nu is much greater than unity for a vigorously convecting system. Equation (5) can then be written in the more succinct form Nu = (Ra/Racr)~.

(8)

This power-law relation between Nusselt number and Rayleigh number is known from laboratory experiments and t h e o r e t -

TERRESTRIAL PLANET COOLING HISTORY ical and numerical calculations to characterize the heat flux from a vigorously convecting system (Ra >> Raor) in a variety of circumstances. Tozer (1967, 1972a, 1972b, 1974) has made extensive use of it throughout his papers discussing the importance of subsolidus creep in regulating the temperatures of the terrestrial planets. Sharpe and Peltier (1978, 1979) have strongly advocated its utility for studying the thermal histories of the planets. The power law relation has also been used by Stevenson and Turner (1979) and Kaula (1979) to investigate thermal history models of the Earth and by Cassen et al. (1979) to model the thermal evolution of the Moon. In a later section we will present a more complete and detailed discussion of the evidence supporting this relation. Early in its thermal evolution, when a planet is quite hot, the mantle will indeed be vigorously eonvecting and (5) or (8) may adequately estimate the heat flux. Certainly, as the planet cools and R a - * Ra~r, the power law relation becomes increasingly unreliable. We have noted that the temperature dependence of mantle viscosity is the single most important factor controlling the cooling history of a planet. This dependence enters our model through the Rayleigh number which is inversely proportional to viscosity. This in turn makes the heat flux from the mantle extremely sensitive to the mantle temperature; from (5) q~ ¢¢ exp(--~A/T¢). Finally, the feedback to the mantle temperature occurs through (3) which says that dTo/dt ~ exp(--BA/T¢). Equations (1) and (3)-(5) constitute a set of first-order, coupled, nonlinear, ordinary differential equations for T¢(t) and l(t). The system can be numerically integrated with respect to time given the initial values To(0) and l(0) and the values of the other parameters. When the lithosphere is thin compared with the thickness of the mantle, i.e., when 1 <~ D, the system decouples and the mantle temperature

197

evolves independently of 1. This is the case throughout almost the entire thermal history of the planet; only after the planet has cooled for very long times and the lithosphere has thickened substantially does T~(t) depend on l(t). Thus, throughout most of the cooling history of the planet To(t) decreases at a rate determined essentially by ~(T~) alone, independent of the growth or even of the very existence of the lithosphere. Calculations of T~(t) with l - 0 yield values of To substantially equivalent to those calculated with l(t)~ 0. Nevertheless, thickening of the lithosphere is an important feature of planetary cooling. In fact, for much of a planet's history, lithospheric growth rather than decreasing mantle temperature is the prime manifestation of the cooling. If there were an Earth similar to our own except for the lack of plate tectonics (Venus perhaps?), then the mantle temperature of that planet would be close to the temperature of our own; the major difference between the two planets would be in the thickness of their lithospheres. In the case 1 ~ D, the mantle temperature can be written as the simple quadrature

rrc dT'c exp(BA/T'¢) ~(o)

(T'o-

Tl) 1+~

- --KI-~ I ag I ~ t.

(9)

D ~-3~ ~ Racr!

HEAT FLUX FROM A VIGOROUSLY CONVECTING FLUID LAYER We have proposed calculating the heat flux from the convecting mantle on the basis of a power-law relation between heat flux or Nusselt number and Rayleigh number. Our purpose in this section is to summarize the theoretical and experimental justifications for such a law. Although these are numerous, it should be stated at the outset that there is presently no experimental data available with which to justify

198

SCHUBERT, CASSEN, AND YOUNG

the power law under the circumstances in which we propose to use it. Rather, we view the application of the power law in our situation as the extension of a principle beyond its rigorously defendable region of validity with the expectation that it will at least be qualitatively correct. Most of the evidence supporting the relation Nu = b Ra~ (10) for the heat transport from a fluid layer at high Rayleigh number exists for the case of a layer heated from below which is transporting a steady quantity of heat. The data from laboratory experiments are fit quite well by Eq. (10). Rossby (1969), Chu and Goldstein (1973), and Garon and Goldstein (1973) report the following results for the heat transport through layers of water. Rossby (1969) : Nu = 0.131Ra °.3, 3 X l0 s ~ Ra 3.4 X 104,

(11)

Chu and Goldstein (1973) : Nu = 0.183Ra °'27s, 1.05 × 10 s >__ Ra >__ 2.76 X 105,

(12)

Garon and Goldstein (1973) : Nu = 0.130Ra °'293, 3.3 X 109 >__ Ra >_ 1.3 X 107•

(13)

Rossby (1969) has also measured the heat transport through layers of silicone oil and mercury. The results are: Nu = 0.184Ra °'~81, 3.5 X 10~ >_ Ra 4 X 10a, silicone oil,

(14)

Nu -- 0.147Ra °-~57, 4.5 X 108 >_ Ra 2 X 104, mercury.

(15)

In (11) to (15) the Rayleigh number is based on the temperature difference across the layer. A number of theories of high Rayleigh number convection provide additional support of relation (10) in the case of heating

from below. Priestley (1954) argued that the heat transport which is controlled by conduction through thin thermal boundary layers at high Ra, should be independent of layer thickness; this requires Nu ¢¢ Ra 1/3. Kraichnan's (1962) mixing length theory also gives Nu ~ Ra 1/3, although a different relation is suggested for extremely high Rayleigh number. The boundary layer theory of Turcotte and Oxburgh (1967) yields Nu = 0.167 Ra 1/3. Howard's (1966) theory of the instability of the thermal boundary layer and Long's (1976) scaling and boundary layer arguments also result in Nu ~ Ra ~/3. There are numerous other experimental, theoretical, and numerical studies which support relation (10) for the fluid layer heated from below; Busse's (1978) recent review provides a detailed discussion of many of these. It is clear from the above that relation (10) adequately describes the steady heat transport through a vigorously convecting layer heated from below. It is equally clear, however, that the values of b and ~ depend somewhat on boundary conditions, on the nature of the fluid (such as the Prandtl number which is the ratio of kinematic viscosity to thermal diffusivity), and weakly on the Rayleigh number itself. In other words, at high Ra and over a limited range of Ra, the heat flux can be well represented by relation (10) whether the convection is two- or three-dimensional, steady or oscillatory, etc. However, the precise values of b and ~ will probably be slowly varying functions of Ra. Because the mantles of the planets are not simply fluid layers heated from below, it is important to expand the justification of relation (10) to more general situations, including internally heated convection and transient convection. Although the models of this paper do not include volumetric heating, such an energy source is undoubtedly important for planetary mantles and the approach of this paper could be readily extended to incorporate it.

TERRESTRIAL PLANET COOLING HISTORY First, consider the application of Eq. (10) to the situation in which steady convection occurs in a layer of fluid that is heated uniformly throughout its interior and has an insulated lower boundary. Kulacki and Nagle (1975) and Kulacki and Emara (1977) have carried out experiments to determine the temperature difference across the layer AT' in this circumstance. They used Joule heating of a dilute aqueous electrolyte to produce energy at the rate per unit volume throughout the layer. They present their results in terms of a Nusselt number Nu' and a Rayleigh number Ra' defined as

Q'

Nu' Ra' -

Q'L'/(k,~T'/L'), ag(L')a(Q'(L')'~. Kp

\

2k

(16) (17)

/

Q'L'

L' is the thickness of the layer, so is the heat flux from the top of the layer in steady state, is the conductive heat flux associated with the temperature difference AT' across the layer. Thus the Nusselt number Nu' is the ratio of the actual to the conductive heat flux. is the temperature difference that would exist across the layer if the internally generated heat were removed by conduction; Ra' is defined by analogy to the ordinary Rayleigh number using this characteristic temperature difference. Their steady-state data are well represented by

k,~T'/L'

Q'L'~/2k

Nu' = 0 . 3 9 6 ( R a ' ) °-2~7, 1.89 X 103 < Ra' < 2.17 X 10 TM.

(18)

Equation (18) can be put into the form of (10) as follows. Rewrite Ra' as Ra'

ag(L')3Kv AT'2(~_7_~L Q'L' _) = ½RaNu',

(19)

199

and substitute (19) into (18) to obtain Nu' = 0.396(½RaNu') °.227.

(20)

Solve (20) for Nu' and find Nu' = 0.228Ra °-294,

(21)

in reasonable agreement with Rossby's (1969) data for water (11). Convection in spherical shells with internal heating and insulated lower boundaries has been studied numerically by Young and Schubert (1974), Schubert and Young (1976), and Schubert (1977). These papers presented thermal models of Mars, Earth, and the Moon from which one can calculate the amounts by which the average temperatures of the convecting mantles exceed the base temperatures of the lithospheres. Log-log plots of these excess temperatures AT against Ra/Raor reveal approximate power law dependences of AT on Ra/Racr which are consistent with (10) and imply/3 ~ 0.3. Next, we consider the application of Eq. (10) to transient convection. Kulacki and Nagle (1975) and Kulacki and Emara (1977) also studied the response of the internal temperature of their convection cell (again, volumetrically heated and insulated at the bottom) to step changes in the heating rate. They measured the time required for a system, initially with a uniform temperature, to come to steady convective equilibrium after a step increase in heating from 0 to and the time for the system in equilibrium with to return to a constant temperature after turning off the heating. For A R a ' > 1.3 )< 105, their results were expressed in the form t'ss = C(/~Ra') m, where t'~8 is the time (measured in units of the conduction time for the temperature to come within 2% of its steady-state value, and 5Ra' is the change in Rayleigh number corresponding to the change in heating 5Q'. The constants C and m were found

et al.

z~Q'

(L')2/K)

AQ'

200

SCHUBERT, CASSEN, AND YOUNG

to be essentially independent of whether AQ' was positive (heating) or negative (cooling) : C = 12 and m = - 0 . 2 1 . We can use Eq. (10) together with the appropriate form of the heat equation to derive an expression for t'~ and compare with the experimental result of Kulacki and E m a r a (1977). T h e time rate of change of AT' is given by

With b = 0.131 and 8 = 0.3, (27) gives

t'~s = 11.4(ARa') -°.~3,

in excellent agreement with the experimental result of Kulacki and E m a r a (1977). For the cooling case, Ra' = 0 and the solution of (23) which satisfies Ra(t' = 0) = (2ARa'/b)~/.+a)

d --

dr'

AT' ----- --b RaaAT ' + L'2Q'/k,

(29)

(22) is

where t' is measured in units of L'2/K. This equation can be written in terms of 1 ~ and R a ' as

/2ARa'\-am+o)

X ((Ra~s/Ra)~-

d -- Ra =

-b

R a ~+~ + 2 R a ' .

(23)

dr' For the heating case Ra' = ARa', and the solution of (23) is t' =

(28)

(2A Ra'b~l~)-al(t+~)

X

f

Ra/Rasj

d u l l -- u 1+~,

(24)

.tO

where we have used the fact t h a t the steady-state Rayleigh number Ra~s is related to ARa' b y bRa~2 +~ = 2ARa'. When Ra/Ras~ = 1 - - ~ t' is t'~,, and (24) yields

and

(25) e is 0.02,

t'~ = (2A Ra'bl/~)-~/(I+~)(I + 8) -~

X B
(1) ,0

,

(26)

1+~ where B(l_,)~+a(1/(1 + 8), 0) is the incomplete beta function (Abramowitz and Stegun, 1964). Since ~ = 0.02 is much less than unity, (26) can be shown to have the approximate value t'ss = - (2ARa'blta) -a/(l+a~ (1 + 8) -1 × ln(,(1 + 8)).

(27)

1).

(30)

T h e time at which Ra = 0.02Rass is (again with b -- 0.131 and 8 = 0.3) tss = 30.3(ARa') -°.~.

(31)

T h e exponent is predicted well, b u t the coefficient is too large in this case, suggesting t h a t (10) overestimates the cooling time. Finally, consider a case in which steady convection is driven b y b o t h internal energy sources and an imposed t e m p e r a t u r e differential (the latter maintained b y fixing the temperatures at b o t h upper and lower boundaries). Cassen and Young (1975) analyzed such a system, with spherical geometry, in connection with a study of the thermal history of the lunar interior. T h e y obtained the heat flux through the lower boundary of the convection zone b y numerically solving the Navier-Stokes equations, and plotted this q u a n t i t y as a function of Q for three values of the viscosity and two values of convection zone thickness. Their Figs. 3 and 4 indicate that, in all cases, the heat flux through the b o t t o m b o u n d a r y is approximately a linear function of Q, with negative slope. This is the expected result for steady state, as long as the flux through the top is independent of Q, as is assumed in (10).

TERRESTRIAL PLANET COOLING HISTORY The observations discussed above indicate that (10) gives a good estimate, if not an exact formula, for the heat flux due to convection in a constant viscosity fluid, regardless of the thermal boundary conditions. Furthermore, the experiments of Booker (1976) and Booker and Stengel (1978) with B~nard convection in a variable viscosity fluid show that temperaturedependent viscosity has a relatively minor effect on the N u - R a relation when the Rayleigh number is based on the viscosity corresponding to the mean temperature of the layer. Also, numerical calculations of convection in non-Newtonian fluids (Parmentier et al., 1976; Parmentier, 1978) suggest that the relation holds for even more complicated theologies. Thus we proceed with some confidence to the applications of expression (10) to the thermal balance of planetary interiors.

TEMPERATURE PROFILE IN THE LITHOSPHERE We have assumed the temperature distribution in the lithosphere to be linear throughout the cooling history of the planet. This can be justified by considering a very closely related problem whose analytic solution has been given by Carslaw and Jaeger (1959). The problem is the freezing of a solid layer at the surface of a liquid. Initially, the liquid fills the region x > 0. The surface x = 0 is maintained at temperature Ts, low enough so that the liquid begins to freeze inward from the surface. At time t, the solid-liquid interface is at x = ~(t). Far from the phase change boundary the liquid temperature is Tf. The interface itself is maintained at the solidification temperature Tm throughout the growth process. The temperature is a solution of the one-dimensional, timedependent heat conduction equation. In

201

the solid, the properties are k, K, p, c, while in the liquid, an enhanced thermal conductivity k* and an enhanced thermal diffusivity K*account for convective heat transport. An energy balance at the phase change interface provides the additional equation for determining its motion ~(t) k(OT/Ox)(x

= ~ - - ~) - - k * ( O T / O x )

(32)

X (x = ~ + ~) = L p ( d ~ / d t ) ,

where e is a small quantity and L is the latent heat. Carslaw and Jaeger (1959) show that the solution for the temperature in the solid layer is (T-

Ts)/(Tm-

T~) = erf[(x/~)~']/erf~,

(33)

= 2~(Kt)1/~,

(34)

where

and ~ is a solution of the transcendental equation (Tin - T~) e-~2/erf~ +

~* (Tin k

--

Tf)

e -~2('/x*)

(K*/~)I/~ erfcE~(~/K*)lt~J =

(i~'/2/c)

~.

(35)

The closeness of this problem to the thickening of the planetary lithosphere is clear. Equations (4) and (32) are identical if we identify c(Tc - Tl) in the thickening lithosphere problem with L in the freezing problem. Further, we can make the correspondences T1-- Ts = Tm-- T8 and T~ Tl = T~-- Tin. Finally, we need to recognize that both k * / k and K*/K arc equivalent to the Nusselt number k * / k = ~*/K = Nu.

(36)

Accordingly, the transcendental equation

202

SCHUBERT, CASSEN, AND YOUNG TABLE I PARAMETER VALUES USED IN COOLING HISTORY CALCULATIONS

Planet

D (km)

O*

640 2030 3000 1740

(~

g (cm/s~) T~ (°C) T¢ (t = 0) (°C) 370 375 1000 160

100 0 0 0

1800 2300 3000 1300

for ~ can be rewritten as (T1-

T,) e-6'

(To - TO erf~ Nul/2e--,~2/Nu = 5 r in.

(37)

erfc(5 Nu -in) If $ << 1 then the temperature in the frozen layer, or the lithosphere, is linear, i.e., (T--

Ts)/(T~- Ts) = x/~.

(38)

T h e temperature profile remains linear as the solid layer thickens. If we assume is small compared with unity we can obtain an explicit expression for ~ from (37) and check on the validity of the assumption a posteriori. We find

"~ --Tr~n(T-~ - 2

T~) N u - In"

(39)

This q u a n t i t y is indeed small throughout most of the thermal histories considered in this paper, thus justifying the use of the linear approximation to the temperature profile in the lithosphere. T h e main difference between the freezing problem discussed in this section and the growing lithosphere problem is t h a t Tf is constant in the freezing problem, while Tc is decreasing with time in the planetary cooling problem. This difference is not expected to alter our conclusion regarding the validity of the linear t e m p e r a t u r e

profile in the lithosphere because T¢ is a very slowly varying function of time throughout almost the entire cooling history of a planet, except in the very early phases when the lithosphere is extremely thin in any case. APPLICATION TO EARTH, MARS, MERCURY, AND MOON We have determined cooling histories for E a r t h ~ , Mars c~, M e r c u r y ~, and Moon ( by integrating (1) and (3)-(5). A number of parameter values are given in Table I in addition to which we chose k = 10 -2 c a l / c m s K, a = 3 X 10 -5 K -I, = 10-2 cm2/s, Tl = 800°C, ~ -- 1.65X 10 ~ cm2/s, A = 5.6 X 104 K, /~ = 0.3, and Racr = 1100. T h e values of k, a, and K are typical of geologic materials. T h e rheological parameter values are representative of those currently thought to be relevant to the E a r t h ' s upper mantle. T h e value of A is based on a formula given b y T u r c o t t e and Oxburgh (1972); this, and the chosen value of ~ give a viscosity of 4.8 X 10 ~ cm2/s at a temperature of 1300°C. T h e choice of TI----800°C is consistent with the minimum temperature generally assumed to be needed for mantle rocks to creep on geologic time scales. There is without question great uncertainty in the rheological parameter values. In a future paper we intend to explore the consequences of varying these values but for our present purposes we consider only a single reasonable set of values. T h e value of /~ was discussed in a previous section. T h e value of Racr is in agreement with the O(103) values of critical Rayleigh number found in a variety of situations. Critical Rayleigh number values depend on the size of the convecting mantle (for spherical geometry), the nature of the b o u n d a r y conditions (isothermal vs adiabatic, for example), whether the mantle is heated from below or from within, etc. Racr = 1100 happens to be nearly the

TERRESTRIAL PLANET COOLING HISTORY critical Rayleigh number for the onset of convection in a spherical shell with dimensions of the Earth's mantle and heated from below (see Schubert and Young (1976) for the pertinent boundary conditions and for other values of Ra¢,). We have assumed that the thermodynamic and rheologic parameters of the mantle rocks have One set of values for all the terrestrial planets. The major differences among the planets occur in the parameters listed in Table I. The value of D for a planet is taken to be the depth to the core-mantle boundary. For Earth, the value is well known. The moment of inertia of the Moon is so close to that of a a homogeneous sphere (Sjogren, 1971; Williams et al., 1974; Kaula et al., 1974), that if the Moon has a pure iron core its radius is limited to 300 to 400 km. Thus we used the entire radius of the Moon for its value of D. The observed mass, radius, and moment of inertia of Mars (Reasenberg, 1977) constrain models of its internal density structure. The value of D for Mars assumes that the planet has a pure iron core whose radius is 0.4 times the planetary radius (Kozlovskaya, 1967; Reynolds and, Summers, 1969 ; Binder, 1969 ; Anderson, 1972; Binder and Davis, 1973; Johnston

203

et al., 1974; Johnston and ToksSz, 1977). Observations of the mass and radius of Mercury place similar constraints on the internal density structure of that planet. The value of D used for Mercury is based o n a value of about 0.75 times t,he planetary radius for the size of a pure iron core (Reynolds and Summers, 1969; Siegfried and Solomon, 1974; Gault et al., 1977). The initial values of mantle temperature To(0) are consistent with estimates of the iron melting temperature (Higgins and Kennedy, 1971) at the pressures corresponding to the core-mantle interfaces in the different planets. They are also consistent with numerical calculations of convective mantle temperatures in the Moon (Cassen and Young, 1975; Schubert et al., 1977; Cassen et al., 1979), Mars (Young and Schubert, 1974), and Earth (Schubert and Young, 1976). Values of gravity and surface temperature are well known. Temperature Tc Figure 2 shows the decrease in mantle temperature with time as the planets cool from their hot initial states. The initial temperature of the Moon is so low, how-

3m

i

, o

FIG. 2. Temperature T of the of the terrestrial planets.

20

40

60 TIME (IO7 yr.)

c o n v e c t i n g m a n t l e as

80

ioo

12o

a function of time for model coolinghistories

204

SCHUBERT, CASSEN, AND YOUNG

ever, that it cools very little during the entire course of its 4.5 billion year evolution. At t = 4.5 Gyr, T¢ is still 1204°C for the Moon, a decrease of only 96°C despite the subsolidus convective cooling of its interior. Thus, although solid state convection m a y be occurring within a planet, if the t e m p e r a t u r e is low enough, or if the viscosity is high enough, the temperature of the deep convecting interior can remain substantially the same over geologic time. As we shall see, cooling produces a thick lithosphere on the M o o n ; it reduces the temperature beneath this lithosphere only slightly. T h e situation is quite different for the other planets which begin cooling from much higher temperatures. There is an extremely rapid reduction in mantle temperature very early in the histories of these planets. This is followed b y a much more gradual decrease in temperature over most of the lifetime of the planet. When To is high, is low, and convection is particularly efficient as exemplified by @ and o~, especially during the first few hundred million years of their evolutions. When Tc is low, v is high and convection cools the planet relatively slowly. Cooling is self-

regulated through the dependence of on T¢. Despite the wide range in initial temperatures and planetary radii, temperatures afrer 4.5 G y r are very much the same for @, o~, and (. At t = 4.5 Gyr, Tc is 1204°C for (, 1283°C for c~, and 1307°C for 0 . For ~, convection cools the interior to Ra = Ra¢,, N u = 1 after only about 1.62 G y r at which time To = 1229°C. Since further cooling of ~ would take place by conduction, we end the convective calculation at 1.62 Gyr. Additional calculations for ~ starting from temperatures between 1500 and 3000°C show (Fig. 3) the extreme insensitivity of Tc after 4.5 G y r to the initial temperature. T h e individual curves in Fig. 3 showing T~ vs t up to 1.2 G y r are for T¢(0) = 3000, 2500, 2000, and 1500°C. T h e To vs t curves for all these initial temperatures lie in the shaded area of the inset which extends the time scale to 4.5 Gyr. For T¢(0) = 1500°C, one finds T~(4.5 Gyr) equal to 1275°C, while for an E a r t h model which starts with T~ = 3000°C, T¢(4.5 Gyr) = 1307°C. Convection reduces the initial 1500°C temperature difference between the models to only 32°C after

33OO

27OO Fi,i n.-

rlhl O. i,l I--

TIME

(I07yr.)

FIG. 3. Cooling histories of the Earth's mantle for different initial temperatures. All cooling curves for starting temperatures between 1500 and 3000°C fall in the shaded region of the inset.

TERRESTRIAL PLANET COOLING HISTORY 4.5 billion years. If a planet is very hot, e.g., after both core formation and solidification of the mantle, subsolidus convection cools it immediately to a temperature level determined essentially by the rheology of the mantle.

IP I I

,

,

205 ,

,

,

(

i

Viscosity v

The dramatic increase in viscosity early in the cooling of ~ , o~, and ~ reflects the rapid decrease in T¢, as shown in Fig. 4. v for ( increases only by about one order of magnitude during 4.5 Gyr because of the small decrease in T¢ over geologic I I I I I I time. At t = 4.5 Gyr, v is substantially the 0 8 16 24 32 40 48 same for ~ , ~ , a n d ( @for @ is4.1 X 1021 TiME (lOeyr.) cm2/s, v for c~ is 7 X 1021 cm2/s, and v FIG. 5. Growth of planetary lithospheres. for ( is 4.8 X 1022 cm2/s) reflecting the approximate equality of temperature cal- Lithosphere Thickness l culated for the mantles of these planets. Lithospheres on all the planets grow from the assumed 100-m initial thickness extremely rapidly during the very earliest stages of cooling (Fig. 5). The growth rate is more modest throughout most of the planet's history. The lithospheres thicken to 1 km at t = 104 yr for (, 2.5 X 104 yr for !~, 6 X 104 yr for c~, and 35 Myr for I021 @. Ten-kilometer-thick lithospheres are formed after only 1 Myr for (, 8 Myr for =, i0 ~ 9, 160 Myr for c~, and 280 Myr for (~. >I-,After 4.5 Gyr, 1 is about 225 km for ~ , o 300 km for c~, and 550 km for (. After 10IS iOiS 1.62 Gyr, l is 246 km for 9. We have, of course, assumed that lithospheric growth IOiS on the planets is unimpeded by other w ,(cm2/s) z processes. During the initial growth of the lithosphere when the planet is hot, vigorous 10=4 convection may preclude the formation of a competent lithosphere. A higher rate of impacting objects during this time period TIME (107y~) than at present probably also slows the accumulation of a competent lithosphere. - 0 8 16 24 32 40 48 Thus, while our calculation shows that TIME (ICP yr.) nearly a third of a billion years is required FIG. 4. K i n e m a t i c viscosity r of t h e e o n v e c t i n g to form a 10-km-thick lithosphere on m a n t l e vs time for cooling history models of t h e terrestrial planets. Earth, the time is probably closer to a

5'° I

.y

206

SCHUBERT, CASSEN, AND YOUNG

"~ 10"4~o, % 8 ~-~

g _J

10-I 5xl(}

l

I

I

I

i

8

16

24

32

40

TIME (lOSy~)

FIG. 6. Decrease with time of surface heat flux (solid curves) and heat flux into the base of the lithosphere (dashed curves) for cooling histories of the terrestrial planets. billion years, in agreement with other lines of evidence suggesting a v e r y thin terrestrial lithosphere during its first billion years of evolution (Wetherill, 1972). Except for the E a r t h , lithosphere growth is so rapid during these early cooling stages t h a t the impediments to lithosphere formation probably make little difference to its thickness after 4.5 Gyr. On Earth, constraints on the growth of a lithosphere continue through the present. Plate tectonics severely limits the lithosphere thickness beneath the ocean basins. However, our estimates of lithosphere thickness should be relevant to the lithosphere beneath continental shields. Since c~ and ( show no evidence of plate tectonics our estimates of lithosphere thickness should be directly applicable. W h y we have plate tectonics on E a r t h with the continual creation and destruction of oceanic lithosphere is still an open question. When we are asked to explain why there is no plate tectonics on Mars,

Mercury, or the Moon we are t e m p t e d to say t h a t these planets have thicker lithospheres which are more resistant to breakup in the style of E a r t h tectonics. However, one may wonder if such tectonics occurred on these smaller terrestrial planets when their lithospheres were thinner only to have the evidence obliterated by meteoritic bombardment. We have no evidence from the present geologic surface features of these planets t h a t plate tectonics ever occurred on them, while very ancient cratered surfaces have survived. Kaula (1975) has suggested that all the terrestrial planets did in fact pass through a stage of plate tectonics in their evolutions, b u t t h a t this stage occurred too early and too rapidly on Mercury and the Moon for a n y traces of it to remain.

Surface Heat Flux qs and Heat Flux into the Base of the Lithosphere q~ Figure 6 shows t h a t both surface heat flux (solid curves) and heat flux into the base of the lithosphere (dashed curves) undergo dramatic decreases early in the cooling of all the planets (the assumed initial lunar temperature is so low t h a t q¢ decays only slowly over the entire thermal evolution of the Moon). T h e initial values of qs for (~, c~, and ( are 800 pcal/cm ~ s since we assumed an initial lithosphere thickness of 100 m and a temperature drop of 800°C across the lithosphere. T h e value of q~ (0) for U is slightly lower since we take only a 700°C temperature difference across its lithosphere. T h e initial values of qo are indicated on the figure for ~ , o~, and ~ by the symbols on the vertical axis (qo(0) for ( is clear from the intersection of the qc(t) curve with the axis). This initial difference between qs and qc is simply a consequence of the model and the arbitrarily chosen initial values of l and T~. However, the lithosphere quickly thickens, T¢ rapidly decreases, and qs and q¢ soon come into a self-consistent relationship. After the al-

TERRESTRIAL PLANET COOLING HISTORY most instantaneous drop in q~ and q~, the heat fluxes decay relatively slowly over most of geologic time. A significant result of our calculations is the lag between the decay of the surface heat flux and the decay of the heat flux into the base of the lithosphere. In our model, this is entirely due to the formation and thickening of the lithosphere. At any time in the cooling history of a planet there is more heat flowing through the surface than is flowing into the base of the lithosphere. Even if the convection in the mantle is in a quasi-steady state, the heat flow from the mantle is not in balance with the heat flow through the surface. Thus the use of present-day surface heat flux observations to infer the total concentration of radiogenic heat sources in a planet on the basis of a presumed steady state thermal balance must be viewed with caution. The discussion of this paragraph assumes that the lithosphere is not itself part of the convecting system. This is true for one-plate planets such as the Moon, Mercury, and Mars. For Earth, the tectonic plates are an integral part of the mantle convection system, and it is not clear if the total volume of lithosphere continues to increase throughout geologic time. Daly and Richter (1978) have recently addressed the issue of whether there is a balance between surface heat flux and instantaneous internal heat sources. On the basis of numerical calculations of convection in a plane layer with decaying radiogenic heat sources, they concluded that the surface heat flux exceeds the instantaneous internal heat production even for convection with initial Rayleigh number as large as 10~. They attribute this to the fact that conduction must still play a role in removing heat from the core region of a convecting system with decaying internal heat sources. This source of nonequilibrium between internal heat production and surface heat flux is distinct from the one in our model. The physical mechanism in our

207

model which accounts for the lag between "internal heat production" (heat from the mantle) and surface heat flux is the thickening of a rigid lithosphere, or alternatively, the tendency of the lithosphere to supplement the decaying heat flux from the mantle by feeding on the internal energy of the mantle. Thus there are at least two reasons to suspect that the estimates of present-day internal heat sources from surface heat flux data may be overestimates. The surface heat flux after 4.5 Gyr of cooling without internal heat sources calculated for the Earth is still 0.35 gcal/cm2 s, a significant fraction of the actual present day mean surface heat flux of 1.5 ucal/cm ~ s (Oxburgh and Turcotte, 1978). The Moon's surface heat flux after 4.5 Gyr of convective cooling from a modest initial temperature is calculated to be 0.15 gcal/cm ~ s. Two in situ surface heat flux measurements for the Moon yielded values of about 0.5 and 0.38 geal/cm ~ s (Langseth et al., 1976). Detailed analyses of topographic effects indicate that 0.33 ucal/cm ~ s, rather than 0.38 ucal/cm 2 s might be more representative of the regional heat flux at the Taurus-Littrow site (Langseth et al., 1976). Thus, a rather substantial fraction of the present day heat flow from a planet could be attributed to primordial heat, still another reason to exercise caution in estimating the present day concentration of radiogenic heat sources in a planet from surface heat flux observations. Stevenson and Turner (1979) have also recently concluded that whole mantle convection in the Earth, driven solely by primordial heat content, could be a significant contribution to the present day terrestrial surface heat flux. The convective lunar thermal history models of Cassen et al. (1979) also exhibit a surface heat flux in excess of that attributable to radioactive heat sources. Caution should be exercised in drawing conclusions about the relative contributions of primordial heat and internal heat

208

SCHUBERT, CASSEN, AND YOUNG Ih12

o

m

Z -1u.l ..J

U

~

Ib

~4

TIME

~Z

(lOSyr.)

4U

46

we present two additional figures of the parameters Ra and N u vs time in Figs. 7 and 8. These dimensionless measures of the vigor of convection and the heat flow from the convecting region are so commonly employed in studies of thermal convection that we feel their presentation is warranted. Recall that the onset of convection occurs for Ra = Racr (assumed here to be 1100) when N u = 1. During their early histories, the planets are highly supercritical; even the present day values of Ra for ~ and indicate vigorously convecting mantles. Primordial heat could still drive a modest supercritical convection in the Moon, whereas it could not maintain present-day convection in the mantle of Mercury. T h e present-day estimates of Nu in ~ and c~ indicate t h a t about 10 times as much heat is being transported to the surface by convection as compared with conduction in the mantles of these planets.

Fro. 7. Decay with time of the Rayleigh number in the mantles of the terrestrial planets. production to the present-day surface heat flux because our model does not include internal heat sources. Since cooling histories are nonlinear, one cannot simply add the internally generated heat to the surface heat flux of the present model to determine the surface heat flux of a model containing b o t h primordial and internal sources of heat. We have, however, carried out calculations for models which contain internal heat sources, and we find t h a t primordial heat still contributes to the present-day surface heat flux in amounts similar to t h a t noted above. Therefore there can be significantly fewer radioactives in a planet's interior t h a n implied b y the heat escaping through its surface. Rayleigh Number Ra and Nusselt Number Nu Although the previous figures completely describe the cooling histories of the planets

6xI03~

103

-~ 1112 O0 Z

b.I

Z

8

16

24

32

40

TIME (IOSyr.) Fro. 8. Decrease of the dimensionless heat flux (Nusselt number) as convection in the mantles of the terrestrial planets decays with time.

TERRESTRIAL PLANET COOLING HISTORY CONCLUDING REMARKS We h a v e not a t t e m p t e d to construct t h e r m a l histories of the terrestrial planets with the completeness undertaken, for example, b y ToksSz et al. (1978) whose models also incorporate effects of solid-state convection. Instead, we have explored some simple cooling models which focus on the main questions of how rapidly and in w h a t m a n n e r planets cool b y subsolidus convection. T h e effects of a n u m b e r of obviously i m p o r t a n t processes including internal heating with decaying radiogenic heat sources, heat input to the mantle from a cooling core, partial melting and the u p w a r d differentiation of radioactives, and the possibility of initially stable t e m p e r a t u r e distributions should be incorporated into future calculations. T h e a d v a n t a g e s of the a p p r o a c h used here, i.e., the use of the power-law relation to predict the heat flux from a convecting mantle, in facilitating such studies are clear. These processes can all be incorporated into the equations without increasing the difficulty of their solution a n d the simplicity of the model makes systematic and extensive investigations practical. A continued effort should be made to determine the limitations of the N u - R a power law. Useful progress in this regard can come from l a b o r a t o r y experiments and rigorous numerical calculations. Although we h a v e strongly emphasized in this p a p e r the a d v a n t a g e s and utility of the N u - R a power law in studying subsolidus convection in the planets, the a p p r o a c h should clearly not be used to the exclusion of other numerical a n d experimental studies. E v e n if we cannot come close to b o t h matching all the relevant dimensionless p a r a m e t e r s of p l a n e t a r y mantles and accounting for all physical p h e n o m e n a in our l a b o r a t o r y and numerical investigations, these m u s t u n d o u b t e d l y be carried on to i m p r o v e our basic knowledge of fundamental processes. B y combining all these

209

approaches, we can hope to m a k e meaningful progress in our understanding of mantle convection. ACKNOWLEDGMENTS We thank David A. Yuen for his interest and helpful comments during the course of this work and Gary Murdock for his programming assistance. This research was supported in part by the Planetology Program, Office of Space Science, NASA, under Grant NGR 05-007-317 and by NASA under Grant NSG 7315. REFERENCES ABRAMOWITZ,M., ANDSTEGUN,I. A. (Eds.) (1964). Handbook of Mathematical Functions, p. 263. U. S. Gov't. Print. Off., Washington, D.C. ANDERSON, D. L. (1972). Internal constitution of Mars. J. Geophys. Res. 77, 789-795. BINDER, A. B. (1969). Internal structure of Mars. J. Geophys. Res. 74, 3110-3118. BINDER, A. B., AND DAVIS, D. R. (1973). Internal structure of Mars. Phys. Earth Planet. Int. 7, 477-485. BIRCH, V. (1965). Energetics of core formation. J. Geophys. Res. 70, 6217-6221. BOOKER, J. R. (1976). Thermal convection with strongly temperature-
210

SCHUBERT, CASSEN, AND YOUNG

O'CONNELL, R. J. (1977). On the scale of mantle IX, pp. 213-214. Lunar and Planetary Inst., convection. T eetonophysics 38, 119-136. Houston, Texas. GARON, A. M., AND GOLDSTEIN~ R. J. (1973). OXBUROH, E. R., AND TURCOTTE, D. L. (1978). Mechanisms of continental drift. Rep. Prog. Phys. Velocity and heat transfer measurements in 41, 1249-1312. thermal convection. Phys. Fluids 16, 1818-1825. GAULT, D. E., BURNS, J. A., CASSEN, P., AND PARMENTIER, E. M. (1978). A study of thermal convection in non-Newtonian fluids. J. Fluid STROM, R. G. (1977). Mercury. Ann. Rev. Astron. Mech. 84, 1-11. Astrophys. 15, 97-126. HEWITT, J. M., McKENZIE, D. P., AND WEISS, N. O. PARMENTIER, E. M., TURCOTTE, D. L., AND TORRANCE, K. E. (1976). Studies of finite amplitude (1975). Dissipative heating in convective flows. non-Newtonian thermal convection with applicaJ. Fluid Mech. 68, 721-738. tion to convection in the Earth's mantle. J. HIGGINS, G., AND KENNEDY, G. C. (1971). The Geophys. Res. 81, 1839-1846. adiabatic gradient and the melting point gradient in the core of the earth. J. Geophys. Res. 76, PELTIER, W. a. (1972). Penetrative convection in the planetary mantle. Geophys. Fluid Dyn. 5, 1870-1878. 47-88. HOWARD, L. N. (1966). Convection at high Rayleigh PELTIER, W. R. (1976). Glacial-isostatic adjustment. number. In "Proc. l l t h Cong. Appl. Mech." II. The inverse problem. Geophys. J. Roy. Astron. (H. GSrtler, Ed.), pp. 1109-1115. Springer-Verlag, Soc. 46, 669-705. New York/Berlin. PRIESTLEY, C. H. B. (1954). Convection from a JARVIS, G. T., AND McKENZIE, D. P. (1979). large horizontal surface. AuNt. J. Phys. 7, 176-201. Infinite Prandtl number compressible convection. REASENBERG, R. D. (1977). The moment of inertia J. Fluid Mech., in press. and isostasy of Mars. J. Geophys. Res. 82,369-375. JOHNSTON, D. H., McGETCHIN, J. R., AND TOKS(iZ, REYNOLDS, R. T., AND SUMMERS, A. L. (1969). M. N. (1974). The thermal state and internal Calculations on the composition of the terrestrial structure of Mars. J. Geophys. Res. 79, 3959-3971. planets. J. Geophys. Re*. 74, 2494-2511. JOHNSTON, D. H., AND TOKS(3Z, M. N. (1977). RossBY, H. T. (1969). A study of B6nard convection Internal structure and properties of Mars. with and without rotation. J. Fluid Mech. 36, Icarus 32, 73-84. 309-335. KAULA, W. M. (1975). The seven ages of a planet. RUNCORN, S. K. (1962). Convection in the Moon, Icarus 26, 1-15. Nature 195, 1150-1151. KAULA, W. M. (1979). Thermal evolution of earth SAMMIS, C. G., SMITH, J. C., SCHUBERT, G., AND and moon growing by planetesimal impacts. YUEN, D. A. (1977). Viscosity-depth profile of J. Geophys. Res., in press. the Earth's. mantle: Effects of polymorphie phase KAULA, W. M., SCHUBERT, G., LINGENFELTER, transitions. J. Geophys. Res. 82, 3747-4761. R. E., SJOGREN, W. L., AND WOLLENHAUPT, W. R. (1974). Apollo laser altimetry and inferences as SCHUBERT, G., TURCOTTE, D. L., AND OXBURGH, E. R. (1969). Stability of planetary interiors. to lunar structure. Proc. Lunar Sci. Conf. 5th, Geophys. J. Roy. Astron. Soc. 18, 441-460. 3049-3058. KOZLOVSKAYA,S. V. (1967). Models for the internal SCHUBERT, G., AND YOUNG, R. E. (1976). Cooling the earth by whole mantle subsolidus convection: structure of the Earth, Venus and Mars. Soviet A constraint on the viscosity of the lower mantle. Astronomy--AJ 10, 865-878. Tectonophysics 35, 201-214. KRAICHNAN, R. H. (1962). Mixing-length analysis SCHUBERT, G., YOUNG, R. E., AND CASSEN, P. of turbulent thermal convection at arbitrary (1977). Solid state convection models of the Prandtl numbers. Phys. Fluids 5, 1374-1389. lunar internal temperature. Philos. Trans. Roy. KUL~CKI, F. A., AND EMARA, A. A. (1977). Steady Soc. London A 285, 523-536. and transient thermal convection in a fluid layer with uniform volumetric energy sources. J. Fluid SHARPE, H. N., AND PELTIER, W. R. (1978). Parameterized mantle convection and the earth's Mech. 83, 375-395. thermal history. Geophys. Res. Lett. 5, 737-740. KULACKI, F. A., AND NAGLE, M. E. (1975). Natural SHARPE, H. N., AND PELTIER, W. R. (1979). A convection in horizontal fluid layers with volumethermal history model for the earth with paramtric energy sources. J. Heat Transfer 97, 204-211. eterized convection. Geophys. J. Roy. Astron. Soc., LANGSETH, M. G., KEIHM, S. J., AND PETERS, K. in press. (1976). Revised lunar heat-flow values. Proc. SIEGFRIED, R. W., II, ANn SOLOMON,S. C. (1974). Lunar Sci. Conf. 7th, 3143 3171. Mercury: Internal structure and thermal evoluLONG, R. R. (1976). Relation between Nusselt tion. Icarus 25, 192-205. number and Rayleigh number in turbulent SJOGREN, W. L. (1971). Lunar gravity estimate: thermal convection. J. Fluid Mech. 73, 445-451.

TERRESTRIAL PLANET COOLING HISTORY Independent confirmation. J. Geophys. Res. 75, 7021-7026. SOLOMON, S. C. AND CHAIKEN, J. (1976). Thermal expansion and thermal stress in the Moon and terrestrial planets : Clues to early thermal history. Proc. Lunar ~Jei. Conf. 7th, 3229-3243. STEVENSON, D. J., AND TURNER, J. S. (1979). Fluid models of mantle convection. In The Earth, Its Or~n, Evolution and Structure (M. W. McElhinney, Ed.). Wiley, New York. (In press) ToRs(iz, M. N., HsuI, A. T., ANn JOHNSTON, D. H. (1978). Thermal evolutions of the terrestrial planets. The Moon and the Planets 18, 265-276. TOZER, D. C. (1965). Thermal history of the earth. 1. The formation of the core. Geophys. J. Roy. Astron. Soc. 9, 95-112. TOZER, D. C. (1967). Towards a theory of thermal convection in the mantle. In The Earth's Mantle (T. F. Gaskell, Ed.), pp. 325-353, Academic Press, London/New York. TOZER, D. C. (1972a). The present thermal state of the terrestrial planets. Phys. Earth Planet. Int. 5, 182-197. TOZER, D. C. (1972b). The moon's thermal state and an interpretation of the lunar electrical conductivity distribution. The Moon 5, 90-105.

211

TOZER, D. C. (1974). The internal evolution of planetary-sized objects. The Moon 9, 167-182. TURCOTTE, D. L., AND OXBURGH, E. R. (1967). Finite amplitude convective cells and continental drift. J. Fluid Mech. 28, 29-42. TURCOTTE, D. L. AND OXBURGH, E. R. (1969). Implications of convection within the Moon. Nature 223, 250-251. TURCOTTE, D. L., AND OXBURGH, E. R. (1972). Mantle convection and the new global tectonics. Ann. ReG. Fluid Mech. 4, 33-68. TURCOTTE, D. L., HsuI, A. T., TORRANCE, K. E., AND SCHUBERT, G. (1974). Influence of viscous dissipation on B~nard convection. J. Fluid Mech. 54, 369-374. WETHERILL, G. W. (1972). The beginning of continental evolution. Tectonophysics 13, 31-45. WILLIAMS, J. G., SINCLAIR, W. S., SLADE, M. A., BENDER, P. L., HAUSSER, J. P., MULHOLLAND, J. D., AND SHELUS, P. J. (1974). Lunar moment of inertia constraints from lunar laser ranging (abstract). In Lunar Science V, p. 845. Lunar Sci. Inst., Houston. YOUNG, R. E., AND SCHUBERT, G. (1974). Temperatures inside Mars: Is the core liquid or solid? Geophys. Res. Left. 1, 157-160.