New temperature and oxygen fugacity data of Martian nakhlite from Northwest Africa (NWA) 5790 and implications for shallow sulphur degassing

Newly analysed titanomagnetite–ilmenite (Tim–Ilm) intergrowths from Martian nakhlite meteorite Northwest Africa (NWA) 5790 yielded crystallisation temperature up to 1032 °C and oxygen fugacity (fO2) up to ΔQFM + 1.6, notably higher than previous estimates for nakhlite magmas (temperature < 950 °C, fO2 = ΔQFM − 0.5 to ΔQFM + 1). To interpret how the magma was reduced from ΔQFM − 0.5 to ΔQFM + 1.6, we used D-Compress to model the sulphur degassing process within a single thick lava pile. For fO2 to significantly decrease in this extended range, a sulphur-rich (S content 4000–7000 ppm) Martian lava flow had to degas all the sulphur species at a certain final degassing pressure, which was 2–4 bar for NWA 988 and Lafayette and < 0.7 bar for Y-000593 and Nakhla. These final degassing pressure data are in good agreement with the Martian nakhlite burial depth estimated by other petrological and geochemical methods. These estimates are also comparable with the excavation depth of ~ 40 m based on the small (6.5 km in diameter) impact crater over the Elysium lava plain. The fO2-controlled sulphur degassing pressure may constitute a method for estimating the burial depth of sulphur-rich lava flows on Mars.


Introduction
Determining the oxygen fugacity (fO 2 ) of Martian shergottite-nakhlite-chassignite meteorites is important to understand the volcanic eruption process in ancient Mars. Both shergottite and chassignite have relatively reduced fO 2 , ranging from ΔQFM − 4 (4 logfO 2 units below quartz-fayalite-magnetite buffer) to ΔQFM − 0.5 (Hewins et al. 2020;Udry et al. 2020). However, fO 2 in this range does not appear to represent magma erupted to the surface of Mars due to the high degree of crystallisation of shergottite and chassignite. Nakhlite has a petrographic texture similar to terrestrial basalts, believed to have been excavated from lava piles buried at the shallow surface of Mars (Mikouchi et al. 2006(Mikouchi et al. , 2012Shuster and Weiss 2005;Righter et al. 2008;Cohen et al. 2017). Therefore, the fO 2 of nakhlite samples can directly provide constraints on lava flows erupted from ancient Martian volcanoes.
There is a considerable variation of fO 2 measured from Fe-Ti oxide phenocrysts in nakhlite samples, which ranges from ΔQFM − 0.5 to ΔQFM + 1 (Szymanski et al. 2010) for those ilmenite-titanomagnetite pairs that passing the Mg/Mn equilibrium test of Bacon and Hirschmann (1988) (for details see Additional file 1: Figure S1 and Table S1). However, the full range of fO 2 in nakhlites has not yet been fully resolved, as Fe-Ti oxides studied likely represent only the relatively low-temperature history of these rocks. The closure temperature (Tc) corresponding to the highest observed fO 2 is < 950 °C, much lower than the crystallisation temperature of clinopyroxene phenocrysts in these melts (1100 °C, Righter et al. 2008). Most studies applied geothermometers and oxybarometers to ilmenite-titanomagnetite pairs that occurred as veinlets and marginal growth along the cracks, which are considered to be late features (e.g., Szymanski et al. 2010;Treiman and Irving 2008), consistent with the low temperatures calculated. Ilmenite exsolution lamellae intergrowths have been reported in Fe-Ti oxide phenocrysts, which may record higher temperature growth (e.g., Balta et al. 2017;Righter et al. 2014;Imae et al. 2005), but none of these intergrowths has been successfully analysed by electron probe microanalysis (EPMA) due to the inadequate thickness of the ilmenite lamellae.
Containing the highest mesostasis (glass phase) proportion, NWA 5790 is the nakhlite sample with the highest cooling rate (Corrigan and Velbel 2015;Mikouchi et al. 2012). Therefore, ilmenite-titanomagnetite pairs in NWA 5790 should record the highest closure temperature. We identified titanomagnetite-ilmenite intergrowths in Fe-Ti oxide phenocrysts in nakhlite NWA 5790, which were wide enough (> 2 μm) to locate a single EPMA spot. Then, we determined the equilibrium temperature and fO 2 of these intergrowths. We conducted sulphur degassing modelling based on the newly obtained fO 2 data and the extended variation range of fO 2 . The modelling results helped us constrain the degassing pressure during the crystallisation of Fe-Ti oxide phenocrysts and the emplacement depth (i.e. burial depth) of different nakhlites.

EPMA analysis
Quantitative chemical analyses of minerals in the polished thin section NWA 5790 were conducted using a JEOL JXA-8230 electron microprobe analyzers (EMPA) equipped with four wavelength dispersive spectrometers at Key Laboratory of Orogenic Belts and Crustal Evolution of School of Earth and Space Sciences, Peking University. No field emission gun was used. The acceleration voltage and the beam current were 15 kV and 10 nA, respectively. The beam diameter was set to "spot" mode (minimum size of the instrument). Counting time of 20 s was used for both peck and background. The SPI 53 minerals standard (U.S.) was utilised for the quantitative analysis: sanidine was employed for K; diopside for Ca and Mg; rutile for Ti; jadeite for Na, Al and Si; chromium oxide for Cr; rhodonite for Mn; hematite for Fe; and nickel silicide for Ni. At the final calibration stage, the PRZ correction was performed. Detection limits of the oxides were 0.01 wt.% for K 2 O, 0.02 wt.% for Al 2 O 3 and Na 2 O, 0.03 wt.% for MgO and CaO, 0.04 wt.% for TiO 2 , Cr 2 O 3 , MnO, FeO, and NiO, and 0.06 wt.% for SiO 2 .

Ilmenite-titanomagnetite geothermometer
The exchange equilibrium of Fe 2+ + Ti 4+ for 2 Fe 3+ between ilmenite and titanomagnetite: has been calibrated in several studies for use as a geothermometer (Andersen and Lindsley 1988;Ghiorso and Sack 1991;Powell and Powell 1977;Sauerzapf et al. 2008). Of these proposed geothermometer models, only the thermometer model of Sauerzapf et al. (2008) was experimentally calibrated at a temperature range (800-1300 °C) that includes the temperature range (900-1100 °C) recorded in Fe-Ti oxides in NWA 5790. Thus, the geothermometer of Sauerzapf et al. (2008) was used in the temperature estimations in this study.

Petrography and P-T estimations
The NWA 5790 sample in this study contains clinopyroxene (55.7%), olivine (6.1%), Fe-Ti oxides (0.4%) and mesostasis (37.8%) ( Fig. 1a and Additional file 1: Figure S2). As shown in Fig. 1b, the Fe-Ti oxide grains are subhedral to anhedral, with a major axis of about 1 mm and a minor axis of about 0.5-1 mm. There are two occurrence domains in the titanomagnetite-ilmenite intergrowths: (1) pervasively distributed exsolution textures characterised by thin ilmenite lamellae occurring in titanomagnetite; and (2) ilmenite veinlets and rim regions along the cracks of titanomagnetite host crystals. EPMA was used to analyse ilmenites from these two domains and their host titanomagnetites. Ilmenite veinlets and rim regions are typically different in composition from lamella, with high TiO 2 , MgO, and MnO contents and low Al 2 O 3 , Cr 2 O 3 , and FeO contents (Table 1).
Titanomagnetite-ilmenite geothermometer and oxybarometer (Sauerzapf et al. 2008) were used to calculate the equilibrium temperature and fO 2 of Tim-Ilm lamellae pairs, which are 1005-1032 °C and ΔQFM + 1.60 to ΔQFM + 1.44, respectively (Table 1). As shown in Fig. 2a, the equilibrium temperature and fO 2 obtained in this study are remarkably higher than that calculated for Tim-Ilm pairs from other nakhlites, including Nakhla, Lafayette, Y-000593, and NWA 998 (Boctor et al. 1976;Bunch and Reid 1975;Sautter et al. 2002;Szymanski et al. 2010;Treiman and Irving 2008). The temperature of 1032 °C narrowed the gap between the crystallisation temperature of Fe-Ti oxide phenocrysts and clinopyroxene phenocrysts, which is estimated to be 1100 °C (Righter et al. 2008). This suggests that the Fe-Ti oxide and clinopyroxene phenocrysts may share a similar crystallisation temperature (> 1032 °C) and fO 2 (~ ΔQFM + 1.60). As shown in Table 1, the veinlet and rim ilmenites in the same sample yield an equilibrium temperature of 848-901 °C and fO 2 of ΔQFM + 1.26 to ΔQFM + 1.49, which are consistent with previous studies (Fig. 2a). This comparison confirms that the high-Ti veinlets and rim regions grew at a later stage, in contrast to the ilmenite lamellae.

Sulphur degassing modelling
During the evolution of silicate magma, any of several processes may change the Fe 3+ /ΣFe ratios in the melt, reflecting the observed fluctuations in fO 2 of nakhlite magma (Kress and Carmichael 1991). Firstly, the crystallisation of pyroxene and olivine increases the Fe 3+ /ΣFe ratio in the residual melt (i.e. oxidises the melt; Brounce et al. 2017). Secondly, degassing of H 2 O and CO 2 can oxidise the melt in open systems (e.g., Mathez (1984); Burgisser and Scaillet (2007); Humphreys et al. (2015)) or has no effect on melt fO 2 in closed systems (e.g., Kelley and Cottrell (2012); Waters and Lange (2016)). These processes are unlikely in nakhlite magmas, which indicate reduction of the magma as crystallisation progresses. Sulphur dissolves in silicate melts and takes the forms of sulphates (S 6+ ) or sulphides (S 2− ), while sulphur is present as SO 2 (S 4+ ) or H 2 S (S 2− ) in volcanic gas (Métrich et al. 2009). Therefore, sulphur degassing can lead to either oxidation or reduction of the residual melt phase, depending on the following reactions (Métrich et al. 2009): and Martian basalts have a sulphur content in the range of 4000-7000 ppm, much more abundant than Earth (1)  (Gaillard and Scaillet, 2009). Thus, we modelled the effect of sulphur degassing on fO 2 changes in nakhlite magma ( Fig. 2b and c). The modelling was based on a progressive degassing model of the C−O−H−S system using the software package D-Compress , which executes the gas−melt equilibrium model of Iacono-Marziano et al. (2012). Suppose that the system has some constraints (i.e. initial H 2 O and CO 2 contents, magmatic fO 2 ). In that case, based on the experimental calibration of the solubility of H 2 , H 2 O, CO, CO 2 , SO 2 , H 2 S, and S 2 in the silicate melt and the calculation of homogeneous equilibria in the gas phase for the same species, D-Compress calculates the concentration and speciation of C, H, O, and S in coexisting gas and silicate melt as a function of pressure and temperature (Brounce et al. 2017).
In this study, the nakhlite parent magma composition "NPM05" recently proposed by Sautter et al. (2012) was used. The initial pCO 2 was 400 mbar (Franz et al. 2019), the H 2 O content was 1.44 wt.% (Weis et al. 2017), and the initial fO 2 was ΔFMQ + 1.6 at 1100 °C. As D-Compress calculates the initial sulphur content of the melt based on implemented solubility laws (Longpré et al. 2017), and sulphur degassing is most sensitive to pressure variation (Gaillard and Scaillet 2009), the initial sulphur content of the melt (4000-7000 ppm) was controlled by  Fig. 2 a Closure temperature and fO 2 (relative to QFM buffer) for NWA 5790 nakhlite titanomagnetite-ilmenite pairs in this study (red dots), compared to closure temperature and fO 2 from titanomagnetite-ilmenite pairs in nakhlites from the literature (Boctor et al. 1976;Bunch and Reid 1975;Sautter et al. 2002;Szymanski et al. 2010;Treiman and Irving 2008). All closure temperatures and fO 2 were recalculated using geothermometer and oxybarometer in Sauerzapf et al. (2008). Error bars for temperature and fO 2 are ± 70 °C and ± 0.5log fO 2 unit, respectively. b Total sulphur content vs. fO 2 (expressed as ΔQFM) degassing simulations. c Pressure vs. fO 2 (expressed as ΔQFM) degassing simulations. Degassing models calculated by D-Compress are in an open C-S-O-H system. For close-system degassing curves, see Additional file 1: Figure S3 adjusting the pressure. Moreover, calculating the "sulphur content at sulfide saturation" (SCSS, O'Neill (2020)) of "NPM05" yielded 4787 ppm, which is in the range of 4000-7000 ppm. Thus, the initial melt sulphur contents used in this study were 4000 ppm, 4787 ppm, 6000 ppm, and 7000 ppm, corresponding to the pressures of 235 bar, 245 bar, 280 bar, and 335 bar, respectively. The main results of degassing modelling are: (1) in all realistic open and closed sulphur degassing systems, a net decrease is observed in the melt fO 2 upon decompression ( Fig. 2c and Additional file 1: Figure S3c); (2) sulphur degassing is predicted to begin after sulphur reaches saturation, but more effective sulphur degassing is achieved only at low pressures (< 50 bar) with a significant drop in fO 2 (Fig. 2b); (3) with increasing initial melt sulphur content, the melt becomes more reduced for a given pressure, and open-system degassing tends to promote this process more efficiently (Fig. 2b and c); (4) the melt fO 2 can decrease to ΔQFM − 0.5 at < 1 bar (Fig. 3), and the fO 2 variation from ΔQFM + 1.6 to ΔQFM − 0.5 can cover all the fO 2 data recorded in equilibrium titanomagnetite−ilmenite intergrowths from nakhlites (Fig. 1). These results are consistent with sulfur degassing controlling the range of observed fO 2 values.

Link between temperature, oxygen fugacity and burial depth
As shown in Fig. 2a, ilmenite-titanomagnetite pairs in different nakhlites recorded variable closure temperatures when the Fe-Ti oxides cooling. When the cooling rate (dT/dt) is low, the closure temperature recorded by the ilmenite-titanomagnetite pairs is also low, due to longer subsolidus re-equilibration process (e.g., Morgado et al. 2019;Hou et al. 2020). As discussed in Hou et al. (2020), the cooling rate (dT/dt) is proportional to Tc 2 / exp(1/Tc), i.e.
When Tc is large, the denominator e (1/Tc) is approximately 1, so (dT/dt) is proportional to the square of Tc, i.e.
In the half-space one-dimensional cooling model, Richter et al. (2016) has shown that the time scale (t) of cooling in any part of a lava flow is proportional to the square of depth (z), where z is the distance to the heat transfer boundary (i.e. burial depth). Therefore, the cooling rate at any height z would be proportional to the reciprocal of z 2 , i.e.
To sum up, we have a relationship between Tc and z via the cooling rate (dT/dt), which is: This means the relationship between fO 2 and Tc in Fig. 2a is governed by the burial depth (z). This relationship was built through the cooling rate (dT/dt) of different nakhlites: the nakhlite from a deeper burial depth has a higher degassing pressure; it also has a lower cooling rate because of the thick overlying lava pile of nakhlites (e.g., Treiman 2005;Daly et al. 2019). Therefore, the minimum fO 2 recorded in the Fe-Ti oxides of nakhlite samples reflects its final burial depth.

Implications on burial depth of nakhlite magma
The nakhlites should have crystallised in an open system, which is supported by the parallel REE patterns of calculated melts (Udry and Day 2018) and the bimodal F/ Cl ratios of two generations of apatite (McCubbin et al. Tc ∝ 1/z.   2013; Birski et al. 2019). Therefore, we can constrain the burial depth of a nakhlite by using its observed fO 2 range and the open-system fO 2 -pressure relationships modelled by D-Compress. For a closed system, the results are shown in Additional file 1: Figure S3. The minimum fO 2 recorded in the Lafayette and NWA 998 is ΔQFM − 0.1 (Fig. 2a), which corresponds to the decompressing pressure range of 2-4 bar (deep blue region in Fig. 3). Likewise, the minimum fO 2 recorded in the Y-000593 and Nakhla is ΔQFM-0.5 (Fig. 2a), which corresponds to the decompressing pressure range of < 0.7 bar (shallow blue region in Fig. 3). The results suggest that nakhlite NWA 998 and Lafayette crystallised at an emplacement depth (i.e. burial depth) of about 20-40 m, while Y-000593 and Nakhla crystallised at less than 7 m (Fig. 3). These emplacement depth constraints are broadly consistent with burial depths reported in previous studies (Mikouchi et al. 2003(Mikouchi et al. , 2006(Mikouchi et al. , 2012. However, nakhlites are improbable to reside at a pressure of less than 0.04 bar (0.4 m), as 0.04 bar is the upper limit of the Martian atmospheric pressure (Ward 1974). Due to insufficient data, the burial depth of NWA 5790 cannot be constrained in this study. However, the pressure range of 235-335 bar obtained by D-Compress reflects the crystallisation depth of phenocrysts (2.4-3.4 km). Thus, a parental melt of nakhlite was formed and began to crystallise ilmenite lamellae in Fe-Ti oxide phenocrysts at a depth of > 2.4-3.4 km. The melt then ascended to a depth of < 40 m or erupted onto the surface, forming late-stage veinlets and rim ilmenites.
The consistency between the burial depth constrained by fO 2 and the modelled sulphur degassing pressure, is consistent with the scenario that the nakhlites originated from the shallow levels of a thick lava flow or a staging lava pile at different depths, which was expressed as the single volcanic event model (e.g., Mikouchi et al. 2003;Treiman 2005;Day et al. 2006;Imae and Ikeda 2007;Udry et al. 2012;McCubbin et al. 2013;Righter et al. 2014;Jambon et al. 2016;Richter et al. 2016;Dottin et al. 2018). This occurrence is also consistent with the excavation depth of ~ 40 m using the small impact crater (6.5 km in diameter) over the Elysium lava plain (Cohen et al. 2017), although their "discrete eruptive events model" is not consistent with the covariation between fO 2 and burial depth.

Error evaluation
At all pressures, the solubility coefficients of volatile can lead to significant relative errors in degassing model (Iacovino et al. 2021). Relative errors induced by extreme solubility coefficients for estimating the gas molar ratio would raise from ~ 30 to 40-50% for basalt, when the pressure drops to < 100 bar . This further raises a potential problem whether the D-Compress modelling can be extrapolated at lower pressure (i.e. a few bars) or not. If this problem does exist, the error bars would possibly be larger than the differences in predicted burial depths for different nakhlites. For example, errors of ± 5 bars or more means that differences in burial depth between NWA 998/Lafayette (~ 2-4 bars) and Y-000593/Nakhla (~ 0.5 bars) would be the same within error. In this case, however, the relative differences between pressure is still sound mathematically. That is, the NWA 998/Lafayette and Y-000593/Nakhla samples were buried at similar, shallow depths and that the former were likely buried deeper than the latter. The fact that the burial depths reported here agree with previous studies suggest that the new values are something close to correct. In other words, it should be clear that the calculated pressure and depth are model-derived and the accuracy of these estimations is model-dependant. A better calibration would yield better model results. We noted that, in current stage, the estimations are semiquantitative rather than perfect quantitative.

Concluding remarks
Based on the newly analysed data of NWA 5790 and the modelling of sulphur degassing, the following conclusions can be reached: (1) The titanomagnetite-ilmenite intergrowth in nakhlite NWA 5790 in this study extended the magma temperature record to 1032 °C and fO 2 record to ΔQFM + 1.60, respectively. (2) The burial depths obtained from the fO 2 -controlled sulphur degassing pressure are consistent with the previously estimated burial depths of various nakhlites, about 20-40 m for NWA 998 and Lafayette and < 7 m for Y-000593 and Nakhla. (3) Although the D-Compress model may be not calibrated sufficiently to resolve differences at low pressure of a few bars, our results suggest that analysis of fO 2 coupled with modelling of sulphur degassing pressure may constitute a method for estimating the burial depth of sulphur-rich lava flows on Mars.