LOS ALAMOS, NM 87544 USA



(Published in Proceedings of the Fifth International Conference on Creationism, R. L. Ivey, Jr., Editor,
Creation Science Fellowship, Pittsburgh, PA, 113-126, 2003; also posted at
http://globalflood.org
.)


KEYWORDS: Genesis Flood, catastrophic plate tectonics, runaway subduction, mantle rheology,
stress-weakening

ABSTRACT

The wealth of new data, mostly from the ocean bottom, that precipitated the acceptance of plate
tectonics during the 1960s simultaneously also opened the door for the first time in more than 200 years
to a technically credible defense of the Genesis Flood. From the mid-1700s through the days of Hutton,
Lyell, and Darwin to the 1960s, it overwhelmed the human mind to imagine a mechanism that could
possibly deliver, in a single brief event, the magnitude and complexity of geological change evident in the
continental rock record above the point where fossils first appear. However, with the new awareness
that the Earths interior could participate in the process and that the stiff layer of rock some 50 miles thick
beneath the oceans could be recycled into the Earth, the stage was set for a breakthrough in regard to
the mechanism for the Flood cataclysm. The crucial final piece of the puzzle has come from laboratory
experiments that have carefully measured the way in which silicate minerals deform under conditions of
high temperature and high stress. These experiments reveal silicate material can weaken dramatically,
by factors of a billion or more, at mantle temperatures and for stress conditions that can exist in the
mantles of planets the size of the Earth. The scenario in which all the Earths ocean lithosphere is
rapidly recycled into the mantle via a runaway process, enabled by this stress-weakening behavior, is
now known as catastrophic plate tectonics [4]. Evidence in the geological record is compelling that such
a cataclysmic episode indeed has occurred in the Earths recent past. A reasonable inference is that this
event corresponds to the Flood described in the Bible and other ancient sources. I report new
computational results from 2D and 3D simulations of this catastrophic plate tectonics process. In
particular, I describe how fundamental advances in computational techniques now make it possible to
advance the numerical solution successfully through the most extreme phase of the runaway regime.

INTRODUCTION

At least as far back as the early 1960s it has been known that the phenomenon of thermal runaway can
potentially occur in materials whose effective viscosity is described by an Arrhenius-like [20] relationship.
The viscosity of such materials varies as e
(E*/RT)
, where T is absolute temperature, E* is the activation
energy, and R is the gas constant. A large variety of materials, including silicate minerals, have
viscosities that vary with temperature in this manner. In 1963 I. Gruntfest showed for a layer subject to
constant applied shear stress and a viscosity with Arrhenius temperature dependence, both the
deformation rate and the temperature within the layer can increase without limit, that is, run away [15].
The criterion for runaway to occur is that the time constant associated with viscous heating be much less
than the characteristic thermal diffusion time of the layer. Several investigators in the late 1960s and
early 1970s explored the possibility of thermal runaway of lithospheric slabs in the mantle. Anderson
and Perkins [2], for example, suggested that the widespread Cenozoic volcanism in the southwestern
U.S. might be a consequence of thermal runaway of chunks of lithosphere in the low-viscosity upper
mantle. They conjectured that surges of melt associated with such runaway events might account for
episodes of volcanism observed at the surface. Lithospheric slabs, because they display an average
temperature some 1000 K or more lower than that of the upper mantle but have a similar bulk chemical composition, are several percent denser than the surrounding upper mantle rock and therefore have a
natural ability to sink. The gravitational body forces acting on a slab lead to high stresses, especially
within the mechanical boundary layer surrounding the slab. As a slab sinks, most of its gravitational
potential energy is released in the form of heat in these regions of high deformation. If conditions are
right, the weakening arising from heating can lead to an increased sinking rate, an increased heating
rate, and greater weakening. This positive feedback associated with thermal weakening can result in
runaway provided the criterion mentioned above is met [5].

Experimental studies of the deformational behavior of silicate minerals over the last several decades
have revealed the strength of such materials not only depends strongly on the temperature but also on
the deformation rate. At shear stresses on the order of 10
-3
times the low-temperature elastic shear
modulus and temperatures on the order of 80% of the melting temperature, silicate minerals deform by a
mechanism known as dislocation creep in which slip occurs along preferred planes in the crystalline
lattice [19]. In this type of solid deformation, material strength depends on the deformation rate in a
strongly nonlinear manner, proportional to the deformation rate to approximately the minus two-thirds
power. At somewhat higher levels of shear stress, these materials display another type of deformational
behavior known as plastic yield, where their strength decreases in an even more nonlinear way, in this
case, inversely with the deformation rate (i.e., proportional to the deformation rate to the minus one
power). When these deformation-rate-weakening mechanisms are combined with the temperature
weakening discussed above, the potential for slab runaway from gravitational body forces is enhanced
dramatically. A point many people fail to grasp is that these weakening mechanisms can reduce the
silicate strength by ten or more orders of magnitude without the material ever reaching its melting
temperature [19].

BREAKTHROUGH IN NUMERICAL MODELING OF THE RUNAWAY MECHANISM

Numerical methods now exist for modeling and investigating this runaway mechanism. Considerable
challenge is involved, however, because of the extreme gradients in material strength that arise [6, 8].
W.-S. Yang, a graduate student with whom I worked closely, focused much of his Ph.D. thesis research
effort at the University of Illinois on finding a robust approach for dealing with such strong gradients in
the framework of the finite element method and an iterative multigrid solver. He showed what is known
as a matrix dependent transfer multigrid approach allows one to treat such problems with a high degree
of success. Although his thesis dealt with applying this method to 3D spherical shell geometry, he
subsequently developed a simplified 2D Cartesian version capable of much higher spatial resolution.
Details of this method together with some sample calculations are provided in a recent paper [27].

This new formulation of the multigrid solver represents a breakthrough in treating large local variations in
rock strength and allows the mantle runaway process to be modeled to completion for the very first time.
Results I have reported in previous ICC papers only tracked the runaway to its earliest stages. Beyond
that point available numerical methods failed. Although the underlying equations themselves indicated
runaway most certainly would occur, computer methods were not available that could handle fully
developed runaway conditions. Moreover, the new solver technique now allows a regime of rock
deformation known as plastic yield that involves an even greater degree of instability. This important
plastic flow regime, because of the increased level of instability it introduces, had not been included in
previous efforts to model the runaway process.

Figure 1 is a plot of the primary deformation regimes as determined by many careful laboratory
experiments for the common mantle mineral olivine. The heavy lines separate the three main regimes:
diffusion creep, dislocation or power-law creep, and plastic yield. Finer lines of constant shear strain
rates are plotted as a function of temperature and shear stress. (For readers unfamiliar with the
terminology, strain has to do with the amount of deformation per unit length and so is dimensionless.
Strain rate is the change in strain per unit time and so has units of inverse time. Stress has units of force
per unit area, the same as pressure.) Note that the rates of strain, or deformation, displayed in this plot
for these solid olivine crystals vary over fourteen orders of magnitude! This range of deformation rate
easily brackets the rates observed in the runaway calculations. (A tectonic plate moving 10 m/s, or 22.4
mph, relative to some substrate below, with a 10 km thick weak zone in between, implies an average
shear strain rate of 10
-3
within the weak zone, for example.) In regard to the three regimes, diffusion
creep involves migration of point defects (extra or missing atoms) through the crystalline lattice in
response to applied stress, while dislocation creep involves planes of atoms moving relative to each
other in a more or less coherent way. In the plastic yield regime, such large numbers of dislocations
emerge that huge increases in deformation rate occur with very little increase in shear stress.

2
Figure 1. De