On the driving mechanism of plate tectonics

Tectonophysics, 38 (1977) 61-88 0 Elsevier Scientific Publishing Company, Amsterdam - Printed in The Netherlands

ON THE DRIVING

MECHANISM

FRANK M. RICHTER

*

Department

of Geodesy and Geophysics,

61

OF PLATE TECTONICS

Cambridge

University, Cambridge (England)

(Received November 5, 1975; revised version received March 25, 1976)

ABSTRACT Richter, F.M., 1977. On the driving mechanism of plate tectonics. In: J. Bonnin and R.S.. Dietz (editors), Present State of Plate Tectonics. Tectonophysics, 38 (l-2): 61-88. A series of idealized models for determining the balance of forces in plate tectonics are presented and the question of the driving mechanism is reduced to evaluating several integrals that represent the forces on the bounding surfaces of the lithospheric plates and the total buoyancy of the downgoing slab. The value of the integrals are plotted as a function of depth along the downgoing slab to graphically represent the complete balance of forces, thus showing the relative importance of the different terms and also giving a measure of the internal stresses within the downgoing slab. The models that best explain the observed velocity of plates and the distribution of stresses within the downgoing slabs as inferred from focal mechanisms require the viscosity in the upper mantle to increase with depth by a factor of about 20. The primary balance of forces in these models is between the negative buoyancy of the downgoing slab and the nonhydrostatic pressure acting on its base. The significant role of the nonhydrostatic pressure is the principal new result of these calculations.

INTRODUCTION

The purpose of this paper is to comment on the driving mechanism and also those forces which resist the motion of lithosperic plates using as a basis for discussion a self-consistent mathematical model of the lithosphere-upper mantle system. Although several specific new cases are presented, in many respects this paper supplements a larger and more complete set of models studied in conjunction with Dan McKenzie, which are now being prepared for joint publication. The single most important concept in plate tectonics is the lack of internal deformation within plates compared to the relative motion between plates. It seems reasonable, perhaps even essential, to retain this concept in formu* Present address: Dept. of the Geophysical Illinois 60637. U.S.A.

Sciences, University of Chicago, Chicago,

62

lating the dynamical problem of what forces maintain the relative motion between the principal spherical caps (plates) that make up the surface of the earth. This is done by considering a two-layer idealization in which undeformable lithospheric plates are overlying a deformable upper mantle. The driving mechanism can then be studied by modelling the interaction of these plates with the upper mantle and with one another. Observations of the actual behaviour of plates are required to test and modify the model until a sufficient degree of understanding and confidence in the results is achieved. The model considered is an idealized two-layer system composed of lithospheric plates moving away from a symmetric ridge and converging at an asymmetric subduction zone where one of the plates plunges back into the upper mantle. The model is ~onceptu~ly similar to an earlier set of models (Richter, 1973), but as will be seen, a more complete analysis is now possible. Given the velocity of the plates, motions in the upper mantle result from viscous coupling and also from the mass flux (return flow) that must exist to balance the mass transport of moving plates of finite thickness. The forces exerted by the upper mantle on the plates can then be calculated given a specific geometry for the plates and a Newtonian constitutive equation governing deformations in the upper mantle. These forces for the most part resist plate motions and the question of the driving mechanism is to find those properties of the plates themselves that produce the forces balancing this resistance. The driving mechanisms considered in this category are the negative bouyancy of the subducted lithosphere (downgoing slab) and the horizontal pressure gradients in the vicinity of the ridge crest (ridge pushing). Convection of the Rayleigh-Benard type in the upper mantle is not considered among the driving mechanisms for reasons discussed in detail later. Within the context of the model, we then have a closed, well posed problem which determines the relative importance of the various terms making up the overall force balance. It is important to keep in mind that a tractable model for the interaction of the plates with the upper mantle depends on assuming the overall geometry of the model (as given by observations) and specifying a Newtonian constitutive relation. Clearly one would like to consider a more general problem in which the geometry of the flow is among the unknowns to be solved for, but for the moment this is still beyond tractable formulation. This eventual goal does not deny the usefulness of the present model, but it does qualify it. The best assessment of the usefulness of the present model derives from comparing the results to observations that can be safely assumed to reflect the interaction of the plates with the upper mantle. The best observations of this type are the stresses within the downgoing slab inferred from focal mechanisms, the mean velocity of plates in relation to their size and type, and the large-scale gravity anomalies. Isacks and Molnar (1969, 1971) have determined a large number of focal mechanisms for intermediate and deep earthquakes within the downgoing slab and conclude that, in general, downdip extensional stresses characterize the shallower events while down-dip

63

compressional stresses predominate below a depth of about 300 km. If, however, the downgoing slabs extend to depths of about 600 km, then the compressive stresses predominate almost to the surface. The cartoon shown in Fig. 1 illustrates qualitatively their principal conclusions regarding stresses within the downgoing slabs. The second type of observations we will compare to the model results involve the average absolute velocities of the plates. By averaged absolute velocity we specifically mean a plate’s velocity with respect to the hot spot frame of reference (Morgan, 1973) averaged over the plate. The need to average arises because it is only a rigid plate’s angular velocity which is constant. Using this definition, the average plate velocities exhibit a bimodal distribution such that plates being subducted over a significant portion of their boundary have velocity in the range 7-9 cm yr-’ while the remaining plates have average velocities always less than 2 cm yr-’ (Forsyth and Uyeda, 1975). These authors also quantify the well known andimportant observation that average plate velocities are essentially independent of the areal extent of plates. This last observation has often been used to suggest that plate tectonics cannot be driven by effects along plate boundaries, but as will be shown, this objection is quite easy to overcome. The last observation used to constrain the models is that the large-scale anomalies of the earth’s gravity field (see for example Kaula, 1972) show no large and systematic increase from ridges towards trenches. This is important in connection with the horizontal pressure gradient driving the return flow in the upper mantle. Too large a pressure gradient would imply an observable increase in gravity across the largest plates and thus the actual gravity field provides an upper bound on the pressure gradient. Schubert and Turcotte (1972) have commented on this point, but while they suggest that the best way to have a small pressure gradient is to have a large depth for the return flow, the implied bound is used here to constrain the magnitude and variation with depth of viscosity in the upper mantle.

... T Fig. 1. Cartoon (after Isacks and Molnar, 1969) showing qualitatively the distribution of stresses within the downgoing slab as a function of depth extent. The filled circles represent down-dip extension, the unfilled circles represent down-dip compression and the size of the circles is a qualitative guide to the activity at that depth.

64

Despite the idealized nature of the models used in the following study, a surprisingly good agreement is found between them and the observations, and while at the outset the observations used may seem limited, they are in fact strong constraints. The models suggest that the primary balance involved in plate tectonics is between the negative buoyancy of downgoing slabs and the nonhydrostatic pressure at their base. The dominant role of the downgoing slab explains the bimodal distribution of the average plate velocities and the relatively minor importance of the drag on the base of plates explains the lack of correlation between plate velocity and plate area. The importance of the nonhydrostatic pressure in resisting plate motions is the essential new feature of the present models and explains how the slabs can both drive plate tectonics and exhibit predominantly compressional stresses. THE MODEL

In this section we consider an idealized model for the interaction of lithospheric plates with the upper mantle. The model, shown in Fig. 2, consists of three regions, trench, interior, and ridge that are calculated separately and later joined to form a complete solution. The interior region consists of a domain where the flow is horizontal and, since the solution is independent of x there, this region can have any desired horizontal extent L depth scales long. The trench region must extend from the point of subduction to a distance such that the flow there is independent of x and can therefore be smoothly joined to the interior solution. Similarly, the ridge region extends from the point of symmetry (ridge crest) to where the solution is independent of 3~.A trench region three depth scales wide and a ridge region extending two depth scales from the ridge crest are sufficient to satisfy these conditions. The lithosphere, shown by ruling in Fig. 2, is 4 of the depth scale thick, a choice which corresponds to 85 km if the entire system is assumed to extend to 700 km, the depth of the deepest recorded earthquakes. The lithosphere is required to move uniformly without internal deformation except for the bending at trenches imposed by the assumed geometry. The moving lithosphere overlies an upper mantle that deforms by a Newtonian constitutive relation and which is subject to the velocity boundary conditions shown in the figure. The velocity is required to vanish at a depth D, which is assumed to be 600 km, a choice that is both very important in terms of its effect on the results and controversial. The region representing the upper mantle is also assumed to be of uniform density, a simplification which allows the neglect of the temperature equation. The model is really a two-layer idealization both in terms of the mechanism of deformation and in terms of density; the essence of the idealization being that continuous but presumably sharp gradients in the earth are replaced by step functions in the model. The usefulness of two-layer systems is well established in geophysics; for example, the heat flow models for ocean ridges (McKenzie, 1967; &later and Franche-

65 TRENCH

REGION

RIDGE -pm-

W.0

u-0

REGION roD

-----c(

WIO.U.0

Fig. 2. An idealized two-layer model for plate tectonics. The top diagram shows the lithosphere (ruled) overlying the upper mantle which is assumed to deform via a Newtonian constitutive relation. The model is made up of three regions, which are joined to form a complete model whose total horizontal extent depends on the choice of L, the extent of the interior region. The boundary conditions are given in terms of the horizontal velocity u and the vertical velocity w, and the subscripts x or z indicate the partial derivative with respect to that variable (i.e. w, E a~/&). The middle diagrams show the subregions in which numerical calculations are carried out given the boundary conditions now shown in terms of stream function $ and vorticity q. Numbers as subscripts indicate the subregion to which the variable belongs. The value of the stream function on the top boundary (F1 or F2) is that required for the upper mantle to balance the mass transported by the moving plates (F1 = VI D/7, Fz = UzD/7). The functions fL, go and fR, gR are the interior solution to the left and right of the downgoing slab as discussed in the text. The lower diagrams are a detail of subregion 3 showing the overlapping boundary conditions and a detail of the outflow region under the ridge crest.

teau, 1970) are essentially two-layer models. These models do of course allow temperature and therefore density variations within the lithosphere, but, as will be seen later, this can be incorporated into the present model. The lithospheric thickness used here is close to that which best explains the heat flow and elevation of ridges relative to older ocean basins (Parsons and Sclater, 1975). To simply say that other two-layer systems are useful does not necessarily imply that the present one is justified and it is therefore worthwhile to consider further the nature of the simplifications. The use of two layers to model the deformation mechanisms can be compared to a

66

more general treatment in the relatively simple interior region where Froidevaux and Schubert (1975) have discussed the effect of a continuous, nonNewtonian rheology. The results suggest that the model used here, especially when the viscosity in the region representing the upper mantle varies with depth, is adequate in the interior region. The effect on the ridge and trench calculation is uncertain but one expects that significant differences will be at most local. The effect of ignoring density variations within the upper mantle is easier to assess. Basically we are ignoring the local cooling (and corresponding density variations) of the upper mantle caused by the subduction of colder lithospheric material and are also not allowing for any other scale of convection besides that involving the plates themselves. The horizontal density gradients outside the downgoing slab are a relatively small source of motion in the upper mantle when compared to the effect of the downgoing slab itself (Richter, 1973) and therefore it seems justified to assume that all the negative buoyancy is contained within the downgoing slab. Ignoring convective scales smaller than the horizontal dimensions of the plates is based on recent experiments on convection under a moving boundary (Richter and Parsons, 1975) which are intended to be models of the interior region and suggest that small-scale convection, that is convection having horizontal scale of the order of the layer depth, cannot be a significant driving mechanism. The experiments show that in the case of RayleighBenard convection under a moving boundary all convective modes which apply stresses to the boundary so as to help or hinder its motion are naturally suppressed. The final form of convection observed was always longitudinal rolls (two dimensional cells with axes aligned in the direction of boundary motion), which only apply stresses at right angles to the direction in which the upper boundary moves. Such modes can be safely ignored if one’s interest is the maintenance of plate motions. Given the above simplifications, the equations governing the system follow. All variables are first nondimensionalized. distance : velocity : [;I density:

= u*[$]

p = p&z) + p’ 1.14

pressure:

p = gD _/ 2’

PO&

po(~)dr + D

p’

where D is the depth extent of the upper mantle, U, is the plate velocity as shown in Fig. 2, pa(z) is the density in the interior region which is used as a reference state, g is the acceleration of gravity and p. is the viscosity of the upper mantle at the base of the lithosphere. The upper limit of integration in the nondimensionalization of pressure is greater than one to account for the

61

thickness of the lithosphere. p’ will be called the perturbation pressure to avoid confusion. The balance of forces in a two-dimensional model of the lithosphere is stated by the momentum equation in integral form. D’SFip’dS

+ c~OU,Jcijnj

(1)

dl’ = 0

where Fi is the body force, the first integral being taken over the crosssectional area of the lithosphere and:

o:j = -p’aij

+ h(,)

is the stress tensor, to be integrated over the boundary, and h(,) depends on the viscosity at depth such that /+) = p. h(,). We are only interested in the force balance in the direction of plate motion, in which case the first integral is simply the total buoyancy of subducted material relative to p,(z). The second integral when expanded consists of: 1.14

P’k’,

s 1

s

Gdl,

W$$

s h(z) ‘2 s P’h’,

*,

dz’,

the integral being taken at the ridge crest

(2)

the integrated length 1

(3)

mean stress on the trench thrust fault of

the drag on the base of the plate integrated crest to trench

from ridge

the drag on the two sides of the downgoing

slab

the perturbed pressure integrated downgoing slab

(4)

(5)

along the base of the (6)

The question of the driving mechanism is thus reduced to evaluating the five integrals above along with the buoyancy integral. It is easiest to evaluate the pressure integral at the ridge crest (2) in terms of dimensional variables using a vertical coordinate t which is zero at the base of the lithosphere in the interior region and d at the top. The calculation of “ridge pushing” has already been carried out by many different authors (see McKenzie, 1972) but is given here for completeness. The dimensional perturbation pressure 5 at the ridge crest is: Etn = Prg, --. PI(#)

(7)

where pr and pI are the total pressure under the ridge crest and in the interior

68

region, respectively. Assuming that pot*) is to a sufficiently good approximation constant, pi becomes: PI = pog(d

-

(8)

8

pr differs from pI due to horizontal density differences and the elevation of

the ridge crest relative to old ocean floor (in model terms, relative to the interior). If the extra height of ridges is E, then: (9)

Pr = P,k?(d + E”-’ El

where pr is the density of material under the ridge assumed to be uniform and E”is the reduced excess elevation E(P, -- 1)/p,, which takes account of the extra elevation being under water. The lack of a large gravity anomaly associated with ridges requires that there be no net mass anomaly between ridge and interior, a condition satisfied if: Pr =

(d

t

f)

(10)

PO

Substituting eq. 10 into eq. 9, then eqs. 8 and 9 into eq. 7, the integral required (2) becomes: d+c

pou*J-p’dz’ =J 0

jklz =

(Pr - 1w 2

= 3 * 10” dyn cm-’

(11)

where terms smaller than the above by a factor e/d have been neglected. pr = 3.3 g cme3, E = 3 km, d = 85 km are used to determine the total force, which is in units of dyn cm-’ rather than dynes because of the twodimensional model. The second integral (3) that has to be evaluated involves the mean stress I? on the thrust fault that constitutes the contact between converging plates. The existence and mechanism of numerous earthquakes along this zone attests to the resistance to plate motion. However, the study of these earthquakes provides not the average stress 0 but the apparent average stress q(T where 77is the poorly known seismic efficiency (the ratio of energy radiated seismically to total energy released). The stress drop Au during an earthquake can also be estimated but Au can be much smaller than the average stress if the final stress remains large. In fact a recent study of a large number of earthquakes around the Mediterranean and in the Middle East (North, 1973) suggests that Au & 0 for all but the largest events. Accordingly the estimate of C’used here will be slightly larger than the typical stress drop associated with the largest shallow earthquakes that have been studied. The best such estimate of Ao is from the 1964 Alaskan earthquake (Au = 28 bar, Kanamori, 1970b) but many studies of other earthquakes (Au = 23 bar, Kanamori, 1970a; Au = 32 bar, Kanamori, 1971; Au = 30 bar, Wu and Kanamori, 1973; Au = 51 bar, Fukao, 1973) show only a small variation from this value. Therefore we will use Z = 100 bar in evaluating integral (3),

69

a value close to that suggested by the apparent average stress of shallow South American earthquakes (Wyss, 1970). The integral is taken over a fault plane dipping at 20” and extending to a depth of 85 km. This choice leads to a larger fault area than that used by Davies and Brune (1971) in estimating the total displacement at convergent plate boundaries by summing seismic moments. Using a larger area seems justified in that we want to include the resistance by aseismic deformation processes likely to be important deeper in the lithosphere. Using?? = 100 bar and dl = 250 km: 6dZ = 2.5 . 10” dyn cm-’ (12) s Given the large uncertainty of this estimate and its similar magnitude to the perturbation pressure integral at the ridge, we assume they cancel each other. This cancellation is independently justified by the estimate of Forsyth and Uyeda (1975) that the ridge pushing term is balanced by the resistance due to plate-plate interactions at trenches and along transform faults. Transform faults cannot be included explicitly in a two-dimensional representation and their effect, most likely to be small, is neglected. In order to evaluate the integrals representing the total viscous drag on the moving plates and the perturbation pressure acting on the bottom of the downgoing slab, the equations governing the model upper mantle must be solved. Assuming a Newtonian viscosity which may depend on depth, these equations are: (13)

v”$ =q

(14)

P = /M(z)

(15)

where 71is the vorticity, 71= (au/az) - (C)w/ax), $J is the stream function such that u = a G/&z and w = --a $/ax, and p is the viscosity which as before varies with depth according to h(z). Given h(z), eqs. 13 and 14 can be solved numerically using the methods described in an earlier paper (Richter, 1973). The trench region has been broken up into three overlapping subregions to allow the computations to be carried out in rectangular domains. Subregion 3 can have different vertical extent depending on how deep we want the downgoing slab to reach. The ridge region includes a portion of boundary where an outflow w = 2Fz/3 is specified. This outflow represents the upwelling of material under the ridge crest having a net mass flux equal to that required by the moving plate and was chosen a linear function of /3 because this is the simplest form that does not cause an artificial singularity in the vorticity at 0 = 0. The only other boundary conditions that require comment are those imposed on the sides of subregions 1, 2, and 4. The condition that the stream function and vorticity are zero under the ridge crest results from assuming the entire

70

model to be symmetric about the ridge axis. The remaining boundary conditions on the sides involve the functions f and g which represent the interior solution and therefore satisfy:

$-%9=0

(16)

CLg

(17)

a2

which for simple choices of p can be integrated directly. The numerical calculations and the results given later are in terms of nondimensionalized variables, the nondimensionalization for stream function and vorticity follows from that previously given for velocity and distance, and P’ = El/PO= h(z) Hereafter all variables are presumed nondimensional unless otherwise stated. Given $ and n, the perturbation pressure at any point can be found by integrating the expressions for the pressure gradient.

(18)

CL az

_f+)

i%

ax

_._

2aiz_ifk a2 axaz

(19)

We are then in a position to evaluate the integrals representing the forces on the plates arising from the perturbation pressure and viscous drag. The viscous drag integrals (4) and (5) can be evaluated using the vorticity along the boundaries.

(20) and:

s

h(z) g

dz =

-jh(z)@z

(21)

These relations hold only on the boundaries and for the boundary conditions used. Two examples of the stream function and vorticity are shown in Figs. 3 and 4. Figure 3 is a uniform viscosity case, h(z) = 1, in which the downgoing slab extends three quarters of the depth, the ridge outflow is over a length equal to half the depth scale (d = 0.5D) and only the plate to the right of the subduction zone is moving (VI = 0). The stream function maps particle trajectories and is very similar to an earlier solution (Richter, 1973), the major difference being that the present solution explicitly includes the

_ __~_ -----Fl Fig. 3. Stream function and vorticity in the trench and ridge region for the model in which the downgoing slab extends to z = 0.25 in Fig. 5A. The interior flow is simply the horizontal continuation of the solutions at the edges of the trench and ridge region.

region under the downgoing slab, which is required to calculate the effect of the perturbation pressure on the plates. The associated vorticity is dominated by local maximas at the corners of the downgoing slab. The concentration of vorticity of the comers where the model lithosphere bends through 90 degrees is in fact a nonintegrable singularity resulting from the unphysical boundary condition that both components of velocity are discontinuous there (Batchelor, 1967, p. 224). Since we propose to integrate the vorticity aIong the boundaries of the lithosphere, this singularity in the vorticity field will be further discussed later. The concentration of vorticity at the bottom of the downgoing slab is physically reasonable and in fact we

Fig. 4. Stream function and vorticity for the model having maximum viscosity ratio of 50 shown in Fig. 512.

I.Jb .

1.5b.

Mb.

0.77

t

- 4.7

B

l P=O

C

l 50.3

(14x3)

BPp=O -15.6 (- 24.0)

6.7

421

u-

,-CG - A60

XT-

U--U

VISCOSITY

1 10

’

J,,dr.-8.25

RATIO

l 15.17

. 15.60

so

l P= 15.17

IBbQ

Fig. 5. The value of the integrals making up the force balance for the different methods computed. The values of perturbation pressure given at three points on the edge of the ridge region and the integrated perturbation pressure in the trench region are with reference to the local zero perturbation pressure shown. In the interior region, the vorticity at the base of the plate and the horizontal pressure gradient are given. A. Uniform viscosity cases with three different depth extents of downgoing slab (to z = 0.5, 0.25, 0.125). Only the plate to the right of the downgoing slab is moving. B. Uniform viscosity case in which both plates move and the downgoing slab extends to a = 0.25. C. Two cases in which the viscosity increases with depth as shown by the viscosity ratio, defined as the viscosity at a given depth divided by the viscosity at the base of the ptate. The downgoing slab extends to z = 0.25, only the plate to the right of the downgoing slab is moving, and the values in parentheses refer to the case with maximum viscosity ratio of 50.

P&i

uro.1u

L.P=0 .P=

;;:

73

can show by comparing an especially constructed model to a laboratory experiment, that the vorticity is accurately resolved there. Figure 4 shows the solution of a problem that differs from the previous one only in that the viscosity varies with depth by a factor of 50. The function h(z) used is shown in the lower right corner of Fig. 5. At first sight it may seem surprising that the vorticity is so little affected by the viscosity increasing with depth, but this behaviour is reasonable since eq. 13 suggests that the vorticity will be significantly different only in regions where there are large viscosity gradients, the local value of the viscosity playing only a minor role. The stream function shows that the principal effect of the variable viscosity is that the plate drags less upper-mantle material along with it towards the trench. Solutions similar to those shown in Figs. 3 and 4 have been calculated for models having different lengths of downgoing slab, different outflow lengths d under ridges, different viscosity depth function h(z) and also cases where the plate to the left of the subduction zone is moving (U, f 0). These solutions are then used to evaluate the integrals 4 to 6, however we postpone the discussion of these in order to first describe an experiment which gives confidence in the computer model’s ability to accurately evaluate these integrals. AN E~ERI~ENT

We want a model problem that can be solved numerically and also be reproduced in a laboratory experiment. Comparing results then provides an estimate of the validity and accuracy of the numerical techniques. Consider an infinitely long horizontal rod of square cross-section rising or falling due to buoyancy forces through a layer of viscous fluid bounded above and below by rigid boundaries. The drag (viscous and pressure) on such a rod can be calculated using the model shown in Fig. 2 if U, and U, are set to zero and a sub-region similar to 3 is inserted between the rod and the upper boundary. For the specific case where the cross-section of the rod is 0.25 D X 0.25 D, the numerical solution estimates a total drag of 20.37 I.IU when the rod is near the mid-level of the layer. Balancing the drag and the buoyancy force, the predicted velocity is: - A/z U= 20.37~

X cross-sectional

area

(22)

A laboratory analogue is a rod of length 1 and cross-section 1 X 1 cm rising at low Reynolds number in a closed, fluid-filled tank 4 cm deep. The length 1 must be sufficiently large for the rod to be effectively infinite, which can be determined experimentally by using rods of increasing length until the observed velocity becomes independent of t. Such an experiment was carried out using glycerol at 21°C and plexiglass (perspex) rods 5,10, and 20 cm long. The lo- and 20-cm rods are found to be effectively infinite to

74

within the accuracy of the velocity measurements (-2%). The velocity of these rods, measured at the mid-level of the tank, was 0.298 cm set-’ (mean of 12 runs) which can be compared to the numerical calculation given the viscosity of glycerol at 21°C and the density cent_rast. The viscosity of glycerol was determined to be 10.57 g cmV1seccl and the density contrast between plexiglass and the glycerol used is -0.066 g cme3 Substituting these values into eq. 22 gives a numerically predicted velocity of 0.300 cm see-‘, the difference between calculation and experiment being well within the estimated experimental error. It may seem redundant to use an experiment to prove the accuracy of the numerical model in such a special case when our greatest concern is rightfully whether the model is a valid analogue of a vastly more complicated system. But of what use is the best analogue if we are not assured that it can be solved correctly? Thus the excellent agreement between the numerical model and the experiment is a necessary assurance that at least the model problems are solved correctly and with sufficient accuracy. THE BALANCE OF FORCES

The solutions shown in Figs. 3 and 4, and those for similar cases, provide the stream function and vorticity from which the perturbation pressure is calculated by inte~ating eqs. 18 and 19. There is then no difficulty in evaluating the resistance to plate motion arising from this pressure term (integral 6). In the case of the vorticity integrals 20 and 21, a decision has to be made regarding the singularities at the two points where the base of the lithosphere meets the downgoing slab. Since these singularities in the vorticity field are clearly unphysical, arising only because of discontinuities in the velocity boundary conditions, they are ignored by stopping the integrals involved a distance -0.1 D from the singularity. One can argue that this choice results in only a small error. In the earth, the radius of bending of the subducting plate appears to be sufficiently large that the vorticity below the bending region will not be much different from the values calculated further along the base of the plate. This implies that the error in ignoring the ~~t-h~d sin~l~ty is of the order of 0.1 L) compared to the total plate length. In the case of the lefthand-side singularity, where the base of the non-subducting plate meets the downgoing slab, local melting, suggested by island-arc volcanism, is likely to significantly reduce the stresses. The stress integrals 4 and 5 are therefore modified so that they now extend only to within 0.1 D of the singularities. These integrals can then be evaluated using relations 20 and 21. Fig. 5 summarizes the models calculated and shows the value of the integrals and the perturbation pressure of selected points. In the ridge region the perturbation pressure, with respect to the zero reference shown, is given at three depths along the edge. The calculated uniformity of pressure with depth there is required by eq. 19 and thus constitutes further evidence of the

75

accuracy of the numerical methods. The particular ridge model shown in Fig. 5A uses d = 0.5 D (see Fig. 2) but models using d = I), d = 0.25 D, and d = 0.125 D were also calculated and have corresponding perturbation pressure at the left edge of 10.00, 24.09, and 43.42. These differences can later be seen to be of little consequence to the total value of the integrated pressure over the base of the downgoing slab. A final point regarding the ridge region is that the zero reference used for perturbation pressure is consistent with the definitions used earlier to estimate ridge pushing. In the interior region the vorticity at the base of the moving plate and the pressure gradient are given. It is interesting to note that the vorticity (the drag) is about five times greater than it would be if there were no return flow. The total drag and pressure change for the interior region are determined by multiplying the vorticity and pressure gradient by the desired length of the region. In the trench region the value of the vorticity integrals along the different segments of the boundary are given and the pressure integral assumes a local zero reference for the perturbation pressure. All the vorticity integrals represent resistances to plate motion except for the integral along the base of the non-subducting plate, which is being driven towards the subduction zone by the flow induced by the downgoing slab. The value of the pressure integral along the base of the downgoing slab shown in Fig. 5 is best regarded as the local contribution in that it uses the perturbation pressure with respect to a local zero reference. The actual pe~urbation pressure at the reference point is found by adding the increase in pressure across the ridge and interior regions. Thus the total pressure integral is: s POX = 1 pdx + O.l4(Ap, total &al

+ Apr)

(23)

where 0.14 is the nondimensional width of the downgoing slab, and Apr and Ap, are the pressure change across the ridge and interior region, respectively. The results in the various regions shown in Fig. 5 are used to construct complete models by joining a given trench region to a desired length of interior flow and this is then joined to the appropriate ridge region. A convenient way of representing the resulting complete force balance of a subducting plate is to plot as a function of z the right- and left-hand sides of: x=R

j x=L

Ph

h-z0+ j+ 20

I~=L+

jk.zIpR + j 20

x=RC

x=R

T@Xlzz1 = j&Apd.z

(24)

20

where the equality holds only when z = 1, z. is the maximum depth reached by the downgoing slab, x = L and x = R are the left and right sides of the downgoing slab, x = RC is the ridge crest, and 1 now represents the width of the downgoing slab. Equation 24 is a restatement of the momentum equation in integral form (1) neglecting both the pressure integral at ridges (2) and the integrated mean stress along the thrust fault between plates (3),

76

as we have already assumed that these terms cancel one another. For convenience we shall call the plot of the righthand-side of eq. 24 the buoyancy curve and the plot of lefthand-side the stress curve. A balanced force system implies that these two curves result in a closed figure. Models with uniform viscosity Figure 6A is the buoyancy-stress curves for the uniform viscosity models shown in 5A assuming a total horizontal plate length of 10,000 km. The buoyancy curves are linear functions of z if a constant density contrast Ap between slab and interior is assumed. The buoyancy curve must join the end points of the corresponding stress curve and its slope is proportional to Ap. The lower part of stress curve shows the separate contributions to the pressure integral (eq. 23) and it is evident that the differences in pressure in the ridge model due to different choices of outflow width d discussed earlier will have little effect on the overall belance. A particularly useful property of buoyancystress plots is that the difference between the buoyancy curve and the stress curve at a given depth divided by the width of the slab 1 is a measure of the internal stress there. When the value of the buoyancy curve is less than the stress curve, the slab will be in down-dip compression, the converse relationship of the curves indicating down-dip extension. The buoyancy~tress plot thus shows both the complete balance of forces and the state of stress in the downgoing slab, which can be compared to the actual stresses inferred from focal mechanisms. The buoyancy-stress curves in Fig. 6A show the effect of increasing the depth extent of the downgoing slab, the principal result being that the greater the depth extent the more significant the local pressure integral becomes in resisting plate motions. The tendency of the models to predict greater downdip compression the deeper the slab extends is in agreement with the stresses inferred from focal mechanisms (Isacks and Molnar, 1971), and had we been able to compute slabs extending even closer to the lower boundary, the cross-over from tension to compression would rapidly rise to shallow depth. In these models the increasing compression with depth extent depends critically on the ~sumption that si~ific~t deformations end at 700 km. In general the cross-over from down-dip compression to down-dip extension in these models is too deep when compared to the observations, which suggest that the tensional region does not extend below 300 km (z = 0.6). Another property of the models is that the slope of the buoyancy curve, which measures the density contrast Ap (in units 1-1 U/gD2) required to maintain a plate velocity U, is only weakly dependent on the depth extent of the slab. Assuming that all slabs have roughly the same Ap, this implies that differences in the depth extent will not cause significant differences in plate velocity. On the other hand, in these uniform viscosity models the plate velocities depend strongly on the horizontal extent of the plate. Figure 6B shows the buoyancy--stress curve for plates having horizontal

77

R Z=O

a

T

I 1

‘

i

1

60

0

z=o I 0

8

1

j

STRESS,

J

60 STRESS, units

units

pU

1

1

1 120

i

I

I

120 ~IJ

Fig. 6. Buoyancy-stress curves (as defined by eq. 24) for the models in Fig. 5A showing the effect of depth extent of the slab when the total horizontal plate length is 10,000 km (A), and the effect of plate lengths of 6000 and 10,000 km when the depth extent is z = 0.25 (B). The buoyancy curves (the diagonal straight lines) represent the buoyancy integrated upwards from the base of the downgoing slab and thus measure the force driving plate motion. The stress curve represents the summing with respect to depth of the resisting terms: the pressure along the base of the slab, the drag along its sides and the drag along the base of the plates; the last term being plotted along z = 1. The curves form a closed figure implying that all the forces are in balance. The units of force are &lUwhere p is the viscosity and U is the plate velocity. In the upper figure the individual contribution to the pressure integral from the ridge, interior and trench region are shown by the letters R, I,T, between tick marks. The difference in a given model between the buoyancy curve and the stress curve at any depth is the internal stress in the downgoing slab in units of /U/l where 1 is the thickness of the slab. The stress is down-dip extension where the buoyancy curve is greater than the stress curve and down-dip compression in the opposite case.

extent of 6,000 and 10,000 km and these cases can be seen to have significantly different buoyancy curve slopes. This difference can be interpreted in two ways. If we require, somewhat artificially, that all plates regardless of extent move at the same velocity then the differences in slope mean that Ap must be different for different plates (assuming p is the same). A more reasonable alternative is to assume that Ap is independent of plate extent, in which case the buoyancy curves must have the same slope. In fig. 6B the buoyancy curves will have the same slope only if the scaling factor I.~Uis different for the two eases. Assuming again that p is the same from plate to plate, the plate velocity U must be different for different horizontal extent of the plates, a result which does nut agree with observed plate velocities (Forsyth and Uyeda, 1975). This d~f~~ulty is later overcome by allowing the viscosity to inerease with depth. Up to this point the value of viscosity p. has not been chosen. This can now be done by requiring that the total buoyancy of the deepest downg~i~~ slab in model af Fig. 6A be of the same order as the estimates from the studies of the thermal structure of the downgoing slab (McKenzie, 1969; Minear and Toksiiz, 1970; Turcotte and Schubert, 1971; Griggs, 1972). We then have: 124~~ U =~Ap&lz

= 10z6 dyn cm-’

and find: p. -

3 *

lP

g cm-’ see-’

(25) where the total buoyancy estimate used is that of McKenzie (1969) when the downgoing slab dips at 45” and U is taken as 8 cm yr-‘, the average velocity of oceanic plates. The value of viscosity found in this way is lower than that estimated using glacial rebound (Cathles, 1971) but it cannot be far off since it already suggests internal stress in the downgoin~ slab of about 500 bar when used to dimensionalize the results in Fig. 6. The value of viscosity is also needed in the estimate of the Dimensions pressure gradient required to drive the return flow in the upper mantle, If we let asterisks denude d~ension~ qu~tities:

using a layer depth D = 630 km, U = 8 cm yr-’ and apfax from the interior region. This pressure gradient requires an uncompensated (in the gravity sense) tilt of the lithosphere of:

where t is the departure of the lithosphere from a level surfaye and Ap is now the density contrast between the l~thos~he~ and sea water. This tilt

79

which is downwards toward the ridge would result in a total elevation change of about 600 m (superimposed on the elevation changes due to thermal subsidence of ridges) for a plate ~0,000 km long (the Pacific plate) and a corresponding gravity anomaly of about 60 mGal. Such a large gravity anomaly is simply not observed across the Pacific plate. This disagreement proves to be the greatest objection to models having uniform viscosity and a return flow above 700 km depth. For this reason, and also because of the unacceptable correlation of plate velocities with their horizontal extent, we will consider in the next section models with variable viscosity. First, however, we will use the results, shown in Fig. 53, of a model in which the non-subducting plate is moving towards the trench at one tenth the velocity of the subducting plate to discuss the balance of forces maintaining the relatively slow motion of non-subducting plates, typically continental plates. Since we have assumed the ridge pushing to be exactly balanced by the mean stress on the contact between plates at the subduction zones, the condition for uniform velocity of a non-subducting plate is simply that the integral of the viscous stresses along the base must vanish. Using the values shown in Fig. 5B, this condition is: 0.82 + 0.48fL

- 3.5) = 1.36

(28)

where L is now the tot& no~d~ension~ plate length, which is found by solving eq. 28 to be L - 4.5. This means that under the present ~sumptions of the model, only plates with total length less than 3000 km will have reasonable velocity. The simplest explanation for this unacceptable result is that the stress on the base of the nun-subducting plate arising from the flow induced by the downgoing slab has been seriously underestima~ by using a slab dipping at 90”. The effect of dip is easy to estimate using the similarity solution for this region discussed by McKenzie (1969). Decreasing the dip to 45” or 22” would increase the stress tending to drive the non-subducting plate towards the trench by a factor of two and four, respectively. Thus using dips in the range 22-45” will result in model velocities for plates less than 10,000 km long in good agreement with the observations. This result implies that the earlier choice of balancing ridge push by the integrated mean stress on the contact between plates is probably quite good. The principal conclusions that result from the uniform viscosity models are listed below. (I) The models explain the tendency of the stresses within the downgoing slab to range from down-dip compression in the deeper parts to down-dip extension at shallower depth as resulting from a primary force balance between the negative buoyancy of the downgoing slab and the combined effects of the nonhydrostatic pressure (acting on the base of the slab) and viscous drag on the base of the lithosphere. The large effect of viscous drag in the overall force balance resu&s in the cross-over from compression to tension occurring at too great a depth. (2) The models suggest that plate velocities should be correlated with

80

plate dimensions, a prediction at variance with the observations. (3) The models, modified to take account of the dip of actual subduction zones, explain the relatively slow velocities of non-subducting plates as resulting from the stresses associated with the local flow induced by the downgoing slab. (4) The models require the viscosity of the upper mantle to be about 3 - 102’ g cm-’ see -‘. This value leads to an estimate of the maximum stresses within the slab of about 500 bar, but also requires a pressure gradient driving the return flow too large by at least a factor of two. MODELS

WITH VARIABLE

VISCOSITY

The main drawback of the uniform viscosity models is that they predict plate velocities which would be correlated with plate area and also require too large an interior pressure gradient. These difficulties can be overcome by including variable viscosity with depth. Two models in which the viscosity increases with depth were calculated and the results are given in Fig. 5C. The corresponding buoyancy-stress curves are shown in Fig. 7A in which a uniform viscosity case is also given for reference. As one might expect the effect of increasing viscosity with depth is to enhance the role of the perturbation pressure and the drag along the sides of the downgoing slab and reduce the importance of the drag on the base of the plate. The immediate consequence is that the cross-over from compression to tension becomes ipcreasingly shallower as the viscosity ratio between the bottom and top of the layer gets large. A viscosity ratio in the range 10-50 would provide reasonable agreement with the inferred stress distribution within actual downgoing slabs. A second consequence of variable viscosity is that, even for a viscosity ratio as small as 10, the plate velocities will no longer be si~ific~tly correlated with plate dimensions (Fig. 7B). Fu~he~ore, by noting in Fig. 5C that the ratio of the integrated stress on the base of the non-subducting plate to the vorticity in the interior region is virtually unchanged by variable viscosity it can be seen that the earlier arguments regarding the velocity of non-subducting plates is unaffected by introducing variable viscosity. The final hurdle is estimating the dimensional pressure gradient driving the return flow in the interior to test if the associated gravity anomaly has been reduced to a reasonable level. As before, we first must determine the viscosity p. at the base of the lithosphere by requiring that the total buoyancy from Fig. 7A be of the order of 1016 dyn cm-‘. The values are p. = 1020 and 3.3 * 101’ g cm-’ see-’ for the viscosity ratios 10 and 50, respectively. These values of p. taken together with a plate velocity of 8 cm yr-l and the interior pressure gradient given in Fig. 5C result in an estimated gravity effect almost exactly the same as in the uniform viscosity case. As we are already using a relatively low estimate of the total buoyancy of the downgoing slab,

81

7=1

NE,

a z=o

I

0

t

t

i

I

I 120

60 NORMALIZED~STRESS

I 0

b 60

120

NORMALIZEDISTRESS

Fig. 7. A. Buoyancy-stress curves for three different maximum viscosity ratios (I?) from Fig. 5A and 5C when the plate is horizontally 10,000 km long and extends to a depth z = 0.25. The units are normalized such that in the cases N = 1, 10, and 50, the units are floU, j&U/S, and ~~~~ where /.+J is the viscosity at the base of the plate in each case. B. Case N = 10 for horizontaI plate lengths of 6,000 and 10,000 km when the downgoing slab extends to z = 0.25. Units &U/3.

there is little hope of using smaller values of p. and thus the solution to the present dilemma must be sought elsewhere. A simple model of the interior region is now posed in order to explicitly explore the effect of variable viscosity on the pressure gradient required to drive the return flow. Consider an interior region composed of two layers of depth dI and d2 (dl represents the upper layer, dl + d2 = 1 and z = 0 is the interface between layers) and having viscosity pl and p2, respectively. To

82

retain the analogy with the interior region used with the numerical calculations, we require the upper boundary of the present system to have velocity U relative to a rigid bottom and the system as a whole must have a return mass flux of 0.14 U. If we call the velocity in the upper layer u1 and in the lower layer u2, and let N be the viscosity ratio pz/llIr then the system has solutions:

and:

Wd2+ 2) Up= (Nd, + d,) + 2pl(%I;

d,)

[-dldz + ($-d:)z

+ (d’+$)z2]

(30)

where the pressure gradient ap,Bx is the same in both layers. The pressure gradient required to drive a mass flux of 0.14U is found from the mass flux condition: dl d2 u,dz = 0.14U u,dz + s f

(311

d2

0

The nondimensional pressure gradient as a function of the viscosity ratio N is shown in Fig. 8 for two choices of dl and d2. We require the dimensional pressure gradient, ap*/ax*, to be less than 0.07 g cm-2sec-2 if the consequent gravity anomaly for a plate 10,000 km long is to be less than 30 mGal. Thus:

aP” -=-ax*

PlU ap ~2

aX

5 0.07 g cmV2 seea

132)

where gl is required by the force balance to be 102’ g cm-r see-’ when N = 10 and 3.3 - 101’ g cm-’ set-’ when N = 50. Points on the lower curve satisfy eq. 32 when U = 8 cm/yr and D = 700 km, but points on the upper curve do not. This suggests that the unacceptably large pressure gradient found earlier in the computer models can easily be brought down to an acceptable value by increasing slightly the thickness of the low-viscosity layer at the base of the plates if the viscosity ratio is greater than 20. Another important result is that once the viscosity ratio is greater than 100 the lower layer viscosity is effectively infinite as far as the return flow and the pressure gradient are concerned. This implies that the pressuregradient and associated gravity anomaly found in the computer models cannot be reduced by simply increasing the model depth D. This is because if we increase D we would then find that we require a large viscosity increase at a depth of about 700 km to account for the distribution of stresses within the downgoing slab, This large viscosity increase would eliminate any reduction of the pressure gradient because we would be in the regime where the lower layer is effectively of infinite viscosity, in other words we would be right

83

viscosity

ratio

Fig. 8. The nondimensional pressure gradient, as a function of viscosity ratio, required to drive a total mass flux at -0.14U in a system composed of two layers of different viscosity (lower viscosity above) and with an upper boundary moving with velocity U. The upper curve is the case where the upper low viscosity layer is 30% of the total depth. The lower curve is the case where both layers have the same depth extent.

back to a model having a lower boundary at about 700 km that we have been considering all along. It therefore seems that the simplest choice is to assume a viscosity ratio greater than 20 and that the low-viscosity layer at the base of the plates extends to a depth of about 350 km. This will result in only minor changes to all computed quantities other than the pressure gradient. It has recently been pointed out to me by Dan McKenzie that an alternative choice of viscosity variation, which would significantly reduce the required pressure gradient and also make plate velocities independent of their horizontal extent, is one in which the viscosity is very low in a very thin layer at the base of the lithosphere. Such a case represents a special limit of the problem considered here, and being a viable alternative, will be explored in some detail in our joint paper. Assuming we no longer have a problem with the pressure gradient by using an appropriate viscosity structure, the other principal conciusions which follow from the variable viscosity models are listed below. (1) The velocity of the non-subducting plate can be explained in the same way as in the uniform viscosity model. (2) The principal balance of forces on a subducting plate is between the buoyancy of the downgoing slab and the combined effect of nonhydrostatic pressure acting on the base of the slab and the drag on its sides. (3) The reduced role of drag on the base of plates when the viscosity ratio is greater than 10 results in plate velocities not significantly correlated with plate dimensions. (4) Viscosity ratios in the range lo-50 will result in a good qualitative agreement between model stresses in the downgoing slab and those inferred from focal mechanism.

84

Effect of variable density contrast and local plate velocity Throughout the discussion of both the uniform and variable viscosity models above, only the simplest form of density contrast between the downgoing slab and the surrounding mantle was used. We assumed a uniform contrast while the actual density contrast in the earth is likely to be more complicated due to the various mechanisms which heat the downgoing slab (conduction, dissipative heating, adiabatic compression) and the local effect of phase changes. Smith and Toksijz (1972) have studied the effect of variable density contrast in connection with the stresses in the downgoing slab and their results support the view that variations in the forces along the boundaries of the slab are more important in controlling the distribution of stresses than the variations in the body force. The relatively minor importance of the density structure can also be seen in Fig. 9 which shows the buoyancystress plot for two different density contrasts in the case where the viscosity varies with depth by a factor of 10 (Fig. 5C). The uniform density contrast-buoyancy curve is as before and the new buoyancy curve assumes a contrast decreasing linearly with depth to zero at the base of the slab. This new curve can be regarded as an extreme case of density variation since the relatively large seismicity (in number of events and magnitude) at the largest depths to which slabs appear to reach suggests that the slab is still distinct from the surrounding upper mantle. The main point is that the difference in the state of stress of the downgoing slab for

z=o

I,1

0

60

120

units pu 5 Fig. 9. Buoyancystress curves showing the effect of two different density contrasts between the downgoing slab and surrounding mantle in the variable viscosity model with maximum viscosity ratio of 10 (Fig. 5C). The dotted buoyancy curve represents a uniform contrast and the dashed curve represents a contrast decreasing linearly with depth to zero at the base of the downgoing slab (at .z = 0.25). Units /KJ/3 where /J is the viscosity at the base of the plate. ISTRESS;

85

these two extreme forms of density contrast is relatively small compared to the effect of depth extent discussed earlier. A potentially much more important effect in terms of the distribution of stresses in the downgoing slab results from the local velocity at a given point on a plate being significantly different than the plate’s average velocity. Since plates are in fact spherical caps, the velocity of any point depends on the distance to the pole of absolute rotation and thus need not be the same as the average velocity. It is most natural to use the average velocity when discussing the overall balance of forces because by averaging we are implicitly integrating over the entire plate, and thus the buoyancy-stress plot must result in a close figure. However in regions where the local plate velocity is different from the average, the local buoyancy~tress curve need not be closed and in fact will not be unless the buoyancy of the downgoing slab depends on the velocity to the same extent as the stresses do. The stresses and the pressure integral along the base of the downgoing slab are proportional to the local velocity while the buoyancy is probably much less sensitive. In the extreme case where the total buoyancy is independent of velocity (in the range 4-16 cm/yr) the effect of local velocity on the stresses is shown in Fig. 10. The effect shown in this figure is undoubtedly exaggerated by the assumption about the buoyancy but the tendency shown by the stresses should exist to some extent in real plates. The most striking feature is the large extensional stresses associated with the slower velocity. The convergence rates near New Zealand are very slow compared to the rest of the Pacific plate and this region may be an example where the above process had lead to a detached piece of slab responsible for the isolated deep earthquakes (Adams, 1963). The effect shown in Fig. 10 can also be a consequence of temporal changes in plate velocity in which case the buoyancy

0

60 f

STRESS;

units

120 ‘8

Fig. 10. Buoyancy--stress curves showing the effect of plate velocity using the model from Fig. 5C having maximum viscosity ratio of 10. Units as in Fig. 9.

86

would be relatively insensitive to the new velocity for a considerable period of time. In general, the Pacific plate can be expected to provide the best evidence of the relationship of local velocity to state of stress of the downgoing slab, a point that has already been discussed by Forsyth and Uyeda (1975). DISCUSSION

In this study we have considered a series of idealized models for plate tectonics that determine the balance of forces involved and the results were compared to those observations which can safely be assumed to reflect the driving mechanism. In the case of subducting plates the terms driving plate motions are the negative buoyancy of the downgoing slab and to a lesser extent the effect of horizontal pressure gradients in the vicinity of the ridge crest (ridge pushing). The primary role of the downgoing slab was almost a foregone conclusion if one accepts that the experiments by Richter and Parsons (1975) rule out significant driving by convective motions having horizontal scales less than that of the plates themselves, The driving terms are balanced by the nonhydrostatic pressure acting on the base of the downgoing slab, the viscous drag on the various bounding surfaces between the plates and the upper mantle, and the stress on the fault zone marking the contact between converging plates. The mean (in time and space) stress on the contact between plates was estimated using the stress drop associated with the largest earthquakes in such regions, and then assumed to cancel the ridge pushing term. The relative importance of the remaining resisting terms depends on the depth extent of the downgoing slab and the variation of viscosity with depth in the upper mantle. The models that best fit the observed velocity of the plates and the distribution of stresses in the downgoing slab inferred from focal mechanisms require the viscosity in the upper mantle to increase with depth by a factor of about 20. In such a model the principal resisting mechanism is the nonhydrostatic pressure followed in importance by the viscous drag on the sides of the downgoing slab. The importance of these terms explains the predominance of down-dip compression in the downgo~g slab even though the slabs are the primary driving mechanism. The relatively minor role of the drag on the base of the plates in variable viscosity models explains the observed lack of correlation between plate velocity and plate area. The relatively slow velocity of non-subducting plates is explained as resulting from the local flow induced under such plates by the faster moving downgoing slab. That this effect is sufficient to explain the velocity of nonsubducting plates is an indication that the earlier assumption that ridge pushing is balanced by the mean stress on the contact between converging plates is reasonable. The most significant new result is the recognition that the nonhydrostatic pressure represents a major resisting term in the balance of forces of plate

87

tectonics. That this is reasonable can be seen if one recalls that even in the simple case of a sphere slowly falling through a fluid layer at least one third of the drag results from the pressure term (Batchelor, 196’7, p. 230). As the sphere approaches a boundary, the role of pressure becomes completely dominant. In the case of the downgoing slab the pressure resistance is even further enhanced by the required return flow in the upper mantle. ACKNOWLEDGEMENTS

I am grateful to Dan P. McKenzie and the Department of Geodesy and Geophysics, Cambridge University, for their always generous hospit~ity during a year’s visit made possible by a Fellowship from the John Simon Guggenheim Foundation. I profited from many discussions with Dan McKenzie on the topics covered in this paper. I also thank Don Forsyth for a preprint of Forsyth and Uyeda (1975), which I found particularly helpful and stimulating. REFERENCES Adams, R.D., 1963. Source characteristics of some deep New Zealand earthquakes. N.Z.J. Geol. Geophys., 6: 209-220. Batchelor, G.K., 196’7. An Introduction to Fluid Dynamics. Cambridge University Press, Cambridge. Cathles, L.M., 19’71. The Viscosity of the Earth’s Mantle. Ph. D thesis, Princeton University, Princeton, N.J. Davies, G. and Brune, J.M., 1971. Regional and global fault slip rates from seismicity. Nature Phys. Sci., 229: 101-107. Forsyth, D. and Uyeda, S., 1975. On the relative importance of driving forces of plate motions. Geophys. J.R. Astron. Sot., 43: 163-200. Froidevaux, C. and Schubert, G., 1975. Plate motion and the structure of the continental asthenosphere: a realistic model of the upper mantle, J. Geophys. Res., 80: 25532564. Fukao, Y., 1973. Thrust faulting at a lithospheric plate boundary: the Portugal earthquake of 1969. Earth Planet. Sci. Lett. 18: 205-216. Griggs, D.T., 1972. The sinking lithosphere and the focal mechanisms of deep earthquakes. In: EC. Robertson (editor), Nature of the Solid Earth. McGraw Hill, New York, pp. 361-384. Isacks, B. and Molnar, P., 1969. Mantle earthquake mechanisms and the sinking of the lithosphere. Nature, 223: 1121-1124 Isacks, B. and Molnar, P., 1971. Distribution of stresses in the descending lithosphere from a global survey of focal-mechanism solutions of mantle earthquakes. Rev. Geophys. Space Phys., 9: 103-174. Kanamori, H., 1970a. Synthesis of long-period surface waves and its application to earthquake source studies - Kurile island earthquake of October 13, 1963. J. Geophys. Res., 75: 5011-5027. Kanamori, H., 1970b. The Alaskan earthquake of 1964 - radiation of long-period surface waves and source mechanism. J. Geophys. Res., 75: 5029-5040. Kanamori, H., 1971. Focal mechanism of Tukachi-Oki earthquake of May 16,1968. Tectonophysies, 12: l-13. Kaula, W.M., 1972. Global gravity and tectonics. In: E.G. Robertson (editor), Nature of the Solid Earth. McGraw Hill, New York, pp. 385-405.

McKenzie, D.P., 1967. Some remarks on heat flow and gravity anomalies. J. Geophys. Res., 72: 6261-6273. McKenzie, D.P., 1969. Speculations on the consequences and causes of plate motions. Geophys. J. R. Astron. Sot., 18: l-32. McKenzie, D.P., 1972. Plate tectonics. In: E.C. Robertson (editor), Nature of the Solid Earth. McGraw Hill, New York, pp. 323-360. Minear, J.W., and Toksiiz, M.N., 1970. Thermal regime of downgoing slab and new global tectonics. J. Geophys. Res., 75: 1397-1419. Morgan, W.J., 1973. Plate motions and deep mantle convection. In: Studies in Earth and Space Sciences. Geol. Sot. Am. Mem., 132: 7-22. North, R.G., 1973. Seismic Source Parameters. Ph. D. thesis, Cambridge University, Cambridge. Parsons, B. and Sclater, J.G., 1977. An analysis of the variation of ocean floor heat flow and bathymetry with age, J. Geophys. Res. Richter, F.M., 1973. Dynamical models for sea-floor spreading. Rev. Geophys. Space Phys., 11: 223-287. Fichter, F.M. and Parsons, B., 1975. On the interaction of two scales of convection in the mantle. J. Geophys. Res., 80: 2529-2541. Schubert, G. and Turcotte, D.L., 1972. One-dimensional model of shallow-mantle convection. J. Geophys. Res., 77: 945-951. &later, J.G. and Francheteau, J., 1970. The implications of terrestial heat-flow observations on current tectonic and geochemical models of the crust and upper mantle of the earth. Geophys. J. R. Astron. Sot., 20: 509-537. Smith, A.T. and Toksiiz, M.N., 1972. Stress distribution beneath island arcs. Geophys. J. R. Astron. Sot., 29: 289-318. Turcotte, D.L.‘and Schubert, G., 1971. Structure of the olivine-spinel phase boundary in the descending lithosphere. J. Geophys. Res., 76: 7980-7987. Wu, F.T. and Kanamori, II., 1973. Source mechanism of February 4, 1965, Rat Island earthquake. J. Geophys. Res., 78: 6082-6092. Wyss, M., 1970. Stress estimates from South America shallow and deep earthquakes. J. Geophys. Res., 75: 1529-1544.