KEYWORDS
Runaway subduction, Genesis flood, powerlaw creep, thermal
runaway, catastrophic plate tectonics
ABSTRACT
Experimental investigation of the solid state deformation
properties of silicates at high temperatures has revealed that the deformation
rate depends on the stress to a power of about 3 to 5 as well as strongly
on the temperature. This highly nonlinear behavior leads to the potential
of thermal runaway of the mantle's cold upper boundary layer as it peels
away from the surface and sinks through the hot mantle. The additional fact
that the mineral phase changes that occur at 660 km depth act as a barrier
to convective flow and lead to a tendency for large episodic avalanche events
compounds the potential for catastrophic dynamics. Twodimensional finite
element calculations are presented that attempt to model these strongly
nonlinear phenomena. It is proposed that such a runaway episode was responsible
for the Flood described in Genesis and resulted in massive global tectonic
change at the earth's surface.
INTRODUCTION
The only event in Scripture since creation capable of the
mass destruction of living organisms evident in the fossil record is the
Genesis Flood. A critical issue in any model for earth history that accepts
the Bible as accurate and true is what was the mechanism for this catastrophe
that so transformed the face of the earth in such a brief span of time.
The correct answer is crucial to understanding the Flood itself and for
interpreting the geological record in a coherent and valid manner. It is
therefore a key element in any comprehensive model of origins from a creationist
perspective. Ideas proposed as candidate mechanisms over the past century
include collapse of a water vapor canopy [5], near collision of a large
comet with the earth [12], rapid earth expansion [11], and violent rupture
of the crust by pressurized subterranean water [4]. There are serious difficulties
with each of these ideas.
Another possibility is that of runaway subduction of the
preFlood ocean lithosphere [2,3]. A compelling logical argument in favor
of this mechanism is the fact that there is presently no ocean floor on
the earth that predates the deposition of the fossiliferous strata. In other
words all the basalt that comprises the upper five kilometers or so of today's
igneous ocean crust has cooled from the molten state since sometime after
the Flood cataclysm began. The age of today's seafloor relative to the fossil
record is based on two decades of deep sea drilling and cataloging of fossils
in the sediments overlying the basalt basement by the Deep Sea Drilling
Program as well as radiometric dating of the basalts themselves [14]. Presumably,
there were oceans and ocean floor before the Flood. If this preFlood seafloor
did not subduct into the mantle, what was its fate? Where are these rocks
today? On the other hand, if the preFlood seafloor did subduct, it must
have done so very rapidly within the year of the Flood. In regard to the
fate of the preFlood seafloor, there is strong observational support in
global seismic tomography models for cold, dense material near the base
of the lower mantle in a belt surrounding the present Pacific ocean [16].
Such a spatial pattern is consistent with subduction of large areas of seafloor
at the edges of a continent configuration commonly known as Pangea.
There are good physical reasons for believing that subduction
can occur in a catastrophic fashion because of the potential for thermal
runaway in silicate rock. This mechanism was first proposed by Gruntfest
[6] in 1963 and was considered by several in the geophysics community in
the early 1970's [1]. Previous ICC papers [2,3] have discussed the process
by which a large cold, relatively more dense, volume of rock in the mantle
generates deformational heating in an envelope surrounding it, which in
turn reduces the viscosity in the envelope because of the sensitivity of
the viscosity to temperature. This decrease in viscosity in turn allows
the deformation rate in the envelope to increase, which leads to more intense
deformational heating, and finally, because of the positive feedback, results
in a sinking rate orders of magnitude higher than would occur otherwise.
It was pointed out that thermal diffusion, or conduction of heat out of
the zone of high deformation, competes with this tendency toward thermal
runaway. It was argued there is a threshold beyond which the deformational
heating is strong enough to overwhelm the thermal diffusion, and some effort
was made to characterize this threshold.
The important new aspect addressed in this paper is the
dependence of the viscosity on the deformation rate itself. Although this
deformation rate dependence of viscosity has been observed experimentally
in the laboratory for several decades, the difficulty of treating it in
numerical models has deterred most investigators from exploring many of
its implications. Results reported in the previous ICC papers did not include
this highly nonlinear phenomenon. Significant improvements in the numerical
techniques that permit large variations in viscosity over small distances
in the computational domain, however, now make such calculations practical.
The result of including this behavior in the analysis of the thermal runaway
mechanism is to discover a much stronger tendency for instability in the
earth's mantle. Moreover, deformation rates orders of magnitude higher than
before throughout large volumes of the mantle now can be credibly accounted
for in terms of this more realistic deformation law. This piece of physics
therefore represents a major advance in understanding how a global tectonic
catastrophe could transform the face of the earth on a time scale of a few
weeks in the manner that Genesis describes Noah's Flood.
Recent papers by several different investigators [10,13,18,19]
have also shown that the mineral phase changes which occur as the pressure
in the mantle increases with depth also leads to episodic dynamics. The
spinel to perovskite plus magnesiowustite transition at about 660 km depth
is endothermic and acts as a barrier to flow at this interface between the
upper and lower mantle. It therefore tends to trap cold material from the
mantle's upper boundary layer as it peels away from the surface and sinks.
Numerical studies show that, with this phase transition present, flow in
the mantle becomes very episodic in character and punctuated with brief
avalanche events that dump the cold material that has accumulated in the
upper mantle into the lower mantle. The episodic behavior occurs without the inclusion of the
physics that leads to thermal runaway. This paper argues that when temperature
and strain rate dependence of the rheology is included, the time scale for
these catastrophic episodes is further reduced by orders of magnitude. In
this light, the Flood of the Bible with its accompanying tectonic expressions
is a phenomenon that is seems to be leaping out of the recent numerical
simulations.
MATHEMATICAL FORMULATION
In this numerical model the silicate mantle is treated
as an infinite Prandtl number, anelastic fluid within a domain with isothermal,
undeformable, tractionfree boundaries. Under these approximations the following
equations describe the local fluid behavior:
Here p denotes pressure, r
density, g gravitational
acceleration, t deviatoric stress, u fluid velocity, T absolute
temperature, g the Grueneisen
parameter, k thermal conductivity, H volume heat production rate, c_{v} specific heat at constant volume, m dynamic shear viscosity, K the
isothermal bulk modulus, and a the volume coefficient of thermal expansion. The quantities
p_{r}, r_{r}, and T_{r} are, respectively, the
pressure, density, and temperature of the reference state. I
is the identity tensor. The superscript T in (4) denotes the tensor transpose.
Equation (1) expresses the conservation of momentum in the infinite Prandtl
number limit. In this limit, the deformational term is so large that the
inertial terms normally needed to describe less viscous fluids may be completely
ignored. The resulting equation (1) then represents the balance among forces
arising from pressure gradients, buoyancy, and deformation. Equation (2)
expresses the conservation of mass under the anelastic approximation. The
anelastic approximation ignores the partial derivative of density with respect
to time in the dynamics and thereby eliminates fast local density oscillations.
It allows the computational time step to be dictated by the much slower
deformational dynamics. Equation (3) expresses the conservation of energy
in terms of absolute temperature. It includes effects of transport of heat
by the flowing material, compressional heating and expansion cooling, thermal
conduction, shear or deformational heating, and local volume (e.g., radiogenic)
heating.
The expression for the deviatoric stress given by equation
(4) assumes the dynamic shear viscosity m
depends on temperature, pressure, and strain rate. The stress therefore
is nonlinear with respect to velocity, and the rheological description is
nonNewtonian. This formulation is appropriate for the deformation regime
in solids known as powerlaw creep to be discussed below. Equation (5) represents
density variations as linearly proportional to pressure and temperature
variations relative to a simple reference state of uniform density, pressure
and temperature. Parameter values used are r_{r }=3400 kg m^{}^{3}, p_{r}=0, T_{r}=1600 K, g=10 m/s, g=1, k=4 W m^{1}K^{1}, H=1.7 x 10^{8}
W m^{3}, c_{v
}=1000 J kg^{1}K^{1}, and K=1 x 10^{12 }Pa.
POWERLAW CREEP
Laboratory experiments to characterize the high temperature
solid state deformation properties of silicates have been carried out by
many investigators over the last three decades [8,9]. These experiments
have established that, for temperatures above about sixty percent of the
melting temperature and strain rates down to the smallest achievable in
the laboratory, silicate materials such as olivine deform according to a
relationship of the form [8]
=A
s^{n} exp[ (E* + pV*)/RT]

(6)

where e is the strain rate, A a material constant, s the
differential stress, n a dimensionless constant on the order of 3 to 5,
E* an activation energy, p is pressure, V* an activation volume, R the universal
gas constant, and T absolute temperature. This relationship implies that
at constant temperature and pressure the deformation rate increases dramatically
more rapidly than the stress. Because the strain rate increases as the stress
to some power greater than one, this type of deformation is known as powerlaw
creep. This relationship may also be expressed in terms of an effective
viscosity m=0.5s/e that depends on the strain rate e as [9, 17, p. 291]
m=B ^{q}
exp[(E* + pV*)/nRT]

(7)

where B=0.5A^{1/n} and q=1  1/n. A value for n of 3.5, appropriate for the
mineral olivine [8,9], yields a q of 0.714. This means that the effective
viscosity m decreases strongly as the strain rate increases. A tenfold increase in the strain rate, for
example, yields an effective viscosity, at fixed temperature and pressure,
a factor of 5.2 smaller! For a 10^{10} increase in strain rate,
the effective viscosity decreases by more than a factor of 10^{7}.
The effect is even more pronounced for larger values of n.
Fig. 1 is a deformation
mechanism map for olivine that shows the region in stresstemperature space
where powerlaw creep is observed. Note that there exists a boundary between
the powerlaw creep regime and that of diffusional creep. Because the strain
rates for diffusional creep are so smalltoo small in fact to be realized
in laboratory experimentsthis boundary is poorly constrained. Kirby [8,
p. 1461] states that the boundary may in actuality be substantially to the
left of where he has drawn it. In any case at a given temperature there
is a threshold value for the strain rate at which point one crosses from
the diffusional regimewhere the strain rate depends linearly on the stressinto
the powerlaw regime. From Fig. 1
this threshold is on the order of 10^{17} to 10^{14} s^{1} for temperatures about 60% of the melting temperature
and stresses of about 1 MPa.
Powerlaw creep is included in the numerical model simply
by using the effective viscosity given by (7) in (4), where the scalar strain
rate e
is obtained by taking the square root of the second invariant of the rate
of strain tensor d=(u
+ u^{T})/2. To remove the singularity
in (7) for zero strain rate and to model explicitly the transition between
diffusion creep and powerlaw creep, a minimum or threshold strain rate
_{o} is incorporated into the
formulation. For regions in the domain where the strain rate
exceeds _{o}, equation (7) applies.
Otherwise the viscosity is strain rate independent. The parameter B is specified
in terms of a reference viscosity m_{o} at reference temperature T_{r} and zero strain rate as
B=m_{o}/{_{o}^{q} exp[(E* + pV*)/nRT_{r}]}. To model the viscosity
contrast between the upper mantle and lower mantle, the reference viscosity
is allowed to vary with depth and increase in a linear fashion by a factor
of 50 between 400 and 700 km. For purposes of numerical stability the threshold
strain rate _{o} is assumed to vary as 1/m_{o}.
PHASE CHANGES
The jumps in seismic quantities observed at depths of about
410 km and 660 km in the earth closely match phase transitions observed
in laboratory experiments at similar temperatures and pressures for olivine
to spinel and from spinel to perovskite silicate structures, respectively.
These phase transitions that occur as the pressure increases and the crystal
structures assume more compact configurations almost certainly play a critical
role in the mantle's dynamical behavior. In a calculation in which silicate
material is transported through these depths and undergoes these phase changes,
two effects need to be taken into account. One is the latent heat released
or absorbed and the other is the deflection of the phase boundary upward
or downward. The latent heat may be accounted for by adding or removing
heat through the volume heating term in equation (3) proportional to the
vertical flux of material through the transition depth. The latent heat
per unit mass is obtained from the Clapeyron equation which expresses that
in a phase transition DH=(dp/dT) T DV, where DH is the enthalpy change, or latent heat, and DV
is the change in specific volume. The Clapeyron slope (dp/dT) is a quantity
that can be determined experimentally for a given transition. The deflection
in the location of a phase boundary occurs because the pressure, and therefore
the depth, at which the phase change occurs depends on the temperature.
The effect of such a deflection enters as a contribution to the buoyancy
term in equation (1). A downward deflection represents positive buoyancy
because the lighter phase now occupies volume normally occupied by the denser
phase. The Clapeyron slope is also a constant of proportionality in the
boundary deflection Dh=(dp/dT) DT/rg that
arises from a deviation DT from the reference temperature. The values for the Clapeyron
slope used here are 1 x 10^{6} Pa/K for the 410 km transition
and 2 x 10^{6 }Pa/K
for the 660 km transition. Note that the exothermic 410 km transition leads
to a positive or upward deflection for a cold slab and hence increased negative
buoyancy, while the endothermic 660 km transition leads to a downward deflection
and reduced negative buoyancy. The 660 km transition therefore acts to inhibit
buoyancy driven flow while the 410 km transition acts to enhance it.
NUMERICAL APPROACH
The set of equations (1)(5) is solved in a discrete manner
on a uniform rectangular mesh with velocities located at the mesh nodes
and temperatures, pressures, and densities at cell centers. Piecewise linear
finite elements are used to represent the velocity field, while the cell
centered variables are treated as piecewise constant over the cells. The
calculational procedure on each time step is first to apply a twolevel
conjugate gradient algorithm [15] to compute the velocity and pressure fields
simultaneously from Eq. (1) and (2). This task involves solving 3n simultaneous
equations for 2n velocity unknowns and n pressure unknowns, where n is the
total number of nodes in the mesh. Key to the procedure is an iterative
multigrid solver that employs an approximate inverse with a 25point stencil.
This large stencil for the approximate inverse enables the method to handle
large variations in viscosity in a stable fashion. The outstanding rate
of convergence in the multigrid solver is responsible for the method's overall
high efficiency. The piecewise linear finite element basis functions provide
secondorder spatial accuracy. The temperature field is updated according
to Eq. (3) with a forwardintime finite difference van Leer limited advection
method.
RESULTS
Two calculations will now be described that illustrate
the effects of power law creep on the stability of a sinking slab. The two
calculations are identical except for the value of the strain rate threshold
above which power law creep occurs. In the first case, the threshold _{o} in the upper 400 km is
3 x 10^{13} s^{1}
which is sufficiently large that no power law creep occurs anywhere in the
domain. In the second case, the threshold is 6.5 x 10^{14}
s^{1}, about a
factor of five smaller. In this case runaway eventually takes place. These
calculations are performed in a rectangular domain 2890 km high and 1280
km wide on a mesh with 64 x 64 cells of uniform size. The viscosity m_{o }at zero strain rate and 1600 K increases in a simple linear
fashion by a factor 50 between 400 km and 700 km depth to represent the
stiffer rheology of the earth's lower mantle compared with the upper mantle.
The phase changes at 410 km and 660 km depth are both included. The endothermic
phase transition at 660 km as well as the higher intrinsic viscosity below
this depth both act to inhibit flow from above. The calculations are initialized
with a uniform temperature of 1600 K except for a slablike anomaly 100 km
wide extending from the top to a depth of 400 km with a central temperature
of 1000 K and a thermal boundary layer at the top such that the temperature
in the topmost layer of cells is initially 1000 K. The upper boundary temperature
is fixed at 700 K and the bottom at 1600 K.
Fig. 2 shows four
snapshots in time spaced roughly 6 x 10^{6} years apart of the calculation
with the larger strain rate threshold. Plots of temperature and effective
viscosity are displayed with velocities superimposed. Note that the initial
maximum velocity drops by a factor of two as the slab encounters increasing
resistance from the higher viscosity and 660 km phase change. The colder
material tends to accumulate and thicken in width in the depth range between
400 and 700 km. When sufficient thickening of the zone of cold material
has occurred, it begins to penetrate slowly into the region below 700 km.
Fig. 3 shows the
effects of a strain rate threshold _{o} sufficiently low that power
law creep is occurring in a significant portion of the problem domain. The
first three snapshots in time for this case resemble those for the previous
case. The main difference are regions of reduced effective viscosity in
the region below 700 km evident in the first and third snapshots due to
the powerlaw rheology. A major change is observed, however, in the fourth
snapshot with an increase in peak velocity and a notable reduction in effective
viscosity below the head of the developing cold plume. In the fifth snapshot
the peak velocity has increased by another 80% and there is more than a
factor of ten reduction in the effective viscosity ahead of the plume. Also
displayed in this snapshot is the viscous heating rate that shows intense
heating surrounding the plume. In the sixth snapshot, the head of the cold
plume is preceded by a belt of high temperature, the velocity has almost
doubled again, the effective viscosities near the plume have dropped even
further, and the heating rate adjacent to the plume has more than doubled.
Shortly after this point in the calculation, runaway occurs and the computation
crashes.
DISCUSSION
What do these calculations have to say about the mantle
and the Flood? First of all, powerlaw rheology dramatically enhances the
potential for thermal runaway. Numerical calculations are not really necessary
to reach this conclusion. Equation (7) indicates an increase in the deformation
rate leads to a reduction in the effective viscosity and reinforces the
reduction in viscosity an increase in temperature provides. These effects
work together in a potent way. An exciting further consequence of the powerlaw
rheology is that high velocities and strain rates can now occur throughout
the mantle. A hint of this can be inferred from the last two snapshots in
Fig. 3. Large and increasing velocities
are not just associated with the sinking plume itself but are observed throughout
the domain. The remaining horizontal sections of the initial cold upper
boundary layer, for example, are also moving at much higher speeds.
In interpreting these numerical experiments it is important
to realize that one is attempting to explore numerically a physically unstable
process. Customary numerical difficulties associated with strong gradients
in the computed quantities are compounded when such a physical instability
occurs. The strategy is to explore the region of parameter space nearby
but not too close to where the instability actually lives. The calculation
of Fig. 3 therefore does not reveal
the true strength of the instability relative to the situation of a moderately
lower value for the threshold strain rate. It is also useful to point out
how various quantities scale relative to one another. The velocities are
inversely proportional to the reference viscosity. A tenfold reduction in
the reference viscosity gives ten times higher velocities. Similarly, the
threshold strain rate for runaway behavior is inversely proportional to
the reference viscosity since strain rate is proportional to velocity. So
reducing the reference viscosity by a factor of ten yields a threshold strain
rate for runaway ten times larger. This neglects the diminished influence
of thermal diffusion at the higher velocities.
How do the parameters used in these calculations compare
with those estimated for the earth? The values used for g, g, k, H, r_{r}, c_{v}, T_{r}, and a in eq. (1)(5) are all reasonable
to within +/30% for the simplified reference state that is employed. The
values used for the Clapeyron slopes for the phase transitions are two to
three times too small and so the effects of the phase changes are underrepresented.
The most important parameters are the reference viscosity and the threshold
strain rate for powerlaw creep. The reference viscosity leads to velocities
prior to runaway that are in accord with current observed plate velocities
of a few centimeters per year. The threshold strain rates used are within
the powerlaw creep region for olivine as given by Kirby (Fig.
1). A large uncertainty is the extrapolation of the creep behavior of
olivine to the minerals of the lower mantle for which there is essentially
no experimental data. The issue is not whether powerlaw creep occurs in
these minerals but what the stress range is in which it occurs. It is likely
the threshold strain rate is not many orders of magnitude different from
olivine. These calculations therefore seem relevant to the earth as we observe
it today.
One difficulty in making a connection between these calculations
and the Flood is their time scale. Some 2 x 10^{7
}years is needed before the instability occurs in the second
calculation. Most of this time is involved with the accumulation of a large
blob of cold, dense material at the barrier created by the phase transition
at 660 km depth. This time span disappears when the initial condition consists
of a large belt of cold material already trapped above this phase transition
in the preFlood mantle. A relatively small amount of additional negative
buoyancy in such a belt can then trigger runaway. One means for providing
a quick pulse of negative buoyancy is by the sudden conversion to spinel
of olivine in a metastable state that resides at depths below the usual
transition depth of about 410 km. Such metastability can arise because the
changes in volume and structure associated with a phase transition do not
necessarily occur spontaneously as transition conditions are reached, especially
if the material is cold. Some means of nucleation of seed crystals of the
new phase is generally required. If such nucleation does not happen, then
substantial amounts of the less dense phase can survive to depths much greater
than what the assumption of a spontaneous transition would imply. Indeed,
there is observational evidence for significant amounts of metastable olivine
in the slab currently beneath Japan [7]. A shock wave passing through such
a volume of metastable material can initiate the nucleation and cause a
sudden conversion to the denser phase. Present day deep focus earthquakes
likely represent manifestations of such a process on a small scale. In the
context of the Flood, it is conceivable that an extraterrestrial impact
of modest size could have triggered a sudden conversion of metastable material
to the denser phase and the resulting earthquakes then propagated in a selfsustaining
manner to convert the metastable material throughout much of the upper mantle
to the denser spinel phase, which in turn initiated the runaway avalanche
of upper mantle rock into the lower mantle. It is also conceivable that
a single large earthquake generated by causes internal to the earth could
have been the event that caused a sudden conversion of the metastable material
and then the runaway avalanche.
CONCLUSIONS
Rapid sinking through the mantle of portions of the mantle's
cold upper boundary facilitated by the process of thermal runaway appears
to be a genuine possibility for the earth. A highly nonlinear deformation
law for silicate minerals at conditions of high temperature known as powerlaw
creep, documented by decades of experimental effort, in which the effective
viscosity decreases strongly with the deformation rate, makes thermal runaway
almost a certainty for a significant suite of conditions. This deformation
law also makes possible strain rates consistent with large scale tectonic
change within the Biblical time frame for the Flood. Mineralogical phase
changes combined with the viscosity contrast between upper and lower mantle
conspire to provide the setting in which a sudden triggering of a runaway
avalanche of material trapped in the upper mantle into the lower mantle
can occur. Calculations by other investigators that include the endothermic
phase transition, but not temperature or strain rate dependent viscosity,
also display the tendency for episodic avalanche events [10,13,18,19]. Such
an episode of catastrophic runaway is here presented as the mechanism responsible
for Noah's Flood.
REFERENCES
[1] O. L. Anderson and P. C. Perkins, Runaway
Temperatures in the Asthenosphere Resulting from Viscous Heating, Journal of Geophysical Research, 79(1974), pp. 21362138.
[2] J. R. Baumgardner, Numerical
Simulation of the LargeScale Tectonic Changes Accompanying the Flood, Proceedings of the International Conference on Creationism,
R. E. Walsh, et al, Editors, 1987, Creation Science Fellowship, Inc.,
Pittsburgh, PA, Vol. II, pp. 1728.
[3] J. R. Baumgardner, 3D
Finite Element Simulation of the Global Tectonic Changes Accompanying
Noah's Flood, Proceedings of the Second International Conference
on Creationism, R. E. Walsh and C. L. Brooks, Editors, 1991, Creation
Science Fellowship, Inc., Pittsburgh, PA, Vol. II, pp. 3545.
[4] W. T. Brown, Jr., In the Beginning, 1989, Center
for Scientific Creation, Phoenix.
[5] J. C. Dillow, The Waters Above, 1981, Moody Press,
Chicago.
[6] I. J. Gruntfest, Thermal
Feedback in Liquid Flow; Plane Shear at Constant Stress,
Transactions of the Society of Rheology, 8(1963), pp. 195207.
[7] T. Iidaka and D. Suetsugu, Seismological Evidence for Metastable Olivine Inside
a Subducting Slab, Nature, 356(1992), pp.
593595.
[8] S. H. Kirby, Rheology
of the Lithosphere, Reviews of Geophysics
and Space Physics, 21(1983), pp. 14581487.
[9] S. H. Kirby and A. K. Kronenberg, Rheology of the Lithosphere: Selected Topics,
Reviews of Geophysics and Space Physics, 25(1987), pp. 12191244.
[10] P. Machetel and P. Weber, Intermittent Layered Convection in a Model Mantle with
an Endothermic Phase Change at 670 km, Nature,
350(1991), pp. 5557. [11] G. R. Morton, The Flood on an Expanding Earth, Creation Research Society Quarterly, 19(1983), pp.
219224.
[12] D. W. Patton, The Biblical Flood and the Ice Epoch,
1966, Pacific Meridian Publishing, Seattle.
[13] W. R. Peltier and L. P. Solheim, Mantle
Phase Transitions and Layered Chaotic Convection, Geophysical Research Letters, 19(1992), pp. 321324.
[14] Proceedings of the Ocean Drilling Program
[15] A. Ramage and A. J. Wathen,
Iterative Solution Techniques for Finite Element Discretisations of Fluid
Flow Problems, Copper Mountain Conference
on Iterative Methods Proceedings, Vol. 1., 1992.
[16] M. A. Richards and D. C. Engebretson, LargeScale
Mantle Convection and the History of Subduction, Nature, 355(1992), pp. 437440. [17] F. D. Stacey,
Physics of the Earth, 2nd ed., 1977, John Wiley & Sons, New York.
[18] P. J. Tackley, D. J. Stevenson, G. A. Glatzmaier,
and G. Schubert, Effects of an Endothermic Phase Transition at 670 km
Depth on Spherical Mantle Convection, Nature,
361(1993), pp. 699704.
[19] S. A. Weinstein, Catastrophic Overturn of the Earth's Mantle Driven by
Multiple Phase Changes and Internal Heat Generation,
Geophysical Research Letters, 20(1993), pp. 101104.

