Mathematical Background¶
Here is a bit of background on the methods used to calculate thermoelastic and thermodynamic properties in BurnMan. More detail can be found in the cited papers.
Endmember Properties¶
Calculating Thermoelastic Properties¶
To calculate the bulk (\(K\)) modulus, shear modulus (\(G\)) and density (\(\rho\)) of a material at a given pressure (\(P\)) and temperature (\(T\)), optionally defined by a geotherm) and determine the seismic velocities (\(V_S, V_P, V_\Phi\)), one uses an Equation of State (EoS). Currently the following EoSs are supported in BurnMan:
 BirchMurnaghan finitestrain EoS (excludes temperature effects, [Poi91]),
 BirchMurnaghan finitestrain EoS with a MieGrüneisenDebye thermal correction, as formulated by [SLB05].
 BirchMurnaghan finitestrain EoS with a MieGrüneisenDebye thermal correction, as formulated by [MBR+07].
 Modified Tait EoS (excludes temperature effects, [HuangChow74]),
 Modified Tait EoS with a pseudoEinstein model for thermal corrections, as formulated by [HollandPowell11].
 CompensatedRedlichKwong for fluids, as formulated by [HP91].
To calculate these thermoelastic parameters, the EoS requires the user to input the pressure, temperature, and the phases and their molar fractions. These inputs and outputs are further discussed in User input.
BirchMurnaghan (isothermal)¶
The BirchMurnaghan equation is an isothermal Eulerian finitestrain EoS relating pressure and volume. The negative finitestrain (or compression) is defined as
where \(V\) is the volume at a given pressure and \(V_0\) is the volume at a reference state (\(P = 10^5\) Pa , \(T\) = 300 K). The pressure and elastic moduli are derived from a thirdorder Taylor expansion of Helmholtz free energy in \(f\) and evaluating the appropriate volume and strain derivatives (e.g., [Poi91]). For an isotropic material one obtains for the pressure, isothermal bulk modulus, and shear modulus:
Here \(K_0\) and \(G_0\) are the reference bulk modulus and shear modulus and \(K_0^\prime\) and \({G}^\prime_{0}\) are the derivative of the respective moduli with respect to pressure.
BurnMan has the option to use the secondorder expansion for shear modulus by dropping the \(f^2\) terms in these equations (as is sometimes done for experimental fits or EoS modeling).
Modified Tait (isothermal)¶
The Modified Tait equation of state was developed by [HuangChow74]. It has the considerable benefit of allowing volume to be expressed as a function of pressure. It performs very well to pressures and temperatures relevant to the deep Earth [HollandPowell11].
MieGrüneisenDebye (thermal correction to BirchMurnaghan)¶
The Debye model for the Helmholtz free energy can be written as follows [MBR+07]
where \(\theta\) is the Debye temperature and \(\gamma\) is the Grüneisen parameter.
Using thermodynamic relations we can derive equations for the thermal pressure and bulk modulus
The thermal shear correction used in BurnMan was developed by [HamaSuito98]
The total pressure, bulk and shear moduli can be calculated from the following sums
This equation of state is substantially the same as that in SLB2005 (see below). The primary differences are in the thermal correction to the shear modulus and in the volume dependences of the Debye temperature and the Gruneisen parameter.
HP2011 (thermal correction to Modified Tait)¶
The thermal pressure can be incorporated into the Modified Tait equation of state, replacing \(P\) with \(P\left(P_{\textrm{th}}  P_{\textrm{th0}}\right)\) in Equation (5) [HollandPowell11]. Thermal pressure is calculated using a MieGrüneisen equation of state and an Einstein model for heat capacity, even though the Einstein model is not actually used for the heat capacity when calculating the enthalpy and entropy (see following section).
\(\Theta\) is the Einstein temperature of the crystal in Kelvin, approximated for a substance \(i\) with \(n_i\) atoms in the unit formula and a molar entropy \(S_i\) using the empirical formula
SLB2005 (for solids, thermal)¶
Thermal corrections for pressure, and isothermal bulk modulus and shear modulus are derived from the MieGrüneisenDebye EoS with the quasiharmonic approximation. Here we adopt the formalism of [SLB05] where these corrections are added to equations (2)–(4):
The \(\Delta\) refers to the difference in the relevant quantity from the reference temperature (300 K). \(\gamma\) is the Grüneisen parameter, \(q\) is the logarithmic volume derivative of the Grüneisen parameter, \(\eta_{S}\) is the shear strain derivative of the Grüneisen parameter, \(C_V\) is the heat capacity at constant volume, and \(\mathcal{U}\) is the internal energy at temperature \(T\). \(C_V\) and \(\mathcal{U}\) are calculated using the Debye model for vibrational energy of a lattice. These quantities are calculated as follows:
where \(\theta\) is the Debye temperature of the mineral, \(\nu\) is the frequency of vibrational modes for the mineral, \(n\) is the number of atoms per formula unit (e.g. 2 for periclase, 5 for perovskite), and \(R\) is the gas constant. Under the approximation that the vibrational frequencies behave the same under strain, we may identify \(\nu/\nu_0 = \theta/\theta_0\). The quantities \(\gamma_0\), \(\eta_{S0}\) \(q_0\), and \(\theta_0\) are the experimentally determined values for those parameters at the reference state.
Due to the fact that a planetary mantle is rarely isothermal along a geotherm, it is more appropriate to use the adiabatic bulk modulus \(K_S\) instead of \(K_T\), which is calculated using
where \(\alpha\) is the coefficient of thermal expansion:
There is no difference between the isothermal and adiabatic shear moduli for an isotropic solid. All together this makes an eleven parameter EoS model, which is summarized in the Table below. For more details on the EoS, we refer readers to [SLB05].
User Input  Symbol  Definition  Units 

V_0  \(V_{0}\) 

m \(^{3}\) mol \(^{1}\) 
K_0  \(K_{0}\)  Isothermal bulk modulus at P=10^5 Pa, T = 300 K  Pa 
Kprime_0  \(K^\prime_0\)  Pressure derivative of \(K_{0}\)  
G_0  \(G_{0}\)  Shear modulus at P = \(10^5\) Pa, T = 300 K  Pa 
Gprime_0  \(G^\prime_0\)  Pressure derivative of \(G_{0}\)  
molar_mass  \(\mu\)  mass per mole formula unit  kg \(\mathrm{mol}^{1}\) 
n  n  number of atoms per formula unit  
Debye_0  \(\theta_{0}\)  Debye Temperature  K 
grueneisen_0  \(\gamma_{0}\)  Grüneisen parameter at P = \(10^5\) Pa, T = 300 K  
q0  \(q_{0}\)  Logarithmic volume derivative of the Grüneisen parameter  
eta_s_0  \(\eta_{S0}\)  Shear strain derivative of the Grüneisen parameter 
This equation of state is substantially the same as that of the MieGruneisenDebye (see above). The primary differences are in the thermal correction to the shear modulus and in the volume dependences of the Debye temperature and the Gruneisen parameter.
CompensatedRedlichKwong (for fluids, thermal)¶
The CORK equation of state [HP91] is a simple virialtype extension to the modified RedlichKwong (MRK) equation of state. It was designed to compensate for the tendency of the MRK equation of state to overestimate volumes at high pressures and accommodate the volume behaviour of coexisting gas and liquid phases along the saturation curve.
Calculating Thermodynamic Properties¶
So far, we have concentrated on the thermoelastic properties of minerals. There are, however, additional thermodynamic properties which are required to describe the thermal properties such as the energy, entropy and heat capacity. These properties are related by the following expressions:
where \(P\) is the pressure, \(T\) is the temperature and \(\mathcal{E}\), \(\mathcal{F}\), \(\mathcal{H}\), \(\mathcal{S}\) and \(V\) are the molar internal energy, Helmholtz free energy, enthalpy, entropy and volume respectively.
HP2011¶
The heat capacity at one bar is given by an empirical polynomial fit to experimental data
The entropy at high pressure and temperature can be calculated by differentiating the expression for \(\mathcal{G}\) with respect to temperature
Finally, the enthalpy at high pressure and temperature can be calculated
SLB2005¶
The Debye model yields the Helmholtz free energy and entropy due to lattice vibrations
Property modifiers¶
The thermodynamic models above consider the effects of strain and quasiharmonic lattice vibrations on the free energies of minerals at given temperatures and pressures. There are a number of additional processes, such as isochemical orderdisorder and magnetic effects which also contribute to the total free energy of a phase. Corrections for these additional processes can be applied in a number of different ways. Burnman currently includes implementations of the following:
 Linear excesses (useful for DQF modifications for [HollandPowell11])
 Tricritical Landau model (two formulations)
 BraggWilliams model
 Magnetic excesses
In all cases, the excess Gibbs free energy \(\mathcal{G}\) and first and second partial derivatives with respect to pressure and temperature are calculated. The thermodynamic properties of each phase are then modified in a consistent manner; specifically:
Subscripts \(_o\) and \(_m\) indicate original properties and modifiers respectively. Importantly, this allows us to stack modifications such as multiple Landau transitions in a simple and straightforward manner. In the burnman code, we add property modifiers as an attribute to each mineral as a list. For example:
from burnman.minerals import SLB_2011
stv = SLB_2011.stishovite()
stv.property_modifiers = [
['landau',
{'Tc_0': 4250.0, 'S_D': 0.012, 'V_D': 1e09}]]
['linear',
{'delta_E': 1.e3, 'delta_S': 0., 'delta_V': 0.}]]
Each modifier is a list with two elements, first the name of the modifier type, and second a dictionary with the required parameters for that model. A list of parameters for each model is given in the following sections.
Linear excesses (linear)¶
A simple linear correction in pressure and temperature. Parameters are ‘delta_E’, ‘delta_S’ and ‘delta_V’.
Tricritical Landau model (landau)¶
Applies a tricritical Landau correction to the properties of an endmember which undergoes a displacive phase transition. These transitions are not associated with an activation energy, and therefore they occur rapidly compared with seismic wave propagation. Parameters are ‘Tc_0’, ‘S_D’ and ‘V_D’.
This correction follows [Putnis92], and is done relative to the completely ordered state (at 0 K). It therefore differs in implementation from both [SLB11] and [HollandPowell11], who compute properties relative to the completely disordered state and standard states respectively. The current implementation is preferred, as the excess entropy (and heat capacity) terms are equal to zero at 0 K.
If the temperature is above the critical temperature, Q (the order parameter) is equal to zero, and the Gibbs free energy is simply that of the disordered phase:
If temperature is below the critical temperature, Q is between 0 and 1. The gibbs free energy can be described thus:
Tricritical Landau model (landau_hp)¶
Applies a tricritical Landau correction similar to that described above. However, this implementation follows [HollandPowell11], who compute properties relative to the standard state. Parameters are ‘P_0’, ‘T_0’, ‘Tc_0’, ‘S_D’ and ‘V_D’.
It is worth noting that the correction described by [HollandPowell11] has been incorrectly used throughout the geological literature, particularly in studies involving magnetite (which includes studies comparing oxygen fugacities to the FMQ buffer (due to an incorrect calculation of the properties of magnetite). Note that even if the implementation is correct, it still allows the order parameter Q to be greater than one, which is physically impossible.
We include this implementation in order to reproduce the dataset of [HollandPowell11]. If you are creating your own minerals, we recommend using the standard implementation.
If the temperature is above the critical temperature, Q (the order parameter) is equal to zero. Otherwise
The second derivatives of the Gibbs free energy are only nonzero if the order parameter exceeds zero. Then
BraggWilliams model (bragg_williams)¶
The BraggWilliams model is essentially a symmetric solid solution model between endmembers, with an excess configurational entropy term predicted on the basis of the specifics of orderdisorder in the mineral, multiplied by some empirical factor. Expressions for the excess Gibbs free energy can be found in [HP96]. Parameters are ‘deltaH’, ‘deltaV’, ‘Wh’, ‘Wv’, ‘n’ and ‘factor’.
Magnetic model (magnetic_chs)¶
This model approximates the excess energy due to magnetic ordering. It was originally described in [CHS87]. The expressions used by BurnMan can be found in [Sun91]. Parameters are ‘structural_parameter’, ‘curie_temperature’[2] (zero pressure value and pressure dependence) and ‘magnetic_moment’[2] (zero pressure value and pressure dependence).
Calculating Solid Solution Properties¶
Many minerals can exist over a finite region of composition space. These spaces are bounded by endmembers (which may themselves not be stable), and each individual mineral can then be described as a solid solution of those endmembers. At an atomic level, different elements substitute for one another on distinct crystallographic sites in the structure. For example, low pressure silicate garnets have two distinct sites on which mixing takes place; a dodecahedral site (of which there are three per unit cell on an eightcation basis) and octahedral site (of which there are two per unit cell). A third tetrahedral cation site (three per unit cell) is usually assumed to be occupied solely by silicon, and therefore can be ignored in solid solution calculations. The chemical formula of many low pressure garnets exist within the solid solution:
We typically calculate solid solution properties by appropriate differentiation of the Gibbs Free energy, where
Implemented models¶
Ideal solid solutions¶
A solid solution is not simply a mechanical mixture of its constituent endmembers. The mixing of different elements on sites results in an excess configurational entropy
where \(s\) is a site in the lattice \(M\), \(c\) are the cations mixing on site \(s\) and \(\nu\) is the number of \(s\) sites in the formula unit. Solid solutions where this configurational entropy is the only deviation from a mechanical mixture are termed ideal. From this expression, we can see that
Symmetric solid solutions¶
Many real minerals are not well approximated as ideal solid solutions. Deviations are the result of elastic and chemical interactions between ions with different physical and chemical characteristics. Regular (symmetric) solid solution models are designed to account for the simplest form of deviations from ideality, by allowing the addition of excess enthalpies, nonconfigurational entropies and volumes to the ideal solution model. These excess terms have the matrix form [DPWH07]
where \(p\) is a vector of molar fractions of each of the \(n\) endmembers and \(W\) is a strictly uppertriangular matrix of interaction terms between endmembers. Excesses within binary systems (\(i\)\(j\)) have a quadratic form and a maximum of \(W_{ij}/4\) halfway between the two endmembers.
Asymmetric solid solutions¶
Some solid solutions exhibit asymmetric excess terms. These can be accounted for with an asymmetric solid solution [DPWH07]
\(\alpha\) is a vector of “van Laar parameters” governing asymmetry in the excess properties.
The \(w_{ij}\) terms are a set of interaction terms between endmembers \(i\) and \(j\). If all the \(\alpha\) terms are equal to unity, a nonzero \(w\) yields an excess with a quadratic form and a maximum of \(w/4\) halfway between the two endmembers.
Subregular solid solutions¶
An alternative way to create asymmetric solution models is to expand each binary term as a cubic expression [HW89]. In this case,
Note the similarity with the symmetric solution model, the primary difference being that there are not two interaction terms for each binary.
Thermodynamic and thermoelastic properties¶
From the preceeding equations, we can define the thermodynamic potentials of solid solutions:
We can also define the derivatives of volume with respect to pressure and temperature
Making the approximation that the excess entropy has no temperature dependence
Including orderdisorder¶
Orderdisorder can be treated trivially with solid solutions. The only difference between mixing between ordered and disordered endmembers is that disordered endmembers have a nonzero configurational entropy, which must be accounted for when calculating the excess entropy within a solid solution.
Including spin transitions¶
The regular solid solution formalism should provide an elegant way to model spin transitions in phases such as periclase and bridgmanite. High and low spin iron can be treated as different elements, providing distinct endmembers and an excess configurational entropy. Further excess terms can be added as necessary.
Calculating Multiphase Composite Properties¶
Averaging schemes¶
After the thermoelastic parameters (\(K_S\), \(G\), \(\rho\)) of each phase are determined at each pressure and/or temperature step, these values must be combined to determine the seismic velocity of a multiphase assemblage. We define the volume fraction of the individual minerals in an assemblage:
where \(V_i\) and \(n_i\) are the molar volume and the molar fractions of the \(i\) th individual phase, and \(V\) is the total molar volume of the assemblage:
The density of the multiphase assemblage is then
where \(\rho_i\) is the density and \(\mu_i\) is the molar mass of the \(i\) th phase.
Unlike density and volume, there is no straightforward way to average the bulk and shear moduli of a multiphase rock, as it depends on the specific distribution and orientation of the constituent minerals. BurnMan allows several schemes for averaging the elastic moduli: the Voigt and Reuss bounds, the HashinShtrikman bounds, the VoigtReussHill average, and the HashinShtrikman average [WDOConnell76].
The Voigt average, assuming constant strain across all phases, is defined as
where \(X_i\) is the bulk or shear modulus for the \(i\) th phase. The Reuss average, assuming constant stress across all phases, is defined as
The VoigtReussHill average is the arithmetic mean of Voigt and Reuss bounds:
The HashinShtrikman bounds make an additional assumption that the distribution of the phases is statistically isotropic and are usually much narrower than the Voigt and Reuss bounds [WDOConnell76]. This may be a poor assumption in regions of Earth with high anisotropy, such as the lowermost mantle, however these bounds are more physically motivated than the commonlyused VoigtReussHill average. In most instances, the VoigtReussHill average and the arithmetic mean of the HashinShtrikman bounds are quite similar with the pure arithmetic mean (linear averaging) being well outside of both.
It is worth noting that each of the above bounding methods are derived from mechanical models of a linear elastic composite. It is thus only appropriate to apply them to elastic moduli, and not to other thermoelastic properties, such as wave speeds or density.
Computing seismic velocities¶
Once the moduli for the multiphase assemblage are computed, the compressional (\(P\)), shear (\(S\)) and bulk sound (\(\Phi\)) velocities are then result from the equations:
To correctly compare to observed seismic velocities one needs to correct for the frequency sensitivity of attenuation. Moduli parameters are obtained from experiments that are done at high frequencies (MHzGHz) compared to seismic frequencies (mHzHz). The frequency sensitivity of attenuation causes slightly lower velocities for seismic waves than they would be for high frequency waves. In BurnMan one can correct the calculated acoustic velocity values to those for long period seismic tomography [MA81]:
Similar to [MBR+07], we use a \(\beta\) value of 0.3, which falls in the range of values of \(0.2\) to \(0.4\) proposed for the lower mantle (e.g. [KS90]). The correction is implemented for \(Q\) values of PREM for the lower mantle. As \(Q_S\) is smaller than \(Q_P\), the correction is more significant for S waves. In both cases, though, the correction is minor compared to, for example, uncertainties in the temperature (corrections) and mineral physical parameters. More involved models of relaxation mechanisms can be implemented, but lead to the inclusion of more poorly constrained parameters, [MB07]. While attenuation can be ignored in many applications [TVV01], it might play a significant role in explaining strong variations in seismic velocities in the lowermost mantle [DGD+12].
User input¶
Mineralogical composition¶
A number of predefined minerals are included in the mineral library and users can create their own. The library includes wrapper functions to include a transition from the highspin mineral to the lowspin mineral [LSMM13] or to combine minerals for a given iron number.
Standard minerals – The ‘standard’ mineral format includes a list of parameters given in the above table. Each mineral includes a suggested EoS with which the mineral parameters are derived. For some minerals the parameters for the thermal corrections are not yet measured or calculated, and therefore the corrections can not be applied. An occasional mineral will not have a measured or calculated shear moduli, and therefore can only be used to compute densities and bulk sound velocities. The mineral library is subdivided by citation. BurnMan includes the option to produce a LaTeX; table of the mineral parameters used. BurnMan can be easily setup to incorporate uncertainties for these parameters.
Minerals with a spin transition – A standard mineral for the high spin and low spin must be defined separately. These minerals are “wrapped,” so as to switch from the high spin to high spin mineral at a give pressure. While not realistic, for the sake of simplicity, the spin transitions are considered to be sharp at a given pressure.
Minerals depending on Fe partitioning – The wrapper function can partition iron, for example between ferropericlase, fp, and perovskite, pv. It requires the input of the iron mol fraction with regards to Mg, \(X_\mathrm{fp}\) and \(X_\mathrm{pv}\), which then defines the chemistry of an MgFe solid solution according to (\(\mathrm{Mg}_{1X_{\mathrm{Fe}}^{\mathrm{fp}}},\mathrm{Fe}_{X_{\mathrm{Fe}}^{\mathrm{fp}}})\mathrm{O}\) or \((\mathrm{Mg}_{1X_{\mathrm{Fe}}^{\mathrm{pv}}},\mathrm{Fe}_{X_{\mathrm{Fe}}^{\mathrm{pv}}})\mathrm{SiO_3}\). The iron mol fractions can be set to be constant or varying with P and T as needed. Alternatively one can calculate the iron mol fraction from the distribution coefficient \(K_D\) defined as
We adopt the formalism of [NFR12] choosing a reference distribution coefficient \(K_{D0}\) and standard state volume change (\(\Delta \upsilon^{0}\)) for the FeMg exchange between perovskite and ferropericlase
where \(R\) is the gas constant and \(P_0\) the reference pressure. As a default, we adopt the average \(\Delta \upsilon^{0}\) of [NFR12] of \(2\cdot10^{7}\) \(m^3 mol^{1}\) and suggest using their \({K_D}_0\) value of \(0.5\).
The multiphase mixture of these minerals can be built by the user in three ways:
1. Molar fractions of an arbitrary number of predefined minerals, for example mixing standard minerals mg_perovskite (\(\mathrm{MgSiO_3}\)), fe_perovskite (\(\mathrm{FeSiO_3}\)), periclase (\(\mathrm{MgO}\)) and wüstite (\(\mathrm{FeO}\)).
2. A twophase mixture with constant or (\(P,T\)) varying Fe partitioning using the minerals that include Fedependency, for example mixing \(\mathrm{(Mg,Fe)SiO_3}\) and \(\mathrm{(Mg,Fe)O}\) with a predefined distribution coefficient.
3. Weight percents (wt%) of (Mg, Si, Fe) and distribution coefficient (includes (P,T)dependent Fe partitioning). This calculation assumes that each element is completely oxidized into its corresponding oxide mineral (\(\mathrm{MgO}\), \(\mathrm{FeO}\), \(\mathrm{SiO_2}\)) and then combined to form ironbearing perovskite and ferropericlase taking into account some Fe partition coefficient.
Geotherm¶
Unlike the pressure, the temperature of the lower mantle is relatively unconstrained. As elsewhere, BurnMan provides a number of builtin geotherms, as well as the ability to use userdefined temperaturedepth relationships. A geotherm in BurnMan is an object that returns temperature as a function of pressure. Alternatively, the user could ignore the geothermal and compute elastic velocities for a range of temperatures at any give pressure.
Currently, we include geotherms published by [BS81] and [And82]. Alternatively one can use an adiabatic gradient defined by the thermoelastic properties of a given mineralogical model. For a homogeneous material, the adiabatic temperature profile is given by integrating the ordinary differential equation (ODE)
This equation can be extended to multiphase composite using the first law of thermodynamics to arrive at
where the subscripts correspond to the \(i\) th phase, \(C_P\) is the heat capacity at constant pressure of a phase, and the other symbols are as defined above. Integrating this ODE requires a choice in anchor temperature (\(T_0\)) at the top of the lower mantle (or including this as a parameter in an inversion). As the adiabatic geotherm is dependent on the thermoelastic parameters at high pressures and temperatures, it is dependent on the equation of state used.
Seismic Models¶
BurnMan allows for direct visual and quantitative comparison with seismic velocity models. Various ways of plotting can be found in the examples. Quantitative misfits between two profiles include an L2norm and a chisquared misfit, but user defined norms can be implemented. A seismic model in BurnMan is an object that provides pressure, density, and seismic velocities (\(V_P, V_\Phi, V_S\)) as a function of depth.
To compare to seismically constrained profiles, BurnMan provides the 1D seismic velocity model PREM [DA81]. One can choose to evaluate \(V_P, V_\Phi, V_S, \rho, K_S\) and/or \(G\). The user can input their own seismic profile, an example of which is included using AK135 [KEB95].
Besides standardized 1D radial profiles, one can also compare to regionalized average profiles for the lower mantle. This option accommodates the observation that the lowermost mantle can be clustered into two regions, a ‘slow’ region, which represents the socalled Large Low Shear Velocity Provinces, and ‘fast’ region, the continuous surrounding region where slabs might subduct [LCDR12]. This clustering as well as the averaging of the 1D model occurs over five tomographic S wave velocity models (SAW24B16: [MegninR00]; HMSLS: [HMSL08]; S362ANI: [KED08]; GyPSuM: [SFBG10]; S40RTS: [RDvHW11]). The strongest deviations from PREM occur in the lowermost 1000 km. Using the ‘fast’ and ‘slow’ S wave velocity profiles is therefore most important when interpreting the lowermost mantle. Suggestion of compositional variation between these regions comes from seismology [TRCT05][HW12] as well as geochemistry [DCT12][JCK+10]. Based on thermochemical convection models, [SDG11] also show that averaging profiles in thermal boundary layers may cause problems for seismic interpretation.
We additionally apply cluster analysis to and provide models for P wave velocity based on two tomographic models (MITP08: [LvdH08]; GyPSuM: [SMJM12]). The clustering results correlate well with the fast and slow regions for S wave velocities; this could well be due to the fact that the initial model for the P wave velocity models is scaled from S wave tomographic velocity models. Additionally, the variations in P wave velocities are a lot smaller than for S waves. For this reason using these adapted models is most important when comparing the S wave velocities.
While interpreting lateral variations of seismic velocity in terms of composition and temperature is a major goal [TDRY04][MCD+12], to determine the bulk composition the current challenge appears to be concurrently fitting absolute P and S wave velocities and incorporate the significant uncertainties in mineral physical parameters).