A thermally conductive Martian core and implications for its dynamo cessation

Thermal conductivity at high pressure and room temperature

We used time-domain thermoreflectance (TDTR) coupled with diamond anvil cell (DAC) (2, 31, 32) (see Materials and Methods) to precisely measure the thermal conductivity of polycrystalline Fe₃S and FeS to ~40 GPa at room temperature. (Although measurements on these materials in liquid phase at high P-T are more relevant to Mars’ liquid core, now, they are technically challenging; we thus performed measurements to provide reliable thermal conductivity values of these solid Fe-S alloys at high P-T and then extrapolated the results into their liquid phase by considering the melting effect on conductivity, see below.) We find that the thermal conductivity of Fe₃S (~16 wt % S) at ambient conditions is ~6.8 W m⁻¹ K⁻¹ (Fig. 1A), which is more than an order of magnitude smaller than that of pure Fe (~76 W m⁻¹ K⁻¹) (2), indicating a strong reduction of thermal conductivity caused by the sulfur impurity. Application of pressure monotonically enhances the thermal conductivity, which reaches ~25 W m⁻¹ K⁻¹ at 40.5 GPa, the pressure at the center of the Mars. Moreover, with ~36 wt % S impurity, the FeS shows drastically distinct values of thermal conductivity and pressure dependence (Fig. 1B). The thermal conductivity of FeS at ambient conditions (FeS I phase) is even lower than that of Fe₃S, down to ~2.9 W m⁻¹ K⁻¹. Upon compression across ~3.4 GPa, the thermal conductivity suddenly increases to 9 W m⁻¹ K⁻¹ at 5.2 GPa (FeS II phase). After the transition at ~6.7 GPa to FeS III phase, considering the data uncertainty, the thermal conductivity remains approximately at 7.5 W m⁻¹ K⁻¹ until ~25 GPa, after which it significantly increases to 12 W m⁻¹ K⁻¹ at 43 GPa, half of the Fe₃S at similar pressures.

It is worth noting that, similar to the effect of silicon alloying (2), the presence of ~16 wt % (Fe₃S) and ~36 wt % S (FeS) in Fe significantly reduces the thermal conductivity, presumably due to the strong inelastic scattering between electrons and impurities (33–35). The impurity effects also alter the pressure dependences of Fe₃S and FeS, which are considerably different from that of pure Fe, where a concave pressure dependence with a minimum at ~40 GPa was observed (2).

Thermal conductivity at high pressure-temperature conditions

To quantify the effects of temperature and sulfur alloying, we further performed high P-T thermal conductivity measurements on the Fe₃S and FeS in externally heated DACs (EHDACs; Fig. 2). The thermal conductivities of Fe₃S and FeS both increase with temperature, similar, again, to the behaviors of Fe-Si alloys (2, 36). Note that along the P-T path of our measurements (i.e., first compressing the sample to a fixed pressure and then raising the temperature), the Fe₃S stayed in the tetragonal structure (37), while the FeS transitioned from FeS III to IV phase at ~730 K (38). Although our measurement temperature was only up to ~1023 K, we observed clear temperature dependences that change with the pressure and S content. We should note that the estimated Martian core temperature was ~2000 to 2400 K (39). If we assume that the thermal conductivity can be phenomenologically modeled as Λ(T) = aTⁿ, where a is a normalization constant, the exponent value n can be determined by the linear slope in the lnΛ-lnT plot. For Fe₃S, n = 0.25 (±0.07) at 15.4 GPa and slightly decreases with pressure to n = 0.2 (±0.04) at 35.5 GPa (Fig. 2A). With higher content of S impurity that results in stronger carrier scattering, the FeS IV phase has a stronger temperature dependence than the Fe₃S at similar pressures: n = 0.36 (±0.07) at 27 GPa and n = 0.33 (±0.05) at 28 GPa, respectively (Fig. 2B). Although the crystal structure of Fe₃S and FeS is different, our results suggest that at a given pressure, the higher content of impurity would have a stronger temperature dependence, which has also been reported in Fe-Si alloys (2, 36) (see Discussion).

Temperature dependence of thermal conductivity

The WF law relates the electronic thermal conductivity Λ_e of a material with its electrical resistivity ρ by Λ_e = L × T/ρ, where L is the Lorenz number that may vary with pressure and temperature, and T is the absolute temperature. Typically, when the temperature is comparable or higher than the Debye temperature of the material, its electrical resistivity can be expressed approximately as ρ = ρ₀ + α + βT, where ρ₀ is the residual resistivity arisen from defect scattering, α is a phenomenological constant associated with the shift in the resistivity due to temperature effect, and β is a temperature coefficient (40, 41). Therefore, Λ_e can be rewritten as Λ_e ~ L × T/(ρ₀ + α + βT). Since our measurement temperature (300 to 1023 K) is higher than one-third of the Debye temperature of Fe-S alloys (~400 K) (42), the high contents of S impurity in Fe makes the sum of (ρ₀ + α) larger than βT (41). The Λ_e of F-S alloys is thus expected to increase with increasing temperature. Furthermore, with increasing S content, the relative contribution of βT to the total electrical resistivity would be further reduced, leading to a stronger temperature dependence. On the other hand, for both Fe₃S and FeS, the temperature exponent n decreases with pressure, in line with the general trend observed for the Fe_0.85Si_0.15 (2). This is presumably caused by the fact that pressure reduces the effective density of impurity in momentum space (less impurity scattering) and enhances the relative contribution of βT, leading to a weaker temperature dependence of thermal conductivity (40).

Thermal conductivity profiles in Martian core and mantle

The recent geophysical observations by the InSight mission (21) suggest that the Martian core is larger than previously thought and confirm that, by contrast to the Earth, the Martian mantle is too thin to allow a transition from ringwoodite to bridgmanite, i.e., ringwoodite is the major mineral in its deep mantle. The mean Martian core density derived from InSight observations spans from ~5.7 to 6.3 g cm⁻³, where a number of potential combinations of S, O, and H contents alloyed in an Fe-Ni core is discussed (21). Geochemically defensible amounts of these light elements with a mean core density >6 g cm⁻³ are at the upper end of such density range. With a lower bulk mantle FeO content, the Insight observations preliminarily suggests that the core is majorly composed of Fe-Ni along with ~10 to 15 wt % S, <5 wt % O, and <1 wt % H and C (21). To understand the thermal evolution history of Mars and its early dynamo, we first model the thermal conductivity of the Martian core and mantle along an estimated Martian areotherm taken from (39) (an estimated present-day radial temperature profile within Mars), based on our present high P-T thermal conductivity of solid Fe-S alloys and previous data for ringwoodite (25). Note that if the core temperature were as high as 3000 K early in Mars’ history, then our modeled thermal conductivity of Martian core would be increased by only ~6%, which is within our data uncertainty. In other words, the temporal changes in the Martian areotherm as it cools down since its early stage have minor effects on our modeled thermal conductivity and numerical simulations on its thermal evolution. For the solid Fe₃S and FeS, we assume their temperature dependences of thermal conductivity at Martian core P-T conditions follow those shown in Fig. 2, i.e., scaling with T^0.21 and T^0.33, respectively. The thermal conductivity of ringwoodite, on the other hand, is assumed to follow a T^−0.5 dependence, typical of Fe-bearing minerals (31, 43–45). The modeled thermal conductivities of ringwoodite (blue curves), solid Fe₃S (black curve), and solid FeS (orange curve) along a Martian areotherm are plotted in Fig. 3.

Literature experimental investigations on the high P-T electrical resistivity ρ of solid Fe-S alloys have been largely limited to ≦10 GPa (28–30). For instance, ρ of Fe-20 wt % S at 4.5 GPa was found to be ~400 microhm cm up to ~1400 K (28), and ρ of FeS at 5 GPa (30) and 8 GPa (28) remained at ~400 and ~700 microhm cm, respectively, up to ~1700 K. Their corresponding thermal conductivities were then further inferred to be ~6 to 11 W m⁻¹ K⁻¹ via the WF law with an ideal Lorenz number. In addition, on the basis of high pressure but room temperature electrical resistivity measurements, ρ of solid Fe-14.2 wt % S at Martian core conditions was estimated to ~100 microhm cm, suggesting a thermal conductivity of 46 W m⁻¹ K⁻¹ at the Martian CMB and 62 W m⁻¹ K⁻¹ at the center (46). Using first-principles calculations, ρ of liquid Fe₃S and Fe₇S were determined to stay at ~107 and ~88 microhm cm, respectively, throughout the Martian core conditions (27). Again, it is important to note that our present study represents the first direct high P-T thermal conductivity measurements of solid Fe-S alloys that do not rely on the WF law and assumption of ideal Lorenz number. Our results for the thermal conductivities of solid Fe-S alloys (solid black and orange curves in Fig. 3) under similar P-T conditions and S content are much higher than those inferred from high P-T experimental electrical resistivity data (28), while only about half of those estimated from room temperature resistivity data (46) and calculated by first-principles (for liquid phase) (27). Table 1 summarizes thermal conductivity of Fe-S alloys at Martian CMB conditions. Here, we only include literature results available at CMB P-T conditions with an S content similar to our present study, as the chemical composition could significantly influence the thermal conductivity of Fe alloys.

Moreover, prior studies have suggested that the effect of melting would reduce the thermal conductivity of Fe and Fe alloys by ~20% or less [see, e.g., (2) and references therein]. Assuming that S is the major light element with ~16 wt % (~25 atomic %) (24) in the Martian core and a 15% conductivity reduction upon melting, our present data suggest that the thermal conductivity of liquid Fe₃S (red dashed curve in Fig. 3) would be ~19 W m⁻¹ K⁻¹ at the top of the Martian core and that it increases by ~70% to ~32 W m⁻¹ K⁻¹ at the center, leading to a much more thermally conductive Martian core than previously inferred. In addition, if the mineralogy of the Martian mantle is similar to that of the Earth’s upper mantle, which is predominantly composed of olivine and its high-pressure polymorph ringwoodite, then the thermal conductivity of ringwoodite in dry phase at the bottom of the mantle would be ~4.5 W m⁻¹ K⁻¹ (blue solid curve in Fig. 3). Note that if the ringwoodite in the Martian mantle contains substantial amounts of water, its thermal conductivity would be lower. For instance, the presence of 1.73 wt % water would lead to a conductivity of ~3 W m⁻¹ K⁻¹ at the bottom of the Martian mantle (see the blue dashed curve). Furthermore, if the ringwoodite in the Martian mantle contains more iron than that in the Earth, then its thermal conductivity would be even lower. Overall, a ~4- to 6-fold thermal conductivity discontinuity is expected across the Martian CMB (liquid Fe₃S versus ringwoodite), which would suppress heat released from the core and promote thermal stratification. Such findings yield distinct behaviors concerning the Q_c (see geodynamic modeling below). In addition, our results for Fe₃S and FeS at Martian core’s P-T conditions provide a platform to model the thermal conductivity of the Fe-S alloys with different S contents, since under different model assumptions the not-well-constrained S content in the Martian core may vary over the range we explored (~16 to 36 wt %) during Mars history. If the Martian core also contains few amounts of other elements (e.g., O, H, C, etc.), then these impurities are expected to further slightly reduce the thermal conductivity of Martian core, while their exact effects require future experimental and computational studies.

Implications for Martian thermal evolution and dynamo cessation

Dynamo operation and duration strongly depend on core thermal conductivity. Higher thermal conductivity promotes faster dynamo decay through ohmic dissipation, and in the absence of gravitational energy released from core crystallization, the dynamo may rapidly switch off. Recent simulations of Mars core-mantle coupled evolution (1) indicate that the cessation time of the Martian magnetic field (11–13) can be explained with a core thermal conductivity in the range of 16 to 35 W m⁻¹ K⁻¹, provided that the mantle reference viscosity, η₀, is between 10¹⁹ and 10²¹ Pa s. Modeling details further show that the pressure dependence of mantle viscosity (with an activation volume fixed to 6 cm³ mol⁻¹) requires a lower range of reference viscosity (10¹⁹ to 10²⁰ Pa s) and thermal conductivity (<20 W m⁻¹ K⁻¹) to explain this cessation time.

Here, we investigate how the radial profile of the core thermal conductivity, Λ_c (red dashed curve in Fig. 3), affects Mars thermal evolution. For this, we performed simulations coupling parameterized models of core and mantle evolution following the method developed in (1) (also see Materials and Methods), with the exception that the thermal conductivity is now allowed to vary with depth. Figure 4 plots evolution of key parameters [(A) heat flow Q_c, (B) entropy arising from ohmic dissipation, E_j, (C) bottom radius of the thermally stratified layer, r_s, and (D) temperature at the CMB] for three different scenarios, including the reference case in (1) (green curves). Thermal stratification is a consequence of a sub-isentropic heat flow (i.e., Q_c is lower than the heat flow along an isentropic temperature gradient, Q_a) and results in higher core temperatures. All cases were performed with the “standard” configuration of (1), i.e., the mantle viscosity does not depend on pressure (activation volume is set to zero), and the mantle abundance in heat-producing elements is that of (47). In addition, the initial mantle temperature and temperature jump at the CMB were set to 2327 and 182 K, respectively, corresponding to the standard reference case in (1) and leading to an initial CMB temperature of 2509 K. Dynamo switches off when E_j is below a small value E_j,min, corresponding to the minimum ohmic dissipation needed to self-sustain feedbacks between magnetohydrodynamic processes. For Earth, and in the case of a poloidal field, this value was estimated to 1.0 MW K⁻¹ (48), which we adopted in our estimation of the Martian dynamo cessation time. Greenwood et al. (1) pointed out that because E_j,min is much smaller than the dissipation provided by secular cooling, its exact value does not substantially influence the dynamo cessation time. The condition E_j < 1.0 MW K⁻¹ occurs around the transition from super- to sub-isentropic heat flow (i.e., Q_c being slightly higher or lower than Q_a). For the reference case (Λ_c is a constant of 24 W m⁻¹ K⁻¹ throughout the core and η₀ = 2.5 × 10²⁰ Pa s; green curves in Fig. 4), this happens ~0.9 Gyr, well within the estimated cessation range (11–13). For the depth-dependent profile of Λ_c based on our experimental data and still with η₀ = 2.5 × 10²⁰ Pa s, dynamo cessation occurs later, ~1.1 Gyr (blue curves in Fig. 4). Within error bars (blue shaded area), however, the cessation time expected from our Λ_c profile agrees with the observed range. In addition, good agreement can be obtained by assuming a slightly higher mantle viscosity (red curves in Fig. 4, the same Λ_c profile as the blue curves, but with η₀ = 3.5 × 10²⁰ Pa s). Thermal stratification starts when Q_c becomes lower than Q_a. Our calculations indicate that Q_a is lower with depth-increasing conductivity, which delays the onset of stratification. For the reference case, stratification starts ~0.3 Gyr before dynamo cessation, while it starts ~0.4 to 0.6 Gyr after dynamo cessation when the depth-dependent Λ_c is accounted for. Moreover, core thermal stratification is much faster with the depth-dependent Λ_c. For instance, the entire core is thermally stratified by 2.8 Gyr with our Λ_c profile and η₀ = 2.5 × 10²⁰ Pa s (blue curve in Fig. 4C), compared to ~4.5 Gyr for the reference case (green curve). A possible explanation for this faster stratification is that, compared to the constant Λ_c case, because of the depth-increasing Λ_c, the increase in cooling rate of the convecting core is greater than that of the stable layer, meaning that the thermal stratification completes earlier.

One aspect that is not accounted for in our simulations is the presence of a layer of molten silicate rocks at the base of the mantle, as recently inferred by InSight seismic data (19, 20). This layer may result from the crystallization of Martian magma ocean and is possibly enriched in iron and heat-producing elements at the base of the mantle (49). Because it is partially molten and enriched in iron, one would further expect the thermal conductivity of this layer to be much lower than that of ringwoodite (blue curve in Fig. 3). In other words, the mantle basal layer may act as a thermal blanket, opposing core cooling and dynamo action. Combined with our present core thermal conductivity, this would, in turn, request a lower mantle viscosity to explain the observed dynamo cessation time. Future direct measurements and computational studies on the thermal conductivity of other Fe-S–light element systems under relevant P-T conditions are required to better constrain the thermal evolution and energy budget of the Mars.

Starting materials and sample preparation

The Fe₃S polycrystals were synthesized by a Kawai apparatus with a 5000-tonne press (run number 5K3187). The starting material is a powder mixture of Fe and S with an atomic ratio of 3:1. The sample was compressed to 22 GPa and first heated to 1000°C, followed by a cooling process from 1000° to 855°C for 9.5 hours and lastly kept at 700°C for 2 hours. The starting FeS troilite powder was extracted from the Cape York IIIAB iron meteorite, also known as the Innaanganeq meteorite, in the Geological Museum of the University of Copenhagen, Denmark. Troilite in iron meteorites is known to be very chemically stoichiometric FeS, which is more suitable for the present study than commercial pyrrhotite samples (Fe_1-XS) with iron deficiency (X = 0 to 0.2). The chemical composition of Fe₃S and FeS was characterized to be Fe_2.97S (with a small amount of Fe₃S₂) and Fe_0.995S, respectively, by electron probe microanalyzer in Academia Sinica.

To prepare samples for thermal conductivity measurements at high pressure and room temperature, each sample was firstly hand-polished down to ~10 μm thick and then thermally evaporated with a 90-nm-thick Al film. The sample along with several ruby spheres was loaded into a symmetric DAC where a pair of 300-μm culet anvils and a Re gasket were used. Because of its relatively low thermal conductivity at high pressures, silicone oil (Chemical Abstracts Service no. 63148-62-9 from Acros Organics) was loaded as the pressure transmitting medium, which remains in a reasonably well quasi-hydrostatic condition over the pressure range we studied (<45 GPa). The pressure within the sample chamber of the DAC was monitored by fluorescence (50) and Raman of the ruby with a typical uncertainty of <5%.

In our high P-T measurements, the ~10-μm-thick sample and ruby spheres were loaded into an EHDAC. Here, polycrystalline NaCl powder (thermally dried at ~120°C for 1 hour before being loaded) was used as the pressure medium. The EHDAC provides a spatially homogeneous temperature distribution over the sample chamber. Since the pressure within the sample chamber may change upon heating, a gas membrane that enables in situ control on the pressure within the EHDAC was also used, allowing us to keep the pressure as close as possible (within ~1 GPa variation range) to our target pressure. The combination of these techniques allows us to investigate the temperature dependence of thermal conductivity of Fe₃S and FeS at a given pressure up to ~1023 K. Details of the EHDAC assemblage and pressure-temperature measurement, as well as the sample geometry and experimental setup, can be found in (43, 51). Note that as expected from previous reports, our measured ruby R1 fluorescence peak does become weak and broad at temperatures higher than ~700 K. At a given high temperature, we determined the peak position with a typical wavelength uncertainty of ~0.05 to 0.2 nm, which could translate a pressure uncertainty of ~0.1 to 0.6 GPa using the temperature correction formula in (52). In addition, the uncertainty of temperature measurement using an R-type thermocouple is only about few Kelvin, which could induce an additional pressure uncertainty of ~0.05 to 0.1 GPa. Thus, in our study, the overall pressure uncertainty at either room temperature or high temperature is typically less than 2 GPa. In our data analysis (Figs. 1 and 2 and fig. S2), we have already taken into account the effect of such pressure uncertainty.

Thermal conductivity measurements

We used TDTR to measure the thermal conductivity of Fe₃S and FeS over a wide range of P-T conditions presented in this work. TDTR is a widely used thermal metrology method that provides thermal conductivity measurements of various materials with high precision. In the past decade, it has also been successfully coupled with high P-T techniques, see, e.g., (2, 31, 43, 53). In short, TDTR is an ultrafast optical pump-probe spectroscopy, which uses a split pump pulse to heat up the Al film coated on the sample and a split probe pulse to detect the heat diffusion dynamics through the sample. By comparing the temporal evolution of the Al’s reflectivity change with numerical calculations based on a bidirectional thermal model, the thermal conductivity of the sample of interest is determined. More details of the principle, experimental setup, and data analysis of the TDTR are described in literatures, e.g., (2, 51, 53–55) and references therein.

Figure S1 shows a set of representative TDTR spectrum for Fe₃S at 40.5 GPa and room temperature along with data fitting by the thermal model calculations. In the thermal model, the volumetric heat capacity of the sample of interest (Fe-S alloys) is an important parameter. At ambient conditions, the volumetric heat capacity of FeS is 2.78 J cm⁻³ K⁻¹ taken from (56) while that of Fe₃S is estimated to be 3.16 J cm⁻³ K⁻¹ by linear interpolation between pure Fe (3.54 J cm⁻³ K⁻¹) and FeS. The heat capacity of both FeS and Fe₃S at high P-T conditions is not known and thus assumed to be a constant as their ambient values. We note that the uncertainty in our data majorly arises from the analysis uncertainty rather than the measurement uncertainty. On the basis of the method (57, 58) described in figs. S1 and S2, the uncertainties in all the parameters in our thermal model calculations would translate <10% error in the derived thermal conductivity of Fe-S alloys before 20 GPa and ~10 to 17% error at 20 to 40 GPa.

Modeling of core evolution

To investigate the evolution of the Martian core, including the growth of a stable layer, and to estimate the cessation time of Martian dynamo, we used the approach in (1), which is coupling models of parameterized convection for the core and for the mantle. For the mantle, this approach uses a parameterization based on stagnant lid thermal convection (59, 60) with a simplified lithosphere. The mantle is heated both from below (heat extracted from the core) and from within (radiogenic heating due to the decay of U, Th, and K). Concentration in heat-producing elements is taken from the compositional model of (47). The convective interior is assumed isothermal, and temperature increases linearly with depth within thermal boundary layers. The crust thickness is considered as constant over time, and melting and crust growth are neglected. The thickness of the stagnant lid (modeling the lithosphere) is set to 300 km. Greenwood et al. (1) showed that the thickness of the stagnant lid does not have a strong impact on the evolution of the mantle convective interior and of the core. Mantle viscosity depends on temperature following an Arrhenius law with activation energy fixed to 300 kJ mol⁻¹ and is specified (reference viscosity η₀) at a temperature of 1600 K. Viscosity is further allowed to increase with temperature, but, for simplicity, here is neglected this dependence (activation volume V_a is set to 0). The time for dynamo’s cessation strongly depends on η₀ and, to a lesser extent, on activation volume (1). For core thermal conductivity in the range of 16 to 35 W m⁻¹ K⁻¹, calculated cessation time agrees with that estimated for Martian dynamo with viscosity in the range of 10²⁰ to 10²¹ Pa s. Accounting for viscosity pressure dependence (V_a = 6 cm³ mol⁻¹) reduces this range by one order of magnitude. In calculations, we tested values of the reference viscosity, η₀ = 2.5 × 10²⁰ Pa s and η₀ = 3.5 × 10²⁰ Pa s.

The core is initially well mixed with an adiabatic temperature profile and is animated with convection. The core is initially assumed super-heated compared to the mantle with a temperature difference dT, which we fixed to 182 K, based on the standard model from (1), and the temperature at the CMB is readjusted according to this treatment. The heat flow at CMB, Q_c, is calculated from scaling relationships built for mantle convection (60) and is further depending on the temperature at the CMB, T_CMB, which is determined from the core evolution. As Q_c becomes lower than the heat flow calculated along an isentropic temperature profile, Q_a, a stable conductive layer starts to grow, and core convection is confined beneath this layer. The evolution of the stable layer follows the treatment in (61). The ohmic dissipation related to magnetic field is measured with the entropy balance E_j between the entropies originating from secular cooling and thermal conduction, E_s and E_k, which are calculated following (61). Because E_k is proportional to the core thermal conductivity, Λ_c, larger conductivity leads to larger E_k and thus lower E_j. In other words, dynamo action is reduced with increasing Λ_c and may stop if conduction is too large.

This PDF file includes: