The onset of convection in fluids with strongly temperature-dependent viscosity cooled from above with implications for planetary lithospheres

Earth and Planetary Science Letters 224 (2004) 371 – 386 www.elsevier.com/locate/epsl

The onset of convection in fluids with strongly temperature-dependent viscosity cooled from above with implications for planetary lithospheres Sarah E. Zaranek *, E.M. Parmentier Department of Geological Sciences, Brown University, Providence, RI 02912, USA Received 17 November 2003; received in revised form 6 May 2004; accepted 10 May 2004

Abstract The convective instability of a viscous fluid with strongly temperature-dependent viscosity cooled from above provides a basis for better understanding the thermal evolution of the Earth and other planets. The role of temperature-dependent viscosity in the initiation of thermal convection is important since thermally activated creep is strongly temperature dependent. This study introduces a criterion for convective instability based on the ratio of the Rayleigh – Taylor (R-T) growth rate of the thermal density stratification due to conductive cooling and the rate of conductive smoothing of convectively generated temperature variations. This definition of a critical boundary layer Rayleigh number that must be reached for the onset of convection is consistent with numerical experiments reported here as well as earlier theoretical and experimental studies. The length and temperature scales appearing in this Rayleigh number provide valuable physical insight into the behavior of the fluid at the onset of convection. A consistently defined viscosity ratio across the convecting thermal boundary layer leads to a characteristic temperature difference (DTc) across this layer. The viscosity ratio inferred from the onset times in numerical experiments reported here (c 10) is similar to that seen in laboratory experiments [Geophys. Res. 99 (1994) 19853] and is consistent with the linearized R-T analysis. Calculation of the onset time of convection and the conductive lid thickness at onset, based on laboratory-derived rheological parameters, provides a better basis to assess the convective instability of the oceanic upper mantle and the possible formation of a thick cratonic lithosphere by the cooling of continental mantle. D 2004 Elsevier B.V. All rights reserved. Keywords: convection; temperature-dependent viscosity; lithosphere

1. Introduction Characterizing the convective cooling of a viscous fluid with strongly temperature-dependent viscosity provides a basis for better understanding the thermal

* Corresponding author. E-mail address: sarah _ [email protected] (S.E. Zaranek). 0012-821X/$ - see front matter D 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.epsl.2004.05.013

evolution of the Earth and other planets. In a fluid cooling from above, the conductively cooled boundary layer thickens with time, eventually becoming convectively unstable. In this study, both 2D numerical experiments and a linearized stability analysis are used to study the effects of temperature-dependent viscosity on the onset of convective motions in an infinite Prandtl number, Boussinesq fluid [2]. The viscosity describing thermally activated solid-state

372

S.E. Zaranek, E.M. Parmentier / Earth and Planetary Science Letters 224 (2004) 371–386

creep in planetary mantles depends strongly on temperature and, over the range of temperatures expected to be present, results in large viscosity variations. In numerical experiments, the magnitude of the viscosity variations that can be well resolved is limited by the change of viscosity between adjacent points in the discretization. The goal of this study is to derive scaling laws that allow extrapolation from numerically resolvable viscosity variations to the much larger variations expected in planetary mantles. The onset of thermal convection is a classical problem dating from the work of Rayleigh and Benard. The classical Rayleigh – Benard problem treating the stability of a fluid layer heated from below has been examined in many previous studies, both theoretical and experimental [3]. Other studies, although fewer in number, have considered a fluid cooled from above, the primary mode of cooling for planetary mantles. Convective instability in a fluid cooling from above differs from the classical case of heating from below in several ways. First, the temperature distribution in the basic state prior to the onset of convection takes the form of a conductively cooled boundary layer above fluid at its initial uniform temperature. A second fundamental difference is that this temperature distribution is changing with time as instability develops. Early theoretical and experimental studies [4– 7] pointed out the importance of this time-dependent basic state on the onset of convection. The earlier studies cited above considered a constant viscosity fluid. For a fluid with viscosity that increases with decreasing temperature, as expected for thermally activated creep, a relatively stagnant, high viscosity lid, develops near the cool top boundary. Heat is transferred through this lid only by conduction. Convective instability in the presence of such a conductive lid was treated theoretically by Yuen and Fleitout [8] with a time-dependent basic state and by Jaupart and Parsons [9] with a quasi-static or ‘‘frozen time’’ analysis. These studies treated a temperatureand pressure-dependent and a purely depth-dependent viscosity, respectively. Jaupart and Parsons indicate that the fluid remains essentially stagnant in a region near the top boundary where the viscosity exceeds that in the deeper, hotter isothermal interior by a factor of about 100. For the temperature distribution due to conductive cooling that exists just prior to the onset of

convection, this viscosity ratio can be used to identify the thickness and temperature difference across the convecting boundary layer at the onset of convection. Laboratory experiments of Davaille and Jaupart [1] studied the effect of the presence of a stagnant lid in fluids with strongly temperature-dependent viscosity on the onset of convection, as well as the convective heat flux subsequent to the onset. These experiments used Golden Syrup or silicon oil as the fluid in a tank with a cooled top boundary and insulated sides and bottom. In contrast to the results of Jaupart and Parsons, the convecting part of the thermal boundary in these experiments was found to be restricted to a viscosity increase of about a factor of 10. Results presented here verify the stagnant lid criterion of Davaille and Jaupart. We develop a new physically based criterion for the onset of convection and compare this prediction with new numerical experiments as well as recent results of Korenaga and Jordan [10], Huang et al. [11], and Zaranek and Parmentier [12]. The results are then applied to better understand the age at which thermal convection initiates beneath the oceanic lithosphere and the possible role of conductive ‘‘top-down’’ cooling in the origin of thick continental lithosphere beneath cratons. In both cases, we identify the rheological parameters that control the age and thickness of the thermal lithosphere or conductive lid at the onset of convective instability.

2. Temperature dependence of viscosity The rate of thermally activated creep in planetary mantles obeys an Arrhenius-type law with a temperature dependence described by an activation energy, Q. The resulting viscosity can be expressed as lðT Þ ¼ li expðeA ððTi =T Þ  1ÞÞ

ð1Þ

where eA ¼

Q ; RTi

ð2Þ

Ti is the initial fluid temperature in Kelvin, the reference viscosity li = l(Ti), and R is the universal gas constant. A typical value Q = 500 kJ/mol for

S.E. Zaranek, E.M. Parmentier / Earth and Planetary Science Letters 224 (2004) 371–386

373

dislocation creep gives eA c 38 so that a small decrease in temperature results in a large increase in viscosity. The viscosity of fluids frequently used in laboratory experiments does not have an exactly Arrheniuslike dependence on temperature. Golden Syrup has an Arrhenius-like variation of viscosity with temperature. But, for example, silicon oil has an exponential-like viscosity – temperature dependence. In numerical experiments, an exponential dependence has often been adopted for simplicity. The exponential viscosity – temperature relationship lðT Þ ¼ li expðeexp ðTi  T Þ=DTi Þ;

ð3Þ

where DTi is the difference between the initial interior temperature and the top boundary temperature, appears to be a good approximation to the Arrhenius law as long as the temperature range is not too large. To extrapolate from laboratory experiments and previous numerical experiments to mantle-like conditions, an understanding of how to relate exponential and Arrhenius viscosity laws is needed. For example, can the onset time for the simpler exponential temperature dependence be applied for Arrhenius-type viscosities? The base of the stagnant or conductive lid, if characterized by a particular viscosity ratio across the convecting part of the thermal boundary layer (see Fig. 1), can be defined by a characteristic temperature difference DTc. eexp and eA can be compared on the basis of the corresponding DTc values that they predict. The eexp value that gives the same DTc as the Arrhenius law, when temperature beneath the convecting boundary layer is still equal to the initial temperature Ti, is given by
 DTi m e* ¼ eexp ¼ eA 1þ ð4Þ eA Ti Here, m is the logarithm of the viscosity ratio l(Ti  DTc)/l(Ti) arising from a temperature difference DTc ¼

m D Ti eexp

ð5Þ

across the convecting thermal boundary. As the temperature dependence of viscosity becomes large (eA ! l), the ratio of eexp to eA approaches the

em

µµ Fig. 1. Temperature (dark) and viscosity (light) as a function of depth just prior to the onset of convection showing conductive lid (dL) and convecting boundary layer (dc). In fluids with strongly temperature-dependent viscosity, convection is limited to a region of low viscosity. The bottom of the stagnant lid (that part of the fluid not participating in convection) is defined by a viscosity increase of about a factor of 10 with respect to the viscosity at depth. This viscosity increase corresponds to a temperature difference, DTc.

constant factor DTi /Ti that is simply a result of how the nondimensional temperature is defined in the two viscosity laws.

3. Physical basis for instability The onset of convection occurs when a suitably defined Rayleigh number reaches a critical value defining the condition of marginal stability in which the amplification of horizontal temperature variations due to buoyant advection balances their damping due to heat conduction. Specifically, the Rayleigh number can be defined as the ratio of the Rayleigh – Taylor (R-T) growth rate (rRT) due to the density variations that drive flow and the rate of conductive smoothing (rCS) of temperature variations. To determine the appropriate form of this Rayleigh number, the R-T growth rate and thermal diffusion rates are calculated individually, as in the recent study of Marquart et al. [13].

374

S.E. Zaranek, E.M. Parmentier / Earth and Planetary Science Letters 224 (2004) 371–386

For temperature-dependent viscosity, we expect rRT ~qo aDTc gdc =li

ð6Þ

where DTc, as defined above, is the temperature difference corresponding to a viscosity increase by a prescribed factor exp(m) relative to li , the viscosity in the uncooled fluid at depth, dc is the thickness of the convecting part of thermal boundary layer as illustrated in Fig. 1. The formula for calculating dc is given in Appendix A. In Eq. (6), the coefficient of thermal expansion a and reference density qo result in a characteristic density difference qoaDTc between the cooled layer and the fluid beneath. Fig. 1 also shows the thickness of the stagnant part of the thermal boundary layer, dL. Conductive smoothing may occur in two ways. In each, the rate of conductive smoothing should be proportional to the thermal diffusivity j divided by the square of an appropriate length scale. First, at a time t when the thermal boundary layer has a thickpffiffiffiffiffiffiffi ness d ¼ pjt , the amplitude of small lateral perturbations in temperature with a wavelength kHd decay exponentially in time at a rate that is proportional to the time for heat to conduct to the top boundary through a distance d rCS ~

j : d2

ð7Þ

rRT qo aDTc gd3c ~ ¼ Racrit : rCS jli

ð9Þ

Davaille and Jaupart [1] used a stability criterion of this form to explain the onset times in their laboratory experiments. However, the results presented in Fig. 2 and discussed below indicate that the conductive cooling to the top boundary through the stagnant lid is dominant in determining the rate of conductive smoothing and thus the formulation of the stability criterion Eq. (8) that matches our numerical experiments. We hypothesize that the appropriate rate of conductive smoothing exists just as instability occurs before plumes have formed. Thus, d rather than dc is the relevant conductive length scale.

4. 2D numerical experiments 4.1. Numerical formulation

An alternative interpretation is that downward propagation of the cooling front at a rate dd/dt overprints growing temperature perturbations. With rCS from Eq. (7) at the time of onset, rRT q0 aDTc gdc d2 ~ ¼ Racrit : rCS jli

the driving force of instability. Once a cool drip has formed, the time for conductive heating depends on its width or thickness which should be proportional to the convecting thermal layer thickness, dc. Since dc < k and shortest length scale should control the rate of conductive heating (see Table 1), dc would then be the appropriate length scale in rCS and

ð8Þ

The timescale for conductive smoothing could alternatively be controlled by horizontal heat conduction over a length scale fk/2. But in our numerical experiments, k>2d (see Table 1), so that d should be the relevant length scale in conductive smoothing rate. Since k is clearly correlated with d(kf3d) in our numerical experiments, the relative importance of the two length scales cannot be directly distinguished. A second way that heat conduction could damp convective motions is by the conductive heating of growing cool downwellings, thus effectively reducing

Conservation equations for mass, momentum, and heat are nondimensionalized using the domain size L as a length scale, a time scale L2/j, a viscosity li, and a temperature scale DTi. For an infinite Prandtl number fluid in which inertia forces are negligible, this nondimensionalization leads to a single dimensionless parameter Rabox. The numerical experiments are formulated with second-order finite difference approximations of these equations on a staggered grid and the Smolarkiewicz [14] method for advection. Solutions for the buoyant viscous flow are obtained using a multigrid iterative solver. The fluid is initially at a uniform temperature Ti with the temperature at the top boundary set to zero (273 K). Numerical experiments are carried out on either 1  1 or 2  2 domains with 256  256 grid points. The bottom boundary is adiabatic and free slip. Horizontally periodic solutions are obtained with

Table 1 Parameter values in numerical experiments Size

Grid

Visc-type

e*

DTc

Onset

d

dc

Pts in dc

k

k/d

k/dc

Racrit

5.00E + 07 5.00E + 07 5.00E + 07 5.00E + 07 5.00E + 07 5.00E + 07 5.00E + 06 5.00E + 06 5.00E + 06 5.00E + 06 5.00E + 06 5.00E + 06 5.00E + 06 5.00E + 06 5.00E + 06 5.00E + 06 5.00E + 06 5.00E + 06 5.00E + 06 5.00E + 06 5.00E + 05 5.00E + 05 5.00E + 05 5.00E + 05 5.00E + 05 5.00E + 05 1.00E + 06 1.00E + 06 1.00E + 06 1.00E + 06 1.00E + 06

1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1

256  256 256  256 256  256 256  256 256  256 256  256 256  256 256  256 256  256 256  256 256  256 256  256 256  256 256  256 256  256 256  256 256  256 256  256 256  256 256  256 256  256 256  256 256  256 256  256 256  256 256  256 256  256 256  256 256  256 256  256 256  256

exp exp exp exp exp exp exp exp exp exp exp exp exp exp exp exp exp exp exp exp exp exp exp exp exp exp a a a a a

3 6 9 12 15 18 3 4 7 8 9 10 11 12 13 14 15 16 17 18 3 6 9 12 15 18 3.636364 4.46281 5.289256 6.115702 6.942149

0.8 0.4 0.266667 0.2 0.16 0.133333 0.8 0.6 0.342857 0.3 0.266667 0.24 0.218182 0.2 0.184615 0.171429 0.16 0.15 0.141176 0.133333 0.8 0.4 0.266667 0.2 0.16 0.133333 0.66 0.537778 0.45375 0.392432 0.345714

0.00087 0.0015 0.002151 0.002718 0.00339 0.00392 0.0038 0.0048 0.00773 0.0091 0.0098 0.0107 0.0125 0.01327 0.01379 0.0149 0.01537 0.01636 0.0169 0.01788 0.01787 0.0317 0.04595 0.05959 0.07362 0.08526 0.0119 0.01589 0.01775 0.020719 0.02195

0.052316 0.068694 0.082265 0.092469 0.10327 0.111049 0.109336 0.122883 0.155942 0.169197 0.175584 0.18347 0.198302 0.204319 0.208283 0.216504 0.219892 0.226863 0.230577 0.237168 0.237102 0.315793 0.380203 0.432971 0.48125 0.517899 0.193484 0.223581 0.236304 0.255304 0.262778

0.043217 0.039155 0.040654 0.04204 0.044339 0.045682 0.090322 0.084598 0.083839 0.086851 0.086771 0.087815 0.092343 0.092891 0.092708 0.094559 0.094412 0.095904 0.0961 0.097563 0.195867 0.18 0.18789 0.196846 0.206627 0.213046 0.140674 0.145373 0.14196 0.144425 0.141703

11 10 10 11 11 12 12 11 11 11 11 11 12 12 12 12 12 12 12 12 25 23 24 25 26 27 36 37 36 37 36

0.166667 0.25 0.25 0.25 0.333333 0.333333 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.666667 0.666667 0.666667 0.666667 0.666667 0.666667 0.666667 1 1 1 1 1 2 0.5 0.5 1 1 1

3.185786 3.639334 3.038975 2.703602 3.227799 3.001672 3.658438 3.255119 2.565062 2.364105 2.278109 2.180197 3.361873 3.262878 3.200768 3.079237 3.031792 2.938628 2.891299 2.810946 4.217599 3.166634 2.630176 2.309622 2.077924 3.861756 2.584188 2.236329 4.231831 3.916906 3.805487

3.856463 6.384852 6.149492 5.946696 7.517764 7.296852 4.428618 4.728257 4.771041 4.605611 4.609847 4.55504 7.219458 7.176848 7.191032 7.050281 7.061251 6.951361 6.937253 6.833211 5.105494 5.555547 5.322269 5.080118 4.839627 9.387656 3.554324 3.439439 7.044238 6.924024 7.057026

4.73E + 03 3.70E + 03 3.67E + 03 3.59E + 03 3.78E + 03 3.76E + 03 4.32E + 03 3.83E + 03 3.50E + 03 3.73E + 03 3.57E + 03 3.55E + 03 3.96E + 03 3.88E + 03 3.71E + 03 3.80E + 03 3.65E + 03 3.70E + 03 3.61E + 03 3.66E + 03 4.40E + 03 3.59E + 03 3.62E + 03 3.69E + 03 3.83E + 03 3.81E + 03 3.48E + 03 3.91E + 03 3.60E + 03 3.69E + 03 3.38E + 03

S.E. Zaranek, E.M. Parmentier / Earth and Planetary Science Letters 224 (2004) 371–386

Rabox

exp: indicates exponential temperature dependence; a: indicates Arrhenius temperature dependence.

375

376

S.E. Zaranek, E.M. Parmentier / Earth and Planetary Science Letters 224 (2004) 371–386

than 10 finite difference grid points are present within the estimated thickness dc at the onset time (as indicated in Table 1). Moresi et al. [15] indicate that reasonable accuracy can be obtained if viscosity variation across a linear element in a finite-element solution cannot be greater than a factor of 2 to 3. In our numerical experiments, the viscosity varies by about a factor of 10 across the convecting thermal boundary layer; so that with more that 10 grid points, the maximum viscosity variation between grid points in this layer is only about 10%. We therefore expect our solutions to be very well resolved. Fig. 2. Onset times of convection from our numerical experiments as a function of the temperature dependence of viscosity (e*). Onset times are scaled by the factor (qoaDTig/lij) 2/3/j predicted by Eq. (8) or Eq. (9). Onset times are well predicted by a constant critical value of the boundary layer Rayleigh number given by Eq. (8) indicating that d is the appropriate length scale for the rate of conductive smoothing.

appropriately defined boundary conditions on the vertical sides of the computational domain. Domains larger than 1  1 at the same Rabox were used to assess the possible effect of domain size on the wavelength of instability. Our numerical experiments are summarized in Table 1. Onset is chosen as the time at which the maximum deviation of the horizontally averaged temperature profile from that due to conductive cooling reaches 1% of the maximum initial temperature Ti. In numerical experiments, the onset time can depend on the form and amplitude of the initial perturbation. Cases with longer onset times allow more time for an initially introduced perturbation to decay by conduction. Simply introducing a perturbation initially could bias the inferred onset times. Therefore, a small perturbation in the temperature field with amplitude 10 6 of the horizontally averaged temperature at each depth and modulated by a random number with zero mean in the range F 0.5 is introduced at each grid point after each 50 time steps in our numerical experiments. Resolution of the numerical experiments is an important consideration. A minimum number of grid points or elements is required within the convecting part of the thermal boundary layer to provide acceptable resolution. In our numerical experiments, no less

4.2. Results from 2D numerical experiments Onset times for 2D numerical experiments for both the Arrhenius and exponential viscosity lie on a single trend as a function of e* (Fig. 2). The solid lines corresponds to a constant value of the boundary layer Rayleigh number as defined in Eqs. (8) and (9). The derivation of the predicted onset time using these criteria is contained in Appendix A. Fig. 2 shows that the onset times in our numerical experiments for a range of Rabox values satisfy the criteria Eq. (8) for a constant value, Racrit. Using the measured onset times and given e* values, values of DTc and of dc at onset can be calculated using equations found in Appendix A. These can then be used to calculate an average Racrit for

Fig. 3. Onset times of convection from the laboratory experiments of Davaille and Jaupart [1] as a function of the temperature dependence of viscosity (e*). Within the range of e* values of experimental data, onset times could be predicted by the thermal boundary Rayleigh number given Eq. (8) or Eq. (9). However, these criteria predict significantly different onset times at higher e* values.

S.E. Zaranek, E.M. Parmentier / Earth and Planetary Science Letters 224 (2004) 371–386

the numerical experiments. Except for one case with a relatively weak viscosity – temperature dependence eexp = 3, our calculated values of Racrit = 3687 with a standard deviation of F129 showing no consistent trend with e* (see Table 1). This Racrit value assumes m = 2.4 as discussed below. 4.3. Comparison with laboratory experiments In the experiments of Davaille and Jaupart [1], the convecting part of the thermal boundary layer (dc) was shown to correspond to a relative viscosity increase of a factor of e2.2f10 or m = 2.2. The convective heat flux after the onset of convection scales with a Rayleigh number of the form Eq. (9) as verified by both laboratory [16] and numerical experiments [17]. The predicted onset time of convection scales with a constant value of a critical Rayleigh number (Racrit) that could be defined by either Eq. (8) or Eq. (9). Either scaling predicts onset times in reasonable agreement with observed onset times over the range of viscosities variations in laboratory experiments (Fig. 3). In fact, the scaling Eq. (9) appears to be a slightly better fit than Eq. (8) over the range of e* in the Davaille and Jaupart laboratory experiments. However, as shown in Fig. 2, this scaling deviates significantly from onset times for numerical experiments with more strongly temperature-dependent viscosity. For thermally activated creep in planetary mantle with e* = eexp = 20 – 35, the criterion Eq. (9) would predict significantly longer onset times than expected from our numerical experiments and the criterion Eq. (8). 4.4. Comparison with other 2D numerical experiments Onset times from numerical experiments with both exponential and Arrhenius viscosities have recently been reported by Korenaga and Jordan [10] and Huang et al. [11]. As shown in Fig. 4, the onset times of Korenaga and Jordan are well predicted by Eq. (8). Their numerical experiments, using a finite element Fig. 4. Onset times of convection from other numerical experiments (symbols) compared with the onset time prediction from Eq. (8). Onset times from Korenaga and Jordan [10] and Huang et al. [11] are in the top two plots and the bottom plot, respectively.

377

method similar to ConMan [18], introduce a random initial temperature perturbations of order 10 3 and 10 5. Changing the initial perturbation from 10 3 to 10 5 increases the onset time by a small and relatively constant factor c1.5. A trend of decreasing normalized onset time with increasing Rabox can be seen in Fig. 4. In other words, with reference to Eq. (8), the Racrit value decreases with increasing Rabox. No such dependence is predicted by Eq. (8) or by the

378

S.E. Zaranek, E.M. Parmentier / Earth and Planetary Science Letters 224 (2004) 371–386

scaling of Korenaga and Jordan. The explanation for this dependence is unclear but could be due to numerical resolution or to numerical diffusion introduced by upwind differencing. Onset times from the numerical experiments of Huang et al. [11] are also well predicted by our stability criterion (Fig. 4). These numerical experiments, using the Citcom finite-element code [15,19], consider an Arrhenius viscosity –temperature dependence and an initial random temperature perturbation of the order of 10 3. Huang et al. describe both the onset times and calculated DTc in terms of power law fits as a function of Rabox and eA tonset ~Ra0:68F0:01 e0:74F0:03 box A

ð10Þ

e0:82F0:06 DTc ~Ra0:01F0:02 box A

ð11Þ

Although our definitions of Racrit and DTc do not intrinsically lead to power laws like Eqs. (10) and (11), power law fits to relationships from (Eqs. (4, 5 and 8) over the same parameter range gives 2=3

tonset ~Rabox e0:70 A

ð12Þ

DTc ~e0:85 A

ð13Þ

in good agreement with Eqs. (10) and (11).

5. Linearized analysis Vertical velocities as a function of viscosity in one of our numerical experiments is shown in Fig. 5a. This figure shows that the relative magnitude of vertical velocities becomes small at any point where l/li>10. Linearized analysis is used to clarify this fundamental observation that convective motion is limited to regions of sufficiently low viscosity and to determine the dependence of rRT on the temperature dependence of viscosity. This approach also shows how the scaling Eq. (6) can be adapted for either Arrhenius or exponential viscosity – temperature dependence. The R-T growth rate of horizontally periodic flow of a given wavelength due to vertical density variations that result from conductive cooling was calculated numerically. Marquart et al. [13] use a similar

approach. However, their results treat a simple twolayer viscosity and density model and therefore cannot be directly compared with our results that have a continuous variation of temperature and viscosity with depth. 5.1. Formulation Equations describing the vertical variation in velocity and the linearized density or temperature perturbations caused by flow are given in Appendix B. Temperature in the basic state is prescribed to follow the conductive cooling distribution with depth (Eq. (A1)). Exponential growth in time of small perturbations from this initially quiescent state are described by eigenvalues and eigenvectors, representing the R-T growth rate of the instability and the depth variation of the velocity, respectively. The eigenvalues and eigenvectors are calculated numerically by a method described in Appendix B. 5.2. Results from R-T growth rate calculations The velocity eigenvectors for a range of eA or eexp values provide a basis to assess the validity of the assumption that flow is confined to regions of sufficiently low viscosity. The velocity eigenvector as a function of depth for a given d (in this case d = 0.4) depend on eA or eexp (Fig. 5c). These eigenvectors, plotted as a function of depth relative to the base of the conductive lid for the appropriate eA or eexp value, are nearly identical (Fig. 5d). Similarly, the velocity eigenvectors plotted as a function of l/li for a range of eA or eexp are nearly identical, and velocities are small in regions where l>10li (Fig. 5b). Velocity as a function of l/li is remarkably similar for a range of e, using both exponential and Arrhenius viscosity laws. The absolute values of vertical velocities from the 2D numerical experiments near the onset of convection are plotted as a function of l/li in Fig. 5a. For both the Arrhenius and exponential viscosities, vertical velocity decays with increasing viscosity, showing very little motion for l>10li and matching the eigenvectors of stability analysis very well (Fig. 5b). For low values of eA or eexp, the exponential and Arrhenius growth rates differ significantly, but this difference decreases with increasing temperature dependence. The growth rates collapse to a single curve

S.E. Zaranek, E.M. Parmentier / Earth and Planetary Science Letters 224 (2004) 371–386

379

Fig. 5. (a) Absolute values of vertical velocities from the 2D numerical experiments as a function of viscosity relative to the interior viscosity (l/ li) at the onset of convection at each grid point in the numerical domain. The velocity distribution is similar to that seen in the vertical velocity eigenvector from the linearized analysis in (c). (b) The same eigenvectors in part (c) as a function of the viscosity normalized by the reference interior viscosity (l/li). Little motion occurs where the ratio of viscosity to interior viscosity is larger than f10. (c) Eigenvectors (vertical velocities) as a function of depth from linearized analysis for a prescribed, constant value of d(d = 0.4) and three different temperature dependences of viscosity (eA = 2,4,6). More strongly temperature-dependent viscosity results in a thicker stagnant lid. (d) The same eigenvectors in part (c) as a function of depth relative the base of the stagnant lid, z  dL. Eigenvalues are almost identical for the range of eA.

when plotted for appropriately rescaled e* values (Fig. 6). In these cases a DTc that corresponds to a viscosity increase of e2.4 (or m = 2.4) provides the best correlation between growth rates. This is very close to the predicted viscosity ratio found by Davaille and Jaupart [1] of e2.24 (f9.3) and is consistent with the vertical velocity eigenvectors describe above. As the temperature dependence of the viscosity increases, both the buoyancy available in the thermal boundary layer qaDTc and dc/d decrease, thus reduc-

ing the R-T growth rate (Fig. 6). For a viscosity contrast exceeding a factor of f10 (e2.4), the growth rate scales with DTc and dc. This scaling is not viable for a viscosity ratio across the thermal boundary layer less than f10 (e2.4). For these viscosity ratios, the growth rate can be scaled by the isoviscous scaling. As e increases from the purely isoviscous case (e = 0) to e = 2.4, the scaled growth rates decrease reflecting higher viscosities within the thermal boundary layer (Fig. 6 inset). When e>2.4, the viscosity variation

380

S.E. Zaranek, E.M. Parmentier / Earth and Planetary Science Letters 224 (2004) 371–386

Fig. 6. Rayleigh Taylor (R-T) growth rate as a function of the temperature dependence of viscosity (e*). Growth rates from Arrhenius (unfilled squares) and exponential temperature-dependent viscosities (solid triangles) are plotted as a function of eexp or eA . When the growth rates are plotted as a function of e* (solid squares and triangles), the curves merge and are well described by the R-T predicted growth rate of Eq. (6). The e* value (see text for definition) represents an appropriate adjusted value of eA or eexp. Inserted graph shows the growth rates a viscosity with exponential temperature-dependent viscosities, before and after scaling by qoaDTgd/li (for e*< 2.4) or qoaDTcgdc/li (for e*z 2.4).

The decrease in rRT by more than a factor of 4 is consistent with a smaller Racrit predicted from our numerical experiments for the isoviscous fluid (1870) than that predicted for cases with temperature-dependent viscosity (3687). Larger rRT would require a smaller Racrit. However, the Racrit increases only by a factor of 2 indicating that the rate of conductive smoothing must also be affected by the presence of a stagnant conductive lid. The rate of conductive smoothing would be expected to decrease in the presence of a conductive lid. Nield [20] calculated the critical Rayleigh number for convection in a constant viscosity fluid layer heated from below beneath a conductive lid of prescribed thickness. In these calculations, the presence of the lid of equal conductivity to the underlying fluid reduces the critical Rayleigh number. This decrease of the critical Rayleigh number is indicative of a decreased rate of conductive smoothing due to heat conduction in the lid consistent with behavior in our numerical experiments described above. In our numerical experiments with a layer cooled only from above and with a conductive lid thickness that changes with time, the form of the critical Rayleigh number is expected to differ from that of Nield.

within the convecting part of the thermal boundary (dc) stays fixed. 6. Implications for planetary lithospheres 5.3. Effect of viscosity distribution and conductive lid on Racrit It has been hypothesized that convection in a fluid with temperature-dependent viscosity can be treated as convection of an isoviscous fluid with a boundary layer of thickness dc and temperature contrast DTc beneath a stagnant, conducting lid of thickness dL. For the isoviscous case, rRT (calculated as described above) scales like qoaDTgd/li. After rRT for the isoviscous and temperature-dependent cases are appropriately scaled, the value of the scaled rRT for the temperature-dependent cases (f0.02) is almost five times lower that of the isoviscous case (f0.094) (Fig. 6 insert). This decrease is because the viscosity distribution in the temperature-dependent case in the convecting thermal boundary layer (dc) is not uniform but varies by a factor of about 10 and with this increased viscosity the rRT decreases.

The stability criterion Eq. (8) agrees with results of numerical experiments to a sufficiently strong viscosity – temperature dependence that extrapolation to the still stronger dependence expected for thermally activated creep in planetary mantles is reasonable. Using this criterion, we explore the conditions that would be needed to explain the thickness of both oceanic and continental lithosphere as the conductive lid due to top-down cooling. For convection beneath cooling oceanic lithosphere, the predicted onset times range from f4 to 80 Myr for value of Q and li from 300 to 550 kJ/mol and 1017 – 1019 Pa s, respectively (Fig. 7). If the 70Ma departure of the oceanic depth – age relationship from that predicted for conductive cooling is caused by the onset of small scale convection (e.g., [1]), the oceanic mantle viscosity (li) derived from the above scaling laws ranges from 4 to 8  1018 Pa s. The

S.E. Zaranek, E.M. Parmentier / Earth and Planetary Science Letters 224 (2004) 371–386

Fig. 7. Conductive lid thickness in km at the onset of convection (top) and the onset time in Myr (bottom) for a range of activation energies ( Q) and reference viscosities (li). Onset time is calculated using stability criterion in Eq. (8) and lid thickness is defined using DTc as defined in the text.

lithosphere or conductive lid thickness (dL) at the onset of convection would be about 110 km (Fig. 7). For values of li c 1017 Pa s, the onset of convection may occur beneath young oceanic lithosphere, as suggested to explain gravity lineations first detected in sea surface altimetry data [21,22]. Oceanic bathymetry may also be explained by models such as the ‘‘Chablis’’ model in which convective heat transfer to the lithosphere occurs at all ages [23]. The results of the current

381

study provide no basis to reject the hypothesis that convective instability may cause both gravity rolls at young ages and a 70-Myr departure from conductive cooling. The presence of melt could be responsible for the low li values required for convective instability beneath young oceanic lithosphere [24]. Perhaps even more importantly, small amounts of water in melt present in the mantle may increase the creep rate of mantle minerals in equilibrium with it [25,26]. The eventual loss of melt and/or water may explain increasing viscosity as the mantle ages. Top-down cooling of initially hot continental mantle may also provide a simple mechanism to create continental lithospheres. However, 150 – 350-kmthick lithosphere beneath old, stable continental interiors [27 – 29] is clearly thicker than the thickest present day oceanic lithosphere. In simple top-down cooling, can the continental lithosphere thickness be reconciled with that of thinner oceanic lithosphere? A lithosphere corresponding to a lid f250 km thick could form by top-down cooling for a mantle viscosity about an order of magnitude greater (3  1019 – 1020 Pa s) than that explaining convective instability beneath 70-Myr-old oceanic lithosphere (Fig. 7). A lid 250 km thick would form after approximately 300 – 400 Myr of conductive cooling. Since the presence of small amounts of water significantly enhance creep rate, higher viscosity is expected for a drier mantle. Analysis of electrical conductivity profiles suggest that Archean continental mantle is in fact drier than oceanic mantle over depths of 150 – 250 km [30,31], possibly reflecting drier residual mantle produced by larger degrees of melting of primitive mantle in the Archean [32,33]. Both the oceanic and continental lithosphere, as residuals of mantle melting during the generation of crust, are expected to be distinct in composition relative to the underlying mantle from which it is derived. Larger degrees of melting in the Archean, than beneath modern oceanic spreading centers, may explain a drier and more creep resistant cratonic mantle. For a strongly temperature-dependent viscosity, no significant thinning of the lid is expected at the onset of convective instability. After convective instability, the lid in our numerical experiments retains a thickness nearly equal to that at the onset of convection (Fig. 8). The lid thickness after the onset of convection is defined as the depth at which the horizontally averaged

382

S.E. Zaranek, E.M. Parmentier / Earth and Planetary Science Letters 224 (2004) 371–386

viscosity increases by a factor exp(m) relative that in the isothermal (adiabatic) region at depth. The lid thickness, thus defined, agrees well with the depth at which vertically advected heat flux becomes small. For a viscosity that is weakly temperature dependent (e* = 9), the lid thins slightly with time after the onset of convection. With more strongly temperature-dependent viscosity (e* = 15), the thermal lid continues to thicken slowly. Slow thickening occurs because the heat flux convected to the base of the lid is smaller than that being conducted through the lid. For less strongly temperature-dependent viscosity, larger DTc results in a larger heat flux to the base of the lid (e.g., [16,38]). For e* = 9, the materially defined lid at onset thins with time. For e* = 15, the materially defined lid retains its thickness at the onset of convection as the thermal lid continues to thicken slowly as cooled mantle accretes to the bottom. Convective instability in the mantle with e* = 20 –36 should reflect this more strongly temperature-dependent behavior. If continental lithosphere were the conductive lid preserved after the onset of convective instability, the lid should retain a composition reflecting the mechanism by which the continental mantle rock formed. Re – Os depletion ages [34] of mantle xenoliths suggests that the minimum age at which continental lithosphere last melted is comparable to the age of overlying crust. This implies the continental mantle sampled by the xenoliths has been stable at least since the crust formed. If processes other than convective instability [35 – 37] did not significantly thin compositionally distinct continental mantle during the 300– 400 Myr of conductive cooling prior to convective instability, the compositionally distinct lithosphere would be equivalent to stagnant lid at the onset convection preserved from that time. Fig. 8. Thickness of conductive lid (shaded), defined from the horizontally averaged viscosity in 2D numerical experiments, as a function of the square root of time for several values of Rabox and e*. Lid thickness and time are normalized by their values at the onset. Contours show the horizontally averaged vertical advected heat flux. Dark dashed straight line corresponds to the lid thickness if purely conductive cooling were to continue. The horizontal light dashed line shows the thickness of the lid at the onset of convection. For cases where the lid continues to thicken, the material cooled and accreted to the bottom of the lid is not the material that was originally beneath it.

S.E. Zaranek, E.M. Parmentier / Earth and Planetary Science Letters 224 (2004) 371–386

On long time scales after the onset of convection, other processes such as heating and erosion by mantle plumes, erosion due to motion of the lithosphere across underlying mantle due to plate motion, or edge-driven convection [35 – 37] may thin an initially thicker compositional lithosphere thickness. If these other mechanisms later thinned the lithosphere, an initially thicker lid at the onset of convective instability could exist. To understand the possible role of top-down cooling in forming both oceanic and continental lithosphere, other aspects of thermally activated creep rheology may be important. For example, the implications of the pressure-dependence viscosity should be considered (e.g., [38,39]). Simple scaling relationships like those discussed above for purely temperature-dependent viscosity could be extended to pressure-dependent viscosity. The possible effect of other, more dynamic aspects of solid-state creep rheology are also interesting to consider. At low stresses, deformation by diffusion creep is strongly grain size dependent. Since equilibrium grain size in laboratory experiments (e.g., [40]) varies inversely with stress or deformation rate, smaller grain size may be present beneath oceanic lithosphere where deformation due to plate spreading reduces grain size and lowers mantle viscosity relative to continental mantle (i.e. [41]). A low deformation rate within continental lithosphere over long periods of time following the last tectonic deformation could allow grain growth leading to an increased viscosity.

7. Conclusion Despite other rheological effects, which may play a still poorly understood role in controlling lithosphere thickness, the effect of temperature-dependent viscosity is clearly fundamental and can be quantified using the simple scaling relationships described in this study. These scaling relationships provide a firm basis to assess convective instability from laboratory creep parameters and available observations or expectations on mantle chemical composition. We find no basis, using laboratory creep laws and the stability criterion Eq. (8), to exclude the possibility of convective instability due

383

to top-down cooling of the oceanic lithosphere at both young and old ages. We further suggest that the thicker cratonic lithosphere can also be created by simple top-down cooling. Thick cratonic lithosphere can be consistent with the thinner lithosphere beneath old oceanic crust if continental lithosphere contains less water than the oceanic upper mantle, as suggested by observations. To account for its greater thickness, the viscosity of the cratonic lithosphere needs to be larger by only a factor of 10.

Acknowledgements This research was funded by grants NAG5-13098 from the NASA Planetary Geology and Geophysics Program and OCE-0242261 from NSF. We thank Anne Davaille for her detailed and useful reviews of our manuscript. We also thank Co-Senior Editor Scott King. We are grateful for the assistance of Jan Hesthaven with the QZ algorithm and for the help of Dayanthie Weeraratne and Lindy Elkins-Tanton in proofreading and editing. [SC]

Appendix A . Derivation of equations for dL, dc, and tonset Definitions of the conductive lid thickness, dL, and the convecting thermal boundary thickness, dc, are shown in Fig. 1. dL is defined as the depth below the cold boundary at which the temperature is Ti  DTc. As discussed in the text DTc ,the temperature difference across the convecting thermal boundary layer is defined as the temperature decrease resulting in a viscosity increase by a factor m. (Eq. (5)). Prior to the onset of convection, conductive cooling leads to the temperature distribution
pffiffiffi  p z T ðzÞ ¼ Ts þ DTi erf ðA1Þ 2 d pffiffiffiffiffiffiffi where d ¼ pjt and Ts is the top boundary temperature. Then
pffiffiffi  p dL DTi  DTc ¼ T ðdL Þ  Ts ¼ DTi erf ðA2Þ 2 d

384

S.E. Zaranek, E.M. Parmentier / Earth and Planetary Science Letters 224 (2004) 371–386

Bux Buz þ ¼0 Bx Bz

and with DTc/DTi = m/e* from Eq. (5) 

m

dL 2 ¼ pffiffiffi erf 1 1  e* d p

ðA3Þ

dL/d is therefore determined by m and e*. The thickness dc is defined by (see Fig. 1) dc ¼

DTc BT Bz jz¼dL

ðA4Þ

dc ¼ 

DT  c  2  DTi p dL d exp  4 d

tonset

ðB4Þ

where q1bq0. We consider solutions of the form ˆ q1 ¼ qðzÞsinkxexpðrtÞ

ðB5Þ

uz ¼ uˆ z ðzÞsinkxexpðrtÞ

ðB6Þ

where r is the R-T growth rate of a horizontally periodic perturbation with wavenumber k = 2p/k. This leads to the following equation


ðA6Þ

which can also be written as 
 1 1 qo aDTi g tonset ¼ pj Racrit li j   m 2
  m 2 2=3 1 exp erf 1 ðA7Þ  e* e* Similarly, an expression for tonset can also be obtained for the criterion Eq. (9).

Appendix B . Linearized Analysis Advection of the density field in a fluid with constant density along with the incompressibility condition and momentum balance for the flow of a very viscous fluid are described by the following equations Bq Bq Bq þ ux þ uz ¼0 Bt Bx Bz

We seek linearized solutions that are a small deviation from a basic or initial state

ðA5Þ

where dL/d is given by Eq. (A3). With the substitution of Eq. (A3) into Eq. (A5), dc/d is seen to be determined by m and e*. Using thepstability ffiffiffiffiffiffiffiffiffiffiffiffiffiffi criterion of Eq. (8), dc/d and DTc, with d ¼ pjtonset and solving for tonset gives 
 1 1 qo aDTc g dc 2=3 ¼ pj Racrit li j d



 Bl 2 Buz B2 l B2 B2  lj uz þ 2 j  uz þ 2 Bz Bz Bx2 Bz2 Bz B2 q ¼g 2 ðB3Þ Bx 4

qðx; z; tÞ ¼ q0 ðzÞ þ q1 ðx; z; tÞ

Solving for dc gives

ðB2Þ

ðB1Þ


  B2 B4 2 Bl B3 2 B k  2k þ k uˆ z uˆ z þ l Bz Bz3 Bz Bz2 Bz4
 1 B2 l B2 ck 2 g 2 ˆ þ þ k ¼  ðB7Þ uˆ z u z l Bz2 Bz2 lr 4

2

where c(z) = Bq0/Bz. Boundary conditions set uˆz = 0 on the top and bottom boundaries at z = 0 and z = L, respectively. The top is a no-slip boundary (Buˆz/ Bz = 0) and the bottom is a free-slip boundary (B2uˆz/ Bz2 = 0). The Eq. (B7) is replaced by a system of four first-order equations in which the derivatives are represented by second-order centered finite difference approximations. This discrete eigenvalue problem is solved numerically using the classical QZ algorithm [42].

References [1] A. Davaille, C. Jaupart, Onset of thermal convection in fluids with temperature-dependent viscosity: application to the oceanic mantle, Geophys. Res. 99 (1994) 19853 – 19866. [2] J.S. Turner, Buoyancy Effects in Fluids, Cambridge Univ. Press, Cambridge, 1973 (368 pp.).

S.E. Zaranek, E.M. Parmentier / Earth and Planetary Science Letters 224 (2004) 371–386 [3] S. Chandrasekar, Hydrodynamic and Hydromagnetic Stability, Clarendon Press, Oxford, 1961 (643 pp.). [4] B.R. Morton, On the equilibrium of a stratified layer of fluid, Q. J. Mech. Appl. Math. X (Pt 4) (1957) 433 – 447. [5] T.D. Foster, Stability of a homogeneous fluid cooled uniformly from above, Phys. Fluids 8 (1965) 1770 – 1774. [6] T.D. Foster, Onset of convection in a layer of fluid cooled from above, Phys. Fluids 11 (1965) 1257 – 1273. [7] P.M. Gresho, R.L. Sani, The stability of a fluid layer subjected to a step change in temperature: transient vs. frozen time analysis, Int. J. Heat Mass Transfer 14 (1971) 207 – 221. [8] D.A. Yuen, L. Fleitout, Stability of the oceanic lithosphere with variable viscosity: an initial-value approach, Phys. Earth Planet Inter. 34 (1984) 173 – 185. [9] C. Jaupart, B. Parsons, Convective instabilities in a variable viscosity fluid cooled from above, Phys. Earth Planet Inter. 39 (1985) 14 – 32. [10] J. Korenaga, T.H. Jordan, Physics of multi-scale convection in the Earth’s mantle, 1, onset of sublithospheric convection, J. Geophys. Res. 108 (2003), doi:10.1029/2002JB001760. [11] J.S. Huang, J. Zhong, J. van Hunen, Controls on sub-lithospheric small-scale convection, J. Geophys. Res. 108 (2003), doi:10.1029/2003JB002456. [12] S.E. Zaranek, E.M. Parmentier, Convective cooling of an initially stably stratified fluid with temperature-dependent viscosity: implications for the role of solid state convection in planetary evolution, J. Geophys. Res. 109 (2003), doi:10.1029/2003JB002462. [13] G. Marquart, H. Schmeling, A. Braun, Small-scale instabilities below the cooling oceanic lithosphere, Geophys. J. Int. 138 (1999) 655 – 666. [14] P.K. Smolarkiewicz, A fully multidimensional positive definite advection scheme with small implicit diffusion, J. Comp. Phys. 54 (1984) 325 – 362. [15] L.N. Moresi, V.S. Solomatov, Numerical investigations of 2D convection withextremely large viscosity variations, Phys. Fluids 7 (1995) 2154 – 2162. [16] A. Davaille, C. Jaupart, Onset of thermal convection in fluids with temperature-dependent viscosity: application to the oceanic mantle, J. Geophys. Res. 99 (1994) 19853 – 19866. [17] O. Grasset, E.M. Parmentier, Thermal convection in a volumetrically heated, infinite Prandtl number fluid with strongly temperature-dependent viscosity: implications for planetary thermal evolution, J. Geophys. Res. 103 (1998) 18171 – 18181. [18] S.D. King, A. Raefsky, H.B. Hager, ConMan: vectorizing a finite element code for incompressible two-dimensional convection in the Earth’s mantle, Phys. Earth Planet Inter. 59 (1990) 195 – 207. [19] S. Zhong, M.T. Zuber, L. Moresi, M. Gurnis, Role of temperature-dependent viscosity and surface plates in spherical shell models of mantle convection, J. Geophys. Res. 105 (2000) 11063 – 11082. [20] D.A. Nield, The Rayleigh – Jeffreys problem with boundary slab of finite conductivity, J. Fluid. Mech. 32 (part 2) (1968) 393 – 398. [21] W.F. Haxby, J.K. Weissel, Evidence for small-scale mantle

[22]

[23]

[24]

[25]

[26]

[27]

[28]

[29]

[30]

[31]

[32]

[33]

[34] [35] [36]

[37]

[38]

385

convection from Seasat altimeter data, J. Geophys. Res. 91 (1986) 3507 – 3520. W.R. Buck, E.M. Parmentier, Convection beneath young oceanic lithosphere: implications for thermal structure and gravity, J. Geophys. Res. 91 (1986) 1961 – 1974. M.P. Doin, L. Fleitout, Thermal evolution of the oceanic lithosphere: an alternative view, Earth Planet. Sci. Lett. 142 (1996) 121 – 136. S. Mei, W. Bai, T. Hiraga, D.L. Kohlstedt, Influence of melt on the creep behavior of olivine – basalt aggregates under hydrous conditions, Earth Planet. Sci. Lett. 201 (2002) 491 – 507. G. Hirth, D.L. Kohlstedt, Water in the oceanic upper mantle: implications for rheology, melt extraction and the evolution of the lithosphere, Earth Planet. Sci. Lett. 144 (1996) 93 – 108. J. Phipps Morgan, The thermodynamics of pressure-release melting of a veined plum-pudding mantle, G3, 2000GC000049 (2001). D.S. Weeraratne, D.W. Forsyth, K.M. Fischer, A. Nyblade, Evidence for an upper mantle plume beneath the Tanzanian craton from Rayleigh wave tomography, J. Geophys. Res. 108 (2003) 1 – 14. R.L. Rudnick, W.F. McDonough, R.J. O’Connell, Thermal structure, thickness, and composition of continental lithosphere, Chem. Geol. 145 (1998) 395 – 411. I.M. Artemievia, W.D. Mooney, Thermal thickness and evolution of Precambrian lithosphere: a global study, J. Geophys. Res. 106 (2001) 16387 – 16414. D. Lizarralde, A. Chave, G. Hirth, A. Schultz, Northeastern Pacific mantle conductivity profile from long-period magnetotelluric sounding using Hawaii to California submarine cable data, J. Geophys. Res. 100 (1995) 17837 – 17854. G. Hirth, R.L. Evans, A.D. Chave, Comparison of continental and oceanic mantle electrical conductivity; is the Archean lithospher dry? G3 2000GC000048 (2000). F.R. Boyd, Compositional distinction between oceanic and cratonic lithosphere, Earth Planet. Sci. Lett. 96 (1989) 15 – 26. C. Herzberg, Phase equilibrium constraints on the formation of cratonic mantle, in: Y. Fei, C.M. Bertka, B.O. Mysen (Eds.), Mantle petrology; field observations and high-pressure experimentation; a tribute to Francis R. (Joe) Boyd, Spec. Publ. Geochem. Soc, 1999 (1989) pp. 241 – 257, Houston, TX. D.G. Pearson, The age of continental roots, Lithos 48 (1999) 171 – 194. N.H. Sleep, Survival of Archean cratonic lithosphere, J. Geophys. Res. 108 (2003) 1 – 29. S.S. Shapiro, B.H. Hager, T.H. Jordan, Stability and dynamics of the continental tectopshere, Lithos 48 (1999) 115 – 133. A. Lenadaric, L.N. Moresi, H. Muhlhaus, The longevity and stability of cratonic lithosphere: insights from numerical simulations of coupled mantle convection and continetal tectonics, J. Geophys. Res. 108 (2003) (in press). M.P. Doin, L. Fleitout, U. Christensen, Mantle convection and stability of depleted and undepleted continental lithosphere, J. Geophys. Res. 102 (1997) 2771 – 2787.

386

S.E. Zaranek, E.M. Parmentier / Earth and Planetary Science Letters 224 (2004) 371–386

[39] C. Dumoulin, M.-P. Dion, L. Fleitout, Numerical simulations of the cooling of an oceanic lithosphere above a convective mantle, Phys. Earth Planet Inter. 125 (2001) 45 – 64. [40] H. Jung, S.-I. Karato, Effects of water on dynamically recrystallized grain-size of olivine, J. Struct. Geol. 23 (1997) 2523 – 2526.

[41] C. E. Hall, and E. M. Parmentier, Influence of grain size evolution on convective instability, G 3 doi:10.1029/ 2002GC000308 (2003). [42] C.B. Moler, G.W. Stewart, An algorithm for generalized matrix eigenvalue problems, SIAM J. Numer. Anal. 10 (2) (1973) 241 – 256.