Density structure of the mantle transition zone and the dynamic geoid

Journal of Geodynamics 59–60 (2012) 183–192

Contents lists available at SciVerse ScienceDirect

Journal of Geodynamics journal homepage: http://www.elsevier.com/locate/jog

Density structure of the mantle transition zone and the dynamic geoid Mikhail K. Kaban a,b,∗ , Valeriy Trubitsyn b a b

Deutsches GeoForschungsZentrum, Telegrafenberg, 14473 Potsdam, Germany Institute of Physics of the Earth, Moscow, Russia

a r t i c l e

i n f o

Article history: Received 3 November 2010 Received in revised form 21 February 2012 Accepted 23 February 2012 Available online 3 March 2012 Keywords: Geoid Mantle transition zone Mantle convection Clapeyron slope Viscosity Tomography

a b s t r a c t Joint inversion of the observed geoid and seismic velocities has been commonly used to constrain mantle properties and convection flow. However, most of the employed tomography models have been obtained without considering the effect of the mantle transition zone (TZ). We use the new-generation tomography model of Gu et al. (2003), in which seismic velocity perturbations are estimated together with variations of the main TZ boundaries. In the inversion of this model with the observed geoid the velocity-to-density scaling factor and density jump at these discontinuities are determined simultaneously. By this, the mantle flow across TZ is defined self-consistently: the undulations of the TZ boundaries suppress or accelerate mantle currents depending on the determined density contrast. For the 410-km discontinuity we obtain the scaling factor and density jump, which are close to mineral physics predictions. Therefore, we conclude that these effects are decoupled in the tomography model. In contrast, the calculated density jump at the 660 discontinuity is approximately 4 times less than the PREM value. We suggest that this is the effect of multiple phase transitions within a depth range of 640–720 km. Under normal thermal conditions, the post-spinel phase transformation is relatively sharp (∼5 km) and clearly visible in seismic models. On the other hand, the transition of majorite garnet to perovskite is much broader (up to ∼50 km) and, therefore hardly detectable. Due to the different sign of the Clapeyron slope, the total gravity effect is drastically decreased. To fit the obtained results, the required value of the Clapeyron slope for the majorite garnet to perovskite transformation should be equal to about +1.7 MPa/K. In the cold zones the same effect might be produced by the transition from ilmenite to perovskite at 610–640 km depth, which is in agreement with multiple reflections revealed in regional seismic studies near the bottom of TZ (Deuss et al., 2006). The estimated amplitude of the mantle flow across TZ is about ±20 mm/year, which corresponds to the whole-mantle convection scheme. The calculated geoid better fits to the observed one than the obtained without considering the TZ effect. © 2012 Elsevier Ltd. All rights reserved.

1. Introduction The observed geoid is one of the most important constraints on mantle parameters, first of all, density distribution and viscosity variations. However, determination of the Earth’s mantle structure has an ambiguous solution when using only surface gravity data. A usual way to cope with such a problem is to combine gravity data with other geophysical data sets to produce the solution, which fits all the data sets and therefore contains fewer degrees of freedom. Seismic tomography models are used for these purposes most often. Since only several parameters (usually depth-dependent velocity-to-density scaling factors and radial viscosity variations) are determined from a large volume of the data, the solution is relatively stable. In global studies this technique is applied starting from the pioneer works by Hager and O’Connell (1981), Ricard et al.

∗ Corresponding author. Tel.: +49 331 288 1172; fax: +49 331 288 1111. E-mail address: [email protected] (M.K. Kaban). 0264-3707/$ – see front matter © 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.jog.2012.02.007

(1984), Richards and Hager (1984), Forte and Peltier (1987), King and Masters (1992), Corrieu et al. (1994) and many others. Despite of the massive efforts, the obtained results are remarkably various and a generalized instantaneous model of the dynamic mantle does not exist at the moment. There might be several explanations for such an indefinite situation, but one of the principal factors is that seismic tomography models, which were analysed in the previous studies, do not include information about the mantle transition zone discontinuities. Up to now, most of the combined gravity-seismic models are based on the tomography data, in which the effects of velocity variations and phase boundaries are mixed (e.g. Forte and Perry, 2000; Steinberger and Calderwood, 2006 and many others). As it is shown in the following chapter, the effect of the phase boundary on seismic velocities depends on many factors, which are not well-defined. The phase boundary might even lead to strong artificial velocity in the vicinity of the phase transition. Inclusion of the phase boundaries in the model significantly influence the convection pattern and by this – the dynamic geoid. Defraigne and Wahr (1991) have developed a numerical technique

184

M.K. Kaban, V. Trubitsyn / Journal of Geodynamics 59–60 (2012) 183–192

to incorporate the effect of phase-boundaries into internal-loading dynamic models. It has been demonstrated that the effect of the phase boundaries on the dynamic geoid might be very significant (Defraigne et al., 1996). Moreover, these authors have found that using a density jump for the 660 discontinuity, which is prescribed in PREM (Dziewonski and Anderson, 1981) would lead to much larger impact to the observed geoid than it might be expected from the observed one. In this paper we use the new-generation tomography model of Gu et al. (2003), in which mantle velocities have been estimated in a joint inversion with the mantle transition zone (TZ) discontinuities (410 and 660 km). Employing joint inversion of all the parameters (seismic velocity anomalies, topography of the TZ discontinuities and observed geoid) we try to understand density structure of the TZ and to determine its impact to the dynamic geoid. A role of the TZ in shaping mantle convection has been discussed for decades. However, up to now the problem of the convection style (a rate of the mass exchange across the TZ) is not resolved completely. In the global studies, dealing with a joint inversion of the geoid and seismic tomography, the style of mantle convection was usually predefined. Most previous works starting from (Hager and O’Connell, 1981) employing the whole-mantle model because this model generally provides better fit to the observed geoid (e.g. Corrieu et al., 1994; Thoraval and Richards, 1997). However, the layered model has been analysed in several studies (Wen and Anderson, 1997; Forte and Woodward, 1997; Cadek et al., 1997; Pari and Peltier, 1998). It has been demonstrated that this model can also provide a reasonable fit to the observed geoid and ˇ and even better explanation for the dynamic topography. Cadek Fleitout (1999) have developed a combined approach, according to which the whole-mantle and layered convection models are mixed in various proportions. Despite this model is artificial and simply depicts a fixed rate of the 660 permeability, which is the same over the whole Earth for all flow patterns, it may help improving the fit to the observed geoid. Steinberger (2007) used another approach to study the effect of the phase boundaries on the geoid; he estimated the topography implied due to temperature variations inferred from seismic velocity variations and Clapeyron slope. The main problem for most of the above models is that the calculated topography of the 660 discontinuity is remarkably different ˇ from observations (e.g. Cadek and Fleitout, 1999; Steinberger and Torsvik, 2007; Steinberger, 2007). Therefore, implementation of seismically determined variations of the main TZ boundaries into geodynamic models is likely the most promising way to overcome the existing problems. Unfortunately existing seismic models of these boundaries are remarkably different (Steinberger, 2007). In contrast to previous models, the TZ topography model of Gu et al. (2003) has been obtained in a joint inversion with the tomography model. Therefore it should be more consistent. These results open innovative possibilities to investigate the role of the TZ boundaries in mantle dynamics. In this case we do not have to identify the convection style explicitly. We apply a joint inversion of the seismic velocities and 410, 660 topography with the observed geoid to determine simultaneously the velocityto-density scaling factor and density jump at the TZ discontinuities. After construction of the mantle density model, the TZ boundaries should control the convection style self-consistently suppressing or intensifying mantle flow through the boundary.

2. Velocity-to-density scaling in the vicinity of a phase boundary The mantle transition zone is extensively studied by seismic methods. Usually two techniques are employed for these purposes. Study of precursors to the SS and PP waves are widely used to

determine global variations of the TZ boundaries (e.g. Shearer, 2000; Deuss, 2009). The 660 discontinuity is the most important for a global dynamic modelling. However, existing studies based on available date often provided controversial results, which could not be explained using a simple mineralogical model of 660 discontinuity as described in PREM (Deuss et al., 2006; Deuss, 2009). Normally, the 660 discontinuity is globally detected by the longperiod S660S precursors, but the long-term P660P precursors are available only for some places (e.g. Estabrook and Kind, 1996; Shearer and Flanagan, 1999). On the other hand, the short-period observations show that the 660-km discontinuity is sharp and prominent in some isolated places, where it is not visible in the long-period PP precursors (Benz and Vidale, 1993). Furthermore, studies of the converted waves often provide the results, which are different from the precursors studies (e.g. Vinnik et al., 1997). All together, these results evidence that the 660 discontinuity is characterized by a complicated structure and is often represented by several boundaries in a depth range from 640 km to 720 km. It has been suggested that this structure might be related by multiple phase transitions in this depth range (e.g. Deuss et al., 2006). In addition to the commonly accepted transition from ringwoodite to perovskite and magnesiowüstite there should be additional boundaries. These features would inevitably influence the mantle flow and dynamic geoid. However, up to now most studies of the global dynamic model of the Earth were based on standard velocity models, which are constructed using body and surface waves tomography methods. It is usually assumed that these models also reflect principal features of the mantle transition zone providing artificial seismic velocity anomalies, which might be properly scaled to density anomalies used in the dynamic modelling. One of the principal parameters, determined in the inversion, is a depth dependant velocity-to-density scaling factor (SF), which defines 3D density structure of the mantle based on seismic models. Let us determine first the “effective” scaling factor (seffect ) for the traditional seismic models, in which the effects of in situ velocity variations and TZ topography are mixed. We assume that temperature variations (T) near any single phase boundary control variations of the material properties: seismic velocity (V), density () and topography of the TZ (h):  = ˛T, 

˛ V = T, V s

h = h T,

(1)

where ˛ is the thermal expansion coefficient, s = ∂ ln()/∂ ln(V) is the in situ velocity-to-density scaling factor (SF) and  h = dh/dT = −−1 g−1 dp/dT the Clapeyron slope for depth. The resolution of seismic tomography models is always limited; therefore it provides velocities, which characterize some depth range (Z) depending on the parameterization, damping and data coverage. Normally it is about 100–150 km at the depths 400–700 km (e.g. Gu et al., 2003). As a result, the variations of seismic velocity in the vicinity of a phase boundary represent a cumulative effect of in situ velocity changes and variations of the boundary dividing layers with basically different velocities. This effect might be described with sufficient accuracy by the simple equation: ˛ Veffect h ˛  T = T + Vtz = T + h Vtz , V s VZ s VZ

(2)

where Vtz = V2 − V1 (Fig. 1). Therefore, the average velocity variations near TZ depend on many factors, which could be variable and not defined precisely. In the case of the negative Clapeyron slope, the effective velocity variations might be of different sign for the same temperature anomaly or close to zero. The situation is even more complicated when seismic velocities and TZ boundary variations are not only controlled by temperature variations but also by other factors (composition, water content, etc.).

M.K. Kaban, V. Trubitsyn / Journal of Geodynamics 59–60 (2012) 183–192

185

Fig. 1. Temperature dependent variations of material properties in the vicinity of a phase boundary.

The corresponding density variations are as follows: total h  T = ˛T + tz = ˛T + h tz  Z Z

(3)

Contrary to seismic velocities, the amplitude of the second term likely exceeds the amplitude of the first one and the effective density variations are basically correlated with temperature anomalies depending on  h . Therefore, velocity anomalies in the “standard” tomography models might not provide reliable information about density variations in the vicinity of the phase boundary with the negative Clapeyron slope. The expressions (2) and (3) define the “effective” scaling factor near the phase boundary: seffective =

total Vtotal /  V

(4)

Consequently, even in this simplified case the effective scaling factor in the vicinity of a phase boundary depends on many factors, which could vary from place to place depending on the model resolution and physical conditions. Several examples for the scaling factor depending on vertical resolution of the tomography model are shown in Fig. 2. The parameters are taken from PREM (Dziewonski and Anderson, 1981), ˛ = 1.2 × 10−5 (solid line) or 2 × 10−5 (dashed line). These curves will be similar when we proportionally reduce the velocity and density jumps.

Fig. 2. Effective velocity-to-density scaling factor (seffect , Eq. (4)) in the vicinity of the phase boundary depending on the depth-resolution of the tomography model.  = 2.3 Mpa/K for the 410 discontinuity and −2 Mpa/K for the 660, the in situ SF s = 0.24, other parameters are taken from PREM (Dziewonski and Anderson, 1981).

Fig. 3. Power spectra of the 660 variations: (Gu et al., 2003, solid line) vs (Houser et al., 2008, dashed line).

For the 410 discontinuity (as well as for any phase boundary with positive Clapeyron slope) the effective scaling factor exceeds the in situ value and gradually decreases while decreasing model resolution. By contrast, the effective scaling factor for the 660 discontinuity might be substantially different (negative or positive). Within the resolution range, which is typical for existing tomography models, it could be also close to zero (Fig. 2), therefore variations of seismic velocities at this depth do not provide information about variations of the phase boundary. The only reliable solution to cope with the problem is to use the tomography models, in which velocity anomalies are determined together with the transition boundaries topography. 3. Data and inversion technique 3.1. Seismic tomography model Modelling of three-dimensional (3D) velocity variations in the mantle and topography of the TZ discontinuities have been considered separately in most previous studies. Gu et al. (2003) have developed an integrated approach; they used a set of S-velocity sensitive data as used before (e.g. Ekström and Dziewonski, 1998) and combined it with a large set of differential travel times of SS-S400S, S-S670S, and direct measurements of S400S–S670S. As a result, S-wave velocity variations (model TOPOS362D1) are obtained consistently with topography of the 410 and 660 boundaries. Although, correlation between this model and that one, obtained without considering the 410 and 660 discontinuities, is good for the whole mantle on the average; within the TZ it is lower (Gu et al., 2003). The new combined model shows much larger velocity perturbation at depth of 660 km than the previous S20a (Ekström and Dziewonski, 1998), which difference might be also related to the presence of the phase boundary with negative Clapeyron slope. A similar approach to determine velocity perturbations and depth variations of 410 and 660 has been used by Houser et al. (2008). However, their interpretation technique, initial data, parameterization and inversion parameters are somewhat different, which leads to different results. The power spectra of the 660 variations, which is the most important parameter for our model, are shown for both models in Fig. 3. The most prominent feature of the model of Houser et al. (2008) is a gap in the spectrum at

186

M.K. Kaban, V. Trubitsyn / Journal of Geodynamics 59–60 (2012) 183–192

the harmonics 3–5. Therefore, regional features dominate in this model, while the global ones are significantly reduced. At the same time, the observed geoid shows maximum power at these wavelengths. Source of the gap in the seismic model at these wavelengths is still unclear. Therefore, in this work we still use the combined model of Gu et al. (2003) in a joint inversion with the geoid anomalies to obtain the most appropriate velocity-to-density scaling factors and density jumps at the 410 and 660 boundaries. A comparative analysis of all existing tomography models will be a subject of a following study.

3.2. Geoid data We use EIGEN-GL04C, which is based on combination of the satellites CHAMP and GRACE and terrestrial data (Förste et al., 2008). The model is complete to degree/order 360 in terms of spherical harmonic coefficients and resolves features of about 55 km width in the geoid and gravity anomaly fields. In the wavelength range used in this study, the GL04C model is very close to the highresolution EGM2008 model (Pavlis et al., 2008). The long-wavelength features of the dynamic geoid contain not only the gravitational signal from deep-seated lateral mass and density inhomogeneities sustained by dynamic Earth mantle processes. To interpret the observed geoid with respect to mantle dynamics and structures, it is desirable first to remove the lithosphere-induced anomalous gravitational potential, which is generated by the topographic surface load and its isostatically compensating masses. It is also important that density variations in the uppermost mantle are to a large extent controlled by compositional variations, which are almost not reflected in the tomography models (e.g. Forte and Perry, 2000; Kaban et al., 2003). Based upon the most recent global compilation of crustal thickness and density data and the age distribution of cooling oceanic lithosphere, several global models of the crust and upper mantle have been developed (Panasyuk and Hager, 2000b; Kaban et al., 1999, 2004). The constructed isostatic models of the lithosphere are supposed to be valid for spatial wavelengths more than 500 km. The field, induced by these models, is subtracted from the observed field to yield the isostatic residual geoid. This residual field is chiefly controlled by the mantle structure below ∼200–250 km and corresponding dynamic topography and, therefore is most suitable for global dynamic modelling (e.g. Panasyuk and Hager, 2000b). In this study we use the model of Kaban et al. (2004). Applying the isostatic correction, the overall pattern of the geoid becomes smoother and the most pronounced features, which are separated in the observed geoid, tend to get connected to larger structures. Despite the overall effect of the isostatic correction is not so significant to revise previous results, in active tectonic areas the isostatic geoid reduction ranges from −21 to +43 m, which should improve the dynamic model. In this study we use the non-isostatic residual geoid, which does not comprise the terms C20 and C40 . Origin of these terms in the observed geoid relative to a hydrostatic spheroid is not yet completely understood (e.g. Nakiboglu, 1982; Mound et al., 2003). Furthermore, the terms C20 and C40 in seismic tomography are determined less reliably than other terms because they depend on the data in the polar areas, which are not well resolved. However, these terms dominate in the observed non-hydrostatic geoid, their globally estimated root mean square (RMS) is equal to 28.5 m, which is almost identical to the RMS of all other terms (31.8 m). Therefore, their uncertainties, especially of C20 , might strongly affect the result. All the models are computed including all the terms, but C20 and C40 were not taken into account when fitting the calculated geoid to the observed one. These terms are also excluded from all figures.

3.3. Inversion technique As starting from Hager and O’Connell (1981), it is assumed that since seismic velocity variations are small, they might be related to density variation via a linear depth-dependent scaling factor (SF). The observed geoid variations are produced by the density variations within the mantle and by the dynamic disturbances of the Earth surface caused by mantle currents. On the other hand, the mantle currents are also controlled by the density variations in the mantle and can be modelled using Stokes equations with a specified viscosity distribution within the mantle. Therefore, assuming some viscosity model, one can determine by inversion the depth-dependent SF, which provides best fit of the observed and calculated geoid. In case of only radial viscosity variant, the Stokes equations (continuity and momentum) together with Poisson’s equation (gravity field flux) including effects of compressibility, self-gravitation and depth-dependent gravity can be solved by a direct method of solving the ODE system for each spherical harmonic mode, which gives a very fast and efficient way to solve the inverse problem. Importance of lateral viscosity variations is still ˇ under debate (e.g. Cadek and Fleitout, 2003; Kaban et al., 2007; Rogozhina, 2008) and we do not consider this effect in the present paper. The inversion procedure is standard. The initial density model is produced by a simple linear conversion of seismic velocity disturbances into density variations using the initial scaling factors ı = 0 a0j (ıVs/Vs). The dynamic geoid is then calculated for separate mantle layers assuming a constant velocity-to-density scaling factor for the given viscosity model. In this study, we chose the initial scaling factor a0j = 0.24, which is consistent with mineral physics results (e.g. Karato, 1993). The initial scaling factors are rescaled in a least-squares adjustment to get the best fit to the observed geoid, therefore to minimize the function 2 .

2 =



⎛ ⎝N l,n − obs

l,n

+



⎞2 l,n l,n ⎠ aj Njl,n − a410 N410 − a660 N660

j



ˇj (aj − 1)2

(5)

j

Nobs is the observed geoid for degree/order l and n; Nj are the geoid variations induced by the layer j; aj are the unknown corrections for the initial scaling factor a0j , and ˇj are the damping factors introduced to stabilize the solution (Kaban and Schwintzer, 2001). The second term is added to stabilize a solution, as usual with a prospect not to wander away too far from the initially chosen scaling factors, which are consistent with the mineral physics predictions. The inversion to solve for unknown scaling factors is performed in the spectral domain by spherical harmonics coefficients. In this study the traditional scheme is modified by adding l,n l,n the terms responsible for the TZ topography a410 N410 and a660 N660 (Eq. (5)). These boundaries are taken directly from the seismic model. In the initial stage we calculate the effects of the 410 and 660 discontinuities assuming the standard density jumps from PREM: 0.18 g/cm3 and 0.39 g/cm3 correspondingly (Dziewonski and Anderson, 1981). Then they are also rescaled in the inversion. The impact of the phase boundaries is fully dynamical; their variations suppress or accelerate mantle flow. Then, we adjust these initial values by calculating the additional scaling factors (a410 and a660 , Eq. (5)) in the inversion together with the scaling factors for velocity perturbations. We do not apply damping for these coefficients; therefore the initial values do not affect the solution. Minimization of the function 2 leads to the normal equation system Bâ = c. The solved-for parameters aˆ j are used then for

M.K. Kaban, V. Trubitsyn / Journal of Geodynamics 59–60 (2012) 183–192

correction of the initial scaling factors:

 d ln   d ln Vs

j

= aˆ j a0j ,

0 410 = a410 410 ,

0 660 = a660 660

(6) The scaling coefficients standard deviations saˆ j resulting from the fit in the least squares adjustment are then:



saˆ j =

1 2  qjj , f

(7)

where qjj are diagonal elements of B−1 , f – the degree of freedom. More details about the inversion techniques is found in Kaban and Schwintzer (2001) and Kaban et al. (2004). 4. Results The whole mantle has been divided into 50 km thick layers starting from the depth 225 km. These layers are used to calculate the dynamic geoid for specified radial viscosity profiles. The crust and upper mantle have been excluded from the inversion since their effect is already removed from the observed geoid as explained in (Section 3.2). The velocity-to-density conversion coefficients are determined for less numbers of depth intervals, which correspond to the basic radial splines used in the tomography model (Gu et al., 2003). Furthermore, in the lower mantle we determine only 4 independent scaling factors to reduce their determination error (Fig. 4). The initial scaling factor is 0.24. Viscosity is specified separately for the layers limited by the following depths: 0, 75, 225, 410, 550, 660, 1600, 2600 and 2900 km. These depth intervals correspond to principal viscosity changes found in previous studies (e.g. Corrieu et al., 1994; Panasyuk and Hager, 2000a). The reference viscosity (1021 Pa s) is set for the depth interval 225–410 km, in other layers it is varied within the limits (Fig. 5), which are in agreement with mineral physics predictions and previous studies (e.g. Panasyuk and Hager, 2000a; Steinberger and Calderwood, 2006). Each viscosity interval has been divided by several points with values changing by a factor of two. Then, the inversion is performed for all viscosity profiles (∼100,000, grid search), and the best results are demonstrated below. 4.1. Velocity-to-density scaling factor and density contrast at the 410 and 660 discontinuities The best solutions for the velocity-to-density scaling factor and density contrast for the 410 and 660 discontinuities are shown in Fig. 4. For comparison we show the SF for the previous S20a model (Ekström and Dziewonski, 1998). This model is based on nearly the same data set but without considering the TZ boundaries (see above). The geoid–tomography inversion is performed without considering the 410 and 660 discontinuities; however other parameters are the same. The scaling factor is less stable in this case; it varies from zero at the bottom of the upper mantle to 0.32 in the layer just below it in the vicinity of the 410 discontinuity. In the lower mantle the SF for the old model systematically exceeds the SF for the new model with the TZ discontinuities. The geoid kernels are very low for the depth interval 1000–1300 km; therefore the determination error for this interval is large. The determined velocity-to-density SF for the combined tomography model is on the average (∼0.2) slightly less than the initial value (0.24). It reliably exceeds this level only inside TZ (500–660 km), where the SF reaches 0.3. This result is consistent with the previous conclusion that the layer containing a phase boundary with the positive Clapeyron slope should be characterized by an increased SF. The transition between wadsleyite and

187

ringwoodite (∼520 km) is not considered explicitly in the tomography model and should affect the SF (Fig. 2). At the same time, the SF, determined for the new combined model, does not differ from the average value at 400 km and 660 km. In contrast, the SF, obtained for the old S20a model, demonstrates a tendency, which is consistent with the previous conclusion. It is larger than the average around 400 km and lesser at 660–1000, which agrees in general with the above conclusions (Section 2, Fig. 2). However, these differences are less than it might be expected from the pure PREM model (Fig. 2). This may be explained by the effect of “damping” and by more complicated structure of the TZ than it is presumed in PREM. The last factor is especially important for the 620–700 km zone (see Section 5). The effective density contrast, determined for the 410 and 660 boundaries, is shown in Fig. 4B. It should be emphasized that these values characterize variations of the discontinuities but not the absolute value of the density change at those depths, since the zero term is not considered. This means that variations of any boundary might be associated with corresponding lateral density changes, which compensate the pure effect of the depth changes but do not affect the absolute density contrast. Similar effect would appear in case of a nearby boundary with opposite variations, which is not visible in the seismic model. As a result, the effective density contrast, which is attributed to a single interface, would be reduced. We have found that the effective density contrast for the 410 discontinuity is close to the PREM value (0.18 g/cm3 ), although the determination error is relatively large. Quite the opposite, the density contrast for the 660 discontinuity is unusually small (0.1 g/cm3 ), and this value is sufficiently reliable. Even, it could represent an upper boundary, because the amplitude of this boundary is likely reduced due to inevitable damping of the seismic inversion. Therefore, the effective density contrast, which characterizes variations of the 660 discontinuity, is four times less than the PREM value (0.39 g/cm3 ). This result agrees with the qualitative estimates of Defraigne et al. (1996) for potential impact of the 660 boundary to the dynamic geoid. Below we discuss possible explanations. 4.2. Viscosity As it is described above, the inversion was performed for different radial viscosity models. We select the model, which is characterized by both minimal residual geoid and perturbations of the scaling factor. The obtained radial viscosity profile (Fig. 5) is generally in agreement with previous studies (e.g. Panasyuk and Hager, 2000a; Corrieu et al., 1994). The average viscosity of the lower mantle is 30 times larger than the viscosity at a depth of 200–300 km, which is taken as a reference value. Different from the results of Forte and Perry (2000) and Forte and Mitrovica (2001) we did not find that strong viscosity variations in the lower mantle may remarkably improve the solution, however we did not investigate this layer in details. One of the principal features of the model is the low-viscosity zone right above the 660 discontinuity, so called “notch”, which is required to provide best fit between the observed and calculated geoid. Before, this feature was typically recognized by postglacial rebound studies (e.g. Milne et al., 1998). In fact, it serves like a barrier between the upper and lower mantle, which significantly reduces the dynamic topography induced by lower mantle heterogeneity and, as a result, affects the dynamic geoid. This effect is somewhat similar to a chemical barrier in the layered model. However, in the studies employing a joint inversion of the geoid and seismic tomography data, identification of the “notch” is uncertain. For example, Panasyuk and Hager (2000a) have found this feature employing a whole-mantle convection model. On the other hand, Lee et al. (2011) have demonstrated that it disappears if one considers the variable uncertainties of tomography models. Therefore,

188

M.K. Kaban, V. Trubitsyn / Journal of Geodynamics 59–60 (2012) 183–192

Fig. 4. (A): Velocity-to-density scaling factors determined in the joint inversion of the non-isostatic geoid and S-wave velocity variations: (solid) for the combined model (Gu et al., 2003) including the TZ boundaries; (dashed) for the previous S20a model (Ekström and Dziewonski, 1998). (B): Determined density contrast across the 410 and 660 discontinuities, which provide best fit to the geoid.

it is important to investigate if this feature appears in the combined model employing seismically determined variations of the main TZ boundaries. The misfit of the calculated geoid and perturbations to the initial scaling factor are shown in Fig. 5B for various viscosity models. It is demonstrated that the best fit to the observed geoid is provided only by the models containing the low-viscosity zone above the 660 discontinuity. Therefore, we may conclude that the low-viscous “notch” above 660 better explains principal features of the observed model. It should be noted that there exists a strong trade-off between thickness and amplitude of the viscosity low (Milne et al., 1998). Therefore, this conclusion should be considered as qualitative, the actual decrease of viscosity might be different depending on the thickness. 4.3. Dynamic geoid The observed and calculated geoid for different tomography models are shown in Fig. 6. All the inversion parameters (damping factor and initial scaling) are the same. Using of the combined model with the TZ boundaries provides better fit of the observed and calculated fields. A variance reduction for geoid is 81% for the combined tomography model versus 78% for the model without the TZ boundaries. Although these values are similar to previously determined (e.g. Steinberger and Calderwood, 2006), it is not correct to compare them, because in the present study we do not consider the terms C20 and C40 , which provide more than a half to

the geoid power. Also we see visible improvements for regional features. The negative anomalies near Antarctica appear only if we take into consideration these boundaries. Shape of the positive anomalies in South Atlantic and Pacific are also reproduced much better in this case. We conclude that considering of the TZ boundaries is essential for modelling of the dynamic geoid. The impact of the 410 discontinuity reaches ±10 m, while the effect of the 660 discontinuity on the geoid is even more pronounced and exceeds ±20 m despite the density contrast is 4 times reduced relative to the PREM value (Fig. 7). This field provides a substantial input to the total geoid (Fig. 6). The dynamic geoid represents a sum of the effects of mantle density heterogeneity and dynamic topography induced by the convecting mantle. The dynamic topography is often compared with the “residual topography” obtained by removing effects of the crust and cooling ocean lithosphere from the observed topography. The residual topography represents that part of the observed topography, which is not compensated or over-compensated by the crustal structure and is not consistent with the standard cooling model of the oceanic plate (Kaban et al., 1999, 2004; Panasyuk and Hager, 2000b). We do not find a good correspondence between the residual and calculated dynamic topography, which agrees with the results of Steinberger (2007). However, the residual topography still depends chiefly on unknown density anomalies within the lithosphere, for example – dense lithospheric roots of the continents; and the dynamic effects likely provide less input to this

Fig. 5. (A) Normalized viscosity profile, which corresponds to the best solution. Thin solid lines depict possible variations of viscosity at different depths, which provide reasonable fit to the geoid and initial scaling factor (circle in the right). Dashed lines show allowable limits of the viscosity changes. (B) Residual geoid vs. perturbations to the initial scaling factor (0.24) for the analysed viscosity models. Different colours correspond to 4 clusters with different viscosity above the 660 discontinuity (550–660 km). The big circle shows the area of the best solutions. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

M.K. Kaban, V. Trubitsyn / Journal of Geodynamics 59–60 (2012) 183–192

189

Fig. 6. (A) Observed geoid with the effect of the crust and upper mantle removed. (B) Geoid obtained with the s20 s model without TZ boundaries (Ekström and Dziewonski, 1998). (C) Geoid based on the joint inversion with the combined model (Gu et al., 2003) including the 410 and 660 boundaries. Terms C20 and C40 are excluded from all fields.

parameter (Kaban et al., 2004). In our opinion, at this stage we cannot use this parameter to constrain global dynamic models of the whole mantle.

the Kuril and Marianas Arcs. These downwellings might be qualitatively associated with the slabs penetrating the 660 discontinuity. For other arc-trench systems we do not see such effect, but these findings are very preliminary.

4.4. Mantle flow across the 660 discontinuity 5. Discussion With the determined density structure of the mantle (via velocity-to-density scaling) and density contrasts at the main boundaries, the resulting mantle flow pattern appears selfconsistently, the observed undulations of the 410 and 660 boundaries control vertical mantle currents (suppress or accelerate them). The calculated mantle flow across the 660 discontinuity is significant (about ±20 mm/year), Fig. 8. It radically exceeds the ˇ and Fleitout (1999) for amplitude of the flow determined by Cadek the mixed model (about ±6 mm/year), where the mantle flow is reduced to one third in comparison with the purely whole-mantle model. Therefore, the combined model estimated in a joint inversion with the TZ boundaries remains essentially whole-mantle. The estimated flow does not correlate in general with the topography of this boundary (Fig. 8B). Therefore, it is controlled by the global convection pattern and we conclude that the 660 zone does not influence the style of convection remarkably. In further discussion we attempt to explain this phenomenon. Resolution of the tomography model is not sufficient to recognize behaviour of subducting slabs, but some assumptions can already be made from Fig. 8. We see strong downwellings under South America and east of Australia behind the Tonga–Kermadec arc. Less but still remarkable descending flow is recognized behind

5.1. Mineralogical density model of the transition zone 410–660 The commonly used pyrolite model of the mantle contains 60% olivine and 40% non-olivine residuum, chiefly garnet in the mantle transition zone. The phase transitions in olivine up to now were considered as dominating in geodynamics (e.g. Schubert et al., 2001). Below a depth of about 410 km olivine transforms to wadsleyite (ı/ ≈ 7%) and, at about 520 km depth, wadsleyite transforms into ringwoodite, which has the spinel structure with a density jump ı/ ≈ 3%. Both these transformations have positive Clapeyron slope,  410 ≈ 1.6–3.0 MPa/K and  520 ≈ 4.3 MPa/K correspondingly and accelerate thermal mantle convection by about 2–3% (Schubert et al., 2001). The transition from olivine to wadsleyite is sharp (5–10 km), therefore, it is clearly detected by seismic methods. By contrast, the transition from wadsleyite to ringwoodite might be spread up to ∼60 km, which leads to difficulties of its detection (e.g. Deuss, 2009). However, the most important transition for mantle convection is the post-spinel transition in olivine from ringwoodite to perovskite and magnesiowüstite at a depth of about 660 km. This transition has been considered for a long time as dominating

Fig. 7. Calculated geoid: (A) effect of the 410 discontinuity; (B) effect of the 660 boundary. Both fields are computed for a fully dynamical model including their effect on the dynamic topography.

190

M.K. Kaban, V. Trubitsyn / Journal of Geodynamics 59–60 (2012) 183–192

Fig. 8. (A) Calculated mantle flow across the 660 discontinuity (red – upwellings, blue – downwellings). (B) Variations of the 660 discontinuity (Gu et al., 2003). The field is truncated after the 20th spherical harmonic. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

at the bottom of the TZ and the corresponding density change ı/ ≈ 9% has been taken for PREM as the only change at this depth (Dziewonski and Anderson, 1981). The transition from ringwoodite to perovskite is characterized by a negative Clapeyron slope and therefore should suppress mantle currents. For a long time it was assumed  660 ≈ −2.5 MPa/K (Schubert et al., 2001). New laboratory studies of the dry olivine evidence that the Clapeyron slope might not be so steep; Katsura et al. (2005) give  660 ≈ −0.6 to 1.2 MPa/K and Litasov and Ohtani (2005)  660 ≈ −0.4 to 1.3 MPa/K. In this case the effect on the mantle convection would be lower. However, further experiments have demonstrated that the material in the transition zone 410–660 can accumulate much more water than the material of the upper and lower mantle. In the wet olivine the Clapeyron slope of the transition Rw–Pv is higher and the former value  660 ≈ −2 MPa/K is still considered as most appropriate (Ohtani and Litasov, 2006). Despite the non-olivine phase change are well known in mineral physics, they are often ignored in seismic and geodynamic studies and the olivine transformations are considered being dominated. However, real structure of the mantle transition zone is more complicated. Phase diagrams for a pyrolitic mantle at 600- to 750km depth after Akaogi et al. (2002) are shown in Fig. 9. In case of the “normal” geotherm (mean mantle temperature ∼1600 ◦ C) ringwoodite transforms to perovskite and magnesiowüstite over a narrow depth interval (∼5 km), the density change corresponding to the 60% content is about ı2 = 0.6·ır-pv = 4.4% and the Clapeyron slope is negative. At the same time the transition of the non-olivine component garnet–perovskite is spread in a wider depth range (∼660–710 km), the density change induced by this transition is about ı3 = 0.4·ıg-pv = 3.8% while the Clapeyron slope is positive (Akaogi et al., 2002; Akaogi, 2007). Therefore, the total density change (4.4 + 3.8 = 8.2%) is close to PREM (9%) but it is distributed between 2 transitions with the opposite sign of the Clapeyron slope (Akaogi et al., 2002).

At lower temperatures, which might relate to a downgoing slab, ringwoodite transforms to perovskite and magnesiowüstite at a greater depth (∼690–700 km). Garnet first transforms to ilmenite (at 610–670 km) and ilmenite transforms to perovskite at a depth range of about 690–730 km, the density changes induced by this transition are 0.4·ıg-il = 3.5% and 0.4·ıil-pv = 2.5% correspondingly, which gives the overall change·ıg-pv = 6% (Akaogi et al., 2002). Then, the total density difference (ringwodite and garnet to perovskite) at low mantle temperatures at the depth range 610–720 km should be equal to: ı660 = 4.4 + 6.0 = 10.4%. Taking into account that all the parameters are not well-defined, the average value of the density change is also close to PREM. The evidences for multiple phase transitions have been already revealed by seismic studies. Shearer and Flanagan (1999) using detailed measurement across the TZ have located a relatively narrow transition from rindwoodite to perovskite (5%). Deuss et al. (2006) have found a complicated structure, showing single and double reflections ranging in depth from 640 to 720 km, which might be explained by the existence of multiple phase transitions at the base of the transition zone 410–660. 5.2. Effect of the multiple phase transitions on the dynamic geoid Defraigne et al. (1996) have already demonstrated that a potential effect of the 660 discontinuity might correspond to a lesser density contrast than it was supposed before. Using the new data it is possible to estimate reliably this value. We found that the effective contrast should be four times less than the PREM value. This result might be explained if we take into account the multiple phase transitions within a depth range 640–720 km. A schematic image of the phase transformations in this depth range is shown in Fig. 10 (normal geotherm). Increased temperature leads to an uplift of the relatively narrow boundary related to the transformation from ringwoodite to perovskite and magnesiowüstite (red).

Fig. 9. Mineral proportions vs. depth in the mantle (Akaogi et al., 2002). (A) The normal mantle geotherm; and (B) the low-temperature geotherm.

M.K. Kaban, V. Trubitsyn / Journal of Geodynamics 59–60 (2012) 183–192

Fig. 10. Variations of the main phase boundaries in case of increased temperature (see text). (For interpretation of the references to colour in the text, the reader is referred to the web version of the article.)

However, the closely located wide zone, which corresponds to the majorite garnet–perovskite transformation, is depressed due to the positive Clapeyron slope (blue). Therefore, the gravity effects of these variations partially compensate each other and the total effect is reduced. This delicate balance leads to much smaller “effective” contribution to the dynamic geoid. Using experimental data on the phase transformations’ properties we can estimate the “effective” density change, which corresponds to the cumulative effect of two boundaries in the 660 zone. Uplift of the Rw–Pv transition leads to the effective mass excess ıM2 = ı2 ıh2 . Its amplitude is controlled by the Clapeyron slope:  = dp/dT = −gdh/dT, therefore ıM2 = − 2 ı2 dT/g. The same is true for the Gar-Pv transition (ıM3 = − 3 ı3 dT/g) but the sign of  3 is opposite. Then we can find the total mass change due to variations of both phase transformations: ıM = ı2 ıh2 − ı3 ıh3 = (−2 ı2 − 3 ı3 )

dT . g

(8)

The relation of the effective mass balance to the value from PREM, which corresponds to only Rw–Pv transformation, is equal to: ıM −2 ı2 − 3 ı3 = ıM2 −2 ıPREM660

(9)

As it was discussed above, ı2 / = 4.4%, ı3 / = 3.8% (Akaogi et al., 2002), ıPREM660 = 9% and  2 = −2 MPa/K for wet ringwoodite (Ohtani and Litasov, 2006). The Clapeyron slope for the Gar-Pv transformation is not well defined and could be within a range of 1–3 MPa/K (e.g. Akaogi and Ito, 1999). Assuming  3 = 1 MPa/K we obtain ıM/ıM2 = 0.35. This means that the total gravity effect is decreased for about 3 time compared to the pure garnet–perovskite transition. However, for  3 = 3 MPa/K, the total gravity effect will be close to zero ıM/ıM2 = 0.07. To fit the obtained result (ıM/ıM2 = 0.25) we should take the effective slope for the Gar-Pv transformation  3 ≈ 1.7 MPa/K. Clearly, this evaluation is simplified. We do not take into account the complicated phase transitions along the cold geotherm (downwelling), where the garnet–perovskite transition splits into the garnet–ilmenite–peroskite transition. The considered example represents a “pure” case of the phase transformations without significant changes, e.g. of minor elements. In reality there exist various factors, which significantly modify the ideal scheme. The variations of the 660 discontinuity determined from S660S precursors do not always correlate with the velocity distribution at the same depth and also do not correlate in many cases with the variations of the 410 discontinuity. We have already mentioned possible effect of the water content, which may vary significantly from place to place. Another important factor is a kinetic effect, which is especially important in the subduction zones (Kubo et al.,

191

2004). Temperature within the slab is less than in the surrounding mantle, which might significantly slow the phase transformations and ringwoodite penetrates below 660 km to a greater depth. This is also visible in the model of Gu et al. (2003). The deepest position of the 660 discontinuity is found under the continental flank of the subduction zones (Fig. 8). Another contribution, which might affect the mantle transition zone, is a latent heat effect (e.g. Steinberger, 2007). Although, in the latter paper it was found that it likely only yields a minor contribution to phase boundary topography. Therefore, the characteristics of the discontinuity can no longer be taken as a global constant. But what is important for the gravity and geodynamic modelling, that independent from all the factors, there always exists a delicate balance between the boundaries of different origin and different properties located within a narrow depth interval. The equivalent density contrast is much lower than for the specified transformation taken alone. Therefore, this small “effective” contrast should be used in a geodynamic modelling. 6. Conclusions The new-generation tomography model of Gu et al. (2003), in which velocity perturbations are obtained simultaneously with undulations of the main TZ boundaries, has been used to investigate density structure of the transition zone 410–660 and to construct an instantaneous dynamic model of the mantle. Velocityto-density scaling factor, “effective” density contrast at the 410 and 660 discontinuities and radial viscosity distribution have been determined in a joint inversion with the observed geoid. The new model provides a possibility to quantify consistently a rate of the mass exchange across the mantle transition zone; and therefore – to define a style of the mantle convection. The 410 and 660 boundaries provide better similarity of the calculated and observed geoid relative to the previous tomography models. The calculated scaling factor and density jump for the 410km discontinuity are very close to the mineral physics prediction; therefore, we can conclude that these effects are really decoupled in the tomography model. The effective density contrast at the 660 boundary, which characterizes the depth variations, is equal to 0.1 g/cm3 , which is 4 times less than the PREM value. This result might be explained by the effect of multiple phase transitions with different Clapeyron slope. Assuming that two phase changes dominate within this depth range (1: ringwoodite to perovskite and magnesiowüstite and 2: garnet to perovskite) and taking into account properties of the first phase transformation ( = −2 MPa/K), we can estimate an average value of the Clapeyron slope for the second one. The corresponding value should be equal to +1.7 MPa/K. In agreement with previous studies, it has been confirmed that viscosity of the layer right above 660 km is remarkably less than of the neighbouring layers. Despite of significant uncertainty of many parameters of the phase transformations we can reliably conclude that the effective density jump at the 660 discontinuity, which corresponds to lateral variations of the seismically determined boundary, is much less then it was expected before. Therefore, the effect of this layer on the global mantle convection should be generally small. The estimated mantle flow across the 660 boundary has the amplitude of about ±20 mm/year, which corresponds to the whole-mantle convection scheme. Acknowledgements We would like to thank H. Schmeling, B. Steinberger and W. Bosch for very constructive reviews and advices, which helped to improve the paper. Funding for this project was kindly provided by

192

M.K. Kaban, V. Trubitsyn / Journal of Geodynamics 59–60 (2012) 183–192

DFG (German Research Foundation) RO-2330/4-II and KA 2669/2I (SPP 1257 “Mass Transport and Mass Distribution in the Earth System”). References Akaogi, M., 2007. Phase transitions of minerals in the transition zone and upper part of the lower mantle. In: Ohtani, E. (Ed.), Geol. Soc. Amer. Memoir. Advances in High Pressure Mineralogy. , pp. 1–13. Akaogi, M., Ito, E., 1999. Calorimetric study on majorite–perovskite transition in the system Mg4 Si4 O12 –Mg3 Al2 Si3 O12 : transition boundaries with positive pressure–temperature slopes. Phys. Earth Planet. Inter. 114, 129–140. Akaogi, M., Tanaka, A., Ito, E., 2002. Garnet–ilmenite–perovskite transitions in the system Mg4 Si4 O12 –Mg3 Al2 Si3 O12 at high pressures and high temperatures: phase equilibria calorimetry and implications for mantle structure. Phys. Earth Planet. Inter. 132, 303–324. Benz, H.M., Vidale, J.E., 1993. Sharpness of upper-mantle discontinuities determined from high-frequence reflections. Nature 365, 147–150. ˇ Cadek, O., Fleitout, L., 1999. A global geoid model with imposed plate velocities and partial latering. J. Geophys. Res. 104, 29055–29075. Cadek, O., Cizkova, H., Yuen, D.A., 1997. Can long-wavelength dynamical signatures be compatible with layered mantle convection? Geophys. Res. Lett. 24, 2091–2094. Corrieu, V., Ricard, Y., Froidevaux, C., 1994. Converting mantle tomography into mass anomalies to predict the Earth’s radial viscosity. Phys. Earth Planet. Inter. 84 (1–4), 3–13. Defraigne, P., Wahr, J.M., 1991. The response of a compressible, non-homogeneous Earth to internal loading theory. J. Geomagn. Geoelectr. 43, 157–178. Defraigne, P., Dehant, V., Wahr, M., 1996. Internal loading of an inhomogeneous compressible Earth with phase boundaries. Geophys. J. Int. 125, 173–192. Deuss, A., 2009. Global observations of mantle discontinuities using SS and PP precursors. Surv. Geophys. 30, 301–326. Deuss, A., Redfern, S., Chambers, K., Woodhouse, J., 2006. The nature of the 660kilometer discontinuity in Earth’s mantle from global seismic observations of PP precursors. Nature 311, 198–201. Dziewonski, A.M., Anderson, D.L., 1981. Preliminary reference Earth model. Phys. Earth Planet. Inter. 25, 297–356. Ekström, G., Dziewonski, A.M., 1998. The unique anisotropy of the Pacific upper mantle. Nature 394, 168–172. Estabrook, C.H., Kind, R., 1996. The nature of the 660-kilometer upper-mantle seismic discontinuity from precursors to the PP phase. Science 274, 1179–1182. Förste, C., Schmidt, R., Stubenvoll, R., Flechtner, F., Meyer, U., Koenig, R., Neumayer, H., Biancale, R., Lemoine, J.-M., Bruinsma, S., Loyer, S., Barthelmes, F., Esselborn, S., 2008. The GeoForschungsZentrum Potsdam/Group de Recherche de Gèodésie Spatiale satellite-only and combined gravity field models: EIGEN-GL04S1 and EIGEN-GL04C. J. Geodesy 82, 331–346. Forte, A.M., Mitrovica, J.X., 2001. Deep-mantle high-viscosity flow and thermochemical structure inferred from seismic and geodynamic data. Nature 410, 1049–1056. Forte, A.M., Peltier, W.R., 1987. Plate tectonics and aspherical Earth structure: the importance of poloidal–toroidal coupling. J. Geophys. Res. 92 (B5), 3645–3679. Forte, A.M., Perry, H.K.C., 2000. Geodynamic evidence for a chemically depleted continental tectonosphere. Science 290, 1940–1944. Forte, A.M., Woodward, R.L., 1997. Seismic–geodynamic constraints on threedimensional structure, vertical flow, and heat transfer in the mantle. J. Geophys. Res. 102, 17981–17994. Gu, Y.J., Dziewonski, A.M., Ekström, G., 2003. Simultaneous inversion for mantle shear velocity and topography. Geophys. J. Int. 154, 559–583. Hager, B.H., O’Connell, R.J., 1981. A simple global model of plate dynamics and mantle convection. J. Geophys. Res. 86, 4843–4867. Houser, C., Masters, G., Flanagan, M., Shearer, P., 2008. Determination and analysis of long-wavelength transition zone structure using SS precursors. Geophys. J. Int. 174, 178–194. Kaban, M.K., Schwintzer, P., 2001. Oceanic upper mantle structure from experimental scaling of Vs and density at different depths. Geophys. J. Int. 147 (1), 199–214. Kaban, M.K., Schwintzer, P., Tikhotsky, S.A., 1999. Global isostatic residual geoid and isostatic gravity anomalies. Geophys. J. Int. 136, 519–536. Kaban, M.K., Schwintzer, P., Artemieva, I.M., Mooney, W.D., 2003. Density of the continental roots: compositional and thermal contributions. Earth Planet. Sci. Lett. 209, 53–69.

Kaban, M.K., Schwintzer, P., Reigber, Ch., 2004. A new isostatic model of the lithosphere and gravity field. J. Geodesy 78 (6), 368–385. Kaban, M.K., Rogozhina, I., Trubitsyn, V., 2007. Importance of lateral viscosity variations in the whole mantle for modelling of the dynamic geoid and surface velocities. J. Geodesy 43 (2), 262–273. Karato, S.I., 1993. Importance of anelasticity in the interpretation of seismic tomography. Geophys. Res. Lett. 20, 1623–1626. Katsura, T., Yamada, H., Shinmei, T., Kubo, A., Ono, Sh., Kanzaki, M., Yoneda, A., Ito, E., Funakoshi, K., Utsumi, W., 2005. Postspinel transition in Mg2 SiO4 . Eos Trans. AGU 86 (52), Fall Meet. Suppl., Abstract MR44A-01. King, S.D., Masters, G., 1992. An inversion for radial viscosity structure using seismic tomography. Geophys. Res. Lett. 19 (15), 1551–1554. Kubo, T., Ohtani, E., Funakoshi, K., 2004. Nucleation and growth kinetics of the ␣␤ transformation in Mg2 SiO4 determined by in situ synchrotron powder X-ray diffraction. American Mineralogist 89 (2–3), 285–293. Lee, C.-K., Hana, S.-C., Steinberger, B., 2011. Influence of variable uncertainties in seismic tomography models on constraining mantle viscosity from geoid observations. Phys. Earth Planet. Inter. 184, 51–62. Litasov, K.D., Ohtani, E., 2005. Influence of water on the phase transitions of olivine and its polymorphs in the Earth’s mantle. Eos Trans. AGU 86 (52), Fall Meet. Suppl., Abstract MR44A-01. Milne, G.A., Mitrovica, J.X., Forte, A.M., 1998. The sensitivity of glacial isostatic adjustment predictions to a low-viscosity layer at the base of the upper mantle. Earth Planet. Sci. Lett. 154, 265–278. Mound, J.E., Mitrovica, J.X., Forte, A.M., 2003. The equilibrium form of a rotating earth with an elastic shell. Geophys. J. Int. 152, 237–241. Nakiboglu, S.M., 1982. Hydrostatic theory of the Earth and its mechanical implications. Phys. Earth Planet. Inter. 28 (4), 302–311. Ohtani, E., Litasov, K., 2006. The effect of water on mantle phase transitions. Rev. Mineral. Geochem. 62, 397–420. Panasyuk, S.V., Hager, B.H., 2000b. Models of isostatic and dynamic topography, geoid anomalies, and their uncertainties. J. Geophys. Res. 105 B, 28199–28209. Pari, G., Peltier, W.R., 1998. Global surface heat flux anomalies from seismic tomography-based models of mantle flow: Implications for mantle convection. J. Geophys. Res. 103, 23743–23780. Pavlis, N.K., Holmes, S.A., Kenyon, S.C., Factor, J.K., 2008. An earth gravitational model to degree 2160. In: EGM2008 European Geosciences Union General Assembly, Vienna, Austria, April 13–18. Ricard, Y., Fleitout, L., Froidevaux, C., 1984. Geoid heights and lithospheric stresses for a dynamic Earth. Ann. Geophys. 2, 267–286. Richards, M.A., Hager, B.H., 1984. Geoid anomalies in a dynamic Earth. J. Geophys. Res. 89 (B7), 5987–6002. Rogozhina I., 2008. Global modeling of the effect of strong lateral viscosity variations on dynamic geoid and mantle flow velocities. PhD Thesis. TU, Berlin, February 4. Schubert, G., Turcotte, D.L., Olson, P., 2001. Mantle Convection in the Earth and Planets. Cambridge Univ. Press, p. 940. Shearer, P.M., 2000. Upper mantle seismic discontinuities. Geophys. Monogr. 117, 115–131. Shearer, P.M., Flanagan, M.P., 1999. Seismic velocity and density jumps across the 410- and 660-kilometer discontinuities. Science 285, 1545–1548. Steinberger, B., 2007. Effects of latent heat release at phase boundaries on flow in the Earth’s mantle, phase boundary topography and dynamic topography at the Earth’s surface. Phys. Earth Planet Inter. 164 (1–2), 2–20. Steinberger, B., Calderwood, A.R., 2006. Models of large-scale viscous flow in the Earth’s mantle with constraints from mineral physics and surface observations. Geophys. J. Int. 167, 1461–1481. Steinberger, B., Torsvik, 2007. Relation between flow in the Earth’s mantle, phase boundary topography and dynamic topography at the Earth’s surface. EGU General Assembly, Vienna, 15–20 April. Geophys. Res. Abst. 9, EGU2007-A-03280. Thoraval, C., Richards, M.A., 1997. The geoid constraint in global geo dynamics: viscosity structure, mantle heterogeneity models and boundary conditions. Geophys. J. Int. 131, 1–8. Vinnik, L., Chevrot, S., Montagner, J.-P., 1997. Evidence for a stagnant plume in the transition zone? Geophys. Res. Lett. 24, 1007–1010. Wen, L., Anderson, D.L., 1997. Layered mantle convection: a model for geoid and topography. Earth Planet. Sci. Lett. 146, 131–143. ˇ O., Fleitout, L., 2003. Effect of lateral viscosity variations in the top 300 km on Cadek, the geoid and dynamic topography. Geophys. J. Int. 152 (3), 566–580. Panasyuk, S.V., Hager, B.H., 2000a. Inversion for mantle viscosity profiles constrained by dynamic topography and the geoid and their estimated errors,. Geophys. J. Int. 143, 821–836.