Time-dependent convection in the Earth's mantle

Physics of the Earth and Planetary Interiors, 36 (1984) 305—327 Elsevier Science Publishers B.V., Amsterdam — Printed in The Netherlands

305

Time-dependent convection in the Earth’s mantle G.T. Jarvis Department of Geology, University of Toronto, Toronto, Ontario, MSS JA1 (Canada) (Received March 2, 1984; revision accepted June 27, 1984) Jarvis, G.T., 1984. Time-dependent convection in the Earth’s mantle. Phys. Earth Planet. Inter., 36: 305—327. The range of Rayleigh numbers accessible to 2-dimensional numerical models of mantle convection has been increased, by almost two orders of magnitude, through the use of high resolution numerical grids with up to 200 intervals in each of the horizontal and vertical directions. The range of Rayleigh numbers studied is believed to encompass conditions which would have existed in a hot Archaean mantle. Although the models are highly idealized, in that they neglect, the compressibility and variations of viscosity in the Earth’s mantle, they suggest that during the Archaean (and earlier), boundary layer instabilities developed more frequently than at present, and did so with a scale and separation comparable to the thickness of the thin thermal boundary layers, rather than to the depth of the convecting layer itself. The rising and sinking plumes which develop from these instabilities drive transient flows which nevertheless have a horizontal scale comparable to the depth of the convecting layer. Thus although surface motions during the Archaean may have been more variable than at present, in terms of speed and direction, they might have occurred on a horizontal scale which is similar to that of present day plate motion. The numerical models indicate that motions at the upper surface of time-dependent convection cells are (1) predominantly controlled by transient sinking plumes, and (2) converge directly towards the source region of the sinking plumes. On the other hand, rising plumes have a more diffuse influence on the upper surface, and surface divergence does not always occur immediately above a rising plume. If oceanic plates on Earth constitute the upper surface of mantle wide convection cells, the more pronounced geophysical signature of subduction zones, compared to oceanic ridges, may be simply due to the fact that buoyancy is most concentrated at the Earth’s surface in the form of sinking plumes. The total surface heat flow is found to fluctuate significantly in the time-dependent convection models. Since estimates of the present-day heat loss from the Earth’s interior are, at best, instantaneous measurements, it is possible that these estimates are not representative (in a time-averaged sense) of the present thermal state of the planetary interior.

1. Introduction Convection within the Earth’s mantle occurs at such high Rayleigh numbers that it is generally believed that the convective motions must be highly time-dependent and irregular (e.g., McKenzie and Richter, 1976). Nevertheless, primarily due to limited numerical resolution, most previous numerical models have been restricted to the study of steady 2-dimensional flow and/or flow at low Rayleigh numbers (e.g., Richter, 1973; Turcotte et aL, 1973; McKenzie et al., 1974; Parmentier et a!., 1976; Lux et al., 1979; Daly, 1980; Houseman, 1983a, b). The appropriate value of the Rayleigh number for convection in the Earth’s mantle depends strongly upon the depth extent of the con0031-9201/84/$03.00

© 1984 Elsevier Science Publishers B.V.

vective circulation, and upon the values of the thermodynamic parameters which are characteristic of this depth extent. Unfortunately, the depth of the mantle circulation which involves plate creation and destruction is not known, and is currently a matter of intense debate. The two main hypotheses are: (1) that convection is confined to the upper 700 km of the mantle (e.g., Turcotte and Oxburgh, 1967, 1972; Richter, 1973; McKenzie et a!., 1974; McKenzie and Weiss, 1975, 1980; Parmentier and Turcotte, 1978; Richter and MeKenzie, 1978; Jacobsen and Wasserburg, 1979; O’Nions et a!., 1979; McKenzie. and O’Nions, 1983); and (2) that convection extends throughout the mantle from the Earth’s surface to the core— mantle boundary (e.g., Hess, 1962; Peltier, 1972;

306

Davies, 1977; O’Connell, 1977; Sharpe and Pe!tier, 1978, 1979; E!sasser et a!., 1979; Peltier, 1980; Olson and Corcos, 1980; O!son, 1981; Peltier and Jarvis, 1982; Allêgre et al., 1984). Geophysical and geochemical data are not yet sufficiently discriminating to rule out, completely, either hypothesis. Consequently convection models based on either of these hypotheses are equally important; comparisons of the predictions of models of these 2 cases may ultimately help to resolve this important issue. The low Ray!eigh number studies, cited above, pertain to convection confined to the upper 700 kiii of the mantle, but not to convection throughout the entire mantle. This paper presents the results of numerical of which mantle are convection at the high Rayleighmodels numbers appropriate to whole mantle convection, both at present, and in the past when mantle temperatures were likely higher than at present. The main purpose here is to emphasize and illustrate the time-dependent nature of convection in the whole mantle context. The convection model employed in this study is, like all previous numerical models of mantle convection, highly idealized with many simplifying assumptions. In particular all of the physical variables which enter the defmition of the Rayleigh number (eq. 9) are assumed to be constant throughout the layer. Consequently the model resuits described below cannot be applied directly to the mantle. As with all such model studies only general aspects of the model’s behaviour under various conditions can be ascribed to convection within the Earth’s mantle. 2. The numerical model The numerical model describes time-dependent 2-dimensional flow in an incompressible constant viscosity fluid with Newtonian rheology and infinite Prandtl number. For such a fluid, the Boussinesq form of the hydrodynarnic equations governing the conservation of mass, momentum and energy, and of the equation of state, may be written as

v u =0 = VP pg + PV2U

0

—

(1) (2)

aT/at =

and

—

v Tu + IV2T+ y

(3)

P = Po(1 aT) (4) where u = (u, 0, w) is the velocity field (with horizontal component u, and vertical component w), P is pressure, p is density, g is the acceleration due to gravity, v is kinematic viscosity, T is ternperature, t is time, K is thermal diffusivity, y is the rate of heating due to internal heat sources, p 0 is the density when T 0°Cand a is the coefficient of thermal expansion (Landau and Lifshitz, 1959). When these equations are nondimensionalized in terms of characteristic units of length, time, ternperature (d,respectively) T = d/V they wheremayV be = 2/v,and i~T,mass and Po’ gthTdin terms of a dimensionless stream function, recast = (0, i~, 0), and temperature as —

V4SP

=

—

aT/ax

(5)

and

aT/at= —v T(v X~)+K’V2T+7’ where ~ = (—

a~p/az,0, a sp/8x) = (v x ~)

(6) (7)

(McKenzie et a!., 1974). All quantities are dimensionless in these last three equations, including the thermal diffusivity K’, and the rate of internal heating y’, which are defined as ic’

=

icv/(gai~Td3)y’

=

yv/(ga~T2d)

(8)

For a given initial temperature field, eqs. 5 and 6 are solved cyclically, using standard space-centred finite difference procedures, to determine the ternporal evolution of T and ‘4’. Zero shear stress (i.e., free) boundary conditions are imposed at the upper and lower surfaces, together with a constant temperature (T = 0) at the upper surface and either a constant temperature (T = i~T)or a constant heat flux at the lower surface. The flow is assumed to be horizontally periodic with a wavelength of A = 2d, where d is the depth of the layer. The numerical scheme employed here is described by Jarvis and Peltier (1982) and similar to that used by McKenzie et a!. (1974) and Jarvis and McKenzie (1980).

307

3. Boundary layer resolution and stability In this section, models are presented for the simplest heating configuration an imposed ternperature difference ~tT across the layer, with no internal heat sources. This is the Bénard configuration for which the Rayleigh number is defined as —

— —

i~

d3

‘~

‘

‘~

19\

RB ~ga T )/~KP) “ ~‘ where the subscript B indicates that this is the Bénard Rayleigh number. The efficiency of convection is measured by the Nusselt number, Nu, which is defined as the ratio of the rate of heat transfer across the convecting layer to that which would occur by conduction alone, When internal heat sources are included the surface heat flux is normally specified as a boundary condition (McKenzie et al., 1974) and the temperature difference ~T becomes a free vanable which has to be determined from the solution. This is often referred to as the Roberts heating configuration (in reference to the initial study with these boundary conditions by Roberts (1967)), and the Rayleigh number in this case is defined in terms of the surface heat flux F as RR

=

(gaPd4)/(Kicv)

(10)

where the subscript R indicates that this is the Roberts Rayleigh number, and K is the thermal conductivity. A new Nusselt number Nu’ is defined for this case by the ratio of the temperature difference across the layer ~7~ond which would be required by a conduction solution to deliver the same surface heat flux .F, to the actual temperature difference ~T which exists when the fluid is convecting. (Since t~Tdecreases as RR increases, the value of Nu’ increases with Rayleigh number.) Jarvis and Peltier (1982) stressed that although similar thermal structures can be produced (in the absence of internal heating) by the imposition of either constant temperature or constant flux boundary conditions, the values of the respective Rayleigh numbers when Nu = Nu’ may be quite different. In particular the value of RR may be one or more orders of magnitude larger than the corresponding value of RB. This distinction is emphasized at the outset here since models with both

types of thermal boundary conditions are included in this paper. Figure 1 shows examples of numerical solutions of the temperature field, obtained with the Bénard heating configuration. Contours of T are shown for various values of RB covering a range of five orders of magnitude. Each set of contours in these figures corresponds to a pair of counter-rotating cells with a cold sinking plume at the centre and hot rising plumes at the outer edges. (Since the outer edges are planes of mirror symmetry only half of each hot plume is seen in these diagrams.) The ratio of the value of RB to the critical value for the onset of convection, R~,is indicated on each frame. In the solution shown in Fig. la, for a Rayleigh number just marginally greater than critical, fluid motion produces only small distortions of the isotherms from horizontal (the pure conduction case). At RB/RC = 10 (Fig. ib) velocities are sufficiently large to advect the cold isotherms downwards and the hot isotherms upwards, generating a broad gravitationally stable central region in which warm fluid overlies cold. As RB is continually increased in passing through frames c to f of Fig. 1, isotherms are swept clear of the interior and concentrated into well defined boundary layers around the edges of each roll. It is evident from this figure that the thickness of the boundary layers decreases with increasing Rayleigh number. It may also be seen, from inspection of Fig. 1, that the horizontal thermal boundary layers are significantly thinner than the vertical boundary layers (i.e., plumes). At the highest Rayleigh number shown here the thickness of the horizontal boundary layers is only a few percent of the layer depth. Since both the stability and heat transfer properties of the entire layer of convecting fluid are determined by the temperatures and velocities within the boundary layers, it is essential that these quantities be determined accurately. In a numerical solution this can only be ensured if the boundary layers are adequately resolved on the numerical mesh. Thus at high Rayleigh numbers, fine meshes must be used; grids with up to 200 vertical intervals were employed for the solutions shown in Fig. 1. As a general rule Moore and Weiss (1973) suggested that at least three grid intervals should span every thermal boundary layer.

308

1031

~O~~IO5

—

C

TEMPERATURE

Fig. 1. Contours of temperature at intervals of i~T/11(where ~T is the difference in temperature between the upper and lower boundaries) for pairs of counter-rotating, steady, 2-dimensional convection rolls. The ratio RB/RC is indicated at the lower left in each frame. Solutions were obtained on the following grids: (a) 24 x 24; (b) 24 x 24; (c) 80 x 80; (d), (e), (1) 200 x 200.

The numerical resolution of the narrow boundary layers seen in Fig. 1 is indicated on the plots shown in Fig. 2 of the mean temperature profiles through the upper boundary layers of model solutions for RB = 1O3R~,104R~and 105R~. On each of these plots the depths of individual grid planes is indicated by the solid circles on the profiles. For the case of RB = iO~ R~,Fig. 2a shows that fourteen vertical grid intervals span the

depth across which the temperature vanes from zero (at the surface) to a maximum (in the temperature overshoot zone). For RB = iO~R~and R~ i05 R~,Fig. 2b and c show that the depth to the maximum temperature in the overshoot zones is spanned by seven and four vertical grid intervals, respectively. Thus all three of these profiles satisfy Moore and Weiss’ (1973) criterion. However, it is clear from Fig. 2 that the boundary layers are

309

3

R8~=10

z

I

3!

8 7)

B c’

I

(RR= 4 X 108)

10 12 I 14

(R~=13 x io

1112

4

Nu~42

16

~16

b

T

118

T~

~-

105’8

z

(RR=8x1o9)

10 14

Nu = 85

~16

C T—~’

Fig. 2. Vertical profiles of horizontally averaged temperature through the upper boundary layer of models for which RB/RC has the values (a) iO~,(b) iO~and (c) i05, as indicated on each frame. The solid circles each indicate the depth of an individual grid plane below the upper surface and the mean temperature at that depth. The corresponding value of Nu, as determined from the numerical solution, and the approximate equivalent value of Rg, are also shown on each frame. Vertical dashed lines indicate the mean temperature of the convecting layer.

better resolved for RB = iO~R~and iO~R~than for the case RB = i0~R~. All the solutions shown in Fig. 1 for RB ~ io~ R~have converged to a steady state. This is only possible if the horizontal boundary layers are stable with respect to all 2-dimensional perturbations, in spite of their steep destabilizing temperature gradients. At high Rayleigh numbers, half of the temperature difference across the convecting fluid

occurs across each of the upper and lower thermal boundary layers. If these each have a thickness of 8, a local Rayleigh number R~can be defined for each boundary layer as R

3)/(2icv) (11) 8 (ga~sT8 which on the basis of boundary layer theory (e.g., Turcotte and Oxburgh, 1967; Roberts, 1979) may

310

be re-written as (12) R8 = (8/d)3RB/2 = cRB/Nu3 where c is a constant (Busse, 1967; McKenzie et al., 1974; Parsons and McKenzie, 1978). Since boundary layer theory also predicts that Nu cr R~’3, it has been concluded in the past (e.g., Busse, 1967; McKenzie et al., 1974) that R~is a constant, independent of RB. Thus if R 8 is sub-critical at low RB, it should remain so at higher values of RB. However, this conclusion depends entirely upon the the boundary layer theory ion validity that Nu of cx R~’3.If the exponent in thispredictpower law relationship differs from 1/3 by a small amount e, the above reasoning leads to a different conclusion. Consider a power law exponent of (1/3 ) so that

The actual value of the power law exponent in eq. 13 may be determined by plotting log(Nu) versus log(RB); the exponent is given by the slope of the best fitting straight line through the points for which RB> 50 R~(Jarvis and Peltier, 1982). Figure 3 is an example of such a plot using the new high resolution solutions presented here. The data plotted in Fig. 3 extend the range of RB one

Nu cx ~

(15).) In other words,

—

(13)

Upon substitution of eq. 13 into (12) we find 3~~ = c’R~ (14) R8 = c’RB/R~ where c’ is a constant. From (14) it is evident that R 8 has a weak dependence on RB. If is positive (i.e., the with powerincreasing law exponent is it <1/3) increase R~and could R8 be will expected that the boundary layers will eventually become unstable. On the other hand if is negative the boundary layers should become more stable with increasing RB. 103

~IL _________________________________ ~Ø2 10_____________________

1)

RPJR~—

Fig. 3. Plots of Nu, and 6’( = 8/d), as functions of RB/RC. ~ values are plotted on logarithmic scales. 8~ represents the shallowest depth at which the mean temperature i~T/2occurs in the mean temperature profile, while S~represents the mean depth of the maximum temperature in the upper thermal boundary layer. The straight lines have been fit through the model predictions for all models with 100 ~ ~ 30000. The data included in these graphs are tabulated in Table I in the Appendix....

order of magnitude beyond that studied by Jarvis and Peltier (1982) and were obtained with twice the resolution. From these data it is found that 318 (is) Nu = 2.24 (RB/R~)° in close agreement with the results of Jarvis and Peltier (1982). (Schubert and Anderson (1984) have recently corroborated this result, finding values of 2.24 and 0.319 for the coefficient and exponent in dependence of R

=

0.015, indicating a weak

8 on RB, and boundary layer instabilities should develop at sufficiently high values of RB. A more direct, although somewhat uncertain, indication of boundary layer instability at higher values of RB, was obtained5 from a this convection R~.For model, model for which RB = 5 x i0 initially steady rolls broke up into a highly disordered and time-dependent flow consisting of hot and cold sheets of fluid frequently breaking away from the lower and upper thermal boundary layers at random times and locations. Figure 4a shows a “snap-shot” of the time-dependent temperature field at one instant, while Fig. 4b indicates the resolution of the mean temperature profile in the upper boundary layer. The horizontal boundary layers are barely discernable in Fig. 4a since the

an grid individual spacing iscontour approximately line. Theequal plumes to the here width have ofa distinctly pulsating character with successive parcels of hot (or cold) fluid breaking away from the same location on the boundary layers. When these parcels of anomalously hot (or cold) material reach the other surface they are stretched out horizontally below (or above) the opposite thermal boundary layer. This produces the temperature overshoots with which we are familiar in steady convection models, but here only in isolated patches. Figure 4b indicates that three grid intervals

311

E~O~J1JJJJJ

__

TIME-

I

0.06

48x48

16

1

8 12

Z

k

____

0

100

TIME Nu —140

14

16

T ~ Fig. 4. Temperature solution for a convection model with RB/RC = 5 x i0~,obtained on a 200x 200 numerical grid. (a) Isotherms for one wavelength of the periodic flow at one instant of time. (b) Vertical profile of the mean temperature in the upper portion of the convecting layer. Solid circles indicate the depths of individual grid planes. The vertical dashed line indicates the mean temperature of the convecting layer. span the thermal boundary layers of the model shown in Fig. 4a. Thus the resolution of the numerical grid appears to satisfy the condition which Moore and Weiss (1973) suggested is necessary to obtain solutions which are accurate to within 1 or 2%. However, this solution is at the limits of numerical resolution on a 200 x 200 urnform numerical mesh, and since the thickness of the boundary layer vanes significantly in such a time-dependent flow the boundary layer must be under-resolved at some locations. Therefore considerable reservation must be attached to the prediction of boundary layer instabilities at high values of RB on the basis of this model. Previous experience with the transition of model solutions from a steady state regime at low Rayleigh numbers to time-dependent flow at high

Fig. 5. Two examples of the dimensionless variations of kinetic energy, E~,with time, for the same numerical model of compressible convection. (a) Model solution obtained on a 24 X 24 finite difference mesh. (b) Model solution obtained on a 48 x 48 numerical grid, after interpolating the solution from 5(a) onto a 48 x 48 grid. One overturn time is indicated by the length of the horizontal bar in 5(b). (After Jarvis and McKenzie, 1980.)

Rayleigh numbers, has shown that the value of the transition Rayleigh number will be underestimated by models which are inadequately resolved. For example, Fig. 5 (from Jarvis and McKenzie, 1980) shows the kinetic energy as a function of time for a model of compressible convection as obtained first on a 24 X 24 numerical grid (Fig. 5a), and then on a 48 X 48 grid (Fig. 5b). On the coarse grid an apparently regular oscillatory solution was obtained. (The amplitude of the oscillation gradually diminished with time.) After 316 dimensionless time units the coarse grid solution was interpolated onto the fine grid and then allowed to evolve as shown in Fig. Sb. The oscillations died out after only a few overturns, and a steady one-roll solution evolved. The difference in behaviour of the same model on the different grids was attnbuted by Jarvis and McKenzie (1980) to the poorer representation on the coarse grid of short wavelength instabilities developing within the thermal

312

boundary layer. When the Rayleigh number was increased by a relatively small factor of 2, Jarvis and McKenzie (1980) obtained well resolved time-dependent solutions on the 48 x 48 grid. Thus although the time dependence indicated in Fig. Sa was found to be artificial, it did appear to signal

contrast, theoretical studies on the stability of 2-dimensional convection between free boundaries (Straus, 1972; Skilbeck and McKenzie, 1979) have found rolls to be stable, within a limited wave number range, up to the maximum value of RB considered (278 R0). This result is consistent with

the imminence of boundary layer instability at slightly higher Rayleigh numbers. By way of analogy with the results shown in Fig. 5, it is speculated here that 5although dependence found R, maythe betime premature in terms of at RB = 5 x i0 the value of R~at which it occurs, it may nevertheless indicate a decreasing stability of the boundary layers as RB increases. This would be consistent with the discussion above with regard to eq. 14. Additional uncertainty about this issue has resulted from a recent model calculation by Schubert and Anderson (1984) which was designed to test the reproducability of the time dependence found here at RB = 5 X i0~R~. Schubert and Anderson (1984) used a newly developed finite element program (as opposed to the finite difference method used here) to solve the convection equations at RB = 5 x iO~R~.A steady state solution was obtained by their numerical scheme. Unfortunately there is no other independent model calculation (at this high value of RB) to corroborate either the time-dependence shown in Fig. 4 or the steady state found by Schubert and Anderson. At face value, however, Schubert and Anderson’s result would support the above speculation that the time dependence shown in Fig. 4 is due to inadequate numerical resolution. The question of whether or not this time dependence signals a decreasing boundary layer stability, and inevitable time dependence at some higher Rayleigh number, remains unanswered. This question is of intrinsic fluid dynamic interest since all of the models discussed above were obtained with zero-shear-stress mechanical boundary conditions (i.e., between free boundanes). Both theoretical (Busse, 1967) and experimental (Krishnamurti, 1970, 1973) studies on the stability of 2-dimensional convection rolls between rigid boundaries have indicated that rolls are only stable within a closed region in Rayleigh

the early boundary layer arguments discussed above and seemed to allow the possibility that steady rolls existed at all values of RB. The new boundary argument that this islayer not the case. presented above suggests Although the stability of ideal 2-dimensional Bénard convection rolls is of fluid dynamic interest, it is not particularly relevant to the nature of convection in the Earth’s mantle. All of the steady 2-dimensional solutions discussed above are artificially stabilized. The steady temperature fields at low Rayleigh numbers are pointwise symmetric about the mid-point of each cell (Moore and Weiss, 1973). Uniformly scaled versions of these were used as the initial fields for the higher Rayleigh number models. This together with the symmetric thermal and mechanical upper and lower boundary conditions, and lack of internal heat sources, ensures that any developing perturbation will also be symmetric about the mid-point of the cell. Thus all disturbances which do not possess this high degree of symmetry are eliminated from the mathematical solution a priori. The simplifying assumption of a constant viscosity also contributes to the stability of the model solutions. If viscosity were allowed to vary with temperature, the upper and lower boundaries would have higher and lower viscosities, respectively, than that of the interior of the cell, thus breaking the symmetry about the midpoint. Moreover, the lower viscosity at the bottom of the layer would have a destabilizing influence on the lower thermal boundary layer. In addition, Jarvis and McKenzie (1980) showed that the assumption of incompressibility has a strongly stabilizing effect on the mathematical solutions. Since none of the restrictive model assumptions discussed here actually pertain to the Earth’s mantle, it is expected that mantle convection is time dependent at Rayleigh numbers much less than i05 R 0. In the following section the symmetry of the boundary conditions is broken by imposing a constant heat flux boundary condition at the lower

number—wave number space. The maximum value of RB for which rolls are stable is RB = 13 R~.In

313

boundary and a constant temperature at the upper boundary. Small amounts of internal heating are also included. At lower equivalent values of R P.’ boundary layers are wider and instabilities develop on a larger scale. By studying such models we may see more clearly the physical processes responsible for the time dependence of the flow,

internal heat sources will be denoted ~t. As either of R R or ~sis increased, the flow will be shown to rapidly increase in complexity. MODEL 1: RR = 5 x 108; ~t = 0.1

4. Time dependent model solutions

A sequence of contour plots of temperature and stream function is shown at three successive times in each of Fig. 6a—e, for a model with RR = 5 x 10~ and p. = 0.1 (i.e., internal heating contributes 10% to the total surface heat flow). The fields shown in

Examples of time-dependent flow are presented in this section, beginning with low values of RR and small amounts of internal heating. The fraction of the total surface heat flux which is due to

each of frames a—e in Fig. 6 follow sequentially from one frame to the next. The contour intervals are the same for each frame and, therefore, the number and spacing of streamlines indicates the relative velocities. Where streamlines are closely

T~I~T~

RR=5.E8

JO.1

Fig. 6. Time-dependent convection, model 1: RR = 5 X i0~ ~t = 0.1; 150 X 150 numerical mesh. Contours of T and ~ are shown at fifteen successive times. Dimensionless times relative to the left-hand frame in 6(a) are: (a) 0.00, 0.054, 0.091, (b) 0.103, 0.114, 0.131, (c) 0.161, 0.202, 0.207, (d) 0.230, 0.252, 0.275, (e) 0.297, 0.320, 0.342. Contour intervals are the same in each frame: dT, 100.5; d~, 0.0837. The labels A, B and C on the temperature field refer to features discussed in the text.

314

~

~ULI

~

I

315

~

D

______________ I_

____________I

_______

_______

_______E

316

spaced, velocities are large where they are widely spaced, velocities are small. When there are more streamlines in total, the overall kinetic energy is higher. Beginning with the left-hand frame in Fig. 6a, attention is drawn to the cold parcel of fluid in the upper left corner (labelled A). This parcel originated in the cold upper boundary layer one convective cycle previously: fluid from the cold upper boundary was swept down at the right and across the bottom immediately above the hot lower boundary. It was then drawn viscously upwards by the hot rising plume, and spread out horizontally below the spreading hot plume material, which —

was itself drawn out horizontally below the cold upper boundary layer. As this cold parcel (A) is swept horizontally, in the time sequence shown in Fig. 6a, it begins to sink somewhat into the interior of the convecting layer. When it reaches the downgoing stream at the right-hand vertical boundary, its negative buoyancy reinforces that of the cold descending plume and velocities (as mdicated by the stream function) increase. The cold parcel sinks quickly to the bottom of the fluid layer and draws some of the hot material (which originated in the hot plume at the left-hand vertical boundary) down with it. The time sequence continues in the left-hand frame of Fig. 6b. Here the cold parcel has reached the bottom of the layer and velocities are a maximum. Moving on, the cold parcel rounds the corner and spreads across the bottom, drawing a hot parcel of fluid (labelled B) behind it. This hot parcel lies above the cold spreading stream of the descending plume which itself overlies the hot lower thermal boundary layer. The streamlines indicate that velocities are decreasing. In terms of potential energy, the hot parcel B in the bottom right corner is analogous to the cold parcel A in the upper left corner in the first frame of Fig. 6a. Continuing the time sequence in Fig. 6c the cold parcel A is beginning to be drawn upwards again and a new cold parcel (labelled C) develops at the top left. Again this

cold parcel was formed from cold descending plume material which, due to insufficient heating by basal heat flow, did not lose its thermal identity before reaching the ascending stream adjacent to the hot plume. Velocities are low where cold material is being lifted but higher where the hot parcel reinforces the buoyancy of the upwelling plume. In the far right hand frame of Fig. 6c the original cold parcel (A) is being stretched as out as the lower half of it begins to sink behind the hot rising parcel (B). As the time sequence continues in Fig. 6d, the original cold parcel A loses its thermal identity, and the second cold parcel, C, begins a similar traverse across the top of the cell followed by the hot parcel, B. In the final sequence

for this model, shown in Fig. 6e, the new cold parcel first sinks into the middle of the cell and then is drawn upwards with the rising stream at the right. In addition instabilities develop in the cold upper boundary layer and are swept downwards with the main flow, increasing velocities once more. The example of time dependent flow presented in Fig. 6 illustrates that the overall flow pattern is determined primarily by the integrated distribution of buoyancy, so that (locally) hot parcels may be drawn downwards and cold parcels may be drawn upwards by the flow. Time dependence in convection cells due to circulating hot and cold parcels, has been found previously in low Prandtl number fluids by Moore and Weiss (1973). These authors attributed this phenomenon to the advection of vorticity out of the vertical plumes (where it is generated) into the horizontal boundary layers, and argued further that a similar phenomenon could not occur in infinite Prandtl number fluids such as those considered in the present study~The reason for this is that vorticity cannot be advected in an infinite Prandtl number fluid since the terms describing this process are eliminated from the momentum equation at the outset. Nevertheless circulating hot and cold parcels are clearly present in the model

Fig. 7. Time-dependent convection, model 2: RR =5 x iOu; ~s= 0.2; 150 x 150 numerical mesh. Contours of T and ~j’are shown at six successive times. Dimensionless times relative to the left-hand frame in ~(a) are: (a) 0.00,0.015, 0.035, (b) 0.047, 0.065, 0.085. Contour intervals are the same in each frame: dT, 167.4; d,Ji, 0.0837. The labels A—I on the isotherms refer to features discussed in the text.

317

~

RR=5.E8 ____________________

~

p=O.2 —

318

solution shown in Fig. 6. In this case these are associated with the pulsating nature of the flow velocity, and the consequent failure of thermal overshoot zones to be eliminated by the time they reach the opposite side of the convection cell. (In other words the balance between horizontal advection and vertical diffusion, which is often assumed for steady convection cells, is not maintained.) MODEL 2: RR

=

5 x 10~p. = 0.2

The time sequence shown in Fig. 7 also mdicates that once a hot parcel of material (e.g., B or C) is drawn into the corner next to a cold sinking plume, it resists being drawn downward because of its own positive buoyancy. Notice in the centre frame of Fig. 7b, that streamlines are diverted around the hot corner material (B). In the final frame, the hot corner material is trapped by the newly developing cold plume and forced downward below it.

In the second model considered here, the fractional amount of internal heating is increased to 20% while the Rayleigh number remains the same as in model 1. As a result, the flow becomes considerably more unsteady than in the previous example. Figure 7 shows a typical time sequence for this model. In the sequence of Fig. 7a, circulating hot and cold parcels are evident: in the upper regions the cold parcel, labelled A, begins to sink into the interior as before, and the hot parcel, labelled B, burrows into the cold upper boundary layer. The uneven heating of the base of the cold boundary layer generates short wavelength horizontal temperature variations which trigger boundary layer instabilities of the same wavelength. This process is indicated more explicitly by the two cold parcels of fluid, labelled D and F, which arrive at the bottom boundary (close t~ the descending plume) and burrow into the hot lower thermal boundary layer. At the bottom left, the base of a young hot plume (G) is swept into the main rising stream and lifts a parcel of trapped cold fluid (H) above it. The sequence continues in Fig. 7b where the two cold parcels (D and, F)~ succeed m tnggenng new instabihties in -the hot boundary layer (e.g., J). The contours of stream function indicate that a considerable change in kinetic energy occurs between the various time frames. Once the hot and cold pulses have reached the opposite boundary layer and begin to move horizontally, the kinetic energy drops. When a new parcel begins to rise or descend, the kinetic energy suddenly increases again. .

9

MODEL 3: RR = 1.5 X 10 ; p. = 0.0 When RR is raised to 1.5 x i09 the flow is time dependent even without internal heat sources (i.e., p. = 0). At this Rayleigh number plumes develop rapidly enough to develop fully before being swept towards pre-existing ascending or descending streams. Typical time sequences for this model are shown in Fig. 8a and b. In the first frame (i.e., the left frame) in Fig. 8a, the head of a mature hot plume (labelled A) has risen almost to the surface while its source region (labelled B) has been swept into the major upwelling at the left-hand vertical boundary. In the next two frames of Fig. 8a, the source region (B) rises rapidly at the left, while a new cold plume develops from the upper thermal boundary layer. Of particular interest here is the manner in which the hot rising parcel (A) is drawn in towards the base of the cold plume (C) by the flow which is generated by the sinking plume. The same process is responsible for the cold parcel at the bottom left and hot parcel at the top right of each frame. The presence of these hot and cold parcels tends to diminish vertical velocity in these areas, and thus facilitates the development of mature plumes at intermediate points in the boundary layers. This drawing in of material towards the site of a descending plume is similar to the flow required in various kinematic models of back arc spreading (e.g., McKenzie, 1969; Andrews and Sleep, 1974) and/ or continental extension above subduction zones (Jowett and Jarvis, 1984). The time sequence

Fig. 8. Time-dependent convection, model 3: RR = 1.5 x i09 ~ = 0.0; 150 x 150 numerical mesh. Contours of T and ~pare shown for three successive times in 8(a), and for another three successive times (somewhat later in the flow solution) in 8(b). Dimensionless times relative to the left.hand frame in 8(a) are: (a) 0.00, 0.0095, 0.0184, (b) 0.0976, 0.1060, 0.1143. Contour intervals are the same in each frame: dT, 669.7; thl, 0.167. The labels A—E on the isotherms refer to features discussed in the text.

319

~

___

A

RR=l.5 E9

~

p=O.O

320

of Fig. 8a is an example of this process in the context of a dynamic model. A second sequence of model solutions for model 3, but at a later time, is shown in Fig. 8b. The main point of interest here is the fact that closely spaced plumes tend to migrate towards each other and merge. By the time of the left-hand frame a major cold sinking plume (labelled A and B) has descended almost to the bottom of the layer, and a secondary cold plume (labelled C) has begun to develop nearby. A hot rising parcel of fluid (labelled D) is sandwiched between the two sinking plumes. By the time of the middle frame in Fig. 8b, the secondary plume (C) is well developed and the streamer (B) from the initial cold sinking

plume can be traced back to a point where it emerges from the side of the secondary plume. The hot parcel, D, has been forced downwards between the two cold plumes. By the final frame in this sequence, there is little indication that the major descending plume, C, was once two separate plumes. Only the cold tail, B, of the original sinking plume and the hot parcel, D, remain as tracers of the merger. This sequence also provides another example of an instability of the cold upper boundary layer being triggered by the arrival of a hot parcel of fluid (labelled E). Comparing the final two frames it is evident that the newly forming cold plume (labelled F) is converging towards the neighbouring plume, C.

~

~

_____ IJ=O.1

RR=l.5E9

Fig. 9. Time-dependent convection, model 4: RR = 1.5 X iou; ~ = 0.1; 150 X 150 numerical mesh. Contours of T are shown for six successive times. Dimensionless times relative to the upper left-hand frame are: 0.00, 0.0056, 0.0111, 0.0166, 0.2221, 0.276. Contour intervals are the same in each frame: dT, 837.1.

321

MODELS 4 and 5:

RR

=

1.5 x

io~

p. = 0 1 and p. = 0 2 Although the previous model did not include internal heat sources the model solutions exhibited time dependence due to the formation of hot and cold plumes in the thermal boundary layers. Similar boundary layer instabilities occur in model 4, in which 10% of the surface heat flow is due to internal heat sources (p. = 0.1). However, as illustrated in Fig. 9, many of these plumes are swept to the vertical boundaries before developing fully. When internal heating is included to the extent that it contributes 20% to the surface heat flow, at the same Rayleigh number (model 5), intermediate plumes mushroom continuously on both

boundaries. This leads to inevitable collisions of opposite buoyancy and consequently to aplumes much of more contorted flow pattern; see Fig. 10.

5. Discussion and conclusions Owing to limited computer resources, this preliminary survey of time dependence in high Rayleigh number convection rolls has been restricted to models with small amounts of internal heating. For larger values of p. the flow will clearly become more irregular and unsteady than that which has been illustrated above. Although larger values of p. are likely to be more appropriate to the mantle, these initial studies at low values of p.

RR=l.5 E9 Time-dependent convection, model 5: RR = 1.5 X iou; p 0.2;

p=O.2

Fig. 10. = 150 X 150 numerical mesh. Contours of T and 4’ are shown for three successive times. Dimensionless times relative the left-hand frame are: 0.00, 0.0115, 0.0230. Contour intervals are the same in each frame: dT, 502.3; d4,, 0.0837.

322

should lead to a better understanding of the physical mechanisms which are responsible for the time dependence. Future calculations at higher values of p. may then be interpreted in the light of the simpler flow structures presented here. The highest value of RR considered (1.5 X 10~) is at least an order of magnitude higher than values appropriate to whole mantle convection today (Jarvis and Peltier, 1982). This value of RR may be more appropriate to the whole mantle during the Archaean when mantle temperatures were higher (Green, 1975; Bickle, 1978; Jarvis and Campbell, 1983). Although the value of p. appropriate to the whole mantle today is uncertain (e.g., Jarvis and Peltier, 1981), it must have been higher in the past when radioelements were more plentiful, and less concentration of these elements into continental crust had taken place. Since the degree of unsteadiness increases rapidly with both R R and p., mantle convection during Archaean times was undoubtedly more highly time dependent than at present. It is expected that this time dependence would be further augmented by spatial variations of viscosity, In an earlier study of time dependence in internally heated fluids at RR = 1.4 x 106, McKenzie et al. (1974) noted that the plumes originating in the boundary layers had spatial scales comparable to the depth of the convecting layer. However, when RR = 1.5 x i0~, Figs. 8—10 indicate that plumes have a horizontal scale which is much smaller than the depth of the layer. Nevertheless, the horizontal scale of the flow induced by the plumes remains comparable with the depth of the convecting layer. A general implication of this is that during the Archaean, surface motions may have been more variable than at present in terms of speed and direction, but may have occurred with a horizontal scale similar to that at present. This implication is corroborated, albeit in a qualitative sense, by the time dependent model solution found for Bénard convection at RB = 5 x i05 R~(which is roughly equivalent to RR = 7 x 10’°).As discussed above it is not clear whether the boundary layer instabilities in this model are of numerical or physical origin. Consequently this model is of no quantitative value. However, regardless of their mode of origin, once instabilities

in the boundary layers develop into macroscopic plumes, these plumes are well resolved and behave physically. Because of the high Rayleigh number, the boundary layers and plumes are thinner than those in the model solutions shown in Figs. 6—10. This model may be thought of as a kinematic model in which plumes are artificially generated at random locations on the boundaries. The qualitative nature of the flow produced by these plumes may then be observed. Figure 11 shows a time sequence for this kinematic model. Since the spatial resolution is equal to the width of the contour lines in this figure, it is evident that the rising and sinking parcels and plumes are well resolved. In this sequence, three successive cold parcels of fluid (labelled A, B and C), which have broken away from previous upper boundary layer plumes, settle above the hot lower boundary layer, spread horizontally and control the spacing of new hot plumes (labelled D, E and F). This behaviour mimics that of genuine plumes of physical origin. The streamlines in Fig. 11 indicate that the plume induced flow organizes itself on a horizontal scale which is comparable to the depth of the convecting layer, even though the plume scale is an order of magnitude smaller. It is important to remember, however, that these flow solutions were obtained for a constant viscosity fluid. In the Earth, hot rising plumes would have a lower viscosity than the bulk of the mantle and would therefore exert less viscous stress on the surrounding mantle than would the cold sinking plumes. Although this asymmetry in the viscous forces driving the circulation would modify the flow field in detail, it is expected that since the highly viscous sinking plumes traverse the entire depth of the layer, the component of flow produced by them would be organized on a horizontal scale comparable to the depth of the layer. In general, all of the time-dependent model solutions suggest that plumes have the most direct effect on the boundary layer from which they originate. This applies both to the velocity and temperature fields. By the time the head of a rising or sinking plume reaches the opposite boundary, the thermal contrast with its surroundings is usually considerably diminished. Consequently velocities at the upper surface are controlled primarily

323

-

~

~

~

*0~•~

~

R8

5 Rc

5x i0 5 x 10~R~200 x 200 numerical

Fig. 11. Time-dependent convection (?): RB = mesh. Contours of T and 4’ are shown for three successive times. Dimensionless times relative the left-hand frame are: 0.00, 0.0016, 0.0042. Contour intervals are the same in each frame: dT, ~T/11; d’~i,0.335. Labels A—F on the isotherms indicate features discussed in the text.

by transient sinking plumes, whereas rising plumes are of primary influence on velocities at greater depths. The influence of rising plumes on the velocity of the upper surface would be further reduced in the Earth’s mantle due to the reduced viscosity of the hot plumes relative to that of the cold plumes. Since sinking plumes begin at the upper surface, surface material converges directly to the location of the sinking plume. In contrast surface divergence above a hot rising plume is not always directly above the plume (see Fig. 8, for example). If the zones of surface convergence and divergence in the simple numerical models are the analogues of subduction zones and mid-ocean ridges on Earth, then the more pronounced geophysical signature of trenches relative to ridges is

understandable. Trenches are produced by concentrated negative buoyancy and ridges by diffuse positive buoyancy. Moreover, the migration along the boundaries of the plumes’ source regions, due to horizontal velocities within the boundary layers, may be analogous to the migration of ridges and trenches on Earth. It must be stressed, however, that since simple numerical models cannot reproduce the complex behaviour of the Earth’s mantle, direct analogies between model predictions and the Earth’s behaviour must be regarded with some skepticism. Nevertheless it is worth noting that the major surface plate on Earth, the north Pacific, has a number of characteristics similar to the time-dependent models. For example, it has both a relatively high surface velocity and a large number of subduction zones. Also, although the transit

324

time from ridge to trench represents only about one-quarter of a convective overturn time, the direction (as indicated by the bend in the motion of this plate experienced a major change in Hawaii—Emperor seamount chain) before this transit was more than half completed. In other words, the convection pattern changed significantly on a time scale which is short compared to the overturn time (i.e., 1/8 of an overturn time

1 H

—

—

—

:~

A 10

F f4o~~

0 80f

in this case). (Evidence summarized by Goodwin (1984) that the Atlantic Ocean has opened and

scale Ma, closed 1350 upon approximately is which a further major indication every shifts 400 ofof Ma the thefor convection short the time past patterns can occur.) Finally, the major area of global seismic energy release on Earth is the irregular shaped set of subduction zones near the Tonga and Fiji Islands (e.g., Gutenburg and Richter, 1954; Bullen, 1963). The geometry of these is not compatible with steady plate motion and subduction for the past 108 y. Non-steady convective flows are associated with variable velocities and heat flow at the upper surface. This variability is illustrated in Fig. 12 where dimensionless values of the surface heat flux, .F~total kinetic energy, E~ maximum horizontal velocity, U~ and maximum vertical velocity, W,~ as predicted by model 3 (discussed above), are each plotted as a function of time. Arrows on the time axis indicate times corresponding to the different frames of Fig. 8a and b. In the relatively short interval shown in Fig. 12, the surface heat flow fluctuates between dimensionless values of 33.2 and 53.1, with a mean value of 40. This range of variability corresponds to 50% of the mean value. The kinetic energy (which may be thought of as a proxy for the square of the mean convective velocity) varies in a pulsating manner from a base level of 100 up to a short lived peak value of 1050. The maximum horizontal and vertical velocities, U,~,ax and W~ax vary (in manner similar to E’k) by a factor of 3 from base level to peak value. The time series for ~ and W~axhave more structure than that for E~since these represent local rather than integral fluctuations. If the motion of surface plates on Earth is related to mantle convective velocities, then Fig. 12 suggests that episodic surface motions

+ I I

80

~

4~

40

I

1

-~

.2

Fig. 12. Time sequence of the dimensionless surface heat flux, F~’;total kinetic energy, Ei,; maximum horizontally velocity, U,~and maximum vertical velocity W,~,a*,for model 3 (see Fig. 8).

are a natural consequence of time dependent flow in the mantle. Comparison of the 4 curves plotted in Fig. 12 reveals little phase correlation between the various parameters of the circulation. The best correlation is between U~ and W~ although at times U~ is a minimum while W,~is a maximum. Similarly, P’ may be at a local minimum while E~is at a maximum, and vice versa. These asynchronous responses occur because the time dependence results from the erratic rising and smIting of plumes. W’max’ and hence E~,is a maximum as a plume (or plumes) leaves its boundary layer of origin. However, at this stage a rising plume would not yet have reached the upper surface to produce the anomalously high (near surface) temperatures responsible for the increase in .I~’.By the time a hot rising plume reaches the upper surface its vertical velocity is significantly reduced. The resultant increase in 1~’may therefore be accompanied by a decrease in E~(unless another plume is concurrently breaking away from either of the boundary layers). Figure 8b indicates that such a sequence of events accounts for the increase in F’s, which is labelled A, and the corresponding decrease in E’k, at the end of the time series shown in Fig. 12. The

325

earlier increase in E~,which is also labelled A,

namic states, and may not deviate from the tem-

corresponds to the release from the lower boundary

poral mean as seriously as Fig. 12 suggests. It is

layer of hot plumes which subsequently produce the increase in F;’ (labelled A). The possibility of a higher than average value of F;’ coexisting with a lower than average value of Ej (and hence of U,~,ax and W~ax)is a significant departure from the predictions of steady state convection models. For steady convection cells, F;’, E~,U,~naxand W~axall increase together according to well defined power law dependencies upon the appropriate Rayleigh number (Turcotte and Oxburgh, 1967; McKenzie et al., 1974; Roberts, 1979; Jarvis and Peltier, 1982). These relationships have been used by Sharpe and Peltier (1979), Elsasser et al. (1979) and Peltier (1980) to estimate the depth of the mantle convective circulation. They are also fundamental to parameterized thermal history models of planetary interiors (McKenzie and Weiss, 1975, 1980; Sharpe and Peltier, 1978, 1979; Schubert et al., 1979; Peltier and Jarvis, 1982). The possibility of a breakdown of these relationships in the case of time-dependent convection raises concerns about the validity of these applications. This is particularly so for inferences based upon present day (i.e., instantaneous) measurements of plate speed and heat flow. Thermal history calculations may not be affected as seriously by these considerations, since they deal with average values of heat flux and velocity on time scales which are long compared to the convective overturn time. However, all thermal history models of the Earth are pinned to one instantaneous datum the present day total heat flux from the Earth’s interior. If convection in the mantle is highly time dependent, the present day heat flow may not coincide with the appropriate mean value of heat flux over the past few convective overturns. Figure 12 indicates that, in model calculations, an instantaneous value of F,,’ may deviate from the mean by as much as

nevertheless important to consider this extra uncertainty when interpreting the predictions of thermal history models. In conclusion it is stressed that the observed present day plate velocities and global heat flow are instantaneous “snap shots” of highly variable time dependent surface dynamics. Numerical studies of time dependent convection in infinite Prandtl number fluids should help to emphasize this aspect of plate tectonics and will hopefully lead to a better understanding of the major dynamic processes.

25%. Actual deviations may be larger than this

Daly, S.F., 1980. Convection with decaying heat sources: con-

since Fig. 12 corresponds to a model with no internal heat sources the most steady case. One

stant viscosity. Geophys. J. R. Astron. Soc., 61: 519—547. Davies, G.F., 1977. Whole-mantle convection and plate tecton. ics. Geophys. J. R. Astron. Soc., 49: 459—486. Elsasser, W.M., Olson, P. and Marsh, B.D., 1979. The depth of mantle convection. J. Geophys. Res., 84: 613. Goodwin, A.M., 1984. Precambrian ring shields: growth, alignment and oscillation. Precambrian Res., in press. Green, D.H., 1975. Genesis of Archean peridotitic magmas and

—

—

mitigating factor in the case of the Earth, is that heat loss from the interior occurs due to the convective re-cycling of several surface plates, each with different characteristics. Thus the total heat loss averages the effects of plates in different dy-

Acknowledgements Financial support for this research from the University of Toronto and the Natural Sciences and Engineering Council of Canada is gratefully acknowledged. References Allegre, CJ., Hart, S.H. and Minster, J.-F., 1984. Chemical structure and evolution of the mantle and continents determined by inversion of Nd and Sr isotopic data. Part 2. Numerical experiments and discussion. Earth Planet. Sci. Lett., in press. Andrews, D.J. and Sleep, N.H., 1974. Numerical modelling of tectonic flow behind island arcs. Geophys. J. R. Astron. Soc., 38: 237—251. Bickle, MJ., 1978. Heat loss from the Earth: a constraint on Archean tectonics from the relation between geothermal gradients and the rate of plate production. Earth Planet. Sci. Lett., 40: 301—315. Bullen, K.E., 1963. An Introduction To the Theory Of Seismology. Cambridge University Press, Cambridge, 3rd edn., 381 pp. Busse, F.H., 1967. On the stability of two-dimensional convection in a layer heated from below. J. Math. Phys., 46:

326 constraints on Archean geothermal gradients and tectonics. Geology, 3: 15—18. Gutenburg, B. and Richter, C.F., 1954. Seismicity Of The Earth and Associated Phenomena. Princeton University Press, Princeton, N.J., 2nd ed., 310 pp. Hess, H.H., 1962. History of the ocean basins. In: A.E.J. Engel, H.L. James and B.F. Leonard (Editors), Petrologic Studies: a volume in honor of A.F. Buddington. Geol. Soc. Am., pp. 599—620. Houseman, G., 1983a. Large aspect ratio convection cells in the upper mantle. Geophys. J. R. Astron. Soc., 75: 309—334. Houseman, GA., 1983b. The deep structure of ocean ridges in a convecting mantle. Earth Planet. Sci. Lett., 64: 283—294. Jacobsen, S.B. and Wasserburg, GJ., 1979. The mean age of mantle and crustal reservoirs. J. Geophys. Res., 84: pp. 7411—7427. Jarvis, G.T. and Campbell, I.H., 1983. Archean komatiites and geotherms: solution to an apparent contradiction. Geophys. Res. Lett., 10: 1133—1136. Jarvis, G.T. and McKenzie, D.P., 1980. Convection in a cornpressible fluid with infinite Prandtl number. J. Fluid Mech. 96: 515—583. Jarvis, G.T. and Peltier, W.R., 1981. Effects of lithospheric rigidity on ocean floor bathymetry and heat flow. Geophys. Res. Lett., 8: 857—860. Jarvis, G.T. and Peltier, W.R., 1982. Mantle convection as a boundary layer phenomenon. Geophys. J. R. Astron. Soc., 68: 389—427. Jowett, E.C. and Jarvis, G.T., 1984. Formation of foreland rifts. Sediment. Geol., 40: 51—72. Krishnamurti, R., 1970. On the transition to turbulent convection. Part 2. The transition to time-dependent flow. J. Fluid Mech., 42: 309—320. Krishnamurti, R., 1973. Some further studies on the transition to turbulent convection. J. Fluid Mech., 60: 285—303. Landau, L.D. and Lifshitz, E.M., 1959. Fluid Mechanics. Pergamon Press, Oxford, 536 pp. Lux, R.A., Davies, G.F. and Thomas, J.H., 1979. Moving lithospheric plates and mantle convection. Geophys. J. R. Astron. Soc., 57: 209—228. McKenzie, D.P., 1969. Speculations on the consequences and causes of plate motions. Geophys. J. R. Astron. Soc., 18: 1—32. McKenzie, D.P. and O’Nions, R.K., 1983. Mantle reservoirs and ocean island basalts. Nature, 301: 229—231. McKenzie, D.P. and Richter, F.M., 1976. Convection currents in the Earth’s mantle. Sci. Am., 235: 72—89. McKenzie, D.P. and Weiss, NO., 1975. Speculations on the thermal and tectonic history of the Earth. Geophys. J. R. Astron. Soc., 42: 131—174. McKenzie, D.P. and Weiss, N.O., 1980. The thermal history of the Earth. In: D.W. Strangway (Editor), The Continental Crust and its Mineral Deposits. Geol. Assoc. Can., Spec. Pap., 20: 575—590. McKenzie, D.P., Roberts, J.M. and Weiss, N.O., 1974. Convection in the Earth’s mantle: towards a numerical simulation. J. Fluid Mech., 62: 465—538.

Moore, D.R. and Weiss, NO., 1973. Two-dimensional Rayleigh—Bénard convection. J. Fluid Mech., 58: 289—312. O’Connell, R.J., 1977. On the scale of mantle convection. Tectonophysics, 38: 119—136. Olson, P., 1981. Mantle convection with spherical effects. J. Geophys. Res., 68: 4881—4890. Olson, P. and Corcos, G.M., 1980. A boundary layer model for mantle convection with surface plates. Geophys. J. R. Astron. Soc., 62: 195—219. O’Nions, R.K., Carter, S.R., Evensen, N.M. and Hamilton, P.J., 1979. Geochemical and cosmochemical applications of Nd isotope analysis. Annu. Rev. Earth Planet. Sci., 7: 11—38. Parmentier, E.M. and Turcotte, D.L., 1978. Two-dimensional mantle flow beneath a rigid accreting lithosphere. Phys. Earth Planet. Inter., 17: 281—289. Parmentier, E.M., Turcotte, D.L. and Torrance, K.E., 1976. Studies of finite amplitude non.Newtoman thermal convection with application to convection in the Earth’s mantle. J. Geophys. Res., 81: 1839—1846. Parsons, B. and McKenzie, D.P., 1978. Mantle convection and the thermal structure of the plates. J. Geophys. Res., 83: 4485—4496. Peltier, W.R., 1972. Penetrative convection in the planetary mantle. Geophys. Fluid Dyn., 5: 47—88. Peltier, W.R., 1980. Mantle convection and viscosity. In: A. Dziewonski and E. Boschi (Editors), Physics of The Earth’s Interior, Proceedings of the Enrico Fermi International School of Physics (Course LXXVIII). North Holland, New York, pp. 362—431. Peltier, W.R. and Jarvis, G.T., 1982. Whole mantle convection and the thermal evolution of the Earth. Phys. Earth Planet. Inter., 29: 281—304. Richter, F.M., 1973. Dynamic models for sea floor spreading. Rev. Geophys. Space Phys., 11: 223—287. Richter, F.M. and McKenzie, D.P., 1978. Simple plate models of mantle convection. J. Geophys., 44: 441—471. Roberts, GO., 1979. Fast viscous Bénard convection. Geophys. Astrophys. Fluid Dyn., 12: 235—272. Roberts, P.H., 1967. Convection in horizontal layers with internal heat generation. Theory. J. Fluid Mccli., 30: 33—49. Schubert, G. and Anderson, CA., 1984. Finite element calculations of very high Rayleigh number thermal convection. Geophys. J. R. Astron. Soc., in press. Schubert, G., Cassen, P. and Young, R.E., 1979. Subsolidus convective cooling histories of terrestrial planets. Icarus, 38: 192—211. Sharpe, H.N. and Peltier, W.R., 1978. Parameterized mantle convection and the Earth’s thermal history. Geophys. Res. Lett., 5: 737—740. Sharpe, H.N. and Peltier, W.R., 1979. A thermal history model for the Earth with parameterized convection. Geophys. J. R. Astron. Soc., 59: 171—205. Skilbeck, J.N. and Md(Cnzie, D.P., 1979. An approximate method for determining the stability of two-scale flow in the mantle. Pageoph, 117: 958—987.

327 Straus, J.M., 1972. Finite amplitude doubly diffusive convection. J. Fluid Mech., 56: 353—374. Turcotte, D.L. and Oxburgh, E.R., 1967. Finite amplitude convective cells and continental drift. J. Fluid Mech., 28: 29—42. Turcotte, D.L. and Oxburgh, E.R., 1972. Mantle convection and the new global tectonics. Annu. Rev. Fluid Mech., 4: 33—68. Turcotte~D.L., Torrance, K.E. and Hsui, AT., 1973. Convection in the Earth’s mantle. Methods Comput. Phys., 13: 431—454.

Appendix Data for power law relationships TABLE 1 Nusselt number, dimensionless boundary layer thickness and size of numerical mesh, for model solutions at different values of RB/RC. (Note: R~= 779.273) _______________________________________________ RS/RC Nu Grid 10 20 50 100 200 500 1000 5000 10000 30000 100000

4.43 5.61 7.69 9.68 12.0 16.2 20.3 33.6 41.8 59.1 85.1

18.3 14.6 10.4 8.33 7.00 5.25 4.16 2.59 2.11 1.50 1.01

2~.2 23.0 16.7 13.8 11.5 8.8 7.0 4.5 3.5 2.5 2.0

24x 24 24x 24 48x 48 80x 80 48x 48 80x 80 200x200 200x200 200x200 200x200 200x200