Genesis Flood

 
COMPUTER MODELING OF
THE LARGE-SCALE TECTONICS
ASSOCIATED WITH THE GENESIS FLOOD
    
 
Presented: Third International Conference on Creationism
July 18-23, 1994
Copyright 1994 by Creation Science Fellowship, Inc.
Pittsburgh, PA  USA - All Rights Reserved


KEYWORDS

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


ABSTRACT

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.


INTRODUCTION

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.


MATHEMATICAL FORMULATION

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= -v (p - pr) + (r - rr) g + v. t
(1)
 
0=v. (r u)
(2)
 
dT/dt=-v. (T u) - (g - 1) T v. u + [v. (k v T) + t : v u + H]/rrcv

(3)
where
t =m [v u + (v u)T - 2 I (v. u)/3]
(4)
and
r=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 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.


THE REFERENCE STATE

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)

 (6)
where
E=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 follows:

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.


PHASE CHANGES

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.


NUMERICAL APPROACH

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.


TREATMENT OF THE RUNAWAY INSTABILITY

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.


INITIAL CONDITIONS

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.


RESULTS

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.


CONCLUSIONS

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.


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, 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., 1992.

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