Lithospheric dynamics

Mike Sandiford

Oct 7, 2000

Contents

I Preliminaries
1  The tectonic agenda
    1.1  A useful description?
    1.2  The tectonic agenda
    1.3  What is the lithosphere?
    1.4  About these notes
chapter 2The geophysical observables11 section 2.1Topography11 subsubsectionMorphology of the ocean floor11 subsubsectionContinental topography13 section 2.2The geoid13 section 2.3Heat flow15 section 2.4Active seismicity16 section 2.5The in-situ stress field16 chapter 3Gravity and the lithosphere17 section 3.1Isostasy17 section 3.2The flexural strength of the lithosphere19 section 3.3Gravitational potential energy 20 section 3.4Lithospheric potential energy and geoid anomalies23 chapter 4The strength of the lithosphere25 section 4.1A rheological primer25 section 4.2Background to lithospheric rheology27 section 4.3A model lithosphere28 section 4.4Uncertainties31 chapter 5The heat of the matter33 section 5.1Conduction and advection33 section 5.2The thermal energy balance34 section 5.3Thermal time constants37 section 5.4Continental geotherms37 section 5.5Natural convection38 chapter 6Isotope Geodynamics45 section 6.1Geochemical reservoirs45 section 6.2Radioactive isotopic decay47 section 6.3Isotopic fractionation during melting50 section 6.4The chondritic earth51 section 6.5The depleted mantle54 partII 1emThe oceans57 chapter 7The ocean lithosphere59 section 7.1Age, bathymetry and heatflow59 subsubsectionThe half-space model60 subsubsectionThe thermal plate model62 section 7.2Force balance on the ocean ridge63 section 7.3Formation of the oceanic crust65 section 7.4Coupling of the -spheres?67 section 7.5Oceanic basalt chemistry68 subsubsectionRare Gas Studies69 chapter 8 Subduction and arc formation73 section 8.1Buoyancy of the ocean lithosphere73 section 8.2Therml structure of subducted slabs75 section 8.3The magnitude of slab pull76 section 8.4Arc dynamics78 section 8.5Arcs and crustal growth80 subsubsectionThe sources of arc magmas81 subsubsectionThe thermal structure of arcs83 subsubsectionMelting in the Subducting Slab84 subsubsectionThe role of water in arc magmatism85 III The continents
9  Continental Deformation
    9.1  Deformation of the lithosphere subject to an end load
    9.2  Deformation within the lithosphere due to basal tractions
10  Stretching continents
    10.1  Isostatic calculations
            Rift Phase Subsidence
            Sag Phase Subsidence
            Heterogeneous stretching of the lithosphere
    10.2  Mechanical consequences of extension
            Inttroduction
            The time scale of thermal sag
             Thermal evolution at finite extensional strain rates
            Mechanical Evolution
            The role of magmas
            The role of mantle plumes
    10.3  Sedimentation in stretched basins
            The Steer's Head
IV Synthesis
10  The driving mechanisms reviewed
    10.1  Torque balance and plate velocity.
    10.2  The African plate
    10.3  Torque balance in the Indo-Australian Plate
chapter 12Precambrian geodynamics117 chapter 13Crustal growth models119 section 13.1Rare gas constraints119 chapter Ai section A.1The equations of motioni section A.2Calculation of ridge-push forceiii


Part 1
Preliminaries


Chapter 1
The tectonic agenda

The Earth is a hot, dense planet in a cold, sparse Universe. At the most basic level the processes that govern the Earth's behaviour are simple physical consequences of this reality, combined with the constraints imposed by the materials that make up the Earth. In particular, heat transfer operates to disperse this thermal energy fluctuation while self-gravitation serves to hold the mass-anomaly together. In order to understand why the Earth is the way it is, and how it evolves in time, it is necessary to understand the consequences of these two, simple, physical processes. The intricate structure of the Earth, at all scales, is a direct consequences of these physical processes, and why it is made of the materials that it is. Moreover, these processes are responsible for the great diversity and beauty of our planet that provides the basic motivation for much of the Earth Science community - a motivation that should not be forgotten.

1.1  A useful description?

Dynamicists are concerned with changes (equivalently, motions) accompanying the passage of time in dynamical systems. The system of concern to the geodynamicist is the Earth and our geological perspective focuses our interest primarily on the lithosphere (on which we reside) and, less directly, the asthenosphere. Consequently, geodynamics is primarily concerned with the basic physical processes that modulate the time evolution of the lithosphere; with one of its principal aims being the description of the behaviour of the lithosphere in the modern Earth.

In order to develop an appropriate geodynamic description we need to understand what constitutes utility. To be useful, the dynamical description must reduce some aspect (optimistically, all aspects) of the behaviour of a complex system to a few general statements. In geodynamics, in which our concern is primarily with the motion of the lithosphere and the forces that drive the motion, we may distinguish between descriptions that are concerned solely with the motion or kinematics without regard to the forces or, alternatively, mechanical descriptions concerned with the interaction between forces and motions.

The kinematic description is far less formidable than the mechanical description and has received considerably more attention as illustrated by plate tectonics. Plate tectonics attempts nothing more than a description of the motion of the lithosphere, and as such is purely kinematic. The fundamental tenets of the plate tectonic description are:

  1. The lithosphere is composed of only a few rigid plates that deform only at their boundaries. That is, plate motions can be described in terms of translation and rotation.
  2. The deforming plate boundaries are very narrow in comparison to the lateral dimensions of the individual plates.
Plate tectonics can only provide a useful description if the number of plates remains small; if large numbers of plates are needed to describe the behaviour of the lithosphere, then the description becomes sufficiently complicated to be rendered useless. Does plate tectonics provide a useful description of the behaviour of the lithosphere? Our prejudice concerning the answer to this question - only in part - provides the key to the way in which we present the subject matter. Let us explain! Since every seismic event represents a deformation of the lithosphere, a logical consequence of the plate tectonic description is that a plate boundary should be drawn through the focal point of every seismic event (or, earthquake). In the ocean basins seismicity is restricted to very narrow zones along the mid ocean ridges, subduction zones and transform zones, that are separated by large, essentially aseismic regions. The lack of seismicity over large regions within the ocean basins implies that, to a first approximation, they do behave as rigid, or elastic, plates. In contrast, seismicity in the continents is widely distributed, even in relatively inactive continents such as Australia, with very large areas of intense, distributed seismicity. The most notable regions of distributed seismicity in the modern earth are the Basin and Range Province in the western USA, the Himalaya - Tibet region in Asia and in the Aegean Sea (which is floored by thin continental crust) in Europe. In these locations, active deformation is taking place over vast areas, which are equally as large as the surrounding aseismic regions. In order to adequately account for such deformation the plate tectonic description would need to invoke hundreds of thousands of individual plates each with lateral dimensions of only a few kilometers.

Clearly, while plate tectonics may provide a useful description of the ocean basins it usefulness for many continental regions is doubtful. Indeed, an alternative description of the continents as essentially continuous ductile media (rather than elastic plates) and which therefore effectively ignores the existence of any discontinuities in the strain field such as faults and, more importantly, plate boundaries is useful for many purposes and will be the description we will emphasize here.

1.2  The tectonic agenda

At the very outset it is useful to establish our agenda in relation to other Earth science studies, such as structural geology, petrology, sedimentology, geomorphology and seismology. The overlap with these fields is substantial, but significant differences remain. We will be concerned primarily with those processes that give rise to the large-scale, regular (or so-called first-order) features that characterize the modern Earth. For example, we will be explicitly interested in understanding the controls on We give due warning that we will not generally be concerning ourselves with the small-scale, higher-order variations in topography, strain, composition, and thermal regimes that are very much the stuff of geomorphology, structural geology and petrology, respectively. For example, we will ignore such things as the detailed lithostratigraphic record of individual sedimentary basins, but we will expend considerable effort in trying to understand the controls on subsidence within basins in general. Similarly, the role of assimilation, fractionation and crystallisation that leads to the great variation of magma types that so excite igneous petrologists will be overlooked in favour of the general nature of the processes that lead to large-scale melting events and the way in which these processes naturally give rise to the characteristic chemical signature of the average melt product.

The great challenge and fun of the geological enterprise lies in deciphering the effects of the first-order processes from the higher-frequency geological noise that make the geological record so complex (and, let it be said, at the same time so rich). From the geodynamic point of view the distinction between the first-order processes and the high-frequency noise is readily appreciated once the significance of the lithosphere as the fundamental tectonic entity is acknowledged. This implies for first-order features a natural length scale of the order of 100 km (or greater), as dictated by the thickness of the lithosphere, and natural timescales of the order of 10-100 Ma (at least), as dictated by the thermal response time of the lithosphere.

Where from comes the noise? At least a large part of the rich and diverse geological noise is modulated by exogenic processes. For example, whereas landscapes are very regular when viewed at the 105 - 106 m-scale (ie., at lithospheric scales), erosion and gullying have natural length-scales that produce very fragmented landscapes on the 101-104 m-scale. Similarly, exogenic forcing, as for example associated with Milankovitch cycles, with natural time scales (102 -106 a) that are very short compared to the thermal time constant of the lithosphere (107 a), influence the third- to fifth-order sequence boundaries that are so obvious in the lithostratigraphic record at the 10-1 - 102 m-scale. It is the rheological anisotropy caused by this short-wavelength lithological variation which is responsible for the structural chaos of orogenic zones when viewed at the map scale, and is the detail which we seek quite deliberately to avoid.

1.3  What is the lithosphere?

As we have emphasized one of our primary concerns is to provide an insight into the nature and behaviour of lithosphere and, particularly, the processes that contribute to its formation. To begin with we need some common idea of what constitutes the lithosphere. We provide two working definitions (Figure 1.1):

  1. The thermal lithosphere is that part of the earth where the dominant mode of heat transfer is conduction (Figure 1.1a).
  2. The kinematic (or mechanical) lithosphere is that part of the earth where there are no substantial vertical gradients in the horizontal component of the velocity field (Figure 1.1b). Of course, large horizontal gradients in the horizontal component of the lithospheric velocity field occur because rotation of the lithosphere around vertical poles may occur without appreciable internal deformation.
Both definitions allow the distinction of the lithosphere from the underlying asthenosphere, where heat transfer is primarily by convection, and which is, as a consequence, characterized by large vertical gradients in the horizontal component of the velocity field.

Figure 1.1: The outer part of the earth can be divided into two distinct zones, the kinematic or mechanical lithosphere (heavy stipple) and the asthenosphere (light stipple) separated by a boundary layer (intermediate stipple) which shows transitional properties. In the lithosphere heat is transported by conduction and thus the geotherm, shown in Figure 1.1a, is a strong function of depth, whereas in the underlying asthenosphere heat is transferred by convection and the geotherm is an adiabat. In the lithosphere there are no significant vertical gradients in the horizontal component of the velocity field (represented by the arrows in Figure 1.1b). Formally the thermal boundary layer is included with the mechanical lithosphere as part of the thermal lithosphere. The Moho, or crust-mantle boundary, lies within the mechanical lithosphere as indicated by the dashed line.

The lithosphere comprises the crust as well as a significant part of the upper mantle, down to a depth of   100 - 150 km. The exact thickness of the lithosphere is not easily determined by conventional geophysical methods such as seismology since the lower boundary of the lithosphere, the lithosphere - asthenosphere boundary, is not defined by a prominent compositional discontinuity. Rather it is defined by changes in the rheology of the mantle-forming rock, peridotite. Since the composition of the mantle is essentially constant throughout, mantle density, rm , is dependent primarily on the temperature (through the volumetric coefficient of thermal expansion, a) and the confining pressure (through the coefficient of compressibility, b). An increase in the temperature causes a decrease in density, while an increase in pressure causes an increase in density. At a constant depth and hence at constant pressure small lateral changes in temperature result in variations in the density. The strength of all materials, including peridotite, is a function of temperature; at asthenospheric temperatures (approximately 1300 oC) peridotite ( for peridotite a is about 10-5 oK-1) is sufficiently weak that it will convect in response to stresses resulting from the density variations associated with small lateral temperature gradients whereas in the lithosphere peridotite is sufficiently strong that it does not undergo convection (although it may subduct as an essentially rigid slab). Since the rheological properties of peridotite are likely to be a continuous function of temperature in the range appropriate to the base of the lithosphere there must be a transition zone separating the lithosphere and asthenosphere. This zone is termed the lower thermal boundary layer and is formally included in the thermal lithosphere, but excluded from the mechanical lithosphere (Figure 1.1).

It is important to realise that these definitions make no direct correspondence between the lithosphere and the compositional structure of the earth. The compositional structure of the earth is constrained in a number of ways, most notably by the seismic velocity structure determined by the way in which seismic energy is radiated following large earthquakes. This approach leads to the distinction of three compositionally distinct shells within the earth bounded by prominent discontinuities in the seismic velocity structure, namely the crust, the mantle and the core. The crust is typically about 6-7 km thick beneath the oceans and 30-50 km thick in the continents with maximum thicknesses of about 70 km beneath the great mountain belts. The detailed velocity structure of the continental crust is complex reflecting a complex compositional structure and history. In comparison the velocity structure of the mantle is more regular with weak discontinuities at about 450 and 650 km, consistent with phase transitions in material with the composition of peridotite. The 650 km discontinuity, corresponding to the maximum pressure stability limit of the silicate spinel structure, is commonly taken as the boundary between the upper and lower mantle. A low seismic velocity region is sometimes present in the upper mantle at depths between 100-200 km. In this scheme, the lower boundary of the lithosphere occurs within the upper mantle probably at depths appropriate to the low velocity zone at around 100-200 km. However, it should be noted that the low velocity zone is not universally present and therefore cannot be used to define the base of the lithosphere.

1.4  About these notes

Already we have suggested that we may need two descriptions to encompass the dynamics of the outer part of the earth; one for the oceans and one for the continents. While this is undoubtedly a simplification, it does suggest a natural subdivision of the subject. In Part 2 we discuss the oceans, about which we have a deeper understanding than the continents (Part 3). Our main aim will be to develop a description of the mechanics of the lithosphere, and so we will be concerned with the nature of the forces driving the motion of the lithosphere (which are reviewed in part 4). To do this we will need to incorporate constraints on the mechanism of formation of the lithosphere which requires an understanding of the geochemical relationships between the crust and mantle as well as lithosphere and asthenosphere. Before we do so it is important to review some fundamental principles which are essential to understanding the dynamics of both the oceanic and the continental lithosphere (Part 1).

A vexing question concerning the preparation of the notes is what level of mathematics should be included. Many of the important physical processes in tectonics admit very concise mathematical definitions, and therefore mathematics provides a natural way to describe the processes. However, we will try to describe the basic physics and have only included the most important of the mathematical definitions in the main body of the text. Peripheral mathematical definitions are included either as marginal notes or in the appendices. For those of you who struggle with the language of mathematics we hope that a qualitative appreciation of the nature of the physical process will provide a new avenue to understanding the mathematical description; certainly this is how we have learned to tackle the mathematics.

Chapter 2
The geophysical observables

The basic constraints on the first order geodynamic processes that have shaped the modern Earth are provided by the observed long-wavelength variations in topography, the geoid (which contains information about the distribution of potential energy), heat flow, seismicity and the in-situ stress field. The main features of each of these geophysical observables at the global scale is summarised below.

2.1  Topography

At the global-scale the distribution of topography is strongly bimodal (Fig. 2.1) with continental topography having a mean elevation of several hundred m and ocean bathymetry having a mean depth of about -4 km.

Figure 2.1: Histograms of topographhy (and bathymetry) at the plate and global scales.

Morphology of the ocean floor

The floors of the ocean basins have a relatively simply morphology (at least when considered at long wavelengths, being dominated by narrow ridges (the mid-ocean ridges), abyssal plains and trenches. The mid-ocean ridges represent submerged volcanic mountain ranges, with an average depth beneath sea level of 2-3 km. Broad swells in the elevation of the ridges occur along their length and locally the ridges are exposed, as is the case in Iceland. The depth of the ocean floor increases systematically away from the ridges, corresponding to time dependent changes in the density structure of the upper mantle beneath the ageing oceanic lithosphere, as reflected for example in the geoid (see Section 7.1). The depth beneath the ridges above 40 Ma oceanic crust is about 2 km, while for 80 Ma old oceanic crust it is about 3 km. The average depth of the ocean waters is therefore clearly related to the age of the oceanic lithosphere, which in turn is related to the rate of generation (spreading rate) of new oceanic lithosphere at ridges (see Section 5.2). In trenches, which mark the sight of subduction of old oceanic lithosphere water depths may reach 10 km.

This regular pattern of oceanic bathymetry is disturbed by:

Continental topography

The continents have a mean elvation of several hundered meters above sea-level with a range from about 300 meters below sea-level (at the break in the continental shelf) up to about 8 km at Mt. Everest. At the 100 km - scale,the highest average elevation is only 5 km.

Continental landscapes provide one of the most dramatic of natural fractal surfaces. The characteristic feature of the fractal is the statitiscal scale-invariance, that is the general character of the surface looks the same independent of the scale of observation. For landscapes this is strictly true only over a limited range of scales (10-1m - 105 m). The invariance of landscapes on this range of scales suggests that the processes operating to sculpt the landscape (i.e., erosion, mass wastage and deposition) are either scale-invariant themselves or that when acting together produce the scale-invariant effect. Erosion acts to roughen landscape at a wide range of scales while both mass wastage and deposition act to smooth the landscape, at the short and long-wavelength, respectively.

2.2  The geoid

Because of the inviscid nature of water (when viewed at geological time-scales) the mean sea-level represents gravitational equipotential surface in the earths gravitational field and therefore its elevation can be related to the density distribution at depth (as well as the the orbital dynamics of the Earth which cause an ellipticity in mass distribution and hence in the equipotential surfaces). The observed geoid height, as determined from the mean sea-level elevation, is expressed in terms of an anomaly relative to some fixed elevation (i.e., the predicted geoid for an idealised sphercially symmetric, rotating mass distribution).

Variations in the gravitational potential energy of the lithosphere, DUl, correlate with the dipole moment of the near-surface density distribution and therefore they can be directly related to the lithospheric component of the observed geoid anomalies, DNl :

DUl = g2
2  p G
 DNl
where g is the gravitational acceleration, and G is the gravitational constant (note that the potential energy varies with the geoid anomaly as approximately 0.23 x 1012 N m-1 per meter). The main problem with using this relationship is resolving the lithospheric contribution of the geoid anomalies from the much larger amplitude anomalies associated with the dynamic processes of plate tectonics and mantle convection.

Positive geoid anomalies of up to 10 - 15 m associated with a number of mid-ocean ridge segments, as well as age-correlated geoid offsets across fracture zones imply that ageing of the ocean lithosphere is accompanied by a decline in potential energy. The geoid anomaly predicted for the cooling half-space model (as well as the thermal plate model) for young ocean lithosphere is about d  (DNo )/d  t = - 0.15 m/Ma, which compares favourably with the observed geoid anomaly over the Mid-Atlantic Ridge at 44.5oN, and elsewhere, as well as with the geoid offsets across fracture zones. The total geoid anomaly of -12.7 m over 84 Ma corresponds decline in Ul of 2.9 x 1012 N m-1.

In comparison with the mid-ocean ridges, the geoid anomalies associated with continental margins and the interior of continents are far less clear. On the basis of averages taken over large areas there appears to be no systematic difference in the geoid height between old ocean basins (older than Cretaceous) and continental masses. Such an interpretation implies that the mean potential energy of the continental lithosphere is equivalent to old ocean basins. However, the data show very substantial differences between continents, with the mean geoid of the African continent some 40 m higher than the North American continent and 10 m higher than the mean for the Atlantic and Pacific ocean basins older than Cretaceous. This observed intercontinental variation far exceeds the plausible lithospheric contributions to geoid anomalies and therefore must reflect long-wavelength sub-lithospheric contributions. Moreover, a number of continental margins are characterized by distinct positive anomalies of the order 6 m across the transition from the ocean basin to sea-level continent and imply that a continental lithospheric column supporting sea level elevation has the potential energy equivalent to ocean lithosphere of age about 44 Ma.

Since the lithospheric contribution to the geoid anomaly reflects the dipole moment of the near-surface density distribution, the observed geoid anomalies across continental margins can also be used to constrain the continental lithospheric density structure. A lithospheric thickness of 125 km and a crustal density of 2750 kg /m3 is consistent with a continental marginal geoid anomaly of + 6 m. Moreover, such a density structure is consistent with the interpretation that an isostatically compensated continental lithospheric column supporting about 1 km of surface elevation above sea level is in potential energy balance with the mid-ocean ridges. While the generally poor resolution of the geoid in mountainous regions precludes definitive correlation between topography and potential energy within the continents, some evidence of the correlation is provided by the lithospheric contribution to the geoid anomaly of 24 - 27 m for the Andean Altiplano. Such inferences are consistent with a geoid that varies with continental topography as 6 - 7 m/km, corresponding to a potential energy variation of about 1.3 x 1012 N m-1. For a continent with an average elevation of 500 m, this correlation suggests a mean continental potential energy of 0.997 x UMOR

2.3  Heat flow

Since the lithosphere represents a conductive lid to mantle convection the variation in heat flow measured at the surface provides insight into the vertical structure of the lithosphere.

In the ocean lithosphere, the highest heat flow regions are associated with mid-ocean ridges, where absolute values are very variable (50-300 mW m-2) due to intense, but localised, hydrothermal activity. In older, deeper lithosphere the heat flow measurements become less variable and gradually decline. The decline in heat flow approximates the same dependence on t0.5 as bathymetry. This age-bathymetry-heatflow law for ocean lithosphere which provides one of the most profound insights into the structure of the lithosphere (see Chapter 7).

2.4  Active seismicity

The distribution of seismicity in the lithosphere provides a direct ßnapshot" of the active deformation.

2.5  The in-situ stress field

The orientation of the stress field in the lithsophere provides an imprtant constraint on mechanisms driving plate motion (see Chapter 10).

Chapter 3
Gravity and the lithosphere

The Earth is self-gravitating in as much as it creates its own gravitational field. The gravitational body force exerted on material of unit volume in this field is given by the product of its density and the acceleration due to gravity (which is function of the mass distribution in the Earth). The way mass is distributed in the lithosphere is central to the notion of isostatic compensation and the consequences of gravity acting on the mass distribution in isostatically compensated lithosphere is central to all of geodynamics.

3.1  Isostasy

It has long been recognized that the surface elevation of the continents is in some way related to the density distribution in the subsurface. The gravity field at the surface of the earth reflects the distribution of mass at depth and gravity measurements across mountain belts show that regions of high elevation generally have a deficiency of mass at depth. That is, there is some specific depth beneath the mountain range where the rocks have a lower density compared with rocks at the same depth beneath the low lying regions flanking the mountain range. The gravity field shows that the deficiency of mass at depth is, to a first approximation, equal to the excess mass in the mountains, implying that at some great depth within the earth, termed the depth of isostatic compensation, the mass of the overlaying rock is equal and independent of the surface elevation. From the gravity data alone, it is not possible to determine exactly how the density is distributed in order to compensate the topography, and two rival isostatic models have been proposed, referred to as Pratt Isostasy and Airy Isostasy (Figure 3.1). These models were proposed long before the concept of the lithosphere was formalized and in the original formulation both models considered that isostatic compensation was achieved entirely within the crust.

Figure 3.1: Pratt (a) and Airy (b) isostasy assume rather different subsurface density distributions to compensate the excess mass associated with the additional topography in mountains. Density distribution is proportional to stipling. Seismic studies beneath mountains such as the Himalayas show that the density structure approaches that of Airy Isostasy, although it is now assumed that isostatic compensation occurs beneath the lithosphere and not at the Moho.

While modern seismic methods have shown that the structure of mountain ranges closely approximates the Airy model, it is important to dispel the notion that compensation is generally achieved at the bottom of the crust. The logical place for compensation to take place is beneath the lithosphere, and in this course we will initially assume, and then attempt to demonstrate, that this is the general case. Indeed the very existence of large lateral temperature gradients in the oceanic mantle lithosphere leads to significant horizontal density gradients. These thermally induced density changes lead to corresponding changes in the surface elevation of the oceanic lithosphere, through a kind of thermal isostasy. The horizontal force resulting from the ocean floor topography, termed ridge push, provides one of the fundamental driving forces for the motion and deformation of the lithosphere. Isostatic compensation can be achieved because the lithosphere essentially floats on a relatively inviscid substrate: the weak peridotite of the asthenosphere. Changes in the buoyancy or elevation of the lithosphere are accommodated by displacement of asthenospheric mantle. However, asthenospheric mantle is not completely inviscid (that is, its viscosity is not negligible), and its displacement in response to lithospheric loading or unloading must take a finite length of time, related to its effective viscosity. An insight into the timescales for the isostatic response of the asthenosphere to loads is provided by the rebound of continental lithosphere following the removal of glacial icecaps. The rebound following the removal of the Pleistocene Laurentide icecap shows that the time scales appropriate to this isostatic response are of the order of 104-105 years. Since mountain belts are built and decay on the time scale of 107-108 years, the isostatic response is effectively instantaneous.

3.2  The flexural strength of the lithosphere

An important question concerning isostasy is the horizontal length scale on which isostatic compensation is achieved. Gravity measurements across mountain belts suggest that, on the scale of several hundred kilometers, the lithosphere often approaches isostatic equilibrium. However the same measurements show that the mass excess associated with small-scale topography, for example, individual mountain peaks within mountain ranges, is generally not compensated. The length-scale of isostatic equilibrium (viz., regional versus local isostatic compensation) relates to the strength of the lithosphere: on small length-scales departures from isostatic equilibrium are supported by the flexural (or elastic) strength of the lithosphere. We illustrate the interaction between the flexural strength of the lithosphere and the horizontal length scale of isostatic compensation by considering a hypothetical infinitely strong lithosphere floating on a completely inviscid substrate. Subject to a substantial, localized load such as a mountain range, such a hypothetical lithosphere will be depressed such that the mass of the displaced asthenosphere is equivalent to the mass of the load. In this case the mass of the mountains is compensated over the regional horizontal dimension of the lithosphere. Of course gravity shows us that this is not the case for the earth, because the density distribution beneath mountains somehow compensates the excess mass of the mountains on a length-scale of comparable dimension to the mountains (Figure 3.2b). However, on the smaller length scale of individual mountain peaks compensation is not achieved; that is the lithosphere is sufficiently strong to distribute the load of an individual mountain over a length scale which is large compared to the lateral dimension of the individual mountain.

Figure 3.2: The response of the lithosphere (a) to applied loads depends on its flexural strength (or rigidity). An infinitely strong lithosphere compensates loads regionally over the lateral dimension of the whole lithosphere (b). At the other end of the spectrum compensation may occur completely beneath the load, so called local isostatic compensation (d). The behaviour of the lithosphere, which has some finite flexural strength, lies somewhere between these two extremes (c), implying that this flexural response of the lithosphere is an important component of deformation associaed with loading. For a visco-elastic lithosphere in which the response to loading is time dependent these figures illustrate schematically the expected temporal evolution following loading as the elastic response gives way to a viscous response through time.

As we shall see, the lithosphere has a finite flexural strength, in as much as it can support a limited amount of loading without permanent deformation. The effective flexural strength of the lithosphere is characterized by the thickness of the elastic lithosphere, which in turn is dependent on the thermal and compositional structure of the lithosphere. Importantly the flexural response of the lithosphere to an applied load may be time dependent (as would be expected for a visco-elastic material); after the emplacement of a load, the effective elastic thickness of the lithosphere may decay with time as the elastic stresses are relieved by permanent deformation.

3.3  Gravitational potential energy

Isostasy is a vertical stress balance! Local isostasy requires that at the level of isostatic compensation the mass of overlaying rock is equalized, independently of the density distribution above this level. Equivalently the normal stress in the vertical direction szz is equalized at the depth of compensation. Necessarily, the vertical density structure at, and below, the level of compensation must everywhere be equal. The vertical stress resulting from gravitational body forces acting on the mass, or body, of the rock. The vertical stress, szz, at depth h is given by:

szz = g 
h

0 
 r dz
(3.1)
where r is the density.

Under hydrostatic conditions (szz = sxx = syy) the lithostatic pressure, P, is equal to szz. Thus for constant density (depth independent) hydrostatic conditions:

P = szz = g  r h
(3.2)

Isostasy is not a complete stress balance! A complete stress balance requires that horizontal stresses as well as vertical stresses must be balanced. In order to achieve this in the presence of a gravitational field all density interfaces must be horizontal. This is clearly not the case for the lithosphere which is characterized by significant lateral variation in the density structure. For example, significant topography characterize the density interfaces at the Earth's surface (rock-air), the Moho (mantle-crust) and within the mantle lithosphere. Variations in lateral density structure within the lithosphere contribute to variations in the lithospheric gravitational potential energy, Ul. The gravitational potential energy of a lithospheric column of unit area is given by the integral of the vertical stress from the surface of the earth to the base of the lithosphere (i.e., depth of isostatic compensation):

Ul
=

h

0 
 szz dz
=
g  

h

0 
 r dz
=
g
h

0 
 (r z ) dz
(3.3)
Because the lithospheric potential energy is defined in terms of a column of unit area it has dimensions J m-2 or equivalently, N m-1. Lateral variations in gravitational potential energy, DUl, induce substantial horizontal forces, termed buoyancy forces. The difference in potential energy between two columns, 1 and 2, integrated over the depth h, is (Figure 3.3):
DUl
=

h

0 
 szz1 dz -
h

0 
szz2 dz
=
g  
h

0 
 (Dr z ) dz
(3.4)
where Dr is (r1 -r2 )z. Equation 3.3 and 3.4 simply state that in the presence of a gravitational field there is a tendency to reduce fluctuations in gravitational potential energy ( note that these equations apply equally to the forces operating on an ageing brie, as it spreads to lower its centre of gravity). While variations in potential energy generally correlate with surface elevation it is important to realise that this is not necessarily so, as illustrated in Figure 3.3. In the absence of an externally applied stress field, regions of high potential energy will experience tension, and regions of low potential energy will experience compression.

Figure 3.3: The magnitude of the buoyancy forces arising from variations in potential energy generally correlate with topographic gradients on density interfaces in the lithosphere is equal to the area between the szz vs. depth curves for both regions. Note that szz is equalized at the level of isostatic compensation, independent of the overlaying density structure.

Figure 3.4: While variations in potential energy generally correlate with variations in elevation, they need not do so as illustrated in this scenario. Clearly column 2 is gravitationally unstable. In a geological context this scenario may have some relevance to old ocean lithosphere (see Chapter 7).

The buoyancy forces arising from variations in gravitational potential energy are large (in fact they provide the fundamental driving forces for the horizontal motion and deformation in the lithosphere). These forces can be sustained because rocks have finite strength just as an immature brie does. An understanding of the rheology of rocks, that is the way in which they respond to applied forces, is fundamental to developing quantitative tectonic models.

Gravity acts on mass and as mass is distributed throughout the volume of a body the forces resulting from the action of gravity are termed body forces (other types of body forces such as magnetic forces are largely irrelevant in tectonics). As we have seen in the previous section the effect of gravitational body forces acting on two columns with different density distributions produces forces acting on the surface between the two columns. Such surface forces are dependent on the area the orientation of the surface. Since the sum of the forces acting on a body are equal to the mass times acceleration (Newton's second law) it is relatively simple to derive the equations of motion, relating displacements and forces (see Appendix A.1). The equations of motion are necessarily couched in terms of the components of the stresses (or, more correctly, the stress tensor) and therefore they must be rendered useful through combination with equations relating stresses to displacements. Such equations, termed constitutive equations, are material dependent and their study is the stuff of rheology.

3.4  Lithospheric potential energy and geoid anomalies

The contributions of lateral variations in the density structure of the lithospheric to observed geoid anomalies can be directly related to the dipole moment of the density distribution:
DN = - 2  p G
g

h

0 
 (Dr z ) dz
(3.5)
where DN is the geoid anomaly in metres. Consequently, variations in lithospheric potential energy are reflected in variations in the geoid, with an anomaly of 5 m corresponding to a variation in potential energy of about 1 x 1012 N m-1.

Chapter 4
The strength of the lithosphere

Constitutive equations specify the relations between stress and strain or the time derivatives of strain and thus can be used in conjunction with the equations of motion to relate stresses and displacements. A number of idealized classes of constitutive relations are recognized. The principal material behaviours are elastic, plastic and viscous. It should be noted that the constitutive relations appropriate to these behaviours are idealized and many materials show more complicated stress-strain relations. For example, a special class of behaviour of relevance to geology is visco-elasticity (and visco-plasticity).

4.1  A rheological primer

In simple, isotropic, elastic materials stress and strain are linearly related. An important aspect of elasticity is that all strain is recovered on relaxation of the stress. Plastic material shows an elastic response up to a critical stress, termed the yield stress, at which the material fails. For perfectly plastic behaviour the material cannot support stresses greater than the yield stress. In viscous materials stresses and strain rates are related. Thus no deformation is recovered with the relaxation of stresses. The general form of the constitutive relation for viscous materials is :

t = k   .
g
 
1/n
 
(4.1)
where k and n are constants, and [(g)\dot] is the shear strain rate. For n = 1 the relationship between stress and strain rate is linear and the material is said to be Newtonian with the viscosity equal to k. Materials with n > 1 are termed power-law or shear thinning fluids. For high values of n the behaviour of the material may approximate plasticity in as much as a dramatic change in the rate of deformation occurs with the relatively small increases in shear stress.

Figure 4.1: Elastic, viscous and visco-elastic response to a step shear strain applied at time t=0.

Visco-elastic materials show time dependent behaviour which is like elastic materials at short time-scales and like viscous materials at longer time-scales as illustrated in the response to the application of a step-like shear strain in Figure 4.1 and in a creep test in Figure 4.1.

Figure 4.2: Elastic, viscous and visco-elastic response to a creep test.

4.2  Background to lithospheric rheology

The distribution of seismicity in the ocean basins suggests that the ocean lithosphere deforms primarily by the rigid body translation and rotation of large "plates"; as we shall see these motions are sustained by the push of the ocean ridges and the pull of the subducting slabs. The lack of internal deformation of the ocean lithosphere reflects its strength; it must be strong enough to sustain the stresses arising from these fundamental driving forces (that is, the ocean lithosphere acts as a stress guide). The behaviour of the ocean lithosphere contrasts with many parts of the continents wherein intense, diffuse, internal deformation is manifested by the wide distribution of seismic events. The existence of diffuse deformation within the continental interiors implies that the continental lithosphere is, in general, considerably weaker than the oceanic lithosphere (or that the magnitude of the stresses sustaining deformation of the continents is larger than for oceans). In addition, differences in the behaviour of the oceanic and continental lithosphere are governed by differences in their buoyancy. Old oceanic lithosphere is negatively buoyant allowing eventual subduction and the effective recycling of the oceanic lithosphere in the convective mantle; one consequence is that there is no oceanic lithosphere older than about 200 Ma on the surface of the modern Earth. In contrast the continental lithosphere is always positively buoyant and cannot be subducted to any significant extent; consequently we have fragments of the continental lithosphere dating back to about 3.9 Ga.

In later chapters we consider the deformation of the continental lithosphere in both compression (which results in the development of mountain belts) and in tension (which results in the development of rifts and the attendant sedimentary basins and flanking uplifted margins). We will not be so much concerned with the internal geometry of deformation (which is very much the realm of structural geology) but the mechanical and thermal consequences of the deformation on the lithospheric scale. Our aim is to develop an understanding of the controls on some simple and regular, yet, as we shall see, fundamentally important features of continental orogenic belts, for example, the first-order controls on elevation. In order to understand the behaviour of the continents during deformation we begin with some simple notions concerning the strength or rheology of the continental lithosphere.

Figure 4.3: The stress needed to deform the lithosphere and the failure mode depend upon the depth, the strain rate, the thermal structure as well as the mineralogical composition of the deforming rocks. At shallow depths (low confining pressures and low temperatures) failure occurs in the brittle mode, at deeper levels (high confining pressures and high temperatures) deformation occurs in a ductile fashion. The strength of rocks is strongly dependent on the temperature and the strain rate. Curve 1 is appropriate to low strain rates or high geotherms, while curve 2 is appropriate to lower geotherms or higher strain rates. The brittle-ductile transition can be viewed as the depth at which failure mode switches, and coincides approximately with the base of the seismogenic zone.

4.3  A model lithosphere

The above discussion implies that the lithosphere has finite strength! That is rocks are able to sustain a finite deviatoric stress (or stress difference) by elastic deformation, all of which is recoverable on the relaxation of the stress. At stress in excess of the critical deviatoric stress rocks will undergo a permanent deformation by processes such as brittle fracture, dislocation creep, pressure solution, the actual deformation mechanism depending on the temperature, strain rate and composition of the rock. The finite elastic strength of the lithosphere implies the general rheological response is plastic, and since the response of the lithosphere can be shown in many cases to be time dependent the lithosphere exhibits a complex form of visco-plasticity. The complete description of the rheological behaviour of such a lithosphere represents a formidable challenge. However it is possible to achieve the general form of a visco-plastic response by specifying a lithospheric rheology governed by a combination of two simple failure mechanisms, as described below.

Figure 4.4: The strength of the lithosphere is given by the area under the t vs depth curves. The strength is a function primarily of the vertical compositional structure of the lithosphere, the geotherm and the strain rate. (a) is appropriate to a low geotherm or high strain rate, i.e. strong lithosphere, where the strongest part of the lithosphere is the brittle upper mantle. (b) is appropriate to a high geotherm or low strain rate (i.e., a weaker lithosphere) where there is no brittle mantle. Note that the proportion of the total strength of the lithosphere concentrated in the crust (especially near the brittle-ductile transition) increases with increasing geotherm.

Active seismicity in the continents is, by and large, restricted to depths less than about 15 kms. Since seismic energy represents the release of elastic energy at failure along a discrete fault plane, the lower limit of seismicity can be thought of as the brittle-ductile transition. The constitutive law describing failure in the brittle mode by a frictional sliding mechanism is frequently called Byrelee's law:

t = co + m sn
(4.2)
where co is the cohesion, t is the shear stress required for failure (t = (s1 - s3 )/2), sn is the normal stress on the failure plane and m is the coefficient of friction. Since the normal stress acting across the failure plane increases with depth the stress needed to cause brittle failure also increases with depth (as a linear function of depth for the case where the density is constant with depth).

The constitutive laws describing deformation in the ductile regime have a power-law form:

(s1 - s3 )
G
=



.
g
 

A




1/n


 
 exp

(Q +P V*)
R T n


(4.3)
where G is the shear modulus (units of Pa), A is a material constant with units s-1, [(g)\dot] is the shear strain rate, Q is a material constant known as the activation energy of units J mol-1, V* is the activation volume in m3 mol-3, and n is a dimensionless material constant known as the power law exponent. The exponential term in Eqn 4.3 expresses the inverse exponential dependence of strength on temperature which obeys an Arhenius relationship. The PV* term is usually small compared to Q and thus Eqn 4.3 can be approximated by:
(s1 - s3 )
G
=



.
g
 

A




1/n


 
 exp

Q
R T n


(4.4)
Equations 4.3 and 4.4 show that the strength of rocks to creep decreases exponentially with increasing temperature and increases with the strain rate. Thus for a given composition, strength must decrease with depth.

In terms of this rheology the brittle-ductile transition can therefore be understood as the depth where the stress required for failure is equal for both ductile creep and brittle failure (Figure 4.2). The material constants for creep vary significantly with the mineralogical makeup of the rock. For example, quartz-rich rocks are much weaker than olivine-bearing rocks. Thus the strength of the lithosphere can be determined only for a specific compositional structure (as well as thermal structure). To a first approximation we can consider that the continental lithosphere comprises two layers: a quartz-rich crustal layer and an olivine-rich mantle layer. For any arbitrary geotherm the strength of such a hypothetical lithosphere depends only on the material constants appropriate to quartz and olivine (Figure 4.3). Permanent deformation of the whole lithosphere at given shear strain rate, [(g)\dot], will only occur when the force applied to the lithosphere exceeds the strength of the lithosphere, Fl, given by (Figure 4.3b):

Fl   =  
zl

0 
(s1 - s3 )  dz
(4.5)

For high strain rates, low geotherms or low applied forces much of the stress in the lithosphere will be supported at least in part elastically (Figure 4.3a). Because of the dependence on strain rate the behaviour of the lithosphere is time dependent. On short time scales the lithosphere is able to support forces in an elastic fashion; on longer timescales it will relieve the applied forces by viscous, permanent deformation. The rheology of the lithosphere may therefore be considered as visco-plastic.

Figure 4.5: The distribution of stress within the lithosphere subject to a horizontal force (in compression or tension) depends on the magnitude of the force, the strain rate and the geotherm (for a given composition). An applied force, Fd, will give instantaneously a homogeneous stress distribution, si, throughout the lithosphere such that si zl = Fd In those regions where si exceeds the critical stress difference for permanent deformation the stresses will be relaxed, with resultant amplification of the stress in the stronger parts of the lithosphere. Permanent deformation of the whole lithosphere can only occur when the whole lithosphere fails (b), that is when the stress concentration everywhere exceeds the stress difference required for permanent deformation. For lower applied forces the stresses will be supported at least in part elastically (a).

4.4  Uncertainties

The simple model for lithospheric rheology outlined in the previous section is highly uncertain for a number of reasons, and consequently must be used with some caution. Firstly, the lithosphere comprises polymineralic aggregates and the extent to which the whole of the crust can be described by quartz-failure and the mantle by olivine-failure is extremely dubious. The strength-depth curve for more realistic compositional models of the lithosphere is likely to be much more complex than shown in models shown in Figs 4.3 & 4.3. Secondly, the material constants are determined in laboratory experiments carried out at strain rates appropriate to human activity (10-5-10-8 s-1), as opposed to geological strain rates (10-13-10-16 s-1). This is particularly important because only small changes in the value of the material constants have a large effect on the calculated strength. Finally, the model excludes a number of deformation mechanisms, such as pressure solution which are known to be operative at least under some conditions.

Such uncertainties render futile the calculation of the absolute strength of the lithosphere using this simplistic model. For example, a 10% uncertainty in the activation energy for creep corresponds to an uncertainty in the strength of the lithosphere of about a factor of 2. However, much more significance can be attached to estimates of the changes in strength accompanying changes in the physical state of the lithosphere using this model, because at least in this case uncertainties in the material constants cancel. For example, a change in the thermal regime of the lithosphere corresponding to a change in the Moho temperature of 100oC causes a change in the strength of about a factor of 2 (Fig. 4.4).

Figure 4.6: Illustration of the dependence of lithospheric strength calculated using our simple rheological model on thermal state of the lithosphere, as reflected by the Moho temperature.

Chapter 5
The heat of the matter

The Earth is a heat engine, and as with all such devices it is doomed eventually to die! Indeed, it is this very heat loss that is responsible for lithospheric growth and ultimately provides the energy source driving the motions of the lithosphere. Heat loss from the Earth's interior is accompanied both by conduction and advection, which can be simply codified in terms of the requirement for the conservation of energy - and which therefore can serve to illustrate one of the great conservation principles of physics.

5.1  Conduction and advection

Thermal energy may be transported in two fundamentally different modes; namely by conduction and by advection. In the former heat diffuses through matter while in the latter it is carried by the physical transport of matter.

Conduction

Everyday experience witnesses the fact that thermal energy is transferred from hot bodies to cool bodies, with the rate of thermal energy transfer, or heat flow, dependent on the nature of the material. Fourier postulated that the heat flow, q, is proportional to the negative of the temperature gradient, because heat flows down temperature gradients, with the constant of proportionality given by a material dependent parameter known as the conductivity, k, a postulate consistent with most subsequent experimental investigation:

q = -k grad  T
(5.1)
Written in this form grad T signifies the gradient of T (often written in the equivalent form T) and as such is a vector operator. That is, grad defines the magnitude and direction of the gradient of a scalar quantity such as temperature.

Advection

The physical motion of material results in advective transport of heat. Using a Eulerian description it is easy to see that the rate of heat advected through a fixed point is dependent on the velocity of material transport and the temperature gradient. Thus, as either the tempertaure gradient or the velocity tends to zero, the rate of heat advection must also tend to zero.

5.2  The thermal energy balance

In order to understand the thermal evolution of material subject to conductive heat transfer it is necessary to derive the thermal energy balance that expresses the fundamental physical constraint provided by the requirement of conservation of energy. Ignoring the advection of heat through the movement of the medium then the thermal energy balance in a control volume can be expressed simply as:

rate of gain of thermal energy
=
rate of thermal energy flowing into the control volume
-
rate of thermal energy flowing out of control volume
+
rate of heat production in control volume
(5.2)
Since the thermal energy of a volume of unit dimensions is given by:
Cp  r  T
where Cp is the heat capacity, the first term in Eqn 5.2, the rate of gain of thermal energy, can be expressed simply as
( Cp r T )
t
which for temperature independent r and Cp is equivalent to
Cp  r  T
t
The difference between the second and third terms in Eqn 5.2 gives the negative of the divergence of the heat flow across the boundary of the control volume:
- div q
Note that the negative sign arises because divergence actually measures the rate of loss whereas in Eqn 5.2 we are really interested in the rate of gain (or convergence). Assuming that k is independent of T then the difference between the second and third terms in Eqn 5.2 can be expressed in terms of T in the following way by substituting the expression for q given in Eqn 5.1:
k div (grad  T )
= k  ·(T )
= k

2T
x2
+ 2T
y2
+ 2T
z2


The final term in Eqn 5.2 is given by:
r H
where H is the heat production per unit mass.

Thus the thermal energy balance given by Eqn 5.2 can be given in the following ways:

T
t
= k(div  ( grad T)) + H
Cp
= k(·T) + H
Cp
= k

2T
x2
+ 2T
y2
+ 2T
z2


+ H
Cp
(5.3)
where k is the thermal diffusivity:
k = k
r Cp
The various forms of the thermal energy balance in Eqn 5.3 form the diffusion equation or, alternatively, the heat equation.

The diffusion-advection equation appropriate to whole rock advection, is obtained by adding an additional term to Eqn 5.3 which expresses the rate at which heat is advected through a point. This advective contribution is given by the product of the velocity and the temperature gradient, giving:

T
t
= k2T - v . T + H
Cp
(5.4)
where the second term on the right is the advective term.

For many purposes, it is useful to recast Eqn 5.4 in a dimensionless form:

T
t
= 2T- PeT  T
(5.5)
where the dimensionless variable T and t are given by
T = T/T0
z = z/l
t = t k/ l2
where l is the appropriate length-scale (for example, the thickness of the lithosphere), and T0 is the characteristic temperature difference at the appropriate length-scale. In Eqn 5.5 PeT is the thermal Peclet number which is a dimensionless variable
PeT = v  l
k
. The Peclet number expresses the ratio of the advective to diffusive terms. For PeT > 10 the advective term dominates the thermal evolution and diffusion can be largely ignored, while for PeT < 1 diffusion dominates and the advective term can be largely ignored. For 1 < PeT < 10 both diffusion and advection contribute to the thermal evolution.

5.3  Thermal time constants

It is useful to define the characteristic diffusive response time to thermal perturbations when heat is transferred by conduction. For example consider a bunsen burner applied to the base of a conducting plate (geologically, this is analagous to the impingement of a mantle plume on the base of the lithosphere). We are interested in how long does it takes for the top surface to feel the effect of the perturbed lower boundary condition. Clearly the response time depends on the length-scale (i.e. the thickness of the plate) and (inversely) on the thermal conductivity (or, equivalently, diffusivity). In this problem we can ignore the effects if internal heat production (which does not change the response time to a change in boundary conditions). From Eqn 5.3 the term on the right, [(T)/( t)], gives the order the temperature divided by the timescale for thermal response, t. Similarly the frist term on the right of Eqn 5.3, k( div  ( grad T)), is of order kT/ l2, where l is the appropriate length scale. Equating the two quantities gives :

t l2
k
The constant of proportionality normally used is 1/p2 giving:
t = l2
p2  k
For the lithosphere typical parameter values are l = 125 km and k = 10-6 m 2 s-1 giving t = 50 Ma.

An important feature of the thermal time constant is its dependence on l2. For the lithosphere this implies that short wavelength fluctuations are damped on timescales which are very short compared to lithospheric-scale perturbations (and, incidently, explains why a 3 m deep cellar is insensitive to seasonal atmospheric temperature fluctuations).

5.4  Continental geotherms

In the steady-state, that is when [(T)/( t)] = 0, and when heat flow is only in one direction, z, Eqn 5.3 reduces to the ordinary differential equation:

0 = k  d2 T
d  z2
+ r H
(5.6)
which can be readily integrated with appropriate boundary conditions to yield a geotherm. For example, assuming a constant distribution of heat producing elements with depth, Eqn 5.6 can be integrated using the boundary conditions appropriate to the continental lithosphere: q = - qs at z = 0 (i.e., the heat flow at the surface of the lithosphere in the direction of increasing depth is the negative of the surface heat flow), and T = Ts at z = 0 (i.e., the temperature at the top of the lithosphere is the surface temperature). The resulting integration yields the expression
Tz = Ts + qs z
k
- r H z2
2 k
(5.7)
relating the temperature at any depth, z, to the conductivity, surface temperature, Ts, surface heat flow, qs, and heat production, H. The base of the lithosphere is defined by the temperature sensitive rheological properties of peridotite and can therefore be treated as an isotherm, Tl = 1280o C. The heat flow at the surface of the lithosphere, qs, therefore reflects the thickness of the lithosphere as well as the distribution of heat sources within the lithosphere.

5.5  Natural convection

Natural convection is the motion induced by thermally created density gradients in viscous materials, otherwise known as fluids. In the presence of gravitational field a compositionally homogeneous fluid heated from below will develop an unstable density stratification. Once convective motion is initiated and hotter, less dense material are juxtaposed against cooler more dense material by slight upward displacement (and vice versa in a downwelling zone) the motion will be maintained by buoyancy forces operating on the individual fluid packets (Figure 3.1). Counteracting the tendency to rise will be the viscous drag of the material and the tendency for the rising material to loose heat to its surroundings by thermal diffusion (Figure 3.2 & 3.3). Whether convective motion is initiated depends the material properties (thermal diffusivity, coefficient of thermal expansion, viscosity). The likelihood of convective motion will be enhanced by the fluid having a high coefficient of thermal expansion, subject to a large temperature gradient over a large vertical extent and diminished by a high thermal diffusivity or by high viscosity. Initiation of the convective motion requires a finite amplitude perturbation of the thermal, and hence density, structure in order to localize upwelling and or downwelling.

Figure 5.1: For the initiation of convection an imbalance of forces is required. A fluid layer subject to an inverted density gradient due, for example, to basal heating and/or surface cooling is capable of undergoing convection. Slight upward displacement of part of the basal layer will produce a horizontal density gradient which will lead to a net upward directed buoyancy force on the displaced fluid parcel, further enhancing the prospect for upward motion and eventually fully developed convective motion. Initiation of the convective motion requires a finite amplitude perturbation of the thermal and hence density structure in order to localize upwelling and or downwelling.

Figure 5.2: The displacement of a parcel of hot light fluid is resisted by viscous dissipation or drag in the surrounding fluid (strictly the rate at which momentum diffuses into the surrounding medium), a measure of which is given by the kinematic viscosity of the fluid.

The Rayleigh number, RaT, is a measure of the relative magnitude of the destabilizing buoyancy forces (Figure 7.1) and the stabilizing effects of momentum and heat diffusion (Figure 7.2 & 7.3) and is given by:

RaT = a DT L3
k n
(5.8)
Where g is the acceleration due to gravity (m s-2), a is the coefficient of thermal expansion (oK-1), dT is the temperature difference across the fluid layer (oK), L is the vertical length of the system (m), k is the thermal diffusivity (m2 s-1) and n is the kinematic viscosity of the fluid (m2 s-1). The kinematic viscosity is the ratio of the viscosity m to the density r.
n = m
r
A measure of the relative thickness of the viscous and thermal boundary layers and is given by the Prandtl number, Pr :
Pr = n
k
(5.9)
Note that both the Rayleigh and Prandtl numbers are dimensionless. When the Prandtl number is large the fluid diffuses momentum relatively faster than it does heat. In such a fluid the slow diffusion of heat from a vertical boundary into the fluid produces a thin layer of hot, low viscosity thermal boundary layer adjacent to the walls. Through the rapid dissipation of momentum, the buoyancy of this layer will drag more fluid into motion with the result that the thickness of the viscous boundary layer is much greater than the thermal boundary layer. In a system in which the temperature gradient is between horizontal surfaces, the critical Rayleigh number for the onset of convection is RaT = 103.

Figure 5.3: Displaced packets of fluids are no longer in thermal equilibrium with the surrounds and thus gain or lose heat at a rate proportional to the thermal diffusivity of the fluid. For large diffusivities the thermal equilibration takes place quickly thereby annihilating the density contrasts that drive the convective motion.

The geometrical form of convection in 3-dimensions may vary considerably with the physical properties of the fluid. However, in all cases continuity requirements dictate a general self-similar organization. At initiation buoyancy forces induce upward motion of the heated fluid and due to the continuity requirement simultaneously cooling parts of the fluid body must begin to sink. These opposite motions obviously cannot occur at the same place and the convective motion must organize itself into convection cells. Upwelling (and downwelling) zones may either be planar or columnar in form.

Two types of fluid flow in convection can be defined; laminar or turbulent. Conditions for these two modes are described by the Reynolds number, Re (again a dimensionless parameter):

Re = w d
n
(5.10)
where w is the mean velocity of the fluid (m s-1), d is the width or diameter of the source (m), and n the kinematic viscosity (m2 s-1). If not confined to a pipe then flow becomes turbulent at values of Re > 100. Confined flow requires much higher values (Re = 2000) to become turbulent. The nature of this flow is clearly important during any sort of process of mixing since turbulence leads to much more efficient exchange of heat and material between two fluids. Convective mixing in laminar flow regimes, as appropriate to the mantle, is ultimately a process of stretching and folding. During this mixing process initial compositional heterogeneities must thin with time, eventually to the extent that diffusional exchange is capable of complete homogenisation and equilibration with the host. At temperatures close to those of the liquidi of magmas, convective flow in magma chambers will be vigorous and turbulent and under these circumstances both efficient mixing and wall-rock exchange is to be expected. As a magma body cools towards its solidus, it will change to laminar flow and eventually cease to convect.

Plate tectonics is in part a consequence of thermal convection in the mantle largely driven by radiogenic heat sources as well as residual heat generated due to the release of gravitational energy during core-mantle separation early in earth history. In contrast, convection does not occur in the interior of the geologically inactive planetary bodies like the Moon and, most probably, Mars. The rate and geometry of mantle convection control the rate of heat flow from the earth's interior to its surface, the rate at which crustal extraction takes place and the rate at which material is recycled via subduction. Any fully developed geodynamic hypothesis must therefore address the nature of mantle convective geometry in the present earth. Convection requires viscous flow or deformation. Clearly liquids or gases flow at rates observable on the scale of human patience however the asthenosphere is also capable of viscous flow albeit with very slow strain rates (of the order of 10-13 - 10-15s-1). The very high viscosity of the asthenosphere (1020 - 1022 poise) leads to the possibility of Rayleigh numbers appropriate for convection only when coupled with the very large length scales of appropriate to the mantle thickness (note that the RaT increases as a function of the cube of the length scale, but decreases only to the first power of viscosity). Convection of at least the upper part of the asthenosphere is probable on the grounds of the relative homogeneity of MORB suggesting a quite well mixed source. This conclusion is also suggested by the timescale and efficiency of extraction of the lithophile elements to the crust from the depleted upper mantle.

Systems such as the Earth with substantial density stratification organise themselves into discrete layered systems. These layers may show internal convection with little material, but significant thermal energy, exchange across layer boundaries. In this case, heat transfer between a lower denser system and the base of an overlying layer will be mediated by conduction through a thermal boundary layer. This situation occurs in the Earth at the core-mantle boundary at the so-called D'' (or double prime) layer. Here, the very dense metallic core heats the base of the silicate mantle resulting in periodic, thermally-induced, gravitational instabilities (Rayleigh-Taylor instabilities). When these instabilities achieve critical positive bouyancy they initiate the ascent of a plume of hot mantle material.

Chapter 6
Isotope Geodynamics

The crust, mantle lithosphere and asthenosphere represent distinct chemical reservoirs. The evolution of these reservoirs through time can be traced by analysis of their radiogenic isotopic signatures. In this section we provide the essential concepts needed to utilize radiogenic isotopes as a constraint on the evolution of these reservoirs. Three decay schemes are used extensively in the study of geodynamics: the Rb-Sr scheme, the Nd-Sm scheme and the U-Pb scheme.

6.1  Geochemical reservoirs

Why are there distinctive geochemical reservoirs? The differentiated state of the present Earth is largely due to the formation and segregation of partial melts in the interior of the planet. Because magmas are less dense than the solid rocks in which they form, they naturally migrate towards regions of lower pressure. The key process in the evolution of distinct geochemical reservoirs formed due to the segregation of partial melts is the way in which elements are partitioned between the melt and the residual solid. A simple model for the variation in concentration of elements in magmas during partial melting based on batch melting in which the melt stays in equilibrium with the solid until the specified fusion proportion is:

Cl
Co
= 1
F+D  (1-F)
(6.1)
where F is the fraction of fusion, Co is the concentration of the element in the initial solid, Cl is the concentration of the element in the melt, D is the bulk distribution coefficient:

D = n

i = 1 
wi KDi
(6.2)
where
KDi = Xia
Xmelta
and where n is the number of mineral phases, w is the weight fraction of mineral i, and Xia is the concentration of element a in phase i. The distribution coefficient of trace elements with respect to specific minerals is a function of ionic charge and radius. Different minerals have different affinities for individual trace elements and therefore different KD values. The D values represent the effective bulk distribution coefficient for the whole mineral assemblage coexisting with the melt. Trace elements that have high D values are referred to as compatible, while those with D < 1 are referred to as incompatible. Note that what constitutes a compatible or incompatible element depends on the specific mineral assemblage. Thus, Sr is an incompatible element for mantle-melting involving olivine and pyroxene assemblages whereas it is a compatible element during crustal melting where feldspars are involved. The crust is the repository of trace elements which are incompatible with respect to the mineralogy of the Earth's mantle.

Equation 6.1 has a number of important special solutions (Figure 6.1). If D = 1 then Cl/Co = 1 and there is no fractionation between source and magma. If D >> 1 then Cl/Co < 1 and the element is preferentially retained in the unmelted residue. If D = 0 then Cl/Co tends to 1/F and the trace element is strongly partitioned into the melt phase and its enrichment is simply a function of the amount of melting. Figure 6.1 also illustrates the effect of extraction of partial melts on the relative concentration of trace elements in the unmelted solid residue. Note that highly incompatible elements are very quickly and efficiently extracted from the melting solid and the residues of such melting are highly depleted in these elements.

Figure 6.1: Solutions to Eqn 4.1 are shown as a function of the proportion of fusion the degree of enrichment or depletion relative to their source (Co), of incompatible and compatible trace elements whose bulk distribution coefficients range from zero to ten.

The solution to Eqn 6.1 for varying amounts of melting of two trace elements both of which have D < 1 but for which DE1 ( = 0.01) < DE2 ( = 0.8) where both elements have Co = 5ppm is shown in Figure 6.1, Figure 6.1 illustrates the variation of the ratio (E1/E2)liquid as a function of F. Clearly, for small to moderate amounts of melting, there is a very significant fractionation of the two trace elements; and the ultimate destination of the melts when segregated (i.e., the outer regions of the Earth) will become preferentially enriched in E1 compared to E2. By contrast, the residue of melting always has a lower ratio E1/E2 than either the original source or the derived melt.

Figure 6.2: Concentration in the melt produced by progressive partial fusion of trace elements whose bulk distribution coefficients are 0.01 and 0.8. Initial concentration in the source was 5 ppm. Also shown is the concentration in the residual solid of the element whose D = 0.01.

It is important to note that Eqn 6.1 describes equilibrium batch partial melting which assumes that the liquid remains in contact and equilibrium with the residue up to the specified degree of fusion, at which stage it is totally extracted. This is a rather simplified model of the processes likely to operate in the earth and we will discuss some more physically realistic mechanisms later.

Figure 6.3: The concentration ratio of a highly incompatible element (Cl (1)) whose D = 0.01 and a moderately incompatible element (Cl (2)) whose D = 0.8 in the melt formed by progressive fusion.

6.2  Radioactive isotopic decay

We can use the models of the type described above to infer something about the compositions of the source regions of magmatic rocks. However, as discussed in the previous section, partial melts of rocks only reflect the composition of those source rocks indirectly. By contrast, it is assumed that at the time of melting, the isotopic compositions of the liquid and it's source rocks are the same and thus the initial isotopic compositions of magmas are the same as those of their source regions at the time of melting.

Naturally occurring radioactive isotopes decay spontaneously to yield stable daughter isotopes. The rate of decay is a physical constant unique to each specific decay scheme. We refer both to the decay constant, l, and to the half-life, t1/2 = (ln  2) /l, of an isotope. The half life being the time taken for half the initial number of parent isotope nucleiides to decay to daughter products.

The equation describing radioactive decay as a function of time is:

N = N0  exp(-lt )
(6.3)
where N0 is the initial number of radioactive parent atoms and N the number left after time t. The numbers of new daughter atoms, d*, produced in time t, is given by:
d* = N0 (1-exp(-lt ) )
(6.4)
In most geological materials, there will already be some of the daughter element present at the start of the time interval of interest which must be added to by the products of radioactive decay. Thus the total radiogenic isotope concentration, d, of a rock or mineral is given by:
d = d0 + N(exp(-lt) -1)
(6.5)
and is clearly dependent on:
  1. the amount of the parent isotope present at the start,
  2. the amount of the daughter isotope present at the start, d0 and,
  3. the time interval during which this system was closed, t and at this stage the amount of the parent isotope left is N.

Isotopes in geological samples are analyzed by use a mass spectrometer as a ratio relative to a reference non-radiogenic isotope of the same element. In the Sr-system, the reference isotope is 86Sr and in the Nd-system it is 144Nd. The so-called isochron equation shows that after the passage of t years (in 109 units) the present day isotope ratio of a closed system depends on its initial ratio and the ratio which is purely a function of the concentration ratio of the parent and daughter elements. For the Rb-Sr scheme the isochron equation is



87Sr
86Sr


=

87Sr
86Sr




i 
+

87Rb
86Sr


(exp(l t) - 1)
(6.6)
(note that in this case 87Sr is the radiogenic isotope and 86Sr the reference and 87Rb is the radioactive parent however we could equally apply the same equation to the other decay schemes with appropriate decay constants).

6.3  Isotopic fractionation during melting

For elements of high atomic mass melting processes do not affect the ratios of two isotopes of the same element because the fractionation is the same for each. However, melting frequently does produce large fractionations between parent and daughter isotopes. The relationship between the isotopic ratio, a, the parent to daughter ratio, m, which is a compositionally dependent term and the time of differentiation, t, in a closed system is given by
a
d t
= l m
(6.7)
where l is the decay constant. The importance of this is considerable not only for the determination of specific rock ages, but also because it provides a connection between the chemistry of specific regions of the earth and their radiogenic isotopic compositions. As discussed in Section 4.1 the changes in concentration of trace elements in both extracted melts and in residues during fusion depends on the D value. Because parent and daughter elements will have different D values partial melting will produce melts and residues whose parent to daughter ratios are different from the source. For instance in general DRb < DSr and DSm > DNd. Therefore melts have higher Rb/Sr and lower Sm/Nd ratios than their sources. Figure 6.3a illustrates the general process of the production of a new separate reservoir, M, by melting and melt extraction of a source, S, at time, TD, leaving a slightly modified source region, R. This specific example is true for a situation in which the parent isotope is more incompatible than the daughter (as is the case for Rb - Sr). If the opposite were true (as in the case of Sm-Nd), the extracted melt trajectory would be less steep than that of the residue. Figures 6.3b and c illustrate two other petrogenetic processes which are also relevant. Of particular importance is isotopic homogenisation which will occur by diffusion during geologically meaningful timescales at temperatures above the so-called blocking temperature. It is an important observation that even at temperatures above the blocking temperature, the rate of diffusional equilibration between adjacent solids or non-convective, non-mixing liquids is only of the order of a meter in 107 years. For this reason the establishment of large separate geochemical reservoirs like the crust or asthenospheric upper mantle, requires both large scale melt extraction and convective heat and mass transfer.

Figure 6.4: Generalized model for three broad types of temporal control on rock radiogenic isotopic composition. a) differentiation by partial melting, S = source, m = melt, TD = time of melting. b) mixing of two compositionally and isotopically distinct melts at time TM to form a homogeneous hybrid. c) Solid-state isotopic homogenisation of two adjacent rocks masses or minerals, A and B, due to metamorphism.

6.4  The chondritic earth

An important process in the chemical differentiation of the earth from an initial compositionally homogeneous globe has been the generation and segregation of partial melts to form several new separated reservoirs (crust, asthenosphere etc). A comprehensive geodynamic model of the earth requires knowledge of the age, life-span and source of these separate regions.

Because the earth has undergone continuous differentiation since its formation at about 4.55 x 109 years ago it is not a simple problem to determine what the bulk composition of the whole, undifferentiated planet was. We only have easy access to the outer crustal portions of the earth and the entire crust only represents a small proportion of the earth's total mass. As we have already seen, this part of the earth has an unrepresentative enrichment of those elements whose distribution coefficients with respect to the mantle minerals are very low. For instance approximately 70% of the earth's entire supply of K, Rb and Cs all now reside in the crust.

By contrast, the Moon which is a considerably smaller planetary body, has undergone much less differentiation than has the earth. It's crust is much more like what we suppose it's whole composition to be and the oldest dates obtained from Moon rocks are much closer to 4.55 Ga. than the oldest terrestrial dates (Amitsoq Gneiss from S.W. Greenland dated at 3.9 Ga.). This illustrates that some form of crustal re-cycling on earth has destroyed most of the original crust, activity absent from the Moon since the earliest Proterozoic. Meteorites are pieces of a fragmented planetary body (probably more than one) which probably occupied an orbit between earth and Mars. Like the Moon, this appears to have had only a short history of geological activity and was then broken-up to form the asteroid belt. Meteorites therefore represent a unique opportunity to directly examine and analyze samples from the interior of a terrestrial planet like earth. They fall into three groups:

The stoney meteorites provide mineral isochrons in both the Rb-Sr and Sm-Nd systems (O'Nions et al, 1979) that indicate the age of the terrestrial planets to be 4.5 -4.55 Ga. and indicate initial 87Rb / 86Sr and 143Nd / 144Nd ratios to be 0.69898 and 0.50682, respectively.

The chondritic meteorites have not been part of a planet that has differentiated core and mantle and their bulk composition is thought to be very similar to the inter-stellar dust from which the planets accreted. The first approximation of the whole-earth's composition is therefore the composition of the chondritic class of meteorites and it is informative to present data on terrestrial samples normalized relative to this as illustrated by Figure 6.4.

Figure 6.5: The dashed line represents the normalized rare earth content of of a chondritic source and the normalized compositions of melts due to 1 and 15% partial fusion are shown. The residue in equilibrium with the 15% melt is the lower most curve.

Given that we are able to define the initial Nd and Sr isotopic composition and age of the Earth and we can use the composition of chondritic meteorites to define its Rb/Sr and Sm/Nd ratios, we are now in a position to predict the variation in the 87Rb / 86Sr and 143Nd / 144Nd ratios for this chondritic uniform reservoir (CHUR) as a function of time (McCulloch and Wasserburg, 1978). The "mantle" lines in the the two lower diagrams in Figure 6.4 are a representation of this reservoir.

It is very common that the Sr and Nd isotopic compositions of rocks are not quoted directly as ratios but are compared to the expected values of the chondritic reservoir. This is referred to as epsilon, e, notation.

eNd(0) =






143Nd
144Nd


(0)

Ichur(0)
- 1




  x   104
(6.8)
where
Ichur(0) = Ichur(T) +

147Sm
144Nd


Chur 


 (el t - 1) \
l = 6.54   x   1012
Ichur(0) = 0.512638        and       147SM
144Ndchur
= 0.1967
eNd(0) (or eSr(0)) refer to the comparison between the present day value of the chondritic reservoir and the present day isotope value of the sample (ie the measured value). If this rock is not of zero age and if the age is known, then its initial isotopic composition at its age, T, can be calculated and compared to the chondritic reservoir at that time, eNd(T) and eSr(T). In the upper diagram of Figure 6.4, the axes could equally be presented as isotope ratio values or e values. The "bulk earth" position in this diagram is the intersection point of eNd(0) = 0 and eSr(0) = 0. This figure shows the field modern day ocean island basalts (OIB) and mid ocean ridge basalt (MORB). These reflect the composition of the present day mantle and show that though some ocean island magmatism is sampling mantle whose composition is close to the predicted chondritic Earth, the modern mantle shows some very strong deviation from this model. This is particularly so for those portions of the mantle sampled by the voluminous magmatic production at the mid ocean ridges. On this same diagram from Figure 6.4, the field of continental crust is also shown. This has high 87Rb / 86Sr (or positive eSr) values and low 143Nd / 144Nd (or negative eNd) values. As the range of Sm/Nd ratios of continental rocks is relatively limited, Eqn 6.8 reduces to express 143Nd / 144Nd as a function of time and this is illustrated as the age contours in the continental crustal quadrant in Figure 6.4. Note that this is not so for 87Rb / 86Sr as there is a very large range of Rb/Sr values for common crustal rocks.

Figure 6.6: Upper: Present day Nd - Sr isotopic covariations for important terrestrial reservoirs. MORB is mid ocean ridge basalts, OIB is ocean island basalts. Lower: Sr- and Nd-time evolutionary diagrams. T model is a model age. Data for the ocean island basalts suggests that in addition to the depleted domain there still exists a mantle reservoir that has not suffered relative Rb and Nd depletion. There is some suggestion that this is a lower mantle or lower upper mantle reservoir which still has a chondrite-like chemistry. This region is also a probable source of rare gas emissions including the now extinct 129Xe suggestive of long term isolation without melt extraction. If this were true it would have important implications in constraining the extent and dimensions of mantle convective cells.

6.5  The depleted mantle

Figure 6.4 shows that for the partial fusion of a chondritic source like that assumed to represent the undifferentiated earth, the degree of total rare earth element enrichment in the melt diminishes with larger amounts of melting and that the residual source is depleted in these incompatible elements. More importantly, the melts have lower Sm/Nd ratios than the original pristine source and the residue has a higher Sm/Nd ratio.

To a simplified first approximation (with respect to the rare earth elements), Figure 6.4 illustrates the evolution of the two important terrestrial geochemical reservoirs:

The lower two diagrams in Figure 6.4 show the present day range of Nd and Sr isotopic compositions of the continental crust with respect to a mantle reservoir. These represent simplifications of the true situation, because (Figure 6.4) the extraction of melts from the mantle bring about complimentary changes to the composition of that mantle. A mantle region that has suffered partial melting and melt extraction will have higher Sm/Nd and lower Rb/Sr and U/Pb ratios than the original chondritic unmelted mantle.

Mid-ocean ridge basalts are melts of the contemporary decompressed asthenospheric mantle beneath the oceanic ridges. Even though these are melts they have rare earth patterns which are like the solid residue in Figure 6.4 and are strongly LREE-depleted. This suggests they are the products of remelting of mantle that has already lost a melt component before. In Figure 6.4, the MORB field is at high eNd and low eSr confirming that their source mantle had achieved its depleted chemistry over a long period of geological time (perhaps 2 Ga). We then view the continents as being complimentary to the depleted mantle now being sampled by mid ocean ridge magmatism. The average age of the continental crust is given by:

< T > =
T

0 
 tm (t)  dt
(6.9)
where tm is the mass of increment of crustal addition aged t.


Part 2
The oceans


Chapter 7
The ocean lithosphere

The age of the oceanic crust, as reflected in the magnetic striping of the ocean floor, increases with distance away from the mid-ocean ridges, indicating that the ridges are the site of generation of new oceanic crust. The volcanic rocks extruded at the surface of the ridges are exclusively basalt (mid-ocean ridge basalt or MORB) which, together with there sub-volcanic intrusive equivalents - gabbros and sheeted dykes, comprise the entire oceanic crust. The total thickness of the oceanic crust generated by mafic igneous activity at the ridges is typically about 5-7 km. The structure of oceanic crust and parts of the subcrustal lithosphere can be directly observed in some ancient orogenic belts where fragments of the oceanic lithosphere have been obducted to form ophiolites during collision processes, for example the Semail ophiolite in Oman.

Figure 7.1: Schematic structure of the ocean lithosphere. The ocean lithosphere consists of about 6-7 km thick crust (heavy stiple), and the mantle lithosphere (intermediate stiple) which thickens with age (and hence distance away from the ridge as indicated by arrows). Light stiple show the region of decompression melting beneath the ridge. Stream lines in the asthenosphere may be largely decoupled from the motion of the overlying lithosphere (see Section 5.5), although asthenosphere must undergo decompression immediately beneath the ridge.

7.1  Age, bathymetry and heatflow

In ocean lithosphere younger than about 80 Ma there is a remarkable correspondence between age of the ocean crust, the depth to the sea floor (bathymetry) and the heat flow through the lithosphere (Fig. 7); with bathymetric depth increasing, and surface heatflow decreasing, with the [age]. This correspondence between age, bathymetry and heatflow is due to time dependent changes in the thickness of the lithosphere. Two models for the ocean lithosphere have been proposed in order to account for this relation: the half-space model and the thermal plate model.

The half-space model

The cooling of ocean lithosphere after formation at a ridge can be treated as a thermal conduction problem (see Chapter 3) in which a non-steady state condition (the situation at the ridge) gradually decays towards a thermally equilibrated state (as the ocean lithOSphere slides away from the ridge). The analogy (Fig 7.1a) can be made with the cooling of a semi-infinite half space, which is given by:
Tz-Tm
Ts-Tm
= erfc



z

k t
 




(7.1)
where Tz is the temperature at depth z, Ts is the temperature at the surface interface of the semi-infinite half space (which in this case is the temperature of ocean water and is taken to be 0oC), Tm is the temperature of the half space in the initial condition and which is maintained at infinite distance for all time (in our case the temperature of the deep mantle, 1280oC), k is the thermal diffusivity and t is time (the error function, erf, and its compliment, erfc, arise commonly in analytical solutions to the heat equation and related differential equations which employ similarity variables).

The behaviour of the error function, and hence Eqn 7.1, is illustrated in Figure 7.1b. As z tends to or t to 0 then:

erfc



z

k t
 




0
and Tz approaches Tm. As z tends to 0 or t to then
erfc



z

k t
 




1
and Tz approaches Ts, providing




z

k t
 




< 2

Figure 7.2: Schematic thermal structure of the ocean lithosphere treated as a problem of the cooling of a semi-infinite half space. (a) shows the thermal structure at the ridge (t0) where asthenosphere at temperature Tm is juxtaposed with ocean waters at temperature Ts. The thermal structure at successive distances away from the ridge where cooling of the initial temperature discontinuity in the semi infinite half space has lead to thickening of the ocean lithosphere is shown by the curves t1 and t2. (b) shows the error function (erf) and complimentary error function (erfc = 1 - erf).

From the half space model, theoretical predictions about the temporal evolution of lithospheric thickness, heat flow and bathymetry can be derived from the basic equations governing the thermal evolution of the lithosphere. Following Turcotte and Schubert (1982, p. 164-165, p. 181-182) the thickness of the lithosphere, zl, as a function of age, t, is given by:

zl = 2.32 
k t
 
(7.2)
The depth of the ocean floor beneath the ridge crest, w, at time, t, after formation is given by:
w = rm a Tm
rm-rw
  
 


k t
p
 
(7.3)
where a is the thermal coefficient of expansion and rm and rw are the density of mantle and water, respectively. The surface heat flow, qs, at time, t, after formation is given by:
qs = k Tm  1

p k t
 
(7.4)
where k in the thermal conductivity.

Some interesting consequences arise from the behaviour described by these equations. For example, Equation 7.3 shows that the average depth of the ocean is proportional to its crustal age. Given a constant volume of sea water, a change in the average age and hence depth of the oceans must result in a change in sea level, which is reflected in the geological record by the extent of ocean onlap on the continents. The average age of the oceans is inversely proportional to the rate of sea floor spreading, and directly proportional to the square root of the mean age of subduction. The mean depth of the oceans, w, as a function of the average age of subduction, t, is given by

w = 1
t
 
t

0 
 w  dt
(7.5)
substituting for w in Eqn 2.2 gives:
w = rm  au Tm
3 (rm-rw)
  
 


k t
p
 
(7.6)
Secular variations in the rate of sea floor spreading, reflected in the mean age of subduction may therefore have important implications to the average height of the oceans. Indeed, this explanation has been used to account for high ocean stands during the Cretaceous (when sea level may have been up to 300 m higher than today) which correspond with periods of fast ocean floor spreading (as indicated by analysis of ocean floor magnetic anomalies).

The thermal plate model

The semi-inifinite half space model predicts continuous cooling (albeit at a rate that gradually decays with time) and therefore thickening of the lithosphere through time (Figure 7.1). While the predictions are in remarkable agreement with the observations on bathymetry and heat flow in young ocean lithosphere these relationships appear to break down in ocean lithosphere older than about 80 Ma, when the thermal structure of the oceanic lithosphere appears to be stabilised.

Figure 7.3: Thermal structure of the ocean lithosphere predicted by the semi-infinite half space model (a) and by the thermal plate model (b).

The semi-infinite half space model assumes that the half space is not convecting. In the earth the deep mantle is convecting, with the consequence that a convective heat flux is provided at the top of the convecting layer (see Chapter 7). The thermal plate model accounts for this apparent time independent behaviour of old oceanic lithosphere by assuming that the convection in the subjacent mantle provides sufficient heat to the base of the cooling lithosphere to stabilise the cooling once a critical thickness is reached, the observations suggest this critical thickness is about 125 km corresponding to the thickness of  80 Ma old lithosphere (Figure 7.1b). Simply stated, the oceanic plate structure is thermally stabilised when the convective heat supply to the base of the lithosphere balances heat lost through the lithosphere.

7.2  Force balance on the ocean ridge

For young ocean lithosphere the cooling of a semi infinite half space provides an acceptable approximation and therefore Eqns 7.1 - 7.3 can be used as the basis to calculate the force balance on the ocean ridge. The isostatic compensation of the oceanic lithosphere causes the youngest ocean to form a high, albeit submerged, mountain range standing out above the abyssal plains. Such profound topographic gradients necessarily lead to substantial horizontal buoyancy forces (Chapter 2), termed the ridge push. In this section we provide the methodology for calculating the magnitude of the ridge push.

Figure 7.4: Ridge push is the force resulting from isostatically induced topographic gradients in the oceanic lithosphere. (a) shows a diagramatic representation of the force balance between a ridge and ocean lithosphere at age t = t1. (b) shows the density distribution appropriate to (a). The dashed line shows the vertical density structure at t = 0. The solid line shows the vertical density structure at t = t1. (c) shows the graphical solution to the force balance.

Refering to Figure 7.2a the ridge push, FR, operating on the oceanic lithosphere of age, t1, and depth below the ridge crest, w, is given by:

FR = F1 - F2 - F3
(7.7)
which is equivalent to solving Eqn 2.3, as shown diagramatically in Figure 7.2c. Note that F1 corresponds to the outward push of the asthenosphere beneath the mid-ocean ridge while F2 and F3 correspond, respectively, to the push of the water column and the old ocean lithosphere inward against the ridge.

Figure 7.5: Ridge push, FR, plotted as function of depth of ocean beneath ridge crest, w.

The quantitative evaluation of Eqn 7.7 is given in the Appendix A.3. The solution of Eqn 7.7 for any depth, w, below the ridge crest is shown in Figure 7.2, assuming the following physical properties a = 5 x 10-5, rm = 3300 kg m-3, rw = 1000 kg m-3, Tm = 1250oC and g = 10 m s-2 .

Figure 7.6: The oceanic lithosphere is characterised by its conductive geothermal gradient. The thermal gradient of the asthenosphere which is relatively well mixed (probably due to convection) is largely the isentropic temperature (adiabatic) gradient of the mantle due to volume and heat capacity (Cp) changes with changes in pressure (depth) The temperature at which this adiabatic temperature gradient extrapolates to the earth's surface is refered to as the potential temperature. Note that with the lithosphere of normal thickness (125 km), the solidus of the mantle peridotite nowhere intersects the geotherm but that on rifting of the lithosphere, the decompressed asthenosphere's adiabat intersects the solidus at about 40 km depth.

7.3  Formation of the oceanic crust

The ridge push resulting from the topography of the ocean floor, and the density structure within the oceanic lithosphere, provides (along with slab pull) one of the primary driving forces for lithospheric motion. The ridge push effectively maintains the constant rupturing of the oceanic lithosphere, and its separation on either side of the ridges. An important result of this rupture of the lithosphere at the ridges relates to the decompression of the underlying asthenosphere.

Figure 7.7: P-T diagram showing melting field of garnet peridotite and adiabatic (isentropic) decompression paths for mantle with potential temperatures of 1280oC, 1380oC, 1480oC and 1580oC, respectively (after McKenzie and Bickle, 1988).

The decompression of asthenosphere beneath spreading ridges is so rapid that there is virtually no loss of heat per unit rock mass; the decompression is therefore isentropic. Since small volume increases occur during isentropic decompression there is necessarily a small decrease in the the heat content per unit volume, and hence temperature. The change in temperature with pressure at constant entropy defines the adiabat. The entropy, S, volume, V, pressure and temperature of a system are related by the Clausius - Clapeyron equation:

DP
DT
= DS
DV
(7.8)
During isentropic (adiabatic) decompression, the decrease in pressure is accompanied by only small volume increases and thus T must decrease. This adiabatic gradient is about 1oC/km in the solid mantle. If the system becomes partly molten, then the change in volume with pressure is larger and T changes more quickly. Since the temperature of the convective mantle is not constant but lies on an adiabat we characterise it by its potential temperature (Tm), which is the projected of the adiabat to the surface of the earth (i.e., at 1 atm).

If sufficient decompression occurs, melting of the asthenosphere will take place when the adiabat passes through its solidus (Figure 7.3 and 7.3). The melt generated by this decompression has the composition of MORB and provides the parental liquid for all igneous rocks that make up the oceanic crust.

The amount of melting generated due to decompression of asthenospheric mantle beneath an active ridge segment depends entirely on the thermal structure of the asthenosphere and the melting properties of the mantle as a function of pressure. For the present day thermal structure (Tm = 1280oC) the amount of melting during complete decompression, amounts to a vertical column some 7 km thick (Figure 7.3). In the past, when the internal temperature may have been considerably hotter than it is today, the column of melt generated beneath the ridges, and hence the oceanic crust, may have been significantly thicker than 7 km.

Figure 7.8: Thickness of melt (measured as a vertical column in kms) present below the given indicated depth produced by the adiabatic decompression of garnet peridotite for different potential temperatures (after McKenzie and Bickle, 1988). For the modern day mantle with a potential temperature of 1280oC melting will not occur at depths les than about 45 km. Adiabatic decompression of the modern day mantle by the complete removal of the overlying lithospheric "lid" (for example at a spreading ridge) will result in the generation of a 7 km pile of MORB-like melt (i.e., the oceanic crust).

7.4  Coupling of the -spheres?

Equation 7.3, derived entirely from theoretical considerations, is in excellent agreement with observed bathymetry of ocean lithosphere younger than about 80 Ma. Indeed, this remarkable agreement between observations and the age-heatflow-bathymetry relationships predicted by Eqns 7.1 - 7.3 provides one of the principal lynchpins of the plate tectonic paradigm and one of the most persuasive lines of argument that the lithosphere is indeed thermally stabilised. Morover, it suggests that the motion of the oceanic lithosphere is largely decoupled from the flow in the underlying asthenosphere. There is as yet no clear understanding of the location, or even the general planform of mantle upwelling in the asthenosphere. Most importantly, we have shown that there is no requirement that the ocean ridges represent the site of upwelling (Figure 7). Wherever asthenospheric upwelling occurs it is likely to modify the thermal structure of the overlying lithosphere, and the suggestion is that the thermal structure of most old oceanic lithosphere has been modified to some degree by upwelling from within the underlying mantle.

7.5  Oceanic basalt chemistry

The oceanic lithosphere is generated entirely by magmatic processes, most of these concentrated at the mid ocean ridges, but with small, but scientifically interesting additions at intraplate hot spots forming the ocean island and sea mount chains (eg Hawaii, the Azores, Reunion, Iceland etc). Seafloor spreading generates about 20 km3a-1 of Mid Ocean Ridge Basalt (MORB) by the decompressional fusion mechanism described above (Figures 7.2-7.3), making these by far the most important volcanic provinces on Earth. From the perspective of this course we are most interested to know what oceanic magmatism tells us about the chemistry and äges" of their mantle source regions.

With respect to important trace element concentrations and the abundance of radiogenic isotopes, the basaltic rocks of the ocean islands (OIB) differ quite significantly from MORB. As oceanic magmatism must be sampling the sub-lithospheric mantle this immediately indicates that there are chemically distinct regions of the oceanic asthenosphere and the differences in Nd, Sr and Pb isotopic compositions (see Figure 7.5) are such that these separate regions must have remained isolated for timescales of the order of 2 Ga or more. This in turn places limits on the type of mantle dynamics and convection that must have occured.

As discussed in the previous chapter, MORB has the Sr, Nd (and Pb) isotopic characteristics, as well as incompatible trace element chemistry, consistent with derivation from a source which underwent the extraction of significant fractions of melt at some period in the past. Broadly speaking the extent of depletion of Rb, U and Nd in the oceanic upper asthenosphere as reflected by present-day MORB is consistent with being the compliment of the continental crust. Clearly, continental crust is not being produced at the present mid oceanic ridges and for the lithophile elements to find their way from the sub-oceanic mantle to the continents we must invoke at least a two-stage process. This second stage process may involve subduction and arc formation. Alternatively, it is possible that the main continental growth occurred in the Precambrian and was unrelated to the Earth's present geodynamic mode.

Figure 7.9: The range of radiogenic isotopic compositions (Nd, Sr, Pb and He) of oceanic basalts. The stippled area is the compositional range of Mid Ocean Ridge Basalts (MORB). The remaining envelope is the region of ocean island basalts (OIB). Clearly the OIB sources of the whole Earth are not all the same (from Allegre, 1987)

Rare Gas Studies

The atmophile inert gases, particularly He, Ar and Xe have recently provided substantial evidence for the existence, nature and life-span of important terrestrial reservoirs. These elements are of course strongly concentrated in the Earth's atmosphere. Perhaps surprisingly however, recent careful measurements of gasses emmitted during volcanic eruptions, included in fresh volcanic rocks (glasses) and even produced from hot springs and wells in continental and oceanic regions, reveal that detectible amounts of these and other gasses are still being released from the earth's interior. Furthermore these often show significant isotopic differences from the atmospheric reservoir indicating that emissions are tapping source regions which have been isolated for significant portions of the Earth's history. 40Ar is the radiogenic product of the decay of 40K and 36Ar is a stable isotope. 4He is the product of various decay series of U and Th and 3He is a stable isotope. The atmosphere has extremely high ratios of 4He/3He (722,000) and 40Ar/36Ar (295) because it is dominated by the accumulation of these products of radioactive decay of lithophile elements.

MORB has high 4He/3He ratios, though values are lower than those of the atmosphere, and extremely high 40Ar/36Ar. This is interpreted to indicate that the source region of MORB was effectively purged of most of its rare gas content early in the earth's history and that most of the He and in particular, the Ar there now is derived from subsequent U-Th and K decay. The existence of significant 3He does however indicate that the early degassing of the earth's interior (probably during core formation soon after accretion) was not total and that some primitive reservoir in the mantle has survived. This conclusion is supported by the 129Xe results which show that oceanic sources (OIB and MORB) are still producing this isotope. 129Xe is the product of the decay of 129I which has a half life of only 16 x 106 years and this decay scheme therefore became extinct very early in the earth's history. By contrast with MORB, some OIB have very low 4He/3He and 40Ar/36Ar ratios, reflecting that a component of their source region is primitive and had not undergone early degassing and melt extraction.

Chapter 8
Subduction and arc formation

The age-bathymetry-heatflow relationship in the ocean basins reflects the time dependent cooling and densification associated with lithospheric thickening (and, moreover, provides strong support for the notion that the oceanic lithosphere is thermally stabilised). Since the great proportion of old oceanic lithosphere is cold mantle peridotite, after some critical time the oceanic lithosphere must become more dense than the hotter peridotite of the subjacent convective asthenosphere thus allowing the possibility of subduction. The gravitational body forces acting on old dense ocean lithosphere give rise to the second fundamental driving force for the motion and deformation of the lithosphere, termed slab pull. The subduction process leads to melting of both the subducting slab and the overlying mantle wedge; the resultant melts segregate to form magmatic arcs (island arcs when they are built on oceanic lithosphere and Andean or Cordilleran arcs when they are built on continental lithosphere).

8.1  Buoyancy of the ocean lithosphere

The condition required for the onset of subduction is the negative buoyancy of the ocean lithosphere, i.e.:
1
zl

zl

0 
 rz  rm
(8.1)
where rz is the density at depth z within the lithosphere of thickness zl and rz is the density of the convective mantle. Assuming an appropriate density structure for the ocean lithosphere we can solve for the critical lithospheric thickness, zlc, for which Eqn 8.1 is true. We begin by assuming the simplest of possible density structures in the ocean lithosphere, as shown in Figure 8.1:
rz = rc
0 < z < zc
rz = rm

1+a Tl  z
zl


zc < z < zl
(8.2)
where rc is the density of the crust. The appropriate equality for zlc is:
1
zlc
 
zlc

0 
 rz dz = rm
(8.3)
Integration of Eqn. 8.3 gives:
1
zlc


rc  zc+(zlc-zc)

rm a (Tl-Tc)
2




= rm
(8.4)

Figure 8.1: Density structure assumed for solution of Eqn 8.1 where Tc is the temperature at the Moho.

For a linear geotherm Tc is given by:

Tc = Tl

zc
zlc


An identical solution is readily obtained graphically by simply making the equivalence between the two regions in Figure 8.1:
zc  (rm -rc) = rm a Tl (zlc-zc)2
zl  2
(8.5)
Substituting in the following values appropriate to the ocean lithosphere : zc = 6000 m, rc = 3000 kg m-3, rm = 3300 kg m-3, Tm - Ts = 1250oC, a = 3 x 10-5 oK-1, gives the critical lithosphere thickness of 40 km. Substituting zlc into Eqn 8.3 yields the critical age of 10 Ma, for the onset of negative buoyancy. In this calculation we have assumed that the only compositional changes associated with the production of oceanic lithosphere are within the oceanic crust due to the basalt extruded at the ridge axis. However, the extraction of basalt from the underlying peridotite must leave a depleted "harzbugitic" residual mantle peridotite, the density of which is somewhat less than the primary mantle peridotite by an amount relating to the proportion of melt extracted (50 kg m-3 for 25% extraction). This "depleted" zone extends beneath the ridges to the base of melting (about 40 km for normal mantle with a potential temperature of 1280oC, Figure 7.3), and is eventually frozen onto the base of the cooling lithosphere. Calculations which account for this "depleted" nature of the upper 40 km of oceanic mantle lithosphere show that the age for onset of negative buoyancy is significantly greater than calculated assuming the density structure as in Eqn 8.2 and Figure 8.1. For example, Oxburgh & Parmentier (1977) calculate a critical age of approximately 50 Ma.

The above analysis shows that much of the old oceanic lithosphere will have the potential for initiating subduction, with the precise age dependent on the average density of the oceanic lithosphere. Negatively buoyant lithosphere may however be held up by the elastic strength of surrounding buoyant lithosphere, and the actual mechanism which triggers the onset of subduction is poorly understood. Moreover, younger positively buoyant oceanic lithosphere may be subducted if it is dragged into pre-existing subduction zones by the pull of the older subducting slab and the push of the existing ridge. Indeed, some workers have suggested it is possible that ocean ridges can be subducted, even though they are necessarily positively buoyant.

8.2  Therml structure of subducted slabs

The penetration of the mantle by the subducting slab must disturb the thermal regime of the mantle, providing the slab does not heat up by diffusion at a rate fast compared to subduction. The question, then, is whether the advective or the diffusive terms in Eqn 4.3 dominate the thermal evolution of the downgoing slab. Simple application of the dimensional analysis outined in section 4.2 allows us to resolve this problem by determining the thermal Peclet number appropriate to subduction (if PeT > 10 then advection dominates weheras if PeT < 1 diffusion dominates). The thermal Peclet number is given by :

PeT = v  l
k
v is given by subduction velocities ( > 1 cm/yr), l by the approrpiate length-scale across which diffusion must act, which is here defined by the thickness of thelithosphere (100 x 103 m) and k by the thermal diffusivity (10-6 m2 s-1. The PeT appropriate to subduction is therefore about 30, thus implying that subduction advects the thermal structure of the ocean lithosphere into the deep mantle.

8.3  The magnitude of slab pull

Subduction occurs due to the gravitational body forces acting on the descending slab, which are proportional to the density contrast between the slab and the surrounding hotter mantle. This density contrast, which drives subduction, arises because of (Figure 8.2): The force balance on the descending slab must also include a term related to the viscous resistance of the mantle to the downward movement of the slab, which will be proportional to subduction velocity. Because the rheological properties of rocks at high temperatures are poorly known, this viscous term in the force balance cannot yet be quantified to any reasonable levels of confidence. An analytical formulation of the slab pull is therefore far more complicated than for ridge push, and here we treat the problem in a qualitative way only.

Figure 8.2: Scematic force balance on a subducting.

That the upper parts of subducting slabs are in general descending under the influence of a gravitational body force is indicated by the existence of extensional stresses in the descending slab as witnessed by focal mechanisms of earthquakes generated in the upper 100 - 200 km of the slab. This implies that for most subduction environments the net gravitational body force minus the viscous resistance (which combine to give slab pull) exceed the magnitude of the ridge push) for normal velocities of subduction (typically in the range 0 -10 cm yr-1). The lack of subduction at rates significantly greater than 10 cm yr-1 may be understood by the increase in the viscous resistence for increased velocities: the viscous resistance providing a negative feedback on the descent velocity.

8.4  Arc dynamics

Kinematics

From the point of view of lithospheric kinematics, subduction represents one of two mechanisms accommodating convergence between plates; the alternative mechanism being the internal deformation of the plates. The relative importance of subduction and internal deformation can be illlstrated by considering the plate convergence velocity Vcon relative to subduction velocity Vsub as in Figures 8.3 & 4.

Figure 8.3: Convergence between plates, Vcon = V2 -V2, may be taken up either by subduction Vsub or by internal deformation of the plates as illustrated in Fig. 8.4.

Figure 8.4: Alternative kinematic scenarios for convergent orogens expressed terms of the velocity of subduction relsative to velocity of plate convergence

Mechanics

While the subduction environment is one of convergence, it may not necessarily involve horizontal compression. We have shown that the driving force for subduction is provided by the negative buoyancy of old oceanic lithosphere which imparts tensional stresses in the descending slab. The state of stress at the trench will depend primarily on the buoyancy of the oceanic lithosphere that is being subducted. We have already commented that it is possible to subduct relatively young buoyant lithosphere. Moreover, there are regions of elevated oceanic lithosphere where the crust is anomalously thick due to, for instance, the accretion of Ocean Island Basalts (OIB). Examples of subduction of anomalously bouyant oceanic lithosphere are provided by the Juan Fernandez and Nazca ridges presently being subducted beneath the amagmatic zones of the South American Andes. Siesmicity from the slab descending beneath these regions shows that this anomalously buoyant oceanic lithosphere is subducting at a much shallower angles than adjacent 'normal' oceanic lithosphere. The state of stress at the trench may therefore be compressional if young or anomalously buoyant lithosphere is being subducted or it may be tensional if old negatively buoyant lithosphere is being subducted.

The force applied at the trench (Ft) provides a boundary condition for deformation in the neibouring arc, which is also subject to buoyancy forces (Fb) arising form the topographic gradients on the arc itself (Chapter 3). The resultant force balance on an arc (Fa) can be formulated as :

Fa = Ft - Fb
(8.6)
where Ft and Fa are positive for compression and Fb is the bouyancy force exerted by the arc on the surrounding lithosphere.

Because of difficulties in formulating Ft we do not provide a complete solution here (one way of treating Ft is in a completely ad hoc fashion by considering it as variable within the bounds of slab pull and ridge push). However, substantial horizontal extensional forces can arise from topography of the magnitude observed in arcs. Extensional failure of an arc gives rise to the development of small ocean basins, termed back-arc basins, between ruptured arc fragments. The foundering of old dense oceanic lithosphere in front of an active arc will lead to a "roll back" effect of the trench, which will greatly enhance the tendency of an arc to spread (Figure 8.4.).

8.5  Arcs and crustal growth

The magmatic island arcs and continental marginal arcs represent the most important sub-aerial terrestrial volcanic provinces. In these settings the visible volcanic activity is volumetrically minor in comparison with contemporary intrusive activity. Compositionally the characteristic lavas of these terrains are andesitic (or quartz diorite, tonalite, granodiorite or quartz monzonite intrusives), though basalts (gabbros) are also common. These zones of convergence between either continental - (e.g. the Andes) or oceanic lithosphere (island arcs) and oceanic lithospheric plates, are sites which are probably critical to the evolution of the continental crust (and subcrustal lithosphere). The ändesite" model of Taylor (1977) embodies this concept, suggesting that the average composition of at least the upper continental crust is andesitic in bulk geochemical terms and that it aquired this composition because the main mode of continental growth has been by marginal accretion of magmatic arcs.

As we have seen in the previous chapters, the oceanic crust has been developed by extensive and repeated fusion of the asthenosphere over a substantial portion of the Earth's history (at least the past 2 Ga). This source of MORB has very large-scale depletions of lithophile (incompatible) trace- and minor elements in comparison with the pristine or primitive (chondritic) mantle. The low 87Sr/86Sr, 207Pb/204Pb, 206Pb/204Pb and high 143Nd/144Nd values of MORB imply that its source region has had low Rb/Sr, Nd/Sm and U/Pb ratios for substantial periods of Earth history. The fact that the continental crust has complementary isotopic and geochemical characteristics with those of MORB suggests that the crust represents the "missing-part" of this asthenospheric source of present-day MORB. The question of what role arcs play in mediating this transition of ocenaic crust to continental crust is fundamental to understanding crustal growth.

Figure 8.5: A possible relationship between MORB and continental growth. If the important events which influence the complementary geochemical /isotopic evolution of the sub-oceanic upper mantle and the continental crust do in deed take place at plate boundaries (and there is some argument against this), then the critical boundary must be the zones of lithospheric plate convergence (subduction zones). It is clear that although the primary extraction of material from the asthenosphere takes place at the ridges the basalt produced here does not resemble the continental crust , the upper portion of which is certainly much more felsic. The following section will discuss how the quite special circumstances that exist in subduction zones combine to produce and emplace relatively felsic magmas in the crust , and why these magmas probably contain a large part of the incompatible element content of the oceanic crust (and hence of the sub-oceanic upper mantle).

Figure 8.6: Nd and Sr isotopic composition of arc magmas and their potential sources.

The sources of arc magmas

One complicating factor in the study of arc magmatism is the question of the source of the magmas. It is easy to show that the covariation of isotope and geochemical data for many suites of arc volcanic rocks must be due to mixing of components from more than one source and in fact many studies can demonstrate up to four or five seperate components. The possible source components which may be available include:
  1. The basaltic (crustal) portion of the subducting slab.
  2. The pelagic sediments which may be carried with the subducting slab.
  3. The depleted peridotite wedge above the subducting slab.
  4. The lower crust of the arc itself.
  5. An enriched (or none-depleted) mantle source (like the OIB source).
Added to this list is the important fact that the oceanic crust undergoing subduction is not pristine, but rather has undergone very extensive alteration, hydration and metamorphism since the time of its formation at the ridge. The upper portion of the subducting slab will therefore undergo dehydration at depth generating fluids which themselves carry considerable concentrations of solutes. These will escape from the subduction zone and their role in fluxing the overlying mantle will have a pivotal role in the generation of arc magmas. The hydrothermal alteration is promoted by the initially high thermal gradient at the ridge crest which promotes circulation of sea water in fracture systems that might extend through the thickness of the whole crust. The oceanic crust is geochemically and mineralogically modified by exchange with sea-water which changes (amongst other things) its O- and Sr-isotopic composition towards those of sea water. The primary igneous minerals and glass are altered to hydrates including talc, actinolite-tremolite , epidote ,chlorite, serpentine and brucite and these minerals carry their bound water to be released at depth from the subducting slab.

The thermal structure of arcs

The heat-flow along an across arc transect from the oceanwards-side of the trench to the back-arc region varies considerably. There is low heat flow above the old oceanic crust, extremely low heat-flow from the trench and significantly higher flow from the arc and backarc regions. The extremely low heat flow at the trenches is an expression of the persistent coldness of the the downgoing slab to depth. This has the importance of delaying the onset of those dehydration reactions that are mainly temperature dependent to much greater depths than under the regime of normal geothermal gradients. Particularly significant amongst these is the the reaction (Figure 8.5):

amphibole - > Na - clinopyroxene (omphacite) + garnet + H2O

which defines the amphibolite - eclogite transition and at pressures less than 25-30 kbars. this takes place at close to 1000oC. Above 25 kbars the reaction is much more pressure dependant (Figure 8.5).

Figure 8.7: The stability field of the hydrate amphibole as a function of temperature and pressure.

Figure 8.8: A simplified view of the thermal structure of the subducting slab and the sub-arc mantle. The stippled region illustrates the stability field of amphibole. The break- down of amphibole when the the slab reaches about 100 km releases H2O and produces the dense rock eclogite consisting of garnet and omphacite (r = 3600 kg m-3).

Melting in the Subducting Slab

The crustal part of the downgoing slab has the potential capacity to produce felsic melts with quite high lithophile element contents, high 87Sr/86Sr, 207Pb/204Pb, 206Pb/204Pb and low 143Nd/144Nd. This is partly due to the pelagic sediment component and partly to the hydrothermal alteration of the basaltic oceanic crust. A popular theory when subduction was first recognised was that arcs are dominated by intermediate (andesitic) magmas because they were derived from melting of mafic sources rather than the ultra-mafic sources that generate basaltic magmas, i.e. the upper part of the down-going slab. However, the common presence of basaltic magmas in arcs implies that the slab is not the major source of these melts . Rather it seems more likely that the wedge of peridotite mantle above the slab is the source of the bulk of the arc magmas and this melts at relatively low temperatures because of the fluxing by water -rich fluids driven off the slab. The mantle wedge is also ënriched" with respect to normal mantle peridotite in the lithophile (incompatible) element component derived from the slab. The hydrous, mafic, primary magmas thus produced then undergo variable fractionation en-route to the upper part of the crust to produce andesite, dacite and rhyolite.

The fate of the slab after it has suffered the extraction of small volumes of melts and fluids is of interest. Ringwood (1974) has determined the density variation of the main earth materials in the subduction system. As illustrated in Figure 8.5, the subduction zone can be considered as a two layerd slab composed of MORB crust and depleted peridotite mantle lithosphere (harzburgite) which is being thrust into asthenospheric mantle. Under conditions of thermal equilibrium, at pressures equivalent to depths up to 650 km, MORB is more dense than both asthenospheric and harzburgititic (lithospheric) peridotite, but is less dense than both types of peridotite in the depth interval 650- 720 km. Depleted, harzburgitic peridotite is less dense than less refractory peridotite (asthenosphere) at all depths except in the 650 - 720 km interval.

Figure 8.9: Compositional structure of the descending slab.

On this basis, Ringwood (1974) suggests that the subducting slab, reaches a neutral buoyancy level at about 650 km where it is preserved as a layer composed of harzburgite with remnant blobs of the basaltic (MORB) portion of the slab. He further suggests that this accumulated layer would provide a physical barrier to promote the development and maintenance of a two - layered convection system in the mantle.

The role of water in arc magmatism

Allegre et al (1987) calculate that over the past 4 Ga. subduction has recycled about 4 times the mass of the present oceans back into the mantle together with considerable amounts of the atmospheric rare gas contents. It is clear that the mantle residence time of this massive flux is essentially zero and that subduction zone volcanism is the obvious site of outgassing of this system. The direct evidence for the role of water in arc magmatism comes from many observations. These include:
  1. Direct discharge of water vapour (with other fugitive species) during volcanic eruption, leading to explosive and pyroclastic volcanism.
  2. High plagioclase content of intermediate magmas resulting from its latent enrichment during fractionation of hydrous magmas. Water delaying the normal appearance of plagioclase to lower temperatures.
  3. The expansion of olivine crystallisation to relatively silicic magmas.
  4. The occurence of hydrous minerals such as amphibole and biotite.
  5. The discharge of magmatic hydrous fluids from decompressing and crystallising magmas at upper crustal pressures to contribute to hydrothermal alteration halos around some plutons

    (porphyry copper deposits).

  6. The delayed, late, or none crystallisation of quartz in many quatrz-normative magma systems.
  7. The complex oscilliatory zonation/ resorbtion patterns of many phenocrysts in arc magmas suggesting growth variable PH2O regimes in crustal magma chambers.
  8. Vapor-bearing fluid inclusions in many phenocrysts.

Figure 8.10: Peridotite solidi for dry, H2O-undersaturated and H2O -saturated (wet) conditions. The stippled area indicates the possible range of geothermal gradients.

It is probable that melting in the peridotite wedge above the subducting slab takes place under conditions like those of the water undersaturated solidus in Figure 8.5. Melting is promoted by the high activity of water from the dehydration of the slab and the initiation of melting will promote the rise of the geotherms in the strip of asthenosphere above the main locus of dehydration in the slab. This in turn will contribute to the rise of the asthenosphere beneath the arc, leading to thinning of the mantle lithosphere beneath arc and further decompressive melting in the asthenosphere. At this stage the situation is not unlike that at the mid ocean ridges and arc splitting and spreading may take place. The combination of the effects of the removal of the mantle lithosphere from beneath the arc and the addition of large volumes magmatic rocks to the crust will tend to depress the base of the arc's crust through rising geotherms, to the effect that the the solidus of mafic rocks may be intersected at lower crustal depths. This will lead to short-timescale magmatic recycling within the arc itself.


Part 3
The continents


Chapter 9
Continental Deformation

The structural geometry of both convergent and extensional continental orogens at the outcrop scale is very (some would say - horribly) complex. More than anything else, this complexity reflects the very strong mechanical anisotropy of crustal rocks; that is, the structural geometry resulting from the deformation is largely a consequence of the inherited structure and is only weakly coupled to the nature of the forces driving the deformation. Indeed this complexity begs the question as to whether it is possible, or even useful to attempt, to evaluate the parameters governing the geodynamic evolution of continental orogenic belts. However, some of the most impressive large scale features of orogenic belts are much more regular than their internal structural geometry. For example, the topography of orogenic belts, while very fragmented (fractal) on small scales, is very regular at the scale of the orogenic belt. Indeed, just as an understanding on the control on topographic variation in the ocean basins provides fundamental insights, understanding the controls on topography provides very important insights into the mechanics of continental orogens.

We begin by examining the controls and some consequences of topography and potential enegrgy using simple calculations based on the assumption of local isostatic equilibrum. This assumption is only likely to be valid for thermally mature orogenic systems which have started to develop plataeus (e.g., Tibet) or, in extension, wide basins, in which the deformation of the lithsophere is induced by forces applied as end loads. The margins of orogenic belts involving wedge shaped thrust belts are certainly not in isostatic equilbrium and so the topographic variation in such circumstances need to be evaluated using a differnet set of boundary conditions. To tackle the mechanics of the frontal parts of mountain belts, and accretionar wedges, we need to investigate the dynamics of critical wedges, where the driving forces are imparted to the deforming crust as basal tractions along some kind of master thrust or decollement.

9.1  Deformation of the lithosphere subject to an end load

Any consideration of the behaviour of the continental lithsophere during deformation needs to recognise the contrasting influence of the crust and the mantle lithosphere in mediating both the thermal and isostatic response: Because of this contrasting influence it is useful to consider the effects of the crust and mantle lithosphere, independently. This can be achieved by describing the deformation of the lithosphere by the ratio of the changes in thickness of the crust, fc, and mantle lithosphere, fm, where these are defined at any stage of the deformation by:
fc = zc
zc0
fm = zm
zm0
where zc and zm are the deformed thickness of the crust and mantle lithosphere, respectively, and zc0 and zm0 are the initial thickness of crust and mantle lithosphere prior to deformation.

Airy Isostasy and crustal thickening

Assuming Airy isostasy then crustal thickening results in both an increase in surface elevation, ch, and the the development of the crustal root, cr.

Figure 9.1: Isostatic effects of crustal thickening.

For a crust of constant density, rc, overlying mantle of density rm, the ratio of the change in surface elevation ch to the thickness of the crustal root cr is given by

ch
cr
= rm - rc
rc
The change in crustal thickness is given by the sum of ch + cr:
fc = ch + cr + zc0
zc0
giving
cr = fc  zc0- ch - zc0
Therefore
ch = zc0(1 -fc )

rc-rm
rc


Note that the change in surface elevation is linear in the change in crustal thickness; that is, it is linear in fc  zc0. The change in potential energy Uc for this scenario is given by:
Uc = g  rc
2
((zc0 fc)2 -zc02 - 2 cr2) - g  rm
2
cr2
Note that potential energy changes with the square of the crustal thickening.

The effects of crustal thinning (i.e., fc < 1) are exactly opposite crustal thickening, that is it results in subsidence and a reduction in potential energy (see Chapter 10).

Airy Isostasy and the mantle lithosphere

See Sandiford & Powell (1990) EPSL

9.2  Deformation within the lithosphere due to basal tractions

Chapter 10
Stretching continents

The most important isostatic consequence of crustal deformation relates to the changes in the density structure of the lithosphere which dictates the isostatically-supported elevation of the lithosphere. For convergent deformations, involving crustal thickening, the isostatic response is (generally) uplift (or mountain building). For extensional deformations the isostatic effect is in part subsidence (or basin formation) and in part uplift, depending on:

Basins formed as the isostatic consequence of extensional deformations are termed stretch basins and include most continental passive margins formed during continental breakup events.

10.1  Isostatic calculations

Rift Phase Subsidence

The treatment of the isostatic consequences of an extensional deformation is similar to te treatment of the uplift caused by convergent deformation, only now we have also to consider the density distribution of the medium filling the basin. We first consider the effects of a homogeneous extensional deformation on the scale of the lithosphere with no magmatic additions (Figure 10.2):

Figure 10.1: Lithospheric-scale extensional geometry used to quantify subsidence.





.
e
 

zz 

z







x 
= 0
(10.1)
Using this condition we can parameterize the deformation at any point within the basin in terms of the one dimensional vertical strain:
b = 1
.
e
 

zz 
(10.2)
Which for plane strain ([(e)\dot]yy = 0) gives b = [(e)\dot]xx or, equivalently, the horizontal stretch across the basin. Note that b is the reciprocal of f, the thickening factor used in Chapter 9.

We can formulate the isostatic effects of active rifting in the following manner. Assuming an intial linear lithospheric geotherm, which amounts to ignoring the effects of any internal heat production, gives the temperature Tz at depth z (Figure 10.1) :

Tz = Ts + (Tl - Ts) z
zl
0 < z < zl
(10.3)
where Tl is the temperature at the base of the lithosphere, Ts is the temperature at the surface of the lithsophere and zl is the thickness of the lithosphere prior to stretching (Figure 10.1). Letting Tm equal the temperature difference across the lithosphere, and setting Ts = 0, gives:
Tz = Tm z
zl
0 < z < zl
(10.4)
After homogeneous stretching by a factor b the thermal structure is given by:

Figure 10.2: Temperature (a) and density (b-d) structure assumed in determining subsidence in a stretched basin. (b) shows denstiy structure prior to stretching, (c) shows density structure immediately after instantaneous stretching, (d) shows density structure after thermal sag.

Tz = Tm b  z
zl
0 < z < zl
b
Tz = Tm
zl
b
< z < zl
(10.5)
The initial vertical density structure is given by (Figure 10.1):
rz = rc

1+a Tm

1- z
zl




0 < z < zc
rz = rm

1+a Tm

1- z
zl




zc < z < zl
(10.6)
where rc and rm are respectively, the density of the crust and mantle at the temperature Tm and a is the volumetric coefficient of thermal expansion. Following stretching by b the density structure is:
rz = rf
0 < z < Si
rz = rc

1+a Tm

1-b
zl




Si < z < Si + zc
b
rz = rm

1+a Tm

1-b
zl




Si + zc
b
< z < Si + zl
b
rz = rm
Si + zl
b
< z < zl
(10.7)
where rf is the density of the medium filling the rift basin of thickness Si. The condition of isostasy amounts to equating the vertical stress, szz, at a common depth, zl, beneath the unstretched and stretched basin. Beneath the unstretched lithosphere szz at zl is given by (allowing that rc  ac = rm  am):
szz \midz = zl
=
g
zl

0 
rz  dz
=
g zc  rc + g (zl - zc) rm + zl
2
 g a rm Tm
(10.8)
Beneath the stretched basin szz at zl is given by:
szz \midz = zl
=
g
zl

0 
rz  dz
=
rf  Si + g zc rc b+ g (zl-zc)rmb
+ g zl a rmTm
2b
+ g rm

zl-Si- zl
b


(10.9)
equating Eqn 10.8 with Eqn 10.9 and solving for Si gives (Figure 10.1):
Si = (b-1)(2 zc(rc-rm)+Tm a zl rm )2 b (rf-rm)
(10.10)

Figure 10.3: Rift phase subsidence as a function of stretching for a basin completely filled by sediment (rf = 2500 kg m-3), and water (rf = 1000 kg m-3).

Sag Phase Subsidence

We assume that in the steady state the thickness of the lithosphere is dictated by balancing the rate of heat flow through the lithosphere with the heat supplied to its base by the convective motion in the interior. Any deformation involving changes in the thickness of the lithosphere therefore induces a departure from this equilibrium condition which in turn influences the ensuing thermal evolution of the lithosphere. Since heat loss is proportional to the temperature gradient in the lithosphere extensional deformation brings the asthenosphere closer to the surface and consequently increases the heat loss through the lithosphere (Chapter 4). This imbalance causes heat to be lost through the lithosphere more quickly than it is supplied. The consequent cooling and thickening of the lithosphere causes it to become more dense, inducing a kind of thermal subsidence or sag. There are two questions of interest here: We can formulate the isostatic effects of the thermal or sag phase subsidence, St, following a stretching by a factor b in the following manner. Following reestablishment of the equilibrium lithospheric thickness, zl, by conductive heat loss the density distribution will be:
rz = rf
0 < z < St + Si
rz = rc

1+a Tm

1- z
zl




St + Si < z < St + Si + zc
b
rz = rm

1+a Tm

1- z
zl




St + Si + zc
b
< z < zl
(10.11)
where St is the subsidence due to thermal sag. Letting the total subsidence, S = St + Si, then szz at zl beneath the saged basin is given by:
szz|z = zl
=
g
zl

0 
 rz  dz
=
g  rf S + g zc rc  b+
g  rm

zl-S- zc
b
+  zl a Tm 
2


(10.12)
Solving for S gives:
S = zc(b-1)( rc-rm))
b(rf-rm)
(10.13)
and St (Figure 10.1):
St = Tm a rm zl (b-1)
b(rf-rm)
(10.14)

Heterogeneous stretching of the lithosphere

The fc-fl diagrams (note that f is the reciprocal of b) shown in Figure 10.1 allows the rift phase subsidence to be calculated for any inhomogeneous stretching deformation, in which the deformation is not equally partitioned between the crust and mantle (although in these diagrams there has been no allowance made of the density of the medium filling the depressions).

Figure 10.4: Sag phase subsidence as a function of stretching for a basin completely filled by sediment (rf = 2500 kg m-3), and water (rf = 1000 kg m-3)

The fc-fl parameterisation allows us to consider the isostatic effects of magma additions to the crust during extension, as for instance has occurred during stretching of the Basin and Range Province in the western USA. Magma additions add mass to the crust, with the result that fc > fl. In the Basin and Range Province magma addition appears to have been sufficent to maintain a crustal thickness in the vicinity of 30 km, while the underlying mantle part of the lithosphere has been severely attenuated (Gans, 1988). The resulting path shown in Figure 10.1 explains the anomalously high elevation of the Basin and Range Province. Finally, the isostatic effects of stretching are dependent to a large extent on the vertical density structure of the lithosphere prior to stretching. For typical continental lithosphere we have shown that homogeneous (or pure shear) thinning results in subsidence. In contrast, Figure 10.1 shows that homogeneous thinning of old oceanic lithosphere results in uplift or Ridge formation. Indeed, this is exactly as we showed in Chapter 7, and this is exactly why the ocean ridges are important engines for lithospheric motion and deformation.

Figure 10.5: fc-fl diagrams for stretching of typical continental lithosphere (zc = 35 km, zl = 125 km, elevation contours in 1 km intervals) and typical oceanic crust (zc = 6 km, zl = 125 km, elevation contours in 0.5 km intervals). Note that homogeneous stretching of continental lithosphere leads to subsidence (elevation contours for subsidence assume no sedimentation or water filling, to change to a filled basin multiply subsidence by rm/(rm - rf)) while homogeneous stretching of oceans leads to uplift (i.e., ridge formations).

10.2  Mechanical consequences of extension

Inttroduction

So far we have treated basin formation in terms of a kinematic model only, that is without consideration of the forces that drive the deformation. In order to develop a mechanical model we must consider the thermal evolution of stretch basins explicitly, because the mechanical strength of the lithosphere is closely allied to its thermal state. We use the mechanical model described in Chapter 2.

The time scale of thermal sag

As discussed in Chapter 4.3 a convenient measure of the response time to thermal perturbations in media in which heat is transported only by conduction is given by the thermal time constant, t:
t = l2
p2 k
(10.15)
where l is the length scale over which the perturbation occurs, k is the thermal diffusivity (units m2 s-1). The thermal time constant provides a measure of time for the decay of approximately 50% of the thermal perturbation, with measurable deacy occuring up until about 3 t. The thermal time constant of the continental lithosphere is of the order of   60 Ma (k = 10-6 m2s-1, l = 105 m), and thus we may expect that subsidence associated with sag of a rifted basin for a period of up to   150-200 Ma (Figure 10.2).

Since the decay of a thermal perturbation is proportional to the second derivative of the temperature gradient with respect to depth (Chapter 4) and is therefore an exponentially decreasing function of time. Consequently, the thermal sag subsidence should also decline exponentially with time.

Figure 10.6: Evolution of the geotherm during sag following instantaneous stretching at time t = 0.

Thermal evolution at finite extensional strain rates

The instantaneous thinning of lithosphere will result in isothermal decompression of all material points in the lithosphere. The resultant increase in the geotherm for an instantaneous stretch is b. However, lithosphere stretched at finite rates will evolve geotherms lower than the theoretical ïnstantaneous" geotherm for two reasons. with important consequence for the evolution of strength within the lithosphere .

Mechanical Evolution

We have shown that convergent deformation of the continental lithosphere induces weakening due largely to the increase of the Moho temperature associated with increasing heat production within the lithosphere. This provides a self localising mechanism for convergent orogens, with the length scales dictated by the rheology, strain rate relationship. Similarly the length scales for divergent deformation will be also controlled by rheology. The lithosphere is, however, weaker in tension than in compression because brittle failure in tension occurs at lower stress difference than in compression (Figure 10.2).

Figure 10.8: See text for discussion

Cconsequenlty, for an equivalent geotherm and strain rate the across strike length scale of an extensional orogen will be significantly less than a compressional orogen.

For small initial finite stretches the lithosphere will undergo essentially isothermal decompression (Figure 10.2). The consequent decrease in the depth of the power law failure envelopes for quartz and olivine will induce a reduction in the strength of the lithosphere; a form of thermal weakening (Figure 10.2a). However, with increasingly large finite strains, the cooling of material points as they are decompressed (Figure 10.2) ultimately leads to an increase in whole lithospheric strength (Figure 10.2). This analysis suggests that early in their evolution rift zones should be self-localising but gradually becoming more distributed. The point at which the switch from weakening to strengthening occurs depends critically on the strain rate with the critical strain decreasing with decreasing strain rate.

The role of magmas

Of course the thermal and mechanical evolution of extensional orogens is depenedent not only on the rate of deformation but also on the involvement of asthenospheric magmatic additions. Decompression of the asthenosphere beneath stretched lithosphere may induce melting in an analogous fashion to melting beneath ridges. For the normal potential temperature of the asthenosphere (1280oC) this requires decompression of asthenosphere beneath depths of approximately 45 km. Assuming initial lithospheric thicknesses of the order of 100 km or more the minimum b value required for melting of the mantle beneath the stretched lithsophere is about 2. Once melting occurs, the rapid advection of heat into the overlying lithosphere will cause dramatic strength reductions, potentially overriding the strain hardening associated with large strains mentioned in the previous section. Moreover, the isostatic consequences of magma addition are to increase the elevation of the lithosphere, as is the case in the Basin and Range Province, potentially giving rise to extensional buoyancy forces which augment the original driving forces for stretching. This magmatically augmented weakening may be an essential requirement for the complete rupturing of stretched continental lithosphere to form new spreading oceans.

The role of mantle plumes

In the discussion above we have not explicitly considered the origin of the forces that drive extension of the continental lithosphere. While it is clear that compressional forces within the continents can be generated by the tractions exerted by subducting slabs it is not at all clear that substantial extensional force may be exerted on normal thickness continetal lithosphere. This is because normal thickness lithosphere (35 km thick crust and 125 km thick lithosphere) is in approximate potential energy balance with the mid-ocean ridges; that is the ridges seem to neither exert compression nor tension on normal continental lithosphere. Of course regions of elevated topography within the continents, such as associated with thickened crust formed in zones of compression, will naturally experience tension when the forces driving convergence are relaxed. Another important way of increasing the potential energy of the lithosphere may be through the impingement of a mantle plume at the base of the lithosphere, which may have the effect of jacking up the lithosphere by as much as 1 km. The associated potential energy gain (for a column of unit area) is given to a first approximation by:
Uplume = g  rc h2
2
(10.16)
where h is the additional elevation caused by the jacking effect of the mantle plume, and for crustal density rc = 2800 kg m-3 is equivalent to a tensional force per unit length of topography of about 1.4 x 1012 N m-1 for every km of additional elevation.

An important additional effect of extension initiated by a plume relates to the higher mantle temperatures which may be up to about 300oC higher than the mean convective mantle temperatures. The implication is that decompression melting in the asthenosphere will begin at a much earlier stage in the rifting process. For example for a potential mantle temperature of 1580oC melting begins at about 120 km depth, implying that plumes may begin generating melts without any appreciable deformation of the overlying lithosphere. Moreover as mentioned in the previous section the transport of such melts into the lithosphere must further enhance the potential energy of the lithosphere and weakens it in such a way to augment any extensional deformation that has already begun.

Figure 10.9: Schematic representation of how lithospheric strength changes as a function of stretching. For low finite strains and/or high strain rates decompression is essentially isothermal and thus the plastic failure envelopes move to shallower depths (since they are dependent on temperature and thus only indirectly on depth through the geotherm - Eqn 5.2) inducing a reduction in lithospheric strength by an amount equal to Area 1 + Area 3. A decrease in the depth of the Moho results in a strength increase by an amount equivalent to Area 2. The net change in the lithospheric strength is given by Area 1 + Area 3 - Area 2. At low strain rates and/or high finite strains cooling of material points occurs such that there is no rise in the plastic failure envelopes with further deformation (i.e., Area 1 and Area 3 tend to 0).The consequent strain hardening is equivalent to Area 2. White and McKenzie (1989) point out that in mantle upwelling zones the temperatures may be in excess of   100oC above the typical potential temperature of the mantle. Stretching above such abnormally hot asthenosphere will induce melting at rather lower b values than for normal mantle, and for a given b value contribute a much greater amount of melt to the overlying lithsophere. The variable magmatic history of extensioanl provibnces and rifted margins may therefore relate to the thermal state of the subjacent convective interior during rifting.

10.3  Sedimentation in stretched basins

We have seen that the evolution of rift basins can be viewed in terms of two distinct phases: the rift phase, and the sag phase. In the the rift phase sedimentation is associated with the isostatic subsidence induced by the active deformation. Rift phase sediments show signs of active tectonism, such as growth faults with rapid facies and thickness changes. Sag phase sedimentation follows rift phase and represents the isostatic effects of freezing asthenosphere onto the base of the thinned lithosphere as the geotherm returns to steady state. Sag-phase sediments are typically laterally extensive without abrupt facies changes. Subsidence rates decrease exponentially through the sag phase with a characteristic time scale indicative of the thermal time constant of the lithosphere.

The Steer's Head

Stretched basin stratigraphy commonly exhibits a Steer's Head or Texas Longhorn shape with the sag phase sediments extending beyond the zone of observable crustal stretching (Figure 8.9). This characteristic stratigraphy may be a function of either the finite flexural rigidity of the lithosphere or by differential stretching of the crust and mantle parts of the lithosphere (White & Mckenzie, 1988).

In the calculations above we made the simplyfying assumption that the lithosphere is everywhere in local isostatic equilibrium. This amounts to treating the lithosphere as infinitely weak (to vertical shearing) and therefore incapable of supporting bending or flexural stresses. The lithosphere is not infinitely weak but rather has a finite flexural rigidity. The consequences of the finite flexural rigidity are that loads on the lithosphere, such as sedimentary or tectonic loads, are generally compensated on a regional scale rather than a local scale, with significant local departures from isostatic equilibrium. In the case of the deposition of sediment in rift basins the finite flexural rigidity of the lithosphere allows compensation of the sedimentary load beyond the zone of stretching allowing downward flexing of the lithosphere beyond the stretched zone. Since the active rift phase is usually associated with a severe thermal perturbation the flexural rigidity of the lithosphere is expected to be very low during active rifting. However cooling attendent with sag pahse subsidence will increase the flexural rigidity causing the loading due to emplacement of sediment to be compensated on length scales which increase with time (Figure 10.3).

Figure 10.10: Schematic cross section showing characteristic steer's head formation of rift basins.

White and McKenzie (1988) argue that the flexural rigidity of continental lithosphere is insufficient to explain the Steer's Head geometry and regard the characteristic geometry as a function of heterogeneous stretching of the lithosphere, with the characteristic length scale for mantle stretching greater than the length scale for crustal stretching. Indeed there is no good reason that stretching should be homogeneous on the scale of the lithsophere, and the resulting values of bc and bl (the reciprocals of fc and fl - fc, respectively) are shown in Figure 10.3. This allows subdivision of the basin into three regions which show characteristic elevation changes. Region a where the rift phase subsidence is greater than for homogeneous stretching; region b where the rift phase elevation change is either subsidence or uplift, but which show substantial sag phase sedimentation with sedimetation onlapping the basin margins progressively through the sag phase, and region c where there is only uplift during rifting and sag restores the initial elevation of the lithosphere (i.e., no net subsidence). The respective fl-fc paths are shown in Figure 10.3.

Figure 10.11: Heterogeneous stretching model for Steer's Head formation (after White and McKenzie, 1988). See text for discussion.

10.4  Topography of normal fault terrains

Rotation of normal fault blocks around horizontal axis is a necessary consequence of extension on a single set of domino faults (Figure ). In the case where the normal faults are planar, rather than listric, it is easy to calculate the finite extension, b, from the dip slope, q, as well as the observed orientation of the faults, a (Figure )
b = li
lo
= sin(q+a)
sin(a)
(10.17)

Figure 10.12: The deformation of planar faults forming a set of dominoes requires rotation about a horizontal axis in order that strain compatibility be maintained. Planar normal faults of this type root in a zone of pervasive aseismic deformation corresponding to the seismic-aseismic transition at around 10-15 km depth.

Figure 10.13: Schematic figure depicting the determination of the stretching factor from orientation of a domino fault system (see text for discussion).

It is possible that the footwall escarpments may be emergent even though the isostatic considerations indicate extension produces overall subsidence. The emergence of the footwall escarpments is due to the horizontal rotation component of tilting normal fault sets. Figure shows schematically how the problem can be conceptualised.The rotation of the dip slopes can be understood as a rotation about a fulcrum. This is the point on the dip slope that in the absence of any isostatically induced subsidence would remain at constant elevation. In fact this fulcrum undergoes the calculated subsidence (Si) for the appropriate extension, b. The condition for an emergent footwall escarpment is clearly that DH > Si. For a given amount of extension (i.e., fixed Si), DH will be proportional to the horizontal spacing of the faults.

Figure 10.14: See text for discussion

Despite the common assertion that normal fault systems are listric (that is they curve to shallower orientations at depth), there is considerable debate amongst seismologists as to whether this is the case, with at least one vociferous group (McKenzie, Jackson etc.) arguing that normal fault sets are generally planar. In modern extensional fault systems normal faults tend to nucleate at around 60o. With increased extension the fault sets rotate to shallower angles and seem to lock at around 30o. Futher extension must be accommodated by the formation of a new set of steeply dipping (60o) faults, which dissect an progressively rotate the older inactive set.


Part 4
Synthesis


Chapter 11
The driving mechanisms reviewed

We have seen that significant forces are associated with the changes in potential energy, accompanying the generation and ageing ocean lithosphere and its subduction; forces which we refer to as ridge push (about 2-3 x 1012 N m-1 as seen by old ocean lithosphere) and slab pull (up to the order of 1013 N m-1). Similarly, potential energy variations associated with continental topography are capable of generating large buoyancy forces. In addition, the passage of the lithosphere over the deeper, hotter convective mantle, or, viceversa, the traction exerted by the convective flow on the base of the lithosphere, may act to either resist or drive plate motion (Fig 11.1). An important constraint on the way these tectonic forces interact is provided by the observation that individual plates are not accelerating. This requires that a basic force balance or, more strictly, a torque balance operates on individual plates. In the following discusion we consider the implications of this reuirement for torque balance for the driving mechanisms a number of the Earth's plates.

Figure 11.1: Mantle flow as a driving mechanism.

11.1  Torque balance and plate velocity.

A force acting on a lithospheric plate produces a torque, T whose magnitude is given by the cross product of the force F and the radius of the Earth (defined by the radius vector,T) :
T = F  x  r
with a torque pole given by the left-hand rule.

Figure 11.2: Torque balance.

The notion that the plates are not appreciably accelerating implies a torque balance, which must reflect the interaction between forces which drive plate motion and those that resist plate motion. It seems reasonable (to me at least) to assume that the torque pole of the force and/or combintaion of forces driving plate motion is correlated with the velocity pole, and thus we should be able to identify the driving forces from the resistive forces, by comparison of torque poles and velocity poles.

There is substantial variation in the absolute velocity of the Earth's major plates which correlates to some degree with the gravitational torques acting on the plate ( see reprint Coblentz et al. On the Potential energy of the Earth's lithosphere)

11.2  The African plate

The African and Antarctic plates are similar in as much as they are both large, slowly moving plates, incorporating significant continental areas bounded amost entirely by passive margins and oceanic ridge systems with only minor lengths of convergent margin boundaries. Moreover, both are characterised by relative slow absolute velocities and by active continental tectonics dominated by rifting (e.g., the African rift system in Africa and the Ross Sea- Trans-Antarctic Mountain rift system).

Traditionally, the rifting that has occurred in the continental regions of both plates has been seen as the consequence of a process actively involving the sub-lithospheric mantle, such as the impingement of a plume on the base of the lithsophere (such as illustrated in Fig. 11.1). It is estimated that the imingement of a plume at the base of the lithdophere can jack the lithsophere upo by as much as 1-2 kms, with a corresponding increase in potential energy of up to 3 x 1012 N m-1. While this is undoubtedly an important process in generating extensional stress regimes in continents, another possibility is simply related to the way in which potential energy is distributed across the plate and how this potential energy distribution evolves in time as the oceans spread outward around the continent. To understand this it is useful idealise the African Plate as a completely circular plate which has grown uniformily with time (Fig.11.3), and imagine how the potential energy of the plate evolves with time. What is growing with time is largely the old ocean lithosphere which represents a potential energy low within the plate. Consequently, the mean plate potential energy must decline with time, and the difference between the continental potential energy and the mean platte potential energy becomes progressively greater producing significant tension in the most elevated parts of the continent.

Figure 11.3: Circular plate analogy to the African and Antarctic Plates.

For those interested more details can be found in the reprints Sandiford & Coblentz :Plate scale potential energy differences and the fragmentationof ageing plates), and Coblentz & Sandiford, Tectonic Stresses in the African Plate: Constraint on the Ambient Stress State.

11.3  Torque balance in the Indo-Australian Plate

The Indo-Australian (IAP), North American and South American plates form a group of fast moving "continental" plates (Minister & Jordan, 1978). In the North American and South American continents the orientation of the maximum horizontal compression (SHmax) is well defined and is clearly aligned with the absolute plate velocity (Richardson, 1992). In contrast, the intraplate stress field within the Australian continent is complex, and thus cannot easily be explained by any single tectonic process. Like North America and South American, but unlike the slower moving plates such as Europe and Africa, the stress field within the interior of the Australian continent is largely compressional. In the northern part of the plate SHmax is aligned N to NNE more or less orthogonal to the collisional boundary in New Guinea. Elsewhere, the orientation of SHmax forms a divergent fan resulting in E-W compression in western Australia, in south-eastern Australia and along the southern margin. While the stress field in the northern and western part of the continent have been relatively successful modelled in terms of plate-scale tectonic processes (Cloetingh & Wortel, 1985,1986; Richardson, 1987; Coblentz et al., 1993a), the sources of the E-W compression in SE Australia remains poorly understood.

In the North and South American plates, the uniform intraplate stress field orientation reflects, in part, the relatively homogeneous boundary configuration of the plates; both having relatively long mid-ocean ridge segments along their trailing (eastern) margins and long continental arc-related mountain tracts along the leading (western) margins. While the IAP is similarly configured with cooling ocean lithosphere dominating the entire southern boundary, the northern and eastern convergent boundaries are heterogeneous consisting variously of continent-continent collisions (Himalaya, New Guinea, New Zealand), continent-arc collisions (Banda Arc), and oceanic-trench segments (Java Trench, Tonga-Kermadec trench). Below we show that the compressive stress pattern in the central and western part of the Australian continent relates to focussing of stresses arising from the plate-scale distribution of gravitational potential energy along the Himalayan and New Guinea collision segments. We then show how similar notions applied to the eastern boundary naturally account for the enigmatic EW compression observed in SE Australia, and, if correct, considerably down-plays the role of subduction and basal traction in the IAP intraplate stress field.

At the outset it is necessary to emphasize that the main constraint we will use here to evaluate the origins of the Australian stress field are the orientations of the in situ stress field. There is now a considerable dataset on stress orientations for the IAP (Zoback,1992). In contrast, our knowledge of stress magnitudes is very poor, and remains the subject of some controversy. This is especially the case for the central Indian Ocean, which is unique in having active deformation of the oceanic lithosphere in the central Indian Ocean (refernces).

The sources of stress that act on plates include: (1) intraplate sources related to lateral variations in gravitational potential energy of the lithosphere (Coblentz et al., 1993b), (2) tractions transmitted across convergent plate margins and (3) tractions at the base of the plate related to its motion over the asthenosphere. The best understood of these are the intraplate variations in potential energy in the ocean lithosphere which give rise to "ridge push" of about 2-3 x 1012 N m-1 of young ocean lithosphere on old (> 80 Ma) ocean lithosphere (note that variations in lithopsheric potential energy can be directly related to geoid anomalies which for ocean ridges are of the order of 10-15 m, Sandwell & Schubert, 1980). Other intraplate sources are associated with continental margins where the difference in potential energy is of the order of 1-2 x 1012 N m-1 (with a corresponding geoid anomaly of about 6 m, Haxby & Turcotte, 1978; Coblentz & Richardson, 1992; Coblentz et al., 1993b) and with areas of high topography developed at continental collisional margins (see discussion below). In comparison with these intraplate stress sources, the magnitudes of stresses imposed at convergent margins and the tractions at the base of the plate are poorly constrained, and estimates very by orders of magnitudes (c.f., Cloetingh & Wortel, 1985 & 1986; and Richardson, 1987). While the negative buoyancy of subducting lithosphere is relatively easy to quantify (of the order of a few times 1013 N m-1 for a fully developed slab), and potentially provides a large "net" tensional force to the trailing plate, the extent to which the stresses arising from the density defect are dissipated in the subduction zone is not known. More controversial are the "trench suction" forces transmitted to the over-riding plate at subduction zones which may apparently be either compressional or tensional. Continent-continent and continent-arc collisions are likely to impose considerable resistance to plate motion because of the buoyancy of continental lithosphere and hence are likely to be a source of intraplate compression. Some idea of the magnitude of the forces associated with collisional processes can be determined by the change of the potential energy associated with the construction of convergent orogens. Precise quantification of these potential energy changes are difficult to establish because of our inadequate knowledge of the deep lithospheric density structure, and because we cannot accurately measure geoid anomalies on the continents. However, the crustal contribution to the excess potential energy of regions of high elevations is proportional to the square of the crustal thickening (e.g., England & McKenzie, 1982), and for the high Himalaya, where crust is approximately double the normal continental thickness may be as much as 5 - 10 x 1012 N m-1 (Zhou & Sandiford, 1992).

Because we have much more confidence in the magnitude of the intraplate sources of stress, than with forces associated with plate boundaries it is useful to model the effect of intraplate sources of stress without applying overly stringent boundary conditions (our modelling strategy is outlined in Appendix 1). We note that the lateral variations in the lithospheric potential energy provides substantial torques in a number of plates (e.g. Richardson, 1992, Coblentz et al., 1993). For plates with the largest gravitational torques (the Pacific, Indo-Australian and South American), there is good correlation between the gravitational torque poles and velocity poles (Coblentz et al., 1993), suggesting gravitational torques may be an important driving mechanism for plate motion. In the IAP the potential energy distribution produces a torque of 6.43 x 1025 N m about a pole at 27o N, 62.4o E which closely matches the velocity pole xxo N, xxo E (the angular misfit is 14o). This torque is mostly reflects potential energy distributions of the ocean ridge systems along the southern and western margin of the plate (which contribute 8.3 x 1025 N m-1 about a pole at 40.7oN, 31.7oE, but also includes contributions from the continental masses.

The assumption that the lithospheric plates are not accelerating (Solomon, 197x; Richardson et al., 1979) and thus are in mechanical equilibrium required balance of the gravitational potential energy torque by some other forces. The close correspondence between gravitational potential energy torque poles and the plate velocity pole requires that the net some of the opposing forces act against in a sense that is resistive to the current plate motion. Such resistance may be achieved either by drag at the base of the plate or by resistance along the convergent margins, and in the following sections we describe the stress fields for each of these scenarios. Figure 2 shows that the stress field produced all resistance provided by basal shear bears little resemblance to the observed stress field. Thus, we consider it is likely that at least some of the torque balance is achieved by resistance imposed along the northern margin and in Figure 3 we show the the stress fields predicted by applying force balance by fixing the whole of the northern boundary and by fixing only those segments involving collision of continental lithosphere of the IAP, that is, Himalaya, New Guinea and New Zealand

The match between the observed and predicted stress fields is far better for Figure 3b than for Figure 3a, especially for the western Australia, suggesting that much of the compressional stress within the IAP is the result of focusing the potential energy torque (mostly arising form ridge push) along collisional boundary segments. The predicted occurrence of EW compression in SE Australia in Figure 3b is of particular importance in that it is in accords with the evidence of a range of in-situ stress indicators (e.g., Zoback, 1992) the origin of which hitherto has been enigmatic. This predicted compression arises from the collisional segment of the plate boundary in New Zealand.

An important aspect of the results indicated in Figure 3b is that the main features of the Australian continental stress field can be reproduced by the interaction of two governing processes, namely gravitational potential energy torques (mainly due to "ridge push") and collisional resistance. Figure 3b clearly shows that the main complexity in the stress field reflects the heterogeneous disposition of the collisional segments along the northern and eastern convergent boundaries of the IAP. If this interpretation of the intraplate stress field is correct it raises the important question of the role of subduction at convergent boundaries in the intraplate stress field, the implication being that subduction processes provide, at best, a second-order control on the IAP stress field. We return to this question following a discussion of the predicted stress magnitudes.

The main conclusions that stem form the work presented here, which may have important global implications are:

Finally, the analysis presented here suggests that the complexity in the Australian stress field, in comparison with other continents such as North America and South America seems simply to reflect the heterogeneous convergent boundary conditions operating on the northern and eastern boundary. While the role of the northern boundary of the IAP has long been suspected, our new interpretation of the origin of the E-W compression in the SE Australia further emphasizes the importance of stress focussing at collisional boundaries.

Also See preprint Sandiford et al. Focussing ridge torques during continental collision in the Indo-Australian plate

Chapter 12
Precambrian geodynamics

Chapter 13
Crustal growth models

13.1  Rare gas constraints

The geochemical/isotopic data discussed above implies that the upper mantle (from the base of the lithosphere down to at least 650 km) has been very efficiently depleted of both its atmophile rare gases and its continental lithophile trace elementswith some 99% of its rare gas content lost within the first 50 Ma of Earth history. Such an interpretation suggests very efficient convection of this early stage may have yielded new crust at 200 times the present rate with the attendant formation of the continental crust from the present MORB source region by the end of the Archaean bringing about its Rb/Sr, U/Pb and Nd/Sm depletion. A component of the OIB source region has not out gassed or undergone significant melt extraction in the Earth's history. This is thought to be a deeper part of the mantle and must either be non-convective, or be part of a separate system unconnected with that of the upper mantle.

Appendix A

A.1  The equations of motion

Consider a small rectangular parallelepiped aligned in a cartesian co-ordinate system, x, y and z, with sides of length dx, dy and dz respectively, as shown in Figure A.1. Remembering that force equals stress by area we begin by decomposing the surfaces forces into stresses acting on each of the faces.

Figure A.1: Co-ordinate system used to derive the equations of motion

Assume that at the point, P, centred within the parallelpiped the normal stresses in the directions x, y and z are given by sxx, syy and szz. The shear stress in the direction normal to x and parallel to y is given by txy and parallel to z by txz. Similarly the shear stresses normal to y are given by tyx and tyz, and normal to z by tzx and tzy. Then the forces acting in the direction x on the faces BB'CC' and AA'DD' are, respectively:



sxx + 1
2
sxx
x
dx

dy dz
and
-

sxx - 1
2
sxx
x
dx

dy dz
(A.1)
where the negative sign is due the fact that stresses are treated as positive in tension and negative in compression. Across the faces A'B'C'D' and ABCD the forces are


txy + 1
2
tyx
y
dy

dx dz
and
-

txy - 1
2
tyx
y
dy

dx dz
(A.2)
Across the faces DCC'D' and ABB'A' the forces are


txz + 1
2
tzx
z
dz

dx dy
and
-

tyz - 1
2
tzx
z
dz

dx dy
(A.3)
The body force acting on the volume with density, r, in the direction x is given by
rX dx dy dz
(A.4)
The total force in the x direction is then


sxx
x
+ tyx
y
+ tzx
z
+ rX

dx dydz
(A.5)
If the component of displacement of point P in the x direction is u then Newton's second law gives
r 2u
2t
= sxx
x
+ tyx
y
+ tzx
z
+ rX
(A.6)
Similar results obtain for the y and z directions. In tectonic settings accelerations can be regarded as neglible, and the only body force is gravity which acts in the vertical direction, taken to be z, and using the convention for summation over repeated indices, the equations of motion can be reduced to:
0 = ti,j
xi
+ ai rg
(A.7)
where ai is the unit vector (0,0,1). Note that in the convention for indices adopted in Eqn A.7 the co-ordinates x, y and z are given by x1, x2 and x3, respectively, while sii = tii. Equation A.7 is general and can be applied to many problems related to tectonic phenomona. However, since it is couched in terms of the components of the stress tensor it must be rendered useful through combination with constitutive equations relating stresses to displacements.

A.2  Calculation of ridge-push force

In order to solve Eqn. 7.7 we need to formulate the density distribution appropriate to Figure 7.2a. The appropriate density distribution is (Figure 7.2b):

rz = rm,
t = 0,
0 < z
rz = rw,
t = t1,
0 < z < w
rz = rm[1 + a (Tm - Tz)],
t = t1,
w < z < w + zl
(A.8)
where rm is the density of mantle at Tm, the temperature of the asthenosphere, rw is the density of water, a is the volumetric coefficient of thermal expansion of peridotite. The density distributions defined by Eqn 7.7 give the following variation (szz)z:
(szz)z = rm g z,
t = 0,
0 < z
(szz)z = rw g z ,
t = t1,
0 < z < w
(szz)z = rw g  w + g 
w+zl

w 
rz dz ,
t = t1,
w < z < w + zl
(A.9)
F1 and F2 are given by:
F1 =
w+zl

0 
(szz)z  dz = rm g (w+zl)2
2
(A.10)
F2 =
w

0 
(szz)z  dz = rw g w2
2
(A.11)
Since in the lithosphere the density is itself a function of depth the third term, F3, in Eqn 7.7 is given by:
F3 =
w+zl

w 
(szz)z  dz
(A.12)
Assuming that the lithospheric geotherm at t1 is linear in depth, the temperature at the surface of the lithosphere Ts = 0oC, and a is independent of temperature then:
Tz = Tm  z
zl
,   rz = rm

1+a Tm

1- z
zl




then

w+zl

w 
rz  dz = g zl rm + g zl rm a Tm
2
= g  zl  rm

1+ a Tm
2


and
F3 = rw  g w zl + g zl2 rm
2


1+ a Tm
2


(A.13)
Thus Eqn 7.7 is given by
FR = rm g (w+zl)2
2
- rw g w2
2
-

rw g w zl+ g zl2 rm
2


1+ a  Tm
2




(A.14)
The condition of isostatic compensation at depth, w +  zl, requires that
szz  (z = w + zl, t = 0) = szz  (z = w +zl,  t = 1)
Solving for zl gives:
zl = w(rm-rw)
a rm Tm


File translated from TEX by TTH, version 2.25.
On 7 Oct 2000, 11:59.