RUNAWAY SUBDUCTION AS THE DRIVING MECHANISM FOR THE GENESIS FLOOD

JOHN R. BAUMGARDNER, Ph.D.
1965 Camino Redondo
Los Alamos, NM 87544

Presented at the Third International Conference on Creationism
Pittsburgh, PA, July 18-23, 1994


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. Two-dimensional 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 pre-Flood
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 pre-Flood seafloor did not subduct into the mantle,
what was its fate? Where are these rocks today? On the other hand,
if the pre-Flood seafloor did subduct, it must have done so very
rapidly --within the year of the Flood. In regard to the fate of the
pre-Flood 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, traction-free boundaries. Under these
approximations the following equations describe the local fluid
behavior:

0=- (p - pr) + (r - rr) g + t

(1)
0= (r u)

(2)
dT/dt=- (T u) - (g - 1) T u + [ (k T) + t : u + H]/rrcv

(3)
wheret =m [ u + ( u)T - 2 I ( u)/3]

(4)
andr=rr + rr(p - pr)/K - a(T - Tr).

(5)

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, cv specific heat at constant volume, m dynamic
shear viscosity, K the isothermal bulk modulus, and a the volume
coefficient of thermal expansion. The quantities pr, rr, and Tr 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
non-Newtonian. This formulation is appropriate for the deformation
regime in solids known as power-law 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 rr =3400 kg m-3, pr=0, Tr=1600 K, g=10
m/s, g=1, k=4 W m-1K-1, H=1.7 x 10-8 W m-3, cv =1000 J kg-1K-1, and
K=1 x 1012 Pa.


POWER-LAW 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 sn 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
power-law 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 1010 increase in strain rate, the
effective viscosity decreases by more than a factor of 107. 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 stress-temperature space where power-law creep is
observed. Note that there exists a boundary between the power-law
creep regime and that of diffusional creep. Because the strain rates
for diffusional creep are so small--too small in fact to be realized
in laboratory experiments--this 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 regime--where the
strain rate depends linearly on the stress--into the power-law
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.

Power-law 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 + uT)/2. To remove the singularity
in (7) for zero strain rate and to model explicitly the transition
between diffusion creep and power-law 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 mo at reference
temperature Tr and zero strain rate as B=mo/{o-q exp[(E* +
pV*)/nRTr]}. 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/mo.

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 106 Pa/K for the 410 km transition
and -2 x 106 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 two-level 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 25-point 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 second-order
spatial accuracy. The temperature field is updated according to Eq.
(3) with a forward-in-time 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 mo 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 106 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 power-law
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, power-law 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 power-law
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, rr, cv, Tr,
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 power-law 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 power-law
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 power-law 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 107 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 pre-Flood 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
self-sustaining 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 power-law 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. 2136-2138.

[2] J. R. Baumgardner, Numerical Simulation of the Large-Scale
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. 17-28.

[3] J. R. Baumgardner, 3-D 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. 35-45.

[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. 195-207.

[7] T. Iidaka and D. Suetsugu, Seismological Evidence for
Metastable Olivine Inside a Subducting Slab, Nature, 356(1992),
pp. 593-595.

[8] S. H. Kirby, Rheology of the Lithosphere, Reviews of
Geophysics and Space Physics, 21(1983), pp. 1458-1487.

[9] S. H. Kirby and A. K. Kronenberg, Rheology of the Lithosphere:
Selected Topics, Reviews of Geophysics and Space Physics,
25(1987), pp. 1219-1244.

[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. 55-57. [11] G. R. Morton, The Flood on an Expanding
Earth, Creation Research Society Quarterly, 19(1983), pp. 219-224.

[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. 321-324.

[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, Large-Scale Mantle
Convection and the History of Subduction, Nature, 355(1992), pp.
437-440. [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.
699-704.

[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. 101-104.