Equation of state and sound wave velocities of fayalite at high pressures and temperatures: implications for the seismic properties of the martian mantle

The elastic properties of a pure, synthetic fayalite aggregate were studied by coupled synchrotron X-ray diffraction and ultrasonic interferometry in a DIA-type multi-anvil press. Measurements at pressures up to about 7 GPa and temperatures up to 873 K yielded an adiabatic bulk modulus, $K_{S\mathrm{0}}=\mathrm{127.2}\pm \mathrm{0.3}$ GPa with $(\partial K_S/\partial P{)}_{T\mathrm{0}}=\mathrm{6.5}\pm \mathrm{0.1}$, and a shear modulus, $G_{\mathrm{0}}=\mathrm{53.3}\pm \mathrm{0.4}$ GPa with $(\partial G/\partial P{)}_{T\mathrm{0}}=\mathrm{1.25}\pm \mathrm{0.05}$. When fixing $(\partial K_S/\partial P{)}_{T\mathrm{0}}=\mathrm{5.3}$ (after $(\partial K_T/\partial P{)}_{T\mathrm{0}}$ from Nestola et al., 2011), K_S0 increases to about 130 GPa. These estimates of ($K_{S\mathrm{0}},(\partial K_S/\partial P{)}_{T\mathrm{0}}$) follow a general linear trend, $K=f(\mathrm{d}K/\mathrm{d}P)$, for fayalite. We define limited ranges for both bulk and shear moduli from previous studies, and we discuss how these variations affect seismic velocities and the determination of a mineralogical model in the context of the Mars InSight SEIS (Seismic Experiment for Interior Structure) experiment.

1 Introduction

The successful landing of the Mars InSight mission and more specifically the installation of the seismometer SEIS (Seismic Experiment for Interior Structure) (Lognonné et al., 2019) on the surface of Mars will bring new constraints on the deep interior of the red planet. A depth profile of the seismic velocities is one of the major goals of this mission because (1) it will provide the first direct picture of the inner structure of the planet, in particular the size of its core (Witze, 2021; Stähler et al., 2021), the presence or not of a lower mantle, and possible discontinuities in its upper mantle (a transition zone?); and (2) it will allow, for the first time, mineralogical models to be fitted to real martian data and therefore the mineralogy of Mars' interior to be constrained.

Building such models requires knowing the thermoelastic properties of candidate minerals, i.e., minerals with compositions close to what is expected for Mars. For a solid solution like (Fe,Mg) olivines, unknown elastic parameters of a given composition are calculated from those of pure endmembers according to their relative amounts (e.g., Eq. 6 in Stixrude and Lithgow-Bertelloni, 2005). As recently proposed, such a relation may not work for all olivine compositions (Béjina et al., 2019; Pamato et al., 2019). Consequently it is important to measure the elastic properties across the entire solid solution in order to clearly define a compositional trend. With this goal in mind we started performing experimental measurements of equations of state (EoS) and sound velocities on silicate minerals with an iron content higher than on Earth. In the text below we will use a simplified notation for the expressions of bulk and shear moduli. Unless mentioned otherwise the elastic parameters at constant temperature and room pressure, K_T0, K_S0, G₀ and their derivatives, $(\partial K_T/\partial P{)}_{T\mathrm{0}}$, $(\partial K_S/\partial P{)}_{T\mathrm{0}}$ and $(\partial G/\partial P{)}_{T\mathrm{0}}$, will be denoted as K_T, K_S, G, $\mathrm{d}{K}_T/\mathrm{d}P$ (or $K_T^{\prime }$), $\mathrm{d}{K}_S/\mathrm{d}P$ (or $K_S^{\prime }$) and $\mathrm{d}G/\mathrm{d}P$ respectively.

The iron content of the martian mantle is estimated to be around 17 wt % (e.g., Dreibus and Wänke, 1985; Taylor, 2013); therefore, olivine, one of its major phases, is thought to have a molar Fe content (X_Fe) around 25 %–30 %. Very few values for bulk modulus (K) and shear modulus (G) exist for olivines with X_Fe>0.2. If the trend of G values with X_Fe seems well defined and decreases linearly with increasing X_Fe (e.g., Liu et al., 2005), the few results on K appear contradictory. Nestola et al. (2011) performed room-temperature static compression experiments and measured an isothermal bulk modulus, K_T=124.7 GPa and $K_T^{\prime }=\mathrm{5.3}$ for olivines with $X_{\mathrm{Fe}}=\mathrm{8}\mathit{\%},\mathrm{20}\mathit{\%},\mathrm{29}\mathit{\%},\mathrm{38}$ %, concluding that Fe has little effect on both K and K^′ in this composition range. By contrast Chung (1971) determined a decreasing linear trend of the adiabatic bulk modulus (K_S) between forsterite (K_S=128.1 GPa) and fayalite (K_S=122.0 GPa) from sound velocity measurements. In addition data on pure fayalite are largely discrepant, with K_T (see Fig. S4 and Table S6 in Béjina et al., 2019) ranging from about 104 GPa (Plymate and Stout, 1990, at 400 ^∘C) to 138 GPa (Smyth, 1975) and K_S values of 138 GPa (Sumino, 1979), 128 GPa (Graham et al., 1988) and 134 GPa (Isaak et al., 1993). First-principles calculations by Stackhouse et al. (2010) predicted an even larger bulk modulus for fayalite, above 139 GPa. Therefore, a Fe-compositional trend of the bulk modulus for Fe-rich olivines is difficult to define and may even be nonlinear (Béjina et al., 2019; Pamato et al., 2019).

This lack of constraints on the effect of iron on the elastic properties of the Fe–Mg olivine solid solution, in particular for X_Fe>0.2, makes it difficult to adjust a mineralogical model to SEIS data and decipher the Fe content of the Mars mantle. We recently published high-pressure room-temperature X-ray diffraction (XRD) data (Béjina et al., 2019) and calculated the isothermal bulk modulus (K_T) for the following olivine compositions: $({\mathrm{Fe}}_x,{\mathrm{Mg}}_{\mathrm{1}-x}{)}_{\mathrm{2}}\mathrm{Si}{\mathrm{O}}_{\mathrm{4}}$ with $x=\mathrm{0.45},\mathrm{0.64},\mathrm{0.82}$ and 1. Because room-temperature measurements at high pressure are subjected to consequent differential stress (t), these were estimated by using the method of Singh et al. (1998) that allows for the determination of d spacing values under hydrostatic pressure and the calculation of t. Subsequently the NaCl pressures were corrected for t. We present here the first results from an experiment at high temperatures and pressures above 1 GPa on a polycrystalline fayalite with coupled XRD and ultrasonic interferometry techniques (with the exception of Liu et al., 2010, but they did not calculate the elastic parameters). XRD measurements show that, as expected, high-temperature data present less differential stress than at room temperature, but t can be either positive or negative. The same correction for stress is therefore applied here.

2 Experimental method

The main details of our experiment, performed at the superconducting wiggler beamline X17B2 at the National Synchrotron Light Source (Brookhaven National Laboratory, Uptown, NY), as well as the preparation and characteristics of our fayalite (Fa) aggregate, can be found in our previous paper (Béjina et al., 2019). The acoustic velocity measurements were obtained in conjunction with X-ray diffraction (XRD) in situ at high pressure and temperature with the setup described by Whitaker et al. (2017). We will therefore only give here a brief description of the experiment.

2.1 Sample characteristics

Our Fe₂SiO₄ sample was synthesized from nanometric powders of Fe₂O₃ and SiO₂. The oxide powder mixture was annealed twice for 2 h at 1000 ^∘C and fO₂ buffered by a gas mixture $\mathrm{Ar}/{\mathrm{H}}_{\mathrm{2}}\mathrm{O}+\mathrm{Ar}/{\mathrm{H}}_{\mathrm{2}}$ close to the $\mathrm{Fe}/\mathrm{Fe}\mathrm{O}$ equilibrium. After checking that the run product was pure fayalite by XRD, a dense pellet was prepared by spark plasma sintering (Béjina and Bystricky, 2009; Guignard et al., 2011). Observations of a fragment of this aggregate by scanning electron microscopy (SEM) using energy dispersive X-ray spectroscopy, orientation contrast and electron backscatter diffraction (EBSD) showed a very weak preferred crystallographic orientation and a homogeneous grain size around 2 µm. Finally our 1.2 mm diameter specimen was cored from this pellet, its length was brought to 1.2 mm, and both ends were carefully polished with Al₂O₃ powders.

2.2 High-pressure and high-temperature experiment

Beamline X17B2 at National Synchrotron Light Source (NSLS) was equipped with a 1000 ton DIA-type multi-anvil apparatus placed in the path of a white X-ray beam and with the DIASCoPE acoustic system for sound velocity measurements (Whitaker et al., 2017). Energy-dispersive XRD patterns (see as an example Fig. 3 in Béjina et al., 2019) were collected at a scattering angle of about 6.5^∘ on the detection system composed of 10 detectors placed in a circle (nine between 0 and 180^∘ azimuth plus one at 270^∘) allowing for the calculation of differential stress (see description in Weidner et al., 2010). We later used the stress values to correct the pressure obtained from NaCl XRD analysis. Our experiment was run using the cell assembly shown in Fig. 1 containing both the fayalite sample and NaCl as a pressure calibrant. A detailed description of the cell is given by Whitaker et al. (2017), the main characteristics being a graphite heater with a WRe3%-WRe25% thermocouple to measure temperature (not corrected for pressure effects); a mixture 10NaCl:1BN at the back of the sample serving as both pressure marker using the NaCl EoS by Brown (1999) and as cushioning to reduce differential stress in the sample; and a BN sleeve along the sample helping to provide a pseudo-hydrostatic environment but also creating low-fO₂ conditions needed to keep fayalite stable. In order to maximize the coupling between the different parts of the cell assembly for better sound wave propagation, 1 µm thick gold foils are placed at each interface along the wave path (Fig. 1). The ones around the sample also serve as markers when measuring sample length on X-radiographic images (see Sect. 2.3). The beamline setup allows us to take images of the full sample at once, these being recorded on a Prosilica charge-coupled device (CCD) camera. A dual-mode 10^∘ Y-cut LiNbO₃ piezoelectric transducer generated both P and S waves at resonance frequencies of 62 and 35 MHz respectively. XRD on NaCl and fayalite, P and S wave signals, and sample images were recorded sequentially at each step of the experiment.

Figure 1 (a) Cell assembly designed for ultrasonic interferometry measurements at high pressure and temperature (for details see Whitaker et al., 2017). The X-ray beam enters the cell via an edge of the boron-epoxy cube at 90^∘ to the thermocouple. The dual-mode 10^∘ Y-cut LiNbO₃ piezoelectric transducer, placed at the back of one of the anvils (grey), generates both P and S waves. Al₂O₃ serves as a buffer rod. To maximize coupling between the different parts on the wave paths, 1 µ thick gold foils are placed at each interface. They also serve as markers to measure the sample length on X-ray images. (b) Examples of P and S wave echoes of the different interfaces: anvil–Al₂O₃, Al₂O₃–fayalite and fayalite–NaCl. The travel time across the sample for each wave is calculated by cross-correlation of the echoes reflected from both ends.

The tonnage-temperature path is presented in Fig. S1 in the Supplement. The Fa aggregate was first cold-pressed to 35 tons and then heated at 1073 K for about 15 min in order to relax the accumulated stress. Because the furnace started to become unstable, the temperature was brought down to 873 K, and this remained the maximum temperature of the experiment. The maximum tonnage-temperature conditions were 50 tons and 873 K. XRD on Fa and NaCl, P and S wave data, and X-ray radiographic images of the sample were collected while decreasing temperature (usually by steps of 100 K). Once at room temperature tonnage was lowered, and temperature increased back to 873 K before collecting a new set of data. Similar steps were repeated until tonnage reached 5 tons. Note that at 5 and 7.5 tons, the maximum temperatures were 573 and 673 K respectively. As discussed in Sect. 3, the maximum pressure estimated from NaCl unit-cell volumes and the EoS of Brown (1999) and corrected for differential stress was about 7 GPa at the highest temperature and about 5 GPa at room temperature because of the relaxation of cell-assembly materials.

2.3 Measurements of sample length and (V_p, V_s)

The length of the Fa sample was determined from X-ray radiographic images taken during the experiment thanks to gold foils placed at each end of our sample that blocked X-rays and acted as markers. The length variation is therefore known as a function of pressure and temperature (Fig. S2). The uncertainty in length is around 2 pixels (e.g., Li et al., 2004) corresponding to about 5 µm.

The round-trip travel times of P and S waves were obtained by cross-correlation of the echoes reflected from the top and bottom of the Fa specimen. This gives a precise superposition of the two echoes, and their cross-correlation presents a unique maximum corresponding to the round-trip travel time across the sample (Fig. 2). To perform these calculations we built a Python script largely based on the ObsPY package (Krischer et al., 2015). Figure 3 shows the P and S velocities obtained with this method as a function of pressure. The main source of uncertainty in the velocities of P and S waves comes from the length measurements described above.

Figure 2 Example of a cross-correlation using our ObsPY-based Python script. (a) Recorded V_p data at 7 GPa and 874 K. The black and red echoes are from both ends of the sample. (b) Superposition of the red and black echoes resulting from the cross-correlation with the indication of the travel time, Δt. (c) Result of the cross-correlation showing the unique maximum that precisely determines the round-trip travel time of the P wave across the specimen.

3 Results

In this section we will discuss both X-ray and ultrasonic interferometry results. It is well accepted that the latter approach is more reliable, in particular because it does not need an external pressure scale such as NaCl EoS (e.g., Kung et al., 2002). Nevertheless, a comparison of the two datasets provides clues on the consistency of our results.

We removed data acquired during the cold-pressing step of our experiment in both X-rays (already presented in Béjina et al., 2019) and sound velocity sets because they present the highest differential stress levels, and these may affect sound velocity measurements.

The microstructure of our sample was also described in our previous publication, and we will only emphasize that our fayalite remained stable throughout the experiment. XRD measurements and postmortem SEM-EBSD observations showed no signs of transformation in our specimen.

3.1 Stress, pressure and unit-cell volume determination

The 10-detector array and other high-pressure equipment we used at the NSLS beamline X17B2 have been described in detail by Weidner et al. (2010). The stress estimation was based on the measurement of lattice strains, a method detailed by Singh et al. (1998). We therefore refer readers to these previous publications (and references therein) for more information on evaluating differential stress in high-pressure experiments. For each fitted hkl peak the measured d spacings, d_m, were determined from the energy peak position on all available detectors (nine detectors subend 180^∘ and 22.5^∘ apart and one at 270^∘). The values d_P(hkl), representing d at hydrostatic pressure and used to estimate all unit-cell volumes, were calculated from the following relation (Singh et al., 1998, their Eqs. 3 and 4):

$$\begin{array}{}\text{(1)}& d_{\mathrm{m}}\left(hkl\right)=d_P\left(hkl\right)\left[\mathrm{1}+\left(\mathrm{1}-\mathrm{3}{\cos }^{\mathrm{2}}\mathit{\psi }\right)Q\left(hkl\right)\right],\end{array}$$

where d_m(hkl) is the measured d spacing of a plane (hkl), ψ is the detector azimuth, and Q(hkl) is a term proportional to differential stress (t):

$$\begin{array}{}\text{(2)}& Q\left(hkl\right)={\displaystyle \frac{t}{\mathrm{3}}}\left(\mathit{\alpha }{\left[\mathrm{2}{G}_{\mathrm{R}}^X\left(hkl\right)\right]}^{-\mathrm{1}}+(\mathrm{1}-\mathit{\alpha }){\left[\mathrm{2}{G}_{\mathrm{V}}\right]}^{-\mathrm{1}}\right).\end{array}$$

In the above equation, $G_{\mathrm{R}}^X$ is the aggregate shear modulus calculated under the Reuss condition (iso-stress), G_V is the shear modulus under the Voigt condition (iso-strain), and their respective expression is to be found in Singh et al. (1998) and Hearmon (1956). The weight factor α=1 (or =0) corresponds to the pure Reuss (or pure Voigt) bound. Our estimation of t was made assuming pure Reuss condition according to Chen et al. (2006).

Examples of plots, $d_{\mathrm{m}}=f(\mathrm{1}-\mathrm{3}{\cos }^{\mathrm{2}}\mathit{\psi })$, are presented in Fig. S3. The value of d_P(hkl) corresponds to d when 3cos ²ψ=1. To estimate the differential stress both in NaCl and Fa we preferred an individual-peak fitting procedure instead of refining the unit-cell parameters by fitting the full measured diffraction patterns with modeled spectra. NaCl peaks were fitted using the software plot85 (Weidner et al., 2010), while fityk (Wojdyr, 2010) was used for Fa peaks. All unit-cell volumes were calculated using both the UnitCell software (Holland and Redfern, 1997) and our own set of Python scripts that make use of the LMFIT minimization package (Newville et al., 2014), also used for the stress calculations. Both procedures gave identical volumes within uncertainties, and UnitCell diagnostics helped filter poorly fitting peaks. In the end, unit-cell volumes of NaCl at high pressure and high temperature were calculated from the d spacings of peaks (200), (220), (222), (400), (420) and (422), and for Fa we kept peaks (061), (112), (130), (131), (132), (133), (134), (152), (211), (222), (241) and (422).

Pressure was first estimated from the NaCl volumes using a third-order Birch–Murnaghan (BM) EoS (Birch, 1947) with the isothermal bulk modulus, K_T, $\mathrm{d}{K}_T/\mathrm{d}P$ and $\mathrm{d}{K}_T/\mathrm{d}T$, from Brown (1999), as recommended by Strässle et al. (2014), and thermal expansion from Spetzler et al. (1972). This pressure was then corrected for the differential stresses, t, measured in both NaCl and Fa (Fig. 4) following the next relation (Béjina et al., 2019):

$$\begin{array}{}\text{(3)}& P_{\mathrm{Fa}}=P_{\mathrm{NaCl}}-{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}(t_{\mathrm{Fa}}-t_{\mathrm{NaCl}}),\end{array}$$

where P_Fa is the pressure felt by the Fa sample, P_NaCl is the pressure calculated from NaCl volumes, and t_Fa and t_NaCl are stresses measured in the Fa sample and in NaCl respectively (detailed data are listed in Table S1). Stress was calculated from d spacings measured at different azimuths of the detector array using Eq. (1) and following the method of Singh et al. (1998). As shown in Fig. 4, differential stress in Fa is mostly compressive at the highest pressures but mostly extensive below about 4 GPa. In comparison NaCl presents less variation in t, and consequently P_Fa is slightly higher than P_NaCl below 4 GPa and lower above (Fig. S4). Table 1 lists unit-cell parameters calculated from the above d_P values as a function of pressure and temperature.

Figure 3 Sound wave velocities as a function of pressure, P_Fa, for V_p (a) and V_s (b).

Figure 4 (a) Differential stress calculated in both NaCl and fayalite using the method of Singh et al. (1998) (compressive stresses are marked as positive). (b) Pressure, P_Fa (left axis and grey squares), determined from the NaCl unit-cell volumes and corrected for the differential stress. Right axis and red dots show temperature. Data from the cold-pressing step were not included in the determination of the elastic parameters, as discussed in the text.

Table 1 Unit-cell parameters of Fa and velocities V_p and V_s measured as a function of P_Fa (P_NaCl corrected for differential stress) and temperature. First column is the Fa XRD pattern numbers. Uncertainties are 1σ.

3.2 Uncertainties

Uncertainties (1σ) in individual peak positions, hence on d_hkl spacing values (Δd_hkl), were directly obtained from the Levenberg–Marquardt fit of the diffractograms in fityk. Because there is no straightforward way to include error bars on data in LMFIT (one can use weighting factors, but there are different ways to estimate them), we used the d−Δd and d+Δd of each hkl peak to calculate the minimum and maximum values of cell parameters and volumes. The uncertainty in V (ΔV) for each (P,T) condition was defined as the maximum value between V−V_min and V_max−V, with V being the nominal volume calculated from d.

Uncertainties in P_NaCl, t_Fa and t_NaCl (Table S1) were obtained from those on unit-cell volumes and d spacings by error propagation. Uncertainties in P_Fa were then determined according to Eq. (3). These were finally included in the weighted fitting routine of EosFit7 (Angel et al., 2014) in order to obtain the fayalite EoS parameters.

This approach leads to larger errors in V and subsequently in all parameters derived from either NaCl or fayalite volumes compared to a purely statistical approach (for example a weighted fit) but could very well be more representative of the real ΔV. Indeed, in many high-pressure studies, underestimated errors are often attributed to pressure and temperature mainly because pressure and thermal gradients are poorly known, the position of the thermocouple with respect to the sample could change with pressure and/or temperature, uncertainties in the EoS used to estimate pressure are not taken into account, etc.

3.3 Isothermal equation of state of fayalite

Unit-cell volumes (Table 1) as a function of stress-corrected pressure and temperature were fitted simultaneously in EosFit7. We used both the second- and third-order Birch–Murnaghan (BM) EoS with a linear expression of $\mathrm{d}{K}_T/\mathrm{d}T$. Because high-temperature data measured at high pressure rarely have the precision of those at 1 atm, the modified Holland–Powell expression was chosen to calculate the thermal expansion (Angel et al., 2014). An acceptable value of a fitted $K_T^{\prime }$, 5.5 for K_T=125.2 GPa, could be obtained after removing the room-temperature data, but uncertainties in all EoS parameters are large. Therefore, for the third-order BM EoS we fixed $K_T^{\prime }$ at 5.3 (Nestola et al., 2011) to accommodate the small pressure range. All results are listed in Table 2, and one can see that K_T values are compatible with our previous value, $K_T=\mathrm{131.9}\pm \mathrm{1.5}$ GPa (Béjina et al., 2019), within the error bars. This K_T was obtained using all room-temperature data, including the pressurization cycle, and therefore with a much higher stress (Fig. 4). The compatibility between room-temperature and high-temperature data, with different stress levels in both, is a hint that our pressure correction leads to reliable results.

Table 2 Birch–Murnaghan EoS parameters determined from our XRD data at all temperatures using EoSFit7 with the modified Holland–Powell equation for thermal expansion (Angel et al., 2014).

^a implied value; ^b fixed value.

As mentioned above high-pressure experiments rarely have precise temperature measurements. Consequently the thermal parameters are not fully reliable and present large uncertainties. Our expression of the thermal expansion tends to overestimate volumes at high temperature compared to Kroll et al. (2012), a point that we discuss in Sect. 4.1. To circumvent this difficulty and test the validity of our values for K_T, we fitted data at each temperature independently using a third-order BM EoS with $K_T^{\prime }=\mathrm{5.3}$ (Fig. 5). Obtained values of V₀ and K_T extrapolated to room pressure are listed in Table 3 as a function of temperature. From these values, an overall $K_{T(\mathrm{0},\mathrm{298})}=\mathrm{129}\pm \mathrm{5}$ GPa was calculated from a linear fit, K_T=f(T), and $V_{(\mathrm{0},\mathrm{298})}=\mathrm{307.9}\pm \mathrm{0.2}$ ų was determined from the modified Holland–Powell equation (Angel et al., 2014). Despite the dispersion in the data we found a K_T(0,298) value compatible with the EoSFit7 results within the uncertainty, therefore validating the simultaneous fit to the whole dataset.

Figure 5 Cell volume of our Fa sample as a function of pressure, P_Fa, at all temperatures. P_Fa is corrected from differential stress as described by Eq. (3). Each fit is an isothermal third-order Birch–Murnaghan equation of state (see results in Table 3).

Table 3 Third-order Birch–Murnaghan EoS parameters obtained from a weighted least-square fit of XRD data at each temperature; $\mathrm{d}{K}_T/\mathrm{d}P$ was fixed at 5.3 (Nestola et al., 2011).

3.4 K_s, $\mathrm{d}{K}_s/\mathrm{d}P$, G and $\mathrm{d}G/\mathrm{d}T$ for fayalite

Following the procedure described, for example, in Kung et al. (2002) and using finite strain equations from Davies and Dziewoński (1975), we calculated K_S0, G₀ and their pressure derivatives at each temperature (Table 4). This method only uses the density of the sample and sound wave velocities and does not need an external pressure gauge such as NaCl (see also Li and Zhang, 2005). The formulation of Kung et al. (2002), thereafter named K02, is a truncated version of Davies and Dziewoński (1975), DD75, this latter giving access to second derivatives of K_S and G with respect to pressure. Our data do not have the resolution necessary to obtain reliable values for these, but they nevertheless allow for the calculation of pressure. A comparison of the pressure calculated by the finite strain equations with P_Fa (Eq. 3) shows a remarkable agreement (Fig. 6). One can also see in Fig. 6 that uncorrected pressure is more discrepant. It is therefore not surprising that fitting $G=\mathit{\rho }{V}_{\mathrm{s}}^{\mathrm{2}}$ and $K_S=\mathit{\rho }{V}_{\mathrm{p}}^{\mathrm{2}}-\mathrm{4}G/\mathrm{3}$ to the data at each temperature using the NaCl pressure scale gives results very close to those from the finite strain method (Table 4).

Figure 6 Pressure calculated using the finite strain equations (P_fs) as a function of P_Fa obtained from unit-cell volumes of NaCl (using EoS by Brown, 1999) and corrected for differential stress in both Fa and NaCl. The line shows the 1:1 ratio (not a fit). Grey crosses show P_NaCl not corrected for stress.

Table 4 Adiabatic bulk modulus, shear modulus and their pressure derivatives from fitting finite strain equations at each temperature using the full form of Davies and Dziewoński (1975), DD75, or their truncated version from Kung et al. (2002), K02, or the stress-corrected NaCl pressure scale, P_Fa.

Table 5 lists the values of all elastic parameters at 298 K calculated after a linear fit of data from Table 4 as a function of temperature. While we list the results from both DD75 and K02 formulations, we only retain the values from the K02 equations because of the large uncertainties around the DD75 results. One can nevertheless note that the two formulations give similar results and that K02 values for K_S, $\mathrm{d}{K}_S/\mathrm{d}P$, G and $\mathrm{d}G/\mathrm{d}P$ are nearly identical to those obtained when fitting the data using P_Fa. On Fig. 7 we show the fit of K02 data (Table 4) as a function of T. Points at the highest temperatures may seem to diverge from the trends, but removing them does not significantly change the results of the weighted least-square fit. The value of K_S we calculated at 298 K is very close to the one found for K_T. The relation $K_S/K_T=\mathrm{1}+\mathit{\gamma }\mathit{\alpha }T$, with γ the Grüneisen parameter, α the thermal expansion and T the temperature, implies a K_T about 1 % smaller than K_S, a difference too small to resolve here because of our uncertainties.

Figure 7 (a, b) Adiabatic bulk modulus and its pressure derivative as a function of temperature. (c, d) Same for shear modulus. Fits are performed with the formulation by Kung et al. (2002). Solid blue lines are the results of a weighted least-square fit to the data, and dashed lines show the 95 % confidence interval. Results (labeled as K02) are listed in Table 5.

Table 5 Results from a linear least-square fit of data (Table 4) with T using the formulations of Davies and Dziewoński (1975), DD75, and Kung et al. (2002), K02, and compared with fitting the data using P_Fa. Uncertainties are 1σ calculated from the covariance of the minimization. K and G are expressed in gigapascals, T in kelvin, and P_Fa in gigapascals.

4 Discussion

4.1 The hard case of fayalite

4.1.1 Oxygen fugacity

Physical properties of Fe-rich minerals are always a challenge to study, in particular at high pressures. For example, variations in the ${\mathrm{Fe}}^{\mathrm{2}+}/{\mathrm{Fe}}^{\mathrm{3}+}$ ratio may cause volume changes that in turn may perturb measurements. A non-exhaustive search of experimentally measured unit-cell volumes (V₀) of Fa at ambient conditions shows values between 306.9 ų (natural sample, Richard and Richet, 1990) and 310.1 ų (synthetic crystal, Graham et al., 1988), a range more than twice as large as that of forsterite, the pure Mg endmember of olivines. In a 1 atm furnace, performing annealing runs at constant oxygen fugacity (fO₂), and therefore at a constant ${\mathrm{Fe}}^{\mathrm{2}+}/{\mathrm{Fe}}^{\mathrm{3}+}$ ratio, is feasible with proper gas mixtures. At high pressures, solid buffers are used to maintain the fO₂ along a chemical equilibrium, often metal/metal oxide, and consequently the fO₂ changes with temperature and pressure. It is therefore possible that the range of V₀ mentioned above reflects the formation conditions and/or the thermal history of the different Fa specimens.

Our volume V₀ fitted from high-temperature data down to room temperature (Table 2) is very close to that of the starting material (Béjina et al., 2019), but, nevertheless, our high-temperature volumes are larger than what is expected from thermal expansion expressions (e.g., Fei, 1995; Kroll et al., 2012). There are very few studies on point defects in minerals in relation to their unit-cell volumes and elastic properties. To our knowledge, among minerals present in planetary mantles, only defects in wüstite, Fe_xO, have been studied in detail (McCammon and Liu, 1984; Zhang, 2000; Zhang and Zhao, 2005). These show that the cell parameter of Fe_xO increases as the mineral tends toward perfect stoichiometry, i.e., as Fe²⁺ stoichiometry tends to 1. However, a comparison between wüstite and fayalite should be made with caution because the point-defect populations are very different in the two minerals. FeO is known for its non-stoichiometry allowing for the accommodation of high levels of metal vacancies (V_Fe) which form clusters of Fe³⁺ surrounded by V_Fe (e.g., Minervini and Grimes, 1999), clusters that are not present in fayalite (e.g., Nakamura and Schmalzried, 1983). On the other hand, in a non-stoichiometric Fe₂SiO₄ wadsleyite, Hazen et al. (2000) proposed that the substitution of 2Fe³⁺ for ${\mathrm{Fe}}^{\mathrm{2}+}+\mathrm{Si}$ increases its unit-cell volume. A parallel could be made with the defect ${\mathrm{Fe}}_{\mathrm{Si}}^{\mathrm{•}}$ (Fe³⁺ on a Si⁴⁺ site) proposed in fayalite (Nakamura and Schmalzried, 1983; Hillert et al., 1996) that probably increases the size of the tetrahedron. Fe³⁺-related defects could very well lead to larger unit cells in Fe-bearing silicates. Like for most Fe-rich minerals, it is difficult to prepare perfectly stoichiometric fayalite. In our high-pressure experiment, our sample was placed in a BN sleeve, a very reducing material, but also in contact with Al₂O₃ and NaCl, making it difficult to estimate the oxygen fugacity in the cell assembly and therefore the nature of the point-defect population. One must also note that our maximum temperature is, so far, the highest at which Fa unit-cell volumes have been measured at high pressure.

4.1.2 Shear modulus

Our value of the shear modulus at room conditions, $G_{(\mathrm{0},\mathrm{298})}=\mathrm{53.3}\pm \mathrm{0.4}$ GPa, is close to previous data (Fig. 8). This is the first estimate of G and G^′ obtained from ultrasonic interferometry and on an aggregate since the study by Chung (1971), all other studies having used Fa single crystals (Graham et al., 1988; Speziale et al., 2004). Except for the G^′ data point of Chung (1971) that is very low, all G^′ are within the range [1.2,1.8]. It is interesting to note that G seems less sensitive to the sample's environment than K. For example, it is accepted that S waves are less sensitive to differential stress than P waves, especially if the wave propagation direction is aligned with the maximum principal stress. Keeping in mind that $G=\mathit{\rho }{V}_{\mathrm{s}}^{\mathrm{2}}$, it is not surprising that G is less affected by differential stress.

Figure 8 Trend of two elastic parameters for fayalite with their respective pressure derivatives. (a) Bulk modulus versus $\mathrm{d}K/\mathrm{d}P$. The blue color is for K_S values, grey for K_T and red for this study. K_S and K_T values are plotted as published, with no conversion from one to the other. The full red square is our data as reported in Table 4 (P_Fa), the open square being K_S fitted to our data with fixed $K_S^{\prime }=\mathrm{5.3}$. Arrows at the end of error bars of the data of Yagi et al. (1975) symbolize that they extend outside the plot range. (b) Shear modulus versus $\mathrm{d}G/\mathrm{d}P$. On both plots, a given symbol refers to the same publication. Dotted lines represent values of K_S and G determined at room pressure; therefore, no pressure derivative was determined.

4.1.3 Bulk modulus

Comparing bulk moduli between studies is always difficult and more so than for G. For Fe–Mg olivines, G decreases linearly with X_Fe by about 30 GPa over the entire solid solution, whereas K varies by only a few gigapascals with no clear trend with X_Fe. As discussed, for example, in Angel (2000), the fitted values of K and K^′ are often strongly correlated, and therefore, fixing either one affects the estimation of the other. This can be seen in Fig. 8 representing experimental values of K of fayalite plotted against $\mathrm{d}K/\mathrm{d}P$. The figure shows a clear trend with a negative slope, meaning that as K^′ increases, K decreases. Many high-pressure experiments, in particular in multi-anvil, only explore a limited pressure range. For such experiments, analyzing data often requires one to fix either K or K^′ and calculate the other from a fit to the data. In our case, for example, Table 2 shows that for $K_T^{\prime }$ fixed at 4 (second-order BM EoS) and 5.3 (Nestola et al., 2011), K_T varies by almost 3 GPa.

Another source of difficulty when comparing EoS parameters comes from the choice of a pressure scale. Historically the NaCl Decker scale (Decker, 1971) has been the most used one, but Strässle et al. (2014) proposed that the scale of Brown (1999) is more accurate. When measuring sound wave velocities, the finite strain approach (Davies and Dziewoński, 1975) allows for the determination of the elastic parameters without prior knowledge of the pressure, which can be calculated afterwards (see their Eq. 1c). These different methods can lead to small changes in pressure estimates with a strong impact on the determination of an EoS, in particular for $\mathrm{d}K/\mathrm{d}P$.

Among other factors affecting the determination of K are, as mentioned above, the sample's stoichiometry and therefore the fO₂ conditions of the experiment, as well as the (often ignored) level of non-hydrostaticity of most high-pressure experiments.

Our result for K_T measured at room temperature, $K_{T(\mathrm{0},\mathrm{302})}=\mathrm{130}\pm \mathrm{1}$ GPa for $K_T^{\prime }=\mathrm{5.3}$ (Table 3), is within the uncertainties of our previous estimate (Béjina et al., 2019). Our choice of $K_T^{\prime }$ comes from the high accuracy of the data of Nestola et al. (2011) measured for olivines with X_Fe between 8 % and 38 %, close to estimates, $K_S^{\prime }=\mathrm{4.9}$ (Speziale et al., 2004) and 5.2 (Graham et al., 1988), for fayalite. When recalculating K_T using EoSFit7 and high-temperature data as well, we obtain $K_{T(\mathrm{0},\mathrm{298})}=\mathrm{129}\pm \mathrm{5}$ GPa for the same $K_T^{\prime }=\mathrm{5.3}$ (Table 2). From the relation between K_S and K_T given in Sect. 3.4, one would expect K_S≃131 GPa with $K_S^{\prime }=\mathrm{5.3}$ at 300 K. The values we deduced from our dataset are $K_S=\mathrm{127.2}\pm \mathrm{0.3}$ GPa and $K_S^{\prime }=\mathrm{6.5}\pm \mathrm{0.1}$ (Table 5). Fixing $K_S^{\prime }$ at 5.3 would increase K_S to about 130 GPa, a value very consistent with our estimate of K_T (Table 2). These values follow the trend $K_S=f(\mathrm{d}{K}_S/\mathrm{d}P)$ from previous studies (Fig. 8) and therefore are in general agreement with these. From this figure, K_S and $\mathrm{d}{K}_S/\mathrm{d}P$ of fayalite seem to fall in a limited range, [125,135] GPa and [4,6.5] respectively.

Figure 9 P wave (top) and S wave (bottom) velocities in fayalite calculated with burnman (Cottaar et al., 2014) using the formalism of EoS of Stixrude and Lithgow-Bertelloni (2011) (SLB2011). Full lines are Mars profiles along the hot areotherm A and dashed lines along the cold areotherm C of Bertka and Fei (1997). Dotted lines are along a geotherm (Priestley and McKenzie, 2013) included in burnman. Pressure and temperature profiles are shown in the inset. Red lines are velocities calculated using the elastic parameters from this work and grey curves those for parameters from Stixrude and Lithgow-Bertelloni (2011). S wave velocities in ahrensite are also shown to visualize the velocity jump when transiting to the spinel structure.

It is in general difficult to pinpoint a main reason for low values of K_S even when correlated with high $K_S^{\prime }$. Among a variety of reasons, we possibly have the wrong temperature reading of the thermocouple, the effect of oxygen fugacity and hence point defects on the elasticity of the olivine structure, the effect of differential stress (Yagi et al., 1975), compositional inhomogeneities or microstructure irregularities, etc. In our case, one could argue that the transition of our Fa specimen to ahrensite (Ahr) may have impacted our estimates of the elastic parameters. However we emphasize that no Ahr was observed by SEM, no distinct Ahr peaks can be seen on the XRD patterns, and sound velocities do not increase when approaching the phase boundary as seen by Liu et al. (2010).

In their study of nanocrystalline forsterite (11 nm grain size prior to the experiment), Couvy et al. (2010) proposed that a small grain size (i.e., many grain boundaries and triple junctions) leads to a lower K_T. This effect has been observed in different materials and for grain sizes usually below 100 nm (e.g., see Marquardt et al., 2011, for MgO and references therein for other materials). To our knowledge no systematic studies have been performed, in particular on olivines, to quantify this effect of grain size on the bulk modulus. The grain size of our Fa sample, prepared from nano-size oxide powders, was about 1 µm prior to the experiment and around 3 µm after. Even though much larger than 100 nm, we cannot rule out that the grain size of our sample could be small enough to start affecting its elastic properties. This assumption is difficult to verify because postmortem microstructures are rarely described in most studies.

4.2 Seismic velocities in Mars' upper mantle

In order to visualize the effect of fayalite elastic parameters on P and S wave velocities, we calculated V_p and V_s using the Python package burnman (Cottaar et al., 2014) with the formalism of Stixrude and Lithgow-Bertelloni (2011) (SLB2011) and data for fayalite (K_T0=135 GPa, $K^{\prime }=\mathrm{4.22}$, G=50.9 GPa, $G^{\prime }=\mathrm{1.46}$) and then with our determination of these parameters (Table 5). We kept the temperature-dependent parameters of SLB2011 because of the large uncertainties in ours. Figure 9 represents the P and S velocities in Fa as a function of pressure along two areotherms for Mars (profiles A and C of Bertka and Fei, 1997) and a geotherm from the CAM2016 model (Priestley and McKenzie, 2013). We chose this geotherm in order to have a very different temperature profile compared to the two areotherms. This figure shows that the differences in P and S velocities with the change in parameters are very small. For the geotherm the effect of the hot temperatures prevails over that of pressure, and we observe a decrease in the velocities up to about 5 GPa.

Figure 10 Velocity profiles in Mars' upper mantle calculated with burnman using the formulation of Stixrude and Lithgow-Bertelloni (2011). The black dots represent the velocities for fayalite using our K_S, $\mathrm{d}{K}_S/\mathrm{d}P$, G and $\mathrm{d}G/\mathrm{d}P$. The red, green and blue domains are the ranges of velocities calculated when varying thermoelastic parameters, as described in the text.

To further test the influence of elastic parameters of Fa on its P and S velocities, we performed the same calculations along the areotherm A from Bertka and Fei (1997), allowing all parameters from the SLB2011 EoS to vary in the intervals: $K_{T\mathrm{0}}=[\mathrm{125},\mathrm{135}]$ GPa, $G_{\mathrm{0}}=[\mathrm{50},\mathrm{53}]$ GPa, $K^{\prime }=[\mathrm{4},\mathrm{6}]$, $G^{\prime }=[\mathrm{1},\mathrm{1.3}]$, the Debye temperature $\mathrm{\Theta }=[\mathrm{550},\mathrm{650}]$ K, the Grüneisen parameter $\mathit{\gamma }=[\mathrm{1},\mathrm{1.3}]$ and exponent $q=[\mathrm{3},\mathrm{4}]$. Figure 10 shows V_p, V_s and $V_{\mathrm{p}}/V_{\mathrm{s}}$ along this areotherm for Fa and two olivines: one with 25 mol % Fe and one with 40 mol % Fe. The EoS for these olivines are calculated by mixing Fe and Mg olivine endmembers according to their molar fraction (Cottaar et al., 2014). In our calculation, only Fa parameters vary, those of forsterite (Fo) being considered as perfectly determined. Consequently mixing Fo and Fa parameters to calculate velocities results in smaller domains of possible velocities compared to pure Fa, but the difference in velocities between these two olivines is rather small. This could be a challenge when fitting data from the SEIS instrument.

Figure 11 Optimum olivine contents determined by adjusting V_p of a simple upper-mantle model to PREM using burnman for three K_T values for fayalite. Large variations in K_T produce small changes in olivine fractions. See text for calculation details.

Another way to look at this is to test the effect of the bulk modulus on the optimized olivine content resulting from a fit of a simple model of Earth's upper mantle to the Preliminary Reference Earth Model (PREM) (Davies and Dziewoński, 1975) using burnman. We only test here the sensitivity with K because its experimental determination presents a much larger dispersion than G. We built a very simple model of the Earth's upper mantle with three mineral phases: olivine, orthopyroxene (Opx) and garnet (Gt), all three being simple Fe–Mg solid solutions. The molar percentage ratio of Opx to Gt was fixed at 2:1, and only the olivine fraction was optimized. These calculations were done for three values of the olivine K_T0, 120, 130 and 140 GPa, which correspond approximately to the range of experimental values determined for this mineral. The properties of this modeled rock were calculated as described by Cottaar et al. (2014) along a geotherm (Anderson, 1982, provided by burnman). Figure 11 shows that the three optimum olivine contents fall in a tight range despite the large variation in K_T0, between 61 mol % and 66 mol %. Therefore, considering the uncertainties in the effect of Fe on the bulk modulus of olivine, it is difficult to say at this point if the SEIS data will have the necessary resolution to determine the Fe concentration of the minerals of Mars' upper mantle with an uncertainty less than 15 %.

5 Conclusions

The bulk and shear moduli of fayalite presented here are the first to be determined by coupled synchrotron X-ray radiation and ultrasonic interferometry. A fit using the third-order Birch–Murnaghan EoS with $K_T^{\prime }=\mathrm{5.3}$ (fixed) to the data gives $K_T=\mathrm{128.6}\pm \mathrm{5.1}$ GPa at room conditions, whereas $K_S=\mathrm{127.2}\pm \mathrm{0.3}$ GPa and $K_S^{\prime }=\mathrm{6.5}$ were obtained from sound wave velocities. Fixing $K_S^{\prime }$ to 5.3 increases K_S to 130 GPa, therefore showing that our K_T and K_S are compatible within uncertainties. These values of K and K^′ follow the linear trend $K=f\left(K^{\prime }\right)$ defined by previously published data. The shear modulus values of fayalite show less dispersion than the bulk modulus; therefore, the latter will be critical when fitting a mineralogical model to Mars InSight SEIS data. The range of K values for the olivine Fe–Mg solid solution is defined between about 125 and 135 GPa. With such variation, a precise compositional and mineralogical model of the Mars mantle will be challenging to obtain. Because experiments on Fe-rich minerals are always difficult, more experiments are needed to better define the trend of K with X_Fe so that more robust values of K can be established for minerals with an Fe content consistent with Mars' mantle composition.