William Paley Institute
Intelligent Design



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

Creation Science Fellowship, Inc.


Any comprehensive model for earth history consistent with the data
from the Scriptures must account for the massive tectonic changes
associated with the Genesis Flood. These tectonic changes include
significant vertical motions of the continental surfaces to allow
for the deposition of up to many thousands of meters of
fossil-bearing sediments, lateral displacements of the continental
blocks themselves by thousands of kilometers, formation of all of
the present day ocean floor basement rocks by igneous processes, and
isostatic adjustments after the catastrophe that produced today's
Himalayas, Alps, Rockies, and Andes. This paper uses 3-D numerical
modeling in spherical geometry of the earth's mantle and lithosphere
to demonstrate that rapid plate tectonics driven by runaway
subduction of the pre-Flood ocean floor is able to account for this
unique pattern of large-scale tectonic change and to do so within
the Biblical time frame.


Many diverse mechanisms have been put forward to explain the
dramatic and rapid geological changes connected with the Genesis
Flood [6,7,13,14]. This event is here conceived to have generated
the portion of the geological record beginning with the initial
abrupt fossil appearance of multicellular organisms and including
all of the so-called Paleozoic and Mesozoic eras and the lower part
of the Cenozoic. In other words, the Flood is understood, in terms
of normal usage of the words in the Genesis account, to be a global
catastrophe that destroyed all the non-aquatic air-breathing life on
the earth except for that preserved in the ark. Since the Scriptures
indicate no large-scale destruction of life between the time of
creation and the Flood, it logically follows that the initial abrupt
appearance of multicellular fossils in the rock record must
represent the onset of this cataclysm. The huge amount of energy
required to accomplish such a vast amount of geological work so
quickly together with the amazing order evident in the stratigraphic
record and the smooth pattern seafloor spreading and continental
drift documented in today's ocean floor obviously impose severe
limitations on candidate mechanisms.

What constraints might one use to discriminate among possible
mechanisms for the Flood? One is the pattern of downwarping and
uplift of the earth's surface that produced the observed patterns of
sedimentation. Broadly speaking, it is possible to divide the
continental regions of today's earth into three general categories
according to the type and amount of sedimentary cover. Cratonic
shield areas such as the Canadian Shield, the African Shield, and
the Scandinavian Shield, represent regions mostly barren of
Phanerozoic, or fossil-bearing, sediment. Surface rocks are instead
pre-Phanerozoic crystalline rocks, frequently displaying strong
metamorphism and often deeply eroded. Cratonic platform areas, a
second category, represent broad regions of continental surface with
generally extensive and uniform Phanerozoic sedimentary deposits
commonly a few kilometers in thickness. The third category includes
Phanerozoic tectonic belts which frequently contain huge thicknesses
of sediments--often up to tens of kilometers--usually with strong
compressive deformations, evidence of large vertical displacements,
and vast amounts of volcanism and metamorphism. These zones are
mostly located along the margins of cratonic shield or platform
regions and usually contain high mountains.

These three categories, in the context of the Flood, respectively
represent broadly uplifted and eroded areas, broadly downwarped
areas that accumulated moderate thicknesses of sediment, and
localized belts where downwarping and deformation were extreme and
where huge thicknesses of sediment accumulated. The evidence
indicates that when the forces responsible for the extreme
downwarping in these tectonic belts abated, high mountains appeared
as the deep, narrow, sediment-filled trenches rebounded
isostatically. The sedimentary patterns therefore suggest that
transient processes, almost certainly operating in the earth's
mantle, caused dynamical subsidence and uplift within craton
interiors and intense localized downwarping at craton edges. In the
context of the Flood, these observational data speak of large and
rapid vertical motions of the earth's surface. Such vertical motions
represent distinctive patterns of internal stress and mechanical
work that must be accounted for by any successful mechanism.

A second major geological constraint concerns the large lateral
displacements of the cratonic blocks that also occurred during the
Flood. From a stress distribution standpoint this requirement of
translating continental blocks by thousands of kilometers in a short
period of time severely constrains candidate mechanisms because it
involves the solid-state deformation of rock in the mantle below.
That craton interiors display so little Phanerozoic deformation
despite the fact the cratons traversed such vast distances so
rapidly means that stress levels within the cratons never approached
the fracture or yield limits and that the forces responsible for
moving these huge bodies of rock were diffuse and relatively uniform
over the area of the block. Mechanisms that move the plates by
applying forces at their edges cannot produce this general absence
of deformation in the craton interiors. The only conceivable
mechanisms able to move plates so far and so rapidly with hardly any
internal deformation are those that involve large scale flow in the
earth's mantle and that apply relatively mild and uniform tractions
on the base of the plates. This constraint as well as the previous
one both point to catastrophic overturning of the mantle driven by
gravitational potential energy in large volumes of cold rock at the
earth's surface and/or in the upper mantle and assisted by a runaway
instability resulting from a temperature and stress dependent
deformation law for silicate rock. The thrust of this paper is to
report advances in numerical modeling of such a mechanism for the
Flood. Results from this effort have been presented in papers at the
two previous ICC meetings in 1986 and 1990 [4,5]. In the 1990 paper
it was shown how subducting ocean floor along the Pangean margins
leads to a pulling apart of the supercontinent in a manner generally
consistent with the pattern of seafloor spreading recorded in the
rocks of today's ocean floor [16]. This paper describes a number of
improvements in the model. One is the use of a much more detailed
reference state for the earth that includes compressibility and
phase changes. Another is the addition of depth variation in the
mantle's viscosity structure that provides for a low viscosity upper
mantle and a higher viscosity lower mantle. Another is a much
improved plate treatment that includes the oceans. The plates are
now tracked using a highly accurate particle-in-cell method. Dynamic
surface topography and sea level are now also computed as part of a
time dependent calculation. This yields maps of the continental
flooding that occurs in response to the mantle's internal dynamics.
In addition there are several numerical improvements that allow
larger time steps and provide increased accuracy.


The earth's mantle in the numerical model is treated as an
irrotational, infinite Prandtl number, anelastic Newtonian fluid
within a spherical shell with isothermal, undeformable,
traction-free boundaries. Under these conditions the following
equations describe the local fluid behavior:

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

0= (r u)

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

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

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

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 radially varying pressure, density, and
temperature of the reference state used for the mantle. 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 (as well as the rotational
terms) 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)

The expression for the deviatoric stress given by equation (4)
assumes a viscosity m that is dependent on the radial temperature
and pressure distribution but independent of the strain rate. The
stress therefore is linear with respect to velocity and represents
the customary description for the deformation of a Newtonian fluid.
This rheological law applies to the type of deformation in solids
known as diffusion creep that is believed to occur in the mantle
under conditions of extremely small strain rate. Equation (5)
represents density variations as linearly proportional to pressure
and temperature variations relative to a reference state. The
compressible reference state is chosen to match observational data
for the earth to a high degree of precision. It includes the density
jumps associated with mineralogical phase changes. In the numerical
model the set of equations (1)-(5) is solved for each grid point in
the computational domain during each time step.


Equations (1)-(5) represent conservation of momentum, mass, and
energy in terms of the local velocity, pressure, and temperature.
The material properties such as thermal conductivity and specific
heat may also vary with position. A much better approach than simply
assuming constant values for these quantities is to rely on a
reference model that provides these material properties as well as
reference values for the temperature, pressure, and density as a
function of depth through the mantle. Substantial effort has been
invested over the last several decades to use seismic and other
geophysical observations to formulate radial seismic earth models
[8]. Such models typically provide density and compressional and
shear wave speeds as a function of depth. It is possible, however,
to construct more comprehensive earth models that give the full
suite of thermodynamic quantities by using an equation of state
together with estimates for material properties of silicate minerals
obtained from experimental measurements. A desirable attribute of
the more comprehensive models is that they reproduce the density
profile of the seismic models.

The reference model used here is based on an equation of state that
represents the density and temperature dependence of pressure as two
independent functions, that is, p(r,T)=p1(r) + p2(T). The Morse
equation of state [2], derived from an atomic potential model of a
crystalline lattice, is employed for the density dependence and
given as follows:

p1(r)=[3Ko/(Ko' - 1)] (r/ro)2/3 E (E - 1)

whereE=exp{ (Ko' - 1) [1 - (r/ro)-1/3] }

Here ro is the uncompressed zero-temperature density, Ko is the
uncompressed zero-temperature isothermal bulk modulus, and Ko' is
the derivative of Ko with respect to pressure. These three material
parameters specify the pressure-density relationship for a given
mineral assemblage. By choosing appropriate values for the upper
mantle, the transition zone between 410 and 660 km depth, and the
lower mantle, one can match the density profile given by the seismic
models quite closely. The values used for these quantities versus
depth are given below .

Depth Range (km) ro (kg/m3) Ko (Pa) Ko'
0-410 3425 1.4 x 1011 5.00
410-510 3695 1.6 x 1011 4.00
510-660 3725 1.6 x 1011 4.00
660-2890 4220 2.6 x 1011 3.85

The thermal component assumed for the equation of state is simply
p2(T)=aKT, where a is the thermal expansivity and K is the
isothermal bulk modulus. Using standard thermodynamic relationships
together with experimentally obtained estimates for quantities such
as the specific heat and Grueneisen parameter, one can integrate the
equation of state with depth through the mantle, starting at the
earth's surface, and obtain a consistent set of thermodynamic
quantities as a function of depth. Because the gravitational
acceleration and the radial density distribution depend on each
other, it is necessary to iterate the calculation to obtain a state
that is in hydrostatic balance as well as in thermodynamic
equilibrium. Depth profiles for rr, Tr, pr, g, K, a, g, and cv
resulting from such a calculation are displayed in Fig. 1.
Temperatures chosen for the top and bottom boundaries were 300 K and
2300 K, respectively. Also shown in Fig. 1 are profiles for the
thermal conductivity and dynamic shear viscosity. The thermal
conductivity is assumed constant with depth except in the bottom
portion of the lower mantle where it is assumed to increase by about
a factor of three because of the elevated temperatures.

Depth variation in the dynamic shear viscosity is modeled using a
temperature and pressure dependent relationship of the form [9]

m=mo exp[ -(E* + prV*)/RTr] (7)

where mo is a depth independent reference viscosity, E* is an
activation energy, V* is an activation volume, and R is the
universal gas constant. As in the case for the parameters of the
Morse equation of state, separate values for E* and V* for the upper
mantle, transition zone, and lower mantle are assumed. These are as

Depth Range (km) E* (kJ/mole) V* (m3/mole)
0-410 500 10.0 x 10-6
410-660 555 6.0 x 10-6
660-2890 640 2.6 x 10-6

Due to limitations in the numerical algorithms, the extremely large
viscosities that arise in the cold upper boundary layer of the
mantle are clipped to a maximum value of 2mo and the values of E*
and V* are scaled by a factor of 0.7 relative to those given above.
The resulting profile showing the depth variation of m is displayed
in Fig. 1.


The jumps in the density profile at 410 and 660 km, respectively,
(Fig. 1) correspond in pressure and temperature to the transitions
observed experimentally between olivine and spinel and between
spinel and perovskite silicate structures. In a dynamical
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 locally 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 -4 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.


The set of equations (1)-(5) is solved in a discrete manner on a
mesh constructed from the regular icosahedron [2,3]. The mesh used
in the calculations (Fig. 2) has 10242 nodes in each of 17 radial
layers for a total of 174,114 nodes. There are 160 nodes around the
equator which implies a horizontal spatial resolution of 250 km at
the earth's surface. Nonuniform spacing of nodes in the radial
direction assists in resolving the boundary layers.

The calculational procedure on each time step is first to apply a
two-level conjugate gradient algorithm [17] to compute the velocity
and pressure fields simultaneously from Eq. (1) and (2). This task
involves solving 4n simultaneous equations for 3n 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 [2]
formulated in terms of a finite element representation of the
continuum equations. The outstanding rate of convergence in the
multigrid solver is responsible for the method's overall high
efficiency. Special piecewise linear spherical finite element basis
functions provide second-order spatial accuracy [2,3]. The
temperature field is updated according to Eq. (3) with a
forward-in-time finite difference interpolated donor cell advection
method. Tectonic plates at the earth's surface are included in this
framework by finding the unique Euler rotation vector w for each
plate such that the net torque y resulting from the surface stress
field acting over the area of the plate is zero. The surface
velocity field corresponding to this piecewise constant set of rigid
plate rotations is then applied as a surface velocity boundary
condition when solving equation (1). An iterative method is employed
to determine the rotation vectors on each time step. For a given
interior velocity field u and an estimate of w for a given plate,
small perturbations in w about the x-, y-, and z-axes are made to
compute torque sensitivities dy/dw. The current estimate for w is
improved by subtracting a correction Dw proportional to y/(dy/dw) in
a manner that drives y to zero.

Because of the need for very accurate treatment of plate boundaries,
a Lagrangian particle-in-cell method is used to define the plates
themselves and to track their motion. Four particles per node have
been found adequate to provide a sufficiently accurate plate
representation. Piecewise linear basis functions are used to map
particle data to the mesh nodes. The particles are moved in a
Lagrangian manner at each time step using the same piecewise linear
basis functions to interpolate the nodal velocities to the
particles. The advantage of this particle method is extremely low
numerical diffusion and hence the ability to minimize the smearing
of the plate edges. When oceanic plate begins to overlap another
plate, the ocean plates particles in the overlap zone are destroyed
to model the disappearance of ocean plate beneath the surface. When
two plates diverge, new oceanic plate is created by generating new
particles. The plate identity of a new particle depends on the
correlation of its velocity with the velocities of the plates on
either side of the gap.


A companion paper describes the consequences of using power law
rheology [10,11] instead of the simpler Newtonian creep law in a
model that also includes phase transitions. The result is to
dramatically increase the potential for episodes of catastrophic
avalanches [12,15,19,20] of cold material from the upper mantle into
the lower mantle. Numerical calculations that include such physics
require high spatial resolution and a robust scheme for treating
strong lateral variations in the effective shear viscosity. Although
this is currently feasible in two dimensions, computational costs
are still prohibitive in three dimensions. An approach used to work
around these limitations has been referred to as the Newtonian
analog method. In this approach the effects of a nonlinear
stress-dependent rheology are partially accounted for by simply
using a Newtonian deformation law and reducing the value of the
viscosity. Although this approach is far from satisfying, it is the
best that can be done from a numerical modeling standpoint at this
time. The results from such a strategy should therefore be
understood as merely suggestive of what the more accurate treatment
would provide. A further consideration is that spatial resolution
currently still restricts the realism even of 3-D global Newtonian
calculations. Some degree of scaling of parameters is usually
necessary for such 3-D calculations to be stable from a numerical
standpoint. One choice for reducing the steepness of the spatial
gradients and thereby achieving the required numerical stability is
to retain the desired value of Newtonian viscosity but to scale the
thermal conductivity and the radiogenic heat production rate to
values larger than those estimated for the real earth. This has the
effect of lowering the overall convective vigor of the system as
measured by the Rayleigh number Ra=agr2cvHd5/mok2 where d is the
depth of the mantle. The strategy then for mimicking the effects of
power law rheology in a Newtonian 3-D model with limited spatial
resolution is to select a reference viscosity mo that yields
appropriate velocities and to scale the thermal conductivity k and
radiogenic heat production H by the amount necessary to yield a
Rayleigh number low enough to be consistent with the available
spatial resolution. For the calculation described below, a reference
viscosity mo of 1 x 1013 Pa-s, a thermal conductivity of 2 x 1010 W
m-1K-1, and a radiogenic heat production rate of 0.02 W/m3 are used.


At a given instant in time the system of equations (1)-(5) can be
solved given only the temperature distribution and boundary
conditions. The only time dependent boundary condition is the plate
configuration. Therefore to initialize a calculation one needs only
an initial temperature distribution and plate configuration. The
calculation shown below assumes an extremely simple initial
temperature field related to an initial plate configuration that
represents a late Paleozoic/early Mesozoic reconstruction of the
supercontinent Pangea (Fig. 3a). This initial state is chosen
because plate motions since this point in earth history are tightly
constrained by observational data in today's ocean floors. The
realism of the calculation in some sense can thus be tested against
these observational constraints. Furthermore, geological evidence is
strong that Gondwanaland --the southern portion of Pangea that
included South America, Africa, India, Australia, and
Antarctica--was intact throughout the Paleozoic era and that North
America, Europe, and Africa also were not far apart during the
Paleozoic. Therefore, the actual pre-Flood continent distribution
may not have been much different from this Pangean configuration.
The plate boundaries in the Pacific hemisphere are chosen to
resemble the present ones. This implies the Pacific spreading ridges
have not migrated significantly since pre-Flood times. The
individual continental blocks represent the present continental
areas mapped to their estimated Pangean locations. Initially the
North American, Greenland, and Eurasian blocks are constrained to
have a common rotation vector. Similarly, the Gondwana blocks
initially rotate as a single unit. Later in the course of the
calculation, these composite blocks are allowed to break into
constituent parts with their own rigid motions.

The initial temperature distribution consists of the reference state
temperatures on which is superimposed a set of slablike
perturbations designed to represent incipient circum-Pangea
subduction. The perturbations have an amplitude of -400 K, a depth
extent of 400 km, and a width that corresponds to a single finite
element basis function (about 250 km). They lie along the Pangean
margin adjacent to South America, North America, and the Pacific and
Tethyan coasts of Asia and along an arc in the ocean from southeast
Asia, through what is now Indonesia and Australia as shown in Fig.
3a. Fig. 3b provides a cross sectional view through the earth in the
plane of the equator and reveals the modest depth extent of the
perturbations. Although they occupy but a tiny fraction of the total
volume of the mantle, these small perturbations are sufficient to
initiate a pattern of motions in the mantle that move the surface
plates by thousands of kilometers. The process, of course, is driven
by the gravitational potential energy existing in the cold upper
boundary at the beginning of the calculation.


Starting with these initial conditions, the numerical model is
advanced in time by solving the momentum, mass, and energy
conservation equations at every mesh point on each time step.
Tractions on the base of the surface plates produced by flow in the
mantle below causes the plates to move and their geometry to change.
Fig. 4 contains a sequence of snapshots at times of 10, 30, 50, and
70 days showing the locations of the continental blocks and the
velocities and temperatures at a depth of 100 km. A notable feature
in the velocity fields of Fig. 4 is the motion of the nonsubducting
continental blocks toward the adjacent zones of downwelling flow.
This motion is primarily a consequence of the drag exerted on a
nonsubducting block by the material below it as this material moves
toward the downwelling zone. Such a general pattern of flow is
evident in the cross-sectional slices of Fig 5. The translation of
the nonsubducting blocks in this manner leads to a backward, or
oceanward, migration of the zones of the downwelling. This oceanward
translation of the continental blocks as well as the subduction
zones therefore acts to pull the supercontinent apart. This behavior
is a basic fluid mechanical result and not the consequence of any
special initial conditions or unusual geometrical specifications
other than the asymmetrical downwelling at the edges of
nonsubducting portions of the surface. That the continental blocks
move apart without colliding and overrunning one another, on the
other hand, depends in a sensitive way on the initial distribution
of thermal perturbations, the shapes of the blocks, and timing of
their breakup. A moderate amount of trial and error was involved in
finding the special set of conditions that leads to the results
shown in Fig. 4 and 5.

An important output from the calculations is the height of the
surface relative to sea level. Fig. 6 shows global topography
relative to sea level at a time of 30 days. Several features are
noteworthy. One is the broad belt of depression and flooding of the
continental surface adjacent to subduction zones, as evident, for
example, along the western margins of North and South America. This
depression of the surface is mostly due to the stresses produced by
the cold slab of lithosphere sinking into the mantle below these
regions. Narrow trenches several kilometers in depth lie inside
these zones. A second feature is the elevation of the topography
above the oceanic spreading ridges. This effect is so strong that
some portions of the ridge are above sea level. Since the volume
occupied by the ridges displaces sea water, a result is to raise the
global sea level and to flood significant portions of the continent
interiors. A third effect is the elevation of continent areas
flanking zones of continental rifting. This is a consequence of the
intrusion of a significant volume of hot buoyant rock from deeper in
the mantle beneath these zones. This produces a belt of mountains
several kilometers high on either side of the rift zone between
North America and Africa, for example. It is worth emphasizing that
the topography dynamically changes with time and that Fig. 6 is but
a snapshot. It illustrates, however, that what is occurring in the
mantle below has a strong and complex effect on the height relative
to sea level of a given point at the earth's surface. Although this
calculation is crude and merely illustrative, it shows that this
mechanism produces the general type of vertical surface motions
required to create key aspects of the global stratigraphic record.
It produces broad scale continental flooding; it creates belts of
thick sediments at the edges of cratons; it uplifts portions of the
continents where broad scale erosion and scouring would be expected
to occur.


This calculation illustrates that with relatively modest initial
perturbations, gravitational potential energy stored in the earth's
upper thermal boundary layer drives an overturning of the mantle
that pulls the Pangean supercontinent apart, moves the continental
blocks by thousands of kilometers, elevates much of the newly formed
seafloor above sea level, floods essential all of the continental
surface, and produces dramatic downwarpings of the continent margins
that lie adjacent to zones of subduction.

The key to the short time scale is the phenomenon of power-law creep
that, for parameter values measured experimentally and for strain
rates observed in the calculation, yields more than eight orders of
magnitude reduction in effective viscosity relative to a condition
of zero strain rate. Indeed maximum strain rates implied by the
calculated velocities are on the order of 10-4 s-1 --precisely in
the range for which laboratory measurements have been made [10,11].
As discussed in more detail in the companion paper, the combination
of the effect of the endothermic phase transition at 660 km depth to
act as a barrier to vertical flow [12,15,19,20] with the tendency of
thermal runaway of regions of cold material from the upper thermal
boundary layer, makes a sudden catastrophic avalanche event a
genuine possibility. Thermal runaway behavior is a direct
consequence of the positive feedback associated with viscous heating
and temperature dependent rheology [1,9] and amplified by an extreme
sensitivity to strain rate. A notable outcome of the recent high
resolution mapping of the surface of Venus by the Magellan
spacecraft is the conclusion that there was a tectonic catastrophe
on Venus that completely resurfaced the planet in a brief span of
time [18]. This event in terms of radiometric time, accounting for
the uncertainties in the cratering rate estimates, coincides almost
precisely with the Flood event on earth. A mechanism internal to
Venus was almost certainly the cause of that catastrophe. It is
reasonable to suspect that simultaneous catastrophes on both the
earth and Venus likely were due to the same phenomenon of runaway
avalanche in their silicate mantles.

This mechanism of runaway subduction then appears to satisfy most of
the critical requirements imposed by the observational data to
successfully account for the Biblical Flood. It leads to a generally
correct pattern of large scale tectonic change; it produces flooding
of the continents; it causes broad uplifts and downwarpings of
craton interiors with intense downwarpings at portions of craton
margins to yield the types of sediment distributions observed. It
also transports huge volumes of marine sediments to craton edges as
ocean floor, in conveyor belt fashion, plunges into the mantle and
most of the sediment is scraped off and left behind. It plausibly
leads to intense global rain as hot magma erupted in zones of plate
divergence, in direct contact with ocean water, creates bubbles of
high pressure steam that emerge from the ocean, rise rapidly through
the atmosphere, radiate their heat to space, and precipitate their
water as rain. That no air-breathing life could survive such a
catastrophe and that most marine life also perished is readily
believable. Finally, numerical modeling appears to be the most
practical means for reconstructing a comprehensive picture of such
an event and for creating a conceptual framework into which the
geological observational data can be correctly integrated and
understood. This calculation, it is hoped, is a modest step in that


[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, A Three-Dimensional Finite Element Model
for Mantle Convection, Ph.D. thesis, 1983, UCLA.

[3] J. R. Baumgardner and P. O. Frederickson, Icosahedral
Discretization of the Two-Sphere, SIAM Journal of Numerical
Analysis, 22(1985), pp. 1107-1115.

[4] 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.

[5] 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.

[6] W. T. Brown, Jr., In the Beginning, 1989, Center for
Scientific Creation, Phoenix. [7] J. C. Dillow, The Waters Above,
1981, Moody Press, Chicago.

[8] A. M. Dziewonski and D. L. Anderson, Preliminary Reference
Earth Model, Physics of Earth and Planetary Interiors, 25(1981),
pp. 297-356.

[9] I. J. Gruntfest, Thermal Feedback in Liquid Flow; Plane Shear
at Constant Stress, Transactions of the Society of Rheology,
8(1963), pp. 195-207.

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

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

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

[14] D. W. Patton, The Biblical Flood and the Ice Epoch, 1966,
Pacific Meridian Publishing, Seattle.

[15] W. R. Peltier and L. P. Solheim, Mantle Phase Transitions and
Layered Chaotic Convection, Geophysical Research Letters,
19(1992), pp. 321-324.

[16] Proceedings of the Ocean Drilling Program

[17] 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.,

[18] R. G. Strom, G. G. Schaber, and D. D. Dawson, The Global
Resurfacing of Venus, Journal of Geophysical Research, 99(1994),
pp. 10899-10926. [19] 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.

[20] 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.

Promoting an Understanding of the Intelligent Design of the Universe