Mantle dynamics with - Semantic Scholar

Physics of the Earth and Planetary Interiors 217 (2013) 48–58

Contents lists available at SciVerse ScienceDirect

Physics of the Earth and Planetary Interiors journal homepage: www.elsevier.com/locate/pepi

Mantle dynamics with pressure- and temperature-dependent thermal expansivity and conductivity Nicola Tosi a,b,⇑, David A. Yuen c,d,g, Nico de Koker e, Renata M. Wentzcovitch f a

Department of Planetary Geodesy, Technische Universität Berlin, Berlin, Germany Department of Planetary Physics, German Aerospace Center (DLR), Berlin, Germany c Department of Earth Sciences, University of Minnesota, Minneapolis, USA d Minnesota Supercomputing Institute, University of Minnesota, Minneapolis, USA e School of Geosciences, University of the Witwatersrand, Wits, South Africa f Department of Chemical Engineering and Materials Science, University of Minnesota, Minneapolis, USA g School of Environmental Sciences, China University of Geosciences, Wuhan, China b

a r t i c l e

i n f o

Article history: Received 30 November 2012 Received in revised form 11 February 2013 Accepted 22 February 2013 Available online 6 March 2013 Edited by K. Hirose Keywords: Thermal expansivity Lattice thermal conductivity Mantle convection

a b s t r a c t In numerical simulations of mantle convection it is commonly assumed that the coefficients of thermal expansion a and thermal conduction k are either constant or pressure-dependent. Pressure changes are generally computed using parametrizations that rely on extrapolations of low-pressure data for a single upper-mantle phase. Here we collect data for both the pressure and temperature dependence of a from a database of first-principles calculations, and of k from recent experimental studies. We use these datasets to construct analytical parametrizations of a and k for the major upper- and lower-mantle phases that can be easily incorporated into exisiting convection codes. We then analyze the impact of such parametrizations on Earth’s mantle dynamics by employing two-dimensional numerical models of thermal convection. When a is the only variable parameter, both its temperature and pressure dependence enhance hot plumes and tend to inhibit the descent of cold downwellings. Taking into account a variable k leads to a strong increase of the bulk mantle temperature, which reduces the buoyancy available to amplify bottom boundary layer instabilities and causes mantle flow to be driven primarily by the instability of cold plates whose surface velocity also tends to rise. When both parameters are considered together, we observe an increased propensity to local layering which favors slab stagnation in the transition zone and subsequent thickening in the lower mantle. Furthermore, the values of k near the coremantle boundary ultimately control the effect of this physical property on convection, which stresses the importance of determining the thermal conductivity of the post-perovskite phase. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction The coefficients of thermal expansion a and thermal conduction k control the way heat is transported in the interior of solid planets by advection and diffusion, thereby affecting the dynamics and thermal evolution of the mantle. From laboratory experiments, it is well known that both parameters exhibit significant variations with pressure and temperature whose magnitude varies among different mineral phases (e.g. Fei, 1995; Xu et al., 2004; Katsura et al., 2009a,b; Manthilake et al., 2011a). Nevertheless, in most numerical simulations of mantle convection in the Earth and terrestrial planets, a and k are either assumed constant or their pressure (i.e. depth) dependence only is taken into account according to simplified parametrizations (e.g. Leitch et al., 1991; Hansen ⇑ Corresponding author at: Department of Planetary Geodesy, Technische Universität Berlin, Berlin, Germany. E-mail address: [email protected] (N. Tosi). 0031-9201/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.pepi.2013.02.004

et al., 1993; Naliboff and L.H., 2006; Zhong, 2006; Foley and Becker, 2009; Schuberth et al., 2009; Nakagawa et al., 2010; Phillips and Coltice, 2010; Tan et al., 2011). A few attempts have been made to describe in a more detailed way the dynamical impact of variable a and k on mantle dynamics. Schmeling et al. (2003) and Ghias and Jarvis (2008) accounted for the pressure and temperature dependence of the thermal expansivity using a Birch–Murnaghan equation of state and the experimental data of Fei (1995) for olivine only in simple models of isoviscous convection. Matyska and Yuen (2005) discussed the effects of the strongly temperaturedependent radiative part of the thermal conductivity on the dynamics of the lower mantle. Tosi et al. (2010) presented global convection models with multiple phase transitions that incorporated experimental constraints on the thermal and transport properties of lower-mantle minerals. Hunt et al. (2012), in conjunction with measurements of the thermal conductivity of a post-perovskite low-pressure analog, showed simple simulations of isoviscous convection including a phase-dependent jump in thermal

N. Tosi et al. / Physics of the Earth and Planetary Interiors 217 (2013) 48–58

conductivity across the perovskite to post-perovskite interface. Maierová et al. (2012) reported results from kinematic models of subduction incorporating variable thermal conductivity associated with a detailed upper mantle phase-diagram. Important efforts have also been made to account for a self-consistent treatment of mantle mineralogy and thermodynamics in numerical models of mantle convection on the base of pre-compiled look-up tables from which physical properties derived from thermodynamic databases can be interpolated (e.g. Nakagawa et al., 2010; Jacobs and van den Berg, 2011). Nevertheless, simple and up-to-date parametrizations for the major mantle minerals of key parameters such as a and k would facilitate their incorporation into existing numerical convection codes. Here we make use of data of thermal expansivity derived from the thermodynamic database of first-principles simulations of the Virtual Laboratory for Earth and Planetary Materials of the University of Minnesota (VLab, www.vlab.msi.umn.edu). In addition, we collect recent thermal conductivity data based on experiments performed on both the upper and lower-mantle major mineral phases (Xu et al., 2004; Manthilake et al., 2011a). We first employ these data sets to construct analytical parametrizations that describe both the pressure and temperature dependence of a and k for forsterite and its high-pressure phases, as well as for perovskite and periclase. After describing the setup of our numerical model, we discuss the impact of the above parametrizations on the dynamics of the Earth’s mantle using 2D numerical convection simulations in Cartesian geometry featuring a rheology which can mimic plate-like behavior and three major phase transitions covering the entire mantle. 2. Parametrization of the thermal expansivity and conductivity 2.1. Thermal expansivity based on VLab data The VLab service-oriented cyberinfrastructure (da Silva et al., 2007) consists of a portal, web services (da Silveira et al., 2008), workflows for first-principles computations of elastic and thermodynamic properties of materials (Wentzcovitch et al., 2010; da Silveira et al., 2012) in distributed environments (da Silveira et al., 2008, 2011), databases of raw inputs/outputs of successful calculations, and internet applications for obtaining thermodynamic properties of minerals (da Silveira et al., 2008, 2011). The latter are the source of the thermodynamic data utilized in this paper and are freely accessible online (VLab-ThoM, www.vlab.msi.umn.edu/resources/thermodynamics). These applications run interactive quasiharmonic (QHA) calculations (Wallace, 1972) using as input pressure-dependent vibrational density of states (VDoS) generated by first-principles calculations for a collection of minerals (Wentzcovitch et al., 2010). All VDoSs were obtained using codes of the Quantum ESPRESSO distribution (Giannozzi et al., 2009). The Thermodynamics of Minerals (ThoM) application produces volume, thermal expansion coefficient, specific heat (at constant pressure or volume), bulk modulus (isothermal or adiabatic), enthalpy, Gibbs free-energy and Grüneisen parameter, all on user-defined pressure–temperature grids, viewed as tables onscreen, to be downloaded, or analyzed visually as 2D plots through an interactive graphical interface (da Silveira et al., 2011). Among these thermodynamic properties, the thermal expansivity exerts a first-order impact on the dynamics of the mantle as it controls the thermal buoyancy that drives convection. QHA expansivities for Mg2SiO4-forsterite (Li et al., 2007; Yu et al., 2008), wadsleyite (Wu and Wentzcovitch, 2007), -ringwoodite (Yu and Wentzcovitch, 2006), MgO-periclase (Karki et al., 2000a; Umemoto et al., 2006) and MgSiO3-perovskite (Karki et al., 2000b; Umemoto et al., 2006) were then regenerated using the ThoM application. As

49

the thermodynamic properties of perovskite and post-perovskite phases of MgSiO3 are very similar at CMB conditions (Tsuchiya et al., 2005), we have used the same expansivity for both. The expansivity values employed in the simulations that we carried out were obtained within the limit of validity of the QHA (Carrier et al., 2007) and for pure Mg end-member compositions. Regarding the latter assumption, it should be mentioned that the presence of Fe impurities is not expected to have a strong effect on the thermal expansivity (Wu et al., 2009; Metsue and Tsuchiya, 2012). Firstprinciples simulations conducted on ferropericlase ((Mg1xFex)O, x = 0.1875, Wu et al., 2009) agree extremely well with experimental measurements, at zero-pressure, of the expansivity of both periclase (x ¼ 0, Touloukian et al., 1977) and ferropericlase (x ¼ 0:36, Westrenen et al., 2005). Nevertheless, spin transitions would affect the thermal expansion and other properties of the aggregate with possible consequences for mantle convection (Bower et al., 2009; Shahnas et al., 2011). According to Metsue and Tsuchiya (2012), the spin state would also influence the expansivity of ferrosilicate perovskite ((Mg1xFex)SiO3, x = 0.0625). The effect however would not be large and would decrease upon compression. Panels a–d in Fig. 1 show the corresponding distributions of a as a function of pressure and temperature. To obtain Fig. 1d, we assumed a standard lower-mantle assemblage consisting of 80% perovskite and 20% periclase. Fig. 1e depicts the radial profile of the thermal expansivity along a reference adiabat with potential temperature of 1600 K and two 100 km thick thermal boundary layers (TBL) (black lines); the blue and red lines show the profile of a along 200 K colder and hotter adiabats, respectively, while the green line is associated with a linear increase of temperature throughout the mantle without boundary layers. The surface and CMB temperatures are set to 300 K and 4000 K, respectively. The variability of the thermal expansivity is mostly evident in the uppermost part of the mantle. Here the strong increase of temperature across the top TBL determines a significant growth of a, which reaches almost twice the reference surface value for the hottest geotherm. Below the TBL, a steady decrease throughout the upper mantle is observed, with a pronounced jump across the olivine–wadsleyite transition where a decreases abruptly. A second, albeit less pronounced, jump is observed at the wadsleyite–ringwoodite transition. At the top of the lower mantle, another discontinuity is evident, with perovskite and periclase having a slightly higher expansivity than ringwoodite. Throughout the lower mantle, the value of a steadily diminishes apart from a minor increase due to temperature across the bottom TBL. Apart from a small increase in the uppermost mantle, the expansivity profile resulting from a linearly increasing temperature (green line) resembles the distribution used in previous studies in which only a depth-dependent and steadily decreasing a was taken into account (e.g. Tosi et al., 2010). Although this assumption appears to be reasonable for the lower mantle, it leads to an underestimation of a by up to  40% in the upper mantle. In order to facilitate the incorporation of the pressure and temperature dependence of the thermal expansivity in mantle convection codes, we used the Levenberg–Marquardt algorithm (e.g. Moré, 1978) to fit the VLab data to the following four-parameters function:

aðT; PÞ ¼ ða0 þ a1 T þ a2 T 2 Þ expða3 PÞ;

ð1Þ

where T is absolute temperature in K, P is hydrostatic pressure in GPa, and ai ; i ¼ 0; . . . ; 3 are phase-dependent coefficients obtained from the inversion of the data. Since the governing equations of thermal convection are generally discretized in terms of depth instead of hydrostatic pressure, we also provide a formula for the temperature and depth dependence of a that we obtained by integrating numerically the PREM density profile together with the associated gravity (Dziewonski and Anderson, 1981):

50

N. Tosi et al. / Physics of the Earth and Planetary Interiors 217 (2013) 48–58

Fig. 1. Pressure and temperature distribution of the thermal expansivity of upper (a–c) and lower mantle phases (d) based on VLab data. Panel e shows the corresponding profiles along a reference (black), hot (red) and cold (blue) adiabat supplemented by top and bottom thermal boundary layers, and along a linear temperature profile (green). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

aðT; zÞ ¼ ðb0 þ b1 T þ b2 T 2 Þ expðb3 zÞ;

ð2Þ

where z indicates the depth in m. Numerical values of the coefficients ai and bi are reported in Table 1. The parametrizations (1) and (2) have been obtained using VLab data from the following

pressure and temperature ranges: P 2 ½0; 25 GPa and T 2 ½300; 2000 K for the upper mantle and P 2 ½25; 140 GPa and T 2 ½300; 4000 for the lower mantle. The functional dependence of temperature employed in Eqs. (1) and (2) is somewhat standard and has already been shown to be

Table 1 Phase-dependent coefficients for the parametrization of aðT; PÞ (ai , Eq. (1)), aðT; zÞ (bi , Eq. (2)), kðT; PÞ (ci , Eq. (3)) and kðT; zÞ (di , Eq. (4)). Forsterite

a0 a1 a2 a3 b0 b1 b2 b3 c0 c1 c2 d0 d1 d2

[K1] [K2] [K] [GPa1] [K1] [K2] [K] [m1] [Wm1 K1] [Wm1 K1 GPa1] [Wm1 K1] [Wm2 K1]

3.15  105 1.02  108 0.76 3.63  102 3.15  105 1.02  108 0.76 1.27  106 2.47 0.33 0.48 2.47 1.15  105 0.48

Wadsleyite

2.84  105 6.49  109 0.88 2.61  102 2.84  105 6.49  109 0.88 9.16  107 3.81 0.34 0.56 3.81 1.19  105 0.56

Ringwoodite

2.68  105 4.82  109 0.93 2.15  102 2.68  105 4.82  109 0.93 7.55  107 3.52 0.36 0.61 3.52 1.26  105 0.61

80% Perovskite

80% Perovskite

20% Periclase (A)

20% Periclase (B)

2.68  105 2.77  109 1.21 8.63  103 2.68  105 2.77  109 1.21 3.76  107 13.6 0.42 0.58 13.6 1.88  105 0.58

3.48 0.12 0.31 3.48 5.17  106 0.31

N. Tosi et al. / Physics of the Earth and Planetary Interiors 217 (2013) 48–58

appropriate in fitting experimental measurements of thermal expansivity (e.g. Fei, 1995; Katsura et al., 2009b; Komabayashi et al., 2008). The use of an exponential dependence for pressure or depth, instead, is different from traditional fits, which employ the density q and the Anderson-Grüneisen parameter dT to describe changes of a upon compression as proportional to ðq0 =qÞdT , where q0 is density at reference pressure. A parametrization that explicitly uses pressure or depth has the advantage over the latter of being readily usable in convection codes, which often assume the Boussinesq approximation and hence constant density. The functional forms that we introduced to parametrize both the thermal expansivity and thermal conductivity (see Eqs. (3) and (4)) reproduce the data within an uncertainty of at most 10% throughout the PT-range considered, which is not expected to have a significant impact on the outcomes of numerical simulations. 2.2. Thermal conductivity based on experimental measurements We compile a model for the lattice thermal conductivity of the entire mantle from experimental measurements of Xu et al. (2004) and Manthilake et al. (2011a), using a modified Debye model to represent the pressure and temperature dependence of k and extrapolate measured values to the conditions present in the lowermost mantle (de Koker, 2010; Manthilake et al., 2011b). In all these experiments k was measured using the Angstrom method in the 500 ton multi-anvil press at the Bayerisches Geoinstitut in Bayreuth (Germany). Xu et al. (2004) measured the lattice thermal conductivity of (Mg0.9Fe0.1)2SiO4 at pressure and temperature conditions of the upper mantle. Manthilake et al. (2011a) reported measurements of the lattice thermal conductivity of both end-member MgSiO3 perovskite and MgO periclase, and of the same two minerals but considering also the presence of Fe. Fig. 2a–c show the PT-dependence of k for upper mantle phases based on the data of Xu et al. (2004). Fig. 2d and e refer to the lower mantle and are based on the data of Manthilake et al. (2011a) assuming 80% perovskite and 20% periclase for pure Mg end-members (hereafter model A) and considering 3 mol% FeSiO3 in perovskite and 20 mol% FeO in periclase (hereafter model B). Fig. 2f illustrates the corresponding profiles along the same temperature distributions used to plot Fig. 1e. The thermal conductivity increases nearly linearly across the upper mantle, apart from small changes in the transition zone, whose complexity however, is not captured by our models and would require more detailed parametrizations (Maierová et al., 2012). As far as the lower mantle is concerned, as discussed extensively by Manthilake et al. (2011a), the incorporation of iron in silicate perovskite and ferropericlase (model B) results in a  50% decrease of the lattice thermal conductivity relative to the end-member compositions (model A). In both cases, as shown in Fig. 2f, k increases steadily below the 660 km discontinuity and reaches values at the core-mantle boundary of  15 W/mK (model A, solid lines) and  8 W/mK (model B, dashed lines). As in what discussed before, the profile resulting from a linearly increasing temperature (green line) resembles the distribution of conductivity used in studies which feature only a depth-dependent and steadily increasing k (e.g. Tosi et al., 2010). While this profile reasonably approximates the conductivity distribution in the lower mantle, it overestimates k by up to  40% in the upper mantle. We parameterized the conductivity data using a three-parameters function of temperature and pressure:

kðT; PÞ ¼ ðc0 þ c1 PÞ

 c 300 2 ; T

and a function of temperature and depth:

ð3Þ

kðT; zÞ ¼ ðd0 þ d1 zÞ

 d 300 2 : T

51

ð4Þ

Numerical values of the coefficients ci and di are reported in Table 1. It should be noted that the above relations were chosen in order to provide an explicit dependence on pressure and depth as these are the variables typically used in convection simulations. Eqs. (3) and (4) differ from those employed to extrapolate the experimental measurements with the Debye model for which a dependence of the lattice conductivity on temperature and volume must be assumed. As in the case of thermal expansivity, they have been obtained using data from the following pressure and temperature ranges: P 2 ½0; 25 GPa and T 2 ½300; 2000 K for the upper mantle and P 2 ½25; 140 GPa and T 2 ½300; 4000 K for the lower mantle. It should be also noted that, differently from the fit to the thermal expansivity data, the above ranges contain extrapolated values. Finally, according to recent experimental evidence, the thermal conductivity of post-perovskite, although probably highly anisotropic, is likely to be up to twice as high as that of perovskite (Cheng et al., 2011; Hunt et al., 2012; Ohta et al., 2012). We explore this behavior by performing a simulation where the conductivity of post-perovskite is represented as double the value for perovskite at the same conditions. 3. Mantle convection model In order to assess the effects of our parametrizations of thermal expansivity and conductivity, we performed numerical simulations of thermal convection in a 2D plane-layer geometry using our finite-volume code YACC (Yet Another Convection Code, code.google.com/p/yacc-convection/, King et al. (2010), Tosi et al. (2010)). To this end, we solved the conservation equations of mass, momentum and energy under the extended Boussinesq approximation (EBA) (Christensen et al., 1985) in a Cartesian box with an aspect-ratio of three with free-slip boundary conditions, isothermal top and bottom boundaries and reflective side-walls. The relatively small aspect-ratio of the box is not expected to have a significant impact on the results as simulations conducted with an aspect-ratio of ten exhibit comparable trends (Tosi et al., 2010). We used a regular computational grid consisting of 300 and 200 points in the horizontal and vertical directions, respectively. As initial condition, we used for all simulations an adiabatic temperature profile with potential temperature of 1600 K supplemented by 100 km thick top and bottom boundary layers and by a small (1%) random perturbation to induce the onset of convection. Although under the EBA, thermal effects due to adiabatic compression and viscous dissipation are taken into account, the compressible anelastic-liquid approximation (ALA) (Jarvis and McKenzie, 1980) would have been generally more appropriate. However, it must be noted that when using the ALA, the presence of a pressure- and temperature-dependent thermal expansivity requires a careful choice of the reference density profile that should be updated in a self-consistent manner during the simulation. The use of the EBA, and hence of constant density everywhere except in the buoyancy term of the momentum equation, allows us to circumvent this additional complexity. As within the temperature and pressure ranges appropriate for the Earth’s mantle, the difference between compressible and EBA formulations is expected to be small, we decided thus to assume our fluid to be incompressible and to analyze the effects of time-variable compressibility in a separate work. The setup of the numerical model and choice of material parameters are essentially the same as those adopted in Tosi et al. (2010) and Tosi and Yuen (2011). All models are isochemical and feature a mixed heating mode, with the ratio between internal heating and

52

N. Tosi et al. / Physics of the Earth and Planetary Interiors 217 (2013) 48–58

Fig. 2. Pressure and temperature distribution of the thermal conductivity of upper mantle phases based on the data of Xu et al. (2004) (a–c) and of lower mantle phases based on the data of Manthilake et al. (2011a) (d, e). Panel f shows the corresponding profiles along a reference (black), hot (red) and cold (blue) adiabat supplemented by top and bottom thermal boundary layers, and along a linear temperature profile (green). Solid and dashed lines refer to the lower mantle conductivity models A and B, respectively (see Section 2.2 for details). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

basal heating Rayleigh numbers set to one in order to avoid an unrealistic increase of the mantle temperature due to the choice of the Cartesian geometry (O’Farrell and Lowman, 2010). In all simulations, we employed a reference thermal Rayleigh number of 2  107 and a dissipation number of 0.6 based on the ambient value of the thermal expansivity of forsterite at 300 K (i.e. 2:6  105 K1). The main difference with respect to our previous studies is in the rheological law. Here, instead of treating the mantle as a linear fluid with a pressure- and temperature-dependent viscosity with a low activation energy to allow for the mobility of the upper surface, we employed a Bayerlee-type plastic yielding. The adopted non-dimensional viscosity g is then a function of depth z, temperature T and of the second invariant of the strain-rate tensor e_ II and takes the following form:



gðz; T; e_ II Þ ¼ min exp ½ lnðcT ÞðT  T 0 Þ þ lnðcz Þz;



rY ðzÞ ; 2e_ II

ð5Þ

where cT ¼ 108 and cz ¼ 102 are the viscosity contrasts due to temperature and depth respectively, T 0 ¼ 1600 K is the temperature at which a reference viscosity of 1022 Pas is set and rY ðzÞ is a depthdependent yield-stress. For the latter we used a surface value of

117 MPa combined with an increase with depth of 162.3 Pa/m (Nakagawa and Tackley, 2005). In addition, we considered a 30-fold phase-dependent viscosity contrast across the transition from ringwoodite to the lower mantle. We took into account four phase transitions: from forsterite to wadsleyite at 400 km depth (Clapeyron slope k ¼ 3 MPa/K, density jump dq ¼ 273 kg/m3), from wadsleyite to ringwoodite at 520 km, from ringwoodite to perovskite and periclase at 670 km (k ¼ 2:5 MPa/K, dq ¼ 342 kg/m3) and from perovskite to postperovskite according to local pressure and temperature conditions (Tosi et al., 2010) (k ¼ 13 MPa/K, dq ¼ 67:5 kg/m3). Because of the complexity of the transition from wadsleyite to ringwoodite and other uncertainties in the density jump (e.g. Yu et al., 2008), we neglected buoyancy effects and latent heat associated with this phase change and took into account only the accompanying, albeit small, variations of a and k across a flat boundary at 520 km depth. Although the rheology of post-perovskite is likely to be weaker than that of perovskite by one order of magnitude or more (Hunt et al., 2009; Ammann et al., 2010; Dobson et al., 2012) and its effects on mantle convection can be significant (Cˇizˇková et al., 2009; Tosi et al., 2010; Samuel and Tosi, 2012), for simplicity, and because uncertainties still exist (Karato, 2010), we assumed that this phase has the same viscosity of the surrounding perovskite.

N. Tosi et al. / Physics of the Earth and Planetary Interiors 217 (2013) 48–58

4. Modeling results In Figs. 3 and 4 we show representative snapshots of selected simulations that were run until they reached statistical steadystate and featuring constant expansivity a and conductivity k (Fig. 3a), variable a and constant k (Figs. 3b, e and 4a), constant a and variable k (Fig. 3c, f and d) and variable a and k (Fig. 3d, g, b, c, e and f). Fig. 3 shows temperature deviations with respect to the laterally averaged profile, i.e. Tðx; zÞ  hTðzÞi, while Fig. 4 shows the non-dimensional distributions of thermal expansivity and/or conductivity for models in which the temperature dependence of these parameters was taken into account. Additionally, in Fig. 5 we present profiles of the temperature (panel a), thermal buoyancy (panel b) and conductivity (panel c) averaged over a non-dimensional time-interval of 0.01, which, scaling the time with our surface value of thermal diffusivity, corresponds to approximately 5 Ga. The thermal buoyancy is defined as aðx; zÞjTðx; zÞ  hTðzÞij. From Fig. 3 it is evident that the influence of variable a and k on mantle convection is dramatic. When both parameters are held

53

constant (Fig. 3a and black lines in Fig. 5), the temperature field reflects that characteristic of convection at relatively high Rayleigh numbers with internal heating. The flow is dominated by cold downwellings, which, despite the viscosity jump between upper and lower mantle, sink straight to the CMB and accumulate there, cooling the deepest part of the mantle and causing the occurrence of locally elevated post-perovskite lenses. Because of the limited temperature jump across the bottom TBL, hot upwellings are generally weak and short-lived. As well known from previous studies (e.g. Leitch et al., 1991; Hansen et al., 1993; Tosi et al., 2010), including variable thermal expansivity (Figs. 3b, e and 4a blue lines in Fig. 5) causes a major change in the convection style: hot, bottom boundary layer instabilities tend to grow during their upward rise, resulting in largescale plumes. Cold downwellings lose buoyancy while sinking, a process which cools the mantle strongly. With respect to the reference run with constant thermodynamic properties, in fact, we observed a temperature drop as high as 0.27 (1000 K) above the bottom TBL. While in the bulk of the mantle the distribution of a

Fig. 3. Non-dimensional temperature anomalies Tðx; zÞ  hTðzÞi from simulations with constant a and k (a), variable a and constant k (b, e), constant a and variable k (c, f), variable a and k (d, g). Arrows indicate flow velocity directions. Thick black lines denote phase boundaries (the boundary between wadsleyite and ringwoodite lies within the transition zone and is not plotted).

54

N. Tosi et al. / Physics of the Earth and Planetary Interiors 217 (2013) 48–58

Fig. 4. Non-dimensional distributions of thermal expansivity (a, b, c) and conductivity (d, e, f) for the simulations shown in Fig. 3 that feature also the temperature dependence of either parameters. Panel a corresponds to Fig. 3e, panel d to Fig. 3f, panels b and e to Fig. 3d, panels c and f to Fig. 3g.

Fig. 5. Time-averaged profiles of non-dimensional temperature (a), thermal buoyancy (b) and conductivity (c) for the simulations shown in Figs. 3 and 4 plus one additional simulation assuming a highly conductive post-perovskite phase (see text for details).

is mostly influenced by pressure, near the boundary layers do the thermal effects become significant (Fig. 4a). The highest values of a are associated with plume heads near the surface, while the lowest are attained near the CMB and coincide with post-perovskite regions which are now much more prominent because of the reduced mantle temperature. The above features are quickly recog-

nized by looking at the blue lines in Fig. 5b showing the profile of thermal buoyancy. Although the convection patterns displayed in Fig. 3b (obtained with depth-dependent a according to the green line in Fig. 1e) and Fig. 3e (obtained with depth- and temperaturedependent a) are similar, a comparison of solid and dashed blue lines of Fig. 5, clearly shows that the temperature dependence of

N. Tosi et al. / Physics of the Earth and Planetary Interiors 217 (2013) 48–58

a is important. A PT-dependent expansivity in fact delivers average mantle temperatures up to 0.06 (220 K) lower than those obtained from the model with P-dependent expansivity only, the largest differences being observed in the transition zone and in the upper half of the lower mantle. In the presence of variable thermal conductivity (Fig. 3c, f, d and red lines in Fig. 5), the mantle temperature increases sufficiently for plumes to be suppressed and post-perovskite to only form where downwellings reach the deep mantle. The distribution of k is dominated by pressure effects apart in the core of cold slabs where it attains its highest values when temperature dependence is included (Fig. 4d). Note however, that despite their high conductivity, slabs do not heat up efficiently via conduction since they sink very rapidly because of the generally low viscosity induced by the surrounding hot mantle. The latter affects in turn also surface plate velocities. On the one hand, using constant properties or variable expansivity, they are in the range of a few cm/a or a fraction of cm/a, respectively. On the other hand, with variable conductivity they can become up to one order of magnitude larger (for example, the peak surface velocity associated with the slab of Fig. 3f is  30 cm/a). The differences between models featuring P-dependent (green lines in Fig. 2f) and PT-dependent k are minor, both qualitatively (compare Fig. 3c and d) and quantitatively (compare solid and dashed red lines in Fig. 5). In terms of temperature distribution, the largest differences are attained in the lower mantle and amount  0:02 (74 K). When both variable a and k are considered together (Figs. 3d, g, 4b, c, e, f and green lines in Fig. 5), yet another convective regime is observed. Because of the temperature increase caused by the elevated thermal conductivity at depth and to the modified buoyancy properties induced by variations of thermal expansivity, up- and downwellings acquire similar importance. In particular, whether we use either A or B models, convection shows an increased propensity towards local layering, with slabs that are often trapped in the transition zone before sinking into the lower mantle where the associated cold thermal anomaly also tends to broaden. Occasionally, slab detachment arising as a consequence of prolonged stagnation is also observed as shown in Fig. 3g. Furthermore, the slowing effect exerted by variable a helps to bring the high surface velocities observed with variable k only to more realistic values of  4  5 cm/a. The profiles of Fig. 5 also include the results of one additional simulation based on conductivity model B but including also a highly conductive post-perovskite phase (dashed green lines). In general, using a milder increase in the thermal conductivity of the lower mantle (i.e. model B instead of model A) leads to a slight reduction of the temperature profile above the CMB (compare solid and dotted green lines). However, when a jump in the conductivity of post-perovskite is considered, both the temperature and buoyancy profiles are hardly distinguishable from those obtained when using an overall higher conductivity throughout the lower mantle (compare solid and dashed green lines).

5. Discussion Still today, numerical simulations of thermal and thermo-chemical convection in the mantle of the Earth and terrestrial planets generally assume that the coefficients of thermal expansion and conduction are constant or at most vary monotonically with depth according to simplified parametrizations that account neither for the temperature nor for the phase dependence of these two important parameters. The possibility offered by the VLab to easily access an online database of results of first-principles simulations prompted us to collect thermal expansivity data for the main min-

55

eral phases of the Earth’s upper and lower mantle over a wide range of pressures and temperatures. For such phases, we also compiled data based on two experimental studies focusing on measurements of the lattice thermal conductivity (Xu et al., 2004; Manthilake et al., 2011a). We made use of the above information to derive two families of parametrizations that describe the pressure (or depth), temperature and phase dependence of a and k. The analytical formulae that we provided (Eqs. (1)–(4)) are simple and can be easily implemented in any mantle convection code. We studied the impact of our parametrizations on the longterm dynamics of the mantle using 2D Cartesian models of purely thermal convection with a visco-plastic rheology in the framework of the extended Boussinesq approximation. Variable a and k influence the dynamics of the mantle in a very significant way. The use of a PT-dependent thermal expansivity lengthens the spatial scale of mantle convection, favoring the importance of plumes over that of downwellings. The former, in fact, gain buoyancy while rising through the mantle due to an increase of a with decreasing pressure. As plumes approach the top boundary layer, the increase of expansivity with temperature becomes relevant and they acquire additional buoyancy. The hottest parts of plume heads, are characterized by expansivity values up to 40% higher than those of the ambient mantle. The opposite argument applies to cold slabs. These are generally associated with low expansivity that obstructs their subduction and causes them to spread slowly, inducing in turn an overall cooling of the lower mantle. It is important to point out that the temperature dependence of a plays a significant role in amplifying these effects. Models in which only its depth dependence is considered deliver average temperatures up to 220 K higher and are characterized by 30–40% reduction of the thermal buoyancy throughout the mantle, but especially in its upper part where thermal effects are more prominent. While the effects of variable thermal conductivity in the upper mantle are minor, the increase of k in the lower mantle determines an overall heating of the system, which can be as high as 500 K when the temperature profiles obtained using variable k are compared against the reference state with constant a and k. This has two major consequences: buoyancy driven plumes tend to disappear and post-perovskite only forms where cold plates reach the lowermost mantle. In addition, cold slabs reaching the CMB heat rapidly due to the inverse relationship between k and T, so that post-perovskite lenses are relatively short-lived. Nevertheless, models with temperature- and pressure-dependent k do not exhibit important changes with respect to their counterparts with pressure-dependent k only. Therefore, as far as the lattice thermal conductivity is concerned, taking into account only its pressure dependence appears to be a good approximation. However, this does not rule out the possibility that the temperature-dependent radiative part of the thermal conductivity (e.g. Goncharov et al., 2008), which is still debated and was here neglected, plays an important role in dictating the dynamics of the lowermost mantle (Matyska and Yuen, 2005; Yuen et al., 2007). We tested two different conductivity models for the lower mantle following Manthilake et al. (2011a). Model A, in which pure Mg end-member perovskite and periclase compositions are assumed, delivers slightly larger temperatures (up to  100 K) in the bottom half of the lower mantle than model B, in which the presence of Fe is also accounted for. Notably, the simulation run in which we used model B in conjunction with a highly conductive post-perovskite is essentially indistinguishable from model A with no conductivity contrast across the perovskite to post-perovskite interface. This demonstrates that the conductivity values above the CMB are those that eventually determine the net effects of this parameter on mantle convection as also recognized by Yanagawa et al. (2005). In this view, recent efforts devoted to estimate the lattice

56

N. Tosi et al. / Physics of the Earth and Planetary Interiors 217 (2013) 48–58

Fig. 6. Selected snapshots of the temperature field obtained from a simulation with temperature- and pressure-dependent thermal expansivity and conductivity (model B) showing different subduction regimes: straight subduction (a), slab detachment (b), slab folding (c) and slab stagnation (d). In all panels only one third of the horizontal domain is shown.

conductivity of post-perovskite are thus particularly important (Cheng et al., 2011; Hunt et al., 2012; Ohta et al., 2012). The dynamical consequences observed when either of the two parameters was held constant while the other varied with pressure and temperature can be also recognized when both a and k are variable. In such a case, however, an intermediate state is found where the increase of temperature induced by the conductivity is reduced in the lowermost mantle by up to  200 K with respect to the case in which k varies with depth and temperature or depth only and a is constant. This, in turn, makes possible the emergence of buoyant upwellings because of the increased temperature jump across the bottom TBL. In addition, irrespective of the conductivity parametrization used (i.e. models A and B, used in Fig. 3d and g, respectively), we observed a tendency to local layering, which was not found in all other simulations we carried out. As a consequence, slabs may stagnate at the bottom of the transition zone and occasionally detach from their upper mantle part. Slab thickening in the lower mantle caused by the depth-decreasing thermal expansivity and occasional folding is another important feature that we found. In Fig. 6, we show selected snapshots of the temperature field obtained using the conductivity model B, which clearly illustrate the variety of subduction regimes that we encountered, namely straight subduction into the lower mantle (a), slab detachment (b), slab folding below the transition zone (c) and slab stagnation (d). This behavior is not characteristic of specific time-spans. It keeps occurring during the course of this simulation and is not observed when either a or k are held constant. Furthermore all the above features are clearly recognized in tomographic models (see e.g. Fukao et al., 2009 for slab stagnation, van der Meer et al. (2010) for detachment and identification of lower mantle slab remnants and Ribe et al. (2007) for slab thickening and buckling in the lower mantle). We can thus argue that the set of models considering both variable a and k is the one that best simulates an actual mantle scenario. It is important to point out that the various convection styles that we discussed are also the result of the assumed rheology, which for simplicity was kept fixed throughout the study. Although the use of different combinations of rheological parameters can have an impact on the predicted convection planforms, the effects induced by the pressure- and temperature dependence of a and k represent a robust feature of our models which can be hardly suppressed by the choice of a different rheology. For example, the use of a visco-plastic rheology instead of a low activation energy to simulate a mobile-lid regime is crucial for observing the variety of subduction styles discussed above (Fig. 6) and and for obtaining a plate-like behavior at the surface and at depth in the mantle. However, it does not affect our conclusions regarding the relative changes of important quantities such as mean temperature or

buoyancy profiles upon variations of thermal expansivity and/or conductivity. 6. Conclusions Despite the inherent simplifications of the numerical model, we expect our findings to be applicable in several geodynamic settings, both at a global and regional scale. For example, the temperature dependence of a could be important for the development of small-scale convection in the upper mantle (e.g. Ballmer et al., 2011), for the generation of plumes which are not directly associated with a deep source (e.g. Conrad et al., 2010), and because of the significant temperature differences it induces in the upper mantle, it could also affect melting processes (e.g. de Smet et al., 1999), the conditions for the onset of plate tectonics (e.g. O’Neill et al., 2007a,b) and the way mantle plumes interact with phase transitions (e.g. Yoshida and Ogawa, 2004). Moreover, the style of subduction, which is often studied assuming constant thermodynamic parameters (e.g. Beˇhounková and Cˇ´ızˇková, 2008; Billen and Hirth, 2007), could be influenced in a significant way, as much as the rheological constraints that can be derived from dynamical subduction models (e.g. Liu and Stegman, 2011; Cˇizˇková et al., 2012). The fact that in our most complex models, which should also be the most realistic among those presented, we observed phenomena like slab stagnation, detachment and thickening suggests that the variability of the two parameters considered here can contribute to explain well established seismic observations with numerical models. Nevertheless, in order to fully assess their robustness, our results will ultimately have to be analyzed in a three-dimensional context. In view of their important effects, we suggest that variable expansivity and conductivity should be routinely included in numerical simulations of mantle convection in the Earth and terrestrial planets. The analytical formulae that we provided should facilitate this task. Acknowledgments We are grateful to Kei Hirose for his editorial handling of the manuscript and to Ikuko Wada and two anonymous reviewers for their constructive comments. We also thank Zhongqing Wu for stimulating discussions about thermal expansivity and Pedro da Silveira for help with the VLab database. N. Tosi acknowledges support form the Deutsche Forschungs Gemeinschaft (Grant No. TO 704/1-1) and the Czech Science Foundation (Project P210/11/ 1366). We also acknowledge support from a NSF grant awarded jointly to R. Wentzcovitch and D.A. Yuen by the Geophysics and Geochemistry programs and to D.A. Yuen from the CMG program.

N. Tosi et al. / Physics of the Earth and Planetary Interiors 217 (2013) 48–58

References Ammann, M.W., Brodholt, J.P., Wookey, J., Dobson, D.P., 2010. First-principles constraints on diffusion in lower-mantle minerals and a weak D00 layer. Nature 465, 462–465. http://dx.doi.org/10.1038/nature09052. Ballmer, M.D., Ito, G., van Hunen, J., Tackley, P.J., 2011. Spatial and temporal variability in hawaiian hotspot volcanism induced by small-scale convection. Nat. Geosci. 4, 457–460. http://dx.doi.org/10.1038/ngeo1187. Billen, M.I., Hirth, G., 2007. Rheologic controls on slab dynamics. Geochem. Geophys. Geosyst. 8, Q08012. http://dx.doi.org/10.1029/2007GC001597. Bower, D.J., Gurnis, M., Jackson, J.M., Sturhahn, W., 2009. Enhanced convection and fast plumes in the lower mantle induced by the spin transition in ferropericlase. Geophys. Res. Lett. 36, L10306. http://dx.doi.org/10.1029/2009GL037706. ˇ´ızˇková, H., 2008. Long-wavelength character of subducted slabs Beˇhounková, M., C in the lower mantle. Earth Planet. Sci. Lett. 275, 43–53. http://dx.doi.org/ 10.1016/j.epsl.2008.07.059. Carrier, P., Wentzcovitch, R.M., Tsuchiya, J., 2007. First principles prediction of crystal structures at high temperatures using the quasiharmonic approximation. Phys. Rev. B. 76, 064116. http://dx.doi.org/10.1103/ PhysRevB.76.064116. Cheng, J.-G., Zhou, J.-S., Goodenough, J.B., Sui, Y., Ren, Y., Suchomel, M.R., 2011. High-pressure synthesis and physical properties of perovskite and postperovskite Ca1xSrxIrO3. Phys. Rev. B 83, 064401. http://dx.doi.org/10.1103/ PhysRevB.83.064401. Christensen, U.R., Yuen, D.A., 1985. Layered convection induced by phase transitions. J. Geophys. Res. 90 (B12), 10,291–10,300. ˇ izˇková, H., C ˇ adek, O., Matyska, C., Yuen, D.A., 2009. Implications of post-perovskite C transport properties for core-mantle dynamics. Phys. Earth Planet. Inter. 180 (3–4), 235–243. http://dx.doi.org/10.1016/j.pepi.2009.08.008. ˇ izˇková, H., van den Berg, A., Spakman, W., Matyska, C., 2012. The viscosity of C Earth’s lower mantle inferred from sinking speed of subducted lithosphere. Phys. Earth Planet. Inter. 200–201, 56–62. http://dx.doi.org/10.1016/ j.pepi.2012.02.010. Conrad, C.P., Wub, B., Smith, E.I., Bianco, T.A., Tibbetts, A., 2010. Shear-driven upwelling induced by lateral viscosity variations and asthenospheric shear: a mechanism for intraplate volcanism. Phys. Earth Planet. Inter. 178, 162–175. http://dx.doi.org/10.1016/j.pepi.2009.10.001. da Silva, C.R.S., da Silveira, P.R.C., Karki, B., Wentzcovitch, R.M., Jensen, P.A., Bollig, E.F., Pierce, M., Erlebacher, G., Yuen, D.A., 2007. Virtual laboratory for planetary materials: system service architecture overview. Phys. Earth Planet. Inter. 163 (321–332), 321–332. http://dx.doi.org/10.1016/j.pepi.2007.04.018. da Silveira, P., da Silva, C.R.S., Wentzcovitch, R.M., 2008. Metadata management for distributed first principles calculations in VLab – a collaborative cyberinfrastructure for materials computation. Comput. Phys. Commun. 178, 186. http://dx.doi.org/10.1016/j.cpc.2007.09.001. da Silveira, P.R.C., nez Valdez, M.N., Wenzcovitch, R.M., Pierce, M., da Silva, C.R.S., Yuen, D.A., 2011. Virtual laboratory for planetary materials (VLab): an updated overview of system service architecture. In: Proceedings of the 2011 TeraGrid Conference: Extreme Digital Discovery. ACM, New York, NY, USA, p. 33, doi:10.1145/2016741.2016777. da Silveira, P.R.C., Holiday, A., nez Valdez, M.N., Wentzcovitch, R.M., 2012. Internet application for first principles elasticity of materials. In: Proc. of 2012e-Science Conference. de Koker, N., 2010. Thermal conductivity of MgO periclase at high pressure: Implications for the D00 region. Earth Planet. Sci. Lett. 292, 392–398. http:// dx.doi.org/10.1016/j.epsl.2010.02.011. de Smet, J.H., van den Berg, A.P., Vlaar, N.J., 1999. The evolution of continental roots in numerical thermo-chemical mantle convection models including differentiation by partial melting. Lithos 48 (48), 153–170. http://dx.doi.org/ 10.1016/S0024-4937(99)00028-6. Dobson, D.P., McCormack, R., Hunt, S.A., Ammann, M.W., Weidner, D., Li, L., Wang, L., 2012. The relative strength of perovskite and post-perovskite NaCoF3. Min. Mag. 76 (4), 925–932. http://dx.doi.org/10.1180/minmag.2012.076.4.09. Dziewonski, A.M., Anderson, D.L., 1981. Preliminary reference Earth model. Phys. Earth Planet. Inter. 25, 297–356. Fei, Y., 1995. Thermal expansion. In: Ahrens, T.J. (Ed.), Mineral Physics and Crystallography – a Handbook of Physical Constants. AGU, Washington, pp. 283–291. Foley, B.J., Becker, T., 2009. Generation of plate-like behavior and mantle heterogeneity from a spherical, viscoplastic convection model. Geochem. Geophys. Geosyst. 10 (8), 4557–4571. http://dx.doi.org/10.1029/2009GC002378. Fukao, Y., Obayashi, M., Nakakuki, T., Utada, H., Suetsugu, D., Irifune, T., Ohtani, E., Yoshioka, S., Shiobara, H., Kanazawa, T., Hirose, K., 2009. Stagnant slab: a review. Annu. Rev. Earth Planet. Sci. 37, 19–46. http://dx.doi.org/10.1146/ annurev.earth.36.031207.124224. Ghias, S.R., Jarvis, G.T., 2008. Mantle convection models with temperature- and depth-dependent thermal expansivity. J. Geophys. Res. 113, B08408. http:// dx.doi.org/10.1029/2007JB005355. Giannozzi, P., Baroni, S., Bonini, N., Calandra, M., Wentzcovitch, R.M., 2009. Quantum ESPRESSO: a modular and open-source software project for quantum simulations of materials. J. Phys.: Condens. Matter 21, 395502. http://dx.doi.org/10.1088/0953-8984/21/39/395502. Goncharov, A., Haugen, B., Struzhkin, V., Beck, P., Jacobsen, S., 2008. Radiative conductivity in the earths lower mantle. Nature 456 (7219), 231–234. http:// dx.doi.org/10.1038/nature07412.

57

Hansen, U., Yuen, D.A., Kroening, S.E., Larsen, T.B., 1993. Dynamical consequences of depth-dependent thermal expansivity and viscosity on mantle circulations and thermal structure. Phys. Earth Planet. Inter. 77, 205–223. Hunt, S.A., Weidner, D.J., Li, L., Wang, L., Walte, N.P., Brodholt, J.P., Dobson, D.P., 2009. Weakening of calcium iridate during its transformation from perovskite to post-perovskite. Nat. Geosci. 2, 794–797. http://dx.doi.org/10.1038/ngeo663. Hunt, S.A., Davies, D.R., Walker, A.M., McCormack, R.J., Wills, A.S., Dobson, D.P., Li, L., 2012. On the increase in thermal diffusivity caused by the perovskite to post perovskite phase transition and its implications for mantle dynamics. Earth Planet. Sci. Lett. 319, 96–103. http://dx.doi.org/10.1016/j.epsl.2011.12.009. Jacobs, M.H.G., van den Berg, A.P., 2011. Complex phase distribution and seismic velocity structure of the transition zone: convection model predictions for a magnesium-endmember olivine–pyroxene mantle. Phys. Earth Planet. Inter. 186 (1–2), 36–48. http://dx.doi.org/10.1016/j.pepi.2011.02.008. Jarvis, G.T., McKenzie, D.P., 1980. Convection in a compressible fluid with infinite Prandtl number. J. Fluid Mech. 96, 515–583. Karato, S., 2010. The influence of anisotropic diffusion on the high-temperature creep of a polycrystalline aggregate. Phys. Earth Planet. Inter. 183, 468–472. http://dx.doi.org/10.1016/j.pepi.2010.09.001. Karki, B.B., Wentzcovitch, R.M., de Gironcoli, S., Baroni, S., 2000a. High-pressure lattice dynamics and thermoelasticity of MgO. Phys. Rev. B 61 (13), 8793–8800. Karki, B.B., Wentzcovitch, R.M., DeGironcoli, S., Baroni, S., 2000b. Ab initio lattice dynamics of MgSiO3 perovskite at high pressure. Phys. Rev. B 62 (22), 14,750– 14,756. Katsura, T., Shatskiy, A., Manthilake, M., Zhai, S., Fukui, H., Yamazaki, D., Matsuzaki, T., Yoneda, A., Ito, E., Kuwata, A., et al., 2009a. Thermal expansion of forsterite at high pressures determined by in situ X-ray diffraction: the adiabatic geotherm in the upper mantle. Phys. Earth Planet. Inter. 174 (1), 86–92. Katsura, T., Yokoshi, S., Kawabe, K., Shatskiy, A., Manthilake, M.A.G.M., Zhai, S., Fukui, H., Hegoda, H.A.C.I., Yoshino, T., Yamazaki, D., Matsuzaki, T., Yoneda, A., Ito, E., Sugita, M., Tomioka, N., Hagiya, K., Nozawa, A., Funakoshi, K., 2009b. P–V– T relations of MgSiO3 perovskite determined by in situ X-ray diffraction using a large-volume high-pressure apparatus. Geophys. Res. Lett. 36, L01305. http:// dx.doi.org/10.1029/2008GL035658. King, S.D., Lee, C., van Keken, P., Leng, W., Zhong, S., Tan, E., Tosi, N., Kameyama, M., 2010. A community benchmark for 2D Cartesian compressible convection in the Earth’s mantle. Geophys. J. Int. 180, 73–87. http://dx.doi.org/10.1111/j.1365246X.2009.04413.x. Komabayashi, T., Hirose, K., Sugimura, E., Sata, N., Ohishi, Y., Dubrovinsky, L.S., 2008. Simultaneous volume measurements of post-perovskite and perovskite in MgSiO3 and their thermal equations of state. Earth Planet. Sci. Lett. 265, 515– 524. http://dx.doi.org/10.1016/j.epsl.2007.10.036. Leitch, A.M., Yuen, D.A., Sewell, G., 1991. Mantle convection with internal heating and pressure-dependent thermal expansivity. Earth Planet. Sci. Lett. 102, 213– 232. Li, L., Wentzcovitch, R.M., Weidner, D.J., da Silva, C.R.S., 2007. Vibrational and thermodynamic properties of forsterite at mantle conditions. J. Geophys. Res. 112, B05206. http://dx.doi.org/10.1029/2006JB004546. Liu, L., Stegman, D.R., 2011. Segmentation of the Farallon slab. Earth Planet. Sci. Lett. 311, 1–10. http://dx.doi.org/10.1016/j.epsl.2011.09.027. ˇ adek, O., C ˇ izkova, H., 2012. The effect Maierová, P., Chust, T., Steinle-Neumann, G., C of variable thermal diffusivity on kinematic models of subduction. J. Geophys. Res. 117, B07202. http://dx.doi.org/10.1029/2011JB009119. Manthilake, G.M., de Koker, N., Frost, D.J., McCammon, C.A., 2011a. Lattice thermal conductivity of lower mantle minerals and heat flux from Earth’s core. Proc. Natl. Acad. Sci. USA 108 (44), 17901–17904. http://dx.doi.org/10.1073/ pnas.1110594108. Manthilake, G.M., de Koker, N., Frost, D.J., 2011b. Thermal conductivity of CaGeO3 perovskite at high pressure. Geophys. Res. Lett. 38, L08301. http://dx.doi.org/ 10.1029/2011GL046882. Matyska, C., Yuen, D.A., 2005. The importance of radiative heat transfer on superplumes in the lower mantle with the new post-perovskite phase change. Earth Planet. Sci. Lett. 234, 71–81. Metsue, A., Tsuchiya, T., 2012. Thermodynamic properties of (Mg, Fe2+)SiO3 perovskite at the lower-mantle pressures and temperatures: an internally consistent LSDA+U study. Geophys. J. Int. 190 (1), 310–322. http://dx.doi.org/ 10.1111/j.1365-246X.2012.05511.x. Moré, J., 1978. The Levenberg–Marquardt algorithm: Implementation and theory. In: Watson, G. (Ed.), Numerical Analysis, Lecture Notes in Mathematics, vol. 630. Springer, Berlin/Heidelberg, pp. 105–116, doi:10.1007/BFb0067700. Nakagawa, T., Tackley, P.J., 2005. The interaction between the post-perovskite phase change and a thermo-chemical boundary layer near the core-mantle boundary. Earth Planet. Sci. Lett. 238, 204–216. Nakagawa, T., Tackley, P.J., Deschamps, F., Connolly, J.A.D., 2010. the influence of MORB and harzburgite composition on thermo-chemical mantle convection in a 3-D spherical shell with self-consistently calculated mineral physics. Earth Planet. Sci. Lett. 296, 403–412. http://dx.doi.org/10.1016/j.epsl.2010.05.026. Naliboff, J.B., Kellogg, L.H., 2006. Dynamic effects of a step-wise increase in thermal conductivity and viscosity in the lowermost mantle. Geophys. Res. Lett. 33, L12S09. http://dx.doi.org/10.1029/2006GL025717. O’Farrell, K.A., Lowman, J.P., 2010. Emulating the thermal structure of spherical shell convection in plane-layer geometry mantle convection models. Phys. Earth Planet. Inter. 182, 73–84. http://dx.doi.org/10.1016/j.pepi.2010.06.010. Ohta, K., Yagi, T., Taketoshi, N., Hirose, K., Komabayashi, T., Baba, T., Ohishi, Y., Hernlund, J., 2012. Lattice thermal conductivity of MgSiO3 perovskite and post-

58

N. Tosi et al. / Physics of the Earth and Planetary Interiors 217 (2013) 48–58

perovskite at the core-mantle boundary. Earth Planet. Sci. Lett. 349–350, 109– 115. http://dx.doi.org/10.1016/j.epsl.2012.06.043. O’Neill, C., Jellinek, A., Lenardic, A., 2007a. Conditions for the onset of plate tectonics on terrestrial planets and moons. Earth Planet. Sci. Lett. 261 (1), 20–32. O’Neill, C., Lenardic, A., Moresi, L., Torsvik, T., Lee, C., 2007b. Episodic precambrian subduction. Earth Planet. Sci. Lett. 262 (3), 552–562. Phillips, B.R., Coltice, N., 2010. Temperature beneath continents as a function of continental cover and convective wavelength. J. Geophys. Res. 115, B04408. http://dx.doi.org/10.1029/2009JB006600. Ribe, N.M., Stutzmann, E., Ren, Y., van der Hilst, R., 2007. Buckling instabilities of subducted lithosphere beneath the transition zone. Earth Planet. Sci. Lett. 254, 173–179. Samuel, H., Tosi, N., 2012. The influence of post-perovskite strength on the Earth’s mantle thermal and chemical evolution. Earth Planet. Sci. Lett. 323, 50–59. http://dx.doi.org/10.1016/j.epsl.2012.01.024. Schmeling, H., Marquart, G., Ruedas, T., 2003. Pressure- and temperature-dependent thermal expansivity and the effect on mantle convection and surface observables. Geophys. J. Int. 154, 2224–2229. http://dx.doi.org/10.1046/ j.1365-246X.2003.01949.x. Schuberth, B.S.A., Bunge, H.P., Ritsema, J., 2009. Tomographic filtering of highresolution mantle circulation models: can seismic heterogeneity be explained by temperature alone? Geophys. Geochem. Geosyst. 10 (5). http://dx.doi.org/ 10.1029/2009GC002401. Shahnas, H., Peltier, W.R., Wu, Z., Wentzcovitch, R.M., 2011. The high pressure electronic spin transition in iron: potential impacts upon mantle mixing. J. Geophys. Res. 116, B08205. http://dx.doi.org/10.1029/2010JB007965. Tan, E., Choi, E., Thoutireddy, P., Gurnis, M., Aivazis, M., 2011. On the location of plumes and lateral movement of thermochemical structures with high bulk modulus in the 3-D compressible mantle. Geochem. Geophys. Geosyst. 12 (7). http://dx.doi.org/10.1029/2011GC003665. Tosi, N., Yuen, D.A., 2011. Bent-shaped plumes and horizontal channel flow beneath the 660 km discontinuity. Earth Planet. Sci. Lett. 312, 348–359. http:// dx.doi.org/10.1016/j.epsl.2011.10.015. ˇ adek, O., 2010. Dynamical consequences in the lower mantle Tosi, N., Yuen, D.A., C with the post-perovskite phase change and strongly depth-dependent thermodynamic and transport properties. Earth Planet. Sci. Lett. 298, 229– 243. http://dx.doi.org/10.1016/j.epsl.2010.08.001. Touloukian, Y., Kirby, R.K., Taylor, R.E., Lee, T.Y.R., 1977. Thermal expansion of nonmeiallic solids. Thermophysical Properties of Matter, vol. 13. Plenum Press, New York. Tsuchiya, J., Tsuchiya, T., Wentzcovitch, R.M., 2005. Vibrational and thermodynamic properties of MgSiO3 postperovskite. J. Geophys. Res. 110, B02204. http:// dx.doi.org/10.1029/2004JB003409.

Umemoto, K., Wentzcovitch, R.M., Allen, P.B., 2006. Dissociation of mgsio3 in the cores of gas giants and terrestrial exoplanets. Science 311 (983). http:// dx.doi.org/10.1126/science.1120865. van der Meer, D.G., Spakman, W., van Hinsbergen, D.J.J., Amaru, M.L., Torsvik, T.H., 2010. Towards absolute plate motions constrained by lower-mantle slab remnants. Nat. Geosci. 3, 36–40. http://dx.doi.org/10.1038/ngeo708. Wallace, D., 1972. Thermodynamics of Crystals. Wiley, New York. Wentzcovitch, R.M., Wu, Z., Carrier, P., 2010. First principles quasiharmonic thermoelasticity of mantle minerals. Rev. Min. Geochem. 71, 99–128. http:// dx.doi.org/10.2138/rmg.2010.71.5. Westrenen, W., Li, J., Fei, Y., Frank, M., Hellwig, H., Komabayashi, T., Mibe, K., Minarik, W., Orman, J., Watson, H., Funakoshij, K., Schmidt, M., 2005. Thermoelastic properties of (Mg0.64Fe0.36)O ferropericlase based on in situ Xray diffraction to 26.7 GPa and 2173 K. Phys. Earth Planet. Int. 151 (1), 163–176. Wu, Z., Wentzcovitch, R.M., 2007. Vibrational and thermodynamic properties of wadsleyite: a density functional study. J. Geophys. Res. 112, B12202. http:// dx.doi.org/10.1029/2007JB005036. Wu, Z., Justo, J., da Silva, C., de Gironcoli, S., Wentzcovitch, R., 2009. Anomalous thermodynamic properties in ferropericlase throughout its spin crossover. Phys. Rev. B 80(1), 014,409. Xu, Y., Shankland, T.J., Linhardt, S., Rubie, D.C., Langenhorst, F., Klasinski, K., 2004. Thermal diffusivity and conductivity of olivine, wadsleyite and ringwoodite to 20 GPa and 1373 K. Phys. Earth Planet. Inter. 143–144, 321–336. http:// dx.doi.org/10.1016/j.pepi.2004.03.005. Yanagawa, T.K.B., Nakada, M., Yuen, D.A., 2005. Influence of lattice thermal conductivity on thermal convection with strongly temperature-dependent viscosity. Earth Planets Space 57, 15–28. Yoshida, M., Ogawa, M., 2004. Influence of two major phase transitions on mantle convection with moving and subducting plates. Earth Planets Space 56 (11), 1019–1033. Yu, Y.G., Wentzcovitch, R.M., 2006. Vibrational and thermodynamic properties of ringwoodite. J. Geophys. Res. 111, B12202. http://dx.doi.org/10.1029/ 2006JB004282. Yu, Y.G., Wu, Z., Wentzcovitch, R.M., 2008. a  b  c transformations in Mg2SiO4 in Earth’s transition zone. Earth Planet. Sci. Lett. 273, 115–122. http://dx.doi.org/ 10.1016/j.epsl.2008.06.023. Yuen, D.A., Monnereau, M., Hansen, U., Kameyama, M., Matyska, C., 2007. Dynamics of superplumes in the lower mantle. In: Yuen, D.A., Maruyama, S., Karato, S.I., Windley, B.F. (Eds.), Superplumes: Beyond Plate Tectonics, Geophysics and Geodesy. Springer, pp. 239–267. Zhong, S., 2006. Constraints on thermochemical convection of the mantle from plume heat flux, plume excess temperature, and upper mantle temperature. J. Geophys. Res. 111, B04409. http://dx.doi.org/10.1029/2005JB003972.