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
Copyright 1994 by Creation Science Fellowship, Inc.
Pittsburgh, PA  USA - All Rights Reserved

Genesis Flood, geological catastrophism, runaway subduction, mantle dynamics, plate tectonics

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:

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 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) heating.

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.
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 direction.