A thermo-chemical regime in the upper mantle in the early Earth inferred from a numerical model of magma-migration in a convecting upper mantle

PHYSICS O F T H E EARTH ANDPLANETARY INTERIORS

ELSEVIER

Physics of the Earth and Planetary Interiors 94 (1996) 187-215

A thermo-chemical regime in the upper mantle in the early Earth inferred from a numerical model of magma-migration in a convecting upper mantle Masanori Kameyama a, *, Hiromi Fujimoto a, Masaki Ogawa

b

a Ocean Research Institute, University of Tokyo, 1-15-1, Minamida~ Nakano-ku, Tokyo, 164, Japan b College of Arts and Sciences, University of Tokyo, 3-8-1, Komaba, Meguro-ku, Tokyo, 153, Japan Received 16 March 1995; revision accepted 15 August 1995

Abstract

A numerical model of mantle magmatism in a convecting upper mantle has been developed to study the thermo-chemical evolution of the upper mantle of the early Earth. The solid-state convection in the upper mantle is modeled by a convection of a binary eutectic material with a Newtonian temperature-dependent rheology in a two-dimensional rectangular box placed on a heat bath as a model of the lower mantle. The density depends on the chemical composition and melt-content as well as temperature of the material. The material contains heat-producing elements that are incompatible and exponentially decay with time. Mantle magmatism is modeled by a permeable flow of melt generated by a pressure-release melting induced by the solid-state convection through matrix. The permeable flow is driven by a buoyancy due to the density difference between the melt and the matrix. The thcrmo-chemical evolution in the box occurs in two stages if the deeper part of the box is not so strongly depleted in heat-producing elements in spite of the upward migration and concentration of heat-producing elements into a crustal layer along the top surface boundary due to magmatism. In the earlier stage, active magmatism occurs because of a strong internal heating due to the heat-producing elements, a chemically stratified structure develops well in the box with dense magmatic products in the deeper part and less dense residual materials in the shallower part, and the temperature distribution becomes strongly superadiabatic over the entire box. The temperature at the base of the box becomes as high as the solidus temperature. The chemically stratified structure is, however, suddenly destroyed by convective mixing and the temperature in the deeper part of the box suddenly drops by several hundred degrees when the internal heat source becomes too weak owing to the decay of heat-producing elements which sustain the active magmatism and hence keep the effect of chemical differentiation due to the magmatism stronger than the effect of convective mixing. In the later stage of the evolution, the box becomes chemically homogeneous and magmatism occurs only mildly. If heat-producing elements are efficiently transported into the crustal layer and the deeper part of the box becomes strongly depleted in heat-producing elements owing to the magmatism, only a mild magmatism occurs even at the beginning, a chemically stratified structure does not develop well, and the temperature in the box rapidly decreases to a stationary value. The regime of hot and chemically stratified upper mantle suggested from the earlier stage of the case with mild depletion of heat-producing elements at depth fits in with many observations of the Archean continental crust.

* Corresponding author. 0031-9201/96/$15.00 © 1996 Elsevier Science B.V. All rights reserved SSDI 0031-9201(95)03102-2

188

M. Kameyarna et al. / Physics of the Earth and Planetary Interiors 94 (1996) 187-215

1. Introduction

One of the most important subjects in the study of the Earth's evolution is to clarify how the Earth's mantle, which was hot just after the formation of the Earth (e.g. Kaula, 1979), has cooled with time and how the chemical structure of the mantle has changed in accordance with the cooling. Some numerical estimates of the cooling history of the mantle suggest that the characteristic time for the cooling of the mantle is several hundred million years to a few billion years depending upon, say, the style of mantle convection and the detail of the rheology of convecting mantle material (e.g. Sharpe and Peltier, 1979; Schubert et al., 1979; Davies, 1980; McKenzie and Richter, 1981; Christensen, 1985; Steinbach et al., 1993; Honda, 1995). The geotherm is most likely to be significantly higher in the early eons of the Earth's history than at present. The high geotherm in the early Earth probably induced a vigorous magmatism in the mantle, particularly in the upper mantle where the solidus temperature is not so high and hence magma can be easily generated. Here, we develop a numerical model of magma-migration in a convecting upper mantle to study how a mantle magmatism influences the thermo-chemical evolution of the hot upper mantle expected in the early Earth. Important constraints on the thermal state of the mantle of the early Earth and hence on the cooling history of the Earth come from the observations of Archean continental crusts. The observations of metamorphic rocks (e.g. Bickle, 1978) and diamonds (Richardson et al., 1984; Boyd and Gurney, 1986) suggest that the continental lithosphere in the Archean is as cold as the present continental lithosphere beneath a shield area. Komatiites, on the other hand, indicate the occurrence of a geotherm much higher than the present geotherm in the Archean (e.g. Bickle et al., 1977; Ohtani, 1984; Takahashi, 1986; Nisbet et al., 1993). The observations suggest that the temperature contrast in the mantle is larger in the Archean than at present (Jeanloz and Morris, 1986). It is not, however, straightforward to clarify the overall picture of the thermal regime of the Archean mantle and hence to constrain the

cooling history of the mantle in the early Earth from the observations. Some of the earlier workers suggest that the average geotherm in the Archean mantle is higher than the average geotherm in the present mantle by several hundred degrees and that the geotherm in the Archean continental lithosphere is much lower than the average because the continental lithosphere is strongly cooled by a vigorous mantle convection beneath oceanic area a n d / o r a root of continents called "tectosphere" (Jordan, 1975) forms beneath the continents (Nisbet and Fowler, 1983; Richter, 1984, 1985). Other workers suggest that the Archean mantle is not so much hotter than the present mantle on average on the basis of the observations of metamorphic rocks and diamonds as well as some of the estimates of the thermal history of early Earth obtained by the use of a simple parameterized convection technique (e.g. Turcotte and Schubert, 1982) and that hot plumes with a large excess temperature uprise through the rather cold Archean mantle to generate a komatiitic magma (e.g. Davies, 1979; Jarvis and Campbell, 1983; Campbell and Griffiths, 1992; Nisbet et al., 1993). The difficulty in clarifying the thermal regime of the Archean mantle from the observations and the simple parameterized convection technique calls for a full numerical mocieling of the thermal regime of the Archean mantle and hence of the thermal evolution of the Earth's early mantle. The thermal evolution of a hot mantle in the early Earth has been suggested to be strongly influenced by mantle magmatism (e.g. Sleep and Windley, 1982; Ogawa, 1988, 1993, 1994; Vlaar et al., 1994). The mantle magmatism in the early Earth is most likely to be much more vigorous than the present mantle magmatism because of the strong dependence of the vigor of mantle magrnatism on geotherm (White and McKenzie, 1989) and the high geotherm expected in the early Earth. A numerical model of a coupled magmatism-mantle convection system (Ogawa, 1993, 1994, hereafter called Papers 1 and 2) suggests that the vigorous mantle magmatism expected in the early Earth makes the upper mantle chemically stratified and makes the geotherm in the early Earth strongly superadiabatic with a

M. Kameyama et al. / Physics of the Earth and Planetary Interiors 94 (1996) 187-215

rather low geotherm in the lithosphere and geotherm much higher than the present geotherm in the deeper part of the upper mantle. Ogawa (1993, 1994) suggests that the temperature distribution obtained in his numerical model is consistent with the observational constraints on the Archean geotherm described above. Mantle magmatism is, however, modeled in a rather a d hoc way and the detail of magma-migration is neglected in Papers 1 and 2. Magma generated in the Earth's mantle migrates as a permeable flow through matrix in partially molten regions (Scott and Stevenson, 1984, 1986, McKenzie, 1984) and / o r as a flow through cracks and dykes formed by magma-intrusion, which is considered to be particularly important for magma-migration in the lithosphere (e.g. Weertman, 1971; Shaw, 1980; Sleep, 1988). The magma migrates upward in the shallower part of the upper mantle but probably migrates downward at greater depths in the upper mantle because of a density inversion between the magma and the coexisting matrix (Stolper et al., 1981; Agee and Walker, 1993). The presumption in Papers 1 and 2 that these details of magma-migration do not affect the overall features of the thermo-chemical evolution of the mantle has not been verified. Here, we develop a numerical model of magma-migration in a convecting upper mantle to study how mantle magmatism influences the thermo-chemical evolution of the hot upper mantle expected in the early Earth; magma-migration is modeled as a permeable flow through matrix. We concentrate on the evolution of the upper mantle only, assuming that the mantle convection in the early Earth occurs as a layered convection. A layered convection is most likely to have occurred in the hot mantle of early Earth even if the mantle is chemically homogenous owing to the transition between the spinel phase and the post spinel phase of (Mg,Fe)2SiO4 (Christensen and Yuen, 1985; Steinbach et al., 1993). (A further discussion on the assumption will be made later.) 2. Numerical model

The solid-state convection in the upper mantle of the Earth is modeled by a convection of a

189

binary eutectic material in a two-dimensional rectangular box of height d = 650 km and aspect ratio (width/depth) w = 1.5 placed on the top of a heat bath, which is a model of the lower mantle. The chemical composition of the eutectic material is denoted by A ¢ B i _ e , where the endmember A stands for olivine and the end-member B stands for a mixture of pyroxene and garnet. The eutectic composition is A0aBo.9, which corresponds to a basaltic composition. The density of the material linearly increases with decreasing temperature T, g, and melt-content ~b except in a thin crustal region along the top surface-boundary; the density is assumed to increase with decreasing ~ in the crustal region to take the effect of basalt-eclogite transition into account. The rheology of the material is Newtonian and the viscosity exponentially depends on temperature. The material contains heat-producing elements that are incompatible and exponentially decay with time. Mantle magmatism, i.e. magma-generation and migration in the mantle, is modeled by a permeable flow of melt (McKenzie, 1984; Scott and Stevenson, 1986, 1989) produced by pressure-release melting. The composition of the melt is calculated from the phase diagram assumed for the binary eutectic material and is equal to the eutectic composition in most cases in the present calculations; a chemical equilibrium is assumed to hold between the melt and the coexisting solid. The melt is, therefore, more enriched in the end-member B as well as heatproducing elements than the coexisting solid. The permeable flow is driven by a buoyancy due to the density difference between the melt and the coexisting matrix. The melt migrates upward in the shallower part of the box but migrates downward in the deeper part of the box in most cases in the present calculation owing to a density inversion between the melt and the coexisting matrix; magma is less dense than the coexisting matrix at shallow depth levels, but is probably denser than the coexisting matrix at greater depths in the Earth's upper mantle (Stolper et al., 1981; Agee and Walker, 1993). The permeability is proportional to the cube of melt-fraction ~b (McKenzie, 1984) at depths greater than 100 km. In order to take into account the effect of magma-migration through cracks and dykes,

190

M. Kameyama et al. /Physics of the Earth and PlanetaryInteriors 94 (1996) 187-215

which is particularly important in the lithosphere (Weertman, 1971; Shaw, 1980; Sleep, 1988), the permeability is assumed to be proportional to ~b at small ~b at depths less than 100 km to make the relative velocity between melt and solid non-zero even when ~b is negligibly small because of the low average temperature in the lithosphere (see below).

1 (

2.1. Material Properties

g(p)=-~

The viscosity r/ of the convecting material depends exponentially on temperature T as (1) where */0 and E are constants taken to be 1.4 x 1024pa s and 6.2 x 10 -3 °C-1, respectively. The values imply that the viscosity becomes O(102°)pa s, a typical value of the viscosity in the upper mantle (Peltier, 1985) at T = 1300°C, a typical temperature of the uppermost mantle (Turcotte and Schubert, 1982). The chemical potential for a solid-phase /~s depends on temperature T and pressure p as ~7 = 770 exp( - E T ) ,

lz'=Cpr+pV,-

T ( S o + C p log r ) ,

(2)

where C,, is specific heat 174 J mo1-1 K -1, a value for olivine calculated from Dulong-Petit's law, S o is reference entropy, and V, is molar volume of solid; C,,, V~ and, as a result, /** are assumed to be common to both of the end-members A and B. The pressure p here and in the following part of Section 2.1 is approximated as a hydrostatic pressure calculated from a reference density P0 = 3300 kg m -3 (see below). The value of V, is varied in the range between 4.24 × 10 -5 m 3 mol-1, a value close to the molar volume of forsterite at zero-pressure, and 2.12 × 10 -5 m 3 mo1-1, a value close to the molar volume of ilmenite phase of MgSiO 3 at the pressure of 650 km-seismic discontinuity (Ohtani, 1983). Boussinesq approximation is employed here and the contribution of thermal expansion to /z" is neglected. The chemical potential of the end-member i (i = A, B) for a melt-phase/,{ is defined as Ix~ = CpT+ Ah + [pV~ + Ah g ( p ) ]

- r So+C, logr+r?

where Ah is latent heat of melting at zero-pressure 9.2 × 104 J mo1-1, the value for fayalite (Robie et al., 1979), R is the gas constant 8.3 J mo1-1 K -1, T,.°, is melting temperature of the end-member i at zero-pressure, and ~[ is the content of the end-member i in melt-phase, i.e. sr~ = 1 - ~ . (Hereafter, ~ is denoted by ~l.) The function g ( p ) depends on pressure p as

nloge , (3)

AV0 - AV® ) pAV®+ 1 / A + l / p '

(4)

where AV0 and AV® are differences in molar volume between melt-phase and solid-phase at the surface and at great depth, respectively, and h is a constant. The functional form of g ( p ) and the values of T~° and h are determined so as to make the solidus temperature calculated from (2) to (4) (see below) approximate the solidus temperature of mantle material well in the pressure range less than 24 GPa (Takahashi, 1986; Ohtani and Sawamoto, 1987; Zhang and Herzberg, 1994). Here we subtracted the effect of adiabatic decompression from the experimentally determined solidus temperature since Boussinesq approximation is employed and, hence, all of the temperatures should be regarded as potential temperatures in the present formulation. The temperatures T° and T° are taken to be 1827°C and 1227°C, respectively, A is taken to be 6.3 GPa, and AV0 and AV® are taken to be 8.27 x 10 -6 m 3 mo1-1 and 2.12 x 10 -7 m 3 mo1-1, respectively. Because of the assumption of chemical equilibrium, /~t. =/z'

(i = A , B)

(5)

holds in a partially molten region. From (2), (3) and (5), the eutectic composition ~e, liquidus t e m p e r a t u r e T l at ~:/> ~:e, and solidus temperature Ts become, exp( A h / R T ° ) ~e = e x p ( A h / R T o ) + e x p ( A h / R T O ) ,

(6)

1 +g(p) T, = T° 1 - ( R T ° / A h ) log ~,'

(7)

and

L = r°[1 + g ( p ) ] ,

(8)

191

M. Kameyama et al. / Physics of the Earth and Planetary Interiors 94 (1996) 187-215

respectively, where

the density of a partially molten material p becomes

Ah

to=

R l o g [ e x p ( A h / R T ° ) + exp(Ah/RT°n)] "

(9) From the given parameter values, we obtain T° --1204°C and Ee = 0.1. (Since Et is always greater than Ee in the present calculations, we do not need to define T t at Et < E,.) The internal energy for the solid-phase e, and the melt-phase e I derived from the chemical potentials defined in (2) and (3) are, e s = CpT

(10)

and

p =po{1 - a T + / 3 ( 1

-~b--~f[1 + / 3 ( I - E , ) ]

,

(15)

where ~b is melt fraction and E = (1 - ~b)Es + ~El

(16)

is bulk composition. The density difference between melt and the coexisting matrix Ap becomes A p =__Ps -- Pl

e t = CpT + Ah[l + g ( p ) ] - p A V ,

(11)

respectively. Here AV is a difference in molar volume between the melt-phase and the solidphase and is calculated from (2), (3) and a thermodynamic relation as,

dg(p) AV= Ah

av0 - av=

d------p- AV= + (1 + p / A ) 2"

(12)

The density of the solid-phase p, and the melt-phase Pt of the binary eutectic material depend on temperature T, El (the content of the end-member A in the melt-phase) and Es (the content of the end-member A in the solid-phase)

as, Ps =Po[ 1 - a T + / 3 ( 1 - E,)]

(13)

= P0(/3(E/- Es) +

AV

[1+/3(1- Et)]

) (17)

Since AV decreases with increasing depth as indicated in (12), Ap becomes negative at great depths for a pair of melt and solid with a sufficiently large compositional contrast. When V~= 4.24 × 10-5 m3 mol-1, Et = Ee and Es = 0.64 (the value given as the initial condition in the present numerical experiments), Ap becomes negative at p > 10 GPa or at depths greater than 300 km. The depth range of greater than 300 km is consistent with the estimate of the depth range for the density inversion in the Earth's mantle suggested by Agee and Walker (1993).

and PI

2.2. Basic equations

=po(1-ar+/3(1- El) g

[l + / 3 ( 1 - E,)] ,

The momentum equation for the convective motion of solid is (14)

where P0 is a reference density 3300 kg m -3, a is thermal expansivity 3 × 10 -5 °C -x, and /3 is a constant. The value of/3 is assumed to be 0.067 at depths greater than 40 km to make Ps equal to the density of eclogite 3500 kg m -3 when E~ is equal to the eutectic composition (--basaltic composition) Ee; /3 = - 0 . 1 is assumed at depths less than 40 km to take the effect of basalt-eclogite transition into account. From (13) and (14),

-pg2-Vp+

V'[rl(vv+tvv)]

=O,

(18)

where the superscript t means transpose of a matrix, V is the velocity of solid, g is gravitational acceleration 9.8 m s -z and £ is the unit vector in z-direction; z-axis is chosen to be the vertical axis with the origin at the bottom of the box and positive upward. The continuity equation is V. V= - V . [~b(v- V)], where v is the velocity of melt.

(19)

192

M. Kameyama et al. / Physics of the Earth and Planetary Interiors 94 (1996) 187-215

We assumed that the relative motion between a melt and the coexisting solid is driven by the density difference between the melt and the solid;

Shaw, 1980; Sleep, 1988). The coefficient k4,/l~4) is assumed to depend on 4) as kg [max(4), 4)e)] 2

k4,

k4,

v - V= -----7(1 - 4))Apg2, /z(p

(20)

where/z is viscosity of melt and k4, is 'effective' permeability. The effect of shear deformation of matrix on magma-migration (McKenzie, 1984) is neglected. The effect of shear deformation of matrix is likely to be negligible except where a very strong shear deformation of matrix occurs, e.g. in a region very close to ridge axis and in the wedge mantle beneath island-arcs (Spiegelman and McKenzie, 1987; Scott and Stevenson, 1989). It is not well constrained how the effective permeability k4, depends on 4) in the Earth. At depths greater than 100 km, we assume that k4, at a sufficiently small 4), i.e. 4) < 4)c where 4)~ is a threshold (see below), is proportional to 4)3 on the basis of a study of a permeable flow of magma in porous media (McKenzie, 1984); the coefficient k4,/lz4) in (20) is assumed to depend on 4) as k4,

(22)

where 4)e is an 'effective porosity' for magmamigration by crack propagation. We assume that 4)e > 0.035 to keep the characteristic time for magma-migration through the lithosphere less than 106 years (McKenzie, 1984). When 4)> 4)c, the coefficient k4,/iz4) is assumed to be independent of 4), since the matrix disintegrates and the relative motion between solid and melt follows Stokes' law at large 4) (Abe, 1993a; Solomatov and Stevenson, 1993); the relative velocity is still proportional to the density difference between melt and solid but does not directly depend on 4). Here, we assume the coefficient k4,/Iz4) at a sufficiently large 4) to be k4,

at 4) > 4)¢.

(23)

We assume 4)c = 0.4. The energy equation is written as 0e +

k,~ 4)2 -

at 6 < 6c,

at 4) < 4)~,

(21)

where 4)o is a reference melt fraction arbitrarily chosen to be 0.05, and k~ is reference permeability. The value of k~/tz appropriate for mantle material is not well known since k~/tz strongly depends on grain size, which is poorly constrained in the Earth's mantle (McKenzie, 1984). Here, the factor k°/tz is varied from 8 6 × 10 -16 4, m 3 s kg-~ to 8.6 × 10-15 m 3 s kg-1. The value of k~/tz is within the range reasonable for the Earth's mantle, though slightly lower than the value suggested by McKenzie (1984). At depths less than 100 km, we estimate the effective permeability based on the observation that magma migrates through cracks and dykes formed by intrusion of magma, i.e. Iv - V I 4: 0, in the lithosphere even when the average temperature is not high enough to cause partial melting and hence 4) is negligibly small (Weertman, 1971;

--+V.(e+V) 0t

if+

' (v = V " ( kTVT) + q - V" 4)-V-fs

-

AV 4) p0Vz

- v)

+ v • (Ke.dVe+),

] (24)

where k v is thermal conductivity 2.5 W m - l K - l, q is the rate of internal heating due to heat-producing elements, Kead is an 'effective eddy diffusivity' (see below), and e ÷=

(1 - 4))es + 4)e +

,

(25)

where

e[ = e ~ + p A V = C p T + A h ( l + g ( p ) ) .

(26)

In (24), the effect of viscous dissipation is ignored since viscous dissipation causes a significant heating only in a highly localized region (Balachandar et al., 1995), while we are interested in the overall

M. Kameyama et al. /Physics o f the Earth and Planetary Interiors 94 (1996) 187-215

thermo-chemical state of the upper mantle here. The third term in the right-hand side of (24) represents the effect of advection by magmamigration, the fourth term represents the effect of gravitational potential energy released upon magma-migration, and the fifth term represents the effect of mixing by the convection of melt in a largely molten region. A mixing by the convection of melt is important when the melt fraction 4) exceeds the threshold 4)c and hence the framework of matrix disintegrates (Abe, 1993b). We model the effect of the convective mixing by diffusion with an effective eddy diffusivity tea d. The estimate of the characteristic time of convective mixing by Hoffman and McKenzie (1985) suggests that the effective diffusivity red a is much larger than thermal diffusivity. We carried out numerical experiments with red d an order of magnitude larger than the thermal diffusivity in most cases. In a region of 4) < ~bc, the value of red a is assumed to gradually increase with increasing 4) to keep numerical stability; we assumed that the effective diffusivity red a depends on melt fraction ~b as ° [min(~b'~be)] 3, red d = red d

(27)

~1)c

where red d ° is the value of eddy diffusivity at 6>6~. The mass transport equation is + v. (ev)

-at

=

- v)]

+ V • (readV~).

(28)

The first term in the right-hand side of (28) represents the effect of advection by magmamigration and the second term represents the effect of the convective mixing in a largely molten region. The internal heating rate q changes with time as,

aq --

+ V. (qV)

at

1 =

Aq exp(-t/r)

r q~ +

Aq exp(-t/r)

q

- V" [dpqt(v - V)] + V ' ( r e a a V q ) .

(29)

193

Here r is a decay time of heat-producing elements taken to be 1.5 Gy, an average of the decay times of 4°K and 235U, q® and Aq are constants (see below), and qt is the internal heating rate in the melt calculated from q as D ql = ( D - 1)4) + 1 q'

(30)

where D is the partition coefficient of heat-producing elements between the melt and the coexisting matrix and is taken to be 100, a value close to the partition coefficient for uranium between melt and pyroxene (Anderson, 1989), in most cases in the present calculations. The first term in the right-hand side of (29) represents the effect of the decay of heat-producing elements. If q is initially uniform and no partial melting and hence no partitioning of heat-producing elements occurs, the solution to (29) becomes q = q= + Aq e x p ( - - t / r ) .

(31)

Namely, q= + Aq is the initial value of internal heating rate and qoo is the internal heating rate at t sufficiently larger than r; q® represents the contribution of the long-lived radioactive elements. The value of q® is determined to be 4.7 × 10 -8 W m -3 from the conclusion that a significant portion of the present surface heat flow (about 60 mW m -2) comes from the internal heat source (for a review, see Christensen, (1985)). Namely, we assumed :q~dz

= 30 mW m -2.

(32)

In (32), the contribution from the second term in the right-hand side of (31) is ignored. The initial value of the internal heating rate q~o+ Aq is several times larger than q® in the Earth (Turcotte and Schubert, 1982) and is taken to be 1.9 × 10 -7 W m-3 in the present calculations. The boundary conditions are (i) impermeable and shear-stress-free conditions at all of the boundaries, (ii) fixed temperature at the top surface-boundary z = d, where d is the height of the box, (iii) adiabatic conditions at the side walls x = O, wd, where w is the aspect ratio of the box, and (iv) uniform temperature equal to the tem-

194

M. Kameyama et al. / Physicsof the Earth and Planetary Interiors 94 (1996) 187-215

perature of the heat bath Thb at the bottomboundary z = 0. Namely,

avx

T-- T~rrc,

aT ax and

ax

az = V ~ = O a t z = d '

(33)

l/x=0 at x = 0 , wd,

(34)

OVx

z=0at z=0.

T=Thb,~=V

(35)

We also assume k,

= Ked d

=

0 at z = 0, d.

(36)

The temperature at the surface Lrfc is assumed to be 27°C. The temperature of the heat bath Thb is calculated from the energy budget of the heat bath as dThb

k

wd

Eq. (37) means that the temperature of the heat bath decreases with time owing to the cooling by the convection in the box. A volume of the heat bath Vhb is larger than the volume of the convecting box by a factor of 2, the ratio of the volume of the lower mantle to the volume of the upper mantle in the Earth. The heat bath does not contain heat-producing elements. The initial conditions for (28) and (29) are uniform composition and uniform internal heating over the entire box, i.e. = ~:0 = 0.64 and q = q0 = q® + Aq,

(38)

respectively. The initial condition for (24) is e += Cp min(T~,

2.3. Numerical method

The basic equations are discretized by a finite difference scheme. A uniform mesh with 78 (vertical direction)× 92 (horizontal direction) mesh points is employed. The flux-corrected transport (FCT) method (Book et al., 1981) with Zalesak's fully two-dimensional flux-limiter (Zalesak, 1979) and Euler-leap-frog time stepping (Barr and Ashurst, 1982) are used to evaluate the contribution of advection by solid-state convection in (28) and (29) as well as the contribution of transport of latent heat by solid-state convection included in the second term in the left-hand side of (24). A third-order interpolation scheme and donner cell scheme with a zeroth order diffusion are used to evaluate the higher and lower order flux in the FCT method, respectively. An upwind scheme called power-law scheme (Patankar, 1980) with Euler time stepping is used to evaluate the contributions of heat transport by solid-state convection and thermal diffusion in (24). Donner cell method is used to evaluate the contribution of advection by migrating melt in (24), (28) and (29). Typically, 500000 time steps are necessary for each run. The momentum equation (18) and the continuity equation (19) are solved for the primitive variables V and p by a direct solver developed in Paper 1. The reliability of the employed numerical code has been already verified by the benchmark tests described in Paper 1. Although the actual calculations are carried out in a nondimensional form, the results are presented here in a dimensional form. The conversion to dimensional quantities is carried out with a length scale of d = 650 km, time scale of d2/KT = 22.4 Gy (IfT= KT~s/C p 6 × 10 -7 m 2 s -1 is thermal diffusivity), and temperature scale of A T = A h / C v = 529 K. =

2173, Tsrfc + 1 × 104× ~ ) .

(39)

Eq. (39) means that ~b is initially zero over the entire box and that the initial temperature at depth is equal to 1900°C (-- 2173 K) and linearly increases with increasing depth near the top surface-boundary. A small perturbation is added to the initial temperature distribution to start convection. The initial temperature in the heat bath is 1900°C.

3. Results

We carried out numerical experiments at various values of the five parameters, molar volume of solid-phase V~, reference permeability k,~, effective porosity in the uppermost part of the box ~be, distribution coefficient of heat-producing ele-

M. Kameyama et al. / Physics of the Earth and Planetary Interiors 94 (1996) 187-215

ments D, and eddy diffusivity in a largely molten region Kea a.° The adopted values of the parameters are shown in Table 1. 3.1. Control experiment We carried out the numerical experiment of Case 1 as the control experiment in this paper. The results of Case 1 are shown in Figs. 1 to 3. In Fig. 1, we show the plots of (a) temperature averaged over the entire box Tar (solid line) and temperature in the heat bath Thb (dotted line), (b) horizontally averaged heat flux across the top surface-boundary h, (c) root-mean-square velocity over the entire box Vrms, and (d) magma-eruption rate rh defined by 1 ,wa rh= w--dJ0 ~bvz dx at z = z c,

(40)

where z c is the height of basalt-eclogite transition versus time; rh is the flux of melt that crosses the depth level of basalt-eelogite transition per unit time. The temperatures Thb and Tar in Fig. l(a) gradually decrease with time during the period of the calculation. Magmatism episodically occurs as can be seen from the spikes in Fig. l(d). The activity of magmatism measured by the height of the spikes is high during the period of t < 1.4 Gy, while negligibly low during the subsequent period of t > 1.4 Gy. There is a sudden transition from a period of active magmatism to a period of inactive magrnatism at t = 1.4 Cry. Spikes occur in the curves of h versus time and Vrms versus time at the times of magmatism too. The root-mean-

195

square velocity V~m~ remains almost constant when magmatism does not occur. Namely, the vigor of the convection in the box remains nearly constant in spite of the cooling of the box shown in Fig. l(a) and the resulting increase in the viscosity of the convecting fluid. Fig. 2 shows the snapshots of (a) isotherms, (b) distribution of internal heating rate q, (c) distribution of ¢, and (d) distribution of melt-content ~b for Case 1 at the elapsed times shown by the numbers in the figure in the unit of Gy, while Fig. 3 shows the plots of horizontally averaged temperature T h against height at several elapsed times. At t = 0.42 Gy, a chemically stratified structure develops in the box as shown in Figs. 2(b) and 2(c). A thin crust of magrnatic products characterized by their low value of ¢ develops along the top surface boundary because of the low density of magrnatic products at depths less than 40 km. (Note that B in (15) is negative in the depth range.) Heat-producing elements strongly concentrate into the crust. Beneath the crust, a thick layer of high-¢ material with ¢ greater than ¢0 = 0.64 (the initial value of ~:) develops as indicated by the shaded area in Fig. 2(c). The high-~: material is a residue of magma and occurs in the upper part of the box because of the chemically induced positive buoyancy of the material (see (15)). The high-¢ material is rather depleted in heat-producing elements. A thick layer of low-¢ material with ¢ < ¢0, which is a mixture of magmatic products and undifferentiated material, develops in the lower part of the box because of the

Table 1 The values of parameters V~ (molar volume of solid-phase), k~/iz (the ratio of the reference permeability at depth to the viscosity of melt), eke (the effective porosity in the uppermost part of the box), K°aa (the eddy diffusivity in a largely molten region) and D (partition coefficient of heat-producing elements between melt and solid)

Case 1 Case 2 Case 3 Case 4 Case 5 Case 6 Case 7 Case 8

Vs (10 -5 m 3 mo1-1)

k~/I.t (10 -15 m 3 s kg -1)

eke

ICed dO (10 -6 m 2 s -1)

D

4.24 2.12 2.83 4.24 4.24 4.24 4.24 4.24

2.9 2.9 2.9 2.9 2.9 2.9 8.6 0.86

0.05 0.05 0.05 0.05 0.1 0.05 0.035 0.05

6 6 6 6 6 3 6 6

100 100 100 1 100 100 100 100

196

M. Kameyama et al. /'Physics of the Earth and Planetary Interiors 94 (1996) 187-215 l

2000

I

I

I

L

I

(a) t e m p e r a t u r e

..................................Thb

t.,,) I-,

Tav

1200 300

, '

, '

, '

I

I

I

I

I

i

I

f

(b) surface heat flux E 1-

0 5

I

I

d

I

i

I

I

i

I

I

I

I

I

I

I

I

(c) r m s v e l o c i t y E ol

E o 0.5

I

I

I

I

I

I

I

I

,

i

~

I

,

I

(d) m a g m a

J

e r u p t i o n rate

E O

"E 0.0 0

2

I

f

i

4

6

8

10

time [Gy] Fig. 1. Plots of (a) temperature averaged over the entire box Tar (solid line) and in the heat bath Thb (dotted line), (b) horizontally averaged heat flux across the top surface boundary h, (c) root-mean-square velocity V,, s, and (d) magma eruption rate rh, versus time for Case 1. chemically induced negative buoyancy of the material. The t e m p e r a t u r e distribution at t = 0.42 Gy is strongly superadiabatic over the entire box as can be seen from Figs. 2(a) and 3. The horizontally averaged t e m p e r a t u r e in Fig. 3 increases with increasing depth even beneath the top cold thermal boundary layer, the depth range shown by the arrow in Fig. 3. The superadiabatic temperature gradient beneath the top cold thermal boundary layer is a result of the chemical stratification. Because of the superadiabatic tempera-

ture gradient, the t e m p e r a t u r e is higher than 1800"C at the base of the box. The hot low-~ material at depth cannot uprise to the surface because of the chemically induced negative buoyancy of the material. Several patches of partially molten materials or diapirs are observed in Fig. 2(d) for t = 0.42 Gy. The melt migrates downward in the diapir at depth, indicated by the arrow in the figure, because of the density inversion between the melt and the coexisting solid described in Subsection 2.1. In the diapirs formed

M. Kameyama et al. / Physics of the Earth and Planetary Interiors 94 (1996) 187-215

at s h a l l o w e r d e p t h levels, t h e m e l t s m i g r a t e u p w a r d s owing to t h e positive b u o y a n c y i n d u c e d by melting. ( A m o r e d e t a i l e d d e s c r i p t i o n o f t h e dia p i r s will b e given in a s e p a r a t e p a p e r . ) The chemically stratified structure further dev e l o p s b e c a u s e o f active m a g m a t i s m a n d t h e t e m p e r a t u r e r e m a i n s high at d e p t h in t h e box d u r i n g t h e s u b s e q u e n t p e r i o d o f t = 0.59 G y to t = 1.21 Gy. Fig. 2 for t = 0.59 G y to 0.71 G y shows how t h e c h e m i c a l l y s t r a t i f i e d s t r u c t u r e f u r t h e r develo p s owing to m a g m a t i s m . A n u p w e l l i n g flow ind u c e d as a r e t u r n flow by t h e cold sinking p l u m e i n d i c a t e d by t h e a r r o w in t h e Fig. 2(a) for t = 0.59

(a) temperature

0.42

0.625

(b) internal heating rate

197

G y causes a p r e s s u r e - r e l e a s e m e l t i n g as i n d i c a t e d by t h e arrows in Fig. 2(d) for t = 0.62 Gy. T h e positive b u o y a n c y d u e to m e l t i n g e n h a n c e s t h e u p w e l l i n g flow a n d c a u s e s f u r t h e r p r e s s u r e - r e l e a s e m e l t i n g d u r i n g t h e p e r i o d o f t = 0.625 G y to 0.634 Gy. A t t = 0.630 Gy, a m e l t p o c k e t is f o r m e d a n d a s t e e p vertical t e m p e r a t u r e g r a d i e n t occurs along t h e t o p s u r f a c e b o u n d a r y owing to a magma-migration through the lithosphere des c r i b e d in (22) as p o i n t e d o u t by t h e arrows in Figs. 2(a) a n d 2(d). A l a r g e b l o c k o f high-~ residu a l m a t e r i a l f o r m s as a result o f t h e c h e m i c a l d i f f e r e n t i a t i o n d u e to t h e p r e s s u r e - r e l e a s e m e l t -

(c) composition

(d) melt content

f

3

'~

Fig. 2. Snapshots of (a) isotherms, (b) the contours of constant internal heating rate q, (c) iso-s¢ lines where ~ is the content of the end-member A, and (d) the contours of constant melt-content ~b for Case 1. The numbers in the figure indicate the elapsed time in the unit of Gy. In (a), the contour interval is 100*C and dashed lines show the isotherms of 1000°C and 15000C. In (b), the contour interval is 3.1 x 10 -8 W m - 3 . The regions occupied by a material with 1.4 x 10 - 7 W m - 3 >q > 4.7 X 10 -s W m - 3 a r e lightly shaded, while the regions with q > 1.4 x 10 -7 W m -3 are heavily shaded. In (c), the contour interval is 0.05. The regions occupied by a material with 0.8 > sc > 0.64 are lightly shaded, where 0.64 is the initial value of ~, while the regions with ~ > 0.8 are heavily shaded. In (d), the contour interval is 0.02. The regions with 0.2 > ~b> 0.05 are lightly shaded, while the regions with 4) > 0.2 are heavily shaded.

198

M. Kameyama et aL / Physics of the Earth and Planetary Interiors 94 (1996) 187-215

ing as indicated by the arrows in Fig. 2(c) for t = 0.71 Gy. Though the cold sinking plume causes a mixing as indicated by the arrows in Fig. 2(c) for t = 0.59 Gy to 0.66 Gy, the amount of the new addition of high-~ material with g > 0.8 is larger than the amount of high-~ material eroded by mixing and, as a whole, the amount of high-~ material with ~ > 0.8 increases during the period of t = 0.59 Gy to 0.71 Gy, i.e. the chemically stratified structure further develops owing to magrnatism. The average temperature at depth slightly drops during the period of t = 0.59 Gy to 0.66 Gy owing to the upward heat transport due to the magma-migration but soon recovers from the effect of the magmatism (see Fig. 2(a) for t = 0.59 Gy to 0.88 Gy) since the elapsed time of t = 0.59 Gy to 0.88 Gy is much shorter than the decay time of heat-producing elements 1.5 Gy and the internal heating rate is still high. Because

(a) temperature

(b) internal heating rate

of the high temperature at depth, a thin layer of partially molten material gradually develops along the bottom boundary with time during the period of t = 0.59 Gy to 1.21 Gy as shown in Fig. 2(d). The melt cannot uprise but stays along the bottom-boundary because of the density inversion. The temperature at depth, however, gradually decrease during the subsequent period of t = 1.21 Gy to 1.30 Gy (see Fig. 2(a)). The uprising diapirs indicated by the arrows in Fig. 2(d) for t = 1.23 Gy and 1.27 Gy transport heat upward and cool the deeper part of the box. The internal heating is not strong enough to keep the deeper part of the box hot any more because of the decay of the heat-producing elements. The chemically stratified structure in the box is catastrophically destroyed during the subsequent short period of t = 1.30 Gy to 1.44 Gy. Diapirs formed at depth uprise to the surface as indi-

(c) composition

(el) melt content

0.630

0.66

0.71 A

Fig. 2 (continued).

M. Kameyama et al./ Physics of the Earth and Planetary Interiors94 (1996) 187-215

cated by the arrows in Fig. 2(d) for t = 1.38 Gy and cause a mixing of low-~ material at depth and high-~ material at the shallower depth levels. Though the magmatism due to the diapirs supplies a new addition of residue and magmatic products as before, the amount of the new addition is not large enough for the chemically stratified structure to survive the mixing; the temperature at depth is no longer high enough to cause a sufficiently vigorous magmatism. The thermal state in the box also significantly changes by t = 1.44 Gy. Both the horizontally averaged temperature and its vertical gradient significantly drop in the height range of z = 100 km to 400 km during the short period as can be seen from the curves of t = 1.21 Gy and 1.44 Gy in Fig. 3 and the regions indicated by I and J in Fig. 2(a) become rather isothermal. The temperature dis-

(a) temperature

(b) internal heating rate

199

tribution shown in Fig. 2(a) for t = 1.44 Gy is characterized by a cold thermal boundary layer along the top boundary, the cold sinking plumes indicated by the arrows in the figure, a hot thermal boundary layer along the bottom boundary, and the isothermal cores I and J, and hence is rather close to a temperature distribution expected in a box where a strictly thermal convection occurs. In spite of the sudden drop in temperature at depth, the horizontally averaged temperature remains nearly constant with time in the top cold thermal boundary layer as can be seen from Fig. 3. By t--1.68 Gy, the box becomes almost chemically homogeneous except in the regions around the crust and the layer of partially molten material along the bottom boundary. Because of the change in the thermal state in the box during the period of t = 1.30 Gy to 1.44 Gy,

(c) composition

0.88 ~

1.23

Fig. 2 (continued).

(d) melt content

200

M. Kameyama et al. / Physics of the Earth and Planetary Interiors 94 (1996) 187-215

the material in the partially molten layer along the bottom boundary is cooled and almost completely solidifies by t = 2.26 Gy. During the subsequent period of t > 3.09 Gy, the low-¢ material in the layer along the bottom boundary is gradually entrained into the shallower part of the box by convection and the box becomes more chemically homogeneous. By t = 9.75 Gy, the box becomes almost completely chemically homogeneous except in the region around the crust. The horizontally averaged geotherm shown in Fig. 3 remains nearly constant with time at height z > 300 km, but gradually decreases with time at z < 300 km during the period. The sudden change in the thermal state in the box at t = 1.4 Gy causes the sudden transition from a period of active magmatism to a period of inactive magmatism observed in Fig. l(d). A large

(a) temperature

(b) internal heating rate

amount of melt is generated at depth and active magmatism occurs during the period of t < 1.4 Gy, because the temperature in the height range of z < 400 km (depths greater than 250 kin) is sufficiently high. The temperature at depth, however, suddenly drops at t = 1.4 Gy. Hot uprising diapirs do not occur at great depth any more and magmatism becomes inactive during the period of t > 1.4 Gy. In summary, the thermo-chemical evolution in the box occurs in two stages in Case 1. In the earlier stage when the internal heating rate is high, an episodic magmatism occurs actively and a chemically stratified structure develops well in the box. The temperature distribution is strongly superadiabatic over the entire box including the region beneath the top cold thermal boundary layer and the temperature at the base of the box is as high as the solidus temoerature. The chemi-

(c) composition

1.27

1.30

1.38

1.44

Fig. 2 (continued)..

(d) melt content

M. Karneyama et al. / Physics of the Earth and Planetary Interiors 94 (1996) 187-215

cally stratified structure is, however, suddenly destroyed and the temperature at depth significantly drops within a short period of time at the beginning of the later stage when the internal heating becomes sufficiently weak owing to the decay of heat-producing elements. The weak internal heating can no longer maintain the high temperature at depth in the box and, hence, can no longer keep magmatism active enough to produce the large amount of residue and magmatic products necessary to make the chemically stratified structure survive a convective mixing. After the sudden change in thermo-chemical state, magmatism becomes inactive, the box becomes chemically homogeneous, and the temperature distribution becomes adiabatic beneath the top cold thermal boundary layer with a rather low temperature at depth. The overall features of thermo-chemical evolution described above are the same as the overall

(a) temperature

(b) internal heating rate

201

features of thermo-chemical evolution obtained in Papers 1 and 2 in spite of the ad hoc modeling of magmatism adopted in the papers. The episodic magmatism observed in Fig. l(d) is also consistent with the results of numerical experiments by Tackley and Stevenson (1993). The result that a chemically stratified structure develops well owing the a mantle magmatism during the period of t < 1.4 Gy does not contradict the results of the earlier numerical experiments: the present ridge volcanism does not cause a significant chemical stratification in the mantle (Gurnis and Davies, 1986; Richards and Davies, 1989; Christensen and Hofmann, 1994). The generation rate of magma and residue obtained in the present numerical model of magrnatism at t < 1.4 Gy is much larger than the generation rate of magma and residue assumed in the earlier numerical experiments because the temperature at depth in the box at t < 1.4 Gy obtained here is

(c) composition

(d) melt content

0 1.68

2.26

ill

3.09 I

9.75

Fig. 2 (continued).

M. Kameyama et al. /Physics of the Earth and Planetary Interiors 94 (1996) 187-215

202

much higher than the geotherm assumed in the earlier numerical models; the generation rate of magma and residue strongly depends on geotherm (White and McKenzie, 1989). The difference in the generation rate of magma and residue causes the difference in degree of chemical stratification in the mantle. Indeed, the weak magmatism during the period of t > 1.4 Gy obtained here (see Fig. l(d)) does not cause a significant chemical stratification as can be seen from Fig. 2(c). We further carried out the calculations of Cases 2 and 3 at smaller ~ , the molar volume of solid phase, than the value assumed in Case 1 (see Table 1) to study how the thermo-chemical evolution in the box depends on the magnitude of positive buoyancy induced by melting. The positive buoyancy induced by melting becomes stronger as ~ becomes smaller as can be seen from (15). We found that the thermo-chemical evolution in the box occurred in two stages in Cases 2 and 3 too, i.e. that the magnitude of the positive buoyancy due to melting did not influence the overall features of the thermo-chemical evolution. We, however, found that the detail of

i

the thermo-chemical evolution depended on the magnitude of the positive buoyancy; (i) the chemical stratification in the first stage was more conspicuous in Cases 2 and 3 than in Case 1 and (ii) the partially molten layer along the bottom boundary observed in Fig. 2(d) did not occur in Cases 2 and 3.

3.2. Effects of distribution of heat-producing elements We carried out the calculations of Cases 4 to 7 (see Table 1) to understand how the thermochemical evolution obtained in Cases 1 to 3 is influenced by the heterogenous distribution of heat-producing elements and, in particular, the strong concentration of heat-producing elements in the crust observed in Cases 1 to 3 (see Figs. 2(b) and 7 below). In Cases 1 to 3, heat-producing elements concentrate in the crust owing to a fractionation of heat-producing elements into melts and the upward migration of the melts. To remove the concentration of heat-producing elements in the crust and make their distribution

i

600

",

\\\

"'" \\ E'

400

....

i~)

--

I~ .

.

.

.

.

.

.

.

solidus

"

t=

.,1', ',

0.42

II

t=,.,,

200 .....

),,,", ,,

t = 1.44

!', I

0

I

500

1000

1500

~\

, t

2000

temperature [C] Fig. 3. Plots of horizontally averaged temperature against height at several elapsed times shown by the numbers in the figure in the units of Gy for Case 1.

M. Kameyama et aL /Physics of the Earth and Planetary Interiors 94 (1996) 187-215

homogeneous, we set the partition coefficient D = 1 in Case 4. In Cases 5, 6, and 7, in contrast, we kept D at the value assumed in Cases 1 to 3, and instead we changed the effective porosity in the uppermost part of the box the, the eddy diffusivity in a largely molten region Kedd, 0 and the reference permeability k~, respectively, from the values assumed in Cases 1 to 3 to see how the thermo-chemical evolution occurs when heat-producing elements concentrate in the crust even more strongly; as will be described below, larger values of the, k~ and a smaller value of Ked 0 d lead to a stronger concentration of heat-producing elements in the crust. Figs. 4 to 6 show the result of Case 5 where the effective porosity in the uppermost part of the box the is twice larger, and hence the 'effective permeability' calculated from (22) with th = the is eight times larger, than the value assumed in Case 1; Fig. 5 shows the plot of the average temperature in the middle part of the box Trod defined as wd

f Trod=

2d/3

dxf 0

d/3 wd

2000

I

I

2d/3

I

........................................

I

I

I

(a) temperature

........................................................................... Thb

.O I--

Tav

1200

300

i

i

i

i

~

i

J

I

i

I

I

I

I

I

(b) surface heat flux ~ • E " 0 5

r

i

I

I

I

i

i

i

i

I

I

I

I

(c) rms velocity l= "~"

>

0 0.5

i

r

i

I

i

~ I

i

i

i

I

(d) magma eruption rate

~'~

dz

I

203

z, t) '

(41)

f o dx f d/3 dz versus time. All of Tar and Thb in Fig. 4(a) and T,, d in Fig. 5 for Case 5 decrease more rapidly with time than Tar, Thb and Trod for Case 1. In particular, Tmd for Case 5 rapidly decreases with time during the period of t < 0.7 Gy while Tma for Case 1 remains nearly constant during the period of t < 1.4 Gy as indicated by the dashed lines in Fig. 5. The magmatism in Case 5 is much less active than the magmatism in Case 1 throughout the period of the calculation. The spikes in Fig. 4(d) are much lower than the spikes in Fig. l(d) and the rather enhanced magrnatism observed at the beginning in Fig. 4(d) continues for only 0.8 Gy, much shorter than the period when an active magmatism occurs in Case 1. Fig. 6 shows that only a mild chemical stratification occurs in Case 5. At t = 0.44 Gy, a thin crust of low-~ material develops along the top surface boundary, a layer of high-~ material de-

.o.

.

0.0 0

2

4 time [Gy]

Fig. 4. The same as Fig. 1 but for Case 5 where the effective porosity in the uppermost part of the box ~be is twice larger than ~be in Case 1.

velops in the upper part of the box beneath the crust, and a layer of low-g material develops in the lower part of the box. The thickness of the layer of high-~ material with ~ > 0.8 is almost the same as the thickness of the layer of high-~ material with ~ > 0.8 at t = 0.42 Gy in Case 1 (see Fig. 2(c)). The temperature distribution shown in Fig. 6(a) for t -- 0.44 Gy is similar to the temperature distribution shown in Fig. 2(a) for t = 0.42 Gy too. The chemically stratified structure, however, does not further develop and, instead, is gradually smeared out by convection in Case 5 in contrast to Case 1. The thickness of the layer of high-~ material with ~ > 0.8 gradually decreases

204

M. Kameyama et aL / Physics of the Earth and Planetary Interiors 94 (1996) 187-215

o

1700

'

'

'

'

'

'

,_

' Case 1

E I,-

1400 1700

. '

.

.

.

.

~

.

.

'

'

'

!

=..,O'-'~E s

e

5

I.1400

0

t

I

I

2

4 time [Gy]

6

8

Fig. 5. Plots of the temperature T,, a averaged over a height range of ~d < z < ~d, where d is the height of the box, against time for the cases shown in the figure.

(a) temperature

0.61 ~

(b) internal heating rate

(c) composition

(d) melt content

:

0.82

1.02

Fig. 6. T h e same as Fig. 2 but for Case 5.

I I

M. Kameyama et al. / Physics of the Earth and Planetary Interiors 94 (1996) 187-215

with time during the entire period of t > 0.44 Gy in Case 5, while significantly increases with time during the period of 0.42 Gy < t < 1.4 Gy in Case 1. The temperature in the depth range shown by the arrows in Fig. 6(a) for t = 0.61 Gy to 1.26 Gy is significantly lower in Case 5 than in Case 1 during the period of 0.59 Gy < t < 1.21 Gy (see Fig. 5 too). Because of the lower temperature in Case 5, the amount of melt generated at depth is much smaller in Case 5 than in Case 1 throughout the period of the calculation as can be seen from Fig. 6(d). To understand why the temperature in the box is low, magmatism is inactive, and a chemically stratified structure does not develop well when t~e is large in Case 5, we plot the horizontally averaged internal heating rate at t = 0.8 Gy against height z in Fig. 7(a) and the degree of concentration of heat-producing elements in the

(a) temperature

(b) internal heating rate

205

crust rq against time in Fig. 7(b) for Cases 1 and 5; rq is defined by

rq= l O O x

wd

d

~a

a

fo dXfz~dzq(x'z't) [%]. fo dxfodzq(x'z't)

(42)

Fig. 7 shows that heat-producing elements concentrate much more strongly in the crust and the amount of heat-producing elements retained at depth is smaller in Case 5 than in Case 1. The value of rq for Case 5 is about 55% and is significantly larger than rq in Case 1, about 30 to 40%, during the period of t < 3 Gy. Because of the smaller amount of heat-producing elements retained at depth, the temperature in the deeper part of the box is lower (see Figs. 5 and 6(a)) and hence magmatism is milder in Case 5 than in

(c) composition

(d) melt content

1.26

2.10

3.70

6.81

ii Fig. 6 (continued).

il)

206

M. Kameyama et al. / Physics of the Earth and Planetary Interiors 94 (1996) 187-215

(a)

650

0

i

600

10

20

30

50

t = 0.8 Gy

..... ................

.1----1-

40

I

I

I

i

-E ,.¢

400

Case 1 Case 4 Case 5

i

4) J¢ 200

o

;" ..... , 1

0

,

,

,'',

2

3

4

5

internal heating rate [10 .7 W/m 3] (b) 60

]

/-~-j

I

j

~

I

_

~

I

I

]

~

I

Case 1

o~

0 60

,,

p

,.j,

I

I

I

I

f

I

I

I

I

I

I

I

I

I

Case 5 o~

0 0

2'

/4

time [Gy]

6'

8

M. Kameyama et aL /Physics of the Earth and Planetary Interiors 94 (1996) 187-215

Case 1. The milder magmatism, i.e. the lower rate of chemical differentiation, is the reason why a chemically stratified structure does not develop well in Case 5. We further carried out the calculation of Case 6 where Kea a° is smaller than the value assumed in Case 1 and the calculation of Case 7 where k~ is larger than the value assumed in Case 1 and obtained results consistent with the result of Case 5. (In Case 7, we also changed the value of 4)e to keep the 'effective permeability' in the uppermost part of the box equal to the 'effective permeability' in Case 1.) We found (i) that heat-producing elements more strongly concentrated in the crust in Cases 6 and 7 than in Case 1, though not so strongly as in Case 5, (ii) that a chemically stratified structure did not develop so well and magmatism was not so active in Cases 6 and 7 as in Case 1, (iii) that the chemically stratified structures in Case 6 and 7 were more conspicuous than the chemically stratified structure observed in Case 5, and (iv) that the magrnatisms in Cases 6 and 7 were more active than the magrnatism in Case 5. We also carried out the calculation of Case 4 where D, the partition coefficient of heat-producing elements, is equal to 1. Since the heat-producing elements homogeneously distribute and do not concentrate in the crust when D = 1, the internal heating source at depth is much stronger in Case 4 than in Case 1 as shown in Fig. 7(a). We found that the strong internal heating at depth induced an active magrnatism in Case 4 for a period much longer than in Case 1, and that a well-developed chemically stratified structure persisted for more than 10 Gy, i.e. a chemically stratified structure developed more strongly in Case 4 than in Case 1 in accordance with the results of Cases 1 and 5 to 7. In summary, the thermo-chemical evolution in the box is strongly influenced by the amount of heat-producing elements that remain at depth in

2000

I

I

I

207

i

I

I

I

(a) temperature

...........................................................

Th b ..........

0 I--

1200 300

i

i

J

i

i

i

i

I

I

I

I

I

I

I

(b) surface heat flux

R"

E E 0 5

i

i

i

i

i

f

i

I

I

I

I

I

I

I

(c) rms velocity E o

>

0 0.5

i

I

i

i

I

,

I

i

i

l

I

,

(d) magma eruption rate E o

"E 0.0 0

,JLtlA, 2

4 time [Gy]

Fig. 8. The same as Fig. 1 but for Case 8 where the reference permeability at depth k~ is smaller than k~ in Case 1 by a factor of 3.

spite of their upward migration into the crust due to melt-migration. The thermo-chemical evolution described in Section 3.1 occurs when the concentration of heat-producing elements in the crust is rather mild and a significant amount of heat-producing elements remain at depth. When heat-producing elements strongly concentrate in the crust, the box rapidly cools, magmatism, which

Fig. 7, (a) Plot of the horizontally averaged internal heating rate q against height at t = 0.8 Gy and (b) plot of the ratio rq of the amount of heat-producing elements contained in the crust to the amount contained in the entire box against time for the cases shown in the figure. The dotted lines in (b) show rq for Case 4, where heat-producing elements are uniformly distributed throughout the period of the calculation. In (a), both the vertical and horizontal scales are changed at the height of 600 km.

208

M. Kameyama et aL / Physics of the Earth and Planetary Interiors 94 (1996) 187-215

actively occurs at the beginning, soon subsides, and a chemically stratified structure does not develop well in the box. The strong influence of the degree of concentration of heat-producing elements in the crust on the thermo-chemical evolution in the box is not observed in the numerical models of Papers 1 and 2. Only the thermochemical evolution described in Section 3.1, occurs in Papers 1 and 2 because of the rather crude modeling of mantle magmatism in the earlier papers. 3.3. Effects o f slow magma-migration

In case 8 we studied how the box thermally and chemically evolves when the permeability is much lower than the permeability assumed in Case 1. We reduced the reference permeability k,~ by a factor of three. The results are shown in Figs. 8 and 9.

(a) temperature

(b) internal heating rate

The overall features of Fig. 8 and the same as the overall features of Fig. 1. Both Thb and T~ gradually decrease with time. The activity of magmatism is high at the beginning (during the period of t < 1.8 Gy) but becomes negligibly low later (during the period of t > 1.8 Gy). The vigor of convection measured by the root-mean-square velocity Vr,,,= shown in Fig. 8(c) remains almost constant when magmatism does not occur in spite of the decreasing Ta~ and Thb. There is, however, a difference between Case 8 and Case 1 in the detail of the activity of magmatism. The magmatism in Case 8 during the period of t < 1.8 Gy is less active than the magmatism in Case 1 during the corresponding period of t < 1.4 Gy as can be seen from the height of the spikes in Figs. l(d) and 8(d). The less active magmatism in Case 8 is a direct result of the lower permeability which reduces the relative velocity between melt and the coexisting matrix in Case 8.

(c) composition

(d) melt content

L

0.48

0.61

1.09

Fig. 9. The same as Fig. 2 but for Case 8.

M. Kameyama et al. / Physics of the Earth and Planetary Interiors 94 (1996) 187-215

(b)internalheating (c)composition

(a) temperature

209

(d)melcontent t

-j

rate

t

1.56

2.30~

~

I

~

1

Fig. 9 (continued).

20

L_

I

I

I

I

I

I

Case I

>

0 20

I

~ o ~ Case 8

0 ~ 0

2

,"~ ~" ~ ' ~ ~ 4 6 time [Gy]

8

Fig. 10. Plots of the ratio r v of the relative velocity between the melt and the coexisting solid ( J v z - Vz J) averaged over the partially molten region with d, > 0.01 to the root-mean-square velocity Vr,~s against time for the cases shown in the figure.

210

M. Kameyama et al. / Physics of the Earth and Planetary Interiors 94 (1996) 187-215

Fig. 9 shows that a chemically stratified structure does not develop so well in Case 8 as in Case 1. At the beginning, the thermo-chemical state in the box in Case 8 is similar to the thermo-chemical state observed in Case 1 as can be seen from Fig. 9 for t = 0.48 Gy. The temperature distribution is superadiabatic over the entire box and a chemically stratified structure in the box develops with a crust of low-~ material along the top surface boundary, a layer of high-~ material in the upper part, and a layer of low-~: material in the deeper part. Heat-producing elements concentrate in the crust. The chemically stratified structure, however, does not further develop during the subsequent period of t > 0.48 Gy as can be seen from the thickness of the layer of high-~ material with ~ > 0.8 in Fig. 9(c). The thickness gradually decreases with time during the entire period of t > 0.48 Gy in Case 8 in spite of the repeated occurrence of pressure-release melting in such big upwelling diapirs as the diapirs indicated by the arrows in Fig. 9(d) for t = 0.48 Gy and 1.09 Gy. The thickness of the layer of low-~ material at depth also decreases with time during the entire period of t > 0.48 Gy. To understand why a chemically stratified structure does not develop so well in Case 8 as in Case 1, we present the plot of the ratio r v --- ( I vz - ~ [)/Vrms, where ( [ v z - V z 1) is the relative velocity between a melt and the coexisting solid averaged over the partially molten region with ~b > 0.01, versus time for both of the cases in Fig. 10. The ratio r~ is a measure of the importance of chemical differentiation by magmatism in the thermo-chemical evolution in the box. As can be seen from the figure, r v in Case 8 is smaller than rv in Case 1. The value of r~ fluctuates around 2 in Case 8, while rv fluctuates around 7 during the period of t < 1.4 Gy in Case 1. Namely, the influence of chemical differentiation on the thermo-chemical evolution in the box is milder in Case 8 than in Case 1. The milder influence of chemical differentiation is the reason why a chemically stratified structure does not develop so well in Case 8 as in Case 1. The value of r~ around 2 in Case 8 suggests that the parameter values assumed in Case 8 are close to the regime of "circulation m o d e " (Scott,

1988) where the effect of mixing by the convection of matrix is stronger than the effect of chemical differentiation by melt-migration. The circulation mode is, however, not likely to be relevant to the mantle magmatism in the Earth. The relative velocity between a melt and the coexisting matrix is most likely to be much larger than the velocity of solid-state convection in the Earth's mantle (McKenzie, 1984).

4. Discussion A presumption in the present numerical models is that mantle convection occurs separately in the upper mantle and the lower mantle. A layered convection occurs in a chemically homogeneous mantle owing to the phase-transition between y-phase and post-spinel-phase of (Mg, Fe)2SiO 4 if the viscosity in the mantle is sufficiently low and, as a result, the Rayleigh number is higher than a threshold on the order of 108 (e.g. Christensen and Yuen, 1985; Steinbach et al., 1993). The condition is most likely to be satisfied in the mantle of the early Earth, which is probably much hotter and much less viscous than the mantle of the present Earth (Steinbach et al., 1993; Yuen et al., 1993). The chemical stratification in the upper mantle obtained in the present numerical experiments (see, for example, Fig. 2(c)) probably enhances the propensity toward a layered convection. The magma generated in the Earth's mantle is more enriched in basaltic components than the primitive mantle material and the magmatic products formed upon cooling and solidification of the magma are highly buoyant in the uppermost part of the lower mantle, i.e. the magmatic products are unlikely to sink deep into the lower mantle (Irifune and Ringwood, 1993; O'Neill and Jeanloz, 1994). The layer of melt along the 650 km boundary observed in some of the present numerical calculations (see Figs. 2(d) and 9(d)) is also likely to enhance the propensity toward a layered convection. Even if a melt formed in the Earth's mantle is denser than the coexisting solid in the deeper part of the upper mantle (Stolper et al., 1981; Agee and Walker,

M. Kameyama et al. / Physics of the Earth and Planetary Interiors 94 (1996) 187-215

1993), the melt is less dense than the coexisting solid in the upper part of the lower mantle (Miller et al., 1991) and does not sink deep into the lower mantle. A layered convection is likely to have occurred in the early Earth owing to the phase transition between y-phase and post-spinel-phase of (Mg, Fe)2SiO 4, the chemical stratification, and the melting at depth. It is not, however, clear if a layered convection occurs in the present mantle where the viscosity is not low enough for the phase transition to completely suppress a convective flow across the phase boundary at the depth of 650 km by itself (Christensen and Yuen, 1985; Machetel and Weber, 1991; Tackley et al., 1993, 1994; Steinbach et al., 1993; Solheim and Peltier, 1994). We apply our numerical results only to the hot upper mantle expected in the early Earth. Our numerical results suggest two possible overall pictures of the thermo-chemical regime in the upper mantle of the early Earth depending upon how strongly heat-producing elements concentrate in the crust and hence how strongly the deeper part of the upper mantle becomes depleted in heat-producing elements owing to a mantle magmatism. If the concentration of heat-producing elements in the crust is mild and the upper mantle is not so strongly depleted in heat-producing elements, our numerical results suggest that the overall picture of t he thermo-chemical regime in the upper mantle of the early Earth becomes as shown in Fig. 2 for t < 1.4 Gy. The upper mantle is chemically stratified with a layer of olivine-rich residual material in the upper part and a layer of material enriched in pyroxene and garnet in the deeper part. The temperature distribution in the upper mantle is strongly superadiabatic beneath the top thermal boundary layer, i.e. lithosphere, as well as within the lithosphere. The temperature in the deeper part of the upper mantle is as high as the solidus temperature and hence is much higher than the present geotherm. The lithosphere in the early Earth is, in contrast, as cold as the present lithosphere. The superadiabatic temperature gradient is maintained by the chemical stratification. Magmatism is induced by upwellings of the hot materials in the deeper part of the upper mantle, and becomes episodic owing

211

to the interaction of magmatism and mantle convection. Our numerical results also suggest that the regime of hot and chemically stratified upper mantle catastrophically ends when the internal heating rate becomes sufficiently low due to a decay of heat-producing elements. The chemically stratified structure is suddenly destroyed and the temperature at depth and the activity of magmatism suddenly drop at the time of the catastrophic end. The overall picture is the same as the picture suggested in Papers 1 and 2. The overall picture suggested from Fig. 2 is consistent with many observations from the Archean continental crusts. The rather cold lithosphere suggested above is consistent with the observations of diamond inclusions (Boyd and Gurney, 1986) and metamorphic rocks (Bickle, 1978); the observations suggest a geotherm as low as the geotherm in the present continental lithosphere beneath the shield area in the Archean continental lithosphere. The very high geotherm at depth obtained in the numerical model is consistent with the observations of komatiites which indicate the occurrence of regions much hotter than the present upper mantle at depth in the Archean upper mantle (e.g. Bickle et al., 1977; Ohtani, 1984; Takahashi, 1986). The episodic nature of magrnatism shown in Fig. l(d) is consistent with the observation that magmatism occurs episodically in the Archean (Kr/Sner, 1985). The catastrophic end of the regime of hot and chemically stratified upper mantle is consistent with an observation of MgO-content in komatiites and picrites of various ages from the Archean to the present (Campbell and Griffiths, 1992). Campbell and Griffiths (1992) plotted the MgO-content in komatiites and picrites against the age of the rocks and found that the maximum of the MgOcontent obtained at a given age suddenly decreases with decreasing age at the end of the Archean, suggesting that the temperature in the source region of the magma significantly dropped within a short period of time at the end of the Archean. The sudden decrease in temperature suggested from the drop in MgO-content is consistent with the sudden decrease in the horizontally averaged temperature at depth at the time of the catastrophic end shown in Figs. 3 and 5.

212

M. Kameyama et al. / Physics of the Earth and Planetary Interiors 94 (1996) 187-215

If heat-producing elements strongly concentrate in the crust and the upper mantle is strongly depleted in heat-producing elements owing to mantle magmatism, our numerical result of Case 5 suggests that a large portion of the upper mantle rapidly cools with a characteristic time of several hundred million years in the earliest stage of the Earth's history, though the temperature at the base of the upper mantle does not decrease so rapidly (see Figs. 4(a) and 6). Magmatism, which actively occurs at the beginning, soon subsides owing to the rapid cooling. A chemically stratified structure develops at the beginning but soon is destroyed, and the upper mantle becomes more chemically homogeneous except at the base of the upper mantle. A rapid cooling of the upper mantle in the earliest stage of the Earth's history has been suggested in the earlier works on the thermal history of the Earth (Sharpe and Peltier, 1979; e.g. Schubert et al., 1979; Davies, 1980). An overall picture of the thermo-chemical regime in the Archean upper mantle suggested from the rapid cooling by earlier workers is that the average temperature in the upper mantle is rather low in the Archean and that hot plumes with a large excess temperature uprise through the rather cold upper mantle (Jarvis and Campbell, 1983; Campbell and Griffiths, 1992; Nisbet et al., 1993). Komatiitic magma is generated in the hot uprising plumes, whereas the rather low geotherm in the continental lithosphere in the Archean suggested from diamond inclusions (Boyd and Gurney, 1986) and metamorphic rocks (Bickle, 1978) is a result of the rather cold upper mantle (Jarvis and Campbell, 1983; Campbell and Griffiths, 1992; Nisbet el al., 1993). The episodic nature of the magmatism in the Archean (KrOner, 1985) may be due to solitary waves in the hot uprising plumes (Schubert et al., 1989) or to the interaction between hot uprising plumes and high-pressure induced phase transitions in the transition layer (Liu et al., 1991). The sudden decrease in the MgO-content in the magmatic products of mantle magmatism at the end of the Archean may imply that the temperature in the source region of the hot uprising plumes, i.e. the hot thermal boundary layer at the bottom, suddenly dropped at the

end of the Archean (Campbell and Griffiths, 1992). Our numerical result of Case 5 where the upper mantle cools rapidly is not, however, consistent with the observations of komatiites. In Case 5, the temperature in the partially molten parts of hot plumes or diapirs that eventually uprise to the surface is only 1500°C at most (see Fig. 6(a)) and is not high enough to produce a komatiitic magma during the period of t > 0.7 Gy, i.e. after the rapid cooling at the earliest stage of the thermal evolution shown in Fig. 5. (Note that the temperature at the bottomboundary Thb shown in Fig. 4(a) does not significantly decrease from the initial value of 1900°C till t = 3 Gy and hence the presumption of layered convection is probably safe during the period of t < 3 Gy even after the rapid cooling.) Furthermore, the diapirs do not extend to a sufficiently deeper part of the box to produce a komatiitic magma as shown in Fig. 6(d); both high temperature and high pressure are necessary for the generation of komatiitic magma (Ohtani, 1984; Takahashi, 1986). Namely, we could not obtain a numerical model where heat-producing elements strongly concentrate in the crust, the box rapidly cools, and still hot plumes or diapirs with an excess temperature large enough to produce komatiitic magma uprise from the deeper part of the box to the surface after the rapid cooling. The absence of hot uprising plumes with a large excess temperature in Case 5 is due to the chemical stratification that is formed by the active magmatism at the beginning of the calculation shown in Fig. 4(d) and persists at the base of the box even after the rapid cooling (see Fig. 6(c)). The chemically induced negative buoyancy of the low-sc material at the base of the box cancels out a large portion of the thermally induced positive buoyancy in the hot thermal boundary layer along the bottom boundary and suppresses hot uprising plumes with a large excess temperature to grow from the hot thermal boundary layer. An active mantle magmatism has probably occurred at the earliest stage of the Earth's history too, because of the hot origin of the Earth (e.g. Kaula, 1979; Stevenson, 1981). A layer of dense magmatic products is, therefore,

M. Kameyama et al. /Physics of the Earth and Planetary Interiors 94 (1996) 187-215

most likely to have formed at the base of the upper mantle in the early Earth to suppress hot uprising plumes with a large excess temperature even if heat-producing elements have strongly concentrated into the crust, the upper mantle has rapidly cooled, and the mantle magmatism has rapidly subsided in the early Earth. Because the earlier stage of the thermo-chemical evolution observed in Cases 1 to 3 meshes with the observations of Archean continental crust while the numerical model of Case 5 does not, the regime of hot and chemically stratified upper mantle together with the catastrophic end suggested from Cases 1 to 3 is more likely as the thermo-chemical regime in the Archean upper mantle.

Acknowledgment We thank the two anonymous reviewers for helpful comments. This work was supported by grant-in-aid for Scientific Research (No. 06222221 and No. 06640548) sponsored by the Ministry of Education, Science, and Culture of Japan. The calculations are carried out on a supercomputer SX3R at the computer center of Tohoku University.

References Abe, Y., 1993a. Thermal evolution and chemical differentiation of the terrestrial magma ocean. In: Evolution of the Earth and Planets, Geophys. Monogr., 74, IUGG Vol. 14, I U G G / A G U , pp. 41-54. Abe, Y., 1993b. Physical state of the very early Earth. Lithos, 30: 223-235. Agee, C.B. and Walker, D., 1993. Olivine flotation in mantle melt. Earth Planet. Sci. Lett., 114: 315-324. Anderson, D.L., 1989. Theory of the Earth, Blackwell, Boston, MA, 366 pp. Balachandar, S., Yuen, D.A., Reuteler, D.M. and Lauer, G.S., 1995. Viscous dissipation in three-dimensional convection with temperature-dependent viscosity. Science, 267: 11501153. Barr, P.IC and Ashurst, W.T., 1982. Evaluation of Zalesak's flux-corrected-transport algorithm for convection and diffusion. Sandia Report SAND81-8223. Sandia National Laboratories, Albuquerque. Bickle, M.J., 1978. Heat loss from the Earth: a constraint on

213

Archaean tectonics from the relation between geothermal gradients and the rate of plate production. Earth Planet. Sci. Lett., 40: 301-315. Bickle, M.J., Ford, C.E. and Nisbet, E.G., 1977. The petrogenesis of peridotitic komatiites: evidence from high-pressure melting experiments. Earth Planet. Sci. Lett., 37: 97-106. Book, D.L., Boris, J.P. and Zalesak, S.T, 1981. Flux-corrected transport. In: D.L. Book (Editor), Finite-Difference Techniques for vectorized Fluid Dynamics Calculations. Springer, New York, pp. 29-55. Boyd, F.R. and Gurney, J.J., 1986. Diamonds and the African lithosphere. Science, 232: 472-477. Campbell, I.H. and Griffiths, R.W., 1992. The changing nature of mantle hotspots through time: implications for the chemical evolution of the mantle. J. Geol., 92: 497-523. Christensen, U.R., 1985. Thermal evolution models for the Earth. J. Geophys. Res., 90: 2995-3007. Christensen, U.R. and Hofmann, A.W., 1994. Segregation of subducted oceanic crust in the convecting mantle. J. Geophys. Res., 99: 19867-19884. Christensen, U.R. and Yuen, D.A., 1985. Layered convection induced by phase transitions. J. Geophys. Res., 90: 1029110300. Davies, G.F., 1979. Thickness and thermal history of continental crust and root zones. Earth Planet. Sci. Lett., 44: 231-238. Davies, G.F., 1980. Thermal histories of convective Earth models and constraints on radiogenic heat production in the Earth. J. Geophys. Res., 85: 2517-2530. Gurnis, M. and Davies, G.F., 1986. Mixing in numerical models of mantle convection incorporating plate kinematics. J. Geophys. Res., 91: 6375-6395. Hoffman, N.R.A. and McKenzie, D.P., 1985. The destruction of geochemical heterogeneities by differential fluid motions during mantle convection. Geophys. J. R. Astron. Soc., 82: 163-206. Honda, S., 1995. A simple parameterized model of Earth's thermal history with the transition from layered to whole mantle convection. Earth Planet. Sci. Lett., 131: 357-369. Irifune, T. and Ringwood, A.E., 1993. Phase transformations in subducted oceanic crust and buoyancy relationships at depths of 600-800 km in the mantle. Earth Planet. Sci. Lett., 117: 101-110. Jarvis, G.T. and Campbell, I.H., 1983. Archean komatiites and geotherms: solution to an apparent contradiction. Geophys. Res. Lett., 10: 1133-1136. Jeanloz, R. and Morris, S., 1986. Temperature distribution in the crust and mantle. Annu. Rev. Earth Planet Sci., 14: 377-415. Jordan, T.H., 1975. The continental tectosphere. Rev. Geophys. Space Phys., 13: 1-12. Kaula, W.M., 1979. Thermal evolution of Earth and Moon growing by planetesimal impacts. J. Geophys. Res., 84: 999-1008. Kr6ner, A., 1985. Evolution of the Archean continental crust. Annu. Rev. Earth Planet. Sci., 13: 49-74.

214

M. Kameyama et al. / Physics of the Earth and Planetary Interiors 94 (1996) 187-215

Liu, M., Yuen, D.A., Zhao, W. and Honda, S., 1991. Development of diapiric structures in the upper mantle due to phase transitions. Science, 252: 1836-1839. Machetel, P. and Weber, P., 1991. Intermittent layered conveetion in a model mantle with an endothermic phase change at 670 km. Nature, 350: 55-57. McKenzie, D., 1984. The generation and compaction of partially molten rock. J. Petrol., 25: 713-765. McKenzie, D. and Richter, F.M., 1981. Parameterized thermal convection in a layered region and the thermal history of the Earth. J. Geophys. Res., 86: 11667-11680. Miller, G.H., Stopler, E.M. and Ahrens, T.J., 1991. The equation of state of a molten komatiite 2: application to komatiite petrogenesis and the Hadean mantle. J. Geophys. Res., 96:11849-11864. Nisbet, E.G. and Fowler, C.M.R., 1983. Model for Archean plate tectonics. Geology, 11: 376-379. Nisbet, E.G., Cheadle, M.J., Arndt, N.T. and Bickle, M.J., 1993. Constraining the potential temperature of the Archaean mantle: a review of the evidence from komatiites. Lithos., 30: 291-307. Ogawa, M., 1988. Numerical Experiments on coupled magmatism-mantle convection system: implications for mantle evolution and Archean continental crusts. J. Geophys. Res., 93:15119-15134. Ogawa, M., 1993. A numerical model of a coupled magmatism-mantle convection system in Venus and the Earth's mantle beneath Archean continental crusts. Icarus, 102: 40-61. Ogawa, M., 1994. Effects of chemical fractionation of heatproducing elements on mantle evolution inferred from a numerical model of coupled magmatism-mantle convection system. Phys. Earth Planet. Inter., 83: 101-127. Ohtani, E., 1983. Melting temperature distribution and fractionation in the lower mantle. Phys. Earth Planet. Inter., 33: 12-25. Ohtani, E., 1984. Generation of komatiite magma and gravitational differentiation in the deep upper mantle. Earth Planet. Sci. Lett., 67: 261-272. Ohtani, E. and Sawamoto, H., 1987. Melting experiment on a model chondritic mantle composition at 25 GPa. Geophys. Res. Lett., 14: 733-736. O'Neill, B. and Jeanloz, R., 1994. MgSiO3-FeSiO3-AI203 in the Earth's lower mantle: perovskite and garnet and 1200 km depth. J. Geophys. Res., 99: 19901-19915. Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow. Hemisphere, Washington, DC, 197 pp. Peltier, W.R., 1985. Mantle convection and viscoelasticity. Annu. Rev. Fluid Mech., 17: 561-608. Richards, M.A. and Davies, G.F., 1989. On the separation of relatively buoyant components from subducted lithosphere. J. Geophys. Res., 16: 831-834. Richardson, S.H., Gurney, J.J., Erlank, A.J. and Harris, J.W., 1984. Origin of diamonds in old enriched mantle. Nature, 310: 198-202. Richter, F.M., 1984. Regionalized models for the thermal

evolution of the Earth. Earth Planet. SCi. Lett., 68: 471484. Richter, F.M., 1985. Models for the Archean thermal regime. Earth Planet. Sci. Lett., 73: 350-360. Robie, R.A., Hemingway, B.S. and Fisher, J.R. 1979. Thermodynamic properties of minerals and related substances at 298.15 K and 1 Bar (10 s Pascals) pressure and at higher temperatures. Geological Survey Bulletin 1452, United States Government Printing Office, Washington. Schubert, G., Cassen, P. and Young, R.E., 1979. Subsolidus convective cooling histories of terrestrial planets. Icarus, 38: 192-211. Schubert, G., Oison, P. Anderson, C. and Goldman, P., 1989. Solitary waves in mantle plumes. J. Geophys. Res., 94: 9523-9532. Scott, D.R., 1988. The competition between percolation and circulation in a deformable porous medium. J. Geophys. Res., 93: 6451-6462. Scott, D.R. and Stevenson, D.J., 1984. Magma solitons. Geophys. Res. Lett., 11: 1161-1164. Scott, D.R. and Stevenson, D.J., 1986. Magma ascent by porous flow. J. Geophys. Res., 91: 9283-9296. Scott, D.R. and Stevenson, D.J., 1989. A self-consistent model of melting, magma migration and buoyancy-driven circulation beneath mid-ocean ridges. J. Geophys. Res., 94: 2973-2988. Sharpe, H.N. and Peltier, W.R., 1979. A thermal history model for the Earth with parameterized convection. Geophys. J. R. Astron. Soc., 59: 171-203. Shaw, H.R., 1980. The fracture mechanisms of magma transport from the mantle to the surface. In: R.B. Hargraves (Editor), Physics of Magmatic Processes. Princeton University Press, Princeton, NJ, pp. 201-264. Sleep, N.H., 1988. Tapping of melt by veins and dikes. J. Geophys. Res., 93: 10255-10272. Sleep, N.H. and Windley, B.F., 1982. Archean plate tectonics: constraints and inferences. J. Geol., 90: 363-379. Solheim, L.P. and Peitier, W.R., 1994. Avalanche effects in phase transition modulated thermal convection: a model of Earth's mantle. J. Geophys. Res., 99: 6997-7018. Solomatov, V.S. and Stevenson, D.J., 1993. Suspension in convective layers and style of differentiation of a terrestrial magma ocean. J. Geophys. Res., 98: 5375-5390. Spiegelman, M. and McKenzie, D., 1987. Simple 2-D models for melt extraction at mid-ocean ridges and island arcs. Earth Planet. Sci. Lett., 83: 137-152. Steinbach, V., Yuen, D.A. and Zhao, W., 1993. Instabilities from phase transitions and the timeseales of mantle thermal evolution. Geophys. Res. Lett., 20: 1119-1122. Stevenson, D.J., 1981. Models of the Earth's core. Science, 214: 611-619. Stolper, E., Walker, D., Hager, B.H. and Hays, J.F., 1981. Melt segregation from partially molten source regions: the importance of melt density and source region size. J. Geophys. Res., 86: 6261-6271. Tackley, P.J. and Stevenson, D.J., 1993. A mechanism for

M. Kameyama et al. / Physics of the Earth and Planetary Interiors 94 (1996) 187-215 spontaneous self-perpetuating volcanism on the terrestrial planets. In: D.B. Stone and S.K. Runcorn (Editors), Flow and Creep in the Solar System: Observations, Modeling and Theory. Kluwer Academic Publishers, Netherlands, pp. 307-321. Tacldey, P.J., Stevenson, D.J., Glatzmaier, G.A. and Schubert, G., 1993. Effects of an endothermic phase transition at 670 km depth in a spherical model of convection in the Earth's mantle. Nature, 361: 699-704. Tacldey, P.J., Stevenson, D.J., Glatzmaier, G.A. and Schubert, G., 1994. Effects of multiple phase transitions in a three-dimensional spherical model of convection in Earth's mantle. J. Geophys. Res., 99: 15877-15901. Takahashi, E., 1986. Melting of a dry peridotite KLB-1 up to 14 GPa: Implications on the origin of periodotitic upper mantle. J. Geophys. Res., 91: 9367-9382. Turcotte, D.L. and Schubert, G., 1982. Geodynamics: Applications of continuum physics to geological problems. Wiley, New York, 450 pp. Viaar, N.J., van Keken, P.E. and van der Berg, A.P., 1994.

215

Cooling of the Earth in the Archaean: consequences of pressure-release melting in a hotter mantle. Earth Planet. Sci. Lett., 121: 1-18. Weertman, J., 1971. Theory of water-fiUed crevasses in glaciers applied to vertical magma transport beneath oceanic ridges. J. Geophys. Res., 76: 1171-1183. White, R. and McKenzie, D., 1989. Magmatism at rift zones: the generation of volcanic continental margins and flood basalts. J. Geophys. Res., 94: 7685-7729. Yuen, D.A., Hansen, U., Zhao, W., Vincent, A.P. and Malevsky, A.V., 1993. Hard turbulent thermal convection and thermal evolution of the mantle. J. Geophys. Res., 98: 5355-5373. Zalesak, S.T., 1979. Fully multidimensional flux-corrected transport algorithms for fluids. J. Comp. Phys., 31: 335362. Zhang, J. and Herzberg, C., 1994. Melting experiments on anhydrous peridotite KLB-1 from 5.0 to 22.5 GPa. J. Geophys. Res., 99: 17729-17742.