Thermoelastic properties of MgSiO 3 -majorite at high temperatures and pressures: A first principles study

As the major component of garnet, the second most abundant phase in Earth's transition zone, MgSiO 3 -majorite plays a fundamental role in controlling the state and dynamics of Earth's mantle. However, due to challenges of experiments and simulations, there are still very limited data on the elastic properties of MgSiO 3 -majorite at simultaneously high temperatures and pressures. In this study, we have carried out extensive first principles calculations to determine the thermoelastic properties of MgSiO 3 -majorite up to 2000 K and 40 GPa. We find that the elastic constants of MgSiO 3 -majorite change significantly over the temperature and pressure range studied, with noticeable non-linearities in their pressure dependences. The seismic anisotropy of MgSiO 3 -ma- jorite is high and generally increases with pressure. It is much higher than that of the other end-members of garnet and ringwoodite, which makes it the most anisotropic mineral in assemblages expected in the lower transition zone. Based on our calculated elastic moduli and with careful elimination of systematic errors, we establish a third-order Birch-Murnaghan-Mie-Grüneisen model for MgSiO 3 -majorite with the parameters: V 114.1 K 0 ′ = 4.44, G 0 ′ = 1.16, γ 0 = 1.08, q 0 = 0.48, η S 0 = 0.76, and θ 0 = 822.5 K. Integrating our results into a thermodynamic model able to predict the properties of mantle assemblages, we find that a pyrolite composition produces velocities that agree with the seismic model AK135 in the upper transition zone. In the lower transition zone, a pyrolite composition fits well with some specific local observations, but a mechanical mixture with 18% basalt and 82% harzburgite is in better agreement with the global seismic model PREM. The much larger abundance of MgSiO 3 -majorite in the garnet phase of harzburgite suggests that the anisotropy in the lower transition zone may not be negligible and would be observable at least in the heterogeneous zones near subducting slabs.


Introduction
From the constraints provided by seismic observations and mineral physics, garnet is recognized to be the major phase in the Earth's upper mantle and transition zone (Stixrude and Lithgow-Bertelloni, 2012). Accounting for up to over 40 mol% in the normal transition zone and over 90 mol% in subducting slabs (Stixrude and Lithgow-Bertelloni, 2012), it largely controls the thermal and chemical evolution of these areas. To strive for a better understanding of the role played by garnet, it is important to know its thermoelastic properties over a wide range of temperatures and pressures, which are mainly determined by those of its constituent components i.e. pyrope, almandite, grossular, majorite, etc.
As a major end-member of garnet, MgSiO 3 -majorite becomes increasingly important from the base of the upper mantle and is present throughout the transition zone (Ringwood, 1967). However, due to challenges in measuring elastic properties at simultaneously high temperature and pressure, most experimental studies of MgSiO 3 -majorite only report values at ambient conditions (Angel et al., 1989;Bass and Kanzaki, 1990;Gwanmesia et al., 1998;Heinemann et al., 1997;Kato and Kumazawa, 1985;Kavner et al., 2000;Liu et al., 2019;Matsubara et al., 1990;Pacalo and Weidner, 1997;Sawamoto, 1987;Sinogeikin and Bass, 2002a;Sinogeikin et al., 1997). Only limited data have been reported at high-pressure and room temperature (Hazen et al., 1994;Yagi et al., 1992). In addition, there are only a small number of simulation studies reporting the elastic properties of MgSiO 3https://doi.org/10.1016/j.pepi.2020.106491 Received 2 December 2019; Received in revised form 27 March 2020; Accepted 8 April 2020 majorite, due to the significant computational expense associated with modeling the relatively large unit cell. Most studies only report values at static conditions (De La Pierre and Belmonte, 2016;Hernandez et al., 2015;Li et al., 2007;Pigott et al., 2015;Vinograd et al., 2006;Yu et al., 2011). One study reports data for the bulk modulus of MgSiO 3 -majorite at simultaneous high temperature and pressure but provides no values for the shear modulus or elastic constants (Yu et al., 2011). There are even fewer studies that report values for seismic anisotropy (Li et al., 2007;Pacalo and Weidner, 1997;Pigott et al., 2015).
In this study we have performed first principles calculations to determine the thermoelastic properties of MgSiO 3 -majorite over a wide range of temperatures and pressures. We determine elastic constants, anisotropies, bulk and shear moduli, seismic velocities and other thermodynamic properties. After carefully eliminating systematic errors in the simulation results, we obtain a third-order Birch-Murnaghan-Mie-Grüneisen model for MgSiO 3 -majorite. Using this model to predict the properties of mineral assemblages in the transition zone, we find that a mechanical mixture of basalt and harzburgite may be a better interpretation than the equilibrium pyrolitic model for global seismic observations in the lower transition zone.

Crystal structure
MgSiO 3 -majorite, Mg 3 (Mg,Si)[SiO 4 ] 3 , is a garnet-structured mineral with tetragonal space group I4 1 /a. The crystal structure of MgSiO 3majorite includes three different types of polyhedra: tetrahedra, octahedra, and dodecahedra. Tetrahedral sites are occupied by Si 4+ , octahedral sites are occupied by Mg 2+ and Si 4+ and dodecahedral sites are occupied by Mg 2+ (Fig. 1). There are two distinct octahedral sites with Wyckoff positions of 8c and 8d. In the ordered I4 1 /a structure, all 8c sites are occupied by Si 4+ cations and all 8d sites by Mg 2+ . At sufficiently high temperatures, a random distribution of Si 4+ and Mg 2+ in these two sites is possible, leading to disordered structures with higher symmetry (Angel et al., 1989;Li et al., 2007).
The starting structure of MgSiO 3 -majorite used in the simulations was the experimentally determined structure (Angel et al., 1989), which has also been used in previous works (Pigott et al., 2015;Vinograd et al., 2006;Yu et al., 2011). It is a tetragonal unit cell with 160 atoms (8 Mg 3 (Mg,Si)[SiO 4 ] 3 units). According to the study of Yu et al. (2011), the order-disorder transition between magnesium and silicon in the octahedral sites occurs at about 3600 K and the long-range order parameter, which is 1.0 for a completely ordered state and 0.0 for a completely disordered state, is always larger than 0.97 when the temperature is lower than 2500 K. Therefore, in line with Yu et al. (2011), we used the lowest energy ordered I4 1 /a structure (as shown in Fig. 1) throughout this study. Simulations of some other cation-ordering arrangements yield comparable results to those described here, as shown in Table S1.

First principles simulations
We determined the effects of pressure and temperature on the structure and elasticity of majorite using first principles simulations. This approach couples density functional theory (DFT; Hohenberg and Kohn, 1964;Kohn and Sham, 1965), used to evaluate the forces between atoms arising from the electronic structure, with Newtonian mechanics, used to simulate the thermal motion of atoms. Simulations were performed using the Vienna ab initio simulation package (VASP; Furthmüller, 1996a, 1996b;Kresse and Hafner, 1993;Kresse and Joubert, 1999), which makes use of a plane-wave expansion of the Kohn-Sham wavefunctions (Payne et al., 1992) and the projector augmented wave (PAW) approach for the description of core electrons (Blöchl, 1994;Kresse and Joubert, 1999). The PBE (Perdew et al., 1996) implementation of GGA (generalized gradient approximation) to the exchange-correlation functional was used throughout. The core radii of the PAW pseudopotentials used were 1.06 Å for Mg (valence configuration with p semi-core valence state is 2p 6 3s 2 ), 0.79 Å for Si (3s 2 3p 2 ) and 0.80 Å for O (2s 2 2p 4 ).
For the static calculations, we used a 2 × 2 × 2 Brillouin zone sampling grid, an energy for the plane-wave cutoff of 800 eV and an energy convergence criterion of 10 −8 eV. Convergence tests indicated that total energies calculated with these parameters were converged to within 0.2 meV/atom. We optimized the structure by minimizing the enthalpy at 0, 10, 20, 30 GPa and then calculated the elastic constants, described below. For the high-temperature simulations, Brillouin zone sampling was restricted to the zone center (Γ point), a plane-wave cutoff of 500 eV was used and the energy convergence criterion was 10 −6 eV. Convergence tests indicated that total energies calculated with these parameters were converged to better than 1.0 meV/atom. The initial configurations for the high-temperature calculations were the optimized structures of MgSiO 3 -majorite from the static simulations. We performed ab initio molecular dynamics (AIMD) simulations within the NVT ensemble, making use of the Nosé thermostat (Nose, 1984). Calculations were run at temperatures of 1000, 1500 and 2000 K and each run for a total of 6 ps, using a time-step of 1 fs. We carefully verified the stress tensors of these undeformed AIMD simulations were hydrostatic, as shown in Fig. S1.

Calculations of elastic properties
To calculate the elastic constants, we follow the method outlined in previous studies (Karki et al., 2001;Oganov et al., 2001;Wallace, 1972;Zhang et al., 2013). Starting from the geometry output from a static relaxation simulation or an equilibrium undeformed molecular dynamics simulation, we apply 4 types of strain: 3 axial (along each crystallographic axis) and 1 triclinic (shear strain) and calculate the resulting stresses. The magnitudes of these strains were ± 0.1%, ± 0.3% for static calculations and ± 1%, ± 2% for the high-temperature simulations. The isothermal elastic constants, C ijkl T , were then calculated by fitting a second-order polynomial to the strain-stress data. The isothermal elastic constants were converted to adiabatic elastic constants, C ijkl S , using standard thermodynamically self-consistent formulae (Davies, 1974;Wallace, 1998): where C V is the isochoric heat capacity and b ij = (∂σ ij /∂T) V is the thermal stress tensor, which were calculated by linear-least-square fitting of internal energy and stress tensor along isochores. Based on these results, we obtained Grüneisen parameters, thermal expansivity, Voigt-Reuss-Hill bulk and shear moduli, and compressional and shear wave velocities. Uncertainties in these quantities were determined with appropriate non-Gaussian statistics via the blocking method (Flyvbjerg and Petersen, 1989). Seismic anisotropies were calculated by using MSAT (Walker and Wookey, 2012) and azimuthal anisotropy for V P and respectively (Kawai and Tsuchiya, 2015).

Lattice parameters and elastic constants
In Table 1, we compare the lattice parameters from this study and those from previous theoretical and experimental studies. Compared with the experiments, our calculated values overestimate the lattice parameters (a and c) by~1% and volume by~4%, which is expected since GGA functionals tend to under-bind and overestimate pressure and moduli. These systematic errors were also found in previous DFT simulations using GGA, such as those from Li et al. (2007), which perfectly match our static calculations over the whole pressure range from 0 to 30 GPa. Compared with other simulations, the maximum differences are~1% in a and c, and~4% in volume, due to differences in the simulation details.

Table 1
Comparisons of the lattice parameters of MgSiO 3 -majorite. Angel et al. (1989) 11.501(1) 11.480(2) 0.998 (2) 1518.6(4) Liu et al. (2017) 11.509 (21) 11.446(43) 0.995 1516.1(60) Sinogeikin et al. (1997) 11.500 (35) 11.445 (50) 11.445 (50) 1513.6(15) Pacalo and Weidner (1997) 11.514 (1)  The raw simulated elastic constants are listed in Table S2. In Fig. 2, we show these data as a function of volume and pressure at various temperatures and compare them with currently available experimental and theoretical data. At static conditions, the results of Li et al. (2007) are in perfect agreement with our own. The elastic constants determined by Pigott et al. (2015) through static lattice-energy calculations (SLEC) generally agree well with ours, but C 23 and C 12 are about 8.8% and~9.5% larger than our predictions. The experimental data of Pacalo and Weidner (1997) are also in accord with our results except for C 12 and C 16 , which strongly deviate from the predictions of all the theoretical simulations. Compared with the results of Pacalo and Weidner (1997), C 12 calculated in this study is~40 GPa larger and C 16 is~16 GPa larger (Table S3).
The temperature and pressure dependence of the elastic constants are shown in Fig. 2(b). C 11 and C 33 increase significantly with decreasing temperature or increasing pressure. They show some nonlinear pressure dependence, especially in the static calculations. This is similar to the pressure dependence of the diagonal component elastic constants in pyrope, as mentioned by Hu et al. (2016). C 12 and C 23 also increase strongly with the pressure, but in a more linear manner. The other elastic constants, C 44 , C 66, and C 16 , vary less with temperature and pressure, although the non-linearity is more evident than C 12 and C 23 .

Seismic anisotropy
The seismic anisotropies of MgSiO 3 -majorite, based on our elastic constants and those from other studies, are shown in Fig. 3. In all cases, the symmetry of the seismic anisotropy is not aligned with the lattice vectors. This is in contrast with crystals with cubic or orthorhombic structured symmetry (Mainprice, 2015), and caused by the finite value in C 16 . The physical reason for this is schematically illustrated in Fig.  S2, where in contrast to cubic pyrope, strain in a-direction (ε 1 ) in the ordered tetragonal structure causes non-zero shear stress (σ 6 ), because the tetrahedra are connected by non-equivalent adjacent octahedra: one connects with four softer Mg octahedra and the other with four harder Si octahedra.
Our calculated seismic anisotropy, at static conditions, is in very good agreement with that of Li et al. (2007) and similar to that of Pigott et al. (2015), but does not agree with the results of Pacalo and Weidner (1997). The orientation of the pattern of anisotropy is the most prominent difference between theoretical studies and the experimental work of Pacalo and Weidner (1997). The C 16 of Pacalo and Weidner (1997) is only 1.4 GPa, which is much smaller than all the theoretical studies. As mentioned in Li et al. (2007), the MgSiO 3 -majorite sample of Pacalo and Weidner (1997) may be somewhat disordered and, although the degree of disorder is unknown, this would plausibly explain the deviation here because disorder would make the two sites more similar on average.
In Fig. 3, the anisotropies from our simulations and those from Li et al. (2007) are almost the same in both A P (8.0%:7.4%) and A Smax (maximum value of A S , 14.70%:13.07%). Compared with those of Pigott et al. (2015), there are some minor differences in both A P (8.0%:5.1%) and A Smax (14.70%:10.44%). Fig. 3(a) and (d) show that our calculated anisotropies at 0 K and 0 GPa are much larger than those measured by Pacalo and Weidner (1997) at ambient conditions. For A P it is about 4 times larger (8.0%:1.8%) and for A Smax it is about 2 times larger (14.70%: 9.06%). In Fig. 4, we show the effects of pressure and temperature on the anisotropy of MgSiO 3 -majorite. At static conditions, A Smax increases, but A P decreases with increasing pressure. At hightemperature conditions (1000, 1500, 2000 K), A P and A Smax vary slightly with temperature and both increase with pressure. At high temperatures and pressures, the values of A P and A Smax become almost invariant, at about 6.0% and 15.0% respectively. Fig. 4 also shows the anisotropies of pyrope (Py), grossular (Gro) and ringwoodite (Rwd) determined in previous studies. They are obviously smaller, making MgSiO 3 -majorite the most anisotropic mineral in assemblages expected in the lower transition zone.

Elastic moduli and velocities
In Fig. 5, we plot calculated pressure and elastic moduli, as a function of volume, and compare them with previous studies. In order to make more precise comparisons with experimental data at ambient conditions, we eliminate the vibrational contribution in these experimental data following the approach proposed by Panero and Stixrude (2004). Compare to these "static" experimental data, the systematic errors from the PBE exchange-correlation functional are found to be up to~10 GPa in pressure,~20 GPa in isothermal bulk modulus and 20 GPa in shear modulus. We applied a rescaling method to eliminate these systematic errors (Zhang and Liu, 2015;Zhang et al., 2018;Zhang et al., 2013) Fig. 3. Seismic anisotropy of MgSiO 3 -majorite at 0 K and 0 GPa calculated using our elastic constants, and those of Pacalo and Weidner (1997), Pigott et al. (2015) and Li et al. (2007). The left legend shows the value range of V P splitting (the leftmost column) and the right legend is for dV S splitting (two columns on the righthand side). The black dots represent maximum and minimum values in the anisotropy spheres.
where M is the property to be rescaled (P, K or G), subscript 0 refers to properties at zero pressure and static conditions, quantities lacking temperature dependence are athermal (static), K is the isothermal bulk modulus, and Δ denotes the change in the quantity from static conditions to the temperature of interest, i. GPa. We report our final corrected data in Table 2. Based on these data, we refitted the third-order Birch-Murnaghan-Mie-Grüneisen model proposed by Lithgow-Bertelloni (2005, 2011). The new parameters for the thermoelastic model of MgSiO 3 -majorite are listed in the caption of Table 2. After eliminating systematic errors, our results are in good agreement with the experimental data within their uncertainties or the scattering of different datasets. Compared with previous simulations, the volumes of Yu et al. (2011) at 300 K are a little higher than our results (Fig. 5(a)) due to small thermal effects. The results from the Lines that connect data points are simple straight lines that are used to show the trend. Solid colored squares are our data for MgSiO 3 -majorite, dark circles and triangles represent theoretical data of pyrope and grossular at static conditions, respectively. Open diamonds represent experimental data for ringwoodite at ambient conditions. Open triangles are experimental data for ringwoodite at 1793 K, which were obtained by extrapolating experimental values (Jackson et al., 2000;Sinogeikin et al., 1998;Sinogeikin et al., 2001;Weidner et al., 1984) to a depth of 550 km. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) SLEC of Vinograd et al. (2006) and Pigott et al. (2015) are consistent with our rescaled data but at high pressures, deviations are noticeable, up to 13.5 GPa (~5.6%) in bulk modulus and 6.7 GPa (~6.3%) in shear modulus. At high temperatures, we find apparent good agreement between our rescaled results and those of Yu et al. (2011) over the whole pressure range, but more careful analysis (Fig. S3) reveals some subtle deviations. These may be due the neglect of anharmonicity in the lattice dynamics simulations of Yu et al. (2011), which is included in our molecular dynamics simulations. Oganov et al. (2000) showed that this can lead to overestimation of the thermal expansion coefficients and Table 2 Elastic moduli of MgSiO 3 -majorite from static conditions to high temperatures and pressures, rescaled from the raw simulation results in this study (listed in Table  S2). Parameters for the thermodynamic model of Stixrude and Lithgow-Bertelloni (2011) refitted from these results are: T 0 = 300 K, V 0 = 114.1 cm 3 /mol, K 0 = 163.6 GPa, G 0 = 86.4 GPa, K 0 ′ = 4.44, G 0 ′ = 1.16, γ 0 = 1.08, q 0 = 0.48, η S0 = 0.76, θ 0 = 822.5 K. V means volumes of unit cell. The uncertainties are within 0.1 GPa for pressure and generally within 3% for bulk and shear moduli.  (Babuska et al., 1978;Bass and Kanzaki, 1990;Chen et al., 1999;Conrad et al., 1999;Gwanmesia et al., 1998;Gwanmesia et al., 2000;Gwanmesia et al., 2009;Gwanmesia et al., 2006;Hazen and Finger, 1989;Isaak and Graham, 1976;Kavner et al., 2000;Leger et al., 1990;Leitner et al., 1980;Levien et al., 1978;Liu et al., 2000;Liu et al., 2019;Liu et al., 2015;Morishima et al., 1999;Pacalo and Weidner, 1997;Pamato et al., 2016;Rigden et al., 1994;Sato et al., 1978;Sinogeikin and Bass, 2000, 2002a, 2002bSinogeikin et al., 1997;Wang et al., 1998;Wang and Ji, 2001;Yagi et al., 1987;Yeganehhaeri et al., 1990;Zhang et al., 1999;Zou et al., 2012). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Y. Lou, et al. Physics of the Earth and Planetary Interiors 303 (2020) 106491 Grüneisen parameters at high temperatures, which we also find ( Fig.  S4). At a fixed volume, overestimation of the thermal expansion coefficient, results in an overstimation of pressure, leading to difference in calculated PeV curves and, in turn bulk moduli (up to over 30 GPa for isothermal bulk modulus and over 20 GPa for the adiabatic bulk modulus), as shown in Fig. S3. In Fig. 6, we show the adiabatic moduli of MgSiO 3 -majorite at various pressures and temperatures. The bulk modulus significantly increases with pressure, doubling when going from 0 to 30 GPa. The shear modulus also increases noticeably by up to 40 GPa over the same pressure range. In common with the elastic constants, both moduli show some non-linear pressure-dependence. The bulk and shear modulus exhibit strong temperature dependence. Increasing the temperature by 2000 K, causes the adiabatic bulk and shear moduli to decrease linearly by 22-26 GPa and 16-18 GPa, respectively.
As an important solid solution in the transition zone, the majoritepyrope join has been the focus of many experimental studies (Bass and Kanzaki, 1990;Gwanmesia et al., 2000;Gwanmesia et al., 2009;Hunt et al., 2010;Liu et al., 2000;Liu et al., 2019;Liu et al., 2015;Morishima et al., 1999;Pamato et al., 2016;Rigden et al., 1994;Bass, 2002a, 2002b;Sinogeikin et al., 1997;Wang et al., 1998;Yagi et al., 1987;Yeganehhaeri et al., 1990). We plot the bulk and shear moduli from these experimental studies against the composition in Fig. 7, together with those calculated from our results for MgSiO 3majorite and those from the predictions of Stixrude and Lithgow-Bertelloni (2011) for Mg 3 Al 2 [SiO 4 ] 3 -pyrope. Assuming linear mixing, as suggested by previous studies (Bass and Kanzaki, 1990;Gwanmesia et al., 2000;Liu et al., 2019), we find the linear equations of K s = 171.2 − 0.065 * X Mj and G = 93.7 − 0.073 * X Mj , determined by our calculations, agree with the experiments within the uncertainties and scattering of the experimental data. We calculate K S , G, V P and V S for majorite-pyrope solid solutions at high temperatures and pressures (Fig. 8), which show good agreement with the measured values of Gwanmesia et al. (2009 and Pamato et al. (2016).

Implications for interpreting the seismic observations in the transition zone
The thermoelastic model of MgSiO 3 -majorite obtained in this study has been used to determine the seismic velocities of pyrolite, basalt, and harzburgite mineral assemblages, with the bulk chemical compositions reported by Workman and Hart (2005), Baker and Beckett (1999) and compiled by Xu et al. (2008), as listed in Table S4. Besides MgSiO 3majorite, we included another 46 end-members and 21 phases, using the parameters listed in Table A1 of Stixrude and Lithgow-Bertelloni (2011), except for MgSiO 3 -perovskite (bridgmanite) for which we used the parameters optimized by Zhang et al. (2013) and cubic CaSiO 3perovskite for which we used the parameters recently re-determined by Thomson et al. (2019). Note that for the garnet phase, in which majorite resides, we have included five typical end-members (pyrope, almandine, grossular, Mg-Majorite and Jd-Majorite, as listed in Table A1 of Stixrude and Lithgow-Bertelloni (2011)). So the effect of iron in the garnet solid solution has only been accounted for on the dodecahedral sites. The adiabatic geotherm with a potential temperature of 1500 K leads to a temperature range of~1690-1820 K in the transition zone, which is within the estimate of 1773 ± 100 K at a depth of 410 km suggested by Frost (2008). With these inputs, we obtained the phase proportions of pyrolite, basalt and harzburgite by minimizing the Gibbs free energy, as proposed by Stixrude and Lithgow-Bertelloni (2011). Using the properties of individual phases, we calculated the seismic velocities of equilibrium mineral assemblages through Voigt-Reuss-Hill averaging.
In Fig. 9, we compare the velocities predicted for a pyrolite mineral assemblage with several seismic models. It shows that the agreement is fair to good for the global seismic models (Dziewonski and Anderson, 1981;Kennett et al., 1995) in the upper transition zone, considering the uncertainties in the seismic observations and the possible uncertainty in temperature (Xu et al., 2008). The agreement extends to the lower transition zone, as compared with the local seismic model of TNA2, which was proposed by Wang et al., 2014 for the transition zone beneath the East Pacific Rise.
Compared with PREM and AK135, the equilibrium pyrolite model seems to predict too low shear-wave velocities with increasing depth from about 520 km. Due to the sluggish chemical diffusivity of mantle materials (10 −14 -10 −16 m 2 s −1 ) in the solid state (Farber et al., 1994;Hofmann and Hart, 1978;Yamazaki et al., 2000), mechanically mixed mantle models are considered by a number of previous studies (Allègre and Turcotte, 1986;Hofmann and Hart, 1978;Ritsema et al., 2009;Xu et al., 2008). Inspired by these studies, we separately calculated the densities and seismic velocities of basalt and harzburgite compositions, and linearly combined them with 18% basalt and 82% harzburgite to get the densities and seismic velocities of the modeled mechanical mixture proposed by Xu et al. (2008). The density of the mechanical mixture is, just like pyrolite, in good agreement with PREM and AK135, but its velocities are found to be in better agreement with PREM than equilibrium pyrolite, as shown in Fig. 9.
Similar to the comparisons in Fig. 9, a number of previous studies support the aggregation of harzburgite in the lower transition zone (Irifune et al., 2008;Nakagawa and Buffett, 2005). This would be important for the anisotropy of the lower transition zone. As shown in Fig.  S5, the amount of MgSiO 3 -majorite in the garnet phase of harzburgite is significantly higher than that in pyrolite and basalt, and according to Sinogeikin et al. (1997), the solid solution of majorite-pyrope maintains the tetragonal structure when MgSiO 3 -majorite accounts for 70% or higher. The noticeable anisotropy, as revealed in Figs. 3 and 4, would, therefore, contribute to the anisotropy in the lower transition zone, which may be very small on a global scale (Huang et al., 2019), but observable, at a local scale, in heterogeneous zones near the subducting slabs (Lynner and Long, 2015).

Conclusions
We have performed first principles calculations to determine the thermoelastic properties of MgSiO 3 -majorite over a wide range of temperatures and pressures. The elastic constants for MgSiO 3 -majorite show noticeable non-linear pressure dependence. The anisotropies in V P and V S of MgSiO 3 -majorite are found to be higher than the other major minerals such as pyrope, grossular, and ringwoodite, making MgSiO 3majorite the most anisotropic mineral in assemblages expected in the lower transition zone. With careful elimination of the systematic errors in the simulation results, the determined elastic properties are shown to be in good agreement with the limited available experimental measurements for MgSiO 3 -majorite and the majorite-pyrope join. Based on the high quality data obtained in this study, we have proposed an optimized thermodynamic model for MgSiO 3 -majorite that is able to accurately predict its elastic properties at simultaneously high temperatures and pressures.
The thermoelastic properties of MgSiO 3 -majorite determined in this study were used to predict the seismic velocities of pyrolite, harzburgite, and basalt mineral assemblages, at the conditions of the transition zone. These showed that, on a global scale, the seismic velocities of the upper transition zone are consistent with a pyrolite-like composition, but this equilibrium assemblage may only extend to the lower transition zone in local regions. Compared with global seismic models, better agreement with PREM than AK135 implies that the lower transition zone is more consistent with an 18% basalt and 82% harzburgite mechanical mixture. The aggregation of harzburgite in the lower transition zone would lead to a much higher MgSiO 3 -majorite content in the garnet phase and cause a cubic-tetragonal phase transition in garnet, which may result in non-negligible anisotropy in the lower transition zone.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Fig. 9. Comparisons of seismic wave velocities and densities among our models (colored solid lines), global models: PREM (dark grey solid line) and AK135 (dark grey dashed line), and local model TNA2 (dark grey dot line). Our models are along the adiabatic geotherm with a potential temperature of 1500 K. Colored shaded zone are Voigt-Reuss bounds of our models. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Y. Lou, et al. Physics of the Earth and Planetary Interiors 303 (2020) 106491 China (#41590620). Simulations were carried out on the computational facilities in the Computer Simulation Lab of IGGCAS and Tianhe-2 at the National Supercomputer Center of China (NSCC) in Guangzhou.
The collaboration leading to this paper was made possible by the U.K.-Sino Centre for Earth and Planetary Science. AMW was partially supported by the Natural Environment Research Council (NE/M000044/ 1).

Data availability
The data for this paper can be found in the supplementary and more extended raw data are available on request.