- Email: [email protected]
Stress history of the moon and of terrestrial planets
m~aus 2, 376-395 (1963)
Stress History of the Moon and of Terrestrial Planets ZI)EN]~K KOPAI. Deparlmenl of .4.~tronomy, U~iver.~ il!/ of Jlran.chesler, E~glanr] 1
lieceived October 2S, 1(.163 The aim of the present investigation has t)een to study the effects of nommiform expansion or contraction of elastic globes of planetary size and mass caused by the secular eseapc of their primodial heat and radiogenic heating. Following a statement of the prot)lem contained in Sec. I, the differential equations are set up (See. II) which control the mechanieM effects of isothermal compression (See. I I l) or thermal expansion (See. IV). The coupling of the mechanical and thermal phenomena is formulated exactly; and first approximation solutions of the problem in the form of series expansions are given in Sec. V. These, in turn, are used as a basis of numerical eomputati(ms to reeonstrnct the thermal and stress history (as well as the future) of the Moon, Mercury, and Mars, to the order of approximation to whieh their radiogenic heat sources can be regarded as uniformly dispersed throughout their mass, and the coefficients {ff thermal diffusivity as well as of volume expansion can be treated as constants. It is shown that, within the scheme of these approximations, the secular changes in external radii of these three celestial [)()dies probably did not exceed ().6 _+ 0.1~:c of their present values throughout their long astronomicM past; but whether or not such changes as did take place were monotonic or oscillatory rem'~ins as yet impossit)le to say. I n a d d i t i o n , if (as could h a r d l y h a v e beer~ otherwise) t h e n e w b o r n p l a n e t c o n t a i n e d in its m a s s a certain f r a c t i o n of r a d i o a c t i v e e l e m e n t s (such as found, for instance, in t h e E a r t h ' s crust or in c h o n d r i t i c m e t e o r i t e s ) , t h e radiogenic h e a t released b y t h e i r d e c a y was b o u n d to p r o v i d e a n a d d i t i o n a l source of h e a t - - i n d e p e n d e n t of, a n d s u p p l e m e n t a r y to, t h a t of tile initial c o n d e n s a t i o n - - w h i c h , like t h e l a t t e r , w a s b o u n d t o s p r e a d t h r o u g h t h e e n t i r e m a s s b y c o n d u c t i o n (or r a d i a t i o n ) . Of t h e r a d i o a c t i v e e l e m e n t s which come in q u e s t i o n in t h i s com~ection, such nuclei as those of 1129, P d ~°7, or A F 6, possessing halflives of t h e o r d e r of 106-107 years, w o u l d h a v e released m o s t of t h e i r h e a t d u r i n g t h e first ten million y e a r s of p l a n e t a r y existence ; a n d if t h e y were sufficiently a b u n d a n t , t h e i r a c t i o n could h a v e m e l t e d t h e entire p l a n e t in a few million years. Needless to say, such m o l t e n globes would h a v e been grossly The major part of this investigation was completed while the writer worked as guest investigator c o n v e c t i v e l y u n s t a b l e ; c o n v e c t i o n c u r r e n t s at the Jet Propulsion Laboratory, California Insti- set u p b y t h e s u p e r - a d i a b a t i c t e m p e r a t u r e g r a d i e n t w o u l d h a v e d r a i n e d a w a y t h e excess tute of Technology. 376
L INTRODUCTION W h e n t h e globe of t h e M o o n a n d of o t h e r t e n ~ s t r i a l p l a n e t s a c c u m u l a t e d f r o m the p r i m o r d i a l cloud f r o m which t h e solar syst e m was formed, some 4.5 billion y e a r s ago, it b r o u g h t w i t h it a c e r t a i n a m o u n t of initial h e a t ( p r o b a b l y n o t large e n o u g h t o cause a large-scale escape of t h e volatiles) which b e c a m e i m p r i s o n e d in t h e n e w l y f o r m e d s p h e r e a n d e n d o w e d it w i t h a c e r t a i n initial d i s t r i b u t i o n of t e m p e r a t u r e . I n t h e course of time, this t e m p e r a t u r e profile was b o u n d t o u n d e r g o a secular change as t h e r e s p e c t i v e p l a n e t a r y b o d y cooled b y c o n d u c t i o n a n d e v e n t u a l escape of h e a t (by r a d i a t i o n ) t h r o u g h its surface; t h e r a t e of c o n d u c t i v e cooling ( r o u g h l y p r o p o r t i o n a l to t h e r a t i o of s u r f a c e - t o - v o l u m e of the r e s p e c t i v e configuration) b e i n g t h e higher, t h e s m a l l e r t h e r a d i u s of t h e b o d y .
STRESS HISTORY
OF T H E
MOON AND
OF T E R R E S T R I A L
PLANETS
377
heat so effectively as to enable the globe to of the corresponding radial strains, has not resolidify in the next few million years. so far been attempted by any investigator. In The second family of radioactive elements the present communication such an analysis of importance for planetary cosmology, will be given, and applied to a study of the K 4°, Th -~32, U 235 or U 2~s, possess half-lives of stress history of the Moon, Mercury, and the order of 109-101° years (i.e., a thousand Mars--the latter two terrestrial planets times longer than the A1+-~family, and com- being of sufficiently small mass to permit us parable with the age of the solar system). to regard their interiors as essentially homo= By their combined action the lunar and geneous without too much risk. The principal planetary interiors have been warming up conclusions arrived at in the course of this secularly and are, moreover, likely to still investigation have already been summarized do so for several thousand million years to in the abstract; and in what follows our task will be to substantiate them in detail. come. These two processes--radiative cooling II. E Q U A T I O N S OF MOTION through the surface, and radiogenic heating As is well known, the fundamental equaof the interior--must, moreover, entail also certain mechanical consequences: namely, tion governing the displacement within an if (or to the extent to which) the planetary elastic self-gravitating body caused by a globes can be regarded as solid elastic spheres nonuniform heating can, in vector form, be characterized by specific coefficients of expressed as thermal expansion (or isothermal compres- 2(1 - ~) grad div u - (1 - 2a) curl curl u sion), any secular change in their tempera2 a(1 + a) grad T ture profiles must bring about a correspond3 ing contraction or dilation of their material. 1 -- 2~ 0211 + 2(1 + a)(1 - 2~) If, moreover, any such change is nonuniform E gP + - -/* P ~ ' throughout the interior, its nonuniformity (1) will be bound to invoke radial strains that may exceed the elastic strength of rocks and where u denotes the displacement vector; g, thus affect the structure of the respective the gravity vector; T, the local (scalar) layers (including the surface). temperature; p, the density; and the quantity Following the pioneer work by Urey (1951, 1952, 1955, 1957, 1962), the purely = ~ 7+ ~ / (2) thermal history of the lunar or planetary past has since been studied by a number of stands for the Poisson's ratio characterizing other investigators (MacDonald, 1959; Levin the respective solid (i.e., ratio of transverse and Maeva, 1960; Levin, 1962; Kopal, 1961, compression to longitudinal expansion); 1962; etc.). On the other hand, a history of whereas the thermomechanical effects associated E 9~ (3) with secular heating or cooling of planetary 3 + t~ bodies has so far scarcely been touched. The only paper devoted to a study of the stress denotes its Young's modulus (i.e., reciprocal history of the Moon arising in this comiec- of linear pressure extension). In these tion appears to be that by MacDonald equations, (1960). In this work MacDonald was, however, concerned only with the secular changes in external radius of the Moon caused by different thermal processes men- signifies the coefficient of isothermal comtioned above, or the total release of strain pression (in cm2/dyne); ~, the coefficient of energy arising from this source, while an rigidity (in dyne/cm0; and investig.ation of the nonuniform expansion or contraction in the interior of the lunar or = - (5) p ~p planetary globes, and of proper stratification
(0o)
378
ZDENI~K KOPAL
represents the coefficient of volume thermal expansion (in d e g - 0 ; that of linear thermal expansion being equal to one-third of a. If we wish to study the displacements u caused by secular heating of the body under consideration, the last term on the righthand side of Eq. (1), involving the second time-derivative of u, becomes utterly neg: ligible. In view of this, as well as of the relations (2) and (3) existing between the pairs of the quantities ~, E and t~, fl, our fundamental equation (1) can be rewritten as [1 + (4/3)t~6] grad div u - (g~) curl curl u = a grad T + Bgp. (6) Assuming, moreover, that the body in question finds itself in hydrostatic equilibrium (which, in view of the finite strength of their material is true of all masses of astronomical magnitude), we can set go = -- g r a d P ,
[1 + (4/3)t~3] grad div u - (g/) curl curl u = a grad T - 3 grad P. (8) Let us assume next that our self-gravitating body in hydrostatic equilibrium is spherical (i.e., that deformation arising froin possible axial rotation is ignorable), and that the nommiform heating which causes distortion is likewise characterized by spherically symmetrical temperature distribution T(r). If so, then = flmction of r only:
which entails that curl u = O; and, in consequence, Eqs. (8) reduce to their radial part asserting that d~"
E 0u,. ] (1 -- a) ~"~ = (1 + ¢)(1 - 2a) ~ Dr+
(11)
u~ = 6 r = r~b; and if so, the validity of Eq. (11) requires that, at r = R, I I --a)r(d6/dr)
+
(1 + z ) 6
= O.
(12)
Moreover, the linear nature of Eqs. (9) permits us to separate the effects of thermal expansion and gravity compression of our elastic sphere as follows: I:et the total
(1:3)
be regarded as a sum of t h e "thermal" displacement Cr and " g r a v i t y " displacement $o, respectively. If so, Eqs. (12) permit us to assert that, whereas the "thermal" component ¢7, of tiw total displacement obeys the e~iuatio,, 1+?~g¢3
r~+4
dr,l = a \ ~ } ,
3(1 --o-) \ d r /
the "gravity" component ~bamust satisfy the
2~)(dP~
parallel relation
-
t +:~,~ -
2~7] =0,
where z,.,. represents the purely radial component of the elastic stress tensor (not to be confused with the suffixless Poisson's ratio a on the right-hand side of the above equation!). Consistent with the assumptions previously made, the relevant parts of the displacement vector are
(7)
where P denotes the pressure; and thus replace (6) by
u = 6r
In order that a solution of this equation should represent the nommiform expansion of a sphere heated within, ~b must (a) remain finite at the origin ( r = 0); and (b) satisfy at the outer boundary (r = R) the condition requiring that
3 + 4n[~ \ d r /
3 + 4~B \ d r / = 6r r.
= -~5~r.;
(l,~)
(lO)
the boundary conditions requiring the finitehess of ~b a t the center, a n d the validity of
where we have at)breviated ~l,
r-~.., + ~ r /
(9)
379
STRESS HISTORY OF THE MOON AND OF TERRESTRIAL PLANETS
(12) at the outer boundary r = R, remain c o m m o n to both cases. The former one of these equations will coutrol the radial dilatation caused b y rising temperature inside our elastic sphere, while the latter should permit us to relate the actual state of such a sphere, under the influence of its own gravitation, with the one which would characterize it at zero pressure, hi what follows, the two cases will be discussed in turn. III.
(7 ~< ½, it
follows t h a t ¢ - ~ r / r will be negative throughout the configuration (as it should, since the gravity acts on the entire sphere to compress it); but the corresponding radial strain
d(~r) dr
d dr (r~)
(22)
will vanish inside our configuration at a fractional radius given b y
GRAVITATIONAL COMPRESSION
?"
3 - - O"
R
3(1 + (7)
< 1,
(23)
Suppose t h a t we consider first the limiting case of a self-gravitating elastic sphere, of an arbitrary mass and size, which would be homogeneous at zero pressure. The Eqs. (7) of hydrostatic equilibrium then reduce to
outside of which d(6r)/dr is positive (the radial strain being extension), and negative inside (corresponding to compression). Therefore, if we introduce a nondimensional constant A defined b y the equation
dP/dr = --go = -- (4/3)TrG~'-'r,
2~GF~C~R2 5(3 + 4ufl) - A,
(16)
where 6: denotes the constant of gravitation, and ~, the constant density of the material. Moreover, if so, Eq. (15) assumes the form
it follows t h a t
¢0(0) dr r ~ + 3 ~
3+4u~r'
¢~;(R) 1 (t (r3@G)
1%- (7 A,
(25)
2(1 -- a) A 1%-(7
(.:2T'(J!P2~ N~ ?,2 %_
= \ 3 %_ 4U/~/
c,,
(18)
--a ¥(7
1 A;
(26)
whereas the average compression
and
27rG~2flr"
Cl
C2
(1())
@e - 5(3 -4- 4u~) %- 3 %- -~' r
where ch2 are integration constants. In order to meet the boundary condition (a) and avoid singularity at the origin, we obviously require t h a t c.~ = 0; while the value of ct follows from the boundary condition (b) and Eq. (12), which combined with (18), yields C1
3 - - (7
(17)
and integrates, successively, into r 2 dr
(24)
15 +
( (20)
by (19). In consequence, the desired particular solution of our problem becomes ~ba(r) = -- ,5(3 + 4ufl) As r ~< R and
%_ (7
~, •
(21)
(for all known materials)
4/o
(~r)r2dr = -
%- a
corresponds to t h a t experienced at the fractional radius r / R = g ~ / 3 . If, therefore, o0 denotes the constant (uncompressed) density which the material of our sphere would possess in the absence of gravitation, while ~ stands (as before) for the mean density of the compressed configuration as a whole, the two are obviously connected by the relation = p0[1 -- 3¢a(R %/2/3)]
(28)
while the actual central (pc) and surface (p,) densities of the material in the gravitational field of the respective body should follow from
380
ZDENI~K K()PAL p~ = p0[1 -
3~a(0)l,
(
(2,~)
I V . THERMAL EXPANSION
Let us consider next the astronomically p,~ = po[1 - 3¢~.(te)]. t more important case of secular expansion In order to apply these results to a given of an elastic solid heated from within, more astronomical b o d y - - s u c h as the M o o n - important because the source of internal we must first specify the value of the Pois- heat, and thus the consequent expansion, son's ratio z characteristic of its material. m a y become functions of the time. For it is ()n the Earth, it can be deduced from the well known since the pioneer work by Urey rate of propagation of seismic waves; and (Urey, 1951) t h a t if the chemical composidown to a depth of some 100 km in the tion of the primordial planetary m a t t e r mantle (where pressures occurring inside contains a similar fraction of long-lived the entire h m a r globe are encountered), the radioactive elements (su<'h as K % Th 2a'-', empirical value of a appears to be between W 2a5, or W2as) a s the average sample of 0.26 and 0.27 (of. Birch, 1952)---close enough <'hondritic meteorites, their spontaneous to a round value of }, which we shall here- disintegration will release heat which will after adopt. If so, however, Eqs. (28) and propagate through a solid configuration (2!0 reduce to from the interior to the surface in accordan<'e with well-known e(luations of heat con0,=~(~+2a+..), ( (30) duction. The thermal conductivity of lunar p,=a(l-~+. • .): ( or planetary globes is, however, known to be so low that time intervals of the order of while 109 years are required for major heat trans= p0(1 + 4.6A + • • • ). (31) port to take place across distances of the order of 10a kin, while the half-lives of the In order to proceed further we must evaluate A. For the Moon (of mass 7.35 radioactive elements mentioned above are X 1025 gm and radius R = 1.738 X 10 s cm), of the same order of magnitude (l()+-10 m 3.34 g m / c m a. If we assume that, for years) as the astronomical ages of the bodies average silicate rocks, /3 = 3 X 10 -12 cm in which they are embedded. Therefore, such bodies must still be secularly warming up, sec2/gxn (cf. Birch, 1952) and (for a = 0.25), and their present internal temperatures ~,/3 = 0.6 while (; = 6.68 X 10-s cma/gm between 1000 2000°K still rising (for most see'-', it follows t h a t A = 0.0157. Accordingly, p~ = 3.34 (1 + 0.0314) = 3.44 gm/cma; p.~ recent mmlerical results, of. Kopal, 1962) = 3.34 (1 - 0.0157) = 3.28 g m / c m 3, while for several billion years to come. The appropriate differential equation (>[ o0 = 3.34/(1 + 4.6 X 0.0157) = 3.10 gin/ heat transfer within a spherically symcm a. If, therefore, the force of gravitation were relaxed, the mean density of the lunar metrical configuration, with due regard to globe at zero pressure would be, almost 7~( thermomechanical effects caused by change lower t h a n the one which we actually ob- of volume accompanying thermal expansion serve, the mean radius of the Moon being or contracti(m, was considered in full genincreased by (1.2A)R = 32.7 k m to 1771 km, erality by l)uhamel (1837) more than a and its surface would be 13430 kin'-' larger. century ago, and extended more recently Needless to say, these hypothetical figures by l.owan (1(,)35) to the case of a radioactively heated sphere. In more specific describe only one part of the whole story the effects of gravitational c.ompression. terms, if ~ denotes the coefficient of heat These are superposed on the linear and conduction of planetary matter (hereafter regarded as constant); Cp.v, its specific heat volume changes due to the se<'.ular heating of the interior of the h m a r globe due to a at constant pressure or volume; ej, the gradual disintegration of radioactive ele- anlount of thermal energy (in e r g / g m . s e c ) inents present there; this other aspect of the lit)erated by spontaneous disintegration of a stress history of astronomical bodies of lunar jth radioactive element; and ),j, the half-life or planetary size will be considered in the of its decay, it follows that e '-' The reader may notice the ch;se analogy l)etwe(,n next section.
381
S T R E S S H I S T O R Y OF T H E MOON A N D OF T E R R E S T R I A L P L A N E T S
OT Ot
K 0 r e &" J
•
(32)
sistent with the boundary conditions (a) and (b) of Sec. II, we shall proceed as follows. If the coefficient a of thermal expansion were zero, the homogeneous part d~v dCr r~-re + 4~ r = 0
where K
denotes the fusivity,"
= K/p('V
coefficient of "thermal
(33) dif-
4~; = ~ / C v
(34)
~' = (b/Cv,
(35)
and
the symbols a and ~br being defined, as before, b y Eqs. (5) and (10). The first t e r m on the right-hand side of Eq. (32) safeguardlug the conservation of energy represents the transfer of heat by conduction; the second, the heat input from radioactive sources; and the third, the loss or gain of mechanical energy expended in contraction or expansion (i.e., heat lost b y expansion while the body is being heated, or set free as a result of contraction which accompanies the cooling). In either case, the conservation of mass requires the function ¢Jr(r,t) to obey the differential equation (14) subject to the boundary condition (12). The boundary conditions to be imposed on T(r,t) are, first, t h a t the initial t e m p e r a t u r e distribution at t = 0 be of the form
T(r,O) = fir),
T(R,t) = constant,
of Eq. (14) would admit of a solution of the form ¢,'T = kr m, where k is a constant of integration, and m stands for a root of the indicial equation m(m + 3) = 0. Consequently, the complete primitive of (38) should be of the form ~bT = k ~ + k e r -'~, where kl.e are suitable constants (or functions of time). Their evaluation (by the method of the variation of parameters) reveals t h a t the solution of the full-dress equation (14), subject to the appropriate boundary conditions, can be expressed in the form
C/Tit,t) -- 1 + (4/3)U~
[fo 1 V
T(r,t)r2dr +
T(r,t)redr ],
(39)
where T(r,t) continues to be given by Eq. (32) of this section. Since, however, (39) discloses t h a t ld
r 2 dr (r~bT) -- 1 + (4-/3)g~ T(r,t) -F 12 foR T(r,t)r2dr],
(36)
where fir) is a given function of r; while, on the surface (r = R),
(38)
(40)
Eq. (32) can assume now the more explicit form
(37)
the value of which can [because of the liuearity of Eq. (32)] be taken as zero without the loss of generality. 3 In order to solve Eq. (14) for CT con-
1 -~- 1 + (4-/g)t~B
f ROT
12(~, -- 1)
-
the thermomechanicM term on the right-hand side of (32) and the term P div v which would arise from convective heat transport in fluids. 3 Should we do so, however, we must remember that the solution of our problem then would, in effect, yield the difference T(r,t) - T(R,t) of the inter]ml and surface temperatures, rather than the absolu~ temperature itself.
Ot
reOr
r20r
+
4a~e×p
(--Xfi).
J
(41) In order to proceed further, let us a d o p t (as in the preceding section) the value of the Poisson's ratio to be a = 0.25, corresponding (by Eq. 2) to ~/~ = 0.6. If so, and provided t h a t both K and j can also be
382
z DENI~K K O P A L
regarded as constants, it can be shown (of. Lowan, 19354) t h a t the particular solution of the foregoing equation (41) which satisfies the prescribed boundary conditions ('an be expressed as a sum of
T(r,t) = T~(r,t) + T._,(r,t),
(42)
where the complementary function
rocks, ('au t)e regarded as approximately constant (i.e., arc unlikely to vary by m o r e than a factor 2 or 3 iu humr of planetary interiors); but as the temperature T m a y v a r y by a factor of the order of 100, the differeu('e "y -- l (being proportional to 7') can scarcely t)e regarded as constant. An inspe('tion of the a(:('ompanyiug Table R
Tdrl) = ~
1 -- 2(sin ~ - +
(sin20~,/20~) exp
.
0,,r
.
(43)
while the particular integral ao
T2(r,t) = rA ~
¢i
1 -- 2(sin ~ - +
ff sin 0 , I/o' exp
× exp l
In these equations, the constant A is defined by 3A = R(4 + 5~) mr-',
(45)
while the 0,,'s are positive roots of the transcendental equation 1 -
0 (tot 0 = 3~,02/4(~ ' - 1).
(46)
If *y = t (i.e., Cr = Cv), then A = R, and 0,, = nlr, in which case Eqs. (14) and (32) ('ease to be simultaneous and ¢~.(r,l) can be evaluated by quadratures from (39) by using for T(r,l) the solution of a simple conductiou prot)lem. To what extent is this likely to be true in planetary cases? Quite generally,
('p -- Cv = (T/p")[(1/p) (Op/OT)11F(1/o)[(1/p) (ap/aP)],r = p T ( ~ ~)
(47)
by (4) and (5), so t h a t Cp - ('v Cv
a ~T -
"y - -
(sin20~,/20,,)J
1 =
p-~(~v"
(48)
Both the coefficieuts a and ~ of thermal expansion and isothermal compression, as listed in Table I for several typical silicate 4 It should, however, be noted that Eq. (1) and Eqs. (7-10) of that paper contain several misprints.
l - - h j t ] dr} rdr. (44)
I listing the thermal properties of siIieates reveals that, on the average, 7
-
1 =
2
X
10-~T
(49)
which, for anticipated temperatures between 1000 ° and 2000 °, becomes a quantity of the order of a few per cent, neither large, nor altogether negligible. An introduction of Eq. (49) in Eq. (41) would, of course, render the latter noulinear; and Eqs. (14) with (32) should then be coupled through this nonlinear term. In such a ease, the Duhamell,owau solution of our problem, as represented by Eqs. (43) and (44), would no longer be exact. The fact t h a t the coefficient of the nonlinear cross term in (32) is numeri¢'ally small suggests, however, an approach to the analytical solution of our generalized thernmmechanical problem by successive approximations: uamely, to insert in the uonlinear terin the known functions (43) and (44) of r and l, and to treat the resulting equation as linear but nonhomogeneous. In the simplest case, it m a y be permissible to ignore the difference "7 -- 1 altogether, and to approximate Tl(r,/) as well as T.,(r,t), as given by Eqs. (43) and (44), b y ~etting A = R and 0,~ = mr on their right-hand sides. Iu the following section this will be done, and mmmrical computations performed to evahmte the magnitude of non-
STRESS
H I S T O R Y OF T H E
MOON A N D OF T E R R E S T R I A L
383
PLANETS
TABLE I THERMAL PROPERTIES OF SILICATESa Mineral
10ca (deg-~)
p (gm/ems)
17 15 20 25 40
2.6 2.76 3.2 3.3 2.65
Albite, NaA1Si308 Anorthite, CaA12Si208 Diopside, CaMgSi:O6 Forsterite, Mg_oSiO4 Quartz, SiO2
10-rCv (erg/gm deg)
106(a~/,oBC~)
1013B (croP/dyne)
0.71 0.7 0.69 0.79 0. 698
(deg-Q
19 l1 11 8 27
8 11 17 30 32
After Birch, 1952. uniform thermal expansion (and of strains arising from this source) due to the radioactive heating of the interiors of the Moon, Mercury, and Mars in the course of their long astronomical past. V.
NUMERICAL
The anMytic solution of this boundaryvalue problem for the arbitrary form of f ( r ) was first given by Lowan (1933). If, in particular, this latter function is regarded as expansible in a series of the form f(r) = T + C ( 1
RESULTS
--x0 +"
• • ,
(51)
If the ratio of specific heats in planetary m a t t e r is taken to be unity, the heat conduction problem as represented by the Eqs. (32), (36), and (37) reduces to
where x = - - r / R stands for the fractional distance from the center, and r, r', etc. are arbitrary constants, the solution of (5) can likewise be expressed as
or K O Ot - r 2 Or \
T(r,t)
Or 07] +
ch~ exp (--
=
rTi(x,t)
+ 6r'T'~(x,t) + . . .
~kjl),
+ To(x/),
J T(r,O) = f(r);
T ( R , t ) = 0.
Tl(X,t) = 2
where
(50)
( _ 1)n+1 sin mrx exp n~X
T'~(x,t) = 2
(52)
(--1) n+l ~ - . ~
--K
t
(53) '
exp
(54)
n =2
and T2(x,t) = 2 ~
~bj ~
(--1)'+1 sin n l r ~ x exp ( - h / ) n~-x
j
(55)
-exp[-K(mr/R)et] K(n~c/R) 2 - ~
n=l
Moreover, once the temperature distribution can thus be regarded as known, the fractional thermal expansion ¢T(r,t) can be obtained by quadratures of the expression (39). If, by analogy with (52), we break it down into individual components (for t~ = 0.6) in accordance with the formula ~bv(r,t)
=
(2/9)a[r~b~(x,t) + C ¢ / l ( x , t ) + • • • + ~bo(x,t)],
(56)
an insertion of the foregoing Eqs. (53-55) in (39) and term-by-term integration reveals that ~l(x,t) =
+ 5 ( - - 1 ) n+l n=l
--K
~
t,
(.57)
384
ZDENI~K
KOPAL
~'l(;r,~) : ~ [(~)2 -{- 5( -- ]) n-klSill l~Tra" --
t~mc cos /z~-:r1
] exp
(.~-.r)~
n =i
(n~) ~
,
(58)
and + 5(-1),,+* j
S i l l ?~r./: - -
t~Jrx COS t~.JrX ]
]
II =1
/((,~/~)~[--K(~rr/R)'-'I] - xi
exp ( -- Xi/) -- exp Ill order to proceed with the actual evahmtion of the functions T~.~ or ¢~,2 as given b y the expansions on the right-hand sides of Eqs. (53-55) and of (57 59), it is necessary first to spe('ify a system of units and the values of the numerical constants involved in the respective formulae. Let us hereafter express the time I ill aeons (AE) or 109 years = 3.154 X 10 ~ see; the radial dislance, r, ill terms of 10 a km = 10 s ('m taken as our unit of length; and the temperature T in units of 10 a degrees. If so, the normalized radii, R, of the p l a n e t a r y bodies under consideration become Moon: Mercury :
1.738 2.48
Mars:
3. 313
and the normalized (nondimensional) value of the coefficient of thermal diffusivity R = [3.154 × 10'V(m~)'-'l K = 3.154K is (for material akin to silicate rocks) likely to lie between 0.010 and 0.030, with K = 0.02 representing p r o b a b l y a fair average. F u r t h e r m o ~ , the values of hi and q~/ characterizing the sources of radiogenic heat (according to Urey, 1962) expressed ill terms of our system of units are likewise summarized in Table II. W i t h all these values duly specified, we call proceed with the s u m m a t i o n of the TABLF II X/ (AE t)
Cj (10 a deg,'AE)
j
Nuclide
1
K 40
0. 545
0. 873
2 3 4
!l'h~a2 ['~aa U ~as
0. 049(.t 0. 972 0.1537
0. 0328 (}. 0942 {1.0533
(5.o)
series on the right-hand sides of Eqs. (53-55) or (57-59). For t = 0, these series yield immediately Tl(x,0) = 1, T't(x,())
¢,(x,0) = :3/'2,
(60)
= (1,/,/,6)(1 --x'-'), ¢'~(.r,()) = (1,/,60)(11 -- 5xe);
(6l)
all(l
T2(x,()) = ¢~(x,0) = 0;
(;2)
h u t for 1 > 0 numerical c o m p u t a t i o n s represent the only avenue of approach. As the series ill question converges rather slowly, a s u m m a t i o n of more t h a n one t h o u s a n d terms in each was found necessary to establish the actual values of these sums correctly to 0.1~7c. This has recently been done with the aid of the I B M 7090 C o m p u t e r <)f the Systems Analysis Division of the ,let Propulsion L a b o r a t o r y , California Institute of T e c h nology. The c o m p u t a t i o n s (programmed b y Mr. P. Poulson) have been carried out for four different values of K, 0.010, 0.015, 0.020, and 0.030 for each <.elestial b o d y under <'onsideration; and the results (corre(.t within one unit of the fourth de<'imal place) tabulated for .c = 0(0.1) 0.90 (0.025) 0.975 and t = 0.5 (0.5) 7.0 AE. These tabulations are too extensive to be reproduced fully in this place, but are accessible to the interested reader elsewhere (of. Kopal, 1962). The essential features of numerical results (obtained for the value of K = 0.02) ('all t)e gathered from the graphs reproduced on tile a c c o m p a n y i n g Figs. 1-12. Ill comparing the families of curves representing the se(-ular cooling or radiogenic heating of p l a n e t a r y interiors with the corresponding curves of ¢~.,(x,t), the reader can ascertain for himself the extent to which the former process leads to a
t~
0
0
0
b
J2 (x, f ) m
..,,-
I
•
.I
)1
Ta (x. t ) I
c
r
c
0 ~m
0 0
o D
0
+1 ( ~ t ) ~
~~i~ i~~ ~~i~iilb
©
~0
o
o
~
=-j
m~
z~
STRESS
HISTORY
OF THE
MOON AND
OF T E R R E S T R I A L
xarlR
Fm. T.
PLANETS
389
0
0
o
r t ix, t ) 0
-
"~ " ~ :iii--~
o
~
::~:~':~-'~
~
o
~
.,!i
~
~
~
o
~'2 (x0 t )
~
~
o
o
o
9
o
~-~
o
o
STRESS HISTORY
OF T H E
MOON AND
OF T E R R E S T R I A L
2,0
1.8
L6
1.4
0.$
0.6
0.4
O.!
0.4
OJ~ x mf~
FIG. 10.
0,0
1,0
PLANETS
391
392
z DENI~K KOPAL
x =r//P
FZG. 11.
STRESS
HISTORY
OF T H E
MOON AND
OF TERRESTRIAL
5
0
0.2
0.4
0.6 x =fir
FIG. 12.
| ,0 ~
PLANETS
393
394
ZDENEK KOPAL
shrinkage of the respective body, or the latter brings about its expansion. It should be stressed that although the curves of T.,(x,l) describe the secular cooling of the respective planet from its initial (arbitrary) temperature, assumed to be constant at I = 0, and the sources of radiogenic heat underlying the computed temperatures T~(x,O have likewise t)een assumed to be unifornfly (listributed throughout the whole interior, the secular expansion or contraction due to both these processes is nonuniform and thus bound to give rise to radial slrains proportional to the derivative d
,:~ 0~,)
(@~,
=
r~
+
~,,,.
collected in Table I) for the coefficient a of thermal expansion (and unless the internal temperatm~s of planetary globes are much higher than, say, 2000°K, or unless the relative abundances of radioactive elements in the material constituting planetary globes are such that tile values of Cj as listed in Table II are grossly unrealistic), the secular changes in external radii of the Moon, Mercury, or Mars did not probably exceed 0.6 :t: 0. ler/v of the present values throughout their entire astronomical past. The change in external radius of a planetary globe characterized by constant coefficients a and Cv is proportional to the difference between the rate at which radiogeni(" heat is being produced in its interior, and the rate at which heat flows out through its surface. A number of factors can influence critically this balance--such as the distribution of radioactivity below the surface, the effective heat conductivity, or the initial temperature of the sut)surface layers. Until more becomes known about them than is availat)le at present, the analytical theory developed in this paper ('amlot, unfortunately, help us even to find out whether such changes which the extenml radii of the Moon or of other terrestrial planets may have undergone in the past 4~ billion years have been monotonic or oscillator)'.
It should be noted in particular that these strains are largely confined to layers near the surface and, moreover, their magnitude grows with the size of the body in question being decidedly larger for Mars than for our Moon. The foregoing results refer, to tie sure, to thermomeehanical action of single specific processes, viz. conductive cooling of an initially isothermal elastic sphere, or radiogenic heating due to the action of sources distributed uniformly throughout its interior. The actual resultant variation of ~b.r(r,l) inside the hmar or planetary globes "it any particular time will, of course, depend t{,EFERENCES vitally on the deviations of the planetary ALLXN, D. W. 11956). T h e solution of a special interiors from these idealized conditions, heat and diffusion equalion. Am. JIalt,. i.e., on the actual form of the initial temMonthly 63, 315-323. perature distribution fir) in the newly B1ncn. F. (1952). Elasticity and ('onstitution of formed planetaLv body, as well as the lhe E a r t h ' s interior. J. Geophtl.~. Res. 57, 227distribution of its radioactive heat sources. '~ 286. As long as these characteristics can be only I)UH..~MnL, J. M. C. (1837). "L('s i)henombnes conjectured at, it seems premature to lhermo-m&:aniques," J. 1;:cole Polytech. 15, 1-57. anticipate the detailed variation of the t(oe,u~, Z. 11961). T h e r m a l history of the M o o n and of ll~(' terrestrial planets. Jet Propulsion external radius (or of internal stresses) of Lab., ('alif. l~l.~l. 7'echnol.. Tech. Rept. 32-108, the hmar of planetary globes throughout pp. 1-24. long time intervals. The data exhibited I(OP.~L, Z. (1962). T h e r m a l history of the M o o n earlier in this section can only tie used to and of the terrestrial planets: Numerical reqssert that if a mean value of 2 X l0 '~deg stills. Jel Propld.~io, Lab., CaliJ. In.~t. Technol.. is adopted (in accordau(,e with the data Teeh. lfcpl. 32-225, l)p. 1-108. F(,r a temper'tture distribution inside a e(m(tu(.ting sphere in which the concentration of heat sources increases exponentially t)lltwards, see L(Jwan (1933) ; for tt cltse itl which stleh SOllrees ~r¢, (.onfine([ to a diser(,te outer shell, see Allan (l,q5t~).
[,EVlN, ]~. Y. (1962). T h e r m a l history of fhe Moon, i~ " T h e M o o n , " I.A.I'. S y m p o s i u m No. 14 (Z. i(opal :m(| Z. K. Mikhailov, ells.), p]). 157-167. Ai'adenfic Press, L o n d o n :m(t New York.
STItESS
HISTORY
OF T H E
MOON
LEVI1'% B. Y., AND MAEVA, S. V. (1960). Certain computations of the thermal history of the Moon. Doklady Akad Nauk. USSR 133, 44-47. Low.a.x, A. N. (1933). On the cooling of a radioactive sphere." Phys. Rev. 44, 769-775. LowA~', A. N. (1935). On the cooling of the Earth. Am. J. Math. 57, 174-182. MAcDo.~aLD, G. J. F. (1959). Calculations of the thermal history of the Earth. J. Geophys. Res. 64, 1967-2000. MACDONALD, G. J. F. (1960). Stress history of the Moon. Platter Space Sci. 2, 249-255. UREv, H. C. (1951). The origin and development of the Earth and other terrestrial planets. Geochim. el Co.~mochim. Acta I. 209-277.
AND
OF T E R R E S T R I A L
PLANETS
395
UREY, It. C. (1952). "The Planets," pp. 49-56 and 168--171. Yale Univ. Press, New Haven, Connecticut. UREV, H. C. (1955). The cosmic abundances of potassium, uranium, and thorium, and heat balances of the Earth, the Moon, and mass. Proc. U. S. Natl. Acad. Sci. 41, 127-144. UREY, It. C. (1957). In "Progress in Physics and Chemistry of the Earth," Vol. 2, pp. 46-76. Pergamon Press, London. UREY, H. C. (1962). Origin and history of the Moon. In "Physics and Astronomy of the Moon" (Z. Kopal, ed.), Chap. 13, pp. 481-523. Academic Press, London and New York.
Stress History of the Moon and of Terrestrial Planets ZI)EN]~K KOPAI. Deparlmenl of .4.~tronomy, U~iver.~ il!/ of Jlran.chesler, E~glanr] 1
lieceived October 2S, 1(.163 The aim of the present investigation has t)een to study the effects of nommiform expansion or contraction of elastic globes of planetary size and mass caused by the secular eseapc of their primodial heat and radiogenic heating. Following a statement of the prot)lem contained in Sec. I, the differential equations are set up (See. II) which control the mechanieM effects of isothermal compression (See. I I l) or thermal expansion (See. IV). The coupling of the mechanical and thermal phenomena is formulated exactly; and first approximation solutions of the problem in the form of series expansions are given in Sec. V. These, in turn, are used as a basis of numerical eomputati(ms to reeonstrnct the thermal and stress history (as well as the future) of the Moon, Mercury, and Mars, to the order of approximation to whieh their radiogenic heat sources can be regarded as uniformly dispersed throughout their mass, and the coefficients {ff thermal diffusivity as well as of volume expansion can be treated as constants. It is shown that, within the scheme of these approximations, the secular changes in external radii of these three celestial [)()dies probably did not exceed ().6 _+ 0.1~:c of their present values throughout their long astronomicM past; but whether or not such changes as did take place were monotonic or oscillatory rem'~ins as yet impossit)le to say. I n a d d i t i o n , if (as could h a r d l y h a v e beer~ otherwise) t h e n e w b o r n p l a n e t c o n t a i n e d in its m a s s a certain f r a c t i o n of r a d i o a c t i v e e l e m e n t s (such as found, for instance, in t h e E a r t h ' s crust or in c h o n d r i t i c m e t e o r i t e s ) , t h e radiogenic h e a t released b y t h e i r d e c a y was b o u n d to p r o v i d e a n a d d i t i o n a l source of h e a t - - i n d e p e n d e n t of, a n d s u p p l e m e n t a r y to, t h a t of tile initial c o n d e n s a t i o n - - w h i c h , like t h e l a t t e r , w a s b o u n d t o s p r e a d t h r o u g h t h e e n t i r e m a s s b y c o n d u c t i o n (or r a d i a t i o n ) . Of t h e r a d i o a c t i v e e l e m e n t s which come in q u e s t i o n in t h i s com~ection, such nuclei as those of 1129, P d ~°7, or A F 6, possessing halflives of t h e o r d e r of 106-107 years, w o u l d h a v e released m o s t of t h e i r h e a t d u r i n g t h e first ten million y e a r s of p l a n e t a r y existence ; a n d if t h e y were sufficiently a b u n d a n t , t h e i r a c t i o n could h a v e m e l t e d t h e entire p l a n e t in a few million years. Needless to say, such m o l t e n globes would h a v e been grossly The major part of this investigation was completed while the writer worked as guest investigator c o n v e c t i v e l y u n s t a b l e ; c o n v e c t i o n c u r r e n t s at the Jet Propulsion Laboratory, California Insti- set u p b y t h e s u p e r - a d i a b a t i c t e m p e r a t u r e g r a d i e n t w o u l d h a v e d r a i n e d a w a y t h e excess tute of Technology. 376
L INTRODUCTION W h e n t h e globe of t h e M o o n a n d of o t h e r t e n ~ s t r i a l p l a n e t s a c c u m u l a t e d f r o m the p r i m o r d i a l cloud f r o m which t h e solar syst e m was formed, some 4.5 billion y e a r s ago, it b r o u g h t w i t h it a c e r t a i n a m o u n t of initial h e a t ( p r o b a b l y n o t large e n o u g h t o cause a large-scale escape of t h e volatiles) which b e c a m e i m p r i s o n e d in t h e n e w l y f o r m e d s p h e r e a n d e n d o w e d it w i t h a c e r t a i n initial d i s t r i b u t i o n of t e m p e r a t u r e . I n t h e course of time, this t e m p e r a t u r e profile was b o u n d t o u n d e r g o a secular change as t h e r e s p e c t i v e p l a n e t a r y b o d y cooled b y c o n d u c t i o n a n d e v e n t u a l escape of h e a t (by r a d i a t i o n ) t h r o u g h its surface; t h e r a t e of c o n d u c t i v e cooling ( r o u g h l y p r o p o r t i o n a l to t h e r a t i o of s u r f a c e - t o - v o l u m e of the r e s p e c t i v e configuration) b e i n g t h e higher, t h e s m a l l e r t h e r a d i u s of t h e b o d y .
STRESS HISTORY
OF T H E
MOON AND
OF T E R R E S T R I A L
PLANETS
377
heat so effectively as to enable the globe to of the corresponding radial strains, has not resolidify in the next few million years. so far been attempted by any investigator. In The second family of radioactive elements the present communication such an analysis of importance for planetary cosmology, will be given, and applied to a study of the K 4°, Th -~32, U 235 or U 2~s, possess half-lives of stress history of the Moon, Mercury, and the order of 109-101° years (i.e., a thousand Mars--the latter two terrestrial planets times longer than the A1+-~family, and com- being of sufficiently small mass to permit us parable with the age of the solar system). to regard their interiors as essentially homo= By their combined action the lunar and geneous without too much risk. The principal planetary interiors have been warming up conclusions arrived at in the course of this secularly and are, moreover, likely to still investigation have already been summarized do so for several thousand million years to in the abstract; and in what follows our task will be to substantiate them in detail. come. These two processes--radiative cooling II. E Q U A T I O N S OF MOTION through the surface, and radiogenic heating As is well known, the fundamental equaof the interior--must, moreover, entail also certain mechanical consequences: namely, tion governing the displacement within an if (or to the extent to which) the planetary elastic self-gravitating body caused by a globes can be regarded as solid elastic spheres nonuniform heating can, in vector form, be characterized by specific coefficients of expressed as thermal expansion (or isothermal compres- 2(1 - ~) grad div u - (1 - 2a) curl curl u sion), any secular change in their tempera2 a(1 + a) grad T ture profiles must bring about a correspond3 ing contraction or dilation of their material. 1 -- 2~ 0211 + 2(1 + a)(1 - 2~) If, moreover, any such change is nonuniform E gP + - -/* P ~ ' throughout the interior, its nonuniformity (1) will be bound to invoke radial strains that may exceed the elastic strength of rocks and where u denotes the displacement vector; g, thus affect the structure of the respective the gravity vector; T, the local (scalar) layers (including the surface). temperature; p, the density; and the quantity Following the pioneer work by Urey (1951, 1952, 1955, 1957, 1962), the purely = ~ 7+ ~ / (2) thermal history of the lunar or planetary past has since been studied by a number of stands for the Poisson's ratio characterizing other investigators (MacDonald, 1959; Levin the respective solid (i.e., ratio of transverse and Maeva, 1960; Levin, 1962; Kopal, 1961, compression to longitudinal expansion); 1962; etc.). On the other hand, a history of whereas the thermomechanical effects associated E 9~ (3) with secular heating or cooling of planetary 3 + t~ bodies has so far scarcely been touched. The only paper devoted to a study of the stress denotes its Young's modulus (i.e., reciprocal history of the Moon arising in this comiec- of linear pressure extension). In these tion appears to be that by MacDonald equations, (1960). In this work MacDonald was, however, concerned only with the secular changes in external radius of the Moon caused by different thermal processes men- signifies the coefficient of isothermal comtioned above, or the total release of strain pression (in cm2/dyne); ~, the coefficient of energy arising from this source, while an rigidity (in dyne/cm0; and investig.ation of the nonuniform expansion or contraction in the interior of the lunar or = - (5) p ~p planetary globes, and of proper stratification
(0o)
378
ZDENI~K KOPAL
represents the coefficient of volume thermal expansion (in d e g - 0 ; that of linear thermal expansion being equal to one-third of a. If we wish to study the displacements u caused by secular heating of the body under consideration, the last term on the righthand side of Eq. (1), involving the second time-derivative of u, becomes utterly neg: ligible. In view of this, as well as of the relations (2) and (3) existing between the pairs of the quantities ~, E and t~, fl, our fundamental equation (1) can be rewritten as [1 + (4/3)t~6] grad div u - (g~) curl curl u = a grad T + Bgp. (6) Assuming, moreover, that the body in question finds itself in hydrostatic equilibrium (which, in view of the finite strength of their material is true of all masses of astronomical magnitude), we can set go = -- g r a d P ,
[1 + (4/3)t~3] grad div u - (g/) curl curl u = a grad T - 3 grad P. (8) Let us assume next that our self-gravitating body in hydrostatic equilibrium is spherical (i.e., that deformation arising froin possible axial rotation is ignorable), and that the nommiform heating which causes distortion is likewise characterized by spherically symmetrical temperature distribution T(r). If so, then = flmction of r only:
which entails that curl u = O; and, in consequence, Eqs. (8) reduce to their radial part asserting that d~"
E 0u,. ] (1 -- a) ~"~ = (1 + ¢)(1 - 2a) ~ Dr+
(11)
u~ = 6 r = r~b; and if so, the validity of Eq. (11) requires that, at r = R, I I --a)r(d6/dr)
+
(1 + z ) 6
= O.
(12)
Moreover, the linear nature of Eqs. (9) permits us to separate the effects of thermal expansion and gravity compression of our elastic sphere as follows: I:et the total
(1:3)
be regarded as a sum of t h e "thermal" displacement Cr and " g r a v i t y " displacement $o, respectively. If so, Eqs. (12) permit us to assert that, whereas the "thermal" component ¢7, of tiw total displacement obeys the e~iuatio,, 1+?~g¢3
r~+4
dr,l = a \ ~ } ,
3(1 --o-) \ d r /
the "gravity" component ~bamust satisfy the
2~)(dP~
parallel relation
-
t +:~,~ -
2~7] =0,
where z,.,. represents the purely radial component of the elastic stress tensor (not to be confused with the suffixless Poisson's ratio a on the right-hand side of the above equation!). Consistent with the assumptions previously made, the relevant parts of the displacement vector are
(7)
where P denotes the pressure; and thus replace (6) by
u = 6r
In order that a solution of this equation should represent the nommiform expansion of a sphere heated within, ~b must (a) remain finite at the origin ( r = 0); and (b) satisfy at the outer boundary (r = R) the condition requiring that
3 + 4n[~ \ d r /
3 + 4~B \ d r / = 6r r.
= -~5~r.;
(l,~)
(lO)
the boundary conditions requiring the finitehess of ~b a t the center, a n d the validity of
where we have at)breviated ~l,
r-~.., + ~ r /
(9)
379
STRESS HISTORY OF THE MOON AND OF TERRESTRIAL PLANETS
(12) at the outer boundary r = R, remain c o m m o n to both cases. The former one of these equations will coutrol the radial dilatation caused b y rising temperature inside our elastic sphere, while the latter should permit us to relate the actual state of such a sphere, under the influence of its own gravitation, with the one which would characterize it at zero pressure, hi what follows, the two cases will be discussed in turn. III.
(7 ~< ½, it
follows t h a t ¢ - ~ r / r will be negative throughout the configuration (as it should, since the gravity acts on the entire sphere to compress it); but the corresponding radial strain
d(~r) dr
d dr (r~)
(22)
will vanish inside our configuration at a fractional radius given b y
GRAVITATIONAL COMPRESSION
?"
3 - - O"
R
3(1 + (7)
< 1,
(23)
Suppose t h a t we consider first the limiting case of a self-gravitating elastic sphere, of an arbitrary mass and size, which would be homogeneous at zero pressure. The Eqs. (7) of hydrostatic equilibrium then reduce to
outside of which d(6r)/dr is positive (the radial strain being extension), and negative inside (corresponding to compression). Therefore, if we introduce a nondimensional constant A defined b y the equation
dP/dr = --go = -- (4/3)TrG~'-'r,
2~GF~C~R2 5(3 + 4ufl) - A,
(16)
where 6: denotes the constant of gravitation, and ~, the constant density of the material. Moreover, if so, Eq. (15) assumes the form
it follows t h a t
¢0(0) dr r ~ + 3 ~
3+4u~r'
¢~;(R) 1 (t (r3@G)
1%- (7 A,
(25)
2(1 -- a) A 1%-(7
(.:2T'(J!P2~ N~ ?,2 %_
= \ 3 %_ 4U/~/
c,,
(18)
--a ¥(7
1 A;
(26)
whereas the average compression
and
27rG~2flr"
Cl
C2
(1())
@e - 5(3 -4- 4u~) %- 3 %- -~' r
where ch2 are integration constants. In order to meet the boundary condition (a) and avoid singularity at the origin, we obviously require t h a t c.~ = 0; while the value of ct follows from the boundary condition (b) and Eq. (12), which combined with (18), yields C1
3 - - (7
(17)
and integrates, successively, into r 2 dr
(24)
15 +
( (20)
by (19). In consequence, the desired particular solution of our problem becomes ~ba(r) = -- ,5(3 + 4ufl) As r ~< R and
%_ (7
~, •
(21)
(for all known materials)
4/o
(~r)r2dr = -
%- a
corresponds to t h a t experienced at the fractional radius r / R = g ~ / 3 . If, therefore, o0 denotes the constant (uncompressed) density which the material of our sphere would possess in the absence of gravitation, while ~ stands (as before) for the mean density of the compressed configuration as a whole, the two are obviously connected by the relation = p0[1 -- 3¢a(R %/2/3)]
(28)
while the actual central (pc) and surface (p,) densities of the material in the gravitational field of the respective body should follow from
380
ZDENI~K K()PAL p~ = p0[1 -
3~a(0)l,
(
(2,~)
I V . THERMAL EXPANSION
Let us consider next the astronomically p,~ = po[1 - 3¢~.(te)]. t more important case of secular expansion In order to apply these results to a given of an elastic solid heated from within, more astronomical b o d y - - s u c h as the M o o n - important because the source of internal we must first specify the value of the Pois- heat, and thus the consequent expansion, son's ratio z characteristic of its material. m a y become functions of the time. For it is ()n the Earth, it can be deduced from the well known since the pioneer work by Urey rate of propagation of seismic waves; and (Urey, 1951) t h a t if the chemical composidown to a depth of some 100 km in the tion of the primordial planetary m a t t e r mantle (where pressures occurring inside contains a similar fraction of long-lived the entire h m a r globe are encountered), the radioactive elements (su<'h as K % Th 2a'-', empirical value of a appears to be between W 2a5, or W2as) a s the average sample of 0.26 and 0.27 (of. Birch, 1952)---close enough <'hondritic meteorites, their spontaneous to a round value of }, which we shall here- disintegration will release heat which will after adopt. If so, however, Eqs. (28) and propagate through a solid configuration (2!0 reduce to from the interior to the surface in accordan<'e with well-known e(luations of heat con0,=~(~+2a+..), ( (30) duction. The thermal conductivity of lunar p,=a(l-~+. • .): ( or planetary globes is, however, known to be so low that time intervals of the order of while 109 years are required for major heat trans= p0(1 + 4.6A + • • • ). (31) port to take place across distances of the order of 10a kin, while the half-lives of the In order to proceed further we must evaluate A. For the Moon (of mass 7.35 radioactive elements mentioned above are X 1025 gm and radius R = 1.738 X 10 s cm), of the same order of magnitude (l()+-10 m 3.34 g m / c m a. If we assume that, for years) as the astronomical ages of the bodies average silicate rocks, /3 = 3 X 10 -12 cm in which they are embedded. Therefore, such bodies must still be secularly warming up, sec2/gxn (cf. Birch, 1952) and (for a = 0.25), and their present internal temperatures ~,/3 = 0.6 while (; = 6.68 X 10-s cma/gm between 1000 2000°K still rising (for most see'-', it follows t h a t A = 0.0157. Accordingly, p~ = 3.34 (1 + 0.0314) = 3.44 gm/cma; p.~ recent mmlerical results, of. Kopal, 1962) = 3.34 (1 - 0.0157) = 3.28 g m / c m 3, while for several billion years to come. The appropriate differential equation (>[ o0 = 3.34/(1 + 4.6 X 0.0157) = 3.10 gin/ heat transfer within a spherically symcm a. If, therefore, the force of gravitation were relaxed, the mean density of the lunar metrical configuration, with due regard to globe at zero pressure would be, almost 7~( thermomechanical effects caused by change lower t h a n the one which we actually ob- of volume accompanying thermal expansion serve, the mean radius of the Moon being or contracti(m, was considered in full genincreased by (1.2A)R = 32.7 k m to 1771 km, erality by l)uhamel (1837) more than a and its surface would be 13430 kin'-' larger. century ago, and extended more recently Needless to say, these hypothetical figures by l.owan (1(,)35) to the case of a radioactively heated sphere. In more specific describe only one part of the whole story the effects of gravitational c.ompression. terms, if ~ denotes the coefficient of heat These are superposed on the linear and conduction of planetary matter (hereafter regarded as constant); Cp.v, its specific heat volume changes due to the se<'.ular heating of the interior of the h m a r globe due to a at constant pressure or volume; ej, the gradual disintegration of radioactive ele- anlount of thermal energy (in e r g / g m . s e c ) inents present there; this other aspect of the lit)erated by spontaneous disintegration of a stress history of astronomical bodies of lunar jth radioactive element; and ),j, the half-life or planetary size will be considered in the of its decay, it follows that e '-' The reader may notice the ch;se analogy l)etwe(,n next section.
381
S T R E S S H I S T O R Y OF T H E MOON A N D OF T E R R E S T R I A L P L A N E T S
OT Ot
K 0 r e &" J
•
(32)
sistent with the boundary conditions (a) and (b) of Sec. II, we shall proceed as follows. If the coefficient a of thermal expansion were zero, the homogeneous part d~v dCr r~-re + 4~ r = 0
where K
denotes the fusivity,"
= K/p('V
coefficient of "thermal
(33) dif-
4~; = ~ / C v
(34)
~' = (b/Cv,
(35)
and
the symbols a and ~br being defined, as before, b y Eqs. (5) and (10). The first t e r m on the right-hand side of Eq. (32) safeguardlug the conservation of energy represents the transfer of heat by conduction; the second, the heat input from radioactive sources; and the third, the loss or gain of mechanical energy expended in contraction or expansion (i.e., heat lost b y expansion while the body is being heated, or set free as a result of contraction which accompanies the cooling). In either case, the conservation of mass requires the function ¢Jr(r,t) to obey the differential equation (14) subject to the boundary condition (12). The boundary conditions to be imposed on T(r,t) are, first, t h a t the initial t e m p e r a t u r e distribution at t = 0 be of the form
T(r,O) = fir),
T(R,t) = constant,
of Eq. (14) would admit of a solution of the form ¢,'T = kr m, where k is a constant of integration, and m stands for a root of the indicial equation m(m + 3) = 0. Consequently, the complete primitive of (38) should be of the form ~bT = k ~ + k e r -'~, where kl.e are suitable constants (or functions of time). Their evaluation (by the method of the variation of parameters) reveals t h a t the solution of the full-dress equation (14), subject to the appropriate boundary conditions, can be expressed in the form
C/Tit,t) -- 1 + (4/3)U~
[fo 1 V
T(r,t)r2dr +
T(r,t)redr ],
(39)
where T(r,t) continues to be given by Eq. (32) of this section. Since, however, (39) discloses t h a t ld
r 2 dr (r~bT) -- 1 + (4-/3)g~ T(r,t) -F 12 foR T(r,t)r2dr],
(36)
where fir) is a given function of r; while, on the surface (r = R),
(38)
(40)
Eq. (32) can assume now the more explicit form
(37)
the value of which can [because of the liuearity of Eq. (32)] be taken as zero without the loss of generality. 3 In order to solve Eq. (14) for CT con-
1 -~- 1 + (4-/g)t~B
f ROT
12(~, -- 1)
-
the thermomechanicM term on the right-hand side of (32) and the term P div v which would arise from convective heat transport in fluids. 3 Should we do so, however, we must remember that the solution of our problem then would, in effect, yield the difference T(r,t) - T(R,t) of the inter]ml and surface temperatures, rather than the absolu~ temperature itself.
Ot
reOr
r20r
+
4a~e×p
(--Xfi).
J
(41) In order to proceed further, let us a d o p t (as in the preceding section) the value of the Poisson's ratio to be a = 0.25, corresponding (by Eq. 2) to ~/~ = 0.6. If so, and provided t h a t both K and j can also be
382
z DENI~K K O P A L
regarded as constants, it can be shown (of. Lowan, 19354) t h a t the particular solution of the foregoing equation (41) which satisfies the prescribed boundary conditions ('an be expressed as a sum of
T(r,t) = T~(r,t) + T._,(r,t),
(42)
where the complementary function
rocks, ('au t)e regarded as approximately constant (i.e., arc unlikely to vary by m o r e than a factor 2 or 3 iu humr of planetary interiors); but as the temperature T m a y v a r y by a factor of the order of 100, the differeu('e "y -- l (being proportional to 7') can scarcely t)e regarded as constant. An inspe('tion of the a(:('ompanyiug Table R
Tdrl) = ~
1 -- 2(sin ~ - +
(sin20~,/20~) exp
.
0,,r
.
(43)
while the particular integral ao
T2(r,t) = rA ~
¢i
1 -- 2(sin ~ - +
ff sin 0 , I/o' exp
× exp l
In these equations, the constant A is defined by 3A = R(4 + 5~) mr-',
(45)
while the 0,,'s are positive roots of the transcendental equation 1 -
0 (tot 0 = 3~,02/4(~ ' - 1).
(46)
If *y = t (i.e., Cr = Cv), then A = R, and 0,, = nlr, in which case Eqs. (14) and (32) ('ease to be simultaneous and ¢~.(r,l) can be evaluated by quadratures from (39) by using for T(r,l) the solution of a simple conductiou prot)lem. To what extent is this likely to be true in planetary cases? Quite generally,
('p -- Cv = (T/p")[(1/p) (Op/OT)11F(1/o)[(1/p) (ap/aP)],r = p T ( ~ ~)
(47)
by (4) and (5), so t h a t Cp - ('v Cv
a ~T -
"y - -
(sin20~,/20,,)J
1 =
p-~(~v"
(48)
Both the coefficieuts a and ~ of thermal expansion and isothermal compression, as listed in Table I for several typical silicate 4 It should, however, be noted that Eq. (1) and Eqs. (7-10) of that paper contain several misprints.
l - - h j t ] dr} rdr. (44)
I listing the thermal properties of siIieates reveals that, on the average, 7
-
1 =
2
X
10-~T
(49)
which, for anticipated temperatures between 1000 ° and 2000 °, becomes a quantity of the order of a few per cent, neither large, nor altogether negligible. An introduction of Eq. (49) in Eq. (41) would, of course, render the latter noulinear; and Eqs. (14) with (32) should then be coupled through this nonlinear term. In such a ease, the Duhamell,owau solution of our problem, as represented by Eqs. (43) and (44), would no longer be exact. The fact t h a t the coefficient of the nonlinear cross term in (32) is numeri¢'ally small suggests, however, an approach to the analytical solution of our generalized thernmmechanical problem by successive approximations: uamely, to insert in the uonlinear terin the known functions (43) and (44) of r and l, and to treat the resulting equation as linear but nonhomogeneous. In the simplest case, it m a y be permissible to ignore the difference "7 -- 1 altogether, and to approximate Tl(r,/) as well as T.,(r,t), as given by Eqs. (43) and (44), b y ~etting A = R and 0,~ = mr on their right-hand sides. Iu the following section this will be done, and mmmrical computations performed to evahmte the magnitude of non-
STRESS
H I S T O R Y OF T H E
MOON A N D OF T E R R E S T R I A L
383
PLANETS
TABLE I THERMAL PROPERTIES OF SILICATESa Mineral
10ca (deg-~)
p (gm/ems)
17 15 20 25 40
2.6 2.76 3.2 3.3 2.65
Albite, NaA1Si308 Anorthite, CaA12Si208 Diopside, CaMgSi:O6 Forsterite, Mg_oSiO4 Quartz, SiO2
10-rCv (erg/gm deg)
106(a~/,oBC~)
1013B (croP/dyne)
0.71 0.7 0.69 0.79 0. 698
(deg-Q
19 l1 11 8 27
8 11 17 30 32
After Birch, 1952. uniform thermal expansion (and of strains arising from this source) due to the radioactive heating of the interiors of the Moon, Mercury, and Mars in the course of their long astronomical past. V.
NUMERICAL
The anMytic solution of this boundaryvalue problem for the arbitrary form of f ( r ) was first given by Lowan (1933). If, in particular, this latter function is regarded as expansible in a series of the form f(r) = T + C ( 1
RESULTS
--x0 +"
• • ,
(51)
If the ratio of specific heats in planetary m a t t e r is taken to be unity, the heat conduction problem as represented by the Eqs. (32), (36), and (37) reduces to
where x = - - r / R stands for the fractional distance from the center, and r, r', etc. are arbitrary constants, the solution of (5) can likewise be expressed as
or K O Ot - r 2 Or \
T(r,t)
Or 07] +
ch~ exp (--
=
rTi(x,t)
+ 6r'T'~(x,t) + . . .
~kjl),
+ To(x/),
J T(r,O) = f(r);
T ( R , t ) = 0.
Tl(X,t) = 2
where
(50)
( _ 1)n+1 sin mrx exp n~X
T'~(x,t) = 2
(52)
(--1) n+l ~ - . ~
--K
t
(53) '
exp
(54)
n =2
and T2(x,t) = 2 ~
~bj ~
(--1)'+1 sin n l r ~ x exp ( - h / ) n~-x
j
(55)
-exp[-K(mr/R)et] K(n~c/R) 2 - ~
n=l
Moreover, once the temperature distribution can thus be regarded as known, the fractional thermal expansion ¢T(r,t) can be obtained by quadratures of the expression (39). If, by analogy with (52), we break it down into individual components (for t~ = 0.6) in accordance with the formula ~bv(r,t)
=
(2/9)a[r~b~(x,t) + C ¢ / l ( x , t ) + • • • + ~bo(x,t)],
(56)
an insertion of the foregoing Eqs. (53-55) in (39) and term-by-term integration reveals that ~l(x,t) =
+ 5 ( - - 1 ) n+l n=l
--K
~
t,
(.57)
384
ZDENI~K
KOPAL
~'l(;r,~) : ~ [(~)2 -{- 5( -- ]) n-klSill l~Tra" --
t~mc cos /z~-:r1
] exp
(.~-.r)~
n =i
(n~) ~
,
(58)
and + 5(-1),,+* j
S i l l ?~r./: - -
t~Jrx COS t~.JrX ]
]
II =1
/((,~/~)~[--K(~rr/R)'-'I] - xi
exp ( -- Xi/) -- exp Ill order to proceed with the actual evahmtion of the functions T~.~ or ¢~,2 as given b y the expansions on the right-hand sides of Eqs. (53-55) and of (57 59), it is necessary first to spe('ify a system of units and the values of the numerical constants involved in the respective formulae. Let us hereafter express the time I ill aeons (AE) or 109 years = 3.154 X 10 ~ see; the radial dislance, r, ill terms of 10 a km = 10 s ('m taken as our unit of length; and the temperature T in units of 10 a degrees. If so, the normalized radii, R, of the p l a n e t a r y bodies under consideration become Moon: Mercury :
1.738 2.48
Mars:
3. 313
and the normalized (nondimensional) value of the coefficient of thermal diffusivity R = [3.154 × 10'V(m~)'-'l K = 3.154K is (for material akin to silicate rocks) likely to lie between 0.010 and 0.030, with K = 0.02 representing p r o b a b l y a fair average. F u r t h e r m o ~ , the values of hi and q~/ characterizing the sources of radiogenic heat (according to Urey, 1962) expressed ill terms of our system of units are likewise summarized in Table II. W i t h all these values duly specified, we call proceed with the s u m m a t i o n of the TABLF II X/ (AE t)
Cj (10 a deg,'AE)
j
Nuclide
1
K 40
0. 545
0. 873
2 3 4
!l'h~a2 ['~aa U ~as
0. 049(.t 0. 972 0.1537
0. 0328 (}. 0942 {1.0533
(5.o)
series on the right-hand sides of Eqs. (53-55) or (57-59). For t = 0, these series yield immediately Tl(x,0) = 1, T't(x,())
¢,(x,0) = :3/'2,
(60)
= (1,/,/,6)(1 --x'-'), ¢'~(.r,()) = (1,/,60)(11 -- 5xe);
(6l)
all(l
T2(x,()) = ¢~(x,0) = 0;
(;2)
h u t for 1 > 0 numerical c o m p u t a t i o n s represent the only avenue of approach. As the series ill question converges rather slowly, a s u m m a t i o n of more t h a n one t h o u s a n d terms in each was found necessary to establish the actual values of these sums correctly to 0.1~7c. This has recently been done with the aid of the I B M 7090 C o m p u t e r <)f the Systems Analysis Division of the ,let Propulsion L a b o r a t o r y , California Institute of T e c h nology. The c o m p u t a t i o n s (programmed b y Mr. P. Poulson) have been carried out for four different values of K, 0.010, 0.015, 0.020, and 0.030 for each <.elestial b o d y under <'onsideration; and the results (corre(.t within one unit of the fourth de<'imal place) tabulated for .c = 0(0.1) 0.90 (0.025) 0.975 and t = 0.5 (0.5) 7.0 AE. These tabulations are too extensive to be reproduced fully in this place, but are accessible to the interested reader elsewhere (of. Kopal, 1962). The essential features of numerical results (obtained for the value of K = 0.02) ('all t)e gathered from the graphs reproduced on tile a c c o m p a n y i n g Figs. 1-12. Ill comparing the families of curves representing the se(-ular cooling or radiogenic heating of p l a n e t a r y interiors with the corresponding curves of ¢~.,(x,t), the reader can ascertain for himself the extent to which the former process leads to a
t~
0
0
0
b
J2 (x, f ) m
..,,-
I
•
.I
)1
Ta (x. t ) I
c
r
c
0 ~m
0 0
o D
0
+1 ( ~ t ) ~
~~i~ i~~ ~~i~iilb
©
~0
o
o
~
=-j
m~
z~
STRESS
HISTORY
OF THE
MOON AND
OF T E R R E S T R I A L
xarlR
Fm. T.
PLANETS
389
0
0
o
r t ix, t ) 0
-
"~ " ~ :iii--~
o
~
::~:~':~-'~
~
o
~
.,!i
~
~
~
o
~'2 (x0 t )
~
~
o
o
o
9
o
~-~
o
o
STRESS HISTORY
OF T H E
MOON AND
OF T E R R E S T R I A L
2,0
1.8
L6
1.4
0.$
0.6
0.4
O.!
0.4
OJ~ x mf~
FIG. 10.
0,0
1,0
PLANETS
391
392
z DENI~K KOPAL
x =r//P
FZG. 11.
STRESS
HISTORY
OF T H E
MOON AND
OF TERRESTRIAL
5
0
0.2
0.4
0.6 x =fir
FIG. 12.
| ,0 ~
PLANETS
393
394
ZDENEK KOPAL
shrinkage of the respective body, or the latter brings about its expansion. It should be stressed that although the curves of T.,(x,l) describe the secular cooling of the respective planet from its initial (arbitrary) temperature, assumed to be constant at I = 0, and the sources of radiogenic heat underlying the computed temperatures T~(x,O have likewise t)een assumed to be unifornfly (listributed throughout the whole interior, the secular expansion or contraction due to both these processes is nonuniform and thus bound to give rise to radial slrains proportional to the derivative d
,:~ 0~,)
(@~,
=
r~
+
~,,,.
collected in Table I) for the coefficient a of thermal expansion (and unless the internal temperatm~s of planetary globes are much higher than, say, 2000°K, or unless the relative abundances of radioactive elements in the material constituting planetary globes are such that tile values of Cj as listed in Table II are grossly unrealistic), the secular changes in external radii of the Moon, Mercury, or Mars did not probably exceed 0.6 :t: 0. ler/v of the present values throughout their entire astronomical past. The change in external radius of a planetary globe characterized by constant coefficients a and Cv is proportional to the difference between the rate at which radiogeni(" heat is being produced in its interior, and the rate at which heat flows out through its surface. A number of factors can influence critically this balance--such as the distribution of radioactivity below the surface, the effective heat conductivity, or the initial temperature of the sut)surface layers. Until more becomes known about them than is availat)le at present, the analytical theory developed in this paper ('amlot, unfortunately, help us even to find out whether such changes which the extenml radii of the Moon or of other terrestrial planets may have undergone in the past 4~ billion years have been monotonic or oscillator)'.
It should be noted in particular that these strains are largely confined to layers near the surface and, moreover, their magnitude grows with the size of the body in question being decidedly larger for Mars than for our Moon. The foregoing results refer, to tie sure, to thermomeehanical action of single specific processes, viz. conductive cooling of an initially isothermal elastic sphere, or radiogenic heating due to the action of sources distributed uniformly throughout its interior. The actual resultant variation of ~b.r(r,l) inside the hmar or planetary globes "it any particular time will, of course, depend t{,EFERENCES vitally on the deviations of the planetary ALLXN, D. W. 11956). T h e solution of a special interiors from these idealized conditions, heat and diffusion equalion. Am. JIalt,. i.e., on the actual form of the initial temMonthly 63, 315-323. perature distribution fir) in the newly B1ncn. F. (1952). Elasticity and ('onstitution of formed planetaLv body, as well as the lhe E a r t h ' s interior. J. Geophtl.~. Res. 57, 227distribution of its radioactive heat sources. '~ 286. As long as these characteristics can be only I)UH..~MnL, J. M. C. (1837). "L('s i)henombnes conjectured at, it seems premature to lhermo-m&:aniques," J. 1;:cole Polytech. 15, 1-57. anticipate the detailed variation of the t(oe,u~, Z. 11961). T h e r m a l history of the M o o n and of ll~(' terrestrial planets. Jet Propulsion external radius (or of internal stresses) of Lab., ('alif. l~l.~l. 7'echnol.. Tech. Rept. 32-108, the hmar of planetary globes throughout pp. 1-24. long time intervals. The data exhibited I(OP.~L, Z. (1962). T h e r m a l history of the M o o n earlier in this section can only tie used to and of the terrestrial planets: Numerical reqssert that if a mean value of 2 X l0 '~deg stills. Jel Propld.~io, Lab., CaliJ. In.~t. Technol.. is adopted (in accordau(,e with the data Teeh. lfcpl. 32-225, l)p. 1-108. F(,r a temper'tture distribution inside a e(m(tu(.ting sphere in which the concentration of heat sources increases exponentially t)lltwards, see L(Jwan (1933) ; for tt cltse itl which stleh SOllrees ~r¢, (.onfine([ to a diser(,te outer shell, see Allan (l,q5t~).
[,EVlN, ]~. Y. (1962). T h e r m a l history of fhe Moon, i~ " T h e M o o n , " I.A.I'. S y m p o s i u m No. 14 (Z. i(opal :m(| Z. K. Mikhailov, ells.), p]). 157-167. Ai'adenfic Press, L o n d o n :m(t New York.
STItESS
HISTORY
OF T H E
MOON
LEVI1'% B. Y., AND MAEVA, S. V. (1960). Certain computations of the thermal history of the Moon. Doklady Akad Nauk. USSR 133, 44-47. Low.a.x, A. N. (1933). On the cooling of a radioactive sphere." Phys. Rev. 44, 769-775. LowA~', A. N. (1935). On the cooling of the Earth. Am. J. Math. 57, 174-182. MAcDo.~aLD, G. J. F. (1959). Calculations of the thermal history of the Earth. J. Geophys. Res. 64, 1967-2000. MACDONALD, G. J. F. (1960). Stress history of the Moon. Platter Space Sci. 2, 249-255. UREv, H. C. (1951). The origin and development of the Earth and other terrestrial planets. Geochim. el Co.~mochim. Acta I. 209-277.
AND
OF T E R R E S T R I A L
PLANETS
395
UREY, It. C. (1952). "The Planets," pp. 49-56 and 168--171. Yale Univ. Press, New Haven, Connecticut. UREV, H. C. (1955). The cosmic abundances of potassium, uranium, and thorium, and heat balances of the Earth, the Moon, and mass. Proc. U. S. Natl. Acad. Sci. 41, 127-144. UREY, It. C. (1957). In "Progress in Physics and Chemistry of the Earth," Vol. 2, pp. 46-76. Pergamon Press, London. UREY, H. C. (1962). Origin and history of the Moon. In "Physics and Astronomy of the Moon" (Z. Kopal, ed.), Chap. 13, pp. 481-523. Academic Press, London and New York.