Physics of the Earth and Planetary Interiors, 31(1983)145—160
Elsevier Scientific Publishing Company, Amsterdam
—
145
Printed in The Netherlands
Earth and Mars: early thermal profiles Angioletta Coradini, Costanzo Federico 1 and Pasquale Lanciano Istituto di Astrofisica Spaziale del CNR, Reparto di Planetologia, dale dell’Università 11. 00185 Roma (Italy)
(Received July 6, 1981; revision accepted August 10. 1982)
Coradini, A., Federico, C. and Lanciano, P., 1983. Earth and Mars: early thermal profiles. Phys. Earth Planet. Inter., 31: 145— 160. The extent of formation heating for the Earth and Mars has been evaluated assuming that the terrestrial planets accumulated from planetesimals. The main result is that, even if a long accumulation time is assumed (i- 100 Ma), it is possible to obtain a planetary structure with a large melted shell taking into account the role played by massive projectiles, which, upon reaching depths of several kilometres, are able to deposit heat significantly below the planetary surface. Internal temperatures, sufficient for the downward migration of the liquid iron alloy, have been obtained.
1. Infroduction One of the most speculative branches of geophysics is the study of thermal processes within the Earth and more generally within a terrestrial planet. However, the subject is important because heat transfer phenomena, directly or indirectly, control most tectonic activity. Moreover igneous activity, regional metamorphism and the response of rocks to stress all depend on thermal control. Sources of heat within a planet are needed to explain both the present heat flux and the internal temperature. Among the possible heat sources are long lived radioactive isotopes, short lived isotopes, energy released during core formation and heat produced during the planet’s formation. In this paper we shall focus on the last process, because several lines of evidence are consistent with the inference of a hot stage soon after the formation of terrestrial planets. We will investigate the extent of accretional heating of a planet formed
Also at Istituto di Geologia deIl’Università degli Studi di Perugia (Italy) 003 l-920l/83/0000—0000/$03.O0
by accumulation of planetesimals. Our goal is to obtain an initial temperature profile that can be used as a starting condition in modelling planetary thermal evolution.
2. Process of planet formation The amount of heat retained in a terrestrial planet during accretion depends on how the formation process took place. In this paper we assume that planets accumulated from planetesimals, as a great deal of earlier work suggests (Safronov, 1969; Weidenschilling, 1976; Greenberg et al., 1978; Vityazev et al., 1978; Wetherill, 1978; Greenberg, 1982; Wetherill, 1980; Coradini et al., 1981). This theory assumes that a swarm of planetesimals are in orbit with low eccentricities and inclinations around the Sun. Due to mutual collisions, these bodies grow, mutual gravitational interactions strengthen and their relative velocities increase. Safronov (1969) has analytically demonstrated that a swarm of planetesimals, differentially rotating around the Sun, would have a relative equi-
© 1983 Elsevier Scientific Publishing Company
146
v
=
m2, 9 3, then from (4) T 130 Ma. This figure is a general result of the theory and(1980). has also been numerically obtained by Wetherill During the growth, the rate of bodies infalling on the planetary embryo is given by ii ( m ) dm Cm ~‘dm (5) =
librium velocity given by (1)
(GM/8R)1~2
where M and R are the mass and the radius of the largest planetesimal in the swarm, 9 is the Safronov number, which depends upon relative velocities, impact dissipation, collision cross sections and mass distribution of planetesimals. Wetherill(l980) has numerically found that the steady state relative velocity varies with the mass of planetesimals in the way predicted by Safronov (1969). In fact, the results of analytical and numerical studies are in agreement, predicting values of the Safronov number in the range 1—S. In determining the values of the Safronov number, a continuous mass distribution has been assumed, with most of the mass in
=
=
where n(m) is the incremental number for density of masses of size m, C is a constant, q, based on analytical and numerical investigations, has a value in the (Wetherill, Duerange to its l.67—2.00 gravitational field, the1980). largest body in the swarm (the planetary embryo) grows more rapidly than the other bodies and becomes more and more massive, collecting mass from a “feeding zone” that enlarges proportionally with the increasing radius of the embryo (Safronov, 1969). By taking into account the enlargement of the feeding zone, the depletion of the feeding zone as the embryo grows and by assuming that a 0, the original surface density of solid material within the feeding zone is constant, the rate of growth of the embryo has been shown by Vityazev et al. (1978) to be 2R2)/B (3) R (1 —A where A2 (2a5”~2a 2 and B 0) l(M©Op/6,r)V (p/a 0)P/(l + 29). Here M© is the mass of the Sun, p the internal density of the growing embryo, a the radial distance of the embryo from the Sun, P the orbital period. In this model the growth time of the embryo to reach up to 97% of the present mass of the corresponding planet is given by =
=
T
=
=
(4)
5.3B/(2A)
Assuming for the Earth, at a
=
1A.U., a0
=
100 kg
-
and it can be shown that C 4~(2 q)pR2R/m~~ (2 =
—
=
—
q)A1/m2” max
with mmax being the mass of the second largest planetesimals in the swarm and i~fis the mass rate of growth of the given embryo. Forthe planetesimal with relative velocities by (1), impact velocity is V 1
the largest bodies. If a mass distribution function in the form of a power law is considered n(m)dm Cm~dm (2)
=
=
V~[
1
+
(1/20)]
1/2
(6)
where Ve is the escape velocity from the embryo. In order to specify all the parameters mentioned above, it is necessary to determine mmax. In a steady state, the size of the second largest planetesimal in a feeding zone of a growing planet has been estimated 3by Safronov (1969) to be given by mmax
=
M/(29)
This size is consistent with Safronov’s independent estimates, based on the observed inclinations of the rotation axes of the planets. However Greenberg (1979) has more recently shown that this is only one of a number of possible estimates and that the size of the second largest planetesimal in the Earth’s zone was in the range of a few km to 2500 km. In spite of the uncertainty in mmax it can be argued that large bodies (r’> 1 km) have played a major role in the formation of the terrestrial planets and thus in initially heating the planets selves. it should be noted that a themcrude estimateIn offact several kilometres (Safronov, 1969; Wetherill, 1976) is given for the radius of a planetesimal able to deposit heat at depth.
3. Impacts as energy sources A large number of theoretical and experimental work has been devoted to impact processes during the past decades; recently critically reviewed by Kieffer and Simonds (1980).
147
From our point of view, two problems need careful attention, the penetration depth of the position. During the impact, a fraction of kinetic projectile and the depth distribution of heat deenergy goes into the target. Half of this energy is transferred into kinetic energy of the particles, that have velocity u~, behind the shock (Gault and Heitowit, 1963). Imposing the condition of momentum continuity at the projectile target interface during the penetration, a simple quadratic equation, with v, and u~,results. If it is assumed, as suggested by Kieffer and Simonds (1980), that the projectile penetrates into the target during the time interval taken by the compression and rarefaction wave to travel through the projectile from the first contact back to the target—projectile interface, then the penetration depth is d
4r’u~/U
(7)
5
where U~is the shock wave front velocity, and r’ the projectile radius. Note that d is different from the crater depth, being a measure of the interface penetration and thus the depth at which the energy is initially deposited. The same chemical corn-
-
dir’
..- - - -
2
/ / -
/ ,~
i
,~
BASALT
—
—----—
DIABASE
-
DIABASE
GRANITE
—
GRANITE
D
_~_,
10, ...
Fig. 2. Geometry of the impact process: d is the projectile penetration depth; r~is the radius of the sphere centered in D. Inside this sphere the pressure is uniform and at its maximum value.
position for the projectile and target is assumed. The penetration depth also depends on the impact velocity. From Fig. 1 it is noted that, in the range of impact velocities characteristic of the accumulation of terrestrial planets, d 2r’. Thus, massive projectiles, reaching depths of several tens of kilometres, begin to deposit heat significantly below the planetary surface. Safronov (1978) and Kaula (1979), neglecting the penetration in their models, have underestimated the different heating efficiency due to small (r’ < 1 km) and large (r’> 1 km) planetesimals. It is assumed, following Kieffer and Simonds (1980), with regard to the attenuation of the shock wave, that the energy transmitted to the target is initially contained in a sphere of radius r0, centered uniformly at depth d shocked (Fig. 2). toInthe thispeak sphere pressure the material F’. After is this stage, the spherical shock front expands symmetrically from its center D and the energy, irreversibly trapped in the target, depends on the processes of compression and rarefaction of the medium. Using the Rankine—Hugoniot equations, a linear shock velocity—particle velocity Hugoniot and the Murnaghan equation of state for the release adiabat, the peak pressure in the shocked material P’ and the specific energy buried. Eu,, can be obtained. In fact subtracting the energy released during the rarefaction stage, from the energy gained by the compressed material, gives
-
/~BASALT
~d
1)
5 10 IS V,(Km s Fig. I. Penetration depth of a projectile normalized by radius r’ versus v. for a target and projectile of the same chemical composition.
E~= l/2P’(~— v’)_f V~P(v) dV (8) It is to be pointed out that F’ = F’ (,~),where r~.is
the radial distance from D; V’ and k~are the
148
specific volume at pressure F’ and P,,, respectively. In order to evaluate E~(i~.),it is necessary to determine the damping of the peak pressure when
jectile E~going into heat is
h
=
l/~[4/3~
d~I
~
r~,> r
0. During the propagation of the shock wave, the energy initially transferred into the target is partly lost by irreversibly heating E~ and partly available for further propagation of the shock through the target EA. Thus an energy balance can be written E~ + E~ constant (9) =
=
where ET is the total energy input. The “residual energy” EA, according to Gault and Heitowit (1963), can be expressed as EA
=
4/31r!~[
—
l/2P’(l
l/VJP( V)
—
V’/V0)
dVl
(10)
where it is assumed that the pressure is constant within the compressed material inside the sphere of radius re.. The “irreversible heating” E~is obtamed integrating eq. 8
(14) where ,~* Sr0 is the distance corresponding to where the buried heat becomes negligible. The value of h ranges from 0.12 to 0.43 (Fig. 3) and depends, on the impact velocity and on the cornposition of material. Such values of h are in agreement with those obtained by O’Keefe and Ahrens (1977), who used a different physical approach. On the other hand the values obtained for h are slightly different from the values adopted by Kaula (1979) in his standard model. The value of h, measuring the heat retained from the impact, is critical and thus its value strongly affects the thermal profile. Our investigation leads to a different law for the distribution of heat with depth, compared to the results of Safronov (1978) and Kaula (1979). In fact Safronov (1978) assumes, for sake of simplicity, a linear decreasing dependence on depth in a plane parallel layer, while Kaula (1979) adopts an exponential decrease.
E~= 4irr~. E’~(r0)/(3V,~) +47T/J’çJ
~
(11)
From eq. 9 it follows
______________________________
H
dEA/dr(.= —dE~/dr~.
Using eqs. 8 and 11, it is possible to obtain (Lanciano, 1981) the relation valid for r,.> r,, 0.4
~=,~expf~l/f(P)dP
(12)
where of the f( material P) depends the theimpact chemical velocity composition through P’. 8and andonon 12 determine, for From every eqs. pressure P’,it is thepossible energytoirreversibly trapped and the corresponding radial distance The value couples (E~,r~.)obtained can be fitted, giving
0 2
~,
E~,(r~) = E~1.~(r0)(r~/r0)”
DIABASE
-‘~‘ / .
—
—
GRANITE
El -
DIABASE
—
GRANITE
(13)
The power law index n is weakly dependent on impact velocity v• and ranges from 4.8—5.1. The fraction of the kinetic energy of the pro-
5 10 IS v, (Km s~) Fig. 3. Fraction of the kinetic projectile energy h irreversibly trapped into the target as a function of the impact velocity v,.
149
We show elsewhere (Appendix) how eq. 13 can be modified to obtain the deposited heat as a function of the distance from the center of the growing embryo, r. Taking into account the mass spectrum of the impacting bodies (2), the rate of heat generation during the accumulation of the planet can be shown as
TABLE I Energy output of long lived radioactive isotopes (after Stacey, 1977) Isotope
C,~, (10-8 g g~’)
X, (10-
238U
1.99
1.55
95
2.04
235U 232Th
0.01 7.00 4.67
9.84 0.49 5.54
562 26 30
92.4 1.25 12.8
‘° y’)
A,
(~W 1) kg
exp(A,T*)
pr
E1(r)
=
~ j=1
J
2, r~, e,(r,
(15)
r’)h(r’) dr’
the summation is carried out on the three subsets in which the planetesimal can be grouped: (I) d>z, (2) d
4. Equation of heat transport for a growing planet To quantitatively determine the thermal profiles of a growing planet, it is necessary to solve the differential equation of heat transport. For a spherically symmetric planet with constant internal density, 2 )( a/ar )[ r 2k( aT/ar )] + ( E/pc~) aT/at ( l/r (16) =
element T* years ago is given by C, C01 exp(A1T*) and the summation is carried out over the four isotopes. The effect of the short lived radioactive elements is ignored because of the uncertainty of their initial concentrations and of the time between nucleogenesis and accumulation of the planets. Another effect that must be included is the rise in temperature caused by adiabatic com=
pression as the internal pressureThe progressively increases during accumulation. temperature profile due to adiabatic compression is T(r)=7exp[(21r/3c~)apG(R2_r2)J
(19)
where T is the temperature, c~,is the specific heat at constant pressure, k is the thermal diffusivity and E the energy generation rate. The energy generation rate due to impact E 1 has already been determined in the previous section. However, there is another source of energy due to the long lived radioactive isotopes so that the total energy rate results
with I~ the temperature at r R(t), and a the volume coefficient of thermal expansion. Obviously the temperature profile can never be lower
E
(17)
where k1 is the “eddy diffusivity” connected with
The rate energy production A, of the four 238U,of 235U, 232Th, 40K is shown together isotopes with the decay constants in Table I. An isotope
the impact (Safronov, 1978), k~ is thermal stirring diffusivity and k~1969, is used to mimic heat transport by convection. k 1 takes into account that during accumulation the dominant mode of heat transfer is through the mixing of previously buried material by the subsequent impacts. Thus k, is a measure of the displacement of material at different temperatures in the complex process of crater formation through ejection and late stage gravitational modifications.
=
E1 +
ER
with decay constant X, was more abundant in the planet by a factor exp(X1i~*),where r* 4600 Ma. Therefore the rate of energy production is =
ER =>A1C1 exp(—X,t)
(18)
where the average concentration of a radioactive
=
than the one given by (19). The diffusivity in eq. 16 can also be expressed as k
=
+
k~+ kV
(20)
150
According to Safronov (1978) k/k
r’
kj(z) =f”/3[H2(r’)/2]{~R~(r’)/ [4~(R_z)2]}.h(r~)
I
(21)
dr’
and we have assumed that the mixing of materials is confined to a cylindrical volume ~rR~H. In eq. 21 r’ is the minimum radius at which a planetesimal is able to mix materials at depth z, r,~axis the radius of the largest planetesimal in the distribution, /3 1/3 and is equal to the square of the ratio of the mixing length to the crater depth H (Safronov, 1969 p. 177). h(r’) is the rate of the infalling bodies and can be obtained from (5). The dimensions of the crater H and R,. have been determined using the experimental data of Gault (1973) and the scaling laws suggested by Safronov (1978). The scaling law depends on the ejection coefficient h e that is the fraction of the kinetic energy of the projectile partitioned into ejecting and moving materials, In eq. 20 the thermal diffusivity k~is the sum ,
.
.
of the diffusivity due to lattice vibration and radiative heat transport, and depends on pressure and temperature. A reasonable mean value, in the range of temperature 2relevant for our ~ (Schatz and numerical Simmons, calculations, is l06disregarded m 1972). We have the poorly known effect of pressure on k~.In the outer layer of the planet, where the impact mixing is efficient, k~is negligible with respect to k 1, contrary to the inner zones, where, before the occurrence of convection, only thermal diffusivity is responsible for heat transfer. Figure 4 depicts the values of the ratio k,/k~ at various depths. Substantial quantities of heat can be transferred upward by the process of thermal convection during the accumulation of the planet. In the heat transport equation, this effect can be simulated by introducing a pseudodiffusivity kV represented by 1/3
k~= Nu k~= (Ra/Ra~) k~ (22) where Nu is the Nusselt number, defined as the
ratio of the convective heat flux to that which would be transported for the same by conduction alone. Ra is the Rayleigh number which entirely governs the onset of convection, being essentially
100
I0
2
I
________________________________________ 400
800
Z(Km)
. . Fig. 4. Ratio of the eddy diffusivity k, on the thermal diffusivity k~ as a function of the depth and at different times during
the growth of the Earth using the standard model. The zone of the planet where k 1 > k~ enlarges as the mean size of the largest bodies increases. In the final phases of the accumulation, when the rate of impacting bodies decreasing, this ratio tends to I. Numbers on the curves give the time in Ma.
the ratio of the driving force (thermal buoyancy) to the retarding force (diffusion of heat and momentum). Ra~is the critical value needed to initiate convection. We have assumed a value equal to 1700. The power law relation (22) between the Nusselt and Rayleigh number is known, from laboratory experiments and theoretical and numerical calculations, to characterize convection in a variety of circumstances. We have calculated the Nusselt number taking into account the strong dependance of viscosity on temperature and of the melting temperature Tm on pressure. The expression used for Tm TmP and v P(Tm, T) are those given in Kaula (1979). The maximum Rayleigh number has been determined for all possible layers where the temperature gradient dT/dr is negative. This method differs from the one adopted by Kaula (1979), who assumes that the convective layer always starts at the level of maximum ternperature. In fact, due to the increase of Tm with depth, the maximum temperature inside the grow=
=
151
ing planet may lie on the plane T z at a point further from the melting curve than the temperature of less hot zones. Consequently the kinematic viscosity, that exponentially depends on can inhibit convection where the maximum temperature is reached. It is possible in (22) to use only k~because the timescale of convective overturn is much longer than the characteristic timescale between two large impacts. In fact the convective overturn characteristic time is given by —
2
2/3
1
\
=
(Stevenson and Turner, 1979). We find that Tv i5 about l0~s, but the mean interval of time ‘r, between two subsequent impacts of planetesimals in the same area is about 1012 s. T, has been evaluated taking into account eq. S and the percent of surface area hit by impacts of kilometresized planetesimals. The surface of a growing planet is also locally melted by impacts, but radiative cooling of these melted zones is efficient enough to drastically reduce the temperature on a timescale shorter than T,. A crude estimate of the time necessary to decrease the local surface temperature by about z.~T 700 K is =
F
T
1
4
4
\1
3
7~=Hcpz1T4a~,T—T)j~xl0H m ~‘
where Hm is the thickness of melted material. For any reasonable value of Hm, ~~7oOis always much shorter than ~ and we assume that kilometre-sized bodies do not affect the surface temperature 1~,. The surface temperature is determined by the infall of small bodies (r’ < I km) and ‘particles, and has a mean value of 300 K (Safronov, 1969, p. 156). From the above discussion, it is clear that eq. 16 must be solved numerically. Using an algorithm of Toksoz et al. (1972), the stability criterion is given by ~t=l(L1r)/k
2
1
2
12m/(2m+l)
2
maxj
where kmax is the maximum value of the diffusivity in each time step ~t, ~r is the radial step and m r/Ltr. To account for the increasing radius of the planet, we have used the method described in =
=
=
=
and when t 0 the initial temperature profile is given by (9), where R 0.O5Rf; Rf is the final radius of the planet. =
L /i,0.2Ra kc) where L is the thickness of the convective layer Tv
Kaula (1979). At the onset of convection, which is the dominant heat transport mechanism during a large part of the life of the planet, .7~t becomes extremely short and, consequently, the computer time increases. However, in order to retain a rigorous treatment of the convective process, numerical procedures were not used to speed up the calculations. The boundary conditions assumed are T(R t) 300 K (dT/dr)r ~ 0
=
5. Results Using the accumulation theory previously described and selecting different meaningful parameters, thermal models of the Earth and Mars, at various stages of their formation, are considered. 5.]. Thermal modelfor the Earth
The “standard” Earth model is characterized by q 1.83, he 0.3, 8 3 and a0 100 kg m2, thus the accumulation time through eq. 4 results T 130 Ma. The internal temperatures, as a func=
=
=
=
=
tion of the radius, are given in Fig. 5 at different .
.
.
.
stages of the embryo growth. In the primordial phases, the energy buried at depth, by impacts, is not enough to melt the body: the peak temperature is at a depth dependent on the penetration of the projectiles and on the conditions of the impact heat transfer. The peak moves inward as the embryo and the impacting bodies Convection 6Rf,grow. 30 Ma after the sets in when the radius is O.5 beginning of the accumulation, and the melting temperature for silicates is reached 2 Ma later at a depth of about 900 km when the radius is O.S8RJ. Convection inhibits a rise of temperature beyond the melting point and the thermal profile lies exactly on the melting curve. At T 130 Ma with the Earth’s present radius, a maximum temperature of 3300 K is reached at 3600 km depth i.e. deeper than the present core—mantle discontinuity. Convection occurs throughout a layer 3300 km .
=
152
T(K)
3000
2000
1OOO~1~
0.2
•
04
06
•
•
0.8
R/Rf
Fig. 5. Early thermal profiles for the Earth standard model. Numbers on the curves give the time in Ma. The curve at 42 Ma is the first to clearly indicate a large melted shell. The dotted line is the melting curve for silicates when the planet reaches its present mass. The distance between two dots corresponds to 50 km.
thick, corresponding approximately to the melted layer. The Rayleigh number at this stage is very high (~6 x 1011) and this vigorous convecting motion seems to preserve the planet homogeneity. The shape of the thermal profile is mainly de-
termined by the external heat source due to the infalling bodies, with the total amount of trapped energy at the end of the accumulation of the Earth at 8% of the total gravitational energy of the planet. Such an amount of conversion of gravita-
T(K)
3000
....
0:2
0.4
•
0.6
0.8
R/R~
Fig. 6. Early thermal profiles for a “cold” Earth. In this model the melting of a large shell is delayed with respect to the standard model.
153
tional to thermal energy in models where no attention has been paid to the physics of the formation process, has been obtained only for unrealistic rapid process (T l0~y) (Hanks and Anderson, 1969; Hiriyama, 1977). A comparison of our final thermal profile with Kaula’s (1979), shows that lower temperatures are reached because the drastic reduction of viscosity at the melting point makes the convection so efficient that the melting temperatures cannot be exceeded. Kaula (1980), reexamining his previous results, had foreseen this behaviour, In order to assess the validity of our standard model, we have drastically modified the set of parameters to obtain “minimum” temperatures. Figure 6 details the early thermal profiles when q 1.9, 9 4, he 0.5 and T 120 Ma. This choice really corresponds to a minimum energy input where, in fact, 50% of the impact energy is lost through ejecting materials (he 0.5). The relative velocities are low (0 4) and q 1.9 implies a mass distribution with 50% of the total mass in bodies less massive than iO~mm,,,,. Otherwise the usual assumption q 1.83 means that only 30% is in bodies less massive than l0~ mm,,,,. Even in this case, temperatures are high enough at the end of the accumulation to reach the melting point in a shell 2000 km wide with convection setting in =
=
=
=
=
=
=
=
before the completion of planet formation and remaining vigorous until the end of the accumulation itself. However, the trapped energy fraction is 7% of the total gravitational energy and does not significantly differ from the value obtained for the standard model. It must be stressed that while the mass distribution of planetesimals is critical in enlarging or narrowing the melted shell, the final thermal profile to a first approximation is quite insensitive to a reduction of the accumulation time, because both the eddy diffusivity and the impact energy generation rate are proportional to the accumulation rate. However, a more detailed treatment of the accumulation and consideration of the link between the accumulation rate and planetesimal velocity distribution, reveals a weak dependence of the thermal profile on the formation time, through the modification of the velocity distribution and the mean mass of the planetesimals. The above mentioned effect is true as long as the accumulation time is shorter than the decay time of the radioactive elements: on a very long time scale the radioactive energy sources also become effective. The importance of convection is considered in Fig. 7 where curve b, in the absence of convection, rapidly overcomes the melting curve. Convection is very important, even with the accumulation of
T(K)
3000
0.2
0.4
0.6
0.8
R/R r
Fig. 7. Standard Earth model: comparison of thermal profiles with (a) and without convection (b). When convection is absent the profile rapidly overcomes the melting curve. In the presence of convection there is an outward enlargement of the molten shell.
.
154
large bodies, that deposit heat at depth, extremely high (supersolidus) temperatures need not be reached. The small differences between the results of the two models affirm that they are not critically dependent on the uncertainties of the adopted parameters. Therefore, the obtained early thermal profiles can be used as initial conditions for further, more realistic, thermal histories. Many studies have been carried out on the subsequent thermal history of the Earth, however, due to the peculiarities of our primordial profile, we sought to verify if, in this case, reasonable evolution could be obtained, We have neglected the process of core formation in our model because of the uncertainties of the role played by light elements in controlling the melting temperature (Brett, 1976). The lowest melting composition for any initial Earth is probably eutectic system Fe—FeS, with a eutectic melting temperature in the range 1000 K at zero pressure, to 1500 K for pressures at the centre of the Earth. So it is possible to argue that, in our model, core formation begins to occur in a very primordial phase when the Earth embryo has a radius 3500 km (Figs. 5, 6). On the other hand for an
alloy with a higher melting point (Fe, Fe—Ni), core differentiation probably occurs during melting of the silicate mantle. Figure 8 depicts the thermal profiles at different times (t T), taking into account convection, heating of radioactive elements, and neglecting any migration of these elements. The temperature and the Rayleigh number decrease with time, because of the decreasing energy released by radioactive elements. The final profile at I 4600 Ma is far from the melting curve of silicates in the mantle and, therefore, an energy release in the range of 2.5 X l0~—4x l0~cal kg ‘(Birch, 1965; Murthy and Hall, 1972) for core formation may be easily accommodated throughout the process. From this figure it is possible to deduce that, even if convection endures for the life of our planet, the primordial heating remains important. However, the thermal profiles in Fig. 8 cannot be used to infer the thermodynamical properties of the material inside the Earth and the timing of the onset of the geomagnetic field, because the core formation, which may have played a dominant role in the subsequent thermal history, has been neglected. Our results seem to indicate that gravitational differentiation began in an intermediate =
T(K) 3OOo~~~:
.,.
02
0.4
0.6
0:8
R/Rf
Fig. 8. Standard Earth model: subsequent thermal evolution neglecting core formation. The central temperature increases due to the energy released by long lived radioactive elements. The Rayleigh number decreases from 6>< 10’ at 130 Ma to 7 x 106 at 4600 Ma. The final profile is far from the silicate melting curve.
155
shell, where the melting temperatures were reached early. Consequently, heating and Rayleigh—Taylor instabilities are involved in the process of core formation. Numerical simulation of core and mantle formation, based on this mechanism, has been performed by Vityazev and Majeva (1980). Their
heat sources (Arvidson et a!., 1980). The foregoing consideration forced a computation of initial thermal profiles for Mars. Unfortunately, there is not a general consensus on the Martian formation time scale as in the case for the Earth. Using eq. 4, assuming 0 2 and a0 7 kg 2, T 1500 Ma. The selected value of the mass m surface density in the vicinity of Mars would have been larger than the one obtained by spreading out the present mass of Mars in its feeding zone. The reason is that the gravitational perturbation of Jupiter would have ejected material out of the Martian feeding zone and increased the planetesimal relative velocities (i.e. low values of 0). On the other hand, as previously stated, numerical studies of simultaneous accumulation of terrestrial planets seem to indicate that in 150 Ma the terrestrial planets reach 99% of their present mass (Wetherill, 1980). Therefore, because of these uncertainties, two different models have been studied. In both models h 0 0.3, 0 2 and q 1.83, but in model I the accumulation rate R has been assumed constant over a period T 100 Ma, while in model II R is given by eq. 3. The thermal profile at the end of accumulation in model I (Fig. 9) shows, at the centre, a tempera=
=
=
results indicate that thecan geochronological andmake geomagnetic constraints be fulfilled. To realistic predictions on the core formation and onset of a geomagnetic field in our model it is also necessary to simulate the process of separation of heavy components from the mantle inside the molten shell. 5.2. Thermal model for Mars
Our present knowledge of the interior of Mars is constrained by lack of data. On the basis of the mean value of density and moment of inertia it can be inferred that this planet is to some degree differentiated, implying that it has been heated
=
during its history. Thermal evolution models have been constructed, considering the formation of the crust, mantle and core, but the time scale can be lengthened or shortened by making various assumptions regarding the initial temperature and
=
=
=
TI K)
2000
100
1000
02
0.4
0.6
0.6
R/Rf
Fig. 9. Early thermal profiles for Mars model I. At the end of accumulation a melted shell 360 km wide is present. The dotted line, corresponding to the melting curve of silicates, has dots equally spaced at 30 km. The convection in the last profile is moderate (Ra l0~).
156
ture of 410 K and this small increase is also due to adiabatic compression. Even if an accumulation time similar to the Earth’s has been assumed, the melting temperature of silicates is reached in a shell 360 km wide with a maximum temperature of 1850 K at a depth of 1300 km. So it is not necessary to assume an unrealistically short time ( iO~y) of formation in order to have a hot beginning for Mars, which is a relatively small planet. In the post-accumulation thermal simulation the melted shell increases rapidly, continuously enlarging (Fig. 10), due to the decay of long lived radioactive elements, until after a few hundred million years (I 460 Ma) it is 1500 km wide. At this time 68% of the planet mass is melted, with the lithosphere thickness reaching a minimum around 1800 Ma, at which stage the process of core formation has probably ended. The thermal consequence of this is a significant increase of central temperatures and a negligible rise in the outer melted zones, because convection can rapidly remove heat from the outer melted zones where roughly 30% of the released energy goes. Only in the latest stages, after t 3200 Ma, does the energy source become insufficient to maintain a melting condition and the lithosphere thickens,
In model I, the amount of gravitational energy stored into the planet during growth is about 20% of the eventual total gravitational energy, which is larger than the fraction trapped in the Earth because for Mars a moderate convection only acts during the final stage of formation. In model II (Fig. 11) the thermal profile depends on ‘r, because the eddy diffusivity k, becomes comparable with k~. At the same time radioactive heat generation rate is comparable with the impact heat generation rate, thus the early thermal history overlaps the subsequent thermal evolution of the planet. The main effect is a higher degree of heating for the whole planet and an enlargement of the melted shell (— 1600 km). Moreover, because of the radioactive heating the central temperature reaches 1200 K. Comparing the results obtained in both models, it is evident that, whatever the accumulation time, after roughly 1500 Ma the thermal profiles evolve in the same way. With reference to the detailed discussion of Arvidson et al. (1980) on thermal evolution and geological constraints, the initially hot Mars we obtained, has had its thermal history moved backward in time with respect to the Solomon (1978) and ToksOz and Hsui (1978) models. However, our
=
=
T( K)
.
•
0.2
04
0.6
0~
R/Rf
Fig. 10. Mars model I: further thermal evolution, neglecting core formation. The final profile moves away from the silicates melting curve and the Rayleigh number is about l0~.
157
T(K) 2000
.
1500
1000
0.2
0.4
0.6
0.8
R/Rf
Fig. II. Mars model II: in the last profile the joint effect of impact energy and radioactive energy produces a larger melted shell.
present thermal profile (Fig. 10) t 4600 Ma is still hot enough, because of the larger amount of accumulation energy stored in the planet. The post Viking view of Mars (Arvidson et al., 1980), suggesting an early crust differentiation, ancient volcanism, extensional and compressional features, seems to confirm our deduction of early and efficient accretional heating. =
sion for the rate of energy per unit volume as a function of the distance r from the centre of the planet (Fig. 12), is now required to insert in the equation of heat transport (16). When a planetesimal with radius r’ and impact velocity v, hits the planet surface, it penetrates a distance d into the ground (Fig. 12). To simplify the calculation, three different regions are considered at different depths with z R r (Fig. 12). =
—
Acknowledgements
A.]. First region: 0
We wish to thank Professors R. Arvidson and D. Tozer for helpful discussions and improvements to this work.
Here the mixing of material is very efficient and the energy is assumed to be uniformly deposited. A symmetrical distribution of the energy, requires that half of the energy is irreversibly trapped in this region. Recalling that the fraction of the impact energy E, irreversibly trapped is h (14), the mean energy per unit volume e1(r’), deposited in this first zone, is
Appendix We have already shown that the specific energy irreversibly trapped as a function of the radial distance t~.can be expressed as follows
z
d
—
e,(r’)
=
(3/8ir)hE,(r’)[R3
—
(R
—
d)]
(2A)
E~~(r~)=E~(r0) if rer,,
~
~
(lA)
where r,, r’ when the projectile and the target have the same chemical composition. The expres=
A.2. Second region.’ d
Where r0 is the radius of the sphere centered in D (Fig. 13) where the whole energy hE1 is initially
158
(
~/ / /
Fig. 12. Adopted geometry in the calculation of the rate of impact heat generation.
trapped. The mean energy per unit volume in this second zone can be further subdivided into two contributions: (1) the energy trapped inside the sphere with radius r0 and (2) the energy deposited outside this sphere, to give 2sin 9 dO e2(r, r’)
=
(p/4~r2)[
+
j
where 0. is the angle between OD and OS (Fig. 12) and S is a point on the intersection between the sphere of radius r,, centered in D, and the sphere of radius r centered in 0. The eq. 3A can be integrated taking into account that r~ (r2 + u2 2ru cos =
—
j°”E~(r)2~r
2E~~(r,,)(ç/r)”2irr2 ~ 0 del (3A)
where u
=
R
—
d. It is then possible to obtain
e 2(r, r’) = (p/2)E~(i~){(I— cos 0~)
~~cond~’
third region
Fig. 13. Sketch of the three different regions as defined in the appendix. d is the penetration depth and r,, is the radius of the sphere where the maximum shock pressure is uniformly reached.
159
2/2
x [(u2
~2_n}
—
where rj2=z/(4u~/U~+ 1) and rç2 < r~,a,,
r~2=zL~,/4u~. If
}
(4A)
+ r2)(
=
E,3(r) =A/(nru)(n
2)j
E~(ro)(re/ro)
2irr2
p/(4irr
=
Xsin 9 dO E~(,~,)i~,”p/{2r(2 — n)u] 2~1V2 — (r —
=
u)2~I] (5A) x [(r2 + u2)( From eq. 14 a useful expression for the constant
3)/(2
—
—
n)J rL, r’~ —
3 0”
r)~”j dr’
(8A) where r 1’3 ~ and r2’3 z/(4u~/lJ~+ 1). The choice of the value of ~ does not affect the numerical results as long as ~ is sufficiently small: r,~jn/r,~,a,, < 10~. =
0
—
x [(u2 + r2)(2~~V2 (
Finally in the third region the energy per unit volume results in 3(r, r’)
then
3”~r
A.3. Third region: z > d + r,,
e
rj’2
=
=
References
E~(r 0) is derived, which appears on the right side of eqs. 4A and 5A. It results from E~(r0)= 3(n — 3)/(4~rn)hE,/(pr,fl h/2(n
=
—
3
3)/n(r’/r0)
phys. Space Phys., 18: 565—603. Birch, F., 1965. Speculations on Earth’s thermal history, Geol. Soc.R., Am.1976. Bull.,The 76: current 133—154.status of speculations on the Brett,
2
V1
The rate of heat generation E1(r’), during the accumulation of the planet, can now be obtained integrating numerically the expressions 2A, 4A and 5A on the impact size spectrum. From eq. 15 2~r’/r’14~’1 max
h(r’) dr
where
p
E,~(r)
3(4 —p)R 3q 2. It follows that
= =
=J
for hypervelocity impact craters formed in rock. Proc. VI Hypervelocity Impact Symp., 2, Firestone Rubber Co., Cleveland, OH, pp. 420—456.
2.kp/r’4~ max
Icarus, 39: Greenberg, R., 141—150. 1979. Growth of large, late-stage planetesimals. Greenberg, R., 1982. Planetesimals to planets. In: A. Brahic (Editor), Formation of Planetary System. Cepaokies-Edilions, Toulouse, France, pp. 5 15—569. Greenberg, R., Hartman, W.K., Chapman, C.R. and Wacker, J.F., 1978. The accretion of planets from planetesimals. In: T. Geherels (Editor), Protostars and Planets. The University of Arizona Press, Tucson, AZ: pp. 599—624. Hanks, T.C. and Anderson, D.L., 1969. The early thermal
—
2r’(3~p){R3
E
—
11(r)
=
Astrophys., 98: 173—185. Gault, D.E., 1973. Displaced depth, diameter and of oblique trajectories formass, impact craters formed in effects dense
1(r, r’)h(r’) dr’
3(4 p)/4hv~R the heat generation in the three different regions results as follows =
Phys., 14: 375-383. Coradini, A., Federico, C. and Magni, G., 1981. Formation of planetesimals in an evolving protoplanetary disk. Astron.
e
and where A
composition of the core of the Earth. Rev. Geophys. Space
cristalline rocks. Moon, 6: 32—44. Gault, D.E. and Heitowit, E.D., 1963. The partition of energy
—
fr~
Arvidson, R.E., Goettel, K.A. and Flohenberg, C.M., 1980. A post Viking view of Martian geologic evolution. Rev. Geo-
J~R
—
d(r’)]
3)
dr’
Af
(6A) where r 1’1 = zL’ç/4u~ and r~’1= r,~,,,,,. If r~’1> E11(r) = 0. E12(r)
=
A(n
—
+ r,,”/[
x[(
3)/nj
(2
—
u2 + r 2
22r1(3_P)/r3{(l r~2
—
cos
ri,,,
history of the Earth. Phys. Earth Planet. Inter., 2: 19—29. Hiriyama, J., 1977. Energy balance in the Earth interior. Tectonophysics, 41: 243—249.
°~)
Kaula, W.K., 1979. Thermal evolution of the Earth and Moon growing by planetesimal impacts. J. Geophys. Res., 84: 999—1008. Kaula, W.K., 1980. The beginning of the Earth’s thermal evolution. In: D.W. Strangway (Editor), The Continental Crust and its Mineral Deposits. Geol. Assoc. Can. Spec. Pap., 20: 25—34.
n )ru]
)(2—n)/2
—
r02
“])
dr’
(7A)
160 Kieffer, SW. and Simonds, C.M., 1980. The role of volatiles and lithology in the impact cratering process. Rev. Geophys. Space Phys., 18: 143—181. Lanciano, P., 1981. Profilo termico iniziale di pianeti terrestri e satelliti gioviani. Laurea thesis, University of Roma, 114 pp. Murthy, W.R. and Hall, H.T., 1972. The origin and chemical composition of the Earth’s core. Phys. Earth Planet. Inter., 6: 123—130. O’Keefe, J.D. and Ahrens, T.S., 1977. Impact induced energy partitioning, melting and vaporization of terrestrial planets. Proc. Lunar Sci. Conf. 8, pp. 3357—3374. Safronov, V.S., 1969. Evolution of protoplanetary cloud and formation of the Earth and the planets. Nauka, Moscow, Translation, 1972, Nasa IT F-67, 206pp. Safronov, VS., 1978. The heating of the Earth during its formation. Icarus, 33: 3—12. Schatz, S.F. and Simmons, G., 1972. Thermal conductivity of the Earth materials at high temperature. J. Geophys. Res., 35: 6966—6983. Solomon, S.C., 1978. On volcanism and thermal tectonics on one plate planets. Geophys. Res. Lett., 5: 461—464. Stacey, F.D., 1977. Physics of the Earth. J. Wiley, New York, 4l4pp. Stevenson, D.J. and Turner, J.S., 1979. Fluid models of the mantle convection. In: M.W. McElhinny (Editor), The
Earth: Its Origin, Structure ahd Evolution. Academic Press, London, pp. 227—264. ToksOz, M.N. and Hsui, A.T., 1978. Thermal history and evolution of Mars. Icarus, 34: 537—547. Toksoz, M.N., Solomon, S.C., Minear, J.W. and Johnston, D.H., 1972. Thermal evolution of the Moon. Moon, 4: 190—213. Vityazev, A.V., Pechernikova, G.V. and Safronov, V.S., 1978. Limiting masses, distances and times for the accumulation of planets of the terrestrial group. Soy. Astron., 22: 60—63. Vityazev, A.V. and Majeva, S.V., 1980. Simulation of the Earth’s core and mantle formation. Phys. Earth Planet. Inter., 22: 296—301. Weidenschilling, S.J., 1976. Accretion of the terrestrial planets II. Icarus, 27: 161—170. Wetherill, G.W., 1976. The role of the large bodies in the formation of the Earth and Moon. Proc. Lunar Sci. Conf. 7, pp. 3245—3257. Wetherill, G.W., 1978. Accumulation of the terrestrial planets. In: T. Geherels (Editor), Protostars and Planets. The University of Arizona Press, Tucson, AZ, pp. 565—598. Wetherill, G.W., 1980. Formation of the terrestrial planets. Ann. Rev. Astron. Astrophys., 18: 77—113.