Constraints on Seismic Models from Other Disciplines – Implications for Mantle Dynamics and Composition

1.23 Constraints on Seismic Models from Other Disciplines – Implications for Mantle Dynamics and Composition A. M. Forte, Universite´ du Que´bec a` Montre´al, Montreal, QC, Canada ª 2007 Elsevier B.V. All rights reserved.

1.23.1 1.23.2 1.23.2.1 1.23.2.2 1.23.2.3 1.23.2.3.1 1.23.2.3.2 1.23.2.3.3 1.23.2.3.4 1.23.2.3.5 1.23.2.3.6 1.23.2.3.7 1.23.2.4 1.23.2.5 1.23.3 1.23.3.1 1.23.3.2 1.23.3.3 1.23.3.4 1.23.3.5 1.23.3.6 1.23.4 1.23.4.1 1.23.4.2 1.23.4.3 1.23.4.4 1.23.5 References

Introduction Geodynamic Observables and Mantle Flow Theory Convection-Related Surface Observations Evidence for Mantle Flow in Correlations between Internal Structure and Surface Gravity Anomalies Fluid Mechanical Modeling of Viscous Mantle Flow Governing equations Spectral treatment of the mantle flow equations Internal boundary conditions Boundary conditions at Earth’s solid surface Boundary conditions at CMB Determining viscous flow Green functions Incorporating tectonic plates as a surface boundary condition Geodynamic Response Functions for the Mantle Depth Dependence of Mantle Viscosity Modeling Geodynamic Observables with Seismic Tomography Seismic Heterogeneity Models Mantle Density Anomalies Predicted Tectonic Plate Motions Predicted Free-Air Gravity Anomalies Predicted Dynamic Surface Topography Predicted CMB Topography Tomography-Based Geodynamic Inferences of Compositional Heterogeneity Constraints from Mineral Physics Compositional Density Anomalies Inferred from Joint Shear- and Bulk-Sound Tomography Compositional Density Anomalies Inferred from ‘Hot’ and ‘Cold’ Mantle Heterogeneity Diffuse Mid-Mantle Compositional Horizon Concluding Remarks

1.23.1 Introduction The earliest models of mantle heterogeneity (e.g., Dziewonski et al., 1977; Dziewonski 1984; Clayton and Comer (1983) – reported in Hager and Clayton (1989); Woodhouse and Dziewonski 1984) exhibited a large-scale three-dimensional (3-D) structure which was reasonably well correlated to the major surface manifestations of mantle convection, namely

805 807 807 809 811 811 814 815 817 819 820 821 824 825 829 830 831 834 836 837 843 846 847 849 850 852 853 854

the large-scale nonhydrostatic geoid (Hager et al. 1985) and the large-scale tectonic plate motions (Forte and Peltier, 1987). These early models of mantle heterogeneity thus helped to show that seismic tomography could indeed resolve the lateral variations in mantle structure which drive the convective flow responsible for the ‘drift’ of the continents and the large-scale perturbations in Earth’s gravitational field. 805

806 Implications for Mantle Dynamics and Composition

Advances in seismic tomographic imaging over the past decade have yielded models of global 3-D mantle structure which poses significantly improved resolution and reliability (e.g., Ekstro¨m and Dziewonski, 1998; Ritsema et al., 1999; Me´gnin and Romanowicz, 2000; Masters et al., 2000; Boschi and Dziewonski, 2000; Grand, 2002; Antolik et al., 2003). These improvements are evident in the remarkable accord among the images of large-scale mantle structure provided by these different models. Such agreement is encouraging, especially when we consider that these models are obtained on the basis of different data sets, different theoretical treatments of seismic wave propagation, and different algorithms for inverting the global seismic data. A review of the recent advances in global seismic tomographic imaging may be found in Romanowicz (2003). While the analyses and discussion presented below will focus on how seismically inferred 3-D mantle structure may be used to explain geodynamic data and to model mantle flow dynamics, it is important to note that much effort has also been dedicated to the converse approach, namely to infer 3-D mantle structure from geodynamic data. This inverse approach based on geodynamic constraints began at almost the same time as the earliest global tomography inversions. Hager (1984) proposed a mantle flow model based on estimated density anomalies associated with seismically active subducted slabs and used this model to obtain a good match to very long wavelength geoid anomalies. Although subsequent studies, for example, by Ricard et al. (1989), attempted a more general inference of mantle heterogeneity using different geodynamic data, the subducted-slab geodynamic models developed by Hager (1984) became prevalent. This led, for example, to a study by Forte and Peltier (1989) which constrained the distribution of subducted slab heterogeneity using other geodynamic observables such as tectonic plate motions and core-mantle boundary topography. The latter study also attempted to use the geodynamic observables to infer the buoyancy forces associated with mantle plumes under midocean ridges. Tectonic plate motions and long-wavelength geoid anomalies have proved to be the most important constraints in developing geodynamic models of 3-D mantle heterogeneity in terms of subducted slabs. The geological history of tectonic plate motions as derived from paleomagnetic data (e.g., Gordon and Jurdy, 1986) have been an especially important ingredient in developing models of both past and present-day subducted slab heterogeneity

(e.g., Richards and Engebretson, 1992; Ricard et al., 1993). These models have been very useful for exploring time-dependent dynamics of the mantle (e.g., Lithgow-Bertelloni and Gurnis, 1997; Lithgow-Bertelloni and Richards, 1998). Direct comparisons between the long-wavelength pattern of subducted slab heterogeneity and the corresponding pattern of heterogeneity derived from seismic tomography have shown good correlations at the longest wavelengths (Richards and Engebretson, 1992). When all wavelengths are considered, however, the global correlations between the slab heterogeneity and different tomography models are relatively poor, with correlation coefficients less than 0.3 throughout the mantle (Figure 7 in LithgowBertelloni and Richards (1998)). There are various factors which might explain the less than satisfactory agreement between seismic models of 3-D mantle structure and the reconstructions of mantle heterogeneity in terms of subducted slabs. First, from a purely technical perspective, the evolution of the slab trajectories were not determined in a fluid mechanically consistent manner by numerically solving the advection–diffusion equations. Progress in this direction has been made by Bunge et al. (1998), McNamara and Zhong (2005), and Que´re´ and Forte (2006) by solving the full set of thermal convection equations with moving tectonic plates as a surface boundary condition. Second, the slab models assume that thermally generated heterogeneity dominates in the mantle and this may not be applicable in regions of the mantle with significant compositional heterogeneity (e.g., Forte and Mitrovica, 2001; McNamara and Zhong, 2005). Perhaps the most important deficiency in these efforts to explain mantle heterogeneity in terms of slab subduction alone is that they do not account for the presence and evolution of hot thermal plumes in the mantle which have been consistently imaged in the global tomography models (e.g., Romanowicz, 2003). A recent appraisal of the origin and importance of mantle plumes in convection models which incorporate Cenozoic plate histories may be found in Que´re´ and Forte (2006). Seismic tomography continues to be the single most important technique for directly inferring the 3-D heterogeneity in the mantle which is associated with the process of thermal convection. The significant advances in global seismic imaging over the past few years provides the underlying motivation for the detailed discussion of the geodynamic implications which will be presented below. A recent detailed

Implications for Mantle Dynamics and Composition

review by Becker and Boschi (2002) of the currently available tomography models focused on a quantitative analysis of similarities and differences between these models. The work presented below will further extend this previous analysis by carring out calculations of the mantle flow predicted on the basis of these tomography models and examining in detail the extent to which the tomography-based flow models provide a satisfactory explanation for the main convection-related surface observables.

1.23.2 Geodynamic Observables and Mantle Flow Theory In the following, the main global geophysical constraints on mantle structure will be presented. These constraints include the free-air gravity anomalies, the dynamic surface, and core–mantle boundary (CMB) topography and the tectonic plate motions. The theoretical relationship between these convectionrelated observables and the internal 3-D structure will be developed in terms of a fluid mechanical model of mantle dynamics.

807

1.23.2.1 Convection-Related Surface Observations An understanding of thermal convection in the mantle is necessary for explaining a multitude of geophysical and geological processes which we can observe and measure at the surface of the Earth, such as continental drift, earthquakes, mountain building, volcanism, perturbations in Earth’s gravitational field, variations in oceanic bathymetry and continental elevation, and long-term changes in global sea-level variations. The principal surface manifestations of mantle convection which have been employed to study the large scale structure and dynamics of the solid Earth are illustrated schematically in Figure 1. The observational constraints on 3-D mantle structure and dynamics which are considered below are the global free-air gravity anomalies, the dynamic surface and CMB topography, and the horizontal divergence of the tectonic plate motions. The Earth’s gravitational potential perturbations are usually represented in terms of geoid anomalies however, as pointed out in Forte et al. (1994), a more detailed and evenly balanced representation

Geoid (gravity) anomalies Dynamic surface topography Tecton

late

nic p

Tecto

ons

dulati

vel un

Sea le

Tectonic plate divergence

Convecting

ic plate

mantle

Dynamic CMB topography

Liquid outer core Figure 1 Convection-related surface observations. The most direct manifestations of thermal convection in the mantle (represented schematically by the upwelling plume colored in red) are the surface motions of tectonic plates – which may be summarized by their horizontal divergence field – and the flow-induced topography of the solid surface and the CMB. The boundary undulations and the density anomalies in the mantle (due, for example, to the hot upwelling plume) give rise to surface gravitational potential perturbations which may be measured in terms of geoid or gravity anomalies

808 Implications for Mantle Dynamics and Composition

(a)

EGM96 free-air gravity anomalies (L = 2–32) 90° 120° 150° 180° –150° –120° –90° –60° –30°

(b)

90°

0° 90°

60°

60°

30°

30°

0°

0°

–30°

–30°

–60°

–60°

–90°

–90°

0°

30°

60°

–60

–40

–20

0

20

40

0°

0° 90°

60°

90° 120° 150° 180° –150° –120° –90° –60° –30°

–1

0

1

0°

2

km

(d)

NUVEL-1 plate divergence (L = 1–32) 30°

60°

–2

60

mGal

(c)

CRUST2.0 corrected topography (L = 1–32) 30°

90° 120° 150° 180° –150°–120° –90° –60° –30°

0°

Mantle

60° 30° 0°

Core

–30° –60° –90°

–1.5 –1.0 –0.5

0.0

0.5

1.0

1.5

rad/(10 Ma)

Excess flattening of CMB = 400 m

Figure 2 Geodynamic observables. (a) The free-air gravity anomalies derived from the joint geopotential model EGM96 (Lemoine et al., 1998). (b) The dynamic surface topography obtained by removing the topography due to isostatically compensated crustal heterogeneity from Earth’s observed topography. The CRUST2.0 (Bassin et al., 2000) crustal heterogeneity model is employed here. (c) The horizontal divergence of the tectonic plate velocities given by the NUVEL-1 model (DeMets et al., 1990). (d) The excess or dynamic CMB ellipticity inferred from core nutation data (Mathews et al., 2002). All fields, with the exception of (d), have been expanded in spherical harmonics up to degree and order 32 (see eqn [1]).

(especially at long wavelengths) of the spectral content is provided by the free-air gravity anomalies shown in Figure 2(a). We will employ the term ‘dynamic topography’ to mean all contributions to Earth’s surface topography which arise from density anomalies in the convecting mantle – including the lithosphere. Observational constraints on dynamic topography therefore require an accurate model of crustal heterogeneity in order to remove all isostatic crustal contributions to Earth’s measured surface topography. The topographic crustal correction, here based on model CRUST2.0 (Bassin et al., 2000), is described in detail by Perry et al. (2003) and the resulting dynamic topography is shown in Figure 2(b).

The tectonic plate velocity field v may be conveniently summarized in terms of two complimentary scalar fields: the horizontal divergence rH ? v and the radial vorticity rˆ ? r  v (Forte and Peltier, 1987). The constraint of plate rigidity imposes a linear dependence between these two scalar fields (Forte and Peltier, 1991) and hence it suffices that we consider only the plate divergence, shown in Figure 2(c). (A more detailed discussion of the implications of plate rigidity, in terms of allowed plate motions, is also presented in Section 1.23.2.3.7.) The most robust global constraint on deep-mantle density heterogeneity and dynamics (Forte et al., 1995) is currently provided by the dynamical ellipticity of the CMB shown in Figure 2(d). The

Implications for Mantle Dynamics and Composition

discrepancy between theoretical predictions of the free core nutation (FCN) period (Wahr, 1981) and the value determined from very long baseline interferometry (Herring et al., 1986), led Gwinn et al. (1986) to conclude that the CMB ellipticity is larger than that implied by theoretical calculations which assumed a rotating Earth in hydrostatic equilibrium. This inference of an ‘excess’ ellipticity, in which the nonhydrostatic radius of the CMB at the poles is 400 m less than at the equator, is supported by recent analyses of the FCN period (Mathews et al., 2002). The most appropriate mathematical basis functions for describing any bounded and continuous function on a spherical surface are the spherical harmonics Ym, (, j), where position on the spherical surface is defined by colatitude  and colongitude j. We may therefore expand the geodynamic observables in Figure 2 in terms of spherical harmonics as follows: f ð; jÞ ¼

1 X þ, X

f,m Y,m ð; jÞ

½1

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u þ, uX , ¼ t f,m f,m

809

½3

m¼ – ,

in which f,m are the spherical harmonic coefficients of the surface field (see expression [1]) and  denotes complex conjugation. A spectral description of the spatial correlation between two fields, as a function of harmonic degree or wavelength, is quantified in terms of ‘degree correlation’ r,, defined as follows: Pþ, m m m¼ – , f1, f2, qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi r, ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Pþ, Pþ, m f m m f m f f m¼ – , 1 , 1 , m¼ – , 2 , 2 ,

½4

in which f1 m, and f2 m, are the spherical harmonic coefficients of the two fields. Degree variances and correlations are very useful spectral characterizations of the amplitude and spatial variation of globally defined surface fields (e.g., O’connell, 1971) and they are used frequently in the discussions presented below.

,¼0 m¼ – ,

where the function f (, j) represents any of the observables (gravity, topography, or plate divergence) and the indices ,, m which characterize each spherical harmonic are called the ‘harmonic degree’ and ‘azimuthal order,’ respectively. A useful introduction to spherical harmonic functions may be found in Jackson (1975). The spatial variation of the spherical harmonics is oscillatory in character; it may be characterized in terms of an equivalent horizontal wavelength. On a spherical surface of radius r, a spherical harmonic Ym, (, j) has a characteristic wavelength , given by the following expression: 2r 2r , ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ; valid for ,  1 ,ð, þ 1Þ , þ 1=2

½2

In practice, the infinite sum over degree , in [1] is truncated at some finite value ,max and for the fields in Figure 2 it is ,max ¼ 32. At the Earth’s surface (r ¼ 6368 km) this is equivalent to a minimum horizontal length scale (or half wavelength) of about 600 km. Spherical harmonic representations of surface fields can be used to quantify their spectral content in terms of an amplitude spectrum. The amplitude spectrum measures the mean (globally averaged) amplitude of a field at a particular wavelength , [2] and it is defined in terms of a ‘degree variance’ , as follows:

1.23.2.2 Evidence for Mantle Flow in Correlations between Internal Structure and Surface Gravity Anomalies The need for a dynamical interpretation of the surface observables (Figure 2) in terms of very deepseated convective flow in the mantle is illustrated in Figure 3 in which long-wavelength gravity anomalies are directly compared to seismically inferred lateral heterogeneity in the lower mantle. The large-scale pattern of faster than average shear velocities which lie below the periphery of the Pacific Ocean (blue colored regions in Figure 3(b)) has long been interpreted in terms of accumulating lithospheric slabs which have subducted below the margins of the Pacific Ocean since Mesozoic times (Chase, 1979; Chase and Sprowl, 1983; Richards and Engebretson, 1992). These seismic anomalies, which presumably represent positive density anomalies (i.e., cold subducted slabs), are closely correlated with a broad ring-like pattern of negative gravity anomalies around the periphery of the Pacific Ocean (red colored regions in Figure 3(a)). The correlation between long-wavelength seismic heterogeneity in the lower mantle and surface gravity anomalies was first identified by Dziewonski et al. (1977), and they postulated that a dynamical interpretation in terms of mantle flow could explain the negative sign of the correlation. In a rigid or elastic

810 Implications for Mantle Dynamics and Composition

(a) 0°

30°

60°

90° 120° 150° 180° –150° –120° –90° –60° –30°

0°

90° 60° 30° 0° –30° –60° –90°

–30

–20

–10

(b) 0°

30°

60°

0 10 mGal

20

30

90° 120° 150° 180° –150° –120° –90° –60° –30°

0°

90° 60° 30° 0° –30° –60° –90°

–1.0

–0.5

0.0 dVs /Vs (%)

0.5

1.0

Figure 3 Surface gravity anomalies and deep-mantle heterogeneity. (a) Long-wavelength free-air gravity anomalies derived from the joint geopotential model EGM96 (Lemoine et al., 1998). (b) Long-wavelength seismic shear velocity anomalies at a depth of 2100 km derived from the tomography model of Grand (2002). Both the gravity and seismic anomaly fields are synthesized from spherical harmonics in the degree range , ¼ 18.

mantle in which flow is absent, a positive correlation between internal density and surface gravity is expected. However, in a convecting mantle, the gravitational signal of internal density anomalies is opposed by the effect of flow-induced topography at the surface and at the CMB, and this dynamical balance can lead to the observed negative correlation (e.g., Pekeris, 1935; Parsons and Daly, 1983; Richards and Hager, 1984). It is instructive to consider a more detailed spectral analysis of the correlation between surface

gravity anomalies and seismically inferred lateral heterogeneity at different depths in the mantle. The degree correlations shown in Figure 4 indicate which depth ranges in the mantle may be important contributors to the different wavelengths of the surface gravity field. The strongest positive correlations appear in the transition-zone region of the mantle (from 400 to 800 km depth) and at the longest wavelength (degree 2) they extend down to about 1400 km depth. The positive sign of these correlations indicates that long-wavelength gravity highs

Implications for Mantle Dynamics and Composition

1.0

50 km

0.0 –1.0 1.0

95% 80% 80% 95%

200 km

0.0

Degree correlation

0.0 –1.0 1.0

1.0

0.0

0.0 –1.0 600 km

1.0

0.0

0.0

–1.0 1.0

–1.0 1.0

800 km

0.0

2000 km

2200 km

–1.0 1000 km

1.0

0.0

0.0

–1.0

–1.0

1.0

1800 km

0.0

–1.0 1.0

1600 km

–1.0 400 km

–1.0 1.0

1400 km

0.0

–1.0 1.0

1.0

1200 km

1.0

0.0

0.0

–1.0

–1.0 0 2 4 6 8 Harmonic degree

2500 km

2800 km

0 2 4 6 8 Harmonic degree

Figure 4 Depth-dependent correlations between 3-D mantle structure and gravity anomalies. Each frame shows the degree correlations between surface gravity anomalies and lateral variations in seismic shear velocity at a particular depth in the mantle. The seismic anomalies are derived from an earlier long-wavelength tomography model SH8/U4L8 (Forte et al., 1993a). Adapted from Forte AM, Dziewonski AM, and Woodward RL (1993a) Aspherical structure of the mantle, tectonic plate motions, nonhydrostatic geoid, and topography of the core-mantle boundary. In: Le Moue¨l J-L, Smylie DE, and Herring T (eds.) Dynamics of the Earth’s Deep Interior and Earth Rotation pp. 135–166. Washington DC: American Geophysical Union.

The correlation between mantle heterogeneity and surface gravity anomalies in a convecting mantle has been shown to be a strong function of the depth dependence of the rheological properties of the mantle (e.g., Hager, 1984; Ricard et al., 1984; Hager et al., 1985; Forte and Peltier, 1987). In particular, as will be discussed further below, the change from positive correlations at the top of the lower mantle to negative correlations at the bottom of lower mantle (Figure 4) can be interpreted in terms of a significant increase in mean viscosity across the mantle. Other surface observables, such as the tectonic plate motions, also provide important constraints on the depth dependence of mantle rheology (e.g., Ricard and Vigny, 1989; Forte et al., 1991). The geodynamic surface observables also provide direct constraints on the 3-D distribution of density anomalies in the mantle and hence these observables provide a fundamentally important and independent means of evaluating the extent to which seismic tomography models successfully resolve the lateral heterogeneity in the mantle which is associated with mantle convection. In the next section, the dynamical link between lateral heterogeneity in the mantle and convectionrelated surface observables is developed in the theoretical framework of a fluid mechanical model of viscous flow in the mantle. The reader wishing to skip over this detailed mathematical treatment of the mantle flow theory is encouraged to jump directly to Sections 1.23.2.4 and 1.23.2.5, where the viscous response of the mantle and current mantle viscosity inferences, respectively, are summarized.

1.23.2.3 Fluid Mechanical Modeling of Viscous Mantle Flow 1.23.2.3.1

overlie similarly large-scale seismically fast regions located at the base of the upper mantle. Negative degree correlations between long-wavelength seismic and gravity anomalies are found in the nearsurface region, down to about 400 km depth, and in lower mantle below 1400 km depth (as in Figure 3) where the strongest negative correlations are found. These depth-dependent degree correlations, obtained on the basis of long-wavelength tomography models derived more than a decade ago, are robust observations which have been verified using the latest seismic tomographic inferences of 3-D mantle structure (Ricard et al., 2006).

811

Governing equations Laboratory and geological evidence of the ability of mantle rocks to creep indefinitely over geological timescales is understood in terms of the existence of atomic-scale defects in the lattice of crystal grains (e.g., Nicolas and Poirier, 1976). If the ambient mantle temperature is sufficiently high, the imposition of nonhydrostatic stresses causes the lattice defects to propagate and thus allows the mantle rocks to creep or ‘flow’ slowly. This process may be characterized in terms of a single parameter, namely an effective viscosity (e.g., Weertman and Weertman, 1975; Weertman, 1978). The characterization of the long-term creep properties of the mantle in terms of an effective viscosity

812

Implications for Mantle Dynamics and Composition

allows us to model the slow flow of the mantle with the conventional hydrodynamic field equations. The hydrodynamic field equations (Landau and Lifshitz, 1959) which express the principles of conservation of mass and momentum are q þ  ? ðuÞ ¼ 0 qt

½5

du ¼  ? s þ g  dt

½6

in which u is the velocity field, d/dt ¼ q/qt þ u ?  is the total (Lagrangian) time derivative, s is the stress tensor, g is the gravitational acceleration, and  is the density. We may represent the gravity field g in terms of a gravitational potential: g ¼ 

ð X  ? Þ

  2 ij ¼ ui; j þ uj ; i – dij uk; k 3

•

½8

An explicit expression for the stress tensor is given by ij ¼ – Pdij þ ij   2 with ij ¼ ui; j þ uj ; i – dij uk; k þ ij uk; k 3

½9 ½10

where  ij is the viscous stress tensor, dij is the identity tensor, P is the total pressure, and  are the ‘isotropic’ viscosity coefficients, and ui,j ¼ qui/qxj represents the derivative of the velocity components i with respect to the coordinate direction j. A number of simplifications are possible when applying the above field equations to the problem of flow in the mantle: mantle rocks creep much slower than the • Since acoustic velocity in the mantle, we can safely ignore the term q/qt in the conservation of mass equation [5] and we thus have the following anelastic-liquid approximation:  ? ðuÞ ¼ 0

•

½11

This simplification will explicitly rule out acoustic waves as solutions of the flow equations. The second viscosity coefficient  in the viscous stress tensor  ij describes the dissipation associated

½12

will further assume a Newtonian (i.e., linear) • We rheology, in which the mantle viscosity is not a

½7

Notice that the sign convention adopted here is opposite to that generally adopted in classical physics, where a ‘negative’ gradient of the potential is used. With this sign convention, Poisson’s equation for the relationship between the gravity potential and density is  ¼ – 4G

with change in fluid volume (density). This volume dissipation may be neglected if the changes in fluid volume occur on timescales which are much longer than for molecular relaxation processes (Landau and Lifshitz, 1959), and this is certainly true for mantle flow. Therefore, the viscous stress tensor  ij will be purely deviatoric:

•

function of stress or strain rate. This assumption is not necessarily appropriate for the mantle (particularly in high stress regions, such as subduction zones) but it will greatly simplify the mathematical resolution of the flow equations which is presented below. A major simplification, from a purely mathematical standpoint, derives from the assumption that mantle viscosity varies only with radius. Although this assumption will significantly simplify the mathematical treatment of the flow theory, it may have important physical implications for the dynamics of buoyancy induced flow in the mantle (e.g., Richards and Hager, 1989). These implications have been examined previously in the context of global-scale flow in 3-D spherical geometry (e.g., Ricard et al., 1988; Martinec et al., 1993; Zhang and Christensen, 1993; Forte and Peltier, ˇ adek and Fleitout, 2003; Moucha et al., 1994; C 2006) and they will be discussed below. A fundamental physical simplification arises from the very large value of the effective viscosity in the mantle. We may nondimensionalize the conservation of momentum equation with the following variable transformations: ðx; y; zÞ ¼ ðdx; dy; dzÞ   t ¼ d 2 = o t

½13

 ¼ o ; g ¼ go g; ¼ o

where the original variables are on the left and the nondimensional ones are on the right. Quantities d, o, o, o, go are characteristic scales for length, thermal diffusivity, density, viscosity, and gravitational acceleration, respectively, in the mantle. The use, in [13], of a thermal diffusion timescale t ¼ d 2 = o

½14

Implications for Mantle Dynamics and Composition

is appropriate since we are dealing with mantle flow arising from thermal convection. The stress tensor and pressure have the same physical units (see expression [9]), and hence their nondimensional transformation is: P ¼ ðo go d ÞP

and

t ¼ ðo go d Þt

½15

where, again, the original quantities are on the left and the nondimensional ones on the right. By virtue of the constitutive relation [12], and expression [15], we obtain the following nondimensional transformation for the flow variable: u¼

o go d 2 u o

½16

Equation [20] shows that in the absence of inertia, there must at all times be a balance between the buoyancy forces g and the forces of viscous disspation described by  ? t. In other words, any changes in internal buoyancy forces will instantly produce changes in fluid flow: this is a consequence of the essentially instantaneous diffusion of momentum in the mantle. We may define an idealized hydrostatic reference state for the mantle, which corresponds to the absence of any internal flow or deformation (i.e., u ¼ 0). In this situation, the deviatoric stress field t vanishes and the momentum conservation equation [6] reduces to: – Po þ o o ¼ 0

or equivalently, u ¼ go t u

½17

2

where t ¼ d / o is the timescale for viscous diffusion of momentum and  o ¼ o/o is the kinematic viscosity. Finally, substitution of expressions [9] and [13]–[17] into [6] yields the following nondimensional momentum conservation equation:

1 qu go  2 t u ? u ¼  ? t – P þ g  þ d  Pr qt

½18

To assess the importance of the terms on the left-hand side of eqn [18], we assume the following values for the scaling variables: d ¼ 3  106 m (for whole-mantle flow), o ¼ 3.3  103 kg m3, go ¼ 10 m s2, o ¼ 1.5  106 m2 s1, o ¼ 1021 Pa s, and we thus obtain: t ¼ 3  10 – 5 s; t ¼ 6  1018 s; go ¼ 3:3  10 – 6 s – 2 ; Pr ¼ 2  1023 d

½19

The key quantity here is the characteristic time for momentum diffusion t which is very small (i.e., mantle flow will come to a complete halt in much less than a millisecond if all buoyancy forces are suddenly removed). The vanishingly small momentum diffusion time, and hence the very large Prandtl number, implies that the inertial forces (du/dt) in the momentum conservation equation are completely negligible. Therefore, using expression [7], the momentum conservation equation [6] simplifies to:  ? t – P þ  ¼ 0

½20

½21

in which Po, o, and o are the pressure, density, and gravity potentials in the hydrostatic state. Poisson’s equation [8] for a hydrostatic planet is: o ¼ – 4Go

½22

We assume that in a dynamic mantle, with a nonvanishing mantle flow u, the pressure, density, and gravity potentials will be perturbed as follows: P ¼ Po þ P1

in which Pr ¼ t /t is the Prandtl number, which characterizes the ratio of temperature and momentum diffusion timescales.

813

 ¼ o þ 1

 ¼ o þ  1

½23

in which all perturbations are assumed to be small, that is,    P1    1 P  o

   1    1   o

   1    1   o

If we now substitute the perturbed variables [23] into the equations of mass and momentum conservation (11–20) and Poisson’s equation [8], and then substract out the hydrostatic reference equations [21] and [22], we finally obtain the following set of first-order accurate, perturbed equations for mantle flow dynamics: mass conservation

 ? ð  o uÞ ¼ 0

momentum conservation  ?  – P1 þ o 1 þ 1 o ¼ 0 gravity 1 ¼ – 4G1

½24

½25 ½26

Notice in eqn [25], that in addition to the driving buoyancy forces (1o), there also exist self-gravitational loads (o1) due to the perturbed gravity field. The above equations must be supplemented by the linear relationship [12] between stress and strain rate, which is valid for an isotropic rheology:

814

Implications for Mantle Dynamics and Composition

Newtonian constitutive equation   ! 2 t ¼ ru þ ur – Ir ? u 3

½27

coordinate system defined by the following complex basis vectors in spherical geometry:  1  eˆ – ¼ pffiffiffi Jˆ – jˆ 2 eˆ0 ¼ rˆ  1  eˆþ ¼ – pffiffiffi Jˆ þ jˆ 2

1.23.2.3.2 Spectral treatment of the mantle flow equations

A classical technique for solving the dynamical flow equations in 3-D spherical geometry is the spectral Green function method, in which all flow variables are expressed in terms of spherical harmonic basis functions introduced in eqn [1] (e.g., Hager and O’Connell, 1981; Richards and Hager, 1984; Ricard et al., 1984; Forte and Peltier, 1987, 1991). The spectral Green functions provide a very convenient mathematical description of the instantaneous relationship between the mantle density anomalies and the viscous flow field, as well as all surface manifestations of the internal flow dynamics. The spectral method is employed here for solving the coupled equations of mass, momentum, and gravity conservation (eqns [24]–[26]), supplemented by the viscous constitutive relation [27]. The method follows that initially developed by Forte and Peltier (1991) for gravitationally consistent compressible flow in spherical geometry. Other treatments of compressible mantle flow in spherical geometry have been presented by Corrieu et al. (1995), Panasyuk et al. (1996), and Defraigne (1997). One begins by rewriting eqns [24]–[27], in the following Cartesian tensor form:

½29

pffiffiffiffiffiffiffi in which l ¼ – 1 and rˆ, #ˆ, jˆ are the unit basis vectors for the standard spherical polar coordinate system. Following Phinney and Burridge (1973), all the tensors appearing in the original system [28] are rotated into the coordinate system defined by [29], thereby yielding the following covariant tensor form of the dynamical equations: u ;  e  ¼ –

 

_o 0 u o

 ;  e þ o ð1 Þ; – 1 go d 0 ¼ 0 ½30  ;   2  d;   

 ; ¼ – P1 e þ u þu – u ed e 3 ð1 Þ;  e  ¼ – 4G1

in which the Greek indices denote the coordinate directions in system [29] and therefore range over the values (1, 0, þ1). The quantities e , e , and d  are the contravariant, covariant, and mixed tensor representations of the Cartesian identity tensor dij. The velocity (u ) and stress ( ) components are expanded in terms of the generalized spherical functions YNm , (Phinney and Burridge, 1973) such that: u – ðr ; ; Þ ¼

X

U,– m ðr ÞY,– 1m ð; Þ

,; m

uk; k ¼ –

u0 ðr ; ; Þ ¼

_o ur o

ij ; j þ o ð1 Þ; i – 1 go rˆ ¼ 0   2 ij ¼ – P1 dij þ ui; j þ uj ; i – dij uk; k 3

X

U,0m ðr ÞY,0m ð; Þ

,; m

uþ ðr ; ; Þ ¼

X

½28

U,þm ðr ÞY,þ1m ð; Þ

,; m

  ðr ; ; Þ ¼

X

ð þÞm

T, m ðr ÞY,

ð; Þ ½32

,; m

ð1 Þ; kk ¼ – 4G1

in which _ o ¼ do/dr, ur ¼ rˆ? u, and (o), i ¼ gorˆ. It should also be noted that in these equations the total stress tensor  is used, rather than the deviatoric stress  used in [25]. The determination of a solution to the system of tensor equations [28] in spherical geometry may be greatly simplified by using an elegant mathematical technique described by Phinney and Burridge (1973). Following their technique, one introduces a new

½31

ð ;  ¼ – 1; 0; þ1Þ

All scalar fields involved in the governing equations [30] are expanded in terms of ordinary spherical harmonics. For example, 1 ðr ; ; Þ ¼

X ð1 Þm, ðr ÞY,m ð; Þ ,; m

X 1 ðr ; ; Þ ¼ ð1 Þm, ðr ÞY,m ð; Þ ,; m

where Y,m XY,0m.

½33

Implications for Mantle Dynamics and Composition

We can simplify subsequent numerical computations by nondimensionalizing all relevant physical variables using the following transformations:

d ¼ 2888 km X radial thickness of mantle go ¼ 9:82 ms – 2 X mean surface gravitational acceleration –3  ¼ 0:1Mg m X characteristic subducted slab density anomaly o ¼ 1021 Pa s X Haskell ð1935Þ reference value

r ¼ dr ; go ðr Þ ¼ go gðr Þ; ðr Þ ¼ o ðr Þ; 1 ¼ ðÞ1     go d 2 U ; T  ¼ go d T  ; U ¼ ½34 o   4GRo d 1 1 ¼ 2, þ 1

Ro ¼ 6371 km X mean surface radius of Earth

Substituting expressions [31]–[33] into the flow equations [30], and using the nondimensionalization scheme [34] as well as the orthogonality properties of the generalized spherical harmonics (as in Phinney and Burridge 1973), yields two independent systems of flow equations. The first system, governing a poloidal geometry of flow, consists of the following six coupled, first-order, ordinary differential equations:

in which the original variables are on the left of each equation and the nondimensional variables are on the right. The scaling quantities we have used are defined as follows:

0 B B B U B B PC B BU C B C B B B 0C B B B d BT C C B C¼B B dr B T P C B C B B C B B – B 1 C B A B @ B B B g1 B @ 0

0

1

1 0 C C0 1 0 C U0 1 0 o C C B 0 0 CB 0 CB P C C B 0 CB U C C B C B CB C C CB 0 C B 4½3 þ r ð_ o =o Þ 61 1 3 _o C gðr Þ T C C 1 ðr Þ B 0 – 0 C B CB C r 2 o r 2, þ 1  r2 o þ C B CB P C B C  0 2

T C C B C B CB C 41 ½3 þ r ð_ o =o Þ 2 2 þ 321 21 3 C B CB C B 0 – – 0 0 C C C  B 2 2 1 A @ r r o o r r C@ A d C ð – 2, þ 1Þ g1 0 0 0 0 0 1 C Ro C ,ð, þ 1Þ 2 A 0 0 0 0 – r2 r ½35 –

½2 þ r ð_o =o Þ r 21 – r

1 r 1 r

0

in which _ o ¼ do/dr, 1 ¼ [,(, þ 1)/2]1/2, 2 ¼ [(,  1)(, þ 2)/2]1/2, and  is Earth’s mean density (5.5143  103 kg m3). The second system, governing a toroidal geometry of flow, consists of the following two coupled, first-order, ordinary differential equations: d dr

815

0 1 1 1 ! ! UT B r C UT C ¼B A @ 2 22 3 TT TT – r2 r

0

0

variables on the harmonic coefficients of flow [31], stress [32], and gravitational potential [33] are U P ðr Þ ¼ U þ ðr Þ þ U – ðr Þ U T ðr Þ ¼ U þ ðr Þ – U – ðr Þ T 0 ¼ T 00 ðr Þ þ

½37

T P ¼ T 0þ ðr Þ þ T 0 – ðr Þ

½36

T T ¼ T 0þ ðr Þ – T 0 – ðr Þ

1.23.2.3.3

The flow, stress, and gravity variables in the poloidal and toroidal flow systems [35] and [36] are all dependent on the same harmonic degree , and order m (this dependence has been dropped for notational convenience). The dependence of these

3 0 ðr Þ  1 ðr Þ 2, þ 1 

Internal boundary conditions The determination of a unique solution of the system of flow equations [35] and [36] derived in the preceding section requires the specification of appropriate boundary conditions at the top and bottom surfaces of the mantle and internal matching

816 Implications for Mantle Dynamics and Composition

conditions at all material interfaces within the mantle. We consider the latter first since the surface boundary conditions can be obtained as a special case of the internal matching conditions. The need for internal boundary conditions in solving for mantle flow is evident upon inspection of the depth variation of mantle density o(r) (see Figure 5) given by the Preliminary Reference Earth Model (PREM) seismic reference model (Dziewonski and Anderson, 1981). The PREM density profile is characterized by two major jumps at depths of 400 km and 670 km. These density jumps have long been interpreted as manifestations of phase-change boundaries (e.g., Jeanloz and Thompson, 1983) which presumably also affect the depth variation of viscosity at these depths (e.g., Sammis et al., 1977). We now derive the matching conditions which must be satisfied when characteristic properties of the mantle, in particular density, viscosity, and chemical composition, change very rapidly across phasechange or chemical horizons in the mantle. We will approximate such rapid changes as mathematical discontinuities. Denoting the mean (i.e., horizontally averaged) location of the internal material boundary as r ¼ ri, geographic variations in the radial location (i.e., deflections) of the boundary are then represented by:   dr i  r ¼ r i þ dr i ; where we assume   1 ri

½38

A discontinuous change in density across an internal boundary implies that       do ¼ o riþ – o ri – dðr – ri Þ dr r ¼r i

½39

By virtue of this last expression, the density perturbation due to the boundary deflection dri is given by the following sheet-mass anomaly:   do dri dr r ¼ri  þ   ¼ – o ri – o ri – dri dðr – ri Þ

di ¼ –

This last expression shows that an internal boundary deflection gives rise to a mantle buoyancy source which may be approximated by a delta-function load, and the corresponding density perturbation is, in non-dimensional terms,  þ   o ri – o ri – dr i  1 ðr Þ ¼ – dðr – r i Þ  d 

        o riþ U 0 riþ ¼ o ri – U 0 ri –

½42

o(r i )

where and are the mantle density immediately above and below the material boundary, respectively. If the internal boundary corresponds to a chemical discontinuity, then we must impose a zero radial velocity condition:

6 Density (g cm–3)

½41

This expression shows that from a purely fluid mechanical perspective, the dynamical effect of a deformed phase-change boundary is indistinguishable from integrated buoyancy forces elsewhere in the mantle. If the material interface instead corresponds to a chemical discontinuity in the mantle, the fluid-mechanical buoyancy effect is also the same. The main distinction between the deflectioninduced buoyancy of phase-change and chemical discontinuities is that the former is controlled by thermodynamics (i.e., the Clapeyron slope) whereas the the latter is due to the change in intrinsic chemical buoyancy across the interface. The matching condition for the radial component of the flow field is obtained by integrating the first þ row in system [35] from r i ¼ ri   to ri ¼ ri þ , taking the limit  ! 0. In terms of the radial flow variable U0, this yields the following mass conservation equation: o(rþ i )

    U 0 riþ ¼ 0 ¼ U 0 ri –

5

4

3

½40

0

500

1500 1000 2000 Depth (km)

2500

Figure 5 Depth-dependent mean density in the mantle from the PREM seismic reference model (Dziewonski and Anderson, 1981).

½43

The matching condition for the tangential mantle flow velocity is obtained by integrating the second equation in [35], which yields the following continuity of the tangential flow variable UP:     U P riþ ¼ U P riþ

½44

This continuity of tangential flow also applies to the case of a internal chemical boundary.

Implications for Mantle Dynamics and Composition

The internal matching condition for the radial stress T0 across a deformed material interface is obtained by integrating the third row in system [35], in which we substitute result [41] into the buoyancy force term 1/ in expression [35] thereby yielding:     T 0 riþ – T 0 ri –

 

 þ o ri – o ri – 3 ¼ u ðri Þ þ 1 ðr i Þ  2, þ 1    

o riþ þ o ri – dr i – gð r i Þ d 

The internal matching conditions for the toroidal components of mantle flow are readily obtained following the same procedure employed for the poloidal flow. Integration of the first row in system [36] yields the following continuity of toroidal flow across any (phase-change or chemical) interface:     U T riþ ¼ U T ri –

½45

    T T riþ ¼ T T ri –

       2 riþ þ ri– 0  þ  ½46 U ri – U 0 ri – u ðri Þ – o ri

    T P riþ – T P ri – ¼ – 1 u ðri Þ

½47

The term u gives rise to an apparent discontinuity in tangential stress and this is again a consequence of our mathematical treatment of a discontinuous change in radial flow across a material interface. The significance of this apparent jump in tangential stress was first noted by Corrieu et al. (1995), and a subsequent analysis by Forte and Woodward (1997a) determined that effect of this discontinuity was negligible. In the case of a chemical discontinuity, the term u vanishes and the tangential stress is then continuous across the interface. Integration of the fifth row in [35] across the material interface located at mean radius ri yields the following matching condition for the perturbed gravitational potential:     1 riþ ¼ 1 ri–

½48

The matching condition for the perturbed gravitational acceleration is obtained by substituting expression [41] into the buoyancy force term 1/ in [35] and then integrating the sixth row of this system, yielding the following nondimensional expression:     g1 riþ – g1 ri –      þ o ri – o ri – dri d ¼ ð2, þ 1Þ d  Ro

in which g1 ¼ d1/dr.

½51

1.23.2.3.4 Boundary conditions at Earth’s solid surface

It will be assumed that the spherically symmetric, hydrostatic reference Earth model is overlain by a global ocean layer which is 3 km thick, as in PREM (Dziewonski and Anderson, 1981). PREM’s crust and seismic lithosphere (LID) are combined into a single mechanical layer containing the same total mass as the two PREM layers assuming, for simplicity, that both layers deform and move together in response the buoyancy driven flow in the mantle. The top surface of the combined crust–lithosphere is located at radius r ¼ 6368 km and the base is located at radius 6291 km (i.e., at a depth of 80 km below the surface of the global ocean layer). This redefined lithosphere has mass density of 3.2 Mg m3. In view of the very complicated mechanical and rheological properties of the crust and underlying lithosphere, in particular the tectonic plates, it is clear that the flow theory developed here with a purely depth-dependent viscosity cannot provide an adequate representation of the near-surface dynamics. An approximate treatment of the effect of surface tectonic plates is presented in Section 1.23.3.7. In the meantime, we consider here two different boundary conditions of relevance at the solid surface: ‘free-slip’ and ‘no-slip’. Earth’s bounding surface at r ¼ a ¼ 6368 km is a chemical (compositional) boundary across which there can be no flow and hence the matching condition (43) is applicable: U 0 ða – Þ ¼ 0; 

½49

½50

Similarly, an integration of the second row in [36] yields the following continuity of tangential stress across any internal material interface:

in which the term

arises from the discontinuous change in radial flow (see expression [42]) due to a density jump across an internal (phase-change) boundary. This term vanishes in the case of a chemical discontinuity since expression [43] then applies. The matching condition for the tangential stress TP is obtained by integrating the fourth equation in [35]:

817

for both free-slip and no-slip

½52

where a denotes the radial location r ¼ a  . We can, to first-order accuracy, ignore surface deflections (i.e., topography) in specifying the boundary conditions on the flow.

818

Implications for Mantle Dynamics and Composition

The condition for the surface tangential flow UP(a) in the lithosphere is ( P

–

U ða Þ ¼

U P ða – Þ to be determined; for free-slip 0

for no-slip

½53

In the case of the surface toroidal flow, the boundary condition is U T ða – Þ ¼ V T

½54

in which VT is determined from the coupling of surface poloidal and toroidal flows due to the presence of rigid surface tectonic plates. The mathematical formulation of this coupling and an explicit expression for VT is presented in Section 1.23.3.7. We can apply the radial stress T 0 matching condition in [45] to the deformed surface boundary which, by virtue of conditions [52], becomes

3 a 1 ða – Þ T 0 ða – Þ – T 0 ða þ Þ ¼ 2, þ 1 

a da – gðaÞ  d

½55

where we have defined the density jump across the solid surface: –

þ

a ¼ o ða Þ – o ða Þ ¼ 3:2 – 1:0 Mg m ¼ 2:2 Mg m – 3

throughout the ocean layer

throughout the ocean layer

½58





3 a a da 1 ða – Þ – gðaÞ  d 2, þ 1  valid for free-slip and no-slip

for free-slip

T P ða – Þ to be determined; for no-slip

[59]

½60

and the condition for the surface toroidal tangential stress is T T ða – Þ ¼ T T ða – Þ to be determined

½61

On the basis of the general result [48], and using result [49], the surface matching conditions for the perturbed gravitational potential and acceleration are: 1 ða – Þ ¼ 1 ðaþ Þ  

d a da g1 ða – Þ ¼ g1 ðaþ Þ þ ð2, þ 1Þ Ro  d

½62 ½63

The ocean-layer potential and gravity fields 1(aþ) and g1(aþ), respectively, are not independent of each other and, as shown in Forte and Peltier (1991), they are both related to the perturbed potential at the surface of the global ocean layer, 1(r ¼ Ro): 1 ðaþ Þ ¼ P, 1 ðRo Þ þ

g1 ða Þ ¼ g, 1 ðRo Þ

½64 ½65

in which the ocean-layer response functions P, and G, are as follows: "   ,þ1  , – 1 # Ro 3 ! Ro ,þ2 a ½66 P, ¼ – – a a 2, þ 1  Ro   d Ro ,þ2 Ro a "  ,þ3  , – 2 # 3 ! d Ro a þ ð, þ 1Þ þ, a 2, þ 1  Ro Ro

G, ¼ – ð, þ 1Þ

½57

Substituting result [58] into expression [55] yields the desired radial-stress-boundary condition at the surface: T 0 ða – Þ ¼

0

½56

and therefore, by virtue of this result, the fourth row in system [35] yields T 0 ðr Þ ¼ 0;

( T P ða – Þ ¼

–3

and where we have also invoked the universally valid condition [48] for the vertical continuity of the perturbed gravitational potential. In expression [55], T 0(aþ) corresponds to the radial stress in the global ocean layer. If we assume that the viscosity in the global ocean layer is negligible (i.e., / o ! 0), then the second row in system [35] yields T P ðr Þ ¼ 0;

in which da is the flow-induced vertical deflection of the solid surface (i.e., dynamic surface topography). The condition for the surface poloidal tangential stress T P(a) is as follows:

½67

in which w X (aþ) ¼ 1Mg m3 is the density of the global ocean layer. Substitution of results [64] and [65] into expressions [62] and [63] yields the complete surface boundary conditions for the gravitational variables: 1 ða – Þ ¼ P, 1 ðRo Þ;

valid for free-slip and no-slip ½68  

d a da g1 ða – Þ ¼G, 1 ðRo Þ þ ð2, þ 1Þ ; Ro  d ½69 valid for free-slip and no-slip

Implications for Mantle Dynamics and Composition

It should be noted that expression [68] should also be substituted into the radial stress condition [59]. The complete set of free-slip and no-slip surface boundary conditions, in terms of the poloidal-flow vector v(r) ¼ [U 0(r), UP(r), T 0(r), T P(r), 1(r), g1(r)]Tr (where Tr denotes transposition) employed in system [35], are: Free-slip. vða – Þ ¼ U P ða – Þy1 þ 1 ðRo Þy2 þ

a da y  d 3

½70

in which the surface basis vectors are 0 1 0 1 0 0 B C B C B C 0 B1C B C B C B 3 a C B C B B0C P, C B C B C 2, þ 1  C y1 ¼ B C y2 ¼ B B C B0C B C B C 0 B C B C B C B0C B C @ A P, @ A 0 G 0 1 , 0 B C 0 B C B C B C – g ða Þ B C B C C y3 ¼ B B C 0 B C B C B C 0 B C   @ d A ð2, þ 1Þ Ro

½71

a da y  d 3

½72

in which the surface basis vectors are 0 1 0 B C B0C B C B C B0C B C y91 ¼ B C and y2 ; y3 are defined in ½71 B1C B C B C B0C @ A

½73

The surface boundary conditions [54] and [61] for the toroidal flow vector u(r) ¼ [UT(r), TT(r)]Tr which is governed by system [36] are ½74

in which the two surface basis vectors are z1 ¼

1 0

! z2 ¼

0 1

Boundary conditions at CMB The derivation of the boundary conditions which apply at the CMB, located at mean radius r ¼ b ¼ 3480 km, is almost identical to the derivation for the surface boundary conditions in the preceding section. The liquid outer core is regarded as having negligible viscosity relative to the mantle and hence the CMB is treated as a purely free-slip boundary. The only difference concerns the application of gravitational matching conditions at r ¼ b, since we must now deal with the interaction between a deformed CMB and a compressible, hydrostatic core. A detailed treatment of the gravitational perturbations maintained in a hydrostatic core is presented in Forte and Peltier (1991), where it is shown that perturbed gravitational acceleration at the top of the core (i.e., immediately below the CMB) is determined by the perturbed potential at the bottom of the mantle (i.e., immediately above the CMB), as follows: g1 ðb – Þ ¼ R ,  1 ðb þ Þ

vðb þ Þ ¼ U P ðb þ Þx1 þ 1 ðb þ Þx2 þ

b db x3  d

½76

in which the CMB basis vectors are:

0

uða – Þ ¼ V T z1 þ T T ða – Þz2

1.23.2.3.5

where b denotes the radial location r ¼ b   (i.e., bottom side of the CMB) and bþ denotes r ¼ b þ  (i.e., top side of the CMB). R, is a numerically determined coefficient which is obtained on the basis of the compressible density profile throughout the core (Forte and Peltier, 1991). The complete set of free-slip CMB conditions in terms of the poloidal-flow vector v(r) ¼ [U0(r), Up(r), T0(r), TP(r), 1(r), g1(r)]Tr in system [35] is

No-slip. vða – Þ ¼ T P ða – Þy91 þ 1 ðRo Þy2 þ

819

! ½75

0 1 0 1 0 0 B C B C B C 0 B1C B C B C B C B C B 3 b C B0C B C B C 2, þ 1  C x1 ¼ B C x2 ¼ B B C B0C B C B C 0 B C B C B C B0C B C @ A 1 @ A 0 R, 0 1 0 B C 0 B C B C B C – gðb Þ B C B C C x3 ¼ B B C 0 B C B C B C 0 B  C @ d A ð2, þ 1Þ Ro

½77

820 Implications for Mantle Dynamics and Composition

where b ¼ o(bþ)  (b) ¼ 4.434 Mg m3 is the density jump across the CMB and db is the deflection (i.e., dynamic topography) of the CMB. The free-slip CMB condition for the toroidalflow vector u(r) ¼ [UT(r), TT(r)]Tr in system [36] is uðb þ Þ ¼ U T ðb þ Þw1

½78

in which the CMB basis vector is 1

w1 ¼

! ½79

0

to r ¼ r9, stopping along the way at all internal material boundaries (r ¼ ri) where we apply the internal matching conditions described previously (Section 1.23.2.3.3). We can similarly propagate each of the CMB basis vectors xi (i ¼ 1, 2, 3) in [77] from the CMB (r ¼ bþ) upward to r ¼ r9. The basic procedure is summarized schematically in Figure 6. At the location r ¼ r9 of the delta-function load we apply the matching conditions [42], [44], [45], and [47]–[49]. For each of these matching conditions, we employ the corresponding components of the poloidal-flow vectors v(r9þ) and v(r9), which are obtained from expressions [70], [72], and [76]: ( vðr 9þ Þ ¼

1.23.2.3.6 Determining viscous flow Green functions

The solution to the system of flow equations [35] may expressed in terms of Green functions which relate the poloidal flow velocity and the stress tensor at an arbitrary radius r ¼ r0 to a delta-function density load 1(r) ¼ d(r  r9) at any other radius r ¼ r9. For arbitrarily complex density (e.g., as in PREM) and viscosity profiles, the poloidal-flow system of eqns (35) must be integrated numerically. When r 6¼ r9, this linear system of equations is homogeneous: d vðr Þ ¼ MP ðr Þvðr Þ dr

ðwhen r 6¼ r 9Þ

½80

Here, MP (r) is the 6  6 matrix appearing in system [35]. We can propagate each of the surface boundary vectors yi (i ¼ 1, 2, 3) in [71] or [73], by numerically integrating [80] from the surface (r ¼ a) downward

Surface basis vectors

)

T P ða – Þy91 ðr þ 9Þ a da y ðr þ 9Þ þ  d 3

þ 1 ðRo Þy2 ðr þ 9Þ ½81

and vðr 9– Þ ¼ U P ðb þ Þx1 ðr 9– Þ þ 1 ðb þ Þx2 ðr 9– Þ þ

b db x3 ðr 9– Þ  d

½82

in which the surface and CMB basis vectors, yi (r9þ) and xi (r9), have been obtained by the numerical integration of [80], as outlined above. The application of the six matching conditions then yields the following system of equations: A p ¼ dp

½83

in which

Poloidal system y1 y 2

U P ða – Þy1 ðr þ 9Þ

Toroidal system

y3

r=a (Surface)

Surface basis vectors

z1

z2

internal { Apply boundary conditions

r = ri Apply matching conditions

internal { Apply boundary conditions

(Material boundary) r = ri r = r’ (Dirac δ-function)

}

r = ri CMB basis x1 vectors

x2

x3

{

Apply internal boundary conditions

CMB basis vector

w1

r=b (CMB)

Figure 6 Numerical integration of poloidal-flow system of eqns [35] (left) and toroidal-flow system of eqns [36] (right). The poloidal-flow basis vectors yi and xi are defined in expressions [71], [73], and [77]. The toroidal-flow basis vectors are defined in expressions [75] and [79].

Implications for Mantle Dynamics and Composition 0 B B B B B B B B p¼B B B B B B B @

U P ða – Þor T P ða – Þ 1 ðRo Þ a da  d U P ðb þ Þ  1 ðb þ Þ b db  d

1

0 1 0 C C B C C B C 0 C B C C B C C B C g ð r 9 Þ C B C C C ½84 Cand dp ¼ B B C 0 C B C C B C C B C 0 C B C C @ C d A A – ð2, – 1Þ Ro

Each column of the 6  6 matrix A in [83] involves the corresponding basis vectors yi (r9þ) and xi (r9) in [81] and [82]. For any given position r ¼ r9 of the delta-function load, we obtain a system of equations given by [83] which can be simply solved to find the unknown vector p. The elements of p define the poloidalflow impulse response (kernel) functions of the mantle which are discussed in Section 1.23.2.4. The toroidal-flow system [36] may be written as the following homogeneous equation: d uðr Þ ¼ MT ðr Þuðr Þ dr

½85

in which MT (r) is the 2  2 matrix appearing in system [36]. We can propagate each of the surface boundary vectors zi (i ¼ 1, 2) in [75], by numerically integrating [85] from the surface (r ¼ a) downward to r ¼ bþ, stopping along the way at all internal material boundaries (r ¼ ri) where we apply the internal matching conditions (50, 51). The basic procedure is summarized schematically in Figure 6. The resulting flow vector u(bþ) is uðb þ Þ ¼ V T z1 ðb þ Þ þ T T ða – Þz2 ðb þ Þ

½86

which must be matched to the toroidal flow in expression [79] uðb þ Þ ¼ U T ðb þ Þw1

½87

We thereby obtain the following simple system: B t ¼ dt

½88

in which t¼

U T ðb þ Þ T T ða – Þ

! and

dt ¼ V T z1 ðb þ Þ

½89

Each column of the 2  2 matrix B in [88] involves the corresponding basis vectors w1 and z2 in [86] and [87]. Expression [88] shows that the toroidal flow throughout the mantle is ‘driven’ by the surface toroidal flow VT generated by rotating tectonic plates,

821

which are themselves driven by the buoyancyinduced poloidal flow in the mantle. This surface coupling of poloidal and toroidal flow, due to the presence of surface tectonic plates, is discussed in the next section.

1.23.2.3.7 Incorporating tectonic plates as a surface boundary condition

The boundary conditions at the top of the mantle are more complex because of the presence of tectonic plates. If the plates are assumed to be an integral part of the underlying mantle flow and they fully participate in the flow, then a simple free-slip surface boundary condition would be appropriate. This free-slip assumption has often been employed in tomography-based flow studies (e.g., Hager et al., 1985; Forte and Peltier, 1987; Hager and Richards, 1989; King and Masters, 1992). If, on the other hand, we recognize that the plates are mechanically and rheologically distinct from the underlying mantle, then a more realistic surface boundary condition that explicitly treats the effective rigidity of plates is required. Essentially rigid tectonic plates are a fundamental aspect of convection dynamics in the Earth and, although their treatment lies beyond classical fluid mechanics theory, their dynamical impact on thermal convection in the mantle must be considered. One approach for modeling the coupling of convection and plates, developed by Ricard and Vigny (1989) and Gable et al. (1991), is based on an explicit treatment of the net vertical torque (or horizontal force) acting on the base of each surface plate. This method is a direct extension of the torque-balance analysis originally employed by Hager and O’Connell (1981) in their modeling of dynamic plate motions. The fundamental underlying assumption in the method developed by Ricard and Vigny is that the plate boundaries are completely stress-free. This method has been employed in a number of mantle-flow models over the past few years (e.g., Corrieu et al., 1994; Lithgow-Bertelloni and Richards, 1998). An alternative method for coupling the motions of rigid surface plates to buoyancy-induced mantle flow was developed by Forte and Peltier (1991, 1994). Only the main aspects of this method are summarized here. It may be shown that for a given geometry of surface plates, the internal density anomalies 1(r, , ) are partitioned into two families: 1 ðr ; ; Þ ¼ ˆ 1 ðr ; ; Þ þ 1 ðr ; ; Þ

½90

822 Implications for Mantle Dynamics and Composition

where the density perturbations ˆ 1 (r, , ) are obtained through a projection operator Pˆ as follows: ˆ 1 ts ðr Þ ¼ Pˆst ; lm ðr Þ1 ml ðr Þ

½91

The calculation of the projection operator Pˆst, lm(r) is given in Forte and Peltier (1994) and it depends on the geometry of the surface plates. The other component of the density anomalies 1 (r, , ) is simply given by the expression, 1 ts ðr Þ ¼ 1 ts ðr Þ – ˆ 1 ts ðr Þ

½92

As Forte and Peltier (1994) show, the poloidal mantle flow field produced by the component ˆ 1 is consistent with the geometry of possible rigid plate motions at the surface whereas the one produced by 1 is orthogonal to any possible plate motion. In other words, the plates participate in the underlying flow driven by ˆ 1, while they resist the flow produced by 1. Hence, free-slip (T p(rs) ¼ 0) and no-slip (V p(rs) ¼ 0) surface boundary conditions are applied to model the internal flows driven by ˆ 1 and 1, respectively. A simple application of the plate-projection operator is illustrated in Figure 7, where we compare the buoyancy driven plate motions predicted on the basis of two hypothetical hot thermal anomalies under the Pacific plate which differ only in their spatial relationship to the nearby East Pacific Ridge (EPR). Both ‘plumes’ are constructed in an ad hoc manner as rectangular-shaped, negative density anomalies with constant amplitude 0.1 g cm3 located in the depth range 400–750 km, equivalent to 2.6% relative density perturbation at 525 km depth. This amplitude is much greater than expected on the basis of the seismic tomography models (see Figure 12) and is only used here for the purpose of illustration. The predicted mantle flow field varies linearly with the amplitude of the driving density anomaly; therefore, a ‘plume’ anomaly which is 10 times smaller (0.01 g cm3) will produce surface motions 10 times smaller than those shown in Figure 7. In the absence of plates, both of the hypothetical plumes (Figures 7(a) and 7(b)) produce the same amplitude and pattern of horizontal flow divergence at the surface of the mantle (compare Figures 7(c) and 7(d)). The situation changes radically when the present-day configuration of rigid plates is imposed at the top of the mantle. The plume which is offset from the nearest plate boundary (Figure 7(a)), the EPR, drives a mantle flow field which produces only a small fraction of the horizontal surface divergence

obtained in the absence of plates (cf. Figures 7(c) and 7(e)). The plate projection operator [91] maps most of the original plume density anomaly 1 into the noslip family of density anomalies 1 which cannot produce observable plate motions. According to the alternative plate-motion theory employed by Ricard and Vigny (1989), this hypothetical plume anomaly produces flow-induced driving torques on the overlying Pacific plate which nearly cancel because the center of the upwelling is too far away from the nearest plate boundary. In contrast, the plume located under the EPR (Figure 7(b)) produces horizontal divergence of the Pacific and Nazca plates with an amplitude comparable to that obtained in the absence of plates (cf. Figures 7(d) and 7(f)). In this case, a substantial fraction of the original density anomaly 1 was mapped into the free-slip family of density anomalies ˆ 1. It should be noted that mutual interactions among the plates will generate non-zero far-field plate divergence and convergence (Figure 7(f), for example in the Indian Ocean) which is a consequence of the assumed rigidity of the plates. From the perspective of the tectonic plates, the existence of a null-space corresponding to the family of internal density anomalies d1 [92] which cannot drive observable motions (as in Figure 7(c)) constitutes a fundamental nonuniqueness in the interpretation of plate tectonics. The present-day plate motions (DeMets et al., 1990), or the geologic reconstructions of the history of plate motions (Gordon and Jurdy, 1986), therefore provide completely nonunique constraints on the mantle density anomalies generated by the thermal convection process. Other convection-related observables (e.g., global gravity anomalies) are required to constrain the density anomalies which fall into the null-space d1. A futher exploration of the implications of the non-unique interpretation of plate motions is presented by Forte and Peltier (1994). The plate motions which are driven by the buoyancy-induced (poloidal) mantle flow will necessarily pose a toroidal component. This toroidal component is directly coupled to the poloidal component of the plate motions as ðrˆ ?   vÞts ¼ Cst ; lm ðH ? vÞml

½93

in which the spherical harmonic coefficients of the radial vorticity of the plate velocity field v are linearly dependent on the harmonic coefficients of the horizontal divergence of the plate motions. The

Implications for Mantle Dynamics and Composition

(a) 0° 90°

(b) 30°

60°

90° 120° 150° 180° –150°–120° –90° –60° –30°

0°

0° 90°

60°

60°

30°

30°

0°

0°

–30°

–30°

–60°

–60°

–90°

–90°

–0.10

(c) 0°

30°

60°

–0.05

0.00 δρ

0.05

0.10

60°

60°

30°

30°

0°

0°

–30°

–30°

–60°

–60°

–90°

–90°

0°

–4 30°

60°

–2

2 0 rad/10 Ma

0°

0°

90°

60°

60°

30°

30°

0°

0°

–30°

–30°

–60°

–60°

–90°

–90°

–1

0 rad/10 Ma

1

30°

60°

(f)

4

90° 120° 150° 180° –150°–120° –90° –60° –30°

60°

90° 120° 150° 180° –150° –120° –90° –60° –30°

–0.10

0° 90°

(e)

90° 120° 150° 180° –150°–120° –90° –60° –30°

30°

(d)

90°

90°

823

–0.05

30°

60°

0.05

0.10

90° 120° 150° 180° –150° –120° –90° –60° –30°

–4 0°

0.00 δρ

–2

2 0 rad/10 Ma

–2

2 0 rad/10 Ma

0°

4

90° 120° 150° 180° –150°–120° –90° –60° –30°

–4

0°

0°

4

Figure 7 Interaction of hypothetical upwelling plumes with rigid surface plates. Maps (a) and (b) show two rectangularshaped, negative density anomalies (the hypothetical ‘plumes’) with amplitude 0.1 g cm3 assumed to be constant in the depth interval 400–750 km. Both plume anomalies have been expanded in terms of spherical harmonics up to degree and order 32. Maps (c) and (d) show the corresponding patterns of predicted horizontal surface divergence for a hypothetical mantle with a free-slip surface boundary (no plates) where the predicted motions are calculated up to harmonic degree 32. The mantle flow calculation was carried out assuming a depth-dependent viscosity variation shown below in Figure 9. Maps (e) and (f) show the corresponding patterns of plate divergence (up to harmonic degree 32) when the present-day configuration of rigid tectonic plates is imposed as a surface boundary condition in the mantle flow calculation.

vorticity and divergence fields are themselves directly dependent on the toroidal and poloidal scalar representations, respectively, of the plate velocity field (Forte and Peltier, 1987). It is important to note

that expression [93] is only valid for toroidal harmonic degree s 2 and that the coupling matrix Cst,lm is dependent only on the geometry of the plates (Forte and Peltier, 1994).

824

Implications for Mantle Dynamics and Composition

As discussed in Section 1.23.2.3.6, the viscous coupling between the plates and underlying mantle will result in the downward propagation of the surface toroidal flow into the mantle. The surface boundary condition [74] requires that the toroidal mantle flow match the toroidal component of plate motions, which is equivalent to their radial vorticity [93]: V t ð rs Þ ¼

rs ðrˆ ?   vÞml 1

½94

pffiffiffiffiffiffiffi where ¼ – 1. The Green function solutions of the poloidal and toroidal flow equations (see Section 1.23.2.3.6) may be used to describe the viscous response of the mantle to an arbitrary distribution of internal density anomalies 1 as follows: U, m ðr Þ T, m ðr Þ

! ¼

Z

rs

U , ðr ; r 9Þ

rCMB

T  , ðr ; r 9Þ

! 1 m, ðr 9Þdr 9

½95

m where U m (r) are the generalized sphe, (r) and T, rical harmonic coefficients of the velocity and stress fields [31] and [32], respectively, and U , (r, r9) and T  , (r, r9) are the corresponding Green functions. Tectonic plates may represent the most extreme manifestation of lateral variations in rheology in the convecting mantle. The plate coupling procedure described above thus provides an explicit means of incorporating the dynamical impact of such extreme lateral heterogeneity in the flow calculations which will be presented below. It is nonetheless true that eqns [35] and [36] governing flow in the mantle do not include lateral variations in viscosity and hence buoyancy forces will only directly excite a poloidal flow. In these calculations, toroidal mantle flow is excited passively from above, via the rotating surface plates. A more general flow theory which models the effects of lateral viscosity variations throughout the mantle (e.g., Martinec et al., 1993; Forte and Peltier, 1994) is required to describe the direct excitation of toroidal flow by mantle buoyancy forces. The dynamical implications of large-amplitude lateral viscosity variations throughout the mantle have recently been studied by Moucha et al. (2006). Despite the considerably more complex theory (relative to equations 35, 36) needed to calculate flow in a 3-D viscosity field, Moucha et al. find the impact on predicted surface observables such as topography and geoid are relatively small. This recent study provides an important verification for the modeling and interpretation of geodynamic observables using the

plate-coupled flow theory with purely depth-dependent viscosity in the mantle. The results obtained on the basis of this simplified treatment of mantle dynamics are presented in the sections that follow. 1.23.2.4 Geodynamic Response Functions for the Mantle The theoretical relationship between the mantle density anomalies and the principal convection related surface observables (i.e., geoid or gravity anomalies, dynamic surface and CMB topography, plate motions) may be summarized in terms of geodynamic response or kernel functions (see Hager and Clayton (1989) and Forte (2000) for reviews). These kernel functions are calculated in the spherical harmonic spectral domain using the viscous flow Green functions defined in expression [95]. The principal assumption we shall make in calculating these geodynamic kernels – in addition to the assumption of purely depth-dependent viscosity below the plates – is that the 670 km seismic discontinuity is not a barrier to radial mantle flow. This whole-mantle flow assumption has been recently subjected to a series of direct tests based on joint inversions of global seismic and geodynamic data sets (Simmons et al., 2005). It is found that the most satisfactory reconciliation of the global seismic and geodynamic constraints on 3-D mantle structure is obtained with whole-mantle flow and significantly poorer fits are obtained with flow models which assume vertical flow barriers at 670 km depth or within the lower mantle (e.g., at 1800 km depth). The geoid kernels G,( ; r) relate the spherical harmonic coefficients of the nonhydrostatic geoid field, dN m, , to the spherical harmonic coefficients of the density perturbations dm, (r) as follows: dN,m ¼

3 ð2, þ 1Þ 

Z

rs

G, ð ; r Þdm, ðr Þdr

½96

rCMB

where  ¼ 5.515 Mg m3 is the mean density of the Earth, rs is the mean radius of the solid surface, and rCMB is the mean radius of the CMB. The geoid kernels are a functional of the nondimensional radial mantle viscosity profile (r)/ o. The free-air gravity anomaly coefficients dGm, are related to the nonhydrostatic geoid coefficients dN m, as follows: dG,m ¼ ð, – 1Þ

go dN m Ro ,

½97

where g0 ¼ 9.82 m s2 (982 000 mGal) is the mean gravitational acceleration at Earth’s surface and

Implications for Mantle Dynamics and Composition

R0 ¼ 6371 km is Earth’s mean radius. Through the combination of equations [96] and [97], we can relate the gravity anomaly coefficients to the internal density perturbations. As discussed in the preceding section, a dynamical flow model which incorporates rigid plate motions must include mixed free-slip and no-slip surface boundary conditions. In the calculation of the flowrelated observables, we must therefore determine the separate contributions provided by the density anomalies dˆ and d  which are convolved with kernel functions derived with free-slip and no-slip surface boundaries, respectively. For example, in the case of the nonhydrostatic geoid anomalies, we rewrite eqn [96] as follows: dN,m ¼

Z rs 3 Gˆ, ð ; r Þdˆ m, ðr Þ ð2, þ 1Þ  rCMB

þ G, ð ; r Þd m, ðr Þ dr

½98

where the geoid kernels Gˆ,( ; r) and G,( ; r) are calculated with free-slip and no-slip surface boundaries, respectively. As shown in Forte and Peltier (1994), the density anomalies dˆ are spatially correlated with the positions of divergent or convergent plate boundaries. Subducting slabs and mid-ocean ridge heterogeneity are therefore resolved as part of the dˆ anomalies. The surface topography kernel functions T,( ; r) relate the spherical harmonic coefficients of the flowinduced surface topography dam, to the spherical harmonic coefficients of the internal density anomalies as follows: da,m ¼

1 mo

Z

rs



m ðr Þ dr Tˆ, ð ; r Þdˆ m þ T, ð ; r Þd ,

,

½99

rCMB

where mo ¼ 2.2 Mg m3 is the density jump across the mantle–ocean boundary and Tˆ,( ; r) and T ,( ; r) are calculated with free-slip and no-slip surface boundaries, respectively. Similarly, the spherical harmonic coefficients of the flow-induced CMB topography dbm, may be expressed in terms of CMB topography kernels B,( ; r) as follows: db,m ¼

1 cm

Z

rs



m, ðr Þ dr Bˆ, ð ; r Þdˆ m, þ B, ð ; r Þd

½100

rCMB

where Bˆ,( ; r) and B,( ; r) are calculated with free-slip and no-slip surface boundaries, respectively, and cm ¼ 4.43 Mg m3 is the density jump across the CMB. As in the case of the nonhydrostatic geoid in

825

eqn [98], the dependence on viscosity of the predicted surface and CMB topography appears explicitly through the nondimensional profile (r)/ o. Surface-plate velocities v may be characterized by their horizontal divergence field H ? v (Forte and Peltier, 1987). The horizontal divergence kernels D,( ; r) relate the spherical harmonic coefficients of the predicted plate divergence (H ? v)m, to the spherical harmonic coefficients of the density perturbations as follows: ðH ? vÞm, ¼

go o

Z

rs

D, ð ; r Þdˆ m, ðr Þdr

½101

rCMB

in which the free-slip component of the mantle denˆ defined in expression [91], is sity anomalies d, required to model the surface divergence of the plates. We also note the presence of the reference scaling value of viscosity 0. This presence implies that plate motions, unlike the other observables in eqns [96]–[100], will be dependent on the absolute value of mantle viscosity.

1.23.2.5 Depth Dependence of Mantle Viscosity Understanding the long timescale rheology of the mantle as represented by its effective viscosity is a central and enduring problem in global geophysics. The diverse methods and data sets which have been employed to constrain mantle viscosity have been a source of ongoing contention and debate. The importance and intensity of this debate are an apt reflection of the fundamental role of mantle viscosity in controlling a wide array of geodynamic processes. For example, millennial timescale glacial isostatic adjustment (GIA) processes such as Pleistocene and Holocene sea-level variations and related anomalies in Earth’s gravitational field and rotational state are known to be strongly dependent on the depth dependence of mantle viscosity. On much longer, million to hundred-million year timescales, viscosity exerts fundamental control on the dynamics of mantle convection and on the corresponding evolution of the thermal and chemical state of Earth’s interior. The very long timescale implications of mantle viscosity also include fundamental surface geological and geophysical processes, such as global scale epeirogeny and associated sea-level changes, global geoid anomalies, and tectonic plate motions. The spatial variation of viscosity in Earth’s mantle can, in principle, be determined from a consideration

826 Implications for Mantle Dynamics and Composition

of the microphysical properties which control the steady-state creep of mantle rocks (e.g., Nicolas and Poirier, 1976; Sammis et al., 1977). In practice, however, the microphysical models of mantle viscosity depend on knowledge of a number of critical physical properties, such as grain size and activation volume and energy, which are poorly known in the mantle and hence require extrapolations of laboratory creep experiments which greatly exceed the physical conditions (i.e., pressure, temperature, or creep rate) under which these experiments were originally performed. It is nonetheless possible to show, on the basis of microphysical considerations, that the global horizontal average viscosity may increase by up to three orders of magnitude across the mantle (e.g., Ranalli 2001). More precise determinations of the radial viscosity profile require additional constraints from geophysical data sets which are sensitive to the long-term rheology of the mantle. The first and most influential geophysical contribution to our understanding of mantle viscosity is Haskell’s (1935) study of Fennoscandian postglacial uplift which was found to require an average viscosity of 1021 Pa s down to a depth approximately equal to the horizontal dimension of the surface load (1000–1500 km). McConnell’s (1968) spectral analysis of the Fennoscandian uplift data led to the inference that mantle viscosity increased significantly with depth while at the same time still satisfying the Haskell average. The Haskell average value of viscosity over the upper 1000–1200 km of the mantle was also supported in subsequent work by O’Connell (1971). These initial studies which supported the compatibility between the Haskell constraint on viscosity and significant increases of viscosity in the lower mantle were displaced by subsequent analyses by Cathles (1975) and Peltier and Andrews (1976) who concluded that viscosity was nearly constant from the base of the lithosphere to the CMB. This view of mantle viscosity characterized by only moderate increases with depth was reinforced in later studies involving larger GIA data sets (e.g., Wu and Peltier 1983; Tushingham and Peltier, 1991). A review of the GIA constraints on viscosity may be found in Mitrovica (1996). The first inferences of the depth variation of mantle viscosity based on surface data sets associated with mantle convection (e.g., Hager, 1984; Richards and Hager, 1984; Ricard et al., 1984) disagreed strongly with the GIA inferences of mantle viscosity which prevailed since the earliest studies by Cathles, Peltier, and collaborators. These convection-based

analyses of viscosity, in particular the first study of Earth’s global geoid anomalies by Hager (1984), suggested a large increase of viscosity in the sublithospheric mantle with an approximately factor of 30 increase in the average viscosity of the lower mantle relative to the upper mantle. This conclusion was supported by subsequent studies by Forte and Peltier (1987) and Forte et al. (1991) who extended the viscous flow modeling to include the plate velocities. The geodynamic inferences for significant (1–2 order of magnitude) increases in the average viscosity of the lower mantle relative to that of the upper mantle were reinforced by subsequent analyses carried out by Ricard et al. (1989), Ricard and Vigny (1989), Hager and Clayton (1989), Forte and Peltier (1991), Forte et al. (1993a, 1994), King and Masters (1992), Corrieu et al. (1994), and Thoraval and Richards (1997). A survey of some of the earliest viscosity inferences derived on the basis of global geoid anomalies are shown in Figure 8. Reviews of mantle viscosity inferences based on either convective flow or GIA data sets may by found in King (1995) and Kaufmann and Lambeck (2000). Additional geodynamic considerations which suggest a substantial increase of viscosity with depth are based on a wide variety of analyses which include the stability of the hot spot reference frame (e.g., Richards, 1991; Que´re´ and Forte, 2006), long-term rates of polar wander (e.g., Sabadini and Yuen, 1989; Spada et al., 1992; Richards et al., 1997; Steinberger and O’Connell, 1997), and the planform of mantle convection (e.g., Zhang and Yuen, 1995; Bunge et al., 1996; Forte and Mitrovica, 2001). The inference of mantle viscosity from convection related surface data, considered in the context of tomography-based mantle flow models, was initially carried out through an effectively trail-and-error fitting procedure (e.g., Forte and Peltier 1987, 1991; Hager and Clayton, 1989; Forte et al., 1993). Subsequently, a number of formal mathematical inversions of the geodynamic data were carried out (e.g., Forte et al., 1991; Ricard and Wuming, 1991; King and Masters, 1992; Corrieu et al., 1994; King, 1995; Forte and Mitrovica, 1996; Kido et al., 1998; Forte, 2000; Panasyuk and Hager, 2000). The viscosity profiles obtained by inverting geodynamic data sets (Figure 8) are all characterized by strong overall increases in viscosity with depth, but it is also clear that there are significant differences in the detailed character of these viscosity inferences. There are many factors which can contribute to these differences and they include, most notably, the use of:

Implications for Mantle Dynamics and Composition

(a)

Hager et al. (1985)

(b) Forte and Peltier (1987, 1991)

(c) Ricard and Wuming (1991)

(e)

(f)

827

0

Depth (km)

500 1000 1500 2000 2500 King and Masters (1992)

(d)

Ricard et al. (1993)

Corrieu et al. (1995)

0

Depth (km)

500 1000 1500 2000 2500 10–1

(g)

101 102 100 Nondimensional viscosity

Forte and Mitrovica (1996)

10–1

(h)

100 101 102 Nondimensional viscosity Lambeck et al. (1996)

10–1

(i)

100 101 102 Nondimensional viscosity Lambeck et al. (1998)

0

Depth (km)

500 1000 1500 2000 2500 1019 1020 1021 1022 1023 Absolute viscosity (Pa s)

1019 1020 1021 1022 1023 Absolute viscosity (Pa s)

1019 1020 1021 1022 1023 Absolute viscosity (Pa s)

Figure 8 A survey of previous geodynamic inferences of depth dependent mantle viscosity. Frames (a) to (f) show relative (nondimensional) variations of viscosity with depth derived on the basis of global long wavelength geoid anomalies. Frames (g) to (i) show radial (absolute) viscosity profiles inferred on the basis of GIA data (in addition to geoid data in the case of frame g). The GIA data sets constrain the absolute value of viscosity whereas the steady-state viscous flow modelling of the geoid data only constrain relative variations of viscosity with depth. (Figure adapted from Kaufmann and Lambeck, 2000).

different tomography-based inferences of mantle density anomalies; different combinations of geodynamic data (e.g., geoid only, geoid plus plate motions, geoid plus dynamic topography, all of the above); different inversion methods (e.g., Monte-Carlo, Bayesian, genetic, Occam algorithms). The formulation of the viscous flow models themselves can also play a critical role in explaining these differences, since it is known that the use of different surface boundary conditions (e.g., free-slip, no-slip, freely

rotating tectonic plates) will lead to rather different viscosity inferences (e.g., Thoraval and Richards, 1997). Finally, the limited resolving power of the geodynamic constraints precludes a direct interpretation of viscosity at any given depth in the mantle. For example, Forte and Mitrovica (2001) found that the combined constraints provided by global free-air gravity anomalies, tectonic plate motions, and excess ellipticity of the CMB could only resolve average values of viscosity in rather broad depth intervals.

828 Implications for Mantle Dynamics and Composition

The discordant viscosity inferences derived from either GIA or convection-related data sets were initially regarded as evidence for the transient nature of mantle viscosity over the timescales that separate these two processes (e.g., Sabadini et al., 1985; Peltier, 1985). The necessity for time-dependent viscosity in the mantle was later weakened by two independent analyses of the GIA data sets. The first argument against transient viscosity was provided by an analysis of differential Late Holocene sea-level high-stands in the Australian region by Nakada and Lambeck (1989), who showed that these GIA data required a nearly two order of magnitude increase of viscosity across the mantle. The second independent argument was provided by Mitrovica (1996), who showed that previous inferences of a nearly constant-viscosity mantle (e.g., Wu and Peltier 1983; Tushingham and Peltier, 1991) were conditioned by a misinterpretation of the original Haskell (1935) constraint. The latter was erroneously assumed to apply only to the top 670 km of the mantle (i.e., the upper mantle), rather than to the average value of viscosity down to depths of about 1000–1500 km (Mitrovica, 1996). The recognition that GIA data are not incompatible with large increases of mantle viscosity with depth spurred new efforts to reconcile GIA and convection data sets with a single profile of mantle viscosity. The first efforts were undertaken by Forte and Mitrovica (1996) and Mitrovica and Forte (1997). Despite the very different time and spatial scales over which these processes operate, these authors found that it was possible to simultaneously explain both decay times associated with the postglacial uplift of Hudson Bay and Fennoscandia, and long-wavelength free-air gravity anomalies associated with mantle convection with a single viscosity profile characterized by an approximately two order of magnitude increase with depth which also satisfied the original Haskell constraint on mantle viscosity. These joint GIA-convection inversions culminated in the recent study by Mitrovica and Forte (2004) in which the range of viscosity profiles consistent with both families of surface data are summarized below in Figure 9. The convection data employed in these viscosity inversions are shown in Figure 2. The GIA data include the Fennoscandian relaxation spectrum (FRS) and a set of decay times determined from the postglacial sea-level history in Hudson Bay and Sweden. The inverted viscosity profiles are characterized by a 3 orders of magnitude increase from the upper mantle (mean value of 4  1020 Pa s) to a highviscosity (>1023 Pa s) peak at 2000 km depth. Below

2000 km, the viscosity shows a significant, 2–3 orders of magnitude, reduction toward the CMB. Similar radial variations in lower-mantle viscosity have been inferred in other independent studies (e.g., Ricard and Wuming 1991). The preferred viscosity profile (solid black line, Figure 9) provides an optimal fit to both (GIA and convection) families of data and, unless stated otherwise, this profile will be employed to calculate all convection related related surface observables presented below. The first efforts to reconcile the constraints on viscosity provided by both glacial isostatic adjustment (GIA) and convection data were undertaken by Forte and Mitrovica (1996) and Mitrovica and Forte (1997). Despite the very different time and spatial scales over which these processes operate, these authors found that it was possible to simultaneously explain both GIA and convection data with a single viscosity profile. These joint GIA–convection inversions culminated in the recent study by Mitrovica and Forte (2004) in which the range of viscosity profiles consistent with both families of surface data is summarized in Figure 8. The convection data employed in these viscosity inversions are shown in Figure 2. The GIA data include the Fennoscandian relaxation spectrum (FRS) and a set of decay times determined from the postglacial sealevel history in Hudson Bay and Sweden. The inverted viscosity profiles are characterized by a 3 orders of magnitude increase from the upper mantle (mean value of 4  1020 Pa s) to a high-viscosity (>1023 Pa s) peak at 2000 km depth. Below 2000 km, the viscosity shows a significant, 2–3 orders of magnitude reduction toward the CMB. Similar radial variations in lower-mantle viscosity have been inferred in other independent studies (e.g., Ricard and Wuming, 1991). The preferred viscosity profile (solid black line, Figure 8) provides an optimal fit to both (GIA and convection) families of data and, unless stated otherwise, this profile will be employed to calculate all convection-related surface observables presented below. The geodynamic response or kernel functions calculated on the basis of the preferred mantle viscosity profile (Figure 9) are illustrated below in Figure 10. Comparing the free-slip and no-slip geoid kernels, we note that the largest contribution from density anomalies in the top half of the mantle is provided by those anomalies dˆ which are efficient in driving observable plate motions. These anomalies correspond to subducting slabs and upper-mantle plumes below the mid-ocean ridges. In the bottom

Implications for Mantle Dynamics and Composition

829

0

500

Depth (km)

1000

1500 Inferred viscosity bounds

2000

2500

3000 18 10

1019

1020 1021 1022 Absolute viscosity (Pa s)

1023

1024

Figure 9 Results of inversions of the GIA and convection data sets. Full details of the inversions are presented in Mitrovica and Forte (2004). The convection data are interpreted in terms of a tomography-based mantle flow model in which the density anomalies are derived from the Grand (2002) tomography model (Figure 3(b)). The thick gray lines illustrate the range of allowable values of mantle viscosity which are consistent with the joint GIA–convection constraints. The solid black line is the preferred viscosity profile, on the basis of the Grand (2002) tomography model, which provides an optimal fit to the entire suite of geodynamic data. The dash-dotted line is the preferred viscosity profile obtained on the basis of the joint shear–bulk-sound tomography model of Masters et al. (2000).

half of the mantle, the density anomalies belonging to the d  family provide the largest contribution to the surface geoid or gravity anomalies. These d  anomalies produce no observable plate motions. At sufficiently short horizontal wavelengths (corresponding to harmonic degrees , 32), we note that the distinction between dˆ and d  anomalies begins to disappear: both provide equal contributions to the surface geoid or gravity anomalies. The free-slip and no-slip surface topography kernels show that the density anomalies d  provide the strongest contribution to dynamic surface topography, especially at long horizontal wavelengths (corresponding to the degree range , ¼ 28). At short horizontal wavelengths, the distinction between the two classes of density anomalies again ceases to be important.

1.23.3 Modeling Geodynamic Observables with Seismic Tomography We now consider the extent to which the 3-D mantle structure resolved in recent seismic tomographic models may be used to explain the global convectionrelated data sets shown in Figure 2. The geodynamic response functions obtained on the basis of the preferred radial viscosity profile (Figure 10) will be used to establish the connection between the seismically inferred heterogeneity in the mantle and the surface geodynamic observables. A broad selection of seismic shear-velocity models, obtained in independent studies using diverse and complementary seismic data and using different seismic modeling and inversion techniques, will be employed here. Seismic shear

830 Implications for Mantle Dynamics and Composition

velocity heterogeneity is the initial focus of this tomography-based modeling because it is expected that S-wave anomalies dVS will be most sensitive to the temperature anomalies maintained by the mantle convection process (e.g., Ro¨hm et al., 2000; Forte and Perry, 2000; Forte et al., 2002).

1.23.3.1

the radial direction and spherical harmonics to equalarea cells in the horizontal direction. To facilitate a direct comparison among these seismic heterogeneity models and the subsequent usage of the models in viscous flow calculations, all models have been projected onto a common set of radial and horizontal basis functions as follows: 32 X 20 X þ, dVS X m m ¼ k c, Tk ðx ÞY, ð; jÞ VS k¼0 ,¼1 m¼ – ,

Seismic Heterogeneity Models

The S-wave tomography models employed here are, in chronological order: model S20A (isotropic version) from Ekstro¨m and Dziewonski (1998), model S20RTS from Ritsema et al. (1999), model SAW24 from Me´gnin and Romanowicz (2000), model SB4_L18 from Masters et al. (2000), model TX2002 from Grand (2002), and model J362D28 from Antolik et al. (2003). The spatial variation of seismic heterogeneity in these S-wave models is parametrized in terms of a wide assortment of basis functions, ranging from B-splines to piecewise-discontinuous layers in Viscosity

(a)

in which Tk(x) are Chebychev polynomials in the normalized radius x ¼ (2r  rsur f  rCMB)/(rsur f  rCMB) and Ym, are ordinary spherical harmonics. Notice that all models are expanded up to order 32 in the radial Chebychev polynomials and up to degree 20 in the surface spherical harmonics. When each of the S-wave heterogeneity models are projected onto the representation (102), the depth variations of the rms amplitudes of the velocity anomalies are shown in Figure 11. Here we note (b)

Divergence kernels

0

Depth (km)

500 1000 1500 l=1 l=2 l=4 l=8 l = 16 l = 32

2000 2500 3000 (c)

1020

1021

1022

1023

Gravity kernels (free-slip)

–0.08 –0.06 –0.04 –0.02 (d)

Gravity kernels (no-slip)

0

Depth (km)

500 1000 1500 2000 2500 3000 Figure 10 (Continued )

–0.2

0

0.2

½102

–0.2

0

0.2

0

Implications for Mantle Dynamics and Composition

(e)

Surf. top. kernels (free-slip)

(f)

831

Surf. top. kernels (no-slip)

0

Depth (km)

500 1000 1500 2000 2500 3000 (g)

l=1 l=2 l=4 l=8 l = 16 l = 32

–1.0 –0.8 –0.6 –0.4 –0.2 0.0 CMB top. kernels (free-slip)

–1.0 –0.8 –0.6 –0.4 –0.2 (h)

0.0

CMB top. kernels (no-slip)

0

Depth (km)

500 1000 1500 2000 2500 3000

0.0 0.2 0.4 0.6 0.8 1.0 1.2

0.0 0.2 0.4 0.6 0.8 1.0 1.2

Figure 10 Geodynamic kernel functions. The geodynamic response functions shown here are calculated on the basis of the preferred GIA–convection-inferred viscosity profile shown in (a) (same as the solid black curve in Figure 7). The kernels corresponding to harmonic degrees , ¼ 1, 2, 4, 8, 16, 32 are identified by the legend in (b) (repeated below in (e)). The horizontal divergence kernels, defined in [101], are shown in (b). The gravity kernels, defined in [98], for free-slip and no-slip surface boundary conditions are shown in (c) and (d), respectively. The surface topography kernels, defined in [99], for freeslip and no-slip conditions are shown in (e) and (f) respectively. The kernels for the dynamic CMB topography, defined in [100], are shown in (g) and (h) for free-slip and no-slip conditions, respectively.

that, with the exception of the top 300 km in the mantle, the average amplitudes of the relative perturbations of seismic shear velocity dVS/VS show significant differences in the mantle, especially in the transition zone region (400–1000 km depth) and at the base of the mantle (2500–2891 km depth). 1.23.3.2

Mantle Density Anomalies

The mantle density anomalies required in modeling the convection-related surface observables may be derived from the seismic tomography models using experimental and theoretical results from mineral ˇ adek et al., 1994; Wang physics (e.g., Karato, 1993; C and Weidner, 1996; Sobolev et al., 1997; Stacey, 1998; Zhang and Weidner, 1999; Trampert et al., 2001;

Stixrude and Lithgow-Bertelloni, 2001; Jackson, 2001; Karato and Karki, 2001; Oganov et al., 2001a). This derivation would, however, require that the amplitudes of seismic velocity anomalies are well constrained and that all relevant mineralogical variables in the mantle (e.g., reference composition, temperature, and equation of state parameters) are sufficiently well known. It is unlikely that all of these requirements can be met at the present time, although previous efforts to simultaneously interpret mantle density and seismic velocity anomalies using mineral physics data have shown some promise (e.g., Forte et al., 1994). A second, alternative method for inferring mantle density anomalies is based on the direct inversion of the convection-related geophysical observables

832 Implications for Mantle Dynamics and Composition

3.0 J 362D 28 SAW 24 S 20A-Iso TX 2002 S 20RTS SB 4_L18

2.5

δVS/ VS (%)

2.0

1.5

1.0

0.5

0.0

0

500

1000

1500 Depth (km)

2000

2500

3000

Figure 11 Amplitude of S-wave heterogeneity in the mantle. Plotted here is the root-mean-square (rms) amplitude of relative perturbations of seismic shear velocity obtained from the following tomography models: J362D28 (Antolik et al., 2003), SAW24 (Me´gnin and Romanowicz, 2000), S20A-Iso (Ekstro¨m and Dziewonski, 1998), TX2002 (Grand, 2002), S20RTS (Ritsema et al., 1999), SB4_L18 (Masters et al., 2000). All tomography models have been projected onto the spatial parametrization given by expression [102].

(e.g., Forte et al., 1993; Corrieu et al., 1994; Karpychev and Fleitout, 2000; Forte and Perry, 2000; Panasyuk and Hager, 2000; Forte and Mitrovica, 2001; Pari, 2001; Deschamps et al., 2001). This approach is useful when contending with significant uncertainties in the amplitudes of the seismic anomalies, as is evident in Figure 11. This figure underscores a key concern, namely that the use of mineral physics data to translate the seismic anomalies can lead to large variations in the estimated mantle density anomalies, depending on the choice of tomography model. To avoid this difficulty, it is preferable to carry out simultaneous inversions of all global geodynamic data (Figure 2) to determine an optimal velocity–density scaling factor d ln /d lnVS (e.g., Forte and Mitrovica, 2001), such that d ¼ 



 d ln  dVS d ln VS VS

½103

where dVS/VS are shear velocity anomalies obtained from a tomography model and d/ are the relative perturbations in density which are constrained by the geodynamic data. This approach has also been employed in other tomography-based mantle flow models (e.g., Hager and Clayton, 1989; Corrieu et al., 1994; Panasyuk and Hager, 2000).

The geodynamic data provide direct, linear constraints on the density anomalies in the mantle (see eqns [98]–[100]), which may be generically expressed as follows: dO,m ¼ f,  f,

Z

a

b N X

K, ð ; r 9Þð1 Þm, ðr 9Þdr 9 K, ð ; ri Þð1 Þm, ðri Þwi

½104

i¼1

in which dOm, are the spherical harmonic coefficients of a geodynamic observable (e.g., geoid or free-air gravity anomalies), K,( ; r9) is the corresponding kernel function, and f, is a factor which depends on the geodynamic observable (see, for example, eqn [98]). The numerical calculation of integrals is usually carried out with equivalent finite sums, such as in [104] where wi is a weighting term which depends on the numerical summation algorithm (e.g., Gauss– Legendre quadrature). On the basis of expression [103], we may further rewrite equation [104] as dO,m ¼ f,

N X i¼1

wi K, ð ; r i Þo ðr i Þ

    dVS m d ln  ðr i Þ ðr i Þ d ln Vs VS , ½105

in which dVS/VS are the shear-velocity anomalies in the mantle and d ln /d ln VS is the corresponding

Implications for Mantle Dynamics and Composition

velocity-to-density scaling which is assumed to vary with depth in the mantle. Expression [105] provides the basis for a discrete, linear inversion of the geodynamic data dOm, to find an optimal velocity-to-density scaling profile d ln / d lnVS. A series of Occam inversions (Constable et al., 1987) of the geodynamic data (Figure 2) were carried out to find the smoothest family of d ln /d lnVS profiles which are consistent with the data. For each Occam inversion, a different seismic shear-wave tomography model (Figure 11) was employed, and in each case the same geodynamic kernel functions (Figure 10) were employed to relate the mantle heterogeneity to the geodynamic data. The results of these inversions are summarized in Figure 12,

833

which shows the optimal profiles of d ln /d lnVS and the resulting depth variation of the rms amplitude of the lateral density anomalies implied by expression [103]. Mineral physics estimates of the velocity–density scaling (Figure 12) show significant departures from the scaling profiles derived from the Occam inversions. As discussed above, this discrepancy may in part arise from the inadequate resolution of seismic heterogeneity in the mantle–especially in the mid-mantle region between 1000 km and 2000 km depth–and it may also arise from inadequacy of the theoretical assumptions inherent in the mineral physics estimates. In the latter case, the chief assumption which may be deficient is the dominance of thermal effects on both the seismic and

(a)

d ln ρ/d ln VS

0.25 0.20 0.15 0.10 0.05 0.00 0

500

1000

1500

2000

2500

3000

2500

3000

(b) J 362D 28 SAW 24B16 S 20A-ISO TX 2002 S 20RTS SB 4_L18 K & K (2001)

δρ/ρ (%)

0.3

0.2

0.1

0

0

500

1000

1500 2000 Depth (km)

Figure 12 Occam inversions for mantle density anomalies. (a) The optimal Occam inferred velocity–density scaling profiles for each of the tomography models illustrated in Figure 11. The legend identifying the results is in (b). For comparison are shown (thick dashed gray lines) a series of theoretical estimates of the depth variation of d ln/d lnVS in the lower mantle obtained by Karato and Karki (2001) on the basis of mineral physics data. These mineral physical estimates assume that both density and seismic anomalies are produced by temperature perturbations in the mantle. (b) The resulting inferences of the depth variation of the rms density anomalies which are calculated on the basis of the seismic shear velocity anomalies using expression [103].

834

Implications for Mantle Dynamics and Composition

density anomalies in the mantle. Compositional effects, especially on mantle density heterogeneity, may be one important explanation for the deviations from theoretical mineral physics estimates based on temperature effects alone (e.g., Stacey, 1998). The possible importance of compositional effects will be discussed below in the context of a further investigation of tomography-based geodynamic constraints on mantle density.

1.23.3.3

Predicted Tectonic Plate Motions

In this and the next few sections, we consider the extent to which the tomography-based inferences of mantle density anomalies presented in Figure 12 provide a successful explanation of the convectionrelated surface observables (Figure 2). We begin here with a consideration of convection-induced tectonic plate motions. The predicted horizontal divergence of the plate motions, calculated on the basis of expression [101] and using the kernels in Figure 10(b), are presented in Figure 13 alongside the observed NUVEL-1 divergence of the plates. A detailed spectral comparison between the predicted and observed plate divergence is shown on Figure 14. A quantitative summary of the overall agreement between the predicted and observed plate divergence, quantified in terms of total rms amplitude, total spatial correlation, and variance reduction, is provided in Table 1. The variance reduction is defined as follows:

systematic bias to low predicted amplitudes is a natural consequence of the least-squares misfit criterion which is employed in almost all inversion algorithms (including the Occam approach) which leads to significantly damped inferences of the d ln /d ln VS scaling profiles. It can be readily shown that the inherent damping of tomography-based inferences of d ln /d ln VS arises as a consequence of poorly resolved 3-D seismic structure in geodynamically important regions of the mantle, such as the mid-mantle region (e.g., Forte et al., 1994). We can illustrate this in an inverse problem in which we seek, for the sake of simplicity, the best-fitting constant d ln /d ln VS value. As shown in expression [104], the predicted convection-related observables are linearly dependent on the internal density perturbations in the mantle. By virtue of expression [105], we may write  d ð; jÞ ¼

where d(, j) represents the surface data and p(, j) is the corresponding mantle-flow prediction assuming d ln /d ln Vs ¼ 1. We can then show that the optimal d ln /d ln VS value that minimizes the least-squares misfit between data and the predictions is given by the following expression: d ln  d ¼ c ðd ; p Þ d lnVS p

where (O  P)m, and Om, are the harmonic coefficients of the difference between the observed and predicted fields and the observed fields, respectively. An important characteristic of the predicted divergence fields shown in Figure 13, which is also confirmed in Table 1, is the subdued amplitude of the predicted fields relative to the data. As will be seen in the next few sections, all tomography-based predictions of geodynamic observables which are calculated on the basis of the optimal Occaminverted density–velocity scaling profiles (Figure 12) are systematically reduced in amplitude relative to the amplitude of the observed fields. This

½107

where d and p are the rms amplitudes of the data and the predictions, ZZ 1 jd ð; jÞj2 d2 S; 4 S1 ZZ 1 2p ¼ jpð; jÞj2 d2 S 4 S1

2d ¼

var: red: ¼ 100% " P P # þ, m m , m¼ – , ðO – P Þ, ðO – P Þ,  1– P Pþ, m m , m¼ – , O, O, ½106

 d ln  pð; jÞ d lnVs

and c ðd ; pÞ ¼

1 4

ZZ

d ð; jÞpð; jÞd2 S

½108

S1

where all integrals are defined over the surface of the unit sphere S1. The quantity c(d, p) defined in [108] is the cross-correlation between the data and the predictions. Equation [107] thus shows that the inferences of d ln/d ln VS will be directly proportional to the spatial correlation between the surface data and the corresponding geodynamic integral of the 3-D structure in the seismic tomography model. As previously shown by Forte et al. (1994), this correlation may be significantly degraded if the 3-D seismic structure is poorly resolved.

Implications for Mantle Dynamics and Composition

NUVEL-1 divergence (L = 1–32) (a)

0°

(b) 0°

30° 60° 90° 120° 150° 180° –150°–120° –90° –60° –30°

0°

90°

J 362D 28 divergence (L = 1–32) 30°

60° 90° 120° 150° 180° –150°–120° –90° –60° –30°

30°

30° 0°

0°

–30°

–30°

–60°

–60° –90°

–90°

SAW 24B 16 divergence (L = 1–32) (c)

0°

30° 60°

90° 120° 150° 180° –150°–120° –90° –60° –30°

0°

(d)

90°

90°

60°

60°

30°

30°

0°

0°

–30°

–30°

–60°

–60°

S 20A -Iso divergence (L = 1–32) 0° 30°

60°

90° 120° 150° 180° –150°–120° –90° –60° –30°

0°

30° 60°

90° 120° 150° 180° –150°–120° –90° –60° –30°

0°

30°

0°

–90°

–90°

TX 2002 divergence (L = 1–32) 0°

S 20RTS divergence (L = 1–32)

30° 60° 90° 120° 150° 180° –150°–120° –90° –60° –30°

0°

(f)

90°

90°

60°

60°

30°

30°

0°

0°

–30°

–30°

–60°

–60°

0°

–90°

–90°

SB 4_L18 divergence (L = 1–32) (g)

0°

90° 60°

60°

(e)

835

0°

30° 60°

0°

90°

90°

60°

60°

30°

30°

0°

0°

–30°

–30°

–60°

–60°

–90°

–90°

–10

‘Mean’ divergence (L = 1–32)

(h)

90° 120° 150° 180° –150° –120° –90° –60° –30°

–5

0 rad/100 Ma

5

60° 90° 120° 150° 180° –150° –120° –90° –60° –30° 0°

10

Figure 13 Tomography-based horizontal divergence of the tectonic plates. (a) Observed NUVEL-1 horizontal plate divergence (DeMets et al., 1990). (b)–(g) The surface divergence predicted on the basis of tomography models J362D28 (Antolik et al., 2003), SAW24 (Me´gnin and Romanowicz, 2000) S20A-Iso (Ekstro¨m and Dziewonski, 1998), TX2002 (Grand, 2002), S20RTS (Ritsema et al., 1999), SB4-L18 (Masters et al., 2000), using the corresponding d ln / d ln VS scaling profiles in Figure 12. All shear-velocity heterogeneity models have been projected onto the common parametrization in exptression [102] prior to calculating the predicted plate motions. (h) The ‘mean’ prediction of plate divergence obtained by calculating the statistical sample average of all the predictions (b)–(g). The observed and predicted plate motions shown here are all synthesized from spherical harmonics up to degree and order 32.

836 Implications for Mantle Dynamics and Composition

(a) 1.0 Data (NUVEL-1) J 362D 28 SAW 24B16 S 20A-ISO TX 2002 S 20RTS S B 4_L18 Mean

Amplitude (rad/100 Ma)

0.8

0.6

0.4

0.2

0.0 (b) 1.0

Correlation

0.8

0.6

0.4

0.2

0.0

99% confidence 95% confidence 0

2

4

6

8

10 12 14 16 18 20 22 24 26 28 30 32 Harmonic degree

Figure 14 Spectral comparison of predicted and observed plate divergence. (a) The degree variance, as defined in eqn [3], of the observed and predicted plate divergence fields shown in Figure 13. (b) The degree correlations, as defined in eqn [4].

Table 1

Comparison of observeda and predictedb horizontal plate divergence

Model

NUVELc

J362D28

SAW24

S20A-Iso

TX2002

S20RTS

SB4-L18

Mean

rms correl var. red. (%)

2.8

1.3 0.77 50

1.1 0.74 44

1.5 0.83 60

1.5 0.78 56

1.6 0.80 59

1.2 0.74 45

1.3 0.82 55

a

Data and predictions (shown in Figure 13) are synthesized from spherical harmonics up to degree 32. Predictions employ the tomography model indicated at the top of each column and they use the corresponding velocity–density scalings in Figure 12. c Horizontal divergence of the plates is derived from the NUVEL-1 (De Mets et al., 1990) description of present-day plate velocities. rms, root-mean-square amplitude, expressed here in units of 108 a1. Correl, global cross-correlation; var. red., variance reduction. b

1.23.3.4 Predicted Free-Air Gravity Anomalies Tomography-based mantle flow models of Earth’s nonhydrostatic gravitational potential have almost exclusively been based on analyses of the equivalent

geoid anomalies, beginning with the earliest studies by Richards and Hager (1984), Hager et al. (1985), and continuing to recent studies (e.g., Panasyuk and Hager, 2000). Since the global geoid anomalies are strongly dominated by horizontal wavelengths

Implications for Mantle Dynamics and Composition

corresponding to harmonic degrees , ¼ 2, 3, they will effectively constrain only the longest wavelength components of 3-D mantle structure. To avoid this very-long-wavelength bias, we consider here a representation of the nonhydrostatic geopotential in terms of equivalent global free-air gravity anomalies (see expression [97]). The free-air gravity anomalies contain a more evenly balanced representation of the different horizontal wavelengths in the nonhydrostatic geopotential (i.e., a ‘flatter’ amplitude spectrum). The observed and predicted gravity anomaly fields, calculated on the basis of expressions [97] and [98] and using the kernels in Figures 10(c) and 10(d), are presented in Figure 15. A comparison of the predicted and observed gravity anomalies in the spectral domain is shown in Figure 16. In Table 2 is presented a detailed summary of the global agreement between the predicted and observed gravity anomalies. In contrast to the predicted plate divergence discussed in the preceding section, we note that the amplitudes of the predicted gravity anomaly fields are somewhat less muted (compare Tables 1 and 2). The overall agreement with the gravity data (in terms of variance reduction) is however poorer than the plate divergence fits because of the decreased spatial correlation between the predictions and the data. 1.23.3.5 Predicted Dynamic Surface Topography The term ‘dynamic topography’ is here defined to include all contributions to the topography of Earth’s solid surface which arise from mantle convection. There has been some controversy as to origin and amplitude of Earth’s dynamic topography, with some confusion as to how this dynamic topography should even be defined (e.g., Gurnis, 1990; Forte et al., 1993b, 1993c; Gurnis, 1993). As will be discussed below, the interpretation of dynamic surface topography employed here, namely the topography arising from all density anomalies in the mantle (including those in the lithosphere), is not universally accepted. A detailed review and analysis of the opposing schools of thought in this debate has been presented by Pari (2001). Numerous studies over the past three decades have focused exclusively on the topography of ocean floor, and the general conclusion is that this topography can be almost completely explained in terms of shallow, thermally induced density

837

anomalies which arise from the age-dependent cooling of the oceanic lithosphere (e.g., Parson and Sclater, 1977; Stein and Stein, 1992). The oceanic bathymetry is modeled in terms of isostatic compensation of density anomalies in the oceanic lithosphere, and hence many studies have regarded this form of thermal isostasy as not being a ‘dynamic’ contribution to surface topography. An adequate review of the extensive literature dealing with this interpretation of seafloor topography is well beyond the scope of this chapter, but recent studies of agedependent oceanic bathymetry with extensive references to past analyses may be found in Doin and Fleitout (2000) and in Crosby et al. (2006). Whether or not the thermal isostatic signal in seafloor bathymetry may be regarded as ‘static’ or ‘dynamic’ depends on the temporal variability of the structure of the upper thermal boundary layer (i.e., the lithosphere in oceanic regions) due to time-dependent mantle convection. A purely mechanical interpretation of the instantaneous, present-day isostatic compensation of thermal anomalies in the lithosphere ignores this background time dependence. The lateral temperature variations in the cooling lithosphere are essentially maintained by a balance between vertical heat conduction and horizontal heat advection (e.g., Jarvis and Peltier, 1982) and the latter is clearly a dynamic effect. In the absence of convection, the upper thermal boundary layer and the corresponding variations in oceanic bathymetry will vanish. The importance of time-dependent changes in the structure of the lithosphere, and hence the corresponding changes in the contribution to surface topography, is highlighted by the significant changes in plate tectonic velocities and plate geometries in the Cenozoic (e.g., Lithgow-Bertelloni and Richards, 1998), which suggest significant departures from steady-state conditions. In the absence of convection, the only contribution to surface topography arises from the isostatic compensation of lateral variations in crustal thickness and crustal density. Earth’s present-day topography is therefore the superposition of the crustal isostatic topography and the dynamic topography – as defined here – due to density anomalies in the mantle which are maintained by mantle convection (Forte et al., 1993b). The crustal isostatic topography is large and it explains most of the observed present-day topography on Earth (e.g., Forte and Perry, 2000; Pari and Peltier, 2000). The dynamic topography is obtained by subtracting the crustal isostatic topography from the observed topography and this therefore requires an accurate model of crustal structure. This

838 Implications for Mantle Dynamics and Composition

EGM 96 gravity anomalies (L = 2–20)

(a) 90°

0°

30°

60°

90° 120° 150° 180° –150°–120° –90° –60° –30°

0°

(b) 90°

J 362D28 gravity anomalies (L = 2–20) 0°

30°

60°

0°

30°

60°

0°

30°

60°

0°

30°

60°

60°

60°

30°

30°

0°

0°

–30°

–30°

–60°

–60°

–90°

–90°

SAW 24B16 gravity anomalies (L = 2–20)

(c)

0°

30°

60°

90° 120° 150° 180° –150°–120° –90° –60° –30°

0°

(d)

90°

90°

60°

60°

30°

30°

0°

0°

–30°

–30°

–60°

–60°

–90°

–90°

TX 2002 gravity anomalies (L = 2–20)

(e)

0°

30°

60°

90° 120° 150° 180° –150°–120° –90° –60° –30°

0°

(f)

90°

90°

60°

60°

30°

30°

0°

0°

–30°

–30°

–60°

–60°

–90°

–90°

SB4_L18 gravity anomalies (L = 2–20)

(g)

0°

30°

60°

90° 120° 150° 180° –150°–120° –90° –60° –30°

0°

(h)

90°

90°

60°

60°

30°

30°

0°

0°

–30°

–30°

–60°

–60°

–90°

–90°

–40 –30 –20 –10

90° 120° 150° 180° –150°–120° –90° –60° –30°

0°

S 20A-Iso gravity anomalies (L = 2–20) 90° 120° 150° 180° –150°–120° –90° –60° –30°

0°

S 20RTS gravity anomalies (L = 2–20) 90° 120° 150° 180° –150°–120° –90° –60° –30°

0°

‘Mean’ gravity anomalies (L = 2–20)

0 10 20 mGal

30

90° 120° 150° 180° –150°–120° –90° –60° –30°

0°

40

Figure 15 Tomography-based global free-air gravity anomalies. (a) Observed EGM96 anomalous free-air gravity field (Lemoine et al., 1998). (b)–(g) The free-air anomalies predicted on the basis of tomography models J362D28 (Antolik et al., 2003), SAW24 (Me´gnin and Romanowicz 2000), S0A-Iso (Ekstro¨m and Dziewonski, 1998), TX2002 (Grand 2002), S20RTS (Ritsema et al. 1999), SB4_L18 (Masters et al. 2000), using the corresponding d ln/d ln VS scaling profiles in Figure 12. All shear-velocity heterogeneity models have been projected onto the common parameterization in expression [102] prior to calculating the gravity anomalies. (h) The ‘mean’ free-air gravity anomaly field obtained by calculating the statistical sample average of all the predictions (b)–(g). The observed and predicted gravity anomalies shown here are all synthesized from spherical harmonics up to degree and order 20.

Implications for Mantle Dynamics and Composition

(a)

6.0

Data (EGM96 ) J 362D 28 SAW 24B16 S 20A-ISO TX 2002 S 20RTS S B 4_L18 Mean

5.0 Amplitude (mGal)

839

4.0 3.0 2.0 1.0 0.0

(b)

1 99% confidence 95% confidence

0.8

Correlation

0.6 0.4 0.2 0 –0.2 –0.4 0

2

4

6

8 10 12 14 Harmonic degree

16

18

20

Figure 16 Spectral comparison of predicted and observed free-air gravity anomalies. (a) The degree variance, as defined in eqn [3], of the observed and predicted gravity anomaly fields shown in Figure 13. (b) The degree correlations, as defined in eqn [4].

Table 2

Comparison of observeda and predictedb free-air gravity anomalies

Model

EGM96c

J362D28

SAW24

S20A-Iso

TX2002

S20RTS

SB4_L18

Mean

rms correl var. red. (%)

14.8

10.4 0.55 28

11.2 0.62 37

10.6 0.54 26

11.4 0.64 39

10.5 0.58 33

8.9 0.42 15

9.3 0.64 41

a

Data and predictions (shown in Figure 15) are synthesized from spherical harmonics up to degree 20. Predictions employ the tomography model indicated at the top of each column and they use the corresponding velocity–density scalings in Figure 12. c The global free-air gravity anomalies are derived from the EGM96 geopotential model (Lemoine et al., 1998). rms, root-mean-square amplitude, expressed here in units of mGal. Correl, global cross-correlation; var. red., variance reduction. b

requirement presents a major challenge because of the uneven and incomplete seismic sampling global crustal thickness and also because of the significant uncertainties in constraining the density of the crust. Currently, the most complete compilation of global crustal heterogeneity is model CRUST2.0 (Bassin

et al., 2000). This crustal model was employed in retrieving the dynamic surface topography presented above in Figure 2(b). The observed and predicted dynamic topography fields, calculated on the basis of expressions [99] and using the kernels in Figures 10(e) and 10(f), are

840 Implications for Mantle Dynamics and Composition

presented in Figure 17. A quantitative spectral comparison of the predicted and observed dynamic topography is shown in Figure 18. In Table 3 is presented a detailed summary of the match between the predicted and observed topography fields. The topography kernels (Figures 10(e) and 10(f)) show that near-surface density anomalies provide the strongest contributions to the predicted surface topography (Figure 17) and that they will be in near-isostatic equilibrium, where perfect isostatic compensation corresponds to a kernel value of 1. We note that in the case of density anomalies which effectively see a no-slip surface boundary (given by expression [92]), the condition of near-isostasy extends to depths of about 200 km (see Figure 10(f )). In view of the debate concerning the magnitude of the thermal isostatic contributions to surface topography in oceanic regions, it is helpful to explore the relative importance of shallow (i.e., lithospheric) and deep-mantle (i.e., sublithospheric) buoyancy using the seismic tomography models. In this regard, it is useful to first evaluate the extent to which the tomography-based topography predictions incorporate the age-dependent cooling of the oceanic lithosphere. Tomography models employing seismic traveltime and/or surface wave data place strong constraints on lithospheric mantle heterogeneity and they contain a clear signature of the oceanic cooling history (e.g., Woodward and Masters, 1991; Su et al., 1992; Trampert and Woodhouse, 1996; Ritzwoller et al., 2004). Here we consider the global tomography models SAW24 (Me´gnin and Romanowicz, 2000) and TX2002 (Grand, 2002), employed in Figures 17(c) and 17(e) respectively. The age-dependent lateral temperature variations in the oceanic lithosphere were calculated using a simple halfspace cooling model with the same thermal parameters as the GDH1 plate-cooling model of Stein and Stein (1992). The seafloor age was determined using the Muller et al. (1997) digital isochrons. The global ocean basins were then sampled on a 1  1 grid and for each grid cell the corresponding values of ocean age, shear velocity perturbation dVS/VS, and temperature perturbation T were determined. These samples were then averaged into 2-million-year age bins and the resulting relationship between seismic shear velocity and temperature is plotted in Figure 19. A linear regression analysis (solid lines in Figure 19) yields an excellent fit to the binned VST variation and the slopes of the regression lines provide estimates of the effective temperature

derivative of S-velocity in the depth range 0–100 km. These estimated derivatives are 8.7  105 and 9.4  105 K1 for models SAW24 and TX2002, respectively. These effective thermal derivatives agree well with independent mineral physics values determined by Stixrude and Lithgow-Bertelloni (2005a) which range between 8  105 and 9  105 K1. These results confirm that the global tomography models considered in this study successfully resolve the pattern and amplitude of agedependent cooling of the oceanic lithosphere. The relative importance of surface topography due to lithospheric versus sublithospheric density anomalies is illustrated in Figure 20. As expected, the topography contributions in oceanic regions from density anomalies in the upper 200 km of the mantle are dominated by the age-dependent cooling of the oceanic lithosphere. Depressions in continental regions exceed observational estimates (Figure 17(a)) because the use of a purely depth-dependent d ln /d ln VS scaling factor (Figure 12) does not resolve the intrinsic chemical buoyancy in the subcontinental mantle. This intrinsic buoyancy can be included in geodynamic models which allow for lateral variations of d ln /d ln VS in the shallow mantle (e.g., Forte and Perry, 2000; Deschamps et al., 2001). The topography generated by sublithospheric buoyancy (Figure 18(b)), with rms amplitude equal to 410 m, is comparable to that generated by lithospheric density anomalies (Figure 20(a)) which has an rms amplitude equal to 580 m. These sublithospheric contributions to topography in the central Pacific and Atlantic Oceans, and in southern Africa, are surface expressions of deep-seated buoyancy with sources in the lower mantle (Figure 3). The depressions below the eastern and western margins of the Pacific (Figure 20(b)) may be interpreted in terms of present-day and Cenozoic subduction history (e.g., Lithgow-Bertelloni and Gurnis, 1997). We note from Figure 17 that the amplitude of the predicted dynamic topography is somewhat larger than the crust-corrected estimate of the observed dynamic topography. Table 3 shows that the rms amplitude of the mean topography prediction is 25% larger than the observed topography. This discrepancy is in marked contrast to that of the predicted gravity and plate divergence fields (Tables 1 and 2) which are, on average, 40–50% smaller than the corresponding observed fields. The mismatched amplitudes of the predicted surface topography has motivated previous efforts (e.g., Le Stunff and Ricard 1995, 1997) to introduce additional buoyancy forces in the mantle arising,

Implications for Mantle Dynamics and Composition

Dynamic topography, CRUST2.0 (L = 1–20)

(a) 0°

30° 60° 90° 120° 150° 180°–150°–120°–90° –60° –30° 0°

J 36D 28 topography (L = 1–20) (b) 0° 30° 60° 90° 120° 150° 180°–150°–120°–90° –60° –30° 0°

90°

90°

60°

60°

30°

30°

0°

0°

–30°

–30°

–60°

–60° –90°

–90°

(c)

SAW 24B16 topography (L = 1–20) 0°

30° 60° 90° 120° 150° 180°–150°–120°–90° –60° –30° 0°

S 20A-Iso topography (L = 1–20) (d) 0° 30° 60° 90° 120° 150° 180°–150°–120°–90° –60° –30° 0°

90°

90°

60°

60°

30°

30°

0°

0°

–30°

–30°

–60°

–60°

–90°

–90°

TX 2002 topography (L = 1–20) (d)

0°

30° 60° 90° 120° 150° 180°–150°–120°–90°

S 20RTS topography (L = 1–20) –60°

–30° 0°

(f)

90°

90°

60°

60°

30°

30°

0°

0°

–30°

–30°

–60°

–60°

–90°

0°

30° 60° 90° 120° 150° 180°–150°–120°–90° –60° –30° 0°

–90°

SB 4_L18 topography (L = 1–20) (g)

841

0°

30° 60° 90° 120° 150° 180°–150°–120°–90° –60° –30° 0°

‘Mean’ dynamic topography (L = 1–20) (h) 0° 30° 60° 90° 120° 150° 180°–150°–120°–90° –60° –30° 0°

90°

90°

60°

60°

30°

30°

0°

0°

–30°

–30°

–60°

–60°

–90°

–90°

–2.0 –1.5 –1.0 –0.5 0.0 0.5 1.0 1.5 2.0 km

Figure 17 Tomography-based global dynamic surface topography. (a) Dynamic surface topography estimated by removing isostatically compensated crustal heterogeneity – described by model CRUST2.0 (Bassin et al., 2000) – from the observed present-day surface topography. (b)–(g) The dynamic surface topography predicted on the basis of tomography models J362D28 (Antolik et al., 2003), SAW24 (Me´gnin and Romanowicz, 2000), S20A-Iso (Ekstro¨m and Dziewonski, 1998), TX2002 (Grand 2002), S20RTS (Ritsema et al., 1999), SB4_L18 (Masters et al., 2000), using the corresponding d ln /d ln VS scaling profiles in Figure 12. All shear-velocity heterogeneity models have been projected onto the common parametrization in expression [102] prior to calculating the dynamic topography. (h) The ‘mean’ dynamic surface topography field obtained by calculating the statistical sample average of all the predictions (b)–(g). The crust-corrected and predicted dynamic topography shown here are all synthesized from spherical harmonics up to degree and order 20.

842 Implications for Mantle Dynamics and Composition

(a)

Data (CRUST2.0) J 362D 28 SAW 24B16 S 20A-Iso TX 2002 S 20RTS S B 4_L18 Mean

Amplitude (km)

0.5 0.4 0.3 0.2 0.1 0.0 (b)

1 99% confidence 95% confidence

0.8

Correlation

0.6 0.4 0.2 0 –0.2

0

2

4

6

8 10 12 14 Harmonic degree

16

18

20

Figure 18 Spectral comparison of predicted and observed (crust-corrected) dynamic surface topography. (a) The degree variance, as defined in eqn [3], of the observed and predicted topography fields shown in Figure 15. (b) The degree correlations, as defined in eqn (4). Table 3

Comparison of observeda and predictedb dynamic surface topography

Model

CRUSTc

J362D28

SAW24

S20A-Iso

TX2002

S20RTS

SB4_L18

Mean

rms correl var. red. (%)

696

908 0.75 25

893 0.76 30

1002 0.80 22

860 0.78 41

972 0.78 24

849 0.79 44

873 0.87 47

a

Data and predictions (shown in Figure 17) are synthesized from spherical harmonics up to degree 20. Predictions employ the tomography model indicated at the top of each column and they use the corresponding velocity–density scalings in Figure 12. c The global dynamic surface topography obtained by subtracting CRUST2.0 (Bassin et al., 2000) isostatic crust from the observed topography. rms, root-mean-square amplitude, expressed here in units of m. correl, global cross-correlation; var. red., variance reduction. b

for example, from undulations of the 670 km seismic discontinuity (e.g., Thoraval et al., 1995; Forte and Woodward, 1997b) which internally compensate the excess surface topography. In assessing the significance of this discrepancy, in terms of mantle dynamics, it is important to recognize

that the crustal correction employed in Figure 17(a) contains significant uncertainties which may in part be responsible for the mismatched amplitudes. Perhaps more important is the question of the validity of a purely depth-dependent velocity–density scaling (Figure 12) employed in the topography predictions.

Implications for Mantle Dynamics and Composition

6 4

Model TX 2002 Model SAW 24

dVS (%)

2 0 –2 –4 –6 –400 –300 –200 –100 0 100 200 300 400 ΔT (K)

500

Figure 19 Comparison of seismic shear velocity and temperature anomalies in the oceanic lithosphere. Seismic shear velocity anomalies, vertically averaged in the depth interval 0–100 km, are derived from tomography models SAW24 (Me´gnin and Romanowicz, 2000) and TX2002 (Grand, 2002). The S-velocity anomalies, measured relative to the global oceanic mean value, were averaged into 2 Ma age bins using the global ocean-age compilation of Muller et al. (1997). The age-dependent temperature variations in the oceanic lithosphere were determined using a simple halfspace cooling model with the same thermal parameters as the GDH1 plate model of Stein and Stein (1992). These temperature variations were vertically averaged down to 100 km depth and subsequently averaged into 2 Ma age bins. The temperature anomalies are measured relative to the global oceanic mean value. The vertical and horizontal error bars represent 1 standard deviation relative to the average shear-velocity and temperature, respectively, in each age bin. The solid blue and red lines are the best-fitting linear regression lines for models SAW24 and TX2002, respectively.

Such depth dependence, which may be justified if one assumes that thermal effects on mantle heterogeneity dominate (e.g., Karato and Karki, 2001), is not a good approximation in the presence of large lateral variations in chemical composition. Such compositional heterogeneity is expected to be important in the shallow subcontinental mantle and it will oppose the local thermal buoyancy, implying significant reductions in the amplitude of continental surface topography (e.g., Forte and Perry, 2000). To understand why the gravity (and plate motion) predictions, in contrast to the surface topography, have substantially smaller amplitudes than the observed fields, we must first appreciate that the global seismic data have difficulty resolving heterogeneity in the deep mantle, particularly in the

843

Southern Hemisphere (e.g., Forte et al., 1994; Simmons et al., 2006). This difficulty does not apply to the shallow mantle where, as shown in Figure 19, the tomography models are able to successfully resolve the cooling oceanic lithosphere. The predicted surface topography is dominated by the density anomalies located in the upper mantle, whereas the gravity predictions contain significant contributions from lower-mantle density anomalies (Figure 10). Joint inversions of global seismic and convection-related data sets have shown that increased amplitudes of the predicted gravity anomalies and substantially improved fits to the global gravity data can be achieved by modifying deep-mantle heterogeneity in a way that is consistent with both the seismic and geodynamic constraints (e.g., Forte et al., 1994; Simmons et al., 2006). Such modifications to lowermantle heterogeneity have little impact on the amplitude of the predicted surface topography. These basic observations suggest that current mismatches between the amplitudes of the predicted and observed dynamic topography are not more significant, or fundamentally different, from the mismatch between the other predictions (i.e., gravity and plate motions) and their corresponding observables. 1.23.3.6

Predicted CMB Topography

Tomography-based studies of the flow-induced deformations of the CMB has been rather limited (e.g., Hager et al., 1985; Forte et al., 1993a, 1995) relative to the much greater number of studies which have focused on the surface geoid or topography. This is in large part a consequence of the uncertain and often contrasting inferences of the CMB topography which have been derived from inversions of seismic phases which are sensitive to this topography (e.g., Morelli and Dziewonski, 1987; Rodgers and Wahr, 1993, Obayashi and Fukao, 1997; Boschi and Dziewonski, 2000; Garcia and Souriau, 2000; Sze and van der Hilst, 2003). A detailed discussion of the seismological uncertainties and difficulties in obtaining reliable global images of the CMB topography may be found in Garcia and Souriau (2000). The most accurate constraint on the topographic undulations of the CMB are derived from spacegeodetic analyses of Earth’s nutations. A brief review of these constraints was presented in Section 1.23.2.1. Unfortunately, these geodetic analyses can only constrain the excess flattening or dynamic ellipticity of the CMB – albeit with high precision – and this corresponds to only one coefficient (, ¼ 2, m ¼ 0) in

844 Implications for Mantle Dynamics and Composition

(a) 90°

0°

30°

60°

90°

120°

150° 180° –150° –120° –90° –60°

–30°

0°

0°

30°

60°

90°

120° 150° 180° –150° –120° –90° –60° –30°

0°

60° 30° 0° –30° –60° –90° (b) 90° 60° 30° 0° –30° –60° –90°

–2.0 –1.5 –1.0 –0.5 0.0 0.5 1.0 1.5 2.0 km Figure 20 Shallow- versus deep-mantle contributions to surface topography. (a) Surface topography predicted on the basis of tomography-derived density anomalies in the upper 200 km of the mantle. (b) Surface topography predicted on the basis of all density anomalies below 200 km depth. In both cases, the density anomalies are derived from tomography model TX2002 (Grand 2002) and using the corresponding d ln /d ln VS scaling profiles in Figure 12. The summation of the topography fields in maps (a) and (b) yields the total field shown above in Figure 17(e). All fields are synthesized from spherical harmonics up to degree and order 20.

a full spherical harmonic expansion (expression [1]) of the spatially varying CMB topography field. This constraint on excess CMB flattening is included in all the Occam-inferred d ln /d ln VS scaling profiles in Figure 12. For the purpose of comparison, we consider a relatively recent seismological inference of the CMB topography (Figure 21(a)) obtained by Boschi and Dziewonski (2000) (henceforth referred to as BD2000) on the basis of seismic PcP and PKP traveltime delays. The flow-induced CMB topography predicted on the basis of the Occam-

inferred d ln /d ln VS profiles (Figure 12), and using the kernels in Figures 10(g) and 10(h), is presented in Figure 21. All predictions of CMB topography shown here exactly reproduce the geodetically constrained 400 m excess CMB ellipticity. A quantitative spectral analysis (Figure 22) of the predicted CMB topography fields shows that they are strongly dominated by the horizontal wavelengths corresponding to harmonic degrees , ¼ 2, 3. In Table 4 is presented a detailed summary of the match between the predicted and observed CMB topography fields.

Implications for Mantle Dynamics and Composition

(a)

BD2000 CMB topography (L = 1–20)

0° 30° 60° 90° 120° 150° 180°–150°–120°–90° –60° –30° 0°

90°

90°

60°

60°

30°

30°

0°

0°

–30°

–30°

–60°

–60°

–90°

–90°

(c)

SAW 24B16 CMB topography (L = 1–20) 0°

(d) 90°

60°

60°

30°

30°

0°

0°

–30°

–30°

–60°

–60°

–90°

(d)

TX 2002 CMB topography (L = 1–20) 0°

S 20A-Iso CMB topography (L = 1–20) 0° 30° 60° 90° 120° 150° 180°–150°–120°–90° –60° –30° 0°

30° 60° 90° 120° 150° 180°–150°–120°–90° –60° –30° 0°

90°

–90° (f)

30° 60° 90° 120° 150° 180°–150°–120°–90° –60° –30° 0°

S 20RTS CMB topography (L = 1–20) 0° 30° 60° 90° 120° 150° 180°–150°–120°–90° –60° –30° 0°

90°

90°

60°

60°

30°

30°

0°

0°

–30°

–30°

–60°

–60° –90°

–90°

SB 4_L14 CMB topography (L = 1–20)

(g) 0°

(h)

30° 60° 90° 120° 150° 180°–150°–120°–90° –60° –30° 0°

‘Mean’ CMB topography (L = 1–20) 0°

90°

90°

60°

60°

30°

30°

0°

0°

–30°

–30°

–60°

–60°

–90°

J 362D 28 CMB topography (L = 1–20)

(b)

0° 30° 60° 90° 120° 150° 180°–150°–120°–90° –60° –30° 0°

845

30° 60° 90° 120° 150° 180°–150°–120°–90° –60° –30° 0°

–90° –4.5 –3.0 –1.5 0.0

1.5

3.0

4.5

–1.5 –1.0 –0.5 0.0 km

0.5

1.0

1.5

Figure 21 Tomography-based global CMB topography. (a) Seismically inferred CMB topography derived by Boschi and Dziewonski (2000). The amplitudes correspond to the top values ( 4.5 km) in the bottom scale bar. (b)–(g) The CMB topography predicted on the basis of tomography models J362D28 (Antolik et al., 2003), SAW24 (Me´gnin and Romanowicz, 2000), S20A-Iso (Ekstro¨m and Dziewonski, 1998), TX2002 (Grand, 2002), S20RTS (Ritsema et al., 1999), SB4_L18 (Masters et al., 2000), using the corresponding d ln /d ln VS scaling profiles in Figure 12. All shear-velocity heterogeneity models have been projected onto the common parametrization in expression [102] prior to calculating the CMB topography. (h) The ‘mean’ CMB topography field obtained by calculating the statistical sample average of all the predictions (b)–(g). The amplitudes for (b)–(h) correspond to the bottom values ( 1.5 km) in the bottom scale bar. The observed and predicted CMB topography shown here are all synthesized from spherical harmonics up to degree and order 20.

846 Implications for Mantle Dynamics and Composition

(a) ‘Data’ (BD2000) J 362D 28 SAW 24B16 S 20A-Iso TX 2002 S 20RTS SB4_L18 Mean

Amplitude (km)

1.5

1.0

0.5

0.0 (b)

1 99% confidence 95% confidence

0.8

Correlation

0.6 0.4 0.2 0 –0.2 –0.4 0

2

4

6

8 10 12 14 Harmonic degree

16

18

20

Figure 22 Spectral comparison of predicted and observed CMB topography. (a) The degree variance, as defined in eqn [3], of the observed and predicted CMB topography fields shown in Figure 19. (b) The degree correlations, as defined in eqn [4].

The overall spatial correlation between the predicted maps of flow-induced CMB topography and the BD2002 model (Figure 21) is quite good; however, it is quite clear that there is a significant difference in the overall amplitude. The CMB topography obtained on the basis of the tomographybased flow calculations is approximately a factor of 3–4 times smaller than that in BD2002. Although, as noted in Section 1.23.3.3, there likely exists an inherent damping in the flow predictions, it is not possible to increase the amplitude of the predicted Table 4

CMB topography (by increasing, for example, the magnitude of the d ln /d ln VS scaling) without at the same time significantly degrading the fit to the other geodynamic observables. The 1.5 km scale undulations in the predicted, flow-induced CMB topography are in the range estimated by Garcia and Souriau (2000) for horizontal wavelengths in excess of 1200 km.

1.23.4 Tomography-Based Geodynamic Inferences of Compositional Heterogeneity Significant lateral variations of chemical composition in the mantle are expected in view of the extensive and repeated episodes of partial melting and subsequent chemical differentiation which have acted throughout the history of our planet. The creation and growth of continental crust is obviously a major manifestation of this chemical differentiation (e.g., Taylor and McLennan, 1995) which has led to large-amplitude compositional heterogeneity in the subcontinental mantle (e.g., Jordan, 1981). In the lower mantle, it is postulated (e.g., Hofmann, 1997) that distinct reservoirs of primordial or primitive composition, which have interacted over geologic time with recycled crustal material entrained in sinking lithospheric slabs, provide an explanation for the highly variable isotopic anomalies which have been observed in ocean island basalts. These island basalts, as well as mid-ocean ridge basalts, have consistently provided evidence of large-scale chemical heterogeneity in the lower mantle (e.g., Hart, 1988). Efforts to resolve and incorporate such large-scale compositional heterogeneity in tomography-based mantle flow models have been relatively recent (e.g., Forte et al., 1995b; Forte and Perry, 2000; Pari and Peltier, 2000; Forte and Mitrovica, 2001). These modeling efforts have depended on significant

Comparison of seismically observeda and predictedb CMB topography

Model

BD2000c

J362D28

SAW24

S20A-Iso

TX2002

S20RTS

SB4_L18

Mean

rms correl var. red. (%)

2.3

0.58 0.56 22

0.42 0.56 18

0.76 0.64 31

0.53 0.59 22

0.61 0.65 27

0.51 0.55 20

0.55 0.62 24

a

Seismic observations and predictions (shown in Figure 21) are synthesized from spherical harmonics up to degree 20. Predictions employ the tomography model indicated at the top of each column and they use the corresponding velocity–density scalings in Figure 12. c CMB topography derived by Boschi and Dziewonski (2000) from global tomographic inversion of seismic delay times. rms, root-mean-square amplitude, expressed here in units of km. correl, global cross-correlation; var. red., variance reduction. b

Implications for Mantle Dynamics and Composition

improvements in the resolution and reliability of global tomography models, largely driven by a major increase in the number of high-quality digital broadband seismic stations deployed worldwide. The most important advance, from the perspective of understanding compositional heterogeneity in the mantle, was the advent of joint inversions for lateral variations in seismic shear velocity dVS/VS and bulksound or acoustic velocity dV/V (Su and Dziewonski, 1997; Kennett et al., 1998; Masters et al., 2000; Antolik et al., 2003). These models have revealed that the perturbations in seismic shear and bulk-sound velocity are anti-correlated in several regions in the lower mantle and, as argued by Stacey (1998), it is difficult to reconcile this anticorrelation with pure thermal heterogeneity. 1.23.4.1

Constraints from Mineral Physics

A quantitative understanding of the implications of the joint VS – V seismic maps in terms of compositional heterogeneity requires knowledge of the thermodynamic properties of mantle minerals. Detailed expositions of the underlying thermodynamic theory and experimental data sets may be found in numerous publications whose discussion is well beyond the scope of the work presented here and the interested reader may find the following studies to be helpful: Stacey (1998), Jackson (1998, 2001), Karato and Karki (2001), Trampert et al., (2001), Cammarano et al. (2003), Mattern et al. (2005), and Stixrude and Lithgow-Bertelloni (2005b). The following discussion is limited to a review of the mineral physics results most relevant to the tomography-based inferences of mantle chemical heterogeneity presented below. In principle, tomography-based inferences of thermal and compositional anomalies can be derived on the basis of temperature and compositional derivatives of density and seismic velocity, assuming a particular background (i.e., reference) composition and temperature for the mantle. A previous analysis by Forte and Mitrovica (2001), henceforth referred to as FM2001, modeled the lower mantle in terms of three components: MgO, FeO, and SiO2, assuming a pyrolite composition for the mantle. This is a simplification which ignores the possible importance of less abundant constituents such as CaSiO3 perovskite (about 5% molar abundance) and solid solution of alumina (Al2O3) in ferromagnesium perovskite. Recent studies have suggested these minor components can have a non-negligible effect on density and seismic-wave velocity (e.g., Jackson, 2001), such that

847

Ca-perovskite increases the shear velocity (Stixrude and Lithgow-Bertelloni, 2005b) and Al-perovskite decreases the shear velocity (Jackson, 2005). Recent efforts to constrain the composition and temperature of the mantle using all five components (MgO–FeO– CaO–Al2O3–SiO2) have revealed considerable nonuniqueness and large tradeoffs (Mattern et al., 2005). Therefore, for the purpose of this discussion, a simplified three-component (MgO–FeO–SiO2) mantle model will again be assumed. The perturbations in density and seismic velocity produced by temperature T and compositional anomalies are modeled as follows: d ln M ¼

q ln M q ln M q ln M dT þ dXFe þ dXPv ½109 qT qXFe qXPv

where M is a generic variable used to denote either density , or shear velocity VS, or bulk-sound velocity V. Following FM2001, the compositional anomalies are parametrized in terms of perturbations in the molar fraction of iron, XFe ¼ [FeO]/ ([FeO] þ [MgO]), and the molar fraction of perovskite, XPv ¼ [Pv]/([Pv] þ [Mw]). The molar fraction of perovskite is equivalent to the silica ratio: XPv ¼ [SiO2]/([FeO] þ [MgO]). Estimates of the thermal and compositional derivatives which are required in expression [109] have been estimated by a number of authors, using different assumptions and methods, and a summary of some of these results is presented in Table 5. The differences between the various estimates presented in this table show that the uncertainties in these derivatives are significant. The most recent determinations by Stixrude and Lithgow-Bertelloni (2005b) also include formal estimates of uncertainties and they are of the same order as the differences in the Table 5. These uncertainties should be borne in mind when considering inferences of thermochemical heterogeneity in the mantle. FM2001 assembled the three equations described by expression [109], involving , VS, and V, into the following linear system: 0

d ln VS

1

0

c11 c12 c13

10

dT

1

B C B CB C B d ln V C ¼ B c21 c22 c23 CB dXPv C @ A @ A@ A d ln  c31 c32 c33 dXFe

½110

in which the matrix elements cij are the thermal or compositional derivatives involved in expression [109]. To extract the maximum possible information from the joint VS–V tomography models, FM2001 introduced the concept of effective thermal and

848

Implications for Mantle Dynamics and Composition

Table 5

Thermal and compositional derivatives Thermal 105  q/qT(K1)

Elastic property Depth ¼ 1970 km (85 GPa) ln  ln VS ln V Depth ¼ 2740 km (127 GPa) ln  ln VS ln V

Compositional

Anharmonic

Anelastic

10  q/qXFe

100  q/qXPv

1.3[1.2](1.3){1.4} 5.0[4.4[(4.7){3.9}h3.5i 0.8[1.5](0.6){0.9}

2.5 0

þ3.2[þ3.0] 2.2[2.5]h1.7i 1.6[1.5]

þ1.1[þ2.8] þ4.9[þ3.5] þ5.5[þ3.4]

1.0[0.8](1.1) 4.7[4.0](4.1)h2.9i 0.7[1.2](0.4)

2.4 0

þ3.2[þ3.0] 2.2[2.5]h1.4i 1.6[1.5]

þ0.4[þ2.2] þ4.5[þ1.9] þ4.8[þ2.0]

All values not enclosed by brackets have been obtained using the approximations and data described in detail by Forte and Mitrovica (2001). Values enclosed in square brackets have been obtained on the basis of the finite-strain and thermoelastic calculations carried out by Forte et al. (2002). Round brackets (parentheses) enclose values estimated by Stacey (1998), while the curly brackets enclose the determinations obtained by Oganov et al. (2001b). Values enclosed in angular brackets are obtained on the basis of results for perovskite in Stixrude and Lithgow-Bertelloni (2005b).

compositional anomalies, dTeff and dXeff, defined as follows: dTeff X dT þ

c12 c23 – c13 c22 dXFe c12 c21 – c11 c22

½111

c21 c13 – c11 c23 dXFe c12 c21 – c11 c22

½112

dXeff X dXPv þ

On the basis of these effective anomalies, the linear system [110] can be reduced to the following simpler 2  2 matrix equation: d ln VS d ln V

! ¼

c11 c12 c21 c22

!

dTeff

!

dXeff

½113

which is supplemented by the following relation, describing the dependence of density anomalies on effective thermal and compositional anomalies and on iron anomalies: d ln  ¼ c31 dTeff þ c32 d Xeff –

 d XFe ½114 c12 c21 – c11 c22

in which  is the determinant of the matrix in eqn [110]. The utility in working with effective thermal and compositional anomalies is that the linear system [113] may be simply inverted, thereby providing a means for directly distilling joint shear/bulk-sound tomography models into maps of mantle thermochemical structure. FM2001 found that the effective compositional heterogeneity dXeff is strongly correlated to the bulk-sound velocity anomalies. This may be understood by inverting system (113), using the FM2001 derivatives in Table 5, which yields dXeff ¼ – 2:3 d ln VS þ 23:0 d ln V

½115

This expression shows that dV provides the dominant contribution to dXeff. The dominant compositional signature of the bulk-sound velocity anomalies is traced to the relatively weak temperature sensitivity of bulk-sound velocity compared to shear velocity (Table 5). This implies that bulk-sound velocity should be especially sensitive to any chemical anomalies in the mantle. In contrast, FM2001 found that the effective thermal anomalies dTeff are strongly correlated to shear velocity anomalies. This may be understood by inverting system [113]:   dTeff ¼ – 15:5 d ln VS þ 14:5 d ln V  103

½116

in which the FM2001 derivatives in Table 5 were used. As will be shown below, the amplitude of the lower-mantle dV anomalies is substantially less than that of the dVS anomalies. This then leads to the prevailing correlation between dVS and dTeff in the deep mantle, as might have been expected by consideration of the contrasting temperature sensitivities of VS and V (see Table 5). Although joint bulk-sound/shear tomography models provide very useful maps of effective temperature and compositional anomalies, they do not provide sufficient information for separating the respective contributions to these effective parametres from temperature, iron, and perovskite (or equivalent silica) variations (expressions [111] and [112]). This separation requires additional independent constraints which can be obtained from geodynamically inferred density anomalies. Previous efforts to constrain lower-mantle density anomalies with tomography-based geodynamic

Implications for Mantle Dynamics and Composition

(a) 2.50

SB 10_L18, VS SB 10_L18, Vφ J 362D 28, VS J 362D 28, Vφ

2.00

δV/V (%)

models (e.g., Forte and Mitrovica, 2001) assumed a purely depth-dependent variation in d ln /d ln VS. As discussed in Section 1.23.3.5, this assumption may not be a good approximation if the lateral variations in mantle composition are sufficiently large. In the next two sections, we therefore consider two approaches which allow lateral variations in the effective d ln /d ln VS throughout the mantle.

849

1.50 1.00 0.50 0.

As noted above, the effective compositional heterogeneity in the lower mantle is correlated to the bulksound velocity anomalies whereas seismic shear velocity anomalies provide a reliable mapping of lateral temperature variations, even in the presence of compositional heterogeneity. We can use these findings to motivate the following heuristic representation for mantle density anomalies which incorporates both thermal and chemical contributions: d ¼ 



   d ln  dVS d ln  dV þ d lnVs VS d lnV V

½117

in which the first term on the right, involving dVS/VS, effectively represents the thermal contribution to mantle heterogeneity while the second term (dV/ V) represents the effective compositional heterogeneity. An important feature of expression [117] is that the relationship between density anomalies and seismic shear/bulk anomalies will vary horizontally because the seismic anomalies dVS/VS and dV/V will not necessarily be correlated in the mantle. Here we consider the joint VS–V tomography model SB10_L18 of Masters et al. (2000) and model J362D28 of Antolik et al. (2003). A quantitative summary of the rms amplitudes and spatial correlations of these models is presented in Figure 23. We first note that the amplitudes of both shear and bulksound velocity anomalies agree fairly well across the mantle, with the exception of the transitionzone region (where the shear velocity anomalies in model J362D28 are much larger) and the shallow mantle (where the bulk sound velocity anomalies in model SB10_L18 are much larger). While the shearvelocity anomalies in models J362D28 and SB10_L18 are fairly well correlated, the bulk-sound velocity anomalies show very little correlation and this underlines the difficulty in adequately resolving V anomalies in the mantle. Model J362D28 shows a

(b)

1

SB 10 –J 362, VS SB 10 –J 362, Vφ VS-Vφ, SB 10_L18 VS-Vφ, J 362D 28

0.5 Correlation

1.23.4.2 Compositional Density Anomalies Inferred from Joint Shear- and Bulk-Sound Tomography

0 –0.5 –1 0

500

1000 1500 2000 Depth (km)

2500

3000

Figure 23 Seismic shear and bulk sound velocity anomalies in the mantle. (a) The depth variation of the rms amplitudes of the relative perturbations of seismic shear velocity (dVS/VS) (red lines) and bulk-sound velocity (dV/V) (black lines) in models SB10_L18 of Masters et al. (2000) and J362D28 of Antolik et al. (2003) are plotted as solid and dashed lines, respectively. (b) The global cross-correlation between the shear-velocity anomalies in models SB10_L18 and J362D28 are plotted as a red solid line and for the bulksound velocity anomalies it is plotted as a black solid line. The global cross-correlation between dVS/VS and dV/V in model SB10_L18 is plotted as a blue solid line and for model J362D28 it is plotted as a blue dashed line.

significant anti-correlation between dVS/VS and dV/V throughout the mantle, whereas model SB10_L18 shows such an anti-correlation only in the lowermost mantle. The density anomalies in expression [117] may again be determined (as in Section 1.23.3.2) by directly inverting the geodynamic observables (Figure 2). Substitution of expression [117] in eqn [104] yields dO,m ¼f,

N X

wi K, ð ; ri Þo ðri Þ

i¼1

    dVS m d ln   ð ri Þ ð ri Þ VS , d lnVS    

dV m d ln  þ ð ri Þ ð ri Þ V , d lnVS

½118

in which the scaling profiles d ln /d ln VS and d ln /d lnV are assumed to vary with depth in

850 Implications for Mantle Dynamics and Composition

the mantle. The radial viscosity profile employed to calculate the kernel functions K, ( ; r) was derived from an inversion based on the mantle heterogeneity in model SB10_L18 (see dash-dotted profile in Figure 9). Occam inversions of expression [118] were then carried out to determine optimally smooth scaling profiles which fit the geodynamic constraints and the results are summarized in Figure 24. The inferred mantle density anomalies (Figure 24(b)) are very similar in amplitude to those determined above on the basis of the seismic shear-velocity models (Figure 12). The major difference is in the mid-mantle region, near 1500 km depth, where we now note that the inferred density anomalies are negatively correlated with the seismic shear-velocity anomalies (a)

0.3

SB 10_L18, δVs SB 10_L18, δVb J 362D 28, δVs J 362D 28, δVb

d ln ρ/d lnV

0.2 0.1 0 –0.1

(Figure 24(c)). Such an anti-correlation cannot be interpreted on the basis of pure thermal effects (see the Karato & Karki scaling profiles in Figure 12 which are positive throughout the lower mantle) and suggest the presence of significant compositional heterogeneity in this region of the lower mantle. A closer inspection of the previous results obtained with the seismic shearvelocity models (Figure 12) also reveals a distinct inflection in the Occam-inferred d ln/d lnVS scaling profiles near 1500 km depth. The fits between the geodynamic data (Figures 2 and 21(a)) and the geodynamic observables predicted on the basis of the Occam-inferred denstity anomalies (Figure 24) are summarized in Table 6. A comparison of these fits with those obtained previously for models J362D28 and SB4_L18 (Tables 1–4) shows that, for model J362D28, there is little improvement obtained by allowing for lateral variations in velocity–density scaling as parametrized by expression [117]. In the case of model SB10_L18, however, there is a significant improvement in fit to all the surface observables, suggesting that the lateral variations in composition implied by expression [117] are supported by the geodynamic data.

–0.2 (b) SB 10_L18 J362D 28

δρ/ρ (%)

0.25

1.23.4.3 Compositional Density Anomalies Inferred from ‘Hot’ and ‘Cold’ Mantle Heterogeneity

0.20 0.15 0.10 0.05 0

Correlation (δρ – δVs)

(c)

1.0 0.5 SB 10_L18 J 363D 28

0 –0.5 –1.0 0

500

1000

1500

2000

2500

3000

Figure 24 Joint shear–bulk-sound inversions for mantle density anomalies. (a) The optimal Occam-inferred velocity– density scaling profiles for the seismic shear (red lines) and bulk-sound (black lines) anomalies in models SB10_L18 (solid lines) and J362D28 (dashed lines). (b) The rms amplitude of the density anomalies, determined on the basis of expression [117] in the main text, for models SB10_L18 (solid lines) and J362D28 (dashed lines). (c) The global cross-correlation between the inferred density anomalies and the seismic shear-velocity anomalies in models SB10_L18 (solid lines) and J362D28 (dashed lines).

Previous interpretations and dynamical models of the large-scale low-velocity anomalies below the central Pacific and below southern Africa (see Figure 3(b)) have identified them as ‘megaplumes’ in the convecting mantle, in which the temperatures are significantly hotter than average mantle and in which the chemical composition may also differ from the overlying mantle (e.g., Becker et al., 1999; Kellogg et al., 1999; Davaille, 1999; Forte and Mitrovica, 2001; Forte et al., 2002; Jellinek and Manga, 2004). These deep-mantle ‘megaplumes’ may be the source regions for the isotopic anomalies which have long been identified in basalts extracted from the Indian and Pacific Oceans (e.g., Hart, 1988; Hofmann, 1997). The hot upwelling plumes in the mantle may therefore have an average chemical composition which differs from that of cold downwellings driven by the descent of dense lithospheric slabs (blue-colored region in Figure 3(b)). We may test the hypothesis of chemically distinct hotter- and colder-than-average regions by parametrizing the mantle density anomalies as follows:

Implications for Mantle Dynamics and Composition Table 6

851

Geodynamic predictions from joint VS/V and ‘hot’/‘cold’ heterogeneity models Gravity

Surf. topography

Horiz. divergence

CMB topography

Model

rms

cor.

fit (%)

rms

cor.

fit (%)

rms

cor.

fit (%)

rms

cor.

fit (%)

J362D28 (dVS/) SB10_L18 (dVS/) TX2002 (dVh/c) SAW24 (dVh/c)

10.4 8.6 11.2 11.6

0.55 0.55 0.70 0.64

28 30 48 39

905 803 753 877

0.75 0.79 0.79 0.75

26 48 54 30

1.3 1.5 1.5 1.4

0.78 0.81 0.84 0.75

51 58 61 51

0.58 0.39 0.53 0.52

0.53 0.37 0.52 0.35

20 10 19 11

rms, root-mean-square amplitude, expressed in units of mGal, m, 108a1, and km for free-air gravity, dynamic surface topography, horizontal divergence, and CMB topography, respectively. cor., global cross-correlation; fit, variance reduction.

d ln d lnVS

      dVS d ln dVS þ d lnVS c VS c h VS h

½119

N X

wi K, ð ; r i Þ o ðr i Þ

i¼1

    dVS d ln ð ri Þ ð ri Þ d lnVS h VS h     m dVS d ln ð ri Þ ð ri Þ þ d lnVS c VS c , 

TX 2002 TX 2002, ‘hot’ TX 2002, ‘cold’ SAW 24B16 SAW 24B16, ‘hot’ SAW 24B16, ‘cold’

2.50

in which (dVS/VS)h (where subscript ‘h’ denotes ‘hot’) represents the regions of the mantle characterized by negative shear-velocity anomalies and presumably hotter-than-average temperatures. Similarly, (dVS/ VS)c (where subscript ‘c’ denotes ‘cold’) represents regions of positive velocity anomalies with presumably colder temperatures. By associating different velocity–density scalings to the ‘hot’ and ‘cold’ regions, expression [119] incorporates a laterally variable d ln/d lnVS profile in the mantle. In the following calculations, we employ the sheartomography models TX2002 (Grand, 2002) and SAW24 (Me´gnin and Romanowicz, 2000). For each model, a spatial separation of the seismic anomalies into faster- (colder-) and slower- (hotter-) than-average regions was carried out and a summary of the rms amplitudes of these anomalies is presented in Figure 25. In both tomography models, we note that the amplitude of the ‘hot’ anomalies is dominant in the lower half of the mantle (below 1500 km depth). Geodynamic constraints on mantle density anomalies are determined, as in Sections 1.23.3.2 and 1.23.4.2, by substituting expression [119] into eqn [104], thereby yielding dO,m ¼ f,

(a) 3.00

½120

The kernel functions K,( ; r) and the associated viscosity profile in this expression are shown in Figure 10. The optimally smooth density–velocity scaling profiles for the ‘hot’, (dVS/VS)h, and ‘cold’, (dVS/VS)c, mantle heterogeneity are obtained through an Occam inversion of all geodynamic

δVS/VS (%)



2.00 1.50 1.00 0.50

(b)

0 1

TX 2002–SAW 24B16

0.8 Correlation

d ¼ 

0.6 0.4 0.2 0

0

500

1000 1500 2000 Depth (km)

2500

3000

Figure 25 ‘Hot’ and ‘cold’ heterogeneity in the mantle. (a) The depth variation of the rms amplitudes of the relative perturbations of seismic shear velocity (dVS/VS) in models TX2002 (Grand 2002) and SAW24 (Me´gnin and Romanowicz 2000) are shown as solid black and dashed black lines, respectively. The rms amplitude of the ‘hot’/‘cold’ heterogeneity in these models (defined as regions where dVS/ VS<0/>0) is shown by the solid and dashed red/blue lines (see legend). (b) The global cross-correlation between the shear-velocity anomalies in models TX2002 and SAW24.

observables in Figure 2 and the results are summarized in Figure 26. For reference and comparison, an inversion was carried out in which no distinction was made between the ‘hot’ and ‘cold’ regions, and the results are also presented in Figure 26. The inverted density–velocity scaling profiles (Figure 26) present a number of interesting features which provide indications of compositional heterogeneity in different depth intervals. In the top 300 km of the

852 Implications for Mantle Dynamics and Composition

(a) 0.6

TX 2002 TX2002, ‘hot’ TX 2002, ‘cold’

d lnρ/d lnVS

0.4 0.2 0

SAW 24B 16 SAW 24B 16, ‘hot’ SAW 24B 16, ‘cold’

–0.2 –0.4

(b) 0.40

TX 2002 TX 2002, ‘hot’ and ‘cold’ SAW 24B 16 SAW 24B 16, ‘hot’ and ‘cold’

δρ/ρ (%)

0.30 0.20 0.10 0

Correlation (δρ−δVS)

(c) 1.0 0.5 0 –0.5 –1.0 0

500

1000

1500

2000

2500

3000

Figure 26 Inversions for ‘hot’ and ‘cold’ mantle density anomalies. (a) The Occam-inferred velocity–density scaling profiles for the ‘hot’/‘cold’ regions of the mantle (as defined in the text and in Figure 25) in models TX2002 (Grand, 2002) and SAW24 (Me´gnin and Romanowicz, 2000) are shown as red/blue solid and dashed lines, respectively (see legend). The velocity–density scaling profiles which are inferred if no distinction is made between the ‘hot’/‘cold’ regions are shown as black solid and dashed lines for models TX2002 and SAW24, respectively. (b) The rms amplitude of the mantle density anomalies, calculated on the basis of expression [119] and the scaling profiles in the top frame, are shown for both models TX2002 and SAW24 as the solid and dashed orange lines, respectively. The amplitude of the density anomalies which is inferred if no distinction is made between ‘hot’/‘cold’ regions are shown as black (solid and dashed) lines. (c) The global cross-correlation between the inferred density anomalies and the seismic shear-velocity anomalies in models TX2002 and SAW24 are shown by the solid and dashed orange lines, respectively. When no distinction is made between ‘hot’/‘cold’ regions, the correlations are shown by the black solid and dashed lines (see legend in (b)).

mantle, the scaling coefficient for the ‘cold’ anomalies is significantly smaller than for the ‘hot’ anomalies and in the top 100 km it is nearly zero. The ‘cold’ anomalies in this depth range are essentially associated with the subcontinental roots and therefore we have another confirmation of previous geodynamic studies (e.g.,

Forte and Perry, 2000) which revealed a mutual cancellation of the thermal and compositional contributions to mantle density below continents. In the lower mantle, the scaling coefficient for the ‘cold’ anomalies is significantly greater than for the ‘hot’ anomalies, where the scaling for the latter actually goes to negative values throughout the mid lower mantle. This inference suggests large amounts of compositional heterogeneity in the lower-mantle ‘megaplumes’ below the Pacific Ocean and below southern Africa. The lower-mantle trend is reversed in the deepest portions of the lower mantle (e.g., in the D0 layer), where we again find that the scaling for the ‘hot’ regions is larger than for the ‘cold’ regions. A consideration of the cross-correlation between the inferred mantle density anomalies and the seismic shear-velocity anomalies (Figure 26) also shows anomalous mid-mantle and deep lower-mantle regions, where the correlations decrease significantly. These decreased correlations explain the inflections in the previously inferred d ln /d ln VS scaling profiles (Figure 12) in the mid- and deep lower mantle. These decorrelations between mantle density and seismic shear-velocity anomalies suggests that thermal contributions are not the sole source of density heterogeneity in these regions and there is a likelihood for significant composition heterogeneity in these portions of the lower mantle. Indeed, a consideration of the fits to the geodynamic data summarized in Table 6 shows a distinct improvement compared to that obtained previously (Tables 1–4). 1.23.4.4 Diffuse Mid-Mantle Compositional Horizon As discussed in the preceding two sections, there appears to be a significant deviation from purely thermal control of mantle heterogeneity in a region of the mid-mantle near 1500 km depth. This was already evident in Figure 12, where the geodynamically inferred values of d ln /d ln VS in the midmantle differ substantially from those inferred by Karato and Karki (2001) on the basis of mineral physics data and assuming that temperature variations are the dominant contributor. The subsequent geodynamic inversions which allowed lateral variations in d ln /d ln VS (Figures 24 and 26) further underlined the anomalous thermochemical behavior of the mid-mantle region. Hints of the anomalous properties of the midmantle region from geodynamic models can also be gleaned from previously published comparisons between subducted slab heterogeneity and seismic

Implications for Mantle Dynamics and Composition

tomography models. For example, LithgowBertelloni and Richards (1998, see their figure 7) found that the lowest correlations between slab heterogeneity and seismically inferred heterogeneity in the lower mantle occurred near 1500 km and that the longest wavelength components of these two fields were anti-correlated at this depth. These authors suggested that these poor correlations could be attributed to insufficient resolution of the global tomography models at these depths. An alternative explanation, if the slabs are viewed as a proxy for thermal heterogeneity, is that the low slab-tomography correlations in the mid-mantle region are a diagnostic for significant compositional heterogeneity at these depths. There have been a number of seismological, geodynamic, and geochemical arguments for the presence of a mid-mantle boundary or interface which segregates a presumably stagnant, intrinsically heavy layer in the bottom half of the lower mantle from the convective flow in the overlying mantle (e.g., Hofmann, 1997; Kellogg et al., 1999; Davaille, 1999). More recently, theoretical and experimental advances in mineral physics have also identified an anomalous region in the mid-mantle (e.g., Badro et al., 2003) and evidence for distinct differences between the elastic and thermodynamic properties in the deep lower mantle, below 1400 km depth, compared to the top of the lower mantle (e.g., Wentzcovich et al., 2004). The geodynamic inversions presented above show that ‘normal’, thermally dominated buoyancy seems to prevail in regions sufficiently far above and below the anomalous mid-mantle region near 1500 km depth. Recent joint seismic–geodynamic inversions (Simmons et al., 2006) have directly tested the hypothesis of a sharp chemical interface at 1800 km depth which inhibits vertical flow across the lower mantle. The fits to the geodynamic data were found to be significantly degraded and this appears to rule out a sharply defined compositional boundary in the lower mantle. These seismic–geodynamic inversions do not, however, rule out the possibility of a diffuse, compositional (or phase-change) horizon at this depth which may have significant vertical undulations, perhaps over distances of a few hundred kilometers, in order to explain the depth interval over which the anomalous values of d ln /d ln VS are inferred (Figures 12, 24, and 26). The verification of this hypothesis will require further investigation, with the greatest insights provided by an integrated approach based on combined seismic,

853

geodynamic, and mineral physical modelling. Recent progress in this direction has been made by Simmons et al. (2007).

1.23.5 Concluding Remarks The current global-scale seismic tomography models have been shown to provide good fits to a wide variety of surface geodynamic observables related to mantle convection. These results suggest that global tomography does indeed provide reliable images of the lateral variations in mantle structure which are produced by the thermal convection process in Earth’s mantle. The tabulated summaries of the results (Tables 1–6) show however that the fits can still be improved and this raises the fundamental question as to what aspects of the global tomography models, or indeed the mantle flow calculations, may be deficient or inadequate. From the perspective of the mantle flow theory, there is a major issue concerning the possible dynamical importance of lateral variations in mantle viscosity which have not been explicitly modeled in the calculations presented above. It is important to recall that the mantle flow calculations include a coupling to rigid tectonic plates at the surface and that the motions of these plates are predicted on the basis of the buoyancy induced flow in the mantle. The surface plates are likely the most important manifestation of lateral rheology variations in the Earth and hence their inclusion in mantle flow modeling is essential. However, the flow model employed here assumes a 1-D radial variation of vis-cosity in the mantle and the question then arises as to the importance of neglecting lateral viscosity variations (LVVs) in the deep mantle. Detailed modeling of the potential bias which may be due to large-scale LVV, estimated on the basis of global tomography models, has recently been carried out by Moucha et al. (2007). One of the main conclusions of this recent study is that the impact of large scale LVV on convectionrelated observables is no larger than the uncertainties arising from differences among the recent tomography models. This relative insensitivity to LVV may well explain why the current tomography-based viscous flow models, which assumes a purely depthdependent mantle viscosity, can yield such a good match to a wide range of convection data. Simultaneous inversions of global seismic and geodynamic data (e.g., Forte et al., 1994; Forte and Woodward, 1997; Simmons et al., 2006) have shown

854 Implications for Mantle Dynamics and Composition

that it is possible to obtain greatly improved fits to the global convection-related data using the mantle flow theory which was presented above. The improved fits are found to be the result of relatively small adjustments (mainly in amplitude) to the seismically inferred heterogeneity in the mantle, especially the mid-mantle region, which are also compatible with the constraints imposed by the geodynamic data. These joint seismic–geodynamic inversions therefore support the basic conclusion of this study, namely that seismic tomography and geodynamics can be reconciled with the large-scale heterogeneity in the mantle that is clearly generated by the mantle convection process.

Acknowledgments This work is supported by grants from the Canadian National Science and Engineering Research Council and the Canada Research Chair Program. Support provided by the Canadian Institute for Advanced Research – Earth Systems Evolution Program is also gratefully acknowledged.

References Antolik M, Gu Yu J, Ekstro¨m G, and Dziewonski AM (2003) J362D28; a new joint model of compressional and shear velocity in the Earth’s mantle. Geophysical Journal International 153: 443–466. Badro J, Fiquet G, Guyot F, et al. (2003) Iron partitioning in Earth’s mantle: Toward a deep lower mantle discontinuity. Science 300: 789–791. Bassin C, Laske G, and Masters G (2000) The current limits of resolution for surface wave tomography in North America. Eos Transactions of the American Geophysical Union 81(48): 897. Becker TW and Boschi L (2002) A comparison of tomographic and geodynamic mantle models. Geochemistry Geophysics Geosystems 3, doi:10.129/2001GC000168. Becker TW, Kellogg JB, and O’Connell RJ (1999) Thermal constraints on the survival of primitive blobs in the lower mantle. Earth and Planetary Science Letters 171: 351–365. Boschi L and Dziewonski AM (2000) Whole Earth tomography from delay times of P, PcP, and PKP phases: Lateral heterogeneities in the outer core or radial anisotropy in the mantle? Journal of Geophysical Research 105: 13675–13696. Bunge H-P, Richards MA, and Baumgardner JR (1996) Effect of depth-dependent viscosity on the planform of mantle convection. Nature 379: 436–438. Bunge H-P, Richards MA, Lithgow-Bertelloni C, Baumgardner JR, Grand SP, and Romanowicz BA (1998) Time scales and heterogeneous structure in geodynamic Earth models. Science 280: 91–95. Cˇadek O and Fleitout L (2003) Effect of lateral viscosity variations in the top 300 km on the geoid and dynamic

topography. Geophysical Journal International 152: 566–580. Cˇadek O, Yuen DA, Steinbach V, Chopelas A, and Matyska C (1994) Lower mantle thermal structure deduced from seismic tomography, mineral physics and numerical modelling. Earth and Planetary Science Letters 121: 385–402. Cammarano F, Goes S, Vacher P, and Giardini D (2003) Inferring upper-mantle temperatures from seismic velocities. Physics of the Earth and Planetary Interiors 138: 197–222. Cathles LM (1975) The Viscosity of the Earth’s Mantle, Princeton: Princeton University Press. Chase CG (1979) Subduction, the geoid, and lower mantle convection. Nature 282: 464–468. Chase CG and Sprowl DR (1983) The modern geoid and ancient plate boundaries. Earth and Planetary Science Letters 62: 314–320. Clayton RW and Comer RP (1983) A tomographic analysis of mantle heterogeneities from body wave travel times. Eos Transactions of the American Geophysical Union 64(45): 776. Constable SC, Parker RL, and Constable CG (1987) Occam’s inversion: A practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics 52: 289–300. Corrieu V, Ricard Y, and Froidevaux C (1994) Converting mantle tomography into mass anomalies to predict the Earth’s radial viscosity. Physics of the Earth and Planetary Interiors 84: 3–13. Corrieu V, Thoraval C, and Ricard Y (1995) Mantle dynamics and geoid Green functions. Geophysical Journal International 120: 516–523. Crosby AG, McKenzie D, and Sclater JG (2006) The relationship between depth, age and gravity in the oceans. Geophysical Journal International 166: 553–573. Davaille A (1999) Simultaneous generation of hotspots and superswells by convection in a heterogeneous planetary mantle. Nature bf402: 756–760. De Mets CR, Gordon RG, Argus DF, and Stein S (1990) Current plate motions. Physics of the Earth and Planetary Interiors 101: 425–478. Defraigne P (1997) Internal loading of a compressible Earth: Effect of a solid lithosphere. Physics of the Earth and Planetary Interiors 101: 303–313. Deschamps F, Snieder R, and Trampert J (2001) The relative density-to-shear velocity scaling in the uppermost mantle. Physics of the Earth and Planetary Interiors 124: 193–212. Doin M-P and Fleitout L (2000) Flattening of the oceanic topography and geoid: Thermal versus dynamic origin. Geophysical Journal International 143: 582–594. Dziewonski AM (1984) Mapping the lower mantle: Determination of lateral heterogeneity in P velocity up to degree and order 6. Journal of Geophysical Research 89: 5929–5952. Dziewonski AM and Anderson DL (1981) Preliminary reference Earth model. Physics of the Earth and Planetary Interiors 25: 297–356. Dziewonski AM, Hager BH, and O’Connell RJ (1977) Large scale heterogeneity in the lower mantle. Journal of Geophysical Research 82: 239–255. Ekstro¨m G and Dziewonski AM (1998) The unique anisotropy of the Pacific upper mantle. Nature 394: 168–172. Forte AM (2000) Seismic-geodynamic constraints on mantle flow: Implications for layered convection, mantle viscosity, and seismic anisotropy in the deep mantle. In: Karato S-i, Forte AM, Liebermann RC, Masters G, and Stixrude L (eds.) Geophysical Monograph Series 117 Earth’s Deep Interior: Mineral Physics and Tomography from the Atomic to the Global Scale, pp. 3–36. Washington DC: American Geophysical Union. Forte AM and Peltier WR (1987) Plate tectonics and aspherical Earth structure: The importance of poloidal–toroidal coupling. Journal of Geophysical Research 92: 3645–3679.

Implications for Mantle Dynamics and Composition Forte AM and Peltier WR (1989) Core-mantle boundary topography and whole-mantle convection. Geophysical Research Letters 16: 621–624. Forte AM and Peltier WR (1991) Viscous flow models of global geophysical observables. Part 1: Forward problems. Journal of Geophysical Research 96: 20131–20159. Forte AM and Peltier WR (1994) The kinematics and dynamics of poloidal-toroidal coupling in mantle flow: The importance of surface plates and lateral viscosity variations. Advances in Geophysics 36: 1–119. Forte AM and Mitrovica JX (1996) New inferences of mantle viscosity from joint inversion of long-wavelength mantle convection and post–glacial rebound data. Geophysical Research Letters 23: 1147–1150. Forte AM and Woodward RL (1997a) Seismic-geodynamic constraints on vertical flow between the upper and lower mantle: The dynamics of the 670 km seismic discontinuity. In: Crossley D (ed.) The Fluid Dynamics of Astrophysics and Geophysics, vol. 7: Earth’s Deep Interior, pp. 337–404. Newark, N.J: Gordon and Breach. Forte AM and Woodward RL (1997b) Seismic-geodynamic constraints on three-dimensional structure, vertical flow, and heat transfer in the mantle. Journal of Geophysical Research 102: 17981–17994. Forte AM and Perry HKC (2000) Geodynamic evidence for a chemically depleted continental tectosphere. Science 290: 1940–1944. Forte AM and Mitrovica JX (2001) Deep-mantle high-viscosity flow and thermochemical structure inferred from seismic and geodynamic data. Nature 410: 1049–1056. Forte AM, Peltier WR, and Dziewonski AM (1991) Inferences of mantle viscosity from tectonic plate velocities. Geophysical Research Letters 18: 1747–1750. Forte AM, Dziewonski AM, and Woodward RL (1993a) Aspherical structure of the mantle, tectonic plate motions, nonhydrostatic geoid, and topography of the core-mantle boundary. In: Le Moue¨l J-L, Smylie DE, and Herring T (eds.) Dynamics of the Earth’s Deep Interior and Earth Rotation, pp. 135–166. Washington DC: American Geophysical Union. Forte AM, Peltier WR, Dziewonski AM, and Woodward RL (1993b) Dynamic surface topography: A new interpretation based upon mantle flow models derived from seismic tomography. Geophysical Research Letters 20: 225–228. Forte AM, Peltier WR, Dziewonski AM, and Woodward RL (1993c) Reply to comment by M. Gurnis. Geophysical Research Letters 20: 1665–1666. Forte AM, Woodward RL, and Dziewonski AM (1994) Joint inversions of seismic and geodynamic data for models of three-dimensional mantle heterogeneity. Journal of Geophysical Research 99: 21857–21877. Forte AM, Mitrovica JX, and Woodward RL (1995a) Seismicgeodynamic determination of the origin of excess ellipticity of the core-mantle boundary. Geophysical Research Letters 22: 1013–1016. Forte AM, Dziewonski AM, and O’Connell RJ (1995b) Continentocean chemical heterogeneity in the mantle based on seismic tomography. Science 268: 386–388. Forte AM, Mitrovica JX, and Espesset A (2002) Geodynamic and seismic constraints on the thermochemical structure and dynamics of convection in the deep mantle. Philosophical Transactions of the Royal Society of London A 360: 2521–2543. Gable CW, O’Connell RJ, and Travis BJ (1991) Convection in three dimensions with surface plates. Journal of Geophysical Research 96: 8391–8405. Garcia R and Souriau A (2000) Amplitude of the core-mantle boundary topography estimated by stochastic analysis of

855

core phases. Physics of the Earth and Planetary Interiors 117: 345–359. Gordon RG and Jurdy DM (1986) Cenozoic global plate motions. Journal of Geophysical Research 91: 12389–12406. Grand SP (2002) Mantle shear-wave tomography and the fate of subducted slabs. Philosophical Transactions of the Royal Society of London A 360: 2475–2491. Gurnis M (1990) Bounds on global dynamic topography from Phanerozoic flooding of continental platforms. Nature 344: 754–756. Gurnis M (1993) Comment. Geophysical Research Letters 20: 1663–1664. Gwinn CR, Herring TA, and Shapiro II (1986) Geodesy by radio interferometry; studies of the forced nutations of the Earth. Part 2: Interpretation. Journal of Geophysical Research 91: 4755–4766. Hager BH (1984) Subducted slabs and the geoid: Constraints on mantle rheology and flow. Journal of Geophysical Research 89: 6003–6015. Hager BH and O’Connell RJ (1981) A simple global model of plate dynamics and mantle convection. Journal of Geophysical Research 86: 4843–4867. Hager BH and Richards MA (1989) Long-wavelength variations in the Earth’s geoid: physical models and dynamical implications. Philosophical Transactions of the Royal Society of London A 328: 309–327. Hager BH and Clayton RW (1989) Constraints on the structure of mantle convection using seismic observations, flow models, and the geoid. In: Peltier WR (ed.) Mantle Convection: Plate Tectonics and Global Dynamics, pp. 657–763. New York: Gordon and Breach. Hager BH, Clayton RW, Richards MA, Comer RP, and Dziewonski AM (1985) Lower mantle heterogeneity, dynamic topography and the geoid. Nature 313: 541–545. Hart SR (1988) Heterogeneous mantle domains; signatures, genesis and mixing chronologies. Earth and Planetary Science Letters 90: 273–296. Haskell NA (1935) The motion of a viscous fluid under a surface load. Physics (NY) 6: 265–269. Herring TA, Gwinn CR, and Shapiro II (1986) Geodesy by radio interferometry; studies of the forced nutations of the Earth. Part 1: Data analysis. Journal of Geophysical Research 91: 4745–4754. Hofmann AW (1997) Mantle geochemistry: The message from oceanic volcanism. Nature 385: 219–229. Jackson JD (1975) Classical Electrodynamics. New York: John Wyley and Sons. Jackson I (1998) Elasticity, composition and temperature of the Earth’s lower mantle: A reappraisal. Geophysical Journal International 134: 291–311. Jackson I (2001) Chemical composition and temperature of the lower mantle from seismological models: Residual uncertainties. Integrated models of Earth structure and evolution, AGU Virtual Meeting - 2001 Spring Meeting, 20 June, [online]. Jackson JM, Zhang J, Shu J, Sinogeikin SV, and Bass JD (2005) High-pressure sound velocities and elasticity of aluminous MgSiO3 perovskite to 45 GPa: Implications for lateral heterogeneity in Earth’s lower mantle. Geophysical Research Letters 32: L21305 (doi:10.1029/2005GL023522). Jarvis GT and Peltier WR (1982) Mantle convection as a boundary layer phenomenon. Geophysical Journal of the Royal Astronomical Society 68: 385–424. Jeanloz R and Thompson AB (1983) Phase transitions and mantle discontinuities. Reviews of Geophysics and Space Physics 21: 51–74. Jellinek AM and Manga M (2004) Links between long-lived hot spots, mantle plumes, D0, and plate tectonics 42: RG3002 (doi:10.1029/2003RG000144).

856 Implications for Mantle Dynamics and Composition Jordan TH (1981) Continents as a chemical boundary layer. Philosophical Transactions of the Royal Society of London A 301: 359–373. Karato S-i (1993) Importance of anelasticity in the interpretation of seismic tomography. Geophysical Research Letters 20: 1623–1626. Karato S-I and Karki BB (2001) Origin of lateral variation of seismic wave velocities and density in the deep mantle. Journal of Geophysical Research 106: 21771–21784. Karpychev M and Fleitout L (2000) Long-wavelength geoid: The effect of continental roots and lithosphere thickness variations. Geophysical Journal International 143: 945–963. Kaufmann G and Lambeck K (2000) Mantle dynamics, postglacial rebound and the radial viscosity profile. Physics of the Earth and Planetary Interiors 121: 301–324. Kellogg LH, Hager BH, and van der Hilst RD (1999) Compositional stratification in the deep mantle. Science 283: 1881–1884. Kennett BLN, Widiyantoro S, and van der Hilst RD (1998) Joint seismic tomography for bulk-sound and shear wavespeed in the Earth’s mantle. Journal of Geophysical Research 103: 12469–12493. Kido M, Yuen DA, Cˇadek O, and Nakakuki T (1998) Mantle viscosity derived by genetic algorithm using oceanic geoid and seismic tomography for whole-mantle versus blockedflow situations. Physics of the Earth and Planetary Interiors 107: 307–326. King S and Masters G (1992) An inversion for the radial viscosity structure using seismic tomography. Geophysical Research Letters 19: 1551–1554. King SD (1995) Models of mantle viscosity. In: Ahrens TJ (ed.) AGU Reference Shelf Series, vol. 2: A Handbook of Physical Constants: Mineral Physics and Crystallography, pp. 227–236. Washington DC: American Geophysical Union. Lambeck K, Johnston P, Smither C, and Nakada M (1996) Glacial rebound of the British Isles. III. Constraints on mantle viscosity. Geophysical Journal International 125: 340–354. Lambeck K, Smither C, and Johnston P (1998) Sea-level change, glacial rebound and mantle viscosity for northern Europe. Geophysical Journal International 134: 102–144. Landau LD and Lifshitz EM (1959) Fluid Mechanics, vol. 6. New York: Pergamon. Lemoine F, Pavlis N, Kenyon S, Rapp R, Pavlis E, and Chao B (1998) New high-resolution model developed for Earth’s gravitational field. Eos Transactions of the American Geophysical Union 79(9): 113. Le Stunff Y and Ricard Y (1995) Topgraphy and geoid due to lithospheric mass anomalies. Geophysical Journal International 122: 982–990. Le Stunff Y and Ricard Y (1997) Partial advection of equidensity surfaces; a solution for the dynamic topography problem? Journal of Geophysical Research 102: 24655–24667. Lithgow-Bertelloni C and Gurnis M (1997) Cenozoic subsidence and uplift of continents from time-varying dynamic topography. Geology 25: 735–738. Lithgow-Bertelloni C and Richards MA (1998) The dynamics of Cenozoic and Mesozoic plate motions. Reviews of Geophysics 36: 27–78. Martinec Z, Matyska C, Cadek O, and Hrdina P (1993) The Stokes problem with 3-D Newtonian rheology in a sphericall shell. Computer Physics Communications 76: 63–79. Masters G, Laske G, Bolton H, and Dziewonski AM (2000) The relative behaviour of shear velocity, bulk sound speed, and compressional velocity in the mantle: Implications for chemical and thermal structure. In: Karato S-i, Forte AM, Libebermann RC, Masters G, and Stixrude (eds.) Geophysical Monograph Series, 117: Earth’s Deep Interior: Mineral Physics and Tomography From the Atomic to the

Global Scale, pp. 63–87. Washington DC: American Geophysical Union, 2000. Mathews PM, Herring TA, and Buffett BA (2002) Modeling of nutation and precession: New nutation series for nonrigid Earth and insights into the Earth’s interior. Journal of Geophysical Research 107: 2068 (doi:10.1029/2001JB000390). McConnell RK (1968) Viscosity of the mantle from relaxation time spectra of isostatic adjustment. Journal of Geophysical Research 73: 7089–7105. McNamara AK and Zhong S-J (2005) Thermochemical structures beneath Africa and the Pacific Ocean. Nature 437: 1136–1139. Me´gnin C and Romanowicz B (2000) The three-dimensional shear velocity structure of the mantle from the inversion of body, surface and higher-mode waveforms. Geophysical Journal International 143: 709–728. Mitrovica JX (1996) Haskell (1935) revisited. Journal of Geophysical Research 101: 555–569. Mitrovica JX and Forte AM (1997) Radial profile of mantle viscosity: Results from the joint inversion of convection and post-glacial rebound observables. Journal of Geophysical Research 102: 2751–2769. Mitrovica JX and Forte AM (2004) A new inference of mantle viscosity based upon joint inversion of convection and glacial isostatic adjustment data. Earth and Planetary Science Letters 225: 177–189. Morelli A and Dziewonski AM (1987) Topography of the core– mantle boundary and lateral homogeneity of the liquid core. Nature 325: 678–683. Moucha R, Forte AM, Mitrovica JX, and Daradich A (2007) Geodynamic implications of lateral variations in mantle rheology on convection related observables and inferred viscosity models. Geophysical Journal International 169: 113–135. Nakada M and Lambeck K (1989) Late Pleistocene and Holocene sea-level change in the Australian region and mantle rheology. Geophysical Journal International 96: 497–517. Nicolas A and Poirier J-P (1976) Crystalline Plasticity and Solid State Flow in Metamorphic Rocks, p444. London: John Wiley and Sons. O’Connell RJ (1971) Pleistocene glaciation and viscosity of the lower mantle. Geophysical Journal of the Royal Astronomical Society 23: 299–327. Obayashi M and Fukao YP (1997) P and PcP travel time tomography for the core–mantle boundary. Journal of Geophysical Research 102: 17825–17841. Oganov AR, Brodholt JP, and Price GD (2001a) Ab initio elasticity and thermal equation of state of MgSiO3 perovskite. Earth and Planetary Science Letters 184: 555–560. Oganov AR, Brodholt JP, and Price GD (2001b) The elastic constants of MgSiO3 perovskite at pressures and temperatures of the Earth’s mantle. Nature 411: 934–937. Mattern E, Matas J, Ricard Y, and Bass J (2005) Lower mantle composition and temperature from mineral physics and thermodynamic modelling. Geophysical Journal International 160: 973–990. Muller RD, Roest WR, Royer J-Y, Gahagan LM, and Sclater JG (1997) Digital isochrons of the world’s ocean floor. Journal of Geophysical Research 102: 3211–3214. Panasyuk SV and Hager BH (2000) Inversion for mantle viscosity profiles constrained by dynamic topography and the geoid, and their estimated errors. Geophysical Journal International 143: 821–836. Panasyuk SV, Hager BH, and Forte AM (1996) Understanding the effects of mantle compressibility on geoid kernels. Geophysical Journal International 124: 121–133.

Implications for Mantle Dynamics and Composition Pari G (2001) Crust 5.1-based inference of the Earth’s dynamic surface topography: Geo-dynamic implications. Geophysical Journal International 144: 501–516. Pari G and Peltier WR (2000) Subcontinental mantle dynamics; a further analysis based on the joint constraints of dynamic surface topography and free-air gravity. Journal of Geophysical Research 105: 5635–5662. Parsons B and McKenzie D Mantle convection and the thermal structure of the plates. Journal of Geophysical Research 83: 4485–4496. Parsons B and Sclater JG (1977) An analysis of the variation of ocean floor bathymetry and heat flow with age. Journal of Geophysical Research 82: 803–827. Parsons B and Daly S (1983) The relationship between surface topography, gravity anomalies, and temperature structure of convection. Journal of Geophysical Research 88: 1129–1144. Peltier WR (1985) New constraint on transient lower mantle rheology and internal mantle buoyancy from glacial rebound data. Nature 318: 614–617. Peltier WR and Andrews JT (1976) Glacial isostatic adjustment, I, The forward problem. Geophysical Journal of the Royal Astronomical Society 46: 605–646. Pekeris CL (1935) Thermal convection in the interior of the Earth. Supplement 3: 343–367. Perry HKC, Forte AM, and Eaton DWS (2003) Upper-mantle thermochemical structure below North America from seismicgeodynamic flow models. Geophysical Journal International 154: 279–299. Phinney RA and Burridge R (1973) Representation of the elasticgravitational excitation of a spherical Earth model by generalized spherical harmonics. Geophysical Journal of the Royal Astronomical Society 34: 451–487. Que´re´ S and Forte AM (2006) Influence of past and present-day plate motions on spherical models of mantle convection: Implications for mantle plumes and hotspots. Geophysical Journal International 165: 1041–1057. Ranalli G (2001) Mantle rheology: Radial and lateral viscosity variations inferred from microphysical creep laws. Journal of Geodynamics 32: 425–444. Ricard Y and Vigny C (1989) Mantle dynamics with induced plate tectonics. Geophysical Journal International 94: 17543–17559. Ricard Y and Wuming B (1991) Inferring viscosity and the 3-D density structure of the mantle from geoid, topography and plate velocities. Geophysical Journal International 105: 561–572. Ricard Y, Fleitout L, and Froidevaux C (1984) Geoid heights and lithospheric stresses for a dynamic Earth. Annales Geophysicae 2: 267–286. Ricard Y, Froidevaux C, and Fleitout L (1988) Global plate motion and the geoid: A physical model. Geophysical Journal International 93: 477–484. Ricard Y, Vigny C, and Froidevaux C (1989) Mantle heterogeneities, geoid, and plane motion – A Monte Carlo inversion. Journal of Geophysical Research 94: 13739–13754. Ricard Y, Richards M, Lithgow-Bertelloni C, and Le Stunff Y (1993) A geodynamic model of mantle density heterogeneity. Journal of Geophysical Research 98: 21895–21909. Ricard Y, Chambat F, and Lithgow-Bertelloni C (2006) Gravity observations and 3D structure of the Earth. Comptes Rendus Geosciences 338: 992–1001. Richards MA (1991) Hot spots and the case for a high viscosity lower mantle. In: Sabadini R, Lambeck K, and Boschi E (eds.) NATO ASI Series C, vol. 334: Glacial Isostasy, Sea-Level and Mantle Rheology, pp. 571–588. Dordrecht: Kluwer Academic.

857

Richards MA and Hager BH (1984) Geoid anomalies in a dynamic Earth. Journal of Geophysical Research 89: 5987–6002. Richards MA and Hager BH (1989) Effects of lateral viscosity variations on long-wavelength geoid anomalies and topography. Journal of Geophysical Research 94: 10299–10313. Richards MA and Engebretson DC (1992) Large-scale mantle convection and the history of subduction. Nature 355: 437–440. Richards MA, Ricard Y, Lithgow-Bertelloni C, Spada G, and Sabadini R (1997) An explanation for the long-term stability of Earth’s rotation axis. Science 275: 372–375. Ritsema J, van Heijst HJ, and Woodhouse JH (1999) Complex shear velocity structure imaged beneath Africa and Iceland. Science 286: 1925–1928. Ritzwoller MH, Shapiro NM, and Zhong S-J (2004) Cooling history of the Pacific lithosphere. Earth and Planetary Science Letters 226: 69–84. Rodgers A and Wahr J (1993) Inference of core-mantle boundary topography from ISC PcP and PKP traveltimes. Geophysical Journal International 115: 991–1011. Ro¨hm AHE, Snieder R, Goes S, and Trampert J (2000) Thermal structure of continental upper mantle inferred from S-wave velocity and surface heat flow. Earth and Planetary Science Letters 181: 395–407. Romanowicz B (2003) Global mantle tomography: Progress status in the past 10 years. Annual Review of Earth and Planetary Sciences 31: 303–328. Sabadini R, Yuen DA, and Gasperini P (1985) The effects of transient rheology on the interpretation of lower mantle rheology. Geophysical Research Letter 12: 361–365. Sabadini R and Yuen DA (1989) Mantle stratification and longterm polar wander. Nature 339: 373–375. Sammis CG, Smith JC, Schubert G, and Yuen DA (1977) Viscosity-depth profile of the Earth’s mantle: Effects of polymorphic phase transitions. Journal of Geophysical Research 82: 3747–3761. Simmons NA, Forte AM, and Grand SP (2006) Constraining modes of mantle flow with seismic and geodynamic data: A joint approach. Earth and Planetary Science Letters 246: 109–124. Simmons NA, Forte AM, and Grand SP (2007) Thermochemical structure and dynamics of the African superplume. Geophysical Research Letter 34: L02301 (doi:10.1029/ 2006GL028009). Sobolev SV, Zeyen H, Granet M, et al. (1997) Upper mantle temperatures and lithosphere-asthenosphere system beneath the French Massif Central constrained by seismic, gravity, petrologic and thermal observations. Tectonophysics 275: 143–164. Spada G, Ricard Y, and Sabadini R (1992) Excitation of true polar wander by subduction. Nature 360: 452–454. Stacey FD (1998) Thermoelasticity of a mineral composite and a reconsideration of lower mantle properties. Physics of the Earth and Planetary Interiors 106: 219–236. Stein CA and Stein S (1992) A model for the global variation in oceanic depth and heat flow with lithospheric age. Nature 359: 123–129. Steinberger B and O’Connell RJ (1997) Change of the Earth’s rotation axis inferred from the advection of mantle density heterogeneities. Nature 387: 169–173. Stixrude L and Lithgow-Bertelloni C (2001) The origin of lateral heterogeneity in the mantle. Integrated models of Earth structure and evolution, AGU Virtual Meeting - 2001 Spring Meeting, 20 June, [online]. Stixrude L and Lithgow-Bertelloni C (2005a) Mineralogy and elasticity of the oceanic upper mantle: Origin of the low-

858 Implications for Mantle Dynamics and Composition velocity zone. Journal of Geophysical Research 110: B03204 (doi:10.1029/2004JB002965). Stixrude L and Lithgow-Bertelloni C (2005b) Thermodynamics of mantle minerals. Part I: Physical properties. Geophysical Journal International 162: 610–632. Su W-J, Woodward RL, and Dziewonski AM (1992) Deep origin of mid-ocean ridge seismic velocity anomalies. Nature 360: 149–152. Su W-J and Dziewonski AM (1997) Simultaneous inversion for 3-D variations in shear and bulk velocity in the mantle. Physics of the Earth and Planetary Interiors 100: 135–156. Sze EKM and van der Hilst RD (2003) Core mantle boundary topography from short period PcP, PKP, and PKKP data. Physics of the Earth and Planetary Interiors 135: 27–46. Taylor SR and McLennan SM (1995) The geochemical evolution of the continental crust. Reviews of Geophysics 33: 241–265. Thoraval C, Machetel P, and Cazenave A (1995) Locally layered convection inferred from dynamic models of the Earth’s mantle. Nature 375: 777–780. Thoraval C and Richards MA (1997) The geoid constraint in global geodynamics: Viscosity structure, mantle heterogeneity models, and boundary conditions. Geophysical Journal International 131: 1–8. Trampert J and Woodhouse JH (1996) High resolution global phase velocity distributions. Geophysical Research Letters 23: 21–24. Trampert J, Vacher P, and Vlaar N (2001) Sensitivities of seismic velocities to temperature, pressure and composition in the lower mantle. Physics of the Earth and Planetary Interiors 124: 255–267. Tushingham AM and Peltier WR (1991) ICE-3G: A new global model of late Pleistocene deglaciation based upon geophysical predictions of postglacial relative sea level change. Journal of Geophysical Research 96: 4497–4523.

Wahr J (1981) The forced mutations of an elliptical, rotating elastic and oceanless Earth. Geophysical Journal of the Royal Astronomical Society 64: 705–727. Wang Y and Weidner DJ (1996) (q/qT)P of the lower mantle. Pure and Applied Geophysics 146: 533–549. Weertman J (1978) Creep laws for the mantle of the Earth. Philosophical Transactions of the Royal Society of London A 288: 9–26. Weertman J and Weertman JR (1975) High temperature creep of rock and mantle viscosity. Annual Review of Earth and Planetary Sciences 3: 293–315. Wentzcovitch RM, Karki BB, Cococcioni M, and de Gironcoll S (2004) Thermoelastic properties of MgSiO3-perovskite: Insights on the nature of the Earth’s lower mantle. Physical Review Letters 92: 018501. Woodhouse JH and Dziewonski AM (1984) Mapping the upper mantle: Three-dimensional modeling of earth structure by inversion of seismic waveforms. Journal of Geophysical Research 89: 5953–5986. Woodward RL and Masters G (1991) Global upper mantle structure from long-period differential travel times. Journal of Geophysical Research 96: 6351–6377. Wu P and Peltier WR (1983) Glacial isostatic adjustment and the free air gravity anomaly as a constraint on deep mantle viscosity. Geophysical Journal of the Royal Astronomical Society 74: 377–449. Zhang S and Christensen U (1993) Some effects of lateral viscosity variations on geoid and surface velocities induced by density anomalies in the mantle. Geophysical Journal International 114: 531–547. Zhang J and Weidner DJ (1999) Thermal equation of state of aluminum-enriched silicate perovskite. Science 284: 782–784. Zhang S and Yuen DA (1995) The influences of lower mantle viscosity strati_cation on 3D spherical-shell mantle convection. Earth and Planetary Science Letters 132: 157–166.