Michael H G Jacobs

Arie P Van Den Berg

Complex phase distribution and seismic velocity structure of the transition zone: convection model predictions for a magnesium-endmember olivine-pyroxene mantle

We have investigated the complex spatial distribution of mineral phases in the Earth's upper mantle, resulting from lateral variations of temperature in a convecting mantle, through numerical modelling of mantle convection. Our convection model includes a self-consistent thermodynamic description for an olivine-pyroxene composition in the SiO 2 , MgO system. The thermodynamic model is based on lattice vibrations and allows for the calculation of thermophysical properties as well as seismic wavespeeds.

Our modelling results show a complex structure in the behavior of physical properties, in particular the seismic shear wavespeed, in a depth range including the mantle transition zone, 400-700 km. We demonstrate that this behavior is related to the distribution of mineral phases in the olivine-pyroxene system. Especially near cold downwelling flows, representing subducting lithopheric plates, our model results show strong lateral variation of mineral phases and associated shear wavespeed. We show that, typically, pockets of contrasting mineral phases smaller than 100 km occur in subduction regions.

In line with current developments in seismic imaging of the mantle transition zone we have computed reflectivity profiles from the shear wavespeed distribution obtained from the convection results. We applied frequency filtering to the raw reflectivity data to investigate the requirements for resolving the heterogenous structure of the transition zone. Our results show that heterogeneous structure from contrasting mineral phase regions is resolved in the reflectivity profiles, for periods below twenty seconds which may be feasible in seismic imaging applications. This opens a perspective for detailed mapping of mineral phase distribution in the mantle offering new constraints on the thermal conditions of the transition zone region, thus providing valuable constraints for geodynamic models of the upper mantle.

Introduction

The mantle transition zone in the depth range 400-700 km is characterized by a large number of mineral phases especially in regions with strong lateral temperature variation, related to subducting lithospheric slabs. Seismological mapping of the existing phase boundaries provides constraints on the geotherm in such regions from identification with phase diagram data for candidate mantle materials. This way Boehler (1996) estimated the temperature at the spinel to post-spinel '660 km' phase boundary at 1900 ± 100 K, based on a pressure of 23.9 GPa from the PREM model [START_REF] Dziewonski | Preliminary reference earth model (PREM)[END_REF]. A recent example of this procedure, extended to lateral variation, is an estimate of temperature distribution in the heterogeneous D" region in the bottom region of the mantle, from high resolution imaging of seismic reflectors, identified as the post perovskite phase boundary (Van der [START_REF] Van Der Hilst | Seismostratigraphy and thermal structure of Earth's core-mantle boundary region[END_REF].

Lateral variation of the upper mantle transition zone has been investigated seismologically using underside reflections, observed as precursors of SS and PP phases [START_REF] Shearer | Experiments in migration processing of SS precursor data to image upper mantle discontinuity structure[END_REF][START_REF] Deuss | Seismic observations of splitting of the mid transition zone in Earths mantle[END_REF][START_REF] Deuss | The nature of the 660-kilometer discontinuity in Earths mantle from global seismic observations of PP precursors[END_REF], resulting in estimates of lateral variation of the depth of the '410' and '660 km' phase boundaries and regional variation in the visibility of the '520 km' phase boundary, associated with the wadsleyite (β) to ringwoodite (γ) transition in the (main) olivine component of the upper mantle. Recently [START_REF] Cao | Imaging the upper mantle transition zone with a generalized Radon transform of SS precursors[END_REF] applied seismic imaging methods, based on the Generalized Radon Transform (GRT) and related to the approach in (van der Hilst et al., 2007) and [START_REF] Wang | Imaging of the lowermost mantle (D") and the core-mantle boundary with SKKS coda waves[END_REF], to SS precursors originating from the transition zone for a tectonically undisturbed region of the pacific. This study shows clear resolution of the topography of the '410' and '660 km' boundaries and laterally varying intermediate structure interpreted as the '520 km' discontinuity that is shown to resolve locally in two separate reflections in the shorter period ranges of the signal. This could be an indication for additional phase boundaries related to A c c e p t e d M a n u s c r i p t mineral phases from the pyroxene system such as majorite and akimotoite that would be stable for possible temperature conditions in the transition zone.

The above examples illustrate that seismic visibility of mineral phase boundaries in complex configurations, as predicted for the mantle transition zone, is of great importance for constraining local temperature, which in turn is of great value for dynamical reconstructions of the upper mantle.

In this paper we investigate this issue, of the seismic expression of the mineral distribution in the transition zone, by means of numerical modelling of mantle convection, taking into account physically self-consistent thermophysical parameters and mineral phase diagram of the olivine-pyroxene system.

To our knowledge no mantle convection studies have adressed the prediction of heterogeneous seismic structure from the distribution of mineral phases in the transition zone, with the exception of [START_REF] Jacobs | The derivation of thermo-physical properties and phase equilibria of silicate materials from lattice vibrations: applications to convection in the earth's mantle[END_REF] restricted to an olivine mantle. Here our mantle convection model is based on the magnesium endmember of the olivine-pyroxene mineralogy [START_REF] Jacobs | Placing constraints on phase equilibria and thermophysical properties in the system MgO-SiO2 by a thermodynamically consistent vibrational method[END_REF]). The same model, further introduced in section 2.2 and the appendix, is also used to compute longitudinal and shear wave speed. Other workers have computed seismic wave speed distributions from mantle convection results without similar consistency between the mineral physics and the convection model [START_REF] Goes | Synthetic seismic signature of thermal mantle plumes[END_REF], (Schuberth et al., 2009a,b). We have chosen to apply the mineral physics model of Jacobs and de [START_REF] Jacobs | Placing constraints on phase equilibria and thermophysical properties in the system MgO-SiO2 by a thermodynamically consistent vibrational method[END_REF], in spite of its restricted composition, because of its degree of internal consistency and accuracy in representing available experimental data, that is superior to alternative mineral physics models that include a wider range of compositions, as is further illustrated below in section 2.2. The extension of our mineral physics model for a broader range of compositions, in particular iron, aluminium and calcium, is work in progress [START_REF] Jacobs | Thermodynamic mixing properties of olivine derived from lattice vibrations[END_REF].

The magnesium endmember olivine-pyroxene model is a simplified model relative to a more realistic, for instance pyrolitic, mantle composition [START_REF] Ringwood | Composition and petrology of the earth's mantle[END_REF]), yet it contains the most significant features of the mantle transition zone including the major phase transitions. The monovariant phase transitions of the magnesium-endmember system produce corresponding sharp transitions in physical properties including the seismic wave velocities. We use this endmember model for an exploratory test [START_REF] Ashida | Thermodynamic aspects of phase boundary among alpha-, beta-, and gamma-Mg 2 SiO 4[END_REF] to investigate the applicability of a self-consistent thermodynamic model in mantle convection computations, [START_REF] Bina | Phase transition Clapeyron slopes and transition zone seismic discontinuity topography[END_REF] to couple the self-consistently calculated sound wave velocities with mineral phase partitioning in the mantle and (3) to investigate the expression of the different mineral phases in seismic reflectivity traces.

In line with the above we have computed model predictions of the complex spatial distribution of mineral phases and corresponding seismic wave velocities resulting from our convection model for the whole mantle, in particular for the transition zone region. We convert the convection results (spatial distribution of pressure and temperature) to shear velocity values for the upper 1000 km of the mantle and correlate the resulting heterogeneous seismic velocity structures with corresponding variations of both temperature and mineral phase distribution. To put the results in perspective with recent developments in seismological reflectivity imaging of mantle structure, we compute 1-D radial reflectivity profiles from our convection results and consider the effect of limited seismic data bandwidth on the representation of the heterogenous structure by reflectivity time traces. This modelling procedure can help in investigating the potential of seismological resolution of the detailed structure of the transition zone. By focusing on the transition zone we avoid strong compositional heterogeneity of other mantle regions such as the D" zone at the bottom of the mantle, not accounted for in our model [START_REF] Xu | The effect of bulk composition and temperature on mantle seismic structure[END_REF].

Seismic expression of temperature variations in mantle convection models has been investigated

by [START_REF] Goes | Synthetic seismic signature of thermal mantle plumes[END_REF] In the next sections we first describe the convection model and underlying parameterization of the mineral physics model. This is followed by a presentation of the results from convection modelling and a discussion of the implications of the results and conclusions.

2 Description of the convection and thermodynamical model

Convection model

We apply a model of a compressible mantle, taking into account the effects of density increase with pressure and solid state phase transitions. The model is defined for a two-dimensional cylindrical domain geometry. The governing equations are,

∂ j (ρu j ) = 0 (1)
expressing mass conservation using a truncated anelastic liquid approximation (Jarvis and McKenzie, 1981[START_REF] Steinbach | Compressible convection in the earth's mantle: a comparison of different approaches[END_REF][START_REF] Steinbach | Numerische experimente zur Konvection in kompressiblen Medien[END_REF][START_REF] Ita | Sensitivity of convection with an endothermic phase change to the form of governing equations, initial conditions, boundary conditions, and equation of state[END_REF]. Conservation of momentum is expressed in the Stokes equation for a highly viscous fluid of infinite Prandtl number,

-∂ i ∆P + ∂ j τ ij + ∆ρg(r)e i = 0, i = 1, 2 (2)
where g(r) and e i are, respectively, the amplitude and direction cosines of the local radial gravity acceleration vector. The density perturbation ∆ρ in (2) is defined in terms of a reference profile ρ r , ρ(P, T ) = ρ r (P ) + ∆ρ(P, T )

where ρ r (P ) = ρ(P, T r (P )) and T r (P ) =< T (P ) > is a reference temperature profile computed from a lateral average. The dynamic pressure gradient in (2) is defined by subtracting ρ r (P )g i from the gradient of pressure ∂ i P in the Stokes equation.

Assuming that the effect of the rate of volume change, that scales with a bulk viscosity is negligable on the convective time scale [START_REF] Jarvis | Convection in a compressible fluid with infinite Prandtl number[END_REF], the deviatoric shear stress tensor for the compressible fluid is defined as,

τ ij = η ∂ i u j + ∂ j u i - 2 3 ∂ k u k δ ij (4)
where the summation convention is implied for repeated indices and η is the temperature and pressure dependent viscosity further specified in Table 1.

The energy transport equation used is,

ρc p DT Dt -αρgwT = ∂ j (k∂jT ) + τ ij ∂ j u i + ρH (5)
D/Dt is the material time derivative and the second term describes the effect of adiabatic heating.

The terms on the right specifiy thermal conduction, heating by viscous dissipation and radiogenic heating respectively. Further symbol definitions are given in Table 1. , 2006). The effect of heat exchange due to solid state phase transitions is not included in the formulation of the energy equation [START_REF] Christensen | Layered convection induced by phase transitions[END_REF]. In the extended boussinesq approximation of [START_REF] Christensen | Layered convection induced by phase transitions[END_REF] this effect is included in terms of parameterizations of the Clapeyron slope and the density contrast across the phase boundary, resulting in modified expressions for effective heat capacity and thermal expansivity in the energy equation. [START_REF] Nakagawa | Incorporating selfconsistently calculated mineral physics into thermochemical mantle convection simulations in a 3-D spherical shell and its influence on seismic anomalies in Earths mantle[END_REF] have extended this approach to a compressible model using tabulated thermophysical properties for a multicomponent model, where effective c p and α are expressed in the temperature-and pressure derivatives of the enthalpy. In our model calculations for the binary system MgO, SiO 2 we have A c c e p t e d M a n u s c r i p t not included such extensions that would require differentiation of the enthalpy across discontinuities at the monovariant phase transitions. We assume that this approximation will only have a minor influence on the overall dynamics and temperature structure of our model. This is further discussed in section 4.

The coupled model equations ( 1), ( 2) and ( 5) are solved using finite element methods. Domain discretization, boundary and initial conditions and computational aspects of the coupling between the model equations and the underlying mineral physics model are discussed in Appendix A.3.

Thermodynamical and mineralogical model

In our convection model the thermophysical properties density ρ, thermal expansivity α and specific heat c P , are derived from a thermodynamic formalism based on lattice dynamics. This formalism by Jacobs and de Jong (2005,2006,2007), further denoted as 'phonon-thermodynamic' method, incorporates the vibrational density of states, for the purely magnesian endmember model of the MgO -SiO 2 system. The overall composition is specified by the mole fractions of the olivine and pyroxene components, (Mg 2 SiO 4 ) 0.56 (MgSiO 3 ) 0.44 .

The absence of at least iron in our model system, relative to a natural system, of pyrolitic composition, prevailing in the mantle has three consequences. The first one is that density in our model system is lower than for a pyrolitic mantle, affecting sound wave velocities and the pressuredepth distribution. This is illustrated in Fig. 1 showing adiabatic profiles of seismic velocities, derived from the phonon-thermodynamic model, for an Earth-like mineral composition within the MgO -SiO 2 system. This figure shows that the main features of the mantle transition zone are reproduced by our mineral physics model based on a restricted composition (not including Al, Ca and Fe). The effect, of adding the major component Fe, on the seismic velocities of olivine, amounts to an overall reduction of the velocity by about 1-2%. Secondly, addition of more elements to the model system displaces phase boundaries. For instance the olivine-wadsleyite phase boundary will shift to an approximately 1 GPa lower pressure when Fe content is increased from zero to 10%.

The wadsleyite-ringwoodite transition shifts in a similar way to lower pressure for increased iron content. Thirdly, specific phase boundaries in the model system split up by the addition of other elements, enclosing multi-phase regions. For instance the olivine-wadsleyite phase boundary splits up in a field in which olivine and wadsleyite coexist when iron is added in our model. Although the restricted model cannot be used to represent the 'real Earth' it is sufficiently realistic for our present purpose of exploring the seismic expression of multiple phase regions in the transition zone generated by the mantle convection proces.

Despite the issue of compositional completeness of our present thermodynamic model we emphasize that the approach we follow here, is not restricted to a particular thermodynamic model. The mantle convection computations are set up in a flexible way such that different thermodynamical models, once available, can be applied, as illustrated further in section A.3 below. In future work this will enable us to refine our present calculations with an extended thermodynamic model containing more elements, and moreover to investigate effects of applying alternative thermodynamic models.

A thermodynamic model incorporating iron, which is constrained by phonon data and applicable in the range from zero Kelvin up to the melting point, is currently under development. Following the work of Jacobs and de Jong (2009), we anticipate that this model is particularly important for iron-bearing materials, which are characterized by antiferromagnetic-paramagnetic transitions at low temperatures, e.g. 67 K for fayalite and 12 K for ringwoodite. This requires a thermodynamic A main advantage of the application of the phonon-thermodynamic method in numerical models of mantle convection is that it enables one to compute the required thermo-physical parameters as well as the phase diagram in a single self-consistent framework. This is important because for many substances crucial experimental data, such as for thermal expansivity and bulk modulus, have been measured at conditions far away from those prevailing in the mantle transition zone.

Therefore, to obtain values for thermodynamic properties at mantle conditions, any physical model is forced to make large extrapolations of available thermodynamic data. For instance Jacobs and de Jong (2005b) showed, by using a calphad-type parameterization technique, that thermal expansivity and volume of gamma-Mg 2 SiO 4 (ringwoodite) could not be constrained with sufficient accuracy to enable a prediction of bulk sound velocity for gamma-(Mg 0.91 Fe 0.09 ) 2 SiO 4 to within 0.8%, twice the tomographic accuracy. Figure 12 of their work demonstrates that bulk sound velocity for this substance modeled by Fabrichnaya (1998), Stixrude and Bukowinsky (1993) and their own work are quite different, calling for a method which more tightly constrains thermodynamic data.

Such accuracy issues are also important because the representation of phase equilibrium boundaries is quite sensitive to variations in the Gibbs energy in pressure-temperature space. For instance a change in the Gibbs energy of akimotoite by about 500 J/mol ( 0.04 %) leads to a temperature shift of about 300 K of the quadruple point where wadsleyite, ringwoodite, akimotoite and stishovite are in equilibrium. [START_REF] Saxena | Earth mineralogical model: Gibbs free energy minimization in the system MgO -FeO -SiO 2[END_REF][START_REF] Saxena | Thermodynamic modelling of the Earths interior[END_REF] gives a detailed analysis of how uncertainties in thermo chemical properties affect the accuracy of the description of phase boundaries.

For all mineral phases present in the model phase diagram of Fig. 5, even for the phases stable at high pressure, experimental data, such as heat capacity, thermal expansivity and bulk modulus, are available at 1 bar pressure. These 1-bar data do not suffer from uncertainties associated with equations of state of pressure reference materials, as is the case for data obtained at high pressure.

It is therefore of great importance to constrain thermodynamic analyses of these phases with 1 bar experimental data. Such data are very well fitted by methods based on polynomial expressions for these properties, (Fabrichnaya et al., 2004), [START_REF] Saxena | Earth mineralogical model: Gibbs free energy minimization in the system MgO -FeO -SiO 2[END_REF], [START_REF] Holland | An internally consistent thermodynamic dataset for phases of petrological interest[END_REF].

However it has been demonstrated by Jacobs and Oonk (2000) that the methods based on polynomial approximations are difficult to use when thermodynamic properties are extrapolated to conditions prevailing in the lower mantle because of physical anomalies that occur in the computed thermal expansivity. Therefore the methods based on polynomials do not perform well, compared to alternative methods, for instance methods based on Mie-Grüneisen-Debye theory, when applied for the full pressure and temperature range of the Earth's mantle.

Significant progress in representing experimental data at lower mantle pressures has been made

by [START_REF] Stixrude | Thermodynamics of mantle minerals-I Physical properties[END_REF]. They used the Mie-Grüneisen-Debye method (MGD) to construct a thermodynamic database for the CFMAS system. However due to the simplified construction of the Vibrational Density of State (VDoS) inherent to the Debye model and the use of the quasi-harmonic approximation, experimental data at 1-bar conditions are less well represented by the MGD method compared to the polynomial approach. This is illustrated for forsterite in Fig. 2, where computed values of the specific heat at constant pressure are compared with experimental data.

Jacobs and de Jong (2005, 2007) introduced a method, based on Kieffers (1979) more elaborate model for the VDoS that also includes the effect of anharmonicity. This method, here labeled MGKa, combines an accurate representation of both 1-bar and high pressure (lower mantle) experimental Figure 2 shows for forsterite, that including anharmonicity combined with a realistic VDoS improves the representation of 1-bar heat capacity significantly. The effect of anharmonicity on vibrational frequencies has been measured for forsterite by [START_REF] Gillet | High-temperature thermodynamic properties of forsterite[END_REF][START_REF] Gillet | Anharmonic properties of Mg 2 SiO 4 -forsterite measured from the volume dependence of the Raman spectrum[END_REF]. Its effect on the heat capacity of forsterite is about 2.5% at 2000 K and 1 bar pressure and it propagates for about 0.4% in the Gibbs energy. Although this effect seems small at first sight, Jacobs and de Jong (2005) showed that it has a significant effect on location and slope of the phase boundary between forsterite and wadsleyite. Neglecting anharmonicity shifts this phase boundary by 2 GPa to lower pressures, whereas the slope decreases with almost a factor two. Not only forsterite, but also other mineral phases such as perovskite and akimotoite are affected by anharmonicity, impacting the position and slopes of phase boundaries as has been demonstrated by Jacobs and de Jong (2007).

In view of the demonstrated impact on the phase diagram, of both a detailed VDoS representation and of anharmonicity, in the lattice dynamics model, a unique characteristic of the MGK-a method, we consider this the most reliable method available to date for the computation of both the phase diagram and the thermo physical properties of the mantle convection model specified in ( 2)-( 5), for an olivine-pyroxene system. varying degree of self-consistency in the representation of thermophysical parameters and seismic wave velocities. The phonon-thermodynamic method is further outlined in the appendix.

Results of numerical modelling

We have modelled thermal convection in a 2-D cylindrical domain with an opening angle of ninety degrees. The ratio of the core radius and the outer radius R c /R is 0.54, similar to the value for Earth. The depth of the mantle is set to 3000 km. The model is both bottom and internally heated.

Uniform and constant radiogenic heating is applied, set to a value H = 2 • 10 -12 W/kg, slightly A newtonian, temperature and pressure dependent viscosity model is used with an exponential parameterization defined in Table 1. This indicates that seismic velocity is related in a simple way to temperature and pressure, in this region without phase transitions. This situation is quite different for the transition zone region between about 400 and 700 km depth as illustrated particularly in Fig. 3e,f. The large scale structure of the smooth temperature field associated with the vertical downwellings is dominant in the lateral variation but much of the small scale structure of the velocity anomalies is not reflected in the smooth variation of the temperature field.

A smooth relation between temperature and seismic wavespeed is underlying interpretations of seismic travel time tomography [START_REF] Goes | Synthetic seismic signature of thermal mantle plumes[END_REF], (Schuberth et al., 2010b). The poor correlation between small scale lateral variation of shear wave speed and tempeperature illustrated in Fig. 3 shows that assumption of such relation is not applicable for small scale variation in the depth range of the transition zone. A clear difference in correlation of seismic velocity and temperature, between regions inside and outside the transition zone is illustrated more directly in Fig. 4, showing discrete samples of wavespeed data obtained from several isobaric cross-sections at different sample depths, taken from the convection results. In our modelling results this transition occurs around 650 km depth as will be shown below. The fact that this transition occurs at depths greater than 520 km in our model is a consequence of the absence of iron in the magnesium endmember composition model used here. This composition related pressure shift of the β → γ transition has been calculated from the thermodynamic database of Saxena (1996) as 1.25 GPa. This effect is amplified in a shift to smaller depth, due to the increased pressure gradient associated with the density increase resulting from adding 10 % iron. Fig. 4 clearly illustrates that interpreting lateral variation in seismic velocity in terms of temperature variation only, will not produce correct results in a multi-phase environment like the mantle transition zone.

That strong lateral variation of physical properties, related to variation in mineral phase, can occur in the mantle transition zone, especially in regions of large lateral variation in temperature, is illustated in Figure 5, by three radial temperature profiles. The profiles are labeled by their angular coordinate, positioned in the zoom-in window around the cold downwellings of Fig. 3b,e. Adiabats with a range of potential surface temperatures of 1400-1800 K are drawn for reference in Figure 5.

This temperature range also includes the red profile, positioned between the two cold downwellings, that corresponds to a background reference profile, more representative for the mantle outside the cold disturbence of the downwelling flows. From the temperature excursion across several phase boundaries of the two low temperature geotherms, associated with the downwelling flows, it is clear that several mineral phases can exist at a given depth/pressure depending on the laterally varying temperature field. This is more directly illustrated in Fig. 6 showing the mineral phase distribution, frames (a,b) and corresponding lateral variation of the shear velocity δv s (c,d) for the two zoom-in windows indicated by black boxes in Fig. 3(a,d). Dashed lines in Fig. 6a correspond to the radial temperature profiles shown in the phase diagram of Fig. 5. Individual pockets of different mineral associations intersected by the radial profiles can be identified with the help of the phase diagram of Fig. 5. Fig. 6a shows a complex configuration of regional domains occupied by contrasting mineral phases, indicated by the same phase labels as in Fig. 5. The large negative temperature perturbations, of up to 1000 K, of the two cold downwellings, shift the local geotherms across several stability fields in the phase diagram in Fig. 5. This causes the strong lateral variation in the phase distribution shown in Fig. 6a.

Frame 6b shows similar lateral variation in the mineral phase distribution for a zoom-in window with a smaller, hot, mantle anomaly. The lateral variation of temperature is mainly reflected here in the variation of the thickness of the various mineral phase layers. The forementioned phase In the modelling results the bottom of the transition zone that represents the boundary between the upper and lower mantle consists of a layer of ringwoodite and perovskite, that transforms into a post-spinel assemblage of periclase and perovskite for pressures higher than the endothermic postspinel transition. The pyroxene component has transformed to perovskite at shallower level at the nearby strongly endothermic transition of akimotoite to perovskite (Chudinovskikh and Boehler, 2004, Jacobs and de Jong, 2007). The magnitude of the Clapeyron slope in (Jacobs and de Jong, 2007) of 4.1 MPa/K is smaller than the value of 6 MPa/K reported in [START_REF] Chudinovskikh | MgSiO 3 phase boundaries measured in the laser-heated diamond cell[END_REF]). The Clapeyron slope of the ringwoodite to postspinel transition of (Jacobs and de Jong,

2007

) is -1.9 MPa/K. The relatively small magnitude of the latter transition and the splitting of the phase transition to perovskite of the olivine and pyroxene reduce the dynamic effect of the composite endothermic phase boundary between the upper and lower mantle, in line with the weak interaction of the cold downwellings with these phase boundaries as shown in Fig. 3.

The radial profile at 37 o , halfway between the cold downwellings, is representative for the global horizontally averaged mantle. The corresponding geotherm in Fig. 5 is close to an adiabat of 1700K. The transition zone has a five layer configuration here with, from top to bottom 1) β + clino-enstatite, 2) β + majorite, 3) β + akimotoite, 4) β + pv 5) γ + perovskite. The transition zone is followed by the lower mantle consisting of a post spinel assemblage of perovskite MgSiO 3 + periclase MgO.

Fig. 6c,d show the lateral variation of the shear velocity for the corresponding zoom-in windows of frames (a) and (b). There is a strong correlation between the velocity anomalies and the mineral phase distribution, to the degree that individual pockets of contrasting mineral content can be identified in the anomaly pattern of the shear velocity. This is in contrast with the corresponding smooth temperature field distributions shown in Fig. 3(b,c).

The '410 km' seismic discontinuity is clearly visible in the velocity frames as the boundary between second and third phase layer from the surface. This phase boundary shows the well-known uplift near the core of the two cold downwellings, in agreement with the positive Clapeyron slope of the α → β phase transition [START_REF] Turcotte | Geodynamics[END_REF]. The phase boundary between layer 1 and 2, related to the transition from ortho-enstatite to clino-enstatie shows an even greater topography related to the high value of the Clapeyron slope, that increases for lower temperatures, as shown in Fig. 5. Table 2: Approximate two-way travel times and phase transitions related to the reflectivity maxima for the background mantle parts of Fig. 9.

The '520 km' seismic discontinuity, associated with the β → γ transition of the olivine component, is represented in our modelling results by the contrast between the β + ak (green in Fig. 6ab) and γ + pv (pink in Fig. 6ab) layers near the boundary between the upper and lower mantle. In the reflectivity traces of Fig. 9 this transition is visible as the second peak from the bottom near 270 seconds. The thin layer of γ + ak, replaced laterally by pv + β, as shown in Fig. 6ab is not resolved in the reflectivity traces.

The '520 km' discontinuity is observed in a peridotite mantle, of pyrolitic composition, roughly halfway between the α → β ('410 km') and γ → M gO + pv ('660 km') phase boundaries (Deus and Woodhouse, 2001, [START_REF] Cao | Imaging the upper mantle transition zone with a generalized Radon transform of SS precursors[END_REF]. Another reason for the relative proximity of the '520' and '660 km' boundaries, besides the forementioned absence of iron in our magnesian endmember model, is the rather high background temperature of the mantle model where the geotherm intersects the spinel to post-spinel boundary (γ + pv → MgO +pv) in Fig. 5 at a temperature of about 2050 K, higher than the experimental value 1900 ± 100 K in [START_REF] Boehler | High-pressure experiments and the phase diagram of lower mantle and core materials[END_REF]. This explanation by a temperature effect is in agreement with the occurence of two thicker pockets of ringwoodite and akimotoite (γ + ak), in the phase distribution shown in Fig. 6a, intersected by dashed profile lines, corresponding to the two colder geotherms labeled 35 o and 39 o in the phase diagram of Fig. 5. From Fig. 9 it can be concluded that the main mineral phase distribution of the transition zone, in the present magnesian endmember model, are clearly expressed in the reflectivity structure, at least for the period range below 20 seconds, well accessible in broadband seismological data, as illustrated in 9b.

Discussion and conclusions

Recent progress in mapping sharp contrasts in seismic shear wavespeed in the transition zone [START_REF] Deuss | The nature of the 660-kilometer discontinuity in Earths mantle from global seismic observations of PP precursors[END_REF][START_REF] Cao | Imaging the upper mantle transition zone with a generalized Radon transform of SS precursors[END_REF] calls for an assessment of phase related heterogeneity generated in mantle dynamical models. In order to investigate such heterogeneous structure of the mantle transition zone, resulting from lateral variation in temperature, we have modelled thermal convection in the Earth's mantle, applying a mineral physics model from which the complex phase diagram for the olivine-pyroxene system can be derived in a self-consistent manner.

Our mantle convection model includes thermophysical properties such as density, thermal expansivity and specific heat, for the magnesium endmembers of a olivine-pyroxene mantle. We have computed lateral variations in seismic shear wave velocities applying the same lattice vibrational model as used in the computation of thermophysical properties (Jacobs and de Jong, 2007).

Our computations are based on a restricted chemical composition of the model material, in particular without the approximately 10% iron content of a pyrolitic mantle composition (Ringwood, (Bina andHelffrich, 1994, Stixrude, 1997). In this sense our model results represents an endmember case for the computation of the seismic expression of phase boundaries in the transition zone.

Our model calculations were done with the simplifying assumption that the effect of latent heat exchange due to phase transitions can be neglected to first order in the large scale dynamics and global temperature distribution. The dynamical effect of the phase boundary on the mantle flow operates mainly through the phase related density contrast in the righthand side of the Stokes equation [START_REF] Bina | Phase transition Clapeyron slopes and transition zone seismic discontinuity topography[END_REF]. The localized temperature effect from the latent heat exchange produces a shift in the transition pressure that is estimated to be smaller than about 0.4 GPa. We note that this value is similar to the uncertainty in the experimental determination of phase boundaries [START_REF] Frost | The structure and sharpness of (Mg, Fe) 2 SiO 4 phase transformations in the transition zone[END_REF]. In view of the forementioned shift of the transition-depth of phase boundaries, resulting from absence of iron in our model, estimated at about 1 to 2 GPa for α → β and β → ringwoodite respectively, the effect of latent heat is acceptable.

In spite of these differences with a more realistic mantle composition, our model reproduces thermophysical properties, including seismic wavespeeds, as well as mineral phase boundaries to within experimental uncertainty, for the magnesium endmember of the olivine-pyroxene system. Selfconsistency of the underlying mineral physics model and high accuracy in representing experimental data, as discussed in section 2.2, is an important characteristic of our convection model, especially when applied to mantle regions with a large number of mineral phases with contrasting physical properties such as the upper mantle transition zone.

We have analysed the distribution of mineral phase and corresponding shear wavespeed, computed from the pressure-temperature fields obtained in the mantle convection calculations. In relatively undisturbed regions, away from cold downwellings and hot mantle plumes, the modelling results show heterogeneity mainly in the variation of layer thicknesses of the different phase associations. Here we observe interesting similarities with seismological results from SS and precursors from underside reflections of the transition zone [START_REF] Deuss | The nature of the 660-kilometer discontinuity in Earths mantle from global seismic observations of PP precursors[END_REF][START_REF] Cao | Imaging the upper mantle transition zone with a generalized Radon transform of SS precursors[END_REF]. The seismological results show lateral variations in amplitude and depth of seismic reflectors interpreted as topography of the '410' and '660 km' phase boundaries. The '520 km' boundary has a more regional character [START_REF] Deuss | Global Observations of Mantle Discontinuities Using SS and PP Precursors[END_REF] and frequency dependent signature that could be related to phase transitions in the pyroxene component [START_REF] Cao | Imaging the upper mantle transition zone with a generalized Radon transform of SS precursors[END_REF]. In our model results ringwoodite occurs in association with stishovite and with akimotoite in the core of a cold downwelling. Away from cold downwellings a thin layer of ringwoodite and akimotoite exists intermittantly replaced by an equally thin layer of wadsleyite and perovskite. This layering can be identified in the reflectivity traces computed from the wavespeed distributions. This configuration is related to the vicinity of the stability fields for Our main objective in this paper has been to investigate the heterogeneity due to the complex phase diagram of a magnesium endmember olivine-pyroxene system, for the mantle depth range including the transition zone, and relate the results to recent developments in seismological reflectivity imaging of the mantle.

We have shown that lateral variations of temperature result in a laterally heterogeneous phase distribution, especially near subduction regions. The results show that the distribution of seismic shear wavespeed is clearly dominated by the phase distribution rather than by the temperature field.

The major phase transitions of both the olivine and pyroxene subsystems show up clearly in the reflectivity traces computed from the shear velocity distribution in the signal bandwith with periods below 20 s, indicating that these phase boundaries can, in principle, be imaged by seismological methods. Geodynamical models that can generate the detailed structure of the transition zone region can therefore be helpful in the interpretaion of such seismological data.

Further progress requires that multi-phase regions, in particular related to a finite iron composition will be modelled. This requires the extension of our thermodynamic data base with iron, which is currently in progress (Jacobs and de Jong, 2009).

A(T, V ) = U ref (V st 0 ) + U st (V ) + A vib (T, V ) (6)
All thermodynamic properties are derived from the Helmholtz energy by applying classical thermodynamic theory. The system MgO -SiO 2 contains insulator materials and for such case the Helmholtz energy is expressed as the sum of a static lattice and a vibrational contribution:

The first term on the right-hand side of eqn. ( 6

P (T, V ) = P st (V ) + P vib (T, V ) (7)
At a selected condition of pressure and temperature the molar volume of the material is solved numerically from [START_REF] Connoly | Computation of phase equilibria by linear programming: A tool for geodynamic modeling and its application to subduction zone decarbonation[END_REF]. Density, ρ, in (3) is then calculated from the molar volume using the molar masses of the endmembers. Once the volume is known, all other thermodynamic properties follow by taking the appropriate differential of eqn. [START_REF] Chudinovskikh | MgSiO 3 phase boundaries measured in the laser-heated diamond cell[END_REF].

The calculation of the heat capacity, c P , in ( 5) is more involved and proceeds in several steps by making use of the well-known thermodynamic relation:

c P = c V + α 2 KV T (8)
The explicit expressions for the heat capacity at constant volume, c V , isothermal bulk modulus, K, and thermal expansivity, α, are derived from eqn (6) by taking the appropriate differentials:

c V = T ∂S ∂T V = -T ∂ 2 A ∂T 2 V (9) K = -V ∂P ∂V T = V ∂ 2 A ∂V 2 T (10) α = 1 K ∂P ∂T V = - 1 K ∂ ∂T ∂A ∂V T V (11)
All differentiations in [START_REF] Daessler | Two-dimensional thermo-kinetic model for the olivine-spinel phase transition in subducting slabs[END_REF][START_REF] Deuss | Seismic observations of splitting of the mid transition zone in Earths mantle[END_REF][START_REF] Deuss | The nature of the 660-kilometer discontinuity in Earths mantle from global seismic observations of PP precursors[END_REF] can be derived analytically from the explicit expression for the Helmholtz energy (Jacobs and de Jong, 2007).

In Kieffers (1979) model the directionally averaged longitudinal and transverse sound velocities are related to the angular vibrational frequency ω i,max by: for the energy equation (factor 4 in ( 15)). n g = 7 gauss integration points are used per element.

u i = 1 4 4πZV 3N A 1/3 ω i,max (12)
These evaluations are repeated two times per integration time step in a predictor/corrector scheme, resulting in the following number of thermo-table evaluation per time step.

N eval = (N stokes + N energy ) × 2 = (N e × n g + 3 × N e × 4 × n g ) × 2 = 26 × N e × n g [START_REF] Frost | The structure and sharpness of (Mg, Fe) 2 SiO 4 phase transformations in the transition zone[END_REF] This amounts to N eval = 1.2 • 10 7 table evaluations per integration time step for the finite element mesh size used in our model calculations where N e = 64000.

  who focus on large scale mantle plumes, and Schubert et al. (2009a,b), who focus on lateral variation in 3-D global convection models in comparison with evidence from global seismic tomography. In contrast we focus on modelling small scale heterogeneity in the transition zone and its expression in high resolution seismic imaging methods.

e = 1 /2e ij e ij 1 / 2

 112 Radial velocity aligned with gravity -e ij = ∂ j u i + ∂ i u j Strain rate tensor -e ij = e ij -1 3 e kk δ ij Strain rate deviator --Second invariant of deviator -η(T, z) = η 0 exp(cz -bT ) P,T dependent viscosity -b = ln(∆η T ), c = ln(∆η P ) --∆η T = 100, ∆η P = 10 -η 0 Viscosity scale value 10 23 Pa s τ ij = ηe ij Viscous stress tensor -α 0 Thermal expansivity scale value 3 • 10 -5

A c c e

  p t e d M a n u s c r i p t framework, which is adequate to reliably describe thermophysical properties from zero Kelvin up to the melting temperature and to pinpoint inconsistencies in existing laboratory data.

  A c c e p t e d M a n u s c r i p t data. As a consequence the MGK-a method is better constrained by, both 1-bar and high pressure, experimental data.

  A c c e p t e d M a n u s c r i p t lower than the chondritic value, in agreement with the smaller surface/volume ratio of the 2-D cylindrical domain, compared to the 3-D spherical case. Temperature is prescribed constant at the outer surface and the core mantle boundary with a temperature contrast ∆T = 3500 K across the mantle. Zero heatflux symmetry conditions are imposed on the radial boundaries. Computations were started from an initial adiabatic temperature field, with a potential temperature T p = 1600 K, truncated at the top-and bottom thermal boundary layers. The model equations were integrated in time for about 1 Gyr.

Fig. 3

 3 Fig. 3 shows results for a representative snapshot of the temperature distribution resulting from the convection model, at a model time of 837 Myr. The snapshot shown is characteristic for the model results with narrow cold downwellings, descending all the way to the core mantle boundary. These downwellings migrate along the outer boundary, similar to subducting lithospheric plates subject to trench migration and slab roll-back. The temperature field is shown in the top row of frames. The bottom row displays corresponding plots of the lateral variation of the seismic shear wavespeed. The temperature field is characterized by two colliding cold downwellings and a region of hot upwelling, in the area indicated by 1500 km wide zoom-in windows, outlined in black in frame (a), that correspond to frames (b) and (c). The cold downwellings show only limited dynamic interaction with the endothermic phase boundaries near 700 km depth, with a slight deflection of the right hand downwelling. This is related to the small value of the Clapeyron slope, -1.9 MPa/K, of the post-spinel phase transition at the bottom of the transition zone in our model. The temperature distribution of the zoom-in windows is smooth and without significant features related to the transition zone. Lateral variation of the seismic shear velocity δv s is shown in frames (d,e,f) that correspond to the temperature frames (a,b,c) respectively. The laterial variation is defined with respect to a smooth background profile, to bring out more clearly the discontinuities in the velocity, related to the different phase boundaries, δv s (x) = (v s (x)-< v s (r) >)/ < v s (r) >. The smooth background profile < v s (r) > is computed from the horizontally averaged values by applying a boxcar moving average filter with a 200 km halfwidth. The velocity and temperature variations shown in Fig.3 are closely correlated in the lower mantle region away from the D" region, where several lens-like structures of post-perovskite show up in the seismic velocity (Hernlund et al., 2005, van der Hilst et al., 2007).

A c c e p t e d M a n u s c r i p t

  Only the data samples taken at 797 km depth(27.5 GPa), in the shallow lower mantle occupied by post-spinel, show a continuous relation between shear velocity and temperature. The velocity-temperature plots of Fig. 4 can be interpreted with the phase-diagram for the magnesium endmember model applied in the convection model shown in Fig. 5. Horizontal dashed lines in Fig. 5 indicate the pressure values of the isobaric velocity-temperature plots in Fig. 4. The intersection of the 9.4 GPa isobar, at 300 km depth, with the phase boundary curve of the ortho-and clino enstatite transition at about 1750 K corresponds to the transition shown in the bottom velocity-temperature curve in Fig. 4. The loop in the sampling curves for 300 and 636 km depth of Fig. 4 is a result of the spread in pressure in the sample points that are located in a 20 km depth range around the indicated depth.In a similar way the intersection of the 21 GPa isobar, at 636 km depth, with the wadsleyiteringwoodite (β → γ) phase boundary near 1800 K corresponds to the velocity transition of the middle curve in Fig.4. The latter phase transition corresponds to the '520 km' discontinuity that has been the target of a growing number of seismological investigations(Shearer et al., Deus et al., 2001,2006[START_REF] Cao | Imaging the upper mantle transition zone with a generalized Radon transform of SS precursors[END_REF].

A c c e p t e d M a n u s c r i p t

  boundary between wadsleyite and ringwoodite, associated with the '520 km' discontinuity, occurs in two different settings here. First as the transition β + ak → γ + ak that is more prominent in the core of the cold downwellings. The local thickening of the β + ak layer is a consequence of the diamond shaped stability field of this mineral assemblage shown in Fig.5. Second, in warmer areas, as the transition β + pv → γ + pv. This results in an intermittent behavior of the β → γ boundary illustrated in 6a,b, which could be related to lateral variations in seismic observations of the '520 km' boundary reported by[START_REF] Deuss | Global Observations of Mantle Discontinuities Using SS and PP Precursors[END_REF] and[START_REF] Cao | Imaging the upper mantle transition zone with a generalized Radon transform of SS precursors[END_REF]. If such detail of the β → γ transition would indeed be observable in the seismic data this would provide strong constraints on lateral variation of the temperature field in the transition zone.

Fig. 7

 7 Fig. 7 shows a compilation of radial distributions of temperature (a,c) and seismic shear velocity (b,d). The different curves in each frame correspond to the profiles indicated by dotted lines in Fig. 6a,b. The black curves in frames (a) and (c) correspond to the horizontally averaged geotherm that is close to the 1700K adiabat in the transition zone.

A c c e p t e d M a n u s c r i p t

  Time in seconds Phase transition 200 α-olivine + clino-enstatite → β-olivine + clino-enstatite 225 β-olivine + clino-enstatite → β-olivine + majorite 240 β-olivine + majorite → β-olivine + akimotoite not-resolved β-olivine + akimotoite → γ-olivine + akimotoite 270 γ-olivine + akimotoite → γ-olivine + perovskite 280 γ-olivine + perovskite → periclase + perovskite

A c c e p t e d M a n u s c r i p t 1975 )

 1975 . We applied this model in our calculations because of the available corresponding thermophysical database (Jacobs and de Jong, 2007) that is unique in accurate prediction of both the phase diagram and material properties specific heat, thermal expansion and density, as discussed in section 2.2. The absence of iron in our model results in a shift to higher pressure of major phase transitions of the model olivine component, olivine-wadsleyite and wadsleyite-ringwoodite, associated with the '410 km' and '520 km' seismic discontinuity respectively (Bina and Helffrich, 1994, Jacobs and de Jong, 2005b). In contrast the spinel-postspinel boundary corresponding to the '660 km' discontinuity is less affected by the model iron content. The restricted composition also implies a lower model density and associated pressure gradient, which also results in a depth increase of the phase boundaries. The monovariant phase transitions of our mineral physics model do not have multi-phase regions characteristic for models containing solid solutions. This implies that phase transitions in our model are associated with discontinuities in physical properties, in contrast to a model containing 10% iron in solid solution that shows smoother transitions between phase domains

A c c e p t e d M a n u s c r i p t

  these phase associations, along an isobar of 22 GPa, in the phase domain of Fig.5. The lateral boundaries between these associations pinpoint the temperature of the 'quadruple point' at 22 GPa in Fig.5and the middle (black) data set shown in Fig.4. This example illustrates the possibility of lateral variation in the signature of seismic reflectors due to phase transitions of the pyroxene component.Other transitions in the pyroxene component that are clearly visible in the reflectivity traces of Fig.9are clinoenstatie to majorite, close to the '410 km' discontinuity, and the majorite to akimotoite transition about halfway in the transition zone.The strongest heterogeneity in our model results is found in regions with cold downwelling flows, simulating lithospheric subduction, where the geotherm deviates by up to 1000 K from the horizontally averaged background value. This heterogeneity shows up as distinct pockets of mineral phase with corresponding, contrasting seismic velocity distribution, where the outlines are clearly visible in the reflectivity traces. It should be noted here that our modelling results have been obtained assuming equilibrium thermodynamics, thereby neglecting possible effects of phase transition kinetics that may become significant in the coldest cores of old subducting lithopheric slabs[START_REF] Daessler | Two-dimensional thermo-kinetic model for the olivine-spinel phase transition in subducting slabs[END_REF][START_REF] Vacher | Numerical tests on the seismic visibility of metastable minerals in subduction zones[END_REF].

  ) is a reference energy contribution at zero Kelvin and zero pressure for a substance in which no vibrational motion is present. The molar volume of this static lattice is denoted by V st 0 . The second term represents the change in energy resulting from the change in volume of the static lattice from V st 0 to V . This term is calculated by applying an equation of state, such as a third order Birch-Murnaghan equation of state or the universal equation of state developed by Vinet et al. (1987). The last term represents the vibrational contribution to the Helmholtz energy and is expressed in vibrational mode frequencies and their mode Grüneisen parameters. Taking the isothermal volume derivative of eqn. (6) results in the expression for the equation of state, which is expressed as the sum of the static lattice and vibrational (thermal) contribution:

Figure 1 :Figure 2 :

 12 Figure 1: Pressure-depth dependence of the seismic, longitudinal (v p ), shear (v s ) and bulk (v b ) wavespeeds, after (Jacobs and de Jong, 2009). Solid red curves are model predictions for mantle adiabats (1420 K) and a similar mineral composition as in our present model. Grey fields represent wavespeeds for compositions ranging from Mg 2 SiO 4 to MgSiO 3 . Green curves represent a model with (Mg 0.9 Fe 0.1 ) 2 SiO 4 . The dashed green curve is from a model based on calphad type parameterizations for the same olivine composition (Jacobs and de Jong, 2005b). Discrete symbols represent the PREM model.

Figure 3 :February 12 , 2011 27 849Figure 4 :

 312274 Figure 3: (a) Full domain temperature snapshot showing two colliding cold downwellings and a hot mantle upwelling. Black boxes represent zoom-in windows corresponding to frames (b,c). (b) zoom-in temperature of the cold downwellings flow. (c) zoom-in temperature of the shallow head of a hot upwelling. (d) Full domain distribution of the lateral variation of the shear wavespeed δv s with respect to the smooth background profile. (e) zoom-in of δv s of the cold downwelling. (f) zoom-in of δv s of the hot mantle plume head.

A c c e p t e d M a n u s c r i p tFebruary 12 , 2011 28 850Figure 5 :February 12 , 2011 29 851Figure 6 :

 1228512296 Figure 5: Phase diagram of the olivine-pyroxene magnesium endmember model used in the convection model. Adiabats of potential temperature T p between 1400 and 1800 K are plotted for reference. Solid curves represent radial temperature profiles, with angular position shown in the legend, and corresponding to the dashed lines in Fig. 6a. Horizontal dashed lines, labeled with a pressure value, correspond to the isobaric cross sections shown in Fig. 4.

A c c e p t e d M a n u s c r i p tFebruary 12 , 2011 31 853Figure 8 :

 12318 Figure 8: Spectral analysis of a single unfiltered reflectivity time trace plotted shown in frame (a), corresponding to the highlighted red trace in Fig. 9(a). (b) powerspectrum of the time trace in (a), showing that the peak power is below the frequency of 0.1 Hz, corresponding to wavelength's from 50 to 60 km in the transition zone.

Figure 9 :

 9 Figure 9: Reflectivity-time traces derived from the s-wavespeed distribution corresponding to the zoom-in window shown in Fig. 6a,c. The corresponding s-wavespeed distribution is shown in Fig. 10. The red colored trace corresponds to Fig. 8. (a) shows the raw data derived from the convection model. (b) and (c) show corresponding time traces after application of low-pass Butterworth filters with corner-period settings of 20 sec. and 40 sec. respectively.

Table 1 :

 1 Physical parameters

	K -1

Table 2

 2

lists reflectivity maxima corresponding to several phase boundaries, mostly related to the pyroxene component, that can be identified in locations, away from the irregular structures related to the two cold downwellings in the mantle region shown.

Acknowledgments

We acknowledge constructive reviews from Takashi Nakagawa and Yuji Higo that were helpful for improving the manuscript. Rob van der Hilst is thanked for making available unpublished results.

Computational resources for this work were made available by the Netherlands Research Center for

Integrated Solid Earth Science (ISES 3.2.5). Collaboration between A. van den Berg and M. Jacobs has been supported through ISES project ME-2.7.

A Thermodynamic formalism and computation of thermophysical parameters and seismic wave velocities A.1 Background of the thermodynamical model

In our modelling experiments thermophysical properties ρ, α, c P and the seismic wave velocities were derived from a formalism based on thermodynamics and lattice dynamics. This phononthermodynamic formalism was developed in (Jacobs andde Jong, 2005, 2007). There are several important differences between this formalism and recently published alternative formulations and corresponding databases by [START_REF] Stixrude | Thermodynamics of mantle minerals-I Physical properties[END_REF], [START_REF] Connoly | Computation of phase equilibria by linear programming: A tool for geodynamic modeling and its application to subduction zone decarbonation[END_REF] and Piazonni et al., (2007). The latter databases employ the Debye model to obtain thermal pressure, but polynomial expressions for heat capacity and thermal expansivity are used in the data base of [START_REF] Piazzoni | A mineralogical model for density and elasticity of the Earths mantle[END_REF].

An important difference is that the formalism used by Jacobs and de Jong (2005,2007) incorporates more details of the phonon spectrum than the Debye formalism allows.

Such details are essential because the Debye model was shown to result in insufficient accuracy in the heat capacity from 0 K up to the melting point of forsterite (Jacobs and de Jong, 2005a).

That compelled them to employ [START_REF] Kieffer | Thermodynamics and lattice vibrations of minerals: 3 Lattice dynamics and an approximation for minerals with application to simple substances and framework silicates[END_REF] model to represent the experimental data associated with the vibrational density of state (VDoS), such as vibrational frequency-temperature-pressure measurements derived from Raman and infrared spectroscopy. Kieffer's model appeared also to be required for other minerals such as MgSiO 3 perovskite.

Secondly, the formalism of Jacobs and de Jong (2005,2007) incorporates intrinsic mode anharmonicity. These data became available for a number of mantle materials, such as forsterite, perovskite and akimotoite through the work of [START_REF] Gillet | High-temperature thermodynamic properties of forsterite[END_REF], [START_REF] Gillet | A thermodynamic model for MgSiO 3 -perovskite derived from pressure, temperature and volume dependence of the Raman mode frequencies[END_REF] and Reynard and [START_REF] Rubie | High-pressure, high-temperature Raman spectroscopic study of ilmenite-type MgSiO 3[END_REF] respectively. Jacobs and de Jong (2005,2007) demonstrated that intrinsic anharmonicity significantly affects slope and position of phase boundaries in the magnesium-olivine and magnesium-pyroxene system.

Thirdly the formalism is constructed in such a way that, at zero Kelvin, static lattice properties are calculated, which can be constrained by or compared with ab initio methods.

The fourth difference is that the phonon-thermodynamic formalism incorporates electronic and magnetic properties, which is shown by Jacobs and de Jong (2009) to be a requirement for accurate representation of Fe bearing mantle materials. They showed that the vibrational formalism is suitable to predicting mixing properties in olivine. Their incorporation of iron in the system MgO -SiO 2 is work in progress. Because the phonon-thermodynamic model incorporates more physical properties relative to polynomial models used by e.g. Fabrichnaya (1998), Holland and Powell (1998), [START_REF] Saxena | Earth mineralogical model: Gibbs free energy minimization in the system MgO -FeO -SiO 2[END_REF], Jacobs et al. (2005b), the thermodynamic description of substances is more unambiguously constrained.

A.2 Computation of thermophysical parameters with the phonon-thermodynamic model

The central equation of the phonon-thermodynamic model is the expression for the Helmholtz energy of a material. For insulating materials such as present in the system MgO -SiO 2 , The Helmholtz In ( 12) u 1 and u 2 represent the sound velocities of the transverse (shear) waves and u 3 that of the Voigt-Reuss-Hill longitudinal wave. The shear wavespeed is calculated according to:

The P-T phase diagram for a specific composition in the MgO -SiO 2 system is derived from the Gibbs energy of each mineral phase. The Gibbs energy determines the stability of the phase and is derived from the Helmholtz energy as:

As stated above, the volume of the mineral phase at a specific pressure and temperature is derived from [START_REF] Connoly | Computation of phase equilibria by linear programming: A tool for geodynamic modeling and its application to subduction zone decarbonation[END_REF] resulting in a value for the P • V term in [START_REF] Fabrichnaya | The assessment of thermodynamic parameters for solid phases in the Fe-Mg-O and Fe-Mg-Si-O systems[END_REF]. This value combined with the value for the Helmholtz energy derived from ( 6) results in the value of the Gibbs energy of each phase. From the Gibbs energy values of each phase and applying general rules of classical thermodynamics it is straightforward to derive which phase assemblage is stable at given (P, T ), [START_REF] Connoly | Computation of phase equilibria by linear programming: A tool for geodynamic modeling and its application to subduction zone decarbonation[END_REF].

A.3 Computational aspects

The model equations ( 1),( 2) and ( 5) are solved on a two dimensional quarter circle cylindrical domain covering the full depth of the mantle, using finite element methods based on the SEPRAN finite element package (Segal and Praagman, 2000).

Free slip impermeable boundaries are prescribed for the Stokes equation ( 2). We apply cartesian coordinates and local transformations to implement the mechanical boundary conditions on the curved boundaries [START_REF] Cuvelier | Finite Element Methods and Navier-Stokes Equations[END_REF]. Temperature is prescribed at the bottom and top horizontal surfaces with a driving temperature contrast ∆T = 3500 K. Zero heatflux symmetry conditions apply on the radial boundaries.

Approximately 2 • 10 5 nodal points are used spanning 64000 quadratic elements for the velocity field using 227 and 567 nodal points in the radial and horizontal direction respectively. A fourfold subdivision is applied on the vertices of the 7-point velocity elements to obtain linear temperature elements for the energy equation (van den [START_REF] Van Den Berg | The Effects of a composite non-Newtonian and Newtonian rheology on Mantle convection[END_REF]). Grid refinement is applied with nodal point spacing reduced to approximately 6 km in the depth interval containing the transition zone.

An integrated finite element method in terms of the velocity and pressure fields is combined with a direct solver for the solution of the Stokes equation where renumbering of the degrees of freedom is applied to reduce the profile of the stiffness matrix.

The energy equation is integrated in time with a predictor-corrector method, applying a Courant type time step control (van den Berg et al., 1993).

Necessary thermo-physical properties density, thermal expansivity and specific heat, are computed through a general interface procedure using equidistant tabulation and bi-linear interpolation of these properties in the P, T domain relevant for the mantle 300 K < T < 4000 K, δT = 2 K, 0 < P < 140 GPa, δP = 0.066 GPa. This represents a flexible and efficient tool that facilitates the use of alternative thermodynamic data bases (van den [START_REF] Van Den Berg | Numerical models of mantle convection based on thermodynamic data for the MgO-SiO2 olivine-pyroxene system[END_REF][START_REF] Jacobs | The derivation of thermo-physical properties and phase equilibria of silicate materials from lattice vibrations: applications to convection in the earth's mantle[END_REF].

P, T domain tabulation is essential since pointwise calculation of thermo-physical properties is not feasible in view of the large number of evaluations of these properties. This becomes clear from the following assessment of the number of evaluations involved per integration time step, for a general model case with N e stokes elements.

For the Stokes equation only ρ needs to be computed. For the energy equation the values of ρ, c P and α must be computed (factor 3 in ( 15)). A four times higher number of elements is used