KEYWORDS

Runaway subduction, Genesis flood, power-law creep, thermal runaway, catastrophic plate tectonics


ABSTRACT

Experimental investigation of the solid state deformation properties of silicates at high temperatures has revealed that the deformation rate depends on the stress to a power of about 3 to 5 as well as strongly on the temperature. This highly nonlinear behavior leads to the potential of thermal runaway of the mantle's cold upper boundary layer as it peels away from the surface and sinks through the hot mantle. The additional fact that the mineral phase changes that occur at 660 km depth act as a barrier to convective flow and lead to a tendency for large episodic avalanche events compounds the potential for catastrophic dynamics. Two-dimensional finite element calculations are presented that attempt to model these strongly nonlinear phenomena. It is proposed that such a runaway episode was responsible for the Flood described in Genesis and resulted in massive global tectonic change at the earth's surface.


INTRODUCTION

The only event in Scripture since creation capable of the mass destruction of living organisms evident in the fossil record is the Genesis Flood. A critical issue in any model for earth history that accepts the Bible as accurate and true is what was the mechanism for this catastrophe that so transformed the face of the earth in such a brief span of time. The correct answer is crucial to understanding the Flood itself and for interpreting the geological record in a coherent and valid manner. It is therefore a key element in any comprehensive model of origins from a creationist perspective. Ideas proposed as candidate mechanisms over the past century include collapse of a water vapor canopy [5], near collision of a large comet with the earth [12], rapid earth expansion [11], and violent rupture of the crust by pressurized subterranean water [4]. There are serious difficulties with each of these ideas.

Another possibility is that of runaway subduction of the pre-Flood ocean lithosphere [2,3]. A compelling logical argument in favor of this mechanism is the fact that there is presently no ocean floor on the earth that predates the deposition of the fossiliferous strata. In other words all the basalt that comprises the upper five kilometers or so of today's igneous ocean crust has cooled from the molten state since sometime after the Flood cataclysm began. The age of today's seafloor relative to the fossil record is based on two decades of deep sea drilling and cataloging of fossils in the sediments overlying the basalt basement by the Deep Sea Drilling Program as well as radiometric dating of the basalts themselves [14]. Presumably, there were oceans and ocean floor before the Flood. If this pre-Flood seafloor did not subduct into the mantle, what was its fate? Where are these rocks today? On the other hand, if the pre-Flood seafloor did subduct, it must have done so very rapidly --within the year of the Flood. In regard to the fate of the pre-Flood seafloor, there is strong observational support in global seismic tomography models for cold, dense material near the base of the lower mantle in a belt surrounding the present Pacific ocean [16]. Such a spatial pattern is consistent with subduction of large areas of seafloor at the edges of a continent configuration commonly known as Pangea.

There are good physical reasons for believing that subduction can occur in a catastrophic fashion because of the potential for thermal runaway in silicate rock. This mechanism was first proposed by Gruntfest [6] in 1963 and was considered by several in the geophysics community in the early 1970's [1]. Previous ICC papers [2,3] have discussed the process by which a large cold, relatively more dense, volume of rock in the mantle generates deformational heating in an envelope surrounding it, which in turn reduces the viscosity in the envelope because of the sensitivity of the viscosity to temperature. This decrease in viscosity in turn allows the deformation rate in the envelope to increase, which leads to more intense deformational heating, and finally, because of the positive feedback, results in a sinking rate orders of magnitude higher than would occur otherwise. It was pointed out that thermal diffusion, or conduction of heat out of the zone of high deformation, competes with this tendency toward thermal runaway. It was argued there is a threshold beyond which the deformational heating is strong enough to overwhelm the thermal diffusion, and some effort was made to characterize this threshold.

The important new aspect addressed in this paper is the dependence of the viscosity on the deformation rate itself. Although this deformation rate dependence of viscosity has been observed experimentally in the laboratory for several decades, the difficulty of treating it in numerical models has deterred most investigators from exploring many of its implications. Results reported in the previous ICC papers did not include this highly nonlinear phenomenon. Significant improvements in the numerical techniques that permit large variations in viscosity over small distances in the computational domain, however, now make such calculations practical. The result of including this behavior in the analysis of the thermal runaway mechanism is to discover a much stronger tendency for instability in the earth's mantle. Moreover, deformation rates orders of magnitude higher than before throughout large volumes of the mantle now can be credibly accounted for in terms of this more realistic deformation law. This piece of physics therefore represents a major advance in understanding how a global tectonic catastrophe could transform the face of the earth on a time scale of a few weeks in the manner that Genesis describes Noah's Flood.

Recent papers by several different investigators [10,13,18,19] have also shown that the mineral phase changes which occur as the pressure in the mantle increases with depth also leads to episodic dynamics. The spinel to perovskite plus magnesiowustite transition at about 660 km depth is endothermic and acts as a barrier to flow at this interface between the upper and lower mantle. It therefore tends to trap cold material from the mantle's upper boundary layer as it peels away from the surface and sinks. Numerical studies show that, with this phase transition present, flow in the mantle becomes very episodic in character and punctuated with brief avalanche events that dump the cold material that has accumulated in the upper mantle into the lower mantle. The episodic behavior occurs without the inclusion of the physics that leads to thermal runaway. This paper argues that when temperature and strain rate dependence of the rheology is included, the time scale for these catastrophic episodes is further reduced by orders of magnitude. In this light, the Flood of the Bible with its accompanying tectonic expressions is a phenomenon that is seems to be leaping out of the recent numerical simulations.
 
 
MATHEMATICAL FORMULATION
 
In this numerical model the silicate mantle is treated as an infinite Prandtl number, anelastic fluid within a domain with isothermal, undeformable, traction-free boundaries. Under these approximations the following equations describe the local fluid behavior:

	0=- (p - pr) + (r - rr) g + t	(1)
	0= (r u)	(2)
	dT/dt=- (T u) - (g - 1) T u + [ (k T) + t : u + H]/rrcv	(3)
where	t =m [ u + ( u)T - 2 I ( 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, 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 pressure, density, and temperature of the reference state. I is the identity tensor. The superscript T in (4) denotes the tensor transpose. Equation (1) expresses the conservation of momentum in the infinite Prandtl number limit. In this limit, the deformational term is so large that the inertial terms normally needed to describe less viscous fluids may be completely ignored. The resulting equation (1) then represents the balance among forces arising from pressure gradients, buoyancy, and deformation. Equation (2) expresses the conservation of mass under the anelastic approximation. The anelastic approximation ignores the partial derivative of density with respect to time in the dynamics and thereby eliminates fast local density oscillations. It allows the computational time step to be dictated by the much slower deformational dynamics. Equation (3) expresses the conservation of energy in terms of absolute temperature. It includes effects of transport of heat by the flowing material, compressional heating and expansion cooling, thermal conduction, shear or deformational heating, and local volume (e.g., radiogenic) heating.

The expression for the deviatoric stress given by equation (4) assumes the dynamic shear viscosity m depends on temperature, pressure, and strain rate. The stress therefore is nonlinear with respect to velocity, and the rheological description is non-Newtonian. This formulation is appropriate for the deformation regime in solids known as power-law creep to be discussed below. Equation (5) represents density variations as linearly proportional to pressure and temperature variations relative to a simple reference state of uniform density, pressure and temperature. Parameter values used are r_r =3400 kg m^-3, p_r=0, T_r=1600 K, g=10 m/s, g=1, k=4 W m^-1K^-1, H=1.7 x 10^-8 W m^-3, c_v =1000 J kg^-1K^-1, and K=1 x 10¹² Pa.
 
 
POWER-LAW CREEP
 
Laboratory experiments to characterize the high temperature solid state deformation properties of silicates have been carried out by many investigators over the last three decades [8,9]. These experiments have established that, for temperatures above about sixty percent of the melting temperature and strain rates down to the smallest achievable in the laboratory, silicate materials such as olivine deform according to a relationship of the form [8]

=A sn exp[ -(E* + pV*)/RT]

(6)

where e is the strain rate, A a material constant, s the differential stress, n a dimensionless constant on the order of 3 to 5, E* an activation energy, p is pressure, V* an activation volume, R the universal gas constant, and T absolute temperature. This relationship implies that at constant temperature and pressure the deformation rate increases dramatically more rapidly than the stress. Because the strain rate increases as the stress to some power greater than one, this type of deformation is known as power-law creep. This relationship may also be expressed in terms of an effective viscosity m=0.5s/e that depends on the strain rate e as [9, 17, p. 291]

m=B -q exp[(E* + pV*)/nRT]

(7)

where B=0.5A^-1/n and q=1 - 1/n. A value for n of 3.5, appropriate for the mineral olivine [8,9], yields a q of 0.714. This means that the effective viscosity m decreases strongly as the strain rate increases. A tenfold increase in the strain rate, for example, yields an effective viscosity, at fixed temperature and pressure, a factor of 5.2 smaller! For a 10¹⁰ increase in strain rate, the effective viscosity decreases by more than a factor of 10⁷. The effect is even more pronounced for larger values of n.

Fig. 1 is a deformation mechanism map for olivine that shows the region in stress-temperature space where power-law creep is observed. Note that there exists a boundary between the power-law creep regime and that of diffusional creep. Because the strain rates for diffusional creep are so small--too small in fact to be realized in laboratory experiments--this boundary is poorly constrained. Kirby [8, p. 1461] states that the boundary may in actuality be substantially to the left of where he has drawn it. In any case at a given temperature there is a threshold value for the strain rate at which point one crosses from the diffusional regime--where the strain rate depends linearly on the stress--into the power-law regime. From Fig. 1 this threshold is on the order of 10^-17 to 10^-14 s^-1 for temperatures about 60% of the melting temperature and stresses of about 1 MPa.

Power-law creep is included in the numerical model simply by using the effective viscosity given by (7) in (4), where the scalar strain rate e is obtained by taking the square root of the second invariant of the rate of strain tensor d=(u + u^T)/2. To remove the singularity in (7) for zero strain rate and to model explicitly the transition between diffusion creep and power-law creep, a minimum or threshold strain rate _o is incorporated into the formulation. For regions in the domain where the strain rate exceeds _o, equation (7) applies. Otherwise the viscosity is strain rate independent. The parameter B is specified in terms of a reference viscosity m_o at reference temperature T_r and zero strain rate as B=m_o/{_o^-q exp[(E* + pV*)/nRT_r]}. To model the viscosity contrast between the upper mantle and lower mantle, the reference viscosity is allowed to vary with depth and increase in a linear fashion by a factor of 50 between 400 and 700 km. For purposes of numerical stability the threshold strain rate _o is assumed to vary as 1/m_o.
 
 
PHASE CHANGES
 
The jumps in seismic quantities observed at depths of about 410 km and 660 km in the earth closely match phase transitions observed in laboratory experiments at similar temperatures and pressures for olivine to spinel and from spinel to perovskite silicate structures, respectively. These phase transitions that occur as the pressure increases and the crystal structures assume more compact configurations almost certainly play a critical role in the mantle's dynamical behavior. In a calculation in which silicate material is transported through these depths and undergoes these phase changes, two effects need to be taken into account. One is the latent heat released or absorbed and the other is the deflection of the phase boundary upward or downward. The latent heat may be accounted for by adding or removing heat through the volume heating term in equation (3) proportional to the vertical flux of material through the transition depth. The latent heat per unit mass is obtained from the Clapeyron equation which expresses that in a phase transition DH=(dp/dT) T DV, where DH is the enthalpy change, or latent heat, and DV is the change in specific volume. The Clapeyron slope (dp/dT) is a quantity that can be determined experimentally for a given transition. The deflection in the location of a phase boundary occurs because the pressure, and therefore the depth, at which the phase change occurs depends on the temperature. The effect of such a deflection enters as a contribution to the buoyancy term in equation (1). A downward deflection represents positive buoyancy because the lighter phase now occupies volume normally occupied by the denser phase. The Clapeyron slope is also a constant of proportionality in the boundary deflection Dh=-(dp/dT) DT/rg that arises from a deviation DT from the reference temperature. The values for the Clapeyron slope used here are 1 x 10⁶ Pa/K for the 410 km transition and -2 x 10⁶ Pa/K for the 660 km transition. Note that the exothermic 410 km transition leads to a positive or upward deflection for a cold slab and hence increased negative buoyancy, while the endothermic 660 km transition leads to a downward deflection and reduced negative buoyancy. The 660 km transition therefore acts to inhibit buoyancy driven flow while the 410 km transition acts to enhance it.
 
 
NUMERICAL APPROACH

The set of equations (1)-(5) is solved in a discrete manner on a uniform rectangular mesh with velocities located at the mesh nodes and temperatures, pressures, and densities at cell centers. Piecewise linear finite elements are used to represent the velocity field, while the cell centered variables are treated as piecewise constant over the cells. The calculational procedure on each time step is first to apply a two-level conjugate gradient algorithm [15] to compute the velocity and pressure fields simultaneously from Eq. (1) and (2). This task involves solving 3n simultaneous equations for 2n velocity unknowns and n pressure unknowns, where n is the total number of nodes in the mesh. Key to the procedure is an iterative multigrid solver that employs an approximate inverse with a 25-point stencil. This large stencil for the approximate inverse enables the method to handle large variations in viscosity in a stable fashion. The outstanding rate of convergence in the multigrid solver is responsible for the method's overall high efficiency. The piecewise linear finite element basis functions provide second-order spatial accuracy. The temperature field is updated according to Eq. (3) with a forward-in-time finite difference van Leer limited advection method.
 
 
RESULTS

Two calculations will now be described that illustrate the effects of power law creep on the stability of a sinking slab. The two calculations are identical except for the value of the strain rate threshold above which power law creep occurs. In the first case, the threshold _o in the upper 400 km is 3 x 10^-13 s^-1 which is sufficiently large that no power law creep occurs anywhere in the domain. In the second case, the threshold is 6.5 x 10^-14 s^-1, about a factor of five smaller. In this case runaway eventually takes place. These calculations are performed in a rectangular domain 2890 km high and 1280 km wide on a mesh with 64 x 64 cells of uniform size. The viscosity m_o at zero strain rate and 1600 K increases in a simple linear fashion by a factor 50 between 400 km and 700 km depth to represent the stiffer rheology of the earth's lower mantle compared with the upper mantle. The phase changes at 410 km and 660 km depth are both included. The endothermic phase transition at 660 km as well as the higher intrinsic viscosity below this depth both act to inhibit flow from above. The calculations are initialized with a uniform temperature of 1600 K except for a slablike anomaly 100 km wide extending from the top to a depth of 400 km with a central temperature of 1000 K and a thermal boundary layer at the top such that the temperature in the topmost layer of cells is initially 1000 K. The upper boundary temperature is fixed at 700 K and the bottom at 1600 K.

Fig. 2 shows four snapshots in time spaced roughly 6 x 10⁶ years apart of the calculation with the larger strain rate threshold. Plots of temperature and effective viscosity are displayed with velocities superimposed. Note that the initial maximum velocity drops by a factor of two as the slab encounters increasing resistance from the higher viscosity and 660 km phase change. The colder material tends to accumulate and thicken in width in the depth range between 400 and 700 km. When sufficient thickening of the zone of cold material has occurred, it begins to penetrate slowly into the region below 700 km.

Fig. 3 shows the effects of a strain rate threshold _o sufficiently low that power law creep is occurring in a significant portion of the problem domain. The first three snapshots in time for this case resemble those for the previous case. The main difference are regions of reduced effective viscosity in the region below 700 km evident in the first and third snapshots due to the power-law rheology. A major change is observed, however, in the fourth snapshot with an increase in peak velocity and a notable reduction in effective viscosity below the head of the developing cold plume. In the fifth snapshot the peak velocity has increased by another 80% and there is more than a factor of ten reduction in the effective viscosity ahead of the plume. Also displayed in this snapshot is the viscous heating rate that shows intense heating surrounding the plume. In the sixth snapshot, the head of the cold plume is preceded by a belt of high temperature, the velocity has almost doubled again, the effective viscosities near the plume have dropped even further, and the heating rate adjacent to the plume has more than doubled. Shortly after this point in the calculation, runaway occurs and the computation crashes.
 
 
DISCUSSION

What do these calculations have to say about the mantle and the Flood? First of all, power-law rheology dramatically enhances the potential for thermal runaway. Numerical calculations are not really necessary to reach this conclusion. Equation (7) indicates an increase in the deformation rate leads to a reduction in the effective viscosity and reinforces the reduction in viscosity an increase in temperature provides. These effects work together in a potent way. An exciting further consequence of the power-law rheology is that high velocities and strain rates can now occur throughout the mantle. A hint of this can be inferred from the last two snapshots in Fig. 3. Large and increasing velocities are not just associated with the sinking plume itself but are observed throughout the domain. The remaining horizontal sections of the initial cold upper boundary layer, for example, are also moving at much higher speeds.

In interpreting these numerical experiments it is important to realize that one is attempting to explore numerically a physically unstable process. Customary numerical difficulties associated with strong gradients in the computed quantities are compounded when such a physical instability occurs. The strategy is to explore the region of parameter space nearby but not too close to where the instability actually lives. The calculation of Fig. 3 therefore does not reveal the true strength of the instability relative to the situation of a moderately lower value for the threshold strain rate. It is also useful to point out how various quantities scale relative to one another. The velocities are inversely proportional to the reference viscosity. A tenfold reduction in the reference viscosity gives ten times higher velocities. Similarly, the threshold strain rate for runaway behavior is inversely proportional to the reference viscosity since strain rate is proportional to velocity. So reducing the reference viscosity by a factor of ten yields a threshold strain rate for runaway ten times larger. This neglects the diminished influence of thermal diffusion at the higher velocities.

How do the parameters used in these calculations compare with those estimated for the earth? The values used for g, g, k, H, r_r, c_v, T_r, and a in eq. (1)-(5) are all reasonable to within +/-30% for the simplified reference state that is employed. The values used for the Clapeyron slopes for the phase transitions are two to three times too small and so the effects of the phase changes are underrepresented. The most important parameters are the reference viscosity and the threshold strain rate for power-law creep. The reference viscosity leads to velocities prior to runaway that are in accord with current observed plate velocities of a few centimeters per year. The threshold strain rates used are within the power-law creep region for olivine as given by Kirby (Fig. 1). A large uncertainty is the extrapolation of the creep behavior of olivine to the minerals of the lower mantle for which there is essentially no experimental data. The issue is not whether power-law creep occurs in these minerals but what the stress range is in which it occurs. It is likely the threshold strain rate is not many orders of magnitude different from olivine. These calculations therefore seem relevant to the earth as we observe it today.

One difficulty in making a connection between these calculations and the Flood is their time scale. Some 2 x 10⁷ years is needed before the instability occurs in the second calculation. Most of this time is involved with the accumulation of a large blob of cold, dense material at the barrier created by the phase transition at 660 km depth. This time span disappears when the initial condition consists of a large belt of cold material already trapped above this phase transition in the pre-Flood mantle. A relatively small amount of additional negative buoyancy in such a belt can then trigger runaway. One means for providing a quick pulse of negative buoyancy is by the sudden conversion to spinel of olivine in a metastable state that resides at depths below the usual transition depth of about 410 km. Such metastability can arise because the changes in volume and structure associated with a phase transition do not necessarily occur spontaneously as transition conditions are reached, especially if the material is cold. Some means of nucleation of seed crystals of the new phase is generally required. If such nucleation does not happen, then substantial amounts of the less dense phase can survive to depths much greater than what the assumption of a spontaneous transition would imply. Indeed, there is observational evidence for significant amounts of metastable olivine in the slab currently beneath Japan [7]. A shock wave passing through such a volume of metastable material can initiate the nucleation and cause a sudden conversion to the denser phase. Present day deep focus earthquakes likely represent manifestations of such a process on a small scale. In the context of the Flood, it is conceivable that an extraterrestrial impact of modest size could have triggered a sudden conversion of metastable material to the denser phase and the resulting earthquakes then propagated in a self-sustaining manner to convert the metastable material throughout much of the upper mantle to the denser spinel phase, which in turn initiated the runaway avalanche of upper mantle rock into the lower mantle. It is also conceivable that a single large earthquake generated by causes internal to the earth could have been the event that caused a sudden conversion of the metastable material and then the runaway avalanche.
 
 
CONCLUSIONS

Rapid sinking through the mantle of portions of the mantle's cold upper boundary facilitated by the process of thermal runaway appears to be a genuine possibility for the earth. A highly nonlinear deformation law for silicate minerals at conditions of high temperature known as power-law creep, documented by decades of experimental effort, in which the effective viscosity decreases strongly with the deformation rate, makes thermal runaway almost a certainty for a significant suite of conditions. This deformation law also makes possible strain rates consistent with large scale tectonic change within the Biblical time frame for the Flood. Mineralogical phase changes combined with the viscosity contrast between upper and lower mantle conspire to provide the setting in which a sudden triggering of a runaway avalanche of material trapped in the upper mantle into the lower mantle can occur. Calculations by other investigators that include the endothermic phase transition, but not temperature or strain rate dependent viscosity, also display the tendency for episodic avalanche events [10,13,18,19]. Such an episode of catastrophic runaway is here presented as the mechanism responsible for Noah's Flood.