Sanloup, C., & Fei, Y. (2004). Closure of the Fe-S-Si liquid miscibility gap at high pressure. *Physics of the Earth and Planetary Interiors, **147**(1), 57–65. https://doi.org/10.1016/j.pepi.2004.06.008

Schubert, G., Chan, K. H., Liao, X., & Zhang, K. (2004). Planetary dynamos: Effects of electrically conducting flows overlying turbulent regions of magnetic field generation. *Icarus*, **172**(2), 305–315. https://doi.org/10.1016/j.icarus.2004.06.007

Schubert, G., Ross, M., Stevenson, D. J., & Spohn, T. (1988). Mercury's thermal history and the generation of its magnetic field. In F. Vilas, C. R. Chapman, & M. Matthews (Eds.), *Mercury*, (pp. 429–460). Tucson: The University of Arizona Press.

Secco, R. A. (2017). Thermal conductivity and Seebeck coefficient of Fe and Fe-Si alloys: Implications for variable Lorenz number. *Physics of the Earth and Planetary Interiors*, **265**, 23–34. https://doi.org/10.1016/j.pepi.2017.01.005

Slack, G. A. (1965). Thermal conductivity of elements with complex lattices: B, P, S. *Physical Review*, **139**(2A), A507–A515. https://doi.org/10.1103/PhysRev.139.A507

Smith, C. S., & Palmer, E. W. (1935). Thermal and electrical conductivities of copper alloys. *Transactions of the American Institute of Mining and Metallurgical Engineers*, **117**, 225–243.

Smith, D. E., Zuber, M. T., Phillips, R. J., Solomon, S. C., Hauck, S. A., Lemoine, F. G., et al. (2012). Gravity field and internal structure of Mercury from MESSENGER. *Science*, **336**(6078), 214–217. https://doi.org/10.1126/science.1218809

Stanley, S., Bloxham, J., Hutchison, W. E., & Zuber, M. T. (2005). Thin shell dynamo models consistent with Mercury's weak observed magnetic field. *Earth and Planetary Science Letters*, **234**(1–2), 27–38. https://doi.org/10.1016/j.epsl.2005.02.040

Stevenson, D. J. (1987). Mercury's magnetic field: A thermoelectric dynamo? *Earth and Planetary Science Letters*, **82**(1–2), 114–120. https://doi.org/10.1016/0012-821X(87)90111-7

Stevenson, D. J. (2003). Planetary magnetic fields. *Earth and Planetary Science Letters*, **6**(3), 213–214. https://doi.org/10.1029/GL006i003p00213

Stevenson, D. J., Spohn, T., & Schubert, G. (1983). Magnetism and thermal evolution of the terrestrial planets. *Icarus*, **54**(3), 466–489. https://doi.org/10.1016/0019-1035(83)90241-5

Tian, Z., Zuber, M. T., & Stanley, S. (2015). Magnetic field modeling for Mercury using dynamo models with a stable layer and laterally variable heat flux. *Icarus*, **260**, 263–268. https://doi.org/10.1016/j.icarus.2015.07.019

Tosi, N., Grott, M., Plesa, A. C., & Breuer, D. (2013). Thermochemical evolution of Mercury's interior. *Journal of Geophysical Research: Planets*, **118**, 2474–2487. https://doi.org/10.1002/jgre.20168

Urakawa, S., Someya, K., Terasaki, H., Katsura, T., Yokoshi, S., Funakoshi, K. I., et al. (2004). Phase relationships and equations of state for FeS at high pressures temperatures and implications for the internal structure of Mars. *Physics of the Earth and Planetary Interiors*, **143-144**(1–2), 469–479. https://doi.org/10.1016/j.pepi.2003.12.015

Veasey, M., & Dumberry, M. (2011). The influence of Mercury's inner core on its physical libration. *Icarus*, **214**(1), 265–274. https://doi.org/10.1016/j.icarus.2011.04.025

Wang, Y., Uchida, T., von Dreele, R., Rivers, M. L., Nishiyama, N., Funakoshi, K. I., et al. (2004). A new technique for angle-dispersive powder diffraction using an energy-dispersive setup and synchrotron radiation. *Journal of Applied Crystallography*, **37**(6), 947–956. https://doi.org/10.1107/S0021889804022502

Wicht, J., Mandea, M., Takahashi, F., Christensen, U. R., Matsushima, M., & Langlais, B. (2007). The origin of Mercury's internal magnetic field. *Space Science Reviews*, **132**(2–4), 261–290. https://doi.org/10.1007/s11214-007-9280-5
<figure><img src="image_258.png" /></figure>

# Thermal conductivities of solid and molten silicates: Implications for dynamos in mercury-like proto-planets

D. Freitas<sup>a,b,*</sup>, J. Monteux<sup>a</sup>, D. Andrault<sup>a</sup>, G. Manthilake<sup>a</sup>, A. Mathieu<sup>a</sup>, F. Schiavi<sup>a</sup>, N. Cluzel<sup>a</sup>

## ARTICLE INFO

**Keywords:**

Thermal diffusivity

Thermal conductivity

Melts

Geodynamo

Mercury

## ABSTRACT

Remanent magnetization and active magnetic fields have been detected for several telluric planetary bodies in the solar system (Earth, Mercury, Moon, Mars) suggesting the presence of core dynamos active at the early stages of the planet formation and variable lifetimes. Among the factors controlling the possibility of core dynamos generation, the dynamics of the surrounding silicate mantle and its associated thermal properties are crucial. The mantle governs the heat evacuation from the core and as a consequence the likeliness of an early thermally driven dynamo. In the case of planets with a thick mantle (associated with supercritical Rayleigh numbers), the core heat is efficiently removed by mantle convection and early thermally-driven dynamos are likely. At the opposite, planets with a thin mantle (associated with subcritical Rayleigh numbers) might evacuate their inner heat by diffusion only, making early thermally-driven dynamos difficult. Within the Solar System, Mercury is a potential example of such a regime. Its small mantle thickness over the planet radius ratio might be inherent to its small orbital semi-axis and hence, might be ubiquitous among the terrestrial objects formed close to their star.

To constrain the likeliness of a thermally driven dynamo on "Mercury-like" planets (i.e. with large Rc/R), we present new thermal diffusivity measurements of various solid, glassy and molten samples. We applied the Angstrom method on cylindrical samples during multi-anvil apparatus experiments at pressures of 2 GPa and temperatures up to 1700 K. Thermal diffusivities and conductivities were estimated for solid and partially molten peridotites, with various melt fractions, and for basaltic and rhyolitic glasses and melts. Our study demonstrates that melts have similar thermal properties despite a broad range of composition investigated. The melts reveal much lower thermal conductivities than the solids with almost an order of magnitude of decrease: 1.70 (±0.19) to 2.29 (±0.26) W/m/K against 0.18 (±0.01) to 0.41 (±0.03) W/m/K for peridotites at high temperatures and various melts respectively. Partially molten samples lie in between and several predictive laws are proposed as a function of the melt fraction and solid/melt texture.

Using our results into forward calculations of heat fluxes for dynamo generation for Mercury-like planets, we quantify the effect of mantle melting on the occurrence of thermally driven dynamos. The presence of a mushy mantle and partial melting could significantly reduce the ability of the mantle to evacuate the heat from the core and can prevent, shut or affect the presence of a planetary magnetic field. The buoyancy and fate of molten material in such bodies can thus influence the magnetic history of the planet. Future observations of Mercury-like planets accreted near their star and the detections of their magnetic signatures could provide constraints on their inner state and partial melting histories.

## 1. Introduction

The presence of internally generated magnetic fields is a variable feature among telluric bodies in the inner solar system. While few are currently active, such as for Mercury and Earth, several are now extinct,

as observed for the Moon, Venus and Mars. On terrestrial planets currently exhibiting a dynamo, the generated magnetic fields characteristics are very different. Earth's magnetic field is very intense 25–65 µT, originating from the core and present for at least 3.5 Ga (Tarduno et al., 2010). Mercury's field strength is much weaker, representing

<footnote><sup>a</sup> Université Clermont Auvergne, CNRS, IRD, OPGC, Laboratoire Magmas et Volcans, F-6300 Clermont-Ferrand, France</footnote>

<footnote><sup>b</sup> School of Geosciences, Grant Institute, The King's Buildings, The University of Edinburgh, West Mains Road, Edinburgh EH9 3JW, United Kingdom</footnote>

<footnote>* Corresponding author at: School of Geosciences, Grant Institute, The King's Buildings, The University of Edinburgh, West Mains Road, Edinburgh EH9 3JW, United Kingdom.</footnote>

<footnote>E-mail address: damien.freitas@ed.ac.uk (D. Freitas).</footnote>
around a 1% of Earth's one (Kabin et al., 2008; Anderson et al., 2011) and its shape is also unique among the different detection in the solar system (Tian et al., 2015). According to the remnant magnetization measured in the crust, it was proposed that such weak magnetic activity occurred during the last 3.9 Ga (Johnson et al., 2015). The source of such a weak-and-prolonged dynamo is still largely debated (Manthilake et al., 2019). In the meantime, there are several evidences for an intrinsic dynamo during the early stages of both Mars and Moon (Acuna et al., 1999; Hood et al., 2010). Their dynamo seems to cease around 4.1–3.9 Ga ago for Mars (Johnson and Phillips, 2005; Lillis et al., 2008; Lillis et al., 2013) and exhibits a somewhat complicated history for the Moon with a strong dynamo between 4.25 and 3.5 Ga, followed by a weak persistence up to 2.5–2 Ga ago (Tikoo et al., 2014; Lawrence et al., 2008; Garrick-Bethell et al., 2009; Mighani et al., 2020). For Venus, the early presence of a dynamo remains yet unconstrained but may be detectable in future explorations (Nimmo, 2002; O'Rourke et al., 2019). More broadly, evidences of paleomagnetic anomalies indicate that the angrite parent bodies, originating from inner regions of the solar system, were subject to an early internally generated dynamo (Weiss et al., 2008). All these elements suggest that transient dynamos might be a somewhat common feature in telluric bodies (Monteux et al., 2011).

The presence of early dynamos is highly conditioned by the internal structure of the planet and its capacity to release the heat accumulated during the accretion processes (accretion, metal/silicate differentiation, core and mantle crystallization) as well as short-lived radiogenic heating. As heat conduction is an inefficient heat transport process in silicates at high temperatures (Hofmeister and Branlund, 2015), the onset of a mantle global convection is a crucial step in planet's thermal history. Convection starts when the Rayleigh number ($Ra$) of the terrestrial mantle is larger than the critical Rayleigh number ($Ra_c$). The higher is $Ra$, the stronger the convection, and the more efficient the heat transport. In contrast, heat is only transported by conduction for a terrestrial mantle with a $Ra < Ra_c$. As $Ra$ scales with mantle thickness ($h_{mantle}$) as $h_{mantle}^3$, convection should take place easily in planets with a thick silicate mantle, even after the solidification of the early magma ocean stage. Consequences are efficient evacuation of the inner heat and the possible occurrence of a dynamo. At the opposite, bodies with a thinner mantle lead to smaller $Ra$ values, making convection unlikely and dynamos more difficult to generate.

Mercury is the most interesting planet for our study. Indeed, its mantle is thin 420 ± 30 km (Hauck et al., 2013) and its core occupies almost 55% of planet's volume and 65% of the planet mass (Strom and Sprague, 2003; Charlier and Namur, 2019). Hence, the mantle of Mercury is controversially at the limit between conductive and convective regimes (Breuer et al., 2007). Different scenarios could explain the small $h_{mantle}$/planet radius ($R$) ratio on Mercury (Charlier and Namur, 2019): (1) primordial nebular processes, yielding to the enrichment of metal over silicate materials in the inner solar system (Ebel and Grossman, 2000; Wurm et al., 2013; Weidenschilling, 1978), (2) highly energetic accretional collisions inducing a major loss of the silicate fraction (Benz et al., 1988; Asphaug and Reufer, 2014), and (3) post-accretion scenarios, with major vaporization of the volatile and silicate elements from the planet during magma ocean stage (Fegley and Cameron, 1987; Boujibar et al., 2015). If scenarios (1) and (3) are dominant, then “Mercury-like” planets with small $h_{mantle}/R$ ratio would be ubiquitous within all planetary systems. Moreover, such a small ratio would be prevailing during the whole accretionary processes. Mercury-like bodies could adopt a wide range of possible compositions depending on their history. For example, bodies accreted from reduced enstatite and/or carbonaceous bencubbinite chondrites (Malavergne et al., 2010) could present a mantle composition similar to terrestrial lherzolite, but with a sulfur content potentially as high as 11 wt% (Namur et al., 2016). Accordingly, the diversity of the silicate samples found on Earth in the forms of rocks, melts and glasses is a good proxy to decipher the properties of a range of mantle-relevant silicate compositions on Mercury-like bodies. Depending on the planet size and the $h_{mantle}/R$ ratio, the

internal pressure ranges from a few MPa to several GPa. We note that extensive mantle melting likely occurred at different stages of the history of such Mercury-like bodies. Major energy incomes are expected from the vicinity to the young Sun (2500 to 3500 K according Charlier and Namur, 2019), internal energy release (chemical and gravitational differentiation, core crystallization), and presence of short-period radioactive elements (Al²⁶, K⁴⁰ etc.). The very high temperatures likely induce extensive melting of the thin mantle up to the core-mantle boundary (CMB).

As mentioned above, dynamo generation could be difficult for Mercury-like planets, due to the subcritical value of $Ra$ possibly disabling mantle convection. In such case, heat transfer by conduction would dominate the planet history and the thermal conductivity of the silicate mantle is a key parameter governing the early core heat flow. Silicates thermal conduction properties are now well characterized at ambient conditions (Hofmeister and Branlund, 2015). Among them, the most common geological minerals were characterized in the forms of single crystal and polycrystalline aggregates: olivine (Osako et al., 2004; Xu et al., 2004; Pertermann and Hofmeister, 2006; Gibert et al., 2005), periclase (Hofmeister and Branlund, 2015), feldspar (Pertermann and Hofmeister, 2006; Hofmeister et al., 2009; Branlund and Hofmeister, 2012) and pyroxenes (Hofmeister et al., 2014, Hofmeister and Pertermann, 2008) as well as peridotite rocks (Gibert et al., 2005; Beck et al., 1978), which have been extensively studied due to their important geological implications. The thermal diffusivities of minerals are linked to the characteristics of their lattice structure and modes of phonon generation and propagation as a function of temperature. For silicates, lattice thermal diffusivities usually decrease with increasing temperature following a 1/T dependence. At the opposite, the diffusivities increase while increasing pressure, however, the temperature dependence is much greater than that of pressure over the considered ranges for small planets. The resulting implication is that mantle rocks and minerals are poor thermal conductors in planetary interiors. Recent measurements of thermal diffusivities of glasses and melts at ambient pressure suggested that the non-crystalline silicates are even more insulating than the minerals (Hofmeister et al., 2009, 2014; Romine et al., 2012).

Measurements of thermal diffusivities at relevant conditions of planetary mantles encompass important difficulties, which were overcome by the use of different techniques (see Hofmeister and Branlund, 2015 for a critical review). Measurements were reported over a wide P and T range for large volume samples of olivine, periclase, bridgmanite using Angström or Pulse method in solid pressure apparatus (Osako et al., 2004; Xu et al., 2004; Manthilake et al., 2011a and Manthilake et al., 2011b, Zhang et al., 2019). Up to now, accurate measurements of geologically relevant silicate glass, partially molten systems and melts at planetary interior conditions remain scarce, if not absent. The available results report a nearly flat evolution of thermal diffusivities with the temperature above 1000 K for silicates, glass and melts (Hofmeister et al., 2014; Hofmeister and Branlund, 2015), suggesting that the measurements above the melting temperature could be safely extrapolated to planetary P-T conditions.

In this study, we aim at better constraining the thermal properties of Mercury-like protoplanets where convection is unlikely and heat is mostly removed by diffusion. We perform HP-HT in situ thermal diffusivities measurements of solid, partially molten and fully molten silicate for various compositions in Multi-anvil apparatus and using Angström method. Then, we constrain the likeliness of a thermally driven dynamo during the early stages of the evolution of a Mercury-like planet. We consider a wide range of planet sizes ($1 \text{ km} < R < R_{\text{Mercury}}$) and thermal states (with solid and partially molten mantles) on the dynamo likeliness.
## 2. Experimental and analytical methods

### 2.1. High-pressure assemblies

High-pressure and high-temperature experiments were performed using a 1500-ton Kawai type Multi-anvil apparatus. All experiments were conducted at 2 GPa, based on a previous press-load vs sample pressure calibration (Boujibar et al., 2014), providing an uncertainty of ~0.1 GPa at pressures <5 GPa. We used octahedral pressure media with full length edges composed of MgO doped with Cr₂O₃ (5 wt%) in a 25/ 17 multi-anvil configuration (octahedron edge length / anvil truncation edge length) (Fig. 1). Our assembly was designed to accommodate the specific requirements for measurements of thermal conductivity of relatively large samples (4 mm long for 4–3.5 mm diameter). All ceramic parts of the cell assembly, including the pressure media, were fired at 1373 K prior assembling in order to remove the absorbed moisture. Oxygen fugacity of the sample was not controlled during the experiments but is expected to be quite reducing due to the presence of the graphite furnace. The use of a steeped graphite furnace helped reducing thermal gradients. Thermal loss from the sample zone was further reduced by the use of a thick zirconia (ZrO₂) sleeve around the furnace. Thermal gradients in our assembly were computed using the software developed by Hernlund et al., 2006. The models show temperature gradients within the sample volume limited to ~7.5 K/mm vertically and even less radially/horizontally (Fig. S1). On the other hand, uncertainties of our thermocouple reading (i.e. where the conductivity measurement is performed, see below) are less than 5 K on absolute and 0.1 K in relative temperatures.

Previous experimental studies described the difficulties to perform a good measurement of thermal diffusivity for a molten sample, due to a potential sample deformation (Hofmeister et al., 2009, 2014; Romine et al., 2012). However, deformation remained minor in our experiments, as evidenced by the good reproducibility of the measurements during the repeated cycles of heating and cooling (Figs. 1 and 2) as well as the shape of the recovered samples (Figs. S6, S7 and S8).

Two tungsten-phenonium (W₉₅Re₅–W₇₄Re₂₆) thermocouples of 75 µm of diameter were used to measure temperature oscillations at the center (in 300 and 600 µm drilled holes in solids and glasses, respectively) and at the edge of the cylindrical sample. Special care was given to ensure that the junction points of the two thermocouples were located in the same sample plane perpendicular to the cylinder axis. Measurements of

<figure><img src="image_259.png" /><figcaption>Fig. 1. Schematic cross-section of the high-pressure cell assembly used for Angstrom method with multi-anvil apparatus. For glass and melt experiments, thermocouples were inserted in alumina tubes of 0.6 mm diameter and 4.5 mm length in the sample zone.</figcaption></figure>

<figure><img src="image_260.png" /><figcaption>Fig. 2. Experimental procedure for the measurement of thermal diffusivity of peridotites, glasses and melts: black, blue and red for cycles 1, 2 and 3 respectively. Heating (H) and cooling (C) parts of the cycles are separated by black dashed lines (for example H1: 1st cycle heating). The expected temperatures for sample dehydration (up to ~500 °C), glass transition and melting are illustrated in shaded colors: green, pink and yellow respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

the thermal conductivity of melts and glasses have always been particularly challenging due to the risk of thermocouple short circuit, because melts, even dry, are good electrical conductors (Tyburczy and Waff, 1983; Ni et al., 2011). To prevent this effect, thermocouples were inserted in alumina tubes of 0.6 mm diameter and 170 µm wall thickness. While these tubes show good resistance to cold compression and almost no reaction with the samples, even at very high temperature, some leakages for low viscosity basaltic melts have been identified. Their occurrence was taken in account for thermal diffusivity estimations (see Supplementary Text S1 and Figs. S2 and S3).

### 2.2. Angström method for thermal diffusivity measurements

We aim at determining thermal diffusivities of geological samples such as peridotites, mafic and felsic glasses and their melts at high pressure and high temperature using a double contact method: the Angström method. Experimental configurations, assemblies and data treatment are similar to several previous studies (Fujisawa et al., 1968; Kanamori et al., 1969; Katsura, 1993; Xu et al., 2004; Manthilake et al., 2011a; Manthilake et al., 2011b). Briefly, a temperature wave is generated radially in a cylindrical sample by the surrounding heater sleeve of the multi-anvil assembly (Fig. 1). Oscillations are generated with controlled frequency and period by modulating the power supply. Periodic temperature signals are recorded by two thermocouples fixed in the center and on the edge of the sample cylinder. At frequencies higher than 1.5 Hz, the signal can get noisier due to a limited time resolution of the recording system.

The recorded signals, for each thermocouple channel, are fitted by a nonlinear least square solver (lsqcurvefit on Matlab© using Levenberg-Marquardt algorithm):

$$ \Delta T = A_0 + A_1 t + A_2 \sin \left( 2\pi A_3 t + \frac{A_4 \pi}{180} \right) \quad (1) $$

Following this method, we obtained the amplitude of the temperature variation (A₂), frequency (A₃), and phase (A₄) of the recorded wave as a function of time (t), for the two thermocouples. Errors are quantified
based on the residue on the non-linear curve fitting (the nlparci function in Matlab©). To infer thermal diffusivity or conductivity from these parameters, the equation of conductive heat transport has to be solved. Here we consider the sample as an infinite uniform cylinder and assume that the heat flow is negligible in the vertical direction, thanks to a relatively long cylindrical heater. The following equation, expressed in cylindrical coordinates by Carlsaw and Jaegger (1959), has to be inverted in order to retrieve the diffusivity of the sample:

$$ \frac{dT}{dt} = D \left( \frac{d^2 T}{dr^2} + \frac{1}{r} \frac{dT}{dr} \right) \qquad (2) $$

where r is the radial distance from the axis, T the temperature, t the time, and D the thermal diffusivity. The boundary condition of our setup is:

$$ \frac{dT}{dr} = 0, \text{ at } r = 0 \qquad (3) $$

If we consider harmonic excitation at a distance r = R from the axis in normal or complex form:

$$ T_R = B_0 + B_1 \cos wt \leftrightarrow T_R = b_0 + b_1 Re (\exp iwt) \qquad (4) $$

where B₀, B₁, b₀ and b₁ are constants and Re is the real part of the exponential. The solution of the radial flow equation with the boundary conditions developed above can be expressed as:

$$ T_R = b_0 + b_1 \left[ J_0 \left( \sqrt{-i} * x \right) * \frac{\exp(iwt)}{J_0 \left( \sqrt{-i} * l \right)} \right] \qquad (5) $$

where

$$ l = (w/\kappa)^{1/2} R \text{ and } x = (w/\kappa)^{1/2} r \qquad (6) $$

for 0 ≤ r ≤ R, where w is the angular frequency and J₀ is the Bessel function of the first kind (integer order n = 0). At r = 0 and r = R we have:

$$ T_0 = b_0 + b_1 \theta \cos(\omega t - \varphi) \qquad (7) $$

$$ T_R = b_0 + b_1 \cos wt \qquad (8) $$

where

$$ \theta = \frac{1}{\sqrt{\text{bei}(u)^2 + \text{ber}(u)^2}} \qquad (9) $$

$$ \varphi = \tan^{-1} \left( \frac{\text{bei}(u)}{\text{ber}(u)} \right) \qquad (10) $$

where θ is the amplitude ratio and φ the phase shift between the two harmonic temperature measurements and bei and ber are imaginary and real parts of the Bessel function of the first kind, respectively. The solution to the Eq. (4) can be written using the dimensionless argument u, from which thermal diffusivity (D) can be directly estimated knowing angular frequency (w) and sample radius (d).

$$ u = d(w/D)^{1/2} \qquad (11) $$

Thanks to the Eqs. (9), (10), and (11), the diffusivity was then retrieved via forward Monte Carlo simulation and neighborhood algorithm (Sambridge and Mosegaard, 2002). In this step, different values of diffusivities are generated and theoretical phase shifts and amplitude ratios are calculated. These values are compared to the values measured in our experiments. When the differences between the computed solution and the experimental determination tend to 0 (minimization step), the correct diffusivity is then obtained if the sample radius (d) and angular frequency (w) are known. We note that mathematical solutions appear every 360° for the phase shift. Such erroneous solutions are checked manually and discarded.

For most of our experiments, we observe a significant variation of the refined raw-diffusivity value as a function of the heat-source frequency. This effect was already reported in the literature. Different equations were proposed to refine a real value of diffusivity, corresponding to the infinite frequency asymptote, based on non-linear equations. While Manthilake et al. (2011b) used:

$$ D = D_{\infty} + A_0 \exp(A_1 * f) \qquad (12) $$

Xu et al. (2004) used:

$$ D = D_{\infty} + A_0 \exp \left( \frac{-f}{f_0} \right) \qquad (13) $$

where D is diffusivity, f source frequency, f₀ asymptote frequency and A₀ and A₁ constants. In these equations, D∞, Aₙ and f₀ are inverted parameters. For a better fit of our experimental data and to minimize the uncertainties on the parameters, we adopt Eq. (13). On the other hand, for experiments presenting no systematic dependence of the raw-diffusivity with frequency, we consider as real value the average value between the raw-diffusivity values measured at all frequencies (Xu et al., 2004; Manthilake et al., 2011b).

## 2.3. Experimental uncertainties

Experimental measurements of thermal diffusivities and their further transformation into thermal conductivities generate uncertainties originating from the estimations of pressure, temperature, sample dimensions and the data fitting itself. Experimental uncertainties on pressure and temperature are presented above (~0.1 GPa and ~5 K, respectively). Sample lengths prior to the sample loading (d₀) and after the melting experiments were determined with a high precision digital gauge (accuracy of ~1 µm) and using the Scanning Electron Microscope (FEG-SEM), respectively (see Table 1). Then, the sample radius during the experiment at high pressure and temperature was calculated using:

$$ d(P, T) = d_0 * \left( 1 - \alpha(T - 298) + \frac{P}{K} \right)^{-1/3} \qquad (14) $$

where α is the thermal expansion and K the bulk modulus of the sample. For these calculations only samples radii before experiments were considered, as post mortem measurements are affected by decompression cracks (highlighted by the larger values measured after the experiments in Table 1). Moreover, the distance between the two thermocouples could not be measured precisely for a few samples. The values of all experimental parameters are provided in Supplementary Materials (Text S4 and Fig. S4). Altogether, final uncertainty on the sample length is between 5 and 10 µm.

There are other uncertainties associated with the procedure of data fitting for the determination of raw thermal diffusivities. Experimental phase shifts and amplitude ratios (Eqs. (9) and (10)) are determined with a precision generally better than 1% and majored by 3% in the worst cases. In the course of the Monte Carlo simulation, the differences between experimental and theoretical diffusivities are recorded and used a posteriori to refine uncertainties within 1σ errors.

Then, real thermal diffusivities are refined from the raw-diffusivities using either (i) the asymptotic non-linear fit for experiments presenting a dependence in frequency (Eq. (13), which yields important uncertainties on the refined parameters) or by averaging (see Methods). The error on average (σAVG) is:

$$ \sigma_{\text{AVG}} = \sqrt{\left(\frac{1}{n}\right)^2 * \sum_{i=n}^{n} \sigma_D^2} \qquad (15) $$

where n is the number of diffusivity measurements performed at different frequencies and σD is the error on each diffusivity. Raw diffusivities recovered from Monte Carlo processing have errors of ~1%, similar to those of phase shift and amplitude ratio. If the fitting step is
<table><caption>**Table 1**
Sum up of our Ansgström experiments with compositions, P-T ranges and associated parameters.</caption><thead><tr><th>Run#</th><th>M662</th><th>M804</th><th>M807</th><th>M808</th><th>M836</th><th>M843</th><th>M844</th><th>M846</th><th>M847</th><th>M848</th><th>M850</th></tr></thead><tbody><tr><td>Sample name</td><td>EPR-MORB</td><td>RP4-S2</td><td>DR07-S2</td><td>GD-S1</td><td>Haplo S1</td><td>ATHO-S1</td><td>DR11-S1</td><td>DR11-S2</td><td>RP4-S3</td><td>ATHO-S2</td><td>GD-S2</td></tr><tr><td>Composition</td><td>MORB</td><td>Peridotite</td><td>MORB</td><td>Rhyolite</td><td>haplobasalt</td><td>Rhyolite</td><td>MORB</td><td>MORB</td><td>Peridotite</td><td>Rhyolite</td><td>Rhyolite</td></tr><tr><td>Pressure (GPa)</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td></tr><tr><td>length before (μm)</td><td>1400</td><td>1771.2</td><td>1061.9</td><td>1098.8</td><td>1124.6</td><td>1541.9</td><td>1141.3</td><td>1048.2</td><td>1652.9</td><td>1045.1</td><td>1336.4</td></tr><tr><td>length after (μm)</td><td>1334.2</td><td>1780.0</td><td>-</td><td>1462.7</td><td>1213.2</td><td>1666.1</td><td>1287.6</td><td>1324.9</td><td>1643.6</td><td>1309.1</td><td>1287.4</td></tr><tr><td>expected Melting T° (K)</td><td>>1600</td><td>>1523</td><td>>1600</td><td>>1173</td><td>>1723</td><td>>1173</td><td>>1600</td><td>>1600</td><td>>1523</td><td>>1173</td><td>>1173</td></tr><tr><td>Max T (K)</td><td>1614</td><td>1681</td><td>1612</td><td>1146</td><td>1769</td><td>1219</td><td>1609</td><td>1284</td><td>1573</td><td>1325</td><td>1073</td></tr><tr><td>Run duration (min)</td><td>544</td><td>1089</td><td>953</td><td>480</td><td>832</td><td>800</td><td>834</td><td>439</td><td>767</td><td>761</td><td>610</td></tr></tbody></table>

realized with Eq. (13), errors of raw data are used as weights in the inversion. The standard deviations are usually between 1 and 10% of the real asymptotic diffusivity. If the averaging method is selected, the standard deviation, estimated via Eq. (15), is usually about 1% of the final value.

A final source of uncertainties come from other technical issues and apparatus reproducibility, which are inherent to such challenging experiments. We considered that final error must be majored by 5% of the value. The relative uncertainties become even higher after conversion into thermal conductivities due to uncertainties and simplifications on the sample density and heat capacity at high pressure and temperature (see Supplementary Text S4).

## 2.4. Experimental procedure

We performed a suite of three heating and cooling cycles to provide an important number of data and maximize the quality of recovered data. A first cycle was run up to 500 °C for moisture removal by 50 °C steps (Fig. 2). The second cycle was run up to glass transition temperature (Tg) by 100 °C steps. Tg varies according to rock composition and its water content (Giordano et al., 2005). For our dry samples, the temperature of 800 °C happened to be above Tg for all our compositions. The temperature was then decreased down to 150 °C by 100 °C steps. The quality of the measurement is usually better during the cooling cycle once (1) the sample has thermally equilibrated with the assembly, (2) moisture has been removed, and (3) a better contact was achieved between the sample and the thermocouples due to local flow in the solid-state. The redistribution of matter under high temperature cancels most of the potential artifacts associated with the presence of pores/voids between the sample and the thermocouples (Hofmeister et al., 2009). Hence, data recorded during the last cycles were considered for final values.

In the last cycle, we heated the sample up to its melting temperature by 100 °C steps or 50 °C near the melting point. The sample was maintained above its melting point for less than 1 h, to avoid sample leakage and chemical reaction with the surrounding parts. The sample was then quenched to 600–550 °C. This step produces a glassy sample from which measurements at low temperature are performed.

Temperature oscillations of the two thermocouples were measured at each temperature step for both heating and cooling cycles, at 12 different oscillating frequencies between 0.1 and 1.5 Hz. The measurements were performed after at least 2 min of thermal equilibration to reach a stable regime (the smaller is the frequency, the longer is the required equilibration time). Then, the recording duration was 1 to 2 min or at least 10 oscillating periods. The measurement at the 12 frequencies took between 20 and 30 min. Thus, the duration of each experiment was more than 12 h.

## 2.5. Chemical and textural analyses

Recovered samples were cut parallel and perpendicular to the cylindrical furnace and each section was polished with great care. We could observe the position of the two thermocouples junctions, measure the sample radius and perform the textural and chemical analyses (see

Supplementary Figs. S6, S7, and S8).

Micro-textures were observed with a Scanning Electron Microscope (SEM) JEOL Jeol JSM-5910 LV using an accelerating voltage of 15 kV and a working distance of 11.4 mm. The 2D phase proportions of our partially molten peridotite samples were obtained from the analyses of qualitative chemical maps obtained by energy-dispersive X-ray spectroscopy (EDX) in the SEM. The images were binarized and phases were individually separated allowing textural analyses with the FOAMS software (Shea et al., 2010). A more detailed description is given in Freitas et al. (2019). On the other hand, quantitative chemical analyses were performed on both our starting materials and experiment products using the electron probe micro analyzer (EPMA). Chemical and textural analyses of starting materials are reported in the Supplementary Text S2, Figs. S6 to S13 and Tables S1, S2, S3 and S4. Analyses of recovered runs are detailed in Supplementary Text S3, supplementary Figs. S6 to S7 and Tables S1, S2, S3 and S4.

Water contents were estimated using the ICP-AES for the peridotite starting materials and via Raman spectroscopy for the recovered samples. Raman spectra were collected with a Renishaw InVia confocal Raman micro spectrometer, equipped with a 532 nm diode laser and a Leica DM 2500 M optical microscope. Measurements were carried out using a 2400 grooves/mm grating, a 100× microscope objective, a slit aperture set to either 20 µm or 65 µm and a laser power of 8 mW for glasses and 16 or 75 mW for olivine. The resulting lateral and axial resolutions were of ~1 and 3 µm, respectively, and the spectral resolution was better than 1 cm⁻¹. Daily calibration of the spectrometer was performed based on the 520.5 cm⁻¹ peak of Si. Spectra were recorded from ~100 to 1300 cm⁻¹ (alumino-silicate network domain) and from ~3000 to 3800 cm⁻¹ (water domain), with variable acquisition times ranging between 5 and 120 s for silicate bands and 120 and 240 s for water domain depending on the water content (Figs. S11 to S17). For water quantification in olivine and glass, we followed the procedures reported by Bolfan-Casanova et al. (2014) and Schiavi et al. (2018). We used both (1) the external calibration procedure, which is based on a set of hydrous olivine standards from (Bolfan-Casanova et al., 2014) and different types of silicate glasses ranging from basaltic to rhyolitic compositions (Schiavi et al., 2018; Médard et al., 2008), and (2) the internal calibration procedure, based on the correlation between the water concentration in olivine or glass and the relative areas of the water and silicate Raman bands (OH/Si integrated intensity ratio). The discrepancy between the two methods is small. Water contents in the standard materials were previously determined using the FTIR technique.

## 3. Results

We performed a total of 11 thermal diffusivity experiments on various chemical compositions, using the Angström method (details in Table 1). In this section, we first describe phase shifts and amplitude ratios between the two thermocouples and their conversion into thermal diffusivities and conductivities. We detail the post-mortem chemical and textural analyses in the Supplementary Text S3, Tables S1 to S4 and Figs. S6 to S13.
### 3.1. Phase shift, amplitude ratio and the refined raw-thermal diffusivities

Signals recorded during the experiments are sinusoidal oscillations with varying frequencies. The thermocouple located at the sample center (labeled TC1 in Fig. 3) presents a phase delay and a smaller amplitude compared to the thermocouple located at the sample edge (labeled TC2 in Fig. 3). Typical raw signals and their fits are presented in Fig. 3, while the refined (absolute) phase shifts and amplitude ratio are

reported in Fig. 4. Both the magnitude of phase shift and amplitude ratio change significantly with the type of sample and the experimental conditions, including the source frequency. At a constant frequency, the phase shift increases with increasing temperature. At a constant temperature, the phase shift increases with increasing the excitation frequency. The amplitude ratio is decreasing with increasing frequency and temperature.

Globally, the refined thermal diffusivities present a comparable

<figure><img src="image_261.png" /><figcaption>**Fig. 3.** Examples of temperature waves recorded at high temperature during M804 experiments. Such data set is acquired after the sample equilibration at target temperature conditions for a few minutes (here T = 250 °C). Acquisitions are performed at an interval ranging from 40 ms to 200 ms. At least 10 periods of the temperature oscillation are recorded to ensure a good extraction of phase shifts and amplitude ratio based on the fitting procedure described in Methods. Upper and lower frames correspond to temperature measurements performed at frequencies of 0.3 and 1 Hz, respectively. Red symbols are measured data and the dash thick lines represent the fitted sine waves. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_262.png" /><figcaption>**Fig. 4.** Example of refined phase shifts (top panels) and amplitude ratio (bottom panels) as a function of frequency (from 0.2 to 1.5 Hz, left panels) and temperature (from 428 K to 1191 K, right panels) for sample M844.</figcaption></figure>

evolution of temperature at all signal frequencies (see an example in Fig. 5). When the frequency dependence is larger than the experimental uncertainty (Fig. 6), we use Eq. (13) to refine the true asymptotic value of the thermal conductivity. Alternatively, when the temperature dependence is below the experimental uncertainty or when no clear frequency trends is visible, we average the different raw-diffusivity values (see Methods). In the wide majority of the cases, diffusivity values inferred from phase shifts appear to be more robust and with a lesser degree of uncertainty, compared to values inferred from difference of amplitude, in agreement with previous studies (Kanamori et al., 1969; Xu et al., 2004). Hence, despite similar values obtained with the two methods, values refined from phase shift were preferred.

### 3.1.1. Results for peridotite

For peridotite, our two successful experiments present a smooth evolution with temperature, yielding diffusivities values decreasing from $1.7(\pm0.1)e^{-6}$ to $7.5(\pm0.4)e^{-7} \text{ m}^2/\text{s}$ for sample M804 and $1.5(\pm0.8)e^{-6}$ to $5.5e^{-7} \text{ m}^2/\text{s}$ for sample M847 (Fig. 6a), with an uncertainty of about 5 to 10% for each sample. We attribute the relative discrepancy to a difference of sample mineralogy between the two samples, due to the coarse grain size of the core drilled peridotite (see Supplementary Texts S2 and S3 and Figs. S8 and S9). Olivine has a higher intrinsic thermal diffusivity than the other phases present in peridotite such as pyroxenes and spinel (Hofmeister and Branlund, 2015). M804 may contain more olivine, inducing higher diffusivities than M847. In this
<figure><img src="image_263.png" /><figcaption>Fig. 5. Typical example of a set of raw values of the thermal diffusivity extracted from our measurements for a basaltic glass, presented as a function of the frequency of the heat source. Experimental temperatures range between ~430 and ~1200 K at a pressure of 2GPa (sample M844). Left and right panels correspond to thermal diffusivities calculated from phase shifts and amplitude ratios, respectively. The colour code is kept similar in both panels. The fit through the data is made following the equation D = D₀ + A exp (−f/f₀) (Xu et al., 2004), where D₀ is the high-frequency asymptote of the diffusivity (see text). Error bars are generally within the size of the symbol.</figcaption></figure>

respect, the results for two samples are thus compatible with each other and the differences are representative of the variability of thermal diffusivity that can be expected among the compositional variability of peridotites (from lherzolites to dunites).

Thermal diffusivities of our peridotite samples scale as a function of ~1/T, as expected from the standard equation provided by Hofmeister and Branlund, 2015:

$$D_{\text{lat}} = a^*T^{-b} + c^*T \quad (16)$$

where T is the temperature (K), and a, b and c are adjustable parameters (Table 2). The quality of the fit is excellent up to a temperature of ~1300 K corresponding to the onset of peridotite melting. We therefore exclude the data points above 1300 K to model the thermal properties of solid peridotite.

### 3.1.2. Results for basalts

For our basaltic samples, the measured thermal diffusivities plot over a broad range of values, from 1.0(±0.1)e⁻⁶ to 3.0(±0.2)e⁻⁷ m²/s (Fig. 6b), with an uncertainty of ~5% for each experiment. The discrepancy is particularly important at low temperatures. Also, a same sample yields thermal diffusivity values significantly different along the different cycles of the experimental procedure (between C2, H3 and the final quench, see Fig. S5). Samples with the highest diffusivities at low temperatures present a rapid decrease of diffusivity with increasing temperature. On the other hand, samples with the lowest diffusivities show very small temperature dependence. It yields to a convergence of all diffusivity measurements at ~1000 K. Based on a study of rhyolitic glasses, Romine et al. (2012) reported a moderate temperature dependence, similar to our samples presenting a low diffusivity, and an increase of 0.0192 mm²/s of the glass diffusivity per percent of crystallinity. For our starting materials with less than 5 vol% of crystals (see Supplementary Text S2 and S3, Fig. S7), the effect of microlites could account for ~0.1 mm²/s of variation in our diffusivity values, which correspond to less than ~10% of observed differences. However, the presence of 35 vol% of crystals in the recovered sample M807 would explain not only its high thermal diffusivity at low temperature, but also its strong temperature dependence that is typical of crystals (see Beck

et al., 1978; Romine et al., 2012; Hofmeister and Branlund, 2015 and our peridotite trends in Fig. 6a). The recrystallization of M807 is not surprising, since the second cycle of annealing was performed above its Tg (Figs. 2 and S5). The crystallinity of other samples depends on the applied cycles of annealing at a temperature eventually close to their Tg. Nonetheless, diffusivities of all samples converge at increasing temperatures, because the conductivity of crystals is not much greater than that of the glass at high temperature, especially if microlites are low thermal diffusivity silicates such as pyroxenes, plagioclase or spinels (see Hofmeister and Branlund, 2015 for mineral thermal diffusivity compilations). Over the 5 basaltic samples investigated in this study, the diffusivity trends indicate either significant recrystallization of M807 (DR07-MORB) and M662 (EPR-MORB), or negligible crystallization of M844 and M846 (DR11-MORB) and M836 (synthetic haplobasalt). The sample crystallization is likely to evolve during the thermal diffusivity measurements in step H3 (up to melting point) of the experiments performed at a temperature significantly above the Tg. For this reason, the crystallinity determined on the recovered samples is only a qualitative measurement of the sample properties at high temperatures.

### 3.1.3. Results for rhyolites

Measured thermal diffusivities of rhyolite samples also plot over a broad range of values from 1.4(±0.8)e⁻⁶ to 4.0(±0.2)e⁻⁷ m²/s, with an uncertainty of about 5% for each experiment (Fig. 6c). The discrepancy appears similar than for the basalt samples, as the sample presenting higher diffusivities also show a major temperature dependence at low temperature. The presence of less than 2 vol% of crystals in the starting material could account for a diffusivity increase of 0.038 mm²/s at maximum (Romine et al., 2012). Still, the diffusivity trends suggest a major recrystallization at high-temperature for M808 (Güney Dag) and M843 (ATHO), some crystallization for M843 (ATHO), and negligible crystallization for M850 (Güney Dag).

### 3.1.4. Properties of melts and partially molten samples

In addition to the evolutions described above, a strong decrease of the thermal diffusivity is observed for most of our samples at the highest temperatures (Fig. 6). The decrease occurs at temperatures around 1300 K for peridotites, 1200 K for basalts and > 1050 K for rhyolites. Such
<figure><img src="image_264.png" /><figcaption>**Fig. 6.** Thermal diffusivities data refined for our peridotites, basalts and rhyolites samples (from top to bottom) as a function of temperature. Associated fits are represented with a solid line of the same colour than the symbols. Error bars are shown when larger than the symbol size. Several basaltic (M807, M662) and rhyolite (M848, M808) samples were partially recrystallized during the annealing and dehydration procedures prior to the conductivity measurement, as indicated by their high conductivity values at low temperatures. Thermal diffusivities values converge at high temperatures and melting.</figcaption></figure>

temperatures are in agreement with the melting or glass transition temperatures, depending if the sample is a peridotite or a glass, recrystallized or not. Similar changes were already reported at temperatures above the glass transition (Hofmeister et al., 2009; Hofmeister et al., 2014; Romine et al., 2012).

The amplitude of the decrease is 45–50% in peridotites, which recovered samples present a degree of partial melting (F) up to 23%, 35–70% for basalts and < 30% for rhyolites samples. The more pronounced decrease in molten peridotites is probably due to a more contrasted change of the local structure at the melting point, between the
<table><caption>**Table 2**
Fitting parameters obtained for the last cycle (C2-H3) of each experiment. Associated errors (1 STD) are given in parenthesis.</caption><thead><tr><th>Run#</th><th>a</th><th>b</th><th>c</th><th>R<sup>2</sup></th><th>R<sup>2</sup> adj</th></tr></thead><tbody><tr><td colspan="6"><strong>Peridotites</strong></td></tr><tr><td>M804</td><td>7.0e-4 (2.2e-4)</td><td>1.06 (0.06)</td><td>2.97e-10 (2.8e-11)</td><td>0.9659</td><td>0.9647</td></tr><tr><td>M847</td><td>8.8e-4 (7e-4)</td><td>1.15 (0.05)</td><td>2.37e-10 (1.8e-11)</td><td>0.9763</td><td>0.9750</td></tr><tr><td colspan="6"><strong>Basalts</strong></td></tr><tr><td>M846</td><td>8.6e-6 (1.5e-5)</td><td>0.5 (0.3)</td><td>9.7e-11 (8.4e-11)</td><td>0.3933</td><td>0.3294</td></tr><tr><td>M662</td><td>7.2e-4 (1.8e-4)</td><td>1.20 (0.04)</td><td>1.01e-10 (4.1e-12)</td><td>0.9964</td><td>0.9962</td></tr><tr><td>M836</td><td>5.6e-5 (7.6e-5)</td><td>0.84 (0.23)</td><td>1.03e-10 (3.5e-11)</td><td>0.5748</td><td>0.5421</td></tr><tr><td>M844</td><td>1.9e-5 (1.1e-5)</td><td>0.63 (0.10)</td><td>8.0e-11 (2.3e-11)</td><td>0.9608</td><td>0.9556</td></tr><tr><td>M807</td><td>8.1e-4 (1.9e-4)</td><td>2.15 (0.04)</td><td>4.1e-11 (9e-12)</td><td>0.9933</td><td>0.9929</td></tr><tr><td colspan="6"><strong>Rhyolites</strong></td></tr><tr><td>M850</td><td>1.9e-5 (3.5e-5)</td><td>0.70 (0.32)</td><td>1.7e-10 (7e-11)</td><td>0.2788</td><td>0.2101</td></tr><tr><td>M808</td><td>1.8e-2 (7e-2)</td><td>1.7 (0.7)</td><td>2.3e-10 (1.4e-10)</td><td>0.8003</td><td>0.7696</td></tr><tr><td>M848</td><td>1.4e-4 (7e-5)</td><td>0.91 (0.09)</td><td>1.3e-10 (2.0e-11)</td><td>0.9851</td><td>0.9826</td></tr><tr><td>M843</td><td>5.1e-3 (3.3e-3)</td><td>1.42 (0.11)</td><td>3.18e-10 (3.6e-11)</td><td>0.9715</td><td>0.9684</td></tr></tbody></table>

minerals and the melt (see discussion). On the other hand, the minor change of diffusivity for rhyolites at Tg could be related to their high SiO₂-content, which preserves a polymerized structure in the melt above the glass transition.

## 3.2. Thermal conductivities

Thermal conductivities ($\kappa$) can then be computed from thermal diffusivities following:

$$ \kappa(P, T) = D(P, T)^* \rho(P, T)^* C_p(P, T) \quad (17) $$

where $\rho(P, T)$ and $C_p(P, T)$ are the sample density and heat capacity, respectively. For our calculations, we considered $\rho$ and $C_p$ values in standard conditions when the P and/or T dependences were not available in the literature (see Supplementary Text S4 and Fig. S4).

Conductivities calculated for peridotites evolve from 4.32 (±0.48) to 1.70 (±0.19) W/m/K with increasing the temperature up to the melting point. Values for basalts range from 1.7 (±8.0e⁻²) to 0.5 (±4.0e⁻²) W/m/K, whereas those for rhyolite lie between 1.1 (±9.0e⁻²) and 0.3 (±3.0e⁻²) W/m/K for low to high temperatures, respectively. At the melting temperature, our partially molten peridotites display values of 1.19 (±0.16) to 0.93 (±0.10) W/m/K, whereas melts present relatively lower values of 0.34 (±0.2) to 0.18 (±0.1) W/m/K for basaltic and 0.41 (±0.3) to 0.31 (±0.3) for rhyolitic compositions (See Fig. 7 and Table 3). We note that the important difference in conductivity between basalts and rhyolites (Fig. 7) is predominantly due to differences between their heat capacities and densities, while their thermal diffusivities are found similar (Fig. 6).

# 4. Interpretation of results

## 4.1. General temperature dependence

For all compositions investigated in this study, thermal diffusivities decrease with increasing temperature until reaching a plateau at temperatures between 700 and 1000 K. Based on experiments performed at room pressure, it was observed that the plateau occurs at about the Debye temperature of the mantle minerals (Hofmeister et al., 2009, 2014). For this reason, it was proposed that thermal diffusivities vary

largely with temperature until the complete activation of the vibration modes (phonons in minerals). The temperature range observed in our study for the occurrence of a plateau is fully compatible with this interpretation.

For basaltic glasses, a comparable but more moderate decrease of thermal diffusivity was reported up to a saturation temperature corresponding well to the glass Tg (Hofmeister et al., 2009, 2014; Romine et al., 2012). In our experiments, the decrease is of ~30% to 60% over the investigated temperature range, depending on the experiment. Such amplitude is compatible with the ~40% decrease observed during the heating of pyroxene glasses (Hofmeister et al., 2009).

For rhyolite samples, the thermal conductivity increases slightly with increasing the temperature (Fig. 7). This is due to the heat capacity that increases more with temperature than the density increases and diffusivity decreases (Eq. (17)). The increase is, however, smaller than reported in Romine et al. (2012), due to the use of a different Cp (Neuville et al., 1993) (see Supplementary Text S4) and a stronger temperature dependence of thermal diffusivities observed in our experiments because of different crystallizations states.

## 4.2. Effect of radiative conduction

Romine et al. (2012) reported an increase in thermal diffusivities of the melts at very high temperatures at ambient pressure. They attributed this feature to an increased role of the radiative component. Such a component can dominate the thermal diffusivity for a sample transparent to the infrared and visible photons at high temperatures. For thin samples, this effect can become problematic if the mean free path of photons is longer than the sample length (ballistic photons, see Hofmeister and Branlund, 2015). No significant increase in thermal diffusivity and conductivity is observed in our high-pressure experiments, except maybe for the rhyolitic samples (Figs. 6 and 7). The difference with the previous work is most probably related to the opacity of our basalts and peridotites samples, hence limiting the radiative transfers.

## 4.3. Effect of glass/melt composition

Overall, our conductivity values are compatible with the values available in the literature (see Fig. 8). Differences in absolute values are nonetheless present. Rhyolite melts (0.31–0.41 W/m/K) are found slightly more conductive than basaltic ones (0.18–0.34 W/m/K). Within the same family of glass, thermal conductivity varies by 0.10 to 0.15 W/m/K for rhyolitic and basaltic melts, respectively. This is slightly larger (of at least 10%) than the experimental uncertainty, an effect possibly due to larger uncertainties on the dimensions of the molten sample. No clear trend can be retrieved from the comparison between our different basaltic or rhyolitic compositions. Among the major elements, iron could be of major importance, due to its critical impact on glass/melt density. Indeed, the thermal diffusivity of glasses was reported to decrease with increasing density (Hofmeister et al., 2014). The comparison between our rhyolites and basalts is coherent with such a trend. However, our haplobasalt (M836) presents diffusivity values comparable with natural basalts (M662, M807, M844), as well as Fe-bearing (M843 and M848) and Fe-free (M808 and M850) rhyolites (Fig. 6), despite a variation of the Fe-content from 0 to 10 wt% in these different samples. Other elements could also impact the melt thermal diffusivity, in particular Si and Al, which favor polymerization of the liquid (and alkali elements for the opposite effect) (Ni et al., 2015). Still, within the experimental uncertainty, we observe no clear trend related to these elements, despite a variation of the SiO₂ content by more than 20%. Additionally, water, with a total content smaller than 1.10 wt%, should have a negligible effect on thermal conductivity (Romine et al., 2012; Ni et al., 2015). It could, however, impact the melt density at a low degree of partial melting (Hofmeister et al., 2014). As the water contents estimated in the recovered samples are similar to the ones obtained in starting materials (see Supplementary Text S2 and S3), water should not
<figure><img src="image_265.png" /><figcaption>**Fig. 7.** Thermal conductivities inferred from diffusivities (Fig. 6), heat capacities and densities (see Fig. S4), based on Eq. (17) of the main text.</figcaption></figure>
<caption>**Table 3**</caption>

Thermal diffusivities and conductivities retrieved for our samples at low and high temperatures for solids, above Tg and/or melting temperature (indicated for each group) for melts and partially molten samples.

<table><thead><tr><th>Sample name</th><th>Temperature range (K)</th><th>Melt fraction (%)</th><th>STD (%)</th><th>K (W/m/K)</th><th>STD</th><th>D (m²/s)</th><th>STD</th></tr></thead><tbody><tr><td colspan="8"><strong>SOLIDS</strong></td></tr><tr><td>M804</td><td>LT &lt; 400 K</td><td>0</td><td>0</td><td>4.32</td><td>0.48</td><td>1.66E-06</td><td>2.05E-09</td></tr><tr><td>M847</td><td>LT &lt; 400 K</td><td>0</td><td>0</td><td>3.28</td><td>0.37</td><td>1.26E-06</td><td>3.34E-09</td></tr><tr><td>M804</td><td>HT &gt;1200 K</td><td>0</td><td>0</td><td>2.29</td><td>0.26</td><td>7.33E-07</td><td>1.80E-09</td></tr><tr><td>M847</td><td>HT &gt;1200 K</td><td>0</td><td>0</td><td>1.70</td><td>0.19</td><td>5.39E-07</td><td>6.13E-10</td></tr><tr><td colspan="8"><strong>Partially molten samples</strong></td></tr><tr><td>M804</td><td>&gt;1300 K</td><td>6.44</td><td>2.98</td><td>1.19</td><td>0.16</td><td>3.74E-07</td><td>2.55E-08</td></tr><tr><td>M847</td><td>&gt;1300 K</td><td>23.31</td><td>2</td><td>0.93</td><td>0.10</td><td>2.92E-07</td><td>1.22E-10</td></tr><tr><td colspan="8"><strong>MELTS</strong></td></tr><tr><td>BASALTS</td><td>MELT &gt;1200 K</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>M807</td><td>&gt;1224 K</td><td>100</td><td>0</td><td>0.34</td><td>0.02</td><td>1.84E-07</td><td>1.51E-09</td></tr><tr><td>M844</td><td>&gt;1326 K</td><td>100</td><td>0</td><td>0.18</td><td>0.01</td><td>9.52E-08</td><td>2.96E-11</td></tr><tr><td>M662</td><td>&gt;1450 K</td><td>100</td><td>0</td><td>0.22</td><td>0.16</td><td>1.20E-07</td><td>8.73E-08</td></tr><tr><td>M846</td><td>&gt;1200 K</td><td>100</td><td>0</td><td>0.22</td><td>0.12</td><td>1.16E-07</td><td>6.47E-08</td></tr><tr><td>M836</td><td>1180 K</td><td>100</td><td>0</td><td>0.28</td><td>0.02</td><td>1.51E-07</td><td>6.94E-10</td></tr><tr><td>rhyolite</td><td>Melt &gt;1000 K</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>M850</td><td>&gt; 1080 K</td><td>100</td><td>0</td><td>0.31</td><td>0.03</td><td>3.05E-07</td><td>2.42E-08</td></tr><tr><td>M808</td><td>&gt;1060 K</td><td>100</td><td>0</td><td>0.35</td><td>0.03</td><td>3.47E-07</td><td>5.95E-10</td></tr><tr><td>M848</td><td>&gt;1260 K</td><td>100</td><td>0</td><td>0.39</td><td>0.03</td><td>3.91E-07</td><td>6.21E-10</td></tr><tr><td>M843</td><td>&gt; 1280 K</td><td>100</td><td>0</td><td>0.41</td><td>0.03</td><td>4.16E-07</td><td>2.36E-10</td></tr><tr><td colspan="8"><strong>Literature</strong></td></tr><tr><td colspan="8"><strong>SOLIDS</strong></td></tr><tr><td>Olivine Xu 0 GPa LT</td><td>294 K</td><td>0</td><td>0</td><td>4.15</td><td>0.42</td><td>4.10E-06</td><td>2.05E-07</td></tr><tr><td>Olivine Xu 0 GPa HT</td><td>1377 K</td><td>0</td><td>0</td><td>1.95</td><td>0.19</td><td>5.00E-07</td><td>2.50E-08</td></tr><tr><td>Olivine Hof 16 AVG</td><td>300 K</td><td>0</td><td>0</td><td>6.37</td><td>0.64</td><td>2.37E-06</td><td>4.75E-08</td></tr><tr><td>Olivine Hof 16 AVG</td><td>&gt;1500 K</td><td>0</td><td>0</td><td>2.76</td><td>0.28</td><td>7.01E-07</td><td>1.40E-08</td></tr><tr><td>Dunite Per 06</td><td>300 K</td><td>0</td><td>0</td><td>4.22</td><td>0.42</td><td>1.57E-06</td><td>3.15E-08</td></tr><tr><td>Dunite 2 Per 06</td><td>300 K</td><td>0</td><td>0</td><td>5.35</td><td>0.53</td><td>1.99E-06</td><td>3.99E-08</td></tr><tr><td>Dunite Per 06</td><td>&gt;1500 K</td><td>0</td><td>0</td><td>1.65</td><td>0.17</td><td>4.20E-07</td><td>8.40E-09</td></tr><tr><td>Dunite 2 Per 06</td><td>&gt;1500 K</td><td>0</td><td>0</td><td>2.23</td><td>0.22</td><td>5.65E-07</td><td>1.13E-08</td></tr><tr><td>Dunite Beck 78 2Gpa</td><td>300 K</td><td>0</td><td>0</td><td>4.15</td><td>0.42</td><td></td><td></td></tr><tr><td>Dunite Beck 78 2Gpa</td><td>300 K</td><td>0</td><td>0</td><td>5.49</td><td>0.55</td><td></td><td></td></tr><tr><td>Dunite 2 GPa Dry Z19</td><td>300 K</td><td>0</td><td>0</td><td>5.21</td><td>0.26</td><td>1.71E-06</td><td>8.55E-08</td></tr><tr><td>Dunite 2 GPa 0.08% wt H₂O Z19</td><td>300 K</td><td>0</td><td>0</td><td>4.86</td><td>0.24</td><td>1.59E-06</td><td>7.95E-08</td></tr><tr><td colspan="8"><strong>Melts</strong></td></tr><tr><td>Di Hof 09</td><td>&gt;1140 K</td><td>100</td><td>0</td><td>1.25</td><td>0.13</td><td>2.89E-07</td><td>5.79E-09</td></tr><tr><td>An Hof 09</td><td>&gt;1290 K</td><td>100</td><td>0</td><td>1.45</td><td>0.15</td><td>3.58E-07</td><td>7.16E-09</td></tr><tr><td>Fk Hof 09</td><td>&gt;1245 K</td><td>100</td><td>0</td><td>1.45</td><td>0.15</td><td>4.75E-07</td><td>9.50E-09</td></tr><tr><td>Ab Hof 09</td><td>&gt;1300 K</td><td>100</td><td>0</td><td>1.60</td><td>0.16</td><td>4.99E-07</td><td>9.99E-09</td></tr><tr><td colspan="8"><strong>Rhyolite/obsidian</strong></td></tr><tr><td>Ncar Ro 12</td><td>1384 K</td><td>100</td><td>0</td><td>1.53</td><td>0.15</td><td>5.17E-07</td><td>1.03E-08</td></tr><tr><td>NCr Ro 12</td><td>1365 K</td><td>100</td><td>0</td><td>1.50</td><td>0.15</td><td>5.08E-07</td><td>1.02E-08</td></tr><tr><td>SCE Ro 12</td><td>1464 K</td><td>100</td><td>0</td><td>1.59</td><td>0.16</td><td>5.31E-07</td><td>1.06E-08</td></tr></tbody></table>

Literature estimates are also provided on the bottom part: “Xu” refers to polycrystalline olivine (dunite) measurements performed at ambient pressure (Xu et al., 2004). “Per06” refers to average values from 3 oriented (100,010,001) single crystal measurements at room pressure (Pertermann and Hofmeister, 2006) and some dunite values are also available in the same study. “Beck78” refers to Dunite measurements performed at 2 GPa (Beck et al., 1978). Regarding melts, literature data are taken from Hofmeister et al. (2009), Di: diopside, An: anorthite, Fk: orthose, Ab: albite. Finally, data from a few rhyolites Ncar, NCr and SCE are selected from Romine et al. (2012). When literature data are available in the form of thermal diffusivities, their conductivities are obtained using either the Cp and ρ provided in the source study or from our own parameters if not available (see Supplementary Materials).

induce any strong diffusivity variation in our data-set. We, therefore, conclude that chemical effects are secondary compared to the structural ones.

### 4.4. Mixing models for thermal conductivities of partially molten peridotite

Our results show that peridotites, and to a lower extent the glasses, present higher thermal diffusivities compared to the melts. This result is not surprising: thermal conductivity is strongly dependent on the local structure and its vibrational properties. The disorder in the atomic structure increases from solid, glasses to molten states (Hofmeister et al., 2014). The composition of the melt appears to be of secondary importance. The final set of thermal diffusivities and conductivities values selected for the applications are given in Tables 3, S5 and Figs. 9 and S14.

To propose a predictive law for thermal conductivity of partially molten rocks several mixing equations are now tested. Such equation generally describes the effect of a small amount of conductive phase into

an insulating matrix. For thermal conduction, the problem is reversed because the melt is less conductive than the solid. In this section, we explore the different predictive models of thermal conductivity of binary mixtures:

1) Linear mixing model consider parallel thermal resistor:

$$ \kappa_{\text{bulk}} = \kappa_s^*(1-F) + F\kappa_m \quad (18) $$

where s corresponds to the solid, m the melt and F the volume fraction of the melt.

2) Tube / Ashbie model consider 1/3 of tubes of melt aligned in the heat flow direction (Grant and West, 1965; Schmeling, 1986):

$$ \kappa_{\text{bulk}} = \frac{1}{3} * \kappa_s^*(1-F) + F\kappa_m \quad (19) $$
<figure><img src="image_266.png" /><figcaption>**Fig. 8.** Comparison with the literature of our thermal diffusivities obtained our peridotite (upper panel), basaltic (center panel) and rhyolitic (lower panel) starting materials. Solids (peridotite and glasses) and melts are represented in thick lines and thick dashed line, respectively. Literature data are reported with thinner dashed lines (grey to black). Diopside (Di) glass and melt, XAlSi₃O₈ glass, Anorthite melt (An), Albite melt (Ab) and orthose melt (FK) from Hofmeister et al., 2009, NCr Rhyolite from Romine et al., 2012 and Enstatite glass from Hofmeister et al., 2014.</figcaption></figure>
<figure><img src="image_267.png" /><figcaption>Fig. 9. Our best model of the evolution of thermal diffusivity as a function of the melt fraction at high temperature (>1200 K), based on all available data sets. Solid peridotites at low and high temperature are represented with green and blue diamonds, respectively. Partially molten peridotites with 6.4% and 23.3% of melt (see Supplementary Text S3) are represented with red circles and melts (basalts and rhyolites) are represented with blue and pink circles. All the errors represented are 1 standard deviation. The fit represented is using a model of thermal resis- tors in series (Eq. (22)), which fits our data set very well (see Fig. S14). (For interpretation of the refer- ences to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

3) Cube model (Waff, 1974) representing cubes of solids into a melt matrix

$$ \kappa_{\text{bulk}} = [1 - F^{2/3}]^{*} \kappa_s \qquad (20) $$

4) Archie's law, an empirical relation developed for electrical conduc- tivity (Watanabe and Kurita, 1993)

$$ \kappa_{\text{bulk}} = C^* (1 - F)^{n^*} \kappa_s \qquad (21) $$

where C and n are constants.

5) Thermal resistors in series:

$$ \kappa_{\text{bulk}} = \kappa_s^* \frac{\kappa_m / \kappa_s}{\kappa_m / \kappa_s + F(1 - \kappa_m / \kappa_s)} \qquad (22) $$

6) Hashin Shtrikman lower bound (HS⁻), representing insulating melt spheres into a conductive solid matrix:

$$ \kappa_{\text{bulk}} = \kappa_s^* \frac{F}{1 / (\kappa_m / \kappa_s) + \frac{1-F}{3\kappa_s}} \qquad (23) $$

7) Maxwell-Eucken relation:

$$ \kappa_{\text{bulk}} = \kappa_s^* \frac{\kappa_m + 2\kappa_s + 2F(\kappa_m - \kappa_s)}{\kappa_m + 2\kappa_s - F(\kappa_m - \kappa_s)} \qquad (24) $$

8) Landauer relation, based on resistors in series:

$$ \kappa_{\text{bulk}} = \frac{1}{4} \left[ \kappa_m (3F - 1) + \kappa_s (2 - 3F) + \left\{ (\kappa_m (3F - 1) + \kappa_s (2 - 3F)) 2 + 8\kappa_s \kappa_m \right\}^{\frac{1}{2}} \right] \qquad (25) $$

9) Russel-Rayleigh relation:

$$ \kappa_{\text{bulk}} = \frac{\kappa_s [\kappa_s + F^{2/3}(\kappa_m - \kappa_s)]}{\kappa_s + (\kappa_m - \kappa_s)(\kappa_m^{2/3} - F)} \qquad (26) $$

These equations provide a different evolution of the thermal con- ductivity with the fraction of melt (F). In Fig. S14, we present the results when either fixing κsolid and κmelt to the average value of our measure- ments, or adjusting their values to minimize the misfit between the mixing models and our results at varying F values (Table S5). Among the variety of fits obtained, it appears that the effect of partial melting is underestimated in most of the cases. Only the thermal resistors in series is capable to reproduce adequately the strong curvature observed experimentally at low F values, as well as the end-member values of solid and melt conductivities. For this reason, we use series model (Eq. (22)) for further discussions.

We note a lack of data points at high melt fractions, preventing to decipher more precisely the mixing trend. At very high temperatures, the experimental measurements become difficult on natural peridotite melting, in particular due to melt escape and chemical reactions with the experimental cell. Complementary data could be acquired working with analog system such as basalt/olivine mixture and could represent a further research direction.

4.5. Thermal conductivity of partially molten peridotite: influence of texture

The 3D solid/melt arrangement is known to influence significantly the geophysical properties of partially molten systems (Von Bargen and Waff, 1986; Laporte et al., 1997; Laporte and Provost, 2000; Minarik and Watson, 1995; Yoshino et al., 2005; Maumus et al., 2005; ten Grotenhuis et al., 2005; Freitas et al., 2019; Freitas and Manthilake, 2019). Their distribution is classically described using dihedral angle value, which translates the ability of a liquid to wet the grain boundaries as a consequence of interfacial energies. The dihedral angle decreases when increasing pressure, temperature, water content and decreasing silica/ alumina content of the melt (Yoshino et al., 2007; Mibe et al., 1998, 1999; Laporte et al., 1997; Watson et al., 1991). Several studies with similar basaltic or peridotite melts (dry) have shown that basalt-like melts at mantle conditions have dihedral angles significantly lower than the interconnection threshold of 60°, with values between 30 and
40° at 2 GPa (Laporte et al., 1997; Yoshino et al., 2005, 2007). Our partially molten samples display a coherent texture and dihedral angles with these observations. Dihedral angles of 23.3° and 19.0° were measured from our samples containing 6.4% and 23.3% of melt, respectively (Fig. S15), in good agreement with previous data given that these mafic melts are moderately hydrous (Table S3). For each melt fraction, a thin layer of melt surround most of the grains, in particular olivines, which in 3D will result into the insulation of the solid gains from their surroundings (Fig. S9). This is very well visible on our low melt fraction sample (F = 6.4%) where the layers of melt are few microns thick (M804). Even if melt is more abundant near clinopyroxenes and spinel sites in M847 (F = 23.3%), the melt pockets are interconnected with similar thin melt layers (<10 µm) (Fig. S9). As a result, thermal conductivity is expected to drop brutally in the first degrees of melting. Still, some grain boundaries should remain un-wetted until the melt fraction rise significantly. For this reason, thermal conductivity should only stabilize at melt fraction corresponding to solid grains completely isolated from each other. This trend is visible in our data and parallel model (Figs. 9 and S14) with a strong decrease in the first 10% of melting highlighted by M804, the change of slope seems to occur around 15% and values decrease more slowly in the 15–50% range as seen in M847, to stabilize and display near-melt values above 50%. The complete isolation of solids should occur at “packing” threshold, which is a function of the solid shapes and size distribution and is expected to occur between 40 and 60% of melting. Thus, the first degrees of melting are very crucial in the case of a wetting liquid and affecting importantly the thermal properties.

### 4.6. Thermal conductivity of peridotite: effect of the grain size

The modelling of thermal conductivity of peridotites and low F molten peridotites should also take in account the effect of grain boundaries thermal resistances as grain size may vary in the different geological contexts (from 100 µm to >1 cm), as seen in natural meteoric examples, (Barrat et al., 1999; Busek, 1997; Keil, 2010; Floran et al., 1978). Indeed, grain boundary scattering could be important when the mean free paths of phonons approach the grain size. This effect, which only concerns solids and low fractions of melt, can be quantified with the following equation (Smith et al., 2003; Smith et al., 2013):

$$ \frac{1}{\kappa_{\text{poly}}} = \frac{1}{\kappa_{\text{single}}} + nR_{\text{boundary}} \quad (27) $$

where $κ_{poly}$ is the thermal conductivity of the polycrystalline sample, $κ_{single}$ the thermal conductivity of a reference single crystal (average from olivine data of (Hofmeister et al., 2007, Table 3), n represents the surface of grain boundaries along the heat flow direction per unit length, and $R_{boundary}$ the thermal resistance of grain boundary plane. The n value should be almost constant with temperature (Smith et al., 2013) and is estimated between $4e^{-4}$ and $7e^{-3}$ m for our two peridotite samples via analyses of SEM images (grain size ranging from 25 to 140 µm, see Table S4 for textural parameters). We calculate $R_{boundary}$ values between $1.9e^{-6}$ and $5.9e^{-6}$ W/m²/K for M804 and $9.3e^{-6}$ and $3.2e^{-5}$ W/m²/K for M847. These values are compatible or slightly higher than hydrous polycrystalline olivine samples (Zhang et al., 2019).

As a result, the thermal conductivities quantified in our experiments are underestimations of natural ones as the grain size is <100 µm, in experiments compared to grain sizes of 100 µm up to >1 cm typical of mantle peridotites and reduced meteorites (which could be relic of bodies interiors, from cumulates (Floran et al., 1978), enstatite chondrite/achondrite (Keil, 2010), diogenite (Barrat et al., 1999) to pallasite (Busek, 1997)). Melts and high F partially molten systems are not affected by such effect, thus the observed decrease of thermal conductivity at the melting temperature is probably smaller in our experiments than in natural cases.

## 5. Implications for geodynamos on Mercury-like proto-planets

### 5.1. Suitable conditions for a dynamo

For a thermally driven dynamo to operate in a terrestrial planet, four conditions were found to be necessary (e.g., Monteux et al., 2011): the core heat flow must be at least adiabatic (1), the thermal convection within the core has to supply enough power to compensate the losses due to ohmic dissipation (2), the Reynolds magnetic number must be supercritical (complex turbulent convection) (3) and the mantle heat flow has to overcome the core heat flow needed to induce a dynamo (4). These conditions can be expressed in terms of heat flow balance (See Fig. 10 for a schematic representation) and are detailed here:

(1) The metallic core has to convect, meaning that the heat flow out of the core needs to overcome the adiabatic heat flow (Stevenson et al., 1983).

For this, the core thermal conductivity ($κ_{core}$) is a dominant parameter. A large $κ_{core}$ value increases the heat flux along the core adiabat and reduces the lifetime of a thermally driven dynamo (Breuer et al., 2015). Several laboratory measurements suggested that the thermal conductivity of polycrystalline iron at Mercury’s core conditions is 113–125 W/m/K (see Deng et al., 2013 and references therein). However, such values for Mercury are recently challenged with several recent studies proposing a much lower conductivity. In a first one, the conductivity of pure Fe and Fe-Si alloys is reported at 30–40 W/m/K and 35–40 W/m/K, respectively (Silber et al., 2019). Then, it is proposed that the thermal conductivity of Fe-S at the P-T conditions of Mercury’s core is as low as ~4 W/m/K, thus 1–2 orders of magnitude lower than that of pure iron (Pommier et al., 2019; Manthilake et al., 2019).

This first condition can be expressed as:

$$ Q_{\text{CMB}} > Q_{\text{Ad}} = \frac{\kappa_c \alpha_c g_c T_{\text{CMB}}}{C_{p,c}} 4\pi R_c^2 \quad (28) $$

To estimate this flux, we assume that $T_{CMB}$ is the melting temperature of pure iron at $P_{CMB}$. This assumption gives a conservative value of the core heat flow in comparison with considering $T_{ICB}$ since the core liquidus is steeper than the core adiabat. We estimate the relation between the melting temperature of pure iron and the pressure using the

<figure><img src="image_268.png" /><figcaption>Fig. 10. Schematic representation of the geometry and of the heat flows computed in our model.</figcaption></figure>
following expression obtained by fitting the experimental results from Anzelli et al., 2013 with a Simon and Glatzel equation:

$$T_{m,Fe} = 1800 \left( \frac{P_{CMB}}{27.9} + 1 \right)^{1/2.08} \quad (29)$$

Such a melting temperature typically lies between the solidus and liquidus of a chondritic mantle for the same pressure conditions (Monteux et al., 2020 and references therein). We also assume that $k_c$, $\alpha_c$ and $C_p$, $c$ are constant (see values in Table 4) and $P_{CMB}$ is calculated as follows (Monteux and Arkani-Hamed, 2014):

$$ \begin{aligned} P_{CMB} &= P(r = R_c) \\ &= \frac{2}{3}\pi G \rho_c^2 (R_c^2 - r^2) + \frac{2}{3}\pi G \rho_{Si}^2 (R^2 - R_c^2) + \frac{4}{3}\pi G \rho_{Si} R_c^3 (\rho_c - \rho_{Si}) \left(\frac{1}{R_c} - \frac{1}{R}\right) \end{aligned} \quad (30) $$

(2) The energy supplied by thermal convection to the geodynamo has to compensate for the loss due to ohmic decay (Buffett, 2002).

This imposes a condition on the core heat flow at the CMB. In fact, core heat flow will need to overcome a critical value. Assuming that dynamo is generated only by thermal convection in the core, we can write:

$$Q_{CMB} > Q_{Ad} + \frac{4\nu_c \bar{B}^2 C_{p,c}}{0.8\mu_c \alpha_c G \rho_c R_c} \quad (31)$$

This heat flux is estimated by considering that the characteristic magnetic length scale equals the radius of the core. The parameter $\bar{B}$ is the average strength of the magnetic field inside the core and was estimated using a scaling from Christensen and Aubert, 2006:

$$\bar{B} = 0.9 \mu_c^{1/2} \rho_c^{1/2} \left( \frac{g_c Q_B (R_c - r_i)}{4\pi R_c r_i} \right)^{1/2} \quad (32)$$

With $Q_B = \alpha_c Q_{CMB}/C_p$, $c$ the buoyancy flux and $r_i$ the radius of the inner core. In the scaling from Christensen and Aubert (2006), the inner core size cannot be set to 0. On Mercury the size of the inner core is currently not well constrained even if recent constraints via geodetic analysis (Genova et al., 2018) suggest its presence and a possible important size ($r_i/R_c$ between 0.3 and 0.7). As we focus here on the effect of thermal cooling on dynamo generation (i.e. we do not consider the effect of compositional convection related to inner core growth), we consider a small inner core with $r_i/R_c = 0.01$. The scaling law used to calculate the average strength of the magnetic field inside the core (Eq. (32)) is valid for the Earth but overestimates $\bar{B}$ in the case of thin shell dynamos such as the one operating within Mercury (Christensen and

<table><caption>**Table 4**

Typical parameter values for heat fluxes calculations displayed in Figs. 11 and 12.</caption><thead><tr><th>κ<sub>c</sub></th><th>Core thermal conductivity</th><th>4–40</th><th>W/m K<sup>-1</sup></th></tr></thead><tbody><tr><td>α<sub>c</sub></td><td>Core thermal expansion</td><td>7.7 × 10<sup>-5</sup></td><td>K<sup>-1</sup></td></tr><tr><td>C<sub>p, c</sub></td><td>Core heat capacity</td><td>800</td><td>J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>ν<sub>c</sub></td><td>Core magnetic diffusivity</td><td>2</td><td>m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>μ<sub>c</sub></td><td>Core magnetic permeability</td><td>4π 10<sup>-7</sup></td><td>H m<sup>-1</sup></td></tr><tr><td>ρ<sub>c</sub></td><td>Core density</td><td>7500</td><td>kg m<sup>-3</sup></td></tr><tr><td>B</td><td>Average magnetic field strength</td><td>Eq. (3)</td><td>T</td></tr><tr><td>R</td><td>Planet radius</td><td></td><td>m</td></tr><tr><td>R<sub>c</sub></td><td>Core radius</td><td>77% R</td><td>m</td></tr><tr><td>r<sub>i</sub></td><td>Inner core radius</td><td>1% R<sub>c</sub></td><td>m</td></tr><tr><td>G</td><td>Gravitational constant</td><td>6.67 × 10<sup>-11</sup></td><td>m<sup>3</sup> kg<sup>-1</sup> s<sup>-2</sup></td></tr><tr><td>κ<sub>Si</sub></td><td>Silicates thermal conductivity</td><td>4.32–1.70</td><td>W/m K<sup>-1</sup></td></tr><tr><td>κ<sub>Si</sub></td><td>Molten silicates thermal conductivity</td><td>0.18–0.41</td><td>W/m K<sup>-1</sup></td></tr><tr><td>α<sub>Si</sub></td><td>Silicates thermal expansion</td><td>1 × 10<sup>-5</sup></td><td>K<sup>-1</sup></td></tr><tr><td>C<sub>p, Si</sub></td><td>Silicates heat capacity</td><td>1000</td><td>J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>T<sub>s</sub></td><td>Surface temperature</td><td>500</td><td>K</td></tr><tr><td>δ<sub>Si</sub></td><td>Mantle thickness</td><td>33% R</td><td>m</td></tr></tbody></table>

Aubert (2006)). Mariner 10 spacecraft measurements showed that Mercury's magnetic field was 100 times weaker than the Earth's one. To account for this discrepancy, we consider that the average strength of the magnetic field is 1% $\bar{B}$ obtained from Eq. (32) when solving Eq. (31). We also consider that $\nu_c$, $\mu_c$ and $\rho_c$ are constants (see values in Table 4). We note that most of the power needed to overcome the criterion related to Eq. (31) can be supplied by thermal core convection (i.e. criterion related to Eq. (28)) especially for large metallic cores.

(3) The magnetic Reynolds number ($Re_m$) must be supercritical in order to have convective motions, inducing a complex structure needed to carry the magnetic field lines (U. R. Christensen and Aubert, 2006).

Reynolds magnetic number is calculated using Christensen and Aubert, 2006 formulation's and assuming that the characteristic magnetic length scale is the radius of the core:

$$Re_m = \left( \frac{Q_{CMB} G \alpha_c}{3 C_{p,c}} \right)^{1/3} \frac{R_c}{\nu_c} > 10 - 100 \quad (33)$$

The value of the critical magnetic Reynolds number ($Re_{m,c}$) is usually constrained for models implying a large inner core. For Mercury-like planets, this value is less documented. Here we make the conservative assumption that this value ranges between 10 and 100.

(4) The mantle heat flow has to overcome the core heat flow needed to generate a dynamo.

The mantle's Rayleigh number ($Ra$) conditions the efficiency of heat evacuation from the mantle:

$$Ra = \frac{\alpha_{Si} \rho_{Si} g \Delta T \delta_{Si}^3}{\eta_{Si} \kappa_{Si}} \quad (34)$$

With $\alpha_{Si}$ the mantle thermal expansion coefficient, $\rho_{Si}$ the mantle density, $\Delta T$ the temperature difference between the core and the surface of the planet, $\delta_{Si}$ the mantle thickness, $\eta_{Si}$ the mantle viscosity, $\kappa_{Si}$ the mantle heat diffusivity ($\kappa_{Si} = k_{Si}/(\rho_{Si}C_p, Si)$), and $C_p, Si$ the mantle heat capacity. If $Ra$ is lower than a critical value ($Ra_C \approx 1000$), the heat is evacuated from the mantle by conduction and the mantle heat flow $Q_{Si}$ is:

$$Q_{Si} = \frac{k_{Si} \Delta T}{\delta_{Si}} 4\pi R^2 \quad (35)$$

Alternatively, if the $Ra > Ra_C$, the heat is evacuated from the mantle by convection and $Q_{Si}$ scales with $Ra^{1/3}$ (Solomatov, 2007):

$$Q_{Si} = 0.089 \frac{k_{Si} \Delta T}{\delta_{Si}} 4\pi R^2 Ra^{1/3} \quad (36)$$

In any case, for a thermally-driven dynamo to operate, $Q_{Si}$ has to be larger than the core heat flow $Q_{Core}$. For these calculations, we consider here that the surface gravity is dominated by the gravity at the CMB ($g = g_c$). On a planet with a core over planet radius ratio ($R_c/R$) of ~0.5, the mantle is thick enough so that $Ra > Ra_c$ and the heat is efficiently evacuated by convection. On such bodies, the dynamo is generated easily and should appear early in the planet history. If the heat evacuation is very efficient, the planet might cool down rapidly, which can potentially lead to a short-lived dynamo. On a Mercury-like body with a large $R_c/R$ ratio, the mantle is thin compared to the core and its Rayleigh number, critical criterion for thermal convection, hardly overcomes $Ra_c$. Hence, conduction should be the main heat evacuation process within the mantle of Mercury-like planets, which can limit the occurrence of an early thermally-driven dynamo.
## 5.2. Simulation results

We now compare the core heat flow needed to generate a dynamo within a Mercury-like body (with e.g. $R_c/R_{planets}$) (Figs. 11 and 12). To do this, we use a maximum possible range of solid conductivities between $\kappa_{Si} = 4.32$ and $1.70$ W/m/K. We consider a maximum planet radius $R = 2440$ km with a thin mantle shell to stay in temperature and pressure conditions equivalent to the experimental conditions detailed above ($P \le 2$GPa and $T \le 1700$ K). We also report the critical size of the protoplanet ($R = 215$ km) above which $Re_m$ is larger than 100 (Eq. (33)). The Fig. 11 shows that a solid conductive mantle is able to evacuate enough heat to induce a thermally driven dynamo if the planet radius is larger than ~215 km, regardless of the value of thermal conductivity in the outer core, at least up to more than 40 W/m/K. In contrast, the change of core conductivities from e.g. 40 W/m/K (Silber et al., 2019) to 4 W/m/K (Manthilake et al., 2019) has a major effect on the maximum size of a planet that can operate a thermally driven dynamo. For a core conductivity $\kappa_{core} = 40$ W/m/K, the maximum planet radius evolves from 1220 to 1880 km, for $\kappa_{Si} = 1.70$ and $4.32$ W/m/K, respectively. When $\kappa_{core} = 4$ W/m/K, the maximum planet radius is more than 2500 km.

We now perform the same type of calculation for a mushy mantle using the silicate-melt conductivity refined in this study (Tables 3 and 4). We consider here neither a complete magma ocean that would evacuate the core heat very efficiently by turbulent convection (Monteux et al., 2016) nor a mushy mantle with a large melt fraction (larger than 40–60%) which would imply a strong decrease of the bulk viscosity

of the mantle (Picard et al., 2013). We consider here a purely conductive mantle where the melt fraction is smaller than 20%. Hence, if this melt fraction is concentrated at the core mantle boundary, the molten layer is thin enough to avoid its convection, while if the liquid is equally distributed within the mantle, we can reasonably assume that the bulk viscosity is weakly affected by the liquid phase and close to the viscosity of subsolidus silicate material. Still, the presence of melt may affect mantle properties. An important parameter is the fate of melts, which is primarily controlled by the solid-melt density contrast. The molten material can either be evacuated at the surface during early volcanism or cumulate at the lowermost mantle. For peridotites, negative melt buoyancy is unlikely to happen at pressures lower than 7 GPa (Sakamaki et al., 2006; Matsukage et al., 2005; Freitas et al., 2017). However, the buoyancy of melts generated from the partial melting of a body with composition largely different from peridotite remains uncertain. Molten reservoirs may be trapped at shallow depths at the end of mushy mantle cooling (Monteux et al., 2020) or during the last stages of fractional crystallization within small bodies (e.g. Frossard et al., 2019). On a Moon-like body, a melt layer may be trapped below an anorthositic crust while on a Mercury-like body, anorthite is denser than the melt and then may not float allowing the melt to form a late shallow magma ocean. We consider here that the melt has neutral buoyancy and is equally distributed within the planetary mantle, with the consequence that the mantle viscosity is primarily controlled by the solid fraction of the mantle.

Due to the lower thermal conductivity of the melt compared to the solid mantle, the mushy mantle should limit the heat flow that can be

<figure><img src="image_269.png" /><figcaption>**Fig. 11.** Heat flows as a function of the planetary radius for a fully solid planetary mantle. The black and grey lines represent the critical heat flow needed for a thermally driven dynamo for $k_c = 40$ and 4 W/m/K, respectively (Conditions 1 and 2, Eq. (31)). The dark and light green lines represent the mantle heat flow in a conductive regime for $\kappa_{Si} = 4.32$ and $1.70$ W/m/K, respectively (Eq. (35)). Green triangles represent the critical planet size at which $Ra = R_{ac}$, therefore corresponding to the minimum heat flow at the CMB, for mantle viscosity ranging from $10^{18}$ to $10^{22}$ Pa.s. Extrapolated values at $R = R_{Mercury} = 2440$ km are also represented with green circles for conductive heat transfers and red and pink squares for convective cases (Eq. (36)) with mantle viscosities indicated nearby. The red dashed line represents the critical value for $Rem = 100$ (Condition 3). At the left side of this line, a dynamo is unlikely. The likely zone for a planet radius capable of powering a dynamo is represented by the pink shaded area in case $\kappa_{Si} = 4.32$ W/m/K and $\kappa_{core} = 40$ W/m/K. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_270.png" /><figcaption>**Fig. 12.** Heat flows as a function of the planetary radius for a stagnant partially molten planetary mantle. The black and grey lines represent the critical heat flow for a thermally driven dynamo for $k_c = 40$ and 4 W/m/K, respectively (Conditions 1 and 2, Eq. (31)). The dark and light green lines represent the mantle heat flow for $κ_{Si} = 0.41$ and 0.18 W/m/K, respectively (Eq. (35)). Green triangles represent the critical planet size at which $Ra = Ra_c$, for mantle viscosity ranging from $10^{18}$ to $10^{22}$ Pa.s. Extrapolated values at $R = R_{Mercury} = 2440$ km are also represented with green circles for conductive heat transfers and red and pink squares for convective cases (Eq. (36)) with viscosities indicated nearby. The red dashed line represents the critical value for $Rem = 100$ (Condition 3). At the left side of this line, a dynamo is unlikely. The likely zone for a planet radius capable of powering a dynamo is represented by the pink shaded area in case $κ_{Si} = 0.41$ W/m/K and $κ_{core} = 40$ W/m/K. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

extracted from the core. In our calculation, we consider for simplicity that the entire mantle has the same conductivity as the melt. As a result, the range of planetary radii for which a dynamo is plausible is much narrower than for a fully solid mantle (Fig. 12). For a core conductivity of 40 W/m/K, a thermally driven dynamo can only operate for planetary sizes ranging between 215 and 350–570 km, depending on the value of $κ_{Si}$ and for $κ_{core} = 40$ W/m/K. The maximum planetary radius becomes 1000–1500 km for a core conductivity of 4 W/m/K. Therefore, the melting of a Mercury-like planetary mantle could limit the establishment of a magnetic field on the planet. This effect could last over different timescales, depending on the size of the planet and of the core, the initial core temperature, the fate of the melt material, etc.

The differences obtained for fully solid (Fig. 11) and molten (Fig. 12) models show that the presence of melt within a “Mercury-like” planet mantle could result into better thermal insulation of the core, which may prevent the occurrence of a thermally driven dynamo. Of course, planets presenting mantle with supercritical Rayleigh number would not be affected by this effect, because mantle convection could extract enough heat from the core at the CMB. As a consequence, the solid mantle viscosity plays a key role because (1) it controls the planet size above which mantle convection is likely, through the Ra value, and (2) it governs the efficiency of mantle heat evacuation once the critical Rayleigh number is overcome (see Eq. (36)). Decreasing the mantle viscosity from $10^{22}$ to $10^{18}$ Pa.s decreases the critical planet size for mantle convection by one order of magnitude (from 1600 to 1900 km to 160–200 km for a solid mantle and from 850 to 1050 km to 80–110 km for a mushy mantle, see Figs. 11 and 12). The figures also show that the lower is mantle viscosity, the easier the planet can meet the criteria for a thermally driven dynamo

when the mantle of the planet is convecting. We note here that considering the influence of melt on the bulk mantle viscosity should lead to a decrease of the critical Rayleigh number and as a consequence of the critical planetary size where heat in only evacuated by thermal conduction (i.e. green triangles would move to the left). However, this effect should be negligible for the small melt fraction considered here.

The impact of mantle partial melting crucially depends on the melt quantity and its location within the mantle. Whether it is distributed in a mushy mantle or accumulated as a pond at a given mantle depth affects both, locally and globally, the mantle viscosity and the thermal diffusivity. As a consequence, the dynamic regime can evolve from conduction to convection, changing the thermal outputs at the CMB and making the heat flux enough, or insufficient, to power a thermally-driven dynamo. In this framework, transient phenomena such as mantle overturns (Elkins-Tanton et al., 2003; Ballmer et al., 2017) implying major melt migration will result in a brutal change in the planet magnetic signature. Similarly, the chemical stratification resulting in the production of a crust (anorthositic, basaltic etc.) may help insulating the planet interior in its early history suggesting weaker but longer sustained dynamos.

## 6. Conclusions

In this study, we constrain the thermal diffusivities of silicate melts with a wide range of composition at high pressure and temperature conditions using the multi-anvil apparatus and the Angström method. We observed that melting induces a significant decrease in both thermal diffusivities and conductivities. We measure melt thermal diffusivities ranging from 0.18 (±0.02) W/m/k to 0.41 (±0.04) W/m/K. The
composition of the melt does not have a significant impact on diffusiv-
ities. We then use our results to address the thermal properties of the
molten silicate reservoirs, and the likeliness of a thermally driven dy-
namo within a Mercury-like planets presenting a large metallic core and
a relatively thin silicate mantle (large Rc/Rplanet). The mantle of such
bodies is not expected to be convecting easily (Ra < Ra_c) leading to an
inefficient heat evacuation and difficulty generating a dynamo
compared to planets with a thicker mantle (where easily Ra > Ra_c). Our
results illustrate that the presence of a mushy mantle above the CMB in
such a Mercury-like planet can significantly reduce the ability of the
mantle to evacuate the heat from the core and limit the likeliness of a
thermally-driven dynamo. As a consequence, the fate and the lifetime of
such a mushy reservoir can have a profound impact on the thermal
history of Mercury-like planets. Future observations of Mercury-like
planets accreted in orbits close to their star and the eventual detection
of their intrinsic magnetic field would in return constrain the state of the
mantle surrounding their metallic core.

Author contributions

**Conceptualization:** DF and GM. **Data curation:** DF, FS, JM, DA and GM. **Resources & Software:** AM, NC, GM, DA and DF. **Writing - original draft:** DF, JM and DA. **Writing - review & editing:** DF, JM, DA, GM and FS. **Funding acquisition:** GM, DA and DF.

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.

Acknowledgments

We acknowledge D. Laporte, A. Gourgeau, O. Sigmarsson, J-L. Froger and S. Jouhannel for providing the natural samples used in this study, F. Pointud, C. Guillot and J.L. Fruquière for their technical assistance with high pressure apparatus and sample/assembly preparation, J-M Henot and E. Voyer for the SEM analyses, J-L Devidal for the electron microscope analyses, C. Fonquernie for flash analyses and M. Benbakkar for ICP-AES measurements. This is Laboratory of Excellence ClerVolc contribution number 452.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.
org/10.1016/j.pepi.2021.106655.

References

Acuna, M.H., Connerney, J.E.P., Lin, R.P., Mitchell, D., Carlson, C.W., McFadden, J., et al., 1999. Global distribution of crustal magnetization discovered by the Mars global surveyor MAG/ER experiment. Science 284 (5415), 790-793.

Anderson, B.J., Johnson, C.L., Korth, H., Purucker, M.E., Winslow, R.M., Slavin, J.A., et al., 2011. The global magnetic field of mercury from MESSENGER orbital observations. Science 333 (6051), 1859-1862. https://doi.org/10.1126/science.1211001.

Anzellini, S., Dewaele, A., Mezouar, M., Loubeyre, P., Morard, G., 2013. Melting of iron at Earth's inner core boundary based on fast X-ray diffraction. Science 340, 464-466.

Asphaug, E., Reufer, A., 2014. Mercury and other iron-rich planetary bodies as relics of inefficient accretion. Nat. Geosci. 7 (8), 564–568. https://doi.org/10.1038/ngeo2189.

Ballmer, M.D., Lourenco, D.L., Hirose, K., Caracas, R., Nomura, R., 2017. Reconciling magma-ocean crystallization models with the present-day structure of the Earth's mantle. Geochem. Geophys. Geosyst. 18 (7), 2785–2806. https://doi.org/10.1002/2017gc006917.

Barrat, J.A., Gillet, P.H., Lesourd, M., Blicher-Toft, J., Poupeau, G.R., 1999. The Tatahouine diogenite: mineralogical and chemical effects of sixty-three years of terrestrial residence. Meteorit. Planet. Sci. 34, 91–97.

Beck, A.E., Darbha, D.M., Schloessin, H.H., 1978. Lattice conductivities of single-crystal and polycrystalline materials at mantle pressures and temperatures. Phys. Earth Planet. Inter. 17 (1), 35–53. https://doi.org/10.1016/0031-9201(78)90008-0.

Benz, W., Slattery, W.L., Cameron, A.G.W., 1988. Collisional stripping of Mercury's mantle. Icarus 74 (3), 516-528. https://doi.org/10.1016/0019-1035(88)90118-2.

Bolfan-Casanova, N., Montagnac, G., Reynard, B., 2014. Measurement of water contents in olivine using Raman spectroscopy. Am. Mineral. 99, 149–156. https://doi.org/10.2138/am.2014.4444.

Boujibar, A., Andrault, D., Ali, M., Bolfan-casanova, N., Devidal, J., Trcera, N., 2014. Metal – silicate partitioning of sulphur, new experimental and thermodynamic constraints on planetary accretion. Earth Planet. Sci. Lett. 391, 42–54. https://doi.org/10.1016/j.epsl.2014.01.021.

Boujibar, A., Andrault, D., Bolfan-Casanova, N., Bouhifd, M.A., Monteux, J., 2015. Cosmochemical fractionation by collisional erosion during the Earth's accretion. Nat. Commun. 6, 1–7. https://doi.org/10.1038/ncomms9295.

Branlund, J.M., Hofmeister, A.M., 2012. Heat transfer in plagioclase feldspars. Am. Mineral. 97, 1145–1154.

Breuer, D., Hauck, S.A., Buske, M., Pauer, M., Spohn, T., 2007. Interior evolution of mercury. Space Sci. Rev. 132, 229–260.

Breuer, D., Rueckriemen, T., Spohn, T., 2015. Iron snow, crystal floats, and inner-core growth: modes of core solidification and implications for dynamos in terrestrial planets and moons. Progr. Earth Planet. Sci. 2 (1), 39. https://doi.org/10.1186/s40645-015-0069-y.

Buffett, B.A., 2002. Estimates of heat flow in the deep mantle based on the power requirements for the geodynamo. Geophys. Res. Lett. 29 (12), 1–4.

Busek, P.R., 1997. Pallasite meteorites-mineralogy, petrology and geochemistry. Meteorit. Cosmochim. Acta 41, 711–740.

Carslaw, H.S., Jaeger, J.C., 1959. Heat conduction in solids. Oxford Univ. Press.

Charlier, B., Namur, O., 2019. The origin and differentiation of planet mercury. Elements 15 (1), 9–14. https://doi.org/10.2138/gselements.15.1.9.

Christensen, U.R., Aubert, J., 2006. Scaling properties of convection-driven dynamos in rotating spherical shells and application to planetary magnetic fields. Geophys. J. Int. 166 (1), 97–114. https://doi.org/10.1111/j.1365-246X.2006.03009.x.

Deng, L., Seagle, C., Fei, Y., Shahar, A., 2013. High pressure and temperature electrical resistivity of iron and implications for planetary cores. Geophys. Res. Lett. 40, 33–37. https://doi.org/10.1029/2012GL054347.

Ebel, D.S., Grossman, L., 2000. Condens. Dust-Enrich. Syst. 64 (2), 339–366.

Elkins-Tanton, L.T., Parmentier, E.M., Hess, P.C., 2003. Magma Ocean fractional crystallization and cumulate overturn in terrestrial planets: implications for Mars. Meteorit. Planet. Sci. 38 (12), 1753–1771. https://doi.org/10.1111/j.1945-5100.2003.tb00013.x.

Fegley, B., Cameron, A.G.W., 1987. A vaporization model for iron/silicate fractionation in the mercury protoplanet. Earth Planet. Sci. Lett. 82 (3-4), 207-222. https://doi.org/10.1016/0012-821X(87)90196-8.

Floran, J.R., Prinz, M., Hlava, P.F., Keil, K., Nehru, C.E., Hinthorne, J.R., 1978. The Chassigny meteorite: a cumulate dunite with hydrous amphibole-bearing melt inclusions. Geochem. Comsochim. Acta 42 (8), 1213–1219. https://doi.org/10.1016/0016-7037(78)90115-1.

Freitas, D., Manthilake, G., 2019. Electrical Conductivity of Hydrous Silicate Melts: Implications for the Bottom-Up Hydration of Earth's Upper Mantle EPSL 523, 1–9. https://doi.org/10.1016/j.epsl.2019.115712.

Freitas, D., Manthilake, G., Schiavi, F., Chantel, J., Bolfan-Casanova, N., Bouhifd, M.A., Andrault, D., 2017. Experimental evidence supporting a global melt layer at the base of the Earth's upper mantle. Nat. Commun. 8, 2186. https://doi.org/10.1038/s41467-017-02275-9.

Freitas, D., Manthilake, G., Chantel, J., Bouhfid, M., Andrault, D., 2019. Simultaneous measurements of electrical conductivity and seismic velocity of partially molten geological materials: implications for melt fraction in the upper mantle. Phys. Chem. Miner. https://doi.org/10.1007/s00269-019-01021-5, 0, 0.

Frossard, P., Boyet, M., Bouvier, A., Hammouda, T., Monteux, J., 2019. Evidence for anorthositic crust formed on an inner solar system planetesimal. Geochem. Perspect. Lett. 11, 28–32. https://doi.org/10.7185/geochemlet.1921.

Fujisawa, H., Fujii, N., Mizutani, H., Kanmori, H., Akimoto, S., 1968. Thermal diffusivity of Mg2SiO4, Fe2SiO4, and NaCl at high pressures and temperatures. J. Geophys. Res. 73, 4727.

Garrick-Bethell, I., Weiss, B.P., Shuster, D.L., Buz, J., 2009. Early lunar magnetism. Science 323 (5912), 356–359. https://doi.org/10.1126/science.1166804.

Genova, A., Goossens, S., Mazarico, E., Lemoine, F.G., Neumann, G.A., Kuang, W., et al., 2018. Geodetic Evidence that Mercury Has a Solid Inner Core, 1–13. https://doi.org/10.1029/2018GL081135.

Gibert, B., Schilling, F.R., Gratz, K., Tommasi, A., 2005. Thermal diffusivity of olivine single crystals and a dunite at high temperature: Evidence for heat transfer by radiation in the upper mantle. Phys. Earth. Planet. Inter. 151, 129–141. https://doi.org/10.1016/j.pepi.2005.02.003.

Giordano, D., Nichols, A.R.L., Dingwell, D.B., 2005. Glass transition temperatures of natural hydrous melts: A relationship with shear viscosity and implications for the welding process. J. Volcanol. Geotherm. Res. 142 (1–2 SPEC ISS), 105–118. https://doi.org/10.1016/j.jvolgeores.2004.10.015.

Grant, F.S., West, G.F., 1965. Interpretation Theory in Applied Geophysics. McGraw-Hill Book Co., NY, p. 583.

Hauck, S.A., Margot, J.L., Solomon, S.C., Phillips, R.J., Johnson, C.L., Lemoine, F.G., et al., 2013. The curious case of Mercury's internal structure. J. Geophys. Res. E: Planets 118 (6), 1204–1220. https://doi.org/10.1002/jgre.20091.

Hernlund, J., Leinenweber, K., Locke, D., Tyburczy, J.A., 2006. A numerical model for steady-state temperature distributions in solid-medium high-pressure cell assemblies. Am. Mineral. 91 (2–3), 295–305. https://doi.org/10.2138/am.2006.1938.
Hofmeister, A.M., Branlund, J.M., 2015. Thermal conductivity of the earth. In: Treatise on Geophysics: Second Edition (Vol. 2). Elsevier B.V. https://doi.org/10.1016/B978-0-444-53802-4-00047-6.

Hofmeister, A.M., Pertermann, M., 2008. Thermal diffusivity of clinopyroxenes at elevated temperature. Eur. J. Mineral. 20, 537–549.

Hofmeister, Anne M., Whittington, A.G., Pertermann, M., 2009. Transport properties of high albite crystals, near-endmember feldspar and pyroxene glasses, and their melts to high temperature. Contrib. Mineral. Petrol. 158 (3), 381–400. https://doi.org/10.1007/s00410-009-0388-3.

Hofmeister, A.M., Branlund, J.M., Pertermann, M., 2007. Properties of rocks and minerals - thermal conductivity of the earth. Treatise Geophys. 2, 543–577. https://doi.org/10.1016/B978-044452748-6-00048-1.

Hofmeister, Anne M., Sehlke, A., Whittington, A.G., 2014. Thermal diffusivity of Fe-rich pyroxene glasses and their melts. Chem. Geol. 384, 1–9. https://doi.org/10.1016/j.chemgeo.2014.06.018.

Hood, L.L., Harrison, K.P., Langlais, B., Lillis, R.J., Poulet, F., Williams, D.A., 2010. Magnetic anomalies near Apollinaris Patera and the Medusae fossae formation in Lucus Planum, Mars. Icarus 208 (1), 118–131. https://doi.org/10.1016/j.icarus.2010.01.009.

Johnson, C.L., Phillips, R.J., Purucker, M.E., Anderson, B.J., Byrne, P.K., Denevi, B.W., et al., 2015. Low-altitude magnetic fieldmeasurements by MESSENGER reveal Mercury's ancient crustal field. Science 348 (6237), 892–895. https://doi.org/10.1126/science.aaa8720.

Johnson, Catherine L., Phillips, R.J., 2005. Evolution of the Tharsis region of Mars: insights from magnetic field observations. Earth Planet. Sci. Lett. 230 (3-4), 241–254. https://doi.org/10.1016/j.epsl.2004.10.038.

Kabin, K., Heimpel, M.H., Rankin, R., Aurnou, J.M., Gómez-Pérez, N., Paral, J., et al., 2008. Global MHD modeling of Mercury's magnetosphere with applications to the MESSENGER mission and dynamo theory. Icarus 195 (1), 1–15. https://doi.org/10.1016/j.icarus.2007.11.028.

Kanamori, H., Mizutani, H., Fujii, N., 1969. Method of thermal diffusivity measurement. J. Phys. Earth 17, 43–57.

Katsura, T., 1993. Thermal diffusivity of silica glass at pressures up to 9 GPa. Phys. Chem. Miner. 20 (3), 201–208. https://doi.org/10.1007/BF00200122.

Keil, K., 2010. Enstatite achondrite meteorites (aubrites) and the histories of their asteroidal parent bodies. Chem. Erde 70 (4), 295–317. https://doi.org/10.1016/j.chemer.2010.02.002.

Laporte, D., Provost, A., 2000. Equilibrium geometry of a fluid phase in a polycrystalline aggregate with anisotropic surface energies: dry grain boundaries. J. Geophys. Res. Solid Earth 105, 25937–25953. https://doi.org/10.1029/2000jb900256.

Laporte, D., Rapaille, C., Provost, A., 1997. Wetting angles, equilibrium melt geometry, and the permeability threshold of partially molten crustal protoliths BT. In: Bouchez, J.L., Hutton, D.H.W., Stephens, W.E. (Eds.), Granite: From Segregation of Melt to Emplacement Fabrics. Springer Netherlands, Dordrecht, pp. 31–54.

Lawrence, K., Johnson, C., Tauxe, L., Gee, J., 2008. Lunar paleointensity measurements: implications for lunar magnetic evolution. Phys. Earth Planet. Inter. 168 (1-2), 71–87.

Lillis, R.J., Frey, H.V., Manga, M., Mitchell, D.L., Lin, R.P., Acuña, M.H., Bougher, S.W., 2008. An improved crustal magnetic field map of Mars from electron reflectometry: Highland volcano magmatic history and the end of the martian dynamo. Icarus 194 (2), 575–596. https://doi.org/10.1016/j.icarus.2007.09.032.

Lillis, R.J., Robbins, S., Manga, M., Halekas, J.S., Frey, H.V., 2013. Time history of the Martian dynamo from crater magnetic field analysis. J. Geophys. Res. E: Planets 118 (7), 1488–1511. https://doi.org/10.1002/jgre.20105.

Malavergne, V., Toplis, M.J., Berthet, S., Jones, J., 2010. Highly reducing conditions during core formation on mercury: implications for internal structure and the origin of a magnetic field. Icarus 206 (1), 199–209. https://doi.org/10.1016/j. icarus.2009.09.001.

Manthilake, G., Chantel, J., Monteux, J., Andrault, D., Bouhifd, M.A., Casanova, N.B., Boulard, E., Guignot, N., King, A., Itie, J.P., 2019. Thermal conductivity of FeS and its implications for Mercury's long sustaining magnetic field. J. Geophys. Res. E Planets. https://doi.org/10.1029/2019JE005979.

Manthilake, G.M., De Koker, N., Frost, D.J., McCammon, C.A., 2011b. Lattice thermal conductivity of lower mantle minerals and heat flux from Earth's core. Proc. Natl. Acad. Sci. U. S. A. 108 (44), 17901–17904. https://doi.org/10.1073/ pnas.1110594108.

Manthilake, M.A.G.M., De Koker, N., Frost, D.J., 2011a. Thermal conductivity of CaGeO3 perovskite at high pressure. Geophys. Res. Lett. 38 (8), 3–6. https://doi.org/10.1029/2011GL046882.

Matsukage, K.N., Jing, Z., Karato, S., 2005. Density of hydrous silicate melt at the conditions of Earth's deep upper mantle. Nature 438, 488–491. https://doi.org/10.1038/nature04241.

Maumus, J., Bagdassarov, N., Schmeling, H., 2005. Electrical conductivity and partial melting of mafic rocks under pressure. Geochim. Cosmochim. Acta 69, 4703–4718. https://doi.org/10.1016/j.gca.2005.05.010.

Médard, E., McCammon, C.A., Barr, J.A., Grove, T.L., 2008. Oxygen fugacity, temperature reproducibility, and H2O contents of nominally anhydrous piston-cylinder experiments using graphite capsules. Am. Mineral. 93 (11–12), 1838–1844. https://doi.org/10.2138/am.2008.2842.

Mibe, K., Fujii, T., Yasuda, A., 1998. Connectivity of aqueous fluid in the Earth's upper mantle. Geophys. Res. Lett. 25, 1233–1236.

Mibe, K., Fujii, T., Yasuda, A., 1999. Control of the location of the volcanic front in island arcs by aqueous fluid connectivity in the mantle wedge. Nature 401, 259–262. https://doi.org/10.1038/45762.

Mighani, S., Wang, H., Shuster, D.L., Borlina, C.S., Nichols, C.I.O., Weiss, B.P., 2020. The end of the lunar dynamo. Science 1–8.

Minarik, W.G., Watson, E.B., 1995. Interconnectivity of carbonate melt at low melt fraction. Earth Planet. Sci. Lett. 133 (3-4), 423-437. https://doi.org/10.1016/0012-821X(95)00085.

Monteux, J., Arkani-Hamed, J., 2014. Consequences of giant impacts in early Mars: Core merging and Martian dynamo evolution. J. Geophys. Res. 119, 480–505.

Monteux, J., Jellinek, A.M., Johnson, C.L., 2011. Why might planets and moons have early dynamos? Earth Planet. Sci. Lett. 310 (3-4), 349–359. https://doi.org/10.1016/j.epsl.2011.08.014.

Monteux, J., Andrault, D., Samuel, H., 2016. On the cooling of a deep terrestrial magma ocean. Earth Planet. Sci. Lett. 448, 140–149. https://doi.org/10.1016/j.epsl.2016.05.010.

Monteux, J., Andrault, D., Guitreau, M., Samuel, H., Demouchy, S., 2020. A mushy Earth's mantle for more than 500 Myr after the magma ocean solidification. Geophys. J. Int. https://doi.org/10.1093/gji/gga064.

Namur, O., Collinet, M., Charlier, B., Grove, T.L., Holtz, F., McCammon, C., 2016. Melting processes and mantle sources of lavas on mercury. Earth Planet. Sci. Lett. 439, 117–128. https://doi.org/10.1016/j.epsl.2016.01.030.

Neuville, D.R., Courtial, P., Dingwell, D.B., Richet, P., 1993. Thermodynamic and rheological properties of silicate melts. Contrib. Mineral. Petrol. 113, 572–581.

Ni, H., Keppler, H., Behrens, H., 2011. Electrical conductivity of hydrous basaltic melts: implications for partial melting in the upper mantle. Contrib. Mineral. Petrol. 162 (3), 637–650. https://doi.org/10.1007/s00410-011-0617-4.

Ni, H., Hui, H., Steinle-Neumann, G., 2015. Transport properties of silicate melts. Rev. Geophys. 53, 715–744. https://doi.org/10.1002/2015RG000485.

Nimmo, F., 2002. Why does Venus lack a magnetic field? Geology 30 (11), 987-990. https://doi.org/10.1130/0091-7613(2002)030-0987:WDVLAM>2.0.CO;2.

O'Rourke, J.G., Buz, J., Fu, R.R., Lillis, R.J., 2010. Detectability of Remanent magnetism in the crust of Venus. Geophys. Res. Lett. 46 (11), 5768–5777. https://doi.org/10.1029/2010GL082725.

Osako, M., Ito, E., Yoneda, A., 2004. Simultaneous measurements of thermal conductivity and thermal diffusivity for garnet and olivine under high pressure. Earth Planet. Inter. 143, 311–320. https://doi.org/10.1016/j.pepi.2003.10.

Pertermann, M., Hofmeister, A.M., 2006. Thermal diffusivity of olivine-group minerals at high temperature. Am. Mineral. 91, 1747–1760.

Picard, D., Arbaret, L., Pichavant, M., Champallier, R., Launeau, P., 2013. The rheological transition in plagioclase-bearing magmas. J. Geophys. Res. Solid Earth 118 (4), 1363–1377.

Pommier, A., Leinenweber, K., Tran, T., 2019. Mercury's thermal evolution controlled by an insulating liquid outermost core? Earth Planet. Sci. Lett. 517, 125–134. https://doi.org/10.1016/j.epsl.2019.04.022.

Romine, W.L., Whittington, A.G., Nabelek, P.I., Hofmeister, A.M., 2012. Thermal diffusivity of rhyolitic glasses and melts: effects of temperature, crystals and dissolved water. Bull. Volcanol. 74 (10), 2273–2287. https://doi.org/10.1007/s00445-012-0661-6.

Sakamaki, T., Suzuki, A., Ohtani, E., 2006. Stability of hydrous melt at the base of the Earth's upper mantle. Nature 439 (7073), 192–194. https://doi.org/10.1038/nature04352.

Sambridge, M., Mosegaard, K., 2002. Monte Carlo methods in geophysical inverse problems. Rev. Geophys. 40 (3), 1–3.

Schiavi, F., Bolfan-Casanova, N., Withers, A.C., Médard, E., Laumonier, M., Laporte, D., et al., 2018. Water quantification in silicate glasses by Raman spectroscopy: correcting for the effects of confocality, density and ferric iron. Chem. Geol. 483 (February), 312–331. https://doi.org/10.1016/j.chemgeo.2018.02.036.

Schmeling, H., 1986. Numerical models on the influence of partial Melton elastic, anelastic and electrical properties of rocks. Part II: electrical conductivity. Phys. Earth Planet. Inter. 43 (2), 123–136. https://doi.org/10.1016/0031-9201(86)90080-4.

Shea, T., Houghton, B.F., Gurioli, L., Cashman, K.V., Hammer, J.E., Hobden, B.J., 2010. Textural studies of vesicles in volcanic rocks: an integrated methodology. J. Volcanol. Geotherm. Res. 190 (3-4), 271–289. https://doi.org/10.1016/j.jvolgeores.2009.12.003.

Smith, D.S., Fayette, S., Grandjean, S., Martin, C., Telle, R., Tonnessen, T., 2003. Thermal resistance of grain boundaries in alumina ceramics and refractories. J. Am. Ceram. Soc. 86 (1), 105–111. https://doi.org/10.1111/j.1151-2916.2003.tb03285.x.

Silber, R.E., Secco, R.A., Yong, W.J., Littleton, J.A.H., 2019. Heat Flow in Earth's Core From Invariant Electrical Resistivity of Fe-Sion the Melting Boundary to 9 GPa: Do Light Elements Matter? J. Geophys. Res.: Solid Earth 124 (6), 5521–5543. https://doi.org/10.1029/2019JB017375.

Smith, D.S., Alzina, A., Bourret, J., Nait-Ali, B., Pennek, F., Tessier-Doyen, N., Otsu, K., Matsubara, H., Elser, P., Gonzenbach, U.T., 2013. Thermal conductivity of porous materials. J. Mater. Res. 28 (17).

Solomatov, V., 2007. Magma oceans and primordial mantle differentiation. In: Schubert, G. (Ed.), Treatise of Geophysics, 9. Elsevier.

Stevenson, D.J., Spohn, T., Schubert, G., 1983. Magnetism and thermal evolution of the terrestrial planets. Icarus 54 (3), 466–489. https://doi.org/10.1016/0019-1035(83)90241-5.

Strom, R.G., Sprague, A.L., 2003. Exploring Mercury: The Iron Planet. Springer Science & Business Media.

Tarduno, J.A., Cottrell, R.D., Watkeys, M.K., Hofmann, A., Doubrovine, P.V., Mamajek, E.E., Liu, D., Sibeck, D.G., Neukirch, L.P., Usui, Y., 2010. Geodynamo, solar wind, and magnetopause 3.4 to 3.45 billion years ago. Science 327 (5970), 1238–1240.

ten Grotenhuis, S.M., Drury, M.R., Spiers, C.J., Peach, C.J., 2005. Melt distribution in olivine rocks based on electrical conductivity measurements. J. Geophys. Res. Solid Earth 110, 1–11. https://doi.org/10.1029/2004JB003462.
Tian, Z., Zuber, M.T., Stanley, S., 2015. Magnetic field modeling for mercury using dynamo models with a stable layer and laterally variable heat flux. *Icarus* 260, 263–268. https://doi.org/10.1016/j.icarus.2015.07.019.

Tikoo, S.M., Weiss, B.P., Cassata, W.S., Shuster, D.L., Gattacceca, J.Ô., Lima, E.A., et al., 2014. Decline of the lunar core dynamo. *Earth Planet. Sci. Lett.* 404, 89–97. https://doi.org/10.1016/j.epsl.2014.07.010.

Tyburczy, J.A., Waff, H.S., 1983. Electrical conductivity of molten basalt and andesite to 25 kilobars pressures: geophysical significance and implications for charge transport and melt structure. *J. Geophys. Res.* 88 (2), 2413–2430.

Von Bargen, N., Waff, H.S., 1986. Permeabilities, interfacial areas and curvatures of partially molten systems: results of numerical computations of equilibrium microstructures. *J. Geophys. Res.* 91, 9261–9276.

Waff, H.S., 1974. Theoretical consideration of electrical conductivity in a partially molten mantle and implications for geothermometry. *J. Geophys. Res.* 79 (26), 4003–4010.

Watanabe, T., Kurita, K., 1993. The relationship between electrical conductivity and melt fraction in a partially molten simple system: Archie's law behavior. *Phys. Earth Planet. Inter.* 78 (1–2), 9–17.

Watson, E.B., Brenan, J.M., Baker, D.R., 1991. Continental Mantle. Oxford University Press.

Weidenschilling, S.J., 1978. Iron/silicate fractionation and the origin of mercury. *Icarus* 35 (1), 99–111. https://doi.org/10.1016/0019-1035(78)90064-7.

Weiss, B.P., Berdal, James S., Elkins-Tanton, L., Stanley, S., Lima, E.A., Carporzen, L., 2008. Magnetism on the Angrite differentiation of Planetesimals. *Sci. Rep.* 322 (October), 713–716.

Wurm, G., Trieloff, M., Rauer, H., 2013. Photophoretic separation of metals and silicates: the formation of mercury-like planets and metal depletion in chondrites. *Astrophys. J.* 769 (1) https://doi.org/10.1088/0004-637X/769/1/78.

Xu, Y., Shankland, T.J., Linhardt, S., Rubie, D.C., Langenhorst, F., Klasinski, K., 2004. Thermal diffusivity and conductivity of olivine, wadsleyite and ringwoodite to 20 GPa and 1373 K. *Phys. Earth Planet. Inter.* 143 (1–2), 321–336. https://doi.org/10.1016/j.pepi.2004.03.005.

Yoshino, T., Takei, Y., Wark, D.A., Watson, E.B., 2005. Grain boundary wetness of texturally equilibrated rocks, with implications for seismic properties of the upper mantle. *J. Geophys. Res. B Solid Earth* 110 (8), 1–16. https://doi.org/10.1029/2004JB003544.

Yoshino, T., Nishihara, Y., Karato, S.I., 2007. Complete wetting of olivine grain boundaries by a hydrous melt near the mantle transition zone. *Earth Planet. Sci. Lett.* 256, 466–472. https://doi.org/10.1016/j.epsl.2007.02.002.

Zhang, B., Ge, J., Xiong, Z., Zhai, S., 2019. Effect of water on the thermal properties of olivine with implications for lunar internal temperature. *J. Geo. Res. Planets* 124 (12), 3469–3481. https://doi.org/10.1029/2019JE006194.
# 6 Dynamique de refroidissement du manteau primitif ter-
restre

Le manteau silicaté constitue l'essentiel du volume de la Terre (83%) et de sa masse (66%). Actuellement, son immense inertie thermique contrôle le lent refroidissement de la planète et maintient la tectonique des plaques depuis au moins 2 milliards d'années et peut-être dès le refroidissement de la surface, quelques centaines de millions d'années après la formation de la Terre, il y a 4.54 milliards d'années. Mais ce refroidissement n'a pas toujours été aussi lent et des épisodes de fusion globale associés à des dynamiques turbulentes ont pu avoir lieu tôt dans l'histoire du manteau terrestre. Les causes de la tectonique des plaques qui n'a été observée sur aucune des planètes et satellites du système solaire, restent encore largement mystérieuses et pourraient bien être le résultat d'une dynamique complexe initiée depuis la formation du système Terre-Lune. Cette section est largement inspirée du chapitre le l'encyclopédie ISTE qui sera publiée prochainement (Monteux and Andrault, in press).

## 6.1 Le manteau au stade océan magmatique

Durant les dernières étapes d'accrétion, le manteau terrestre a probablement subi plusieurs épisodes de fusion plus ou moins complète appelés "stade océan magmatique". En effet, la formation du noyau terrestre s'est déroulée très tôt dans son évolution (entre 30 et 50 Ma) (Kleine et al., 2002) ce qui nécessite un important taux de fusion pour accommoder une séparation fer/silicates rapide. De plus, les données géochimiques montrent qu'une croûte a pu se former très tôt (entre 30 et 200 Ma) comme on peut l'attendre à la surface d'un océan magmatique (Boyet and Carlson, 2005). S'il existe un consensus sur l'existence de ces épisodes d'océans magmatiques, de nombreuses incertitudes demeurent sur leur durée de vie, l'évolution de la fraction de matériel fondu au cours du refroidissement, le degré de fractionnement chimique lors de leur solidification ainsi que sur les couplages avec le noyau métallique sous-jacent.

### 6.1.1 Un système en convection vigoureuse

La dynamique d'un océan magmatique avec une épaisseur supérieure à 1000 km est extrêmement turbulente (Solomatov, 2007). Afin de quantifier la vigueur de la convection dans un réservoir impliquant potentiellement toute l'épaisseur du manteau, il est possible de calculer le nombre de Rayleigh qui est le rapport entre le temps caractéristique de transfert de chaleur par conduction et le temps caractéristique de transfert de chaleur par advection (voir Eq.2) Pour entrer en convection, un reservoir de fluide doit avoir un Ra supérieur à une valeur critique $Ra_c$. Pour le manteau terrestre, $Ra_c \approx 1000$. En dessous de cette valeur critique, le refroidissement a lieu uniquement par conduction alors que pour des $Ra > Ra_c$, le transfert de chaleur se fait par advection et conduction : la convection. Le nombre de Rayleigh est
estimé dans le manteau terrestre actuel à $10^6 – 10^8$ ce qui indique qu'il doit convecter. En considérant que l'océan magmatique concerne tout le manteau et que la viscosité d'un magma liquide est comprise entre $10^{-2}$ et $10^2$ Pa.s, les nombres de Rayleigh obtenus sont compris entre $10^{20}$ et $10^{30}$ (Solomatov, 2007) ce qui indique que la convection est très vigoureuse dans cet environnement. Les tailles caractéristiques des écoulements sont ainsi d'autant plus petites avec de petits tourbillons infra-métriques. Afin de caractériser une telle dynamique, plusieurs approches peuvent être envisagées.

### 6.1.2 Modélisation d'un océan magmatique

Pour comprendre la dynamique de systèmes à très hauts nombres de Rayleigh, on peut utiliser des modèles analogiques en laboratoire impliquant un reservoir rempli de gaz (He ou SF₆ par exemple) chauffé par le bas et refroidi par le haut (Roche, 2020). Couplées à des modèles analytiques théoriques (Siggia, 1994), ces expériences permettent de contraindre le refroidissement d'un fluide très turbulent et notamment de proposer des relations entre le flux de chaleur observé en surface du réservoir (caractérisé par le nombre de Nusselt, Nu) et la vigueur de la convection (caractérisé par le nombre de Rayleigh, Ra). Différents régimes de convection peuvent ainsi être mis en évidence allant de turbulent (avec $Nu \propto Ra^{1/3}$) jusqu'à des régimes ultra-turbulents (avec $Nu \propto Ra^{2/7}$) dans lesquels se trouvent probablement les océans magmatiques (Solomatov, 2007). D'autres modèles analogiques permettent aussi de quantifier la stabilité d'une couche de cristaux dans un réservoir en convection laminaire et turbulente (Solomatov et al., 1993).

Une autre approche de modélisation est cette fois-ci numérique. Elle consiste à résoudre les équations de conservation de la quantité de mouvement, de la masse et de la chaleur par un code numérique. De nombreux modèles d'évolution uni-dimensionnels d'océan magmatiques ont été développés (Solomatov and Stevenson, 1993; Abe, 1997; Monteux et al., 2016a) et présentent les avantages suivants : (1) ils peuvent caractériser des échelles de temps courtes et leur résolution spatiale peut atteindre celle des longueurs caractéristiques d'un milieu en convection turbulente (2) les modèles 1D peuvent modéliser des contrastes de viscosité extrêmes tel que ceux existants entre fraction solide et liquide dans un océan magmatique (3) il est assez facile d'implémenter dans ces modèles des contraintes issues des données récentes de la pétrologie expérimentale tels que les courbes de fusion (solidus et liquidus) ou des paramètres dont la valeur pourrait changer en fonction de la profondeur comme la densité ou la viscosité. Cependant, si les modèles 1D calculent des profils de température moyens sur toute la profondeur du réservoir, ils ne peuvent pas caractériser d'hétérogénéités latérales de température ou de densité. Il est donc impossible pour ces modèles de suivre l'initiation de la convection solide qui peut se mettre en place avant même la fin de la cristallisation de l'océan magmatique ou encore de mettre en évidence des phases d'overturn engendrées par des gradients de densité ou de température ce qu'il est possible de faire avec des outils de modélisation numérique (2D ou 3D (Maurice et al., 2017)). De plus, les modèles 1D ne permettant pas de calculer le champ de vitesse, ils doivent donc faire des hypothèses fortes sur le calcul des épaisseurs des couches
limites qui se mettent en place dans les zones où la vitesse d'écoulement diminue fortement (au fond ou à la surface de l'océan magmatique ou à proximité des zones cristallisées et stables).

### 6.1.3 Profils adiabatiques

Dans un système tel que le manteau terrestre, la pression est hydrostatique et augmente avec la profondeur. Si ce milieu est à l'équilibre, la température augmente avec la pression et donc avec la profondeur en suivant un profil adiabatique (ou adiabat). Des hétérogénéités thermiques, liées par exemple au fait que le noyau soit très chaud, provoquent des instabilités dans ce système qui peuvent rendre le gradient de température supérieur au gradient adiabatique et se traduire par de la convection thermique. Quand cette convection est vigoureuse, le gradient de température est légèrement supérieur au gradient adiabatique. Afin d'estimer le profil de température dans un océan magmatique, il est donc possible de se baser sur le calcul du profil de température adiabatique en suivant la relation :

$$ \frac{dT}{dP} = \frac{\alpha T}{\rho C_p} \qquad (3) $$

avec $T$ la température, $P$ la pression, $\alpha$ le coefficient de dilatation thermique, $C_p$ la capacité thermique et $\rho$ la masse volumique dans le manteau. Ainsi, en faisant un hypothèse sur la température potentielle $T_p$ qui est la température de l'adiabat à la surface de la Terre à $P = 0$, il est possible de calculer au premier ordre des profils de température dans le manteau (Figure 21). Actuellement les données pétrologiques permettent d'estimer que le manteau profond terrestre suit un adiabat avec une température potentielle d'environ 1600 K. Or ce profil de température était probablement bien plus élevé juste après l'impact géant qui a donné naissance au système Terre-Lune. Si l'on fait l'hypothèse que tout le manteau était fondu, la température potentielle caractéristique du manteau post-impact était donc nécessairement supérieure à 3000 K (Figure 21).

### 6.1.4 Dynamique de refroidissement du manteau au stade "océan magmatique"

Suite à un impact géant, l'augmentation de température a probablement permis de faire fondre l'intégralité du manteau terrestre. Après l'impact, le noyau était potentiellement très chaud alors que la surface de la Terre était à une température d'équilibre pouvant dépasser les 3000 K selon les caractéristiques de l'atmosphère primitive (Lebrun et al., 2013). Rapidement après cet impact géant, l'océan magmatique s'est donc mis à convecter afin d'évacuer sa chauffer interne par la surface de la Terre où la température était beaucoup plus faible alors que le noyau était quant à lui toujours chaud. Après un refroidissement rapide à partir de la surface, la température décroit brutalement et passe en dessous de la température du liquidus au niveau du contact avec l'atmosphère où une croûte très fine et très instable se forme comme cela se produit à la surface d'un lac de lave. Deux couches limites thermiques extrêmement fines se sont
<figure><img src="image_271.png" /><figcaption>FIGURE 21 – Profils de température adiabatiques à l'intérieur du manteau pour différentes température potentielles en surface. Le solidus et liquidus (chondritique et péridotitique) sont représentés respectivement en vert et rouge (à partir de Monteux et al. (2016a)).</figcaption></figure>

mises en place (1) à la base du manteau où s'opère une contraste entre du matériau silicaté et du matériau métallique (2) à la surface de l'océan magmatique à la frontière entre l'atmosphère et le manteau. Ces couches limites sont associées à une discontinuité du champ de vitesse et donc à un transfert de chaleur uniquement par conduction. Le profil de température à l'intérieur de l'océan magmatique s'est ajusté rapidement sur un profil adiabatique (Figure 22). Le refroidissement se traduit par un abaissement progressif de ce profil interne adiabatique avec une température potentielle de plus en plus faible.

Le manteau poursuit son refroidissement global jusqu'à que le profil de température croise le liquidus (Figure 22). Dans un modèle de manteau chondritique, ce croisement se produit au niveau de la frontière manteau noyau et la solidification de l'océan magmatique va donc s'établir du bas du manteau en direction de la surface (Monteux et al., 2016a). Cependant, si le croisement se fait au milieu du manteau, deux fronts de solidification peuvent se mettre en
<figure><img src="image_272.png" /></figure>

FIGURE 22 – Évolution temporelle des profils de température à l'intérieur du manteau ter-
restre depuis une fusion globale (profil le plus à droite) jusqu'au profil actuel (le plus à
gauche). Le solidus et liquidus chondritiques sont illustrés respectivement en vert et rouge.
Le profil de température en bleu représente le profil à 40% de taux de fusion (i.e. 60% de
cristaux). Arrivé à ce profil, le manteau passe du stade océan magmatique au stade mush qui
est bien plus visqueux (à partir de Monteux et al. (2016a, 2020)).

place : un front qui part du milieu du manteau vers la surface et un front qui part du milieu en direction de la frontière manteau-noyau. Dans ce cas-là, un océan magmatique peut se retrouver isolé à la base du manteau solide (Labrosse et al., 2007).

Le refroidissement de l'océan magmatique va subir un ralentissement brutal lorsque son
profil de température va atteindre des valeurs correspondant à 40-50% de taux de fusion (i.e. un
profil de température compris entre liquidus et solidus correspondant à 50-60% de cristallisa-
tion) (Figure 22). Ce ralentissement pourrait intervenir moins d'un million d'année après l'im-
pact à l'origine du système Terre/Lune selon la viscosité du magma (Monteux et al., 2016a).
À ce moment, la viscosité du matériel augmente fortement car les cristaux commencent à
interagir mécaniquement les uns avec les autres et à empêcher la déformation du fluide. Le mélange magma/minéraux appelé "mush" a un comportement rhéologique proche du comportement d'un solide. La fraction de matériel solide augmente jusqu'à ce qu'un seuil séparant le régime turbulent du régime visqueux à l'état solide soit atteint (Abe, 1997). Aussitôt que ce seuil est atteint, l'efficacité du refroidissement du manteau primitif devient limitée. Ceci a des conséquences importantes sur l'état thermique du manteau profond mais aussi sur la capacité de refroidissement du noyau et ainsi sur l'histoire magnétique de la Terre.

Une deuxième phase de refroidissement, beaucoup plus lente, va ainsi se mettre en place une fois que tout le manteau se trouve à environ 40% de taux de fusion. Durant cette phase, la viscosité de la phase solide devient critique car c'est elle qui va gouverner le taux de refroidissement mais aussi les profondeurs auxquelles vont se maintenir les dernières poches de fusion (Monteux et al., 2020).

## 6.2 Des océans magmatiques à la dynamique mantellique actuelle

Au cours du stade océan magmatique, le manteau terrestre primitif a subi un refroidissement global rapide jusqu'à ce que sa fraction fondue diminue jusqu'à une valeur critique $\phi_{crit}$ (~40%) associée à une augmentation importante de sa viscosité (Solomatov, 2000; Monteux and Arkani-Hamed, 2016). Cette étape a été suivie d'une étape de refroidissement lent qui a provoqué une cristallisation quasi complète du manteau silicaté (Solomatov, 2000). Cette deuxième étape de refroidissement aurait pu durer des centaines de Ma, voire quelques Ga car sa dynamique a été gouvernée par la rhéologie du manteau solide à déformation lente, contrainrement à la première qui était gouvernée par la viscosité du magma (Solomatov, 2000; Ulvrová et al., 2012; Monteux et al., 2016a, 2020).

La profondeur à laquelle la cristallisation complète est atteinte est probablement régie par (1) la température de solidus qui contrôle la profondeur et la température de solidification, (2) l'épaisseur de la couche limite thermique supérieure où la chaleur convective du manteau profond est transférée à la surface par conduction, (3) le contraste de température des deux côtés de la couche limite thermiques qui contrôle l'efficacité de l'évacuation de la chaleur et (4) l'atmosphère primitive qui pourrait ralentir le refroidissement du manteau (Lebrun et al., 2013; Salvador et al., 2017). Dans le régime partiellement fondu, la viscosité du manteau solide est un paramètre clé, qui régit la dynamique du manteau et la formation des couches limites thermiques. Aujourd'hui, la viscosité du manteau terrestre est difficile à contraindre, en particulier dans le manteau inférieur (Cížková et al., 2012) et la valeur de la viscosité du manteau terrestre avant l'apparition d'événements majeurs de différenciation est encore plus hypothétique. Contraindre la valeur de cette grandeur physique est de première importance, car la viscosité du manteau solide profond régit probablement l'initiation de processus géologiques majeurs tels que la tectonique des plaques et la remontée de panaches mantelliques depuis la CMB (Sleep et al., 2014; Foley et al., 2014).
### 6.2.1 Taux de fusion et modèles de viscosité

La fraction de matériel fondu $\phi_{melt}$ joue un rôle fondamental dans la viscosité du manteau primitif ($\eta$) et donc sur sa dynamique de refroidissement. Au cours de son refroidissement, le manteau partiellement fondu cristallise et la fraction de matériel fondu diminue. Lorsque la fraction de matériel solide est suffisamment grande pour que les cristaux ne puissent plus entrer en suspension dans le fluide se produit une transition rhéologique majeure pour laquelle la viscosité augmente de plusieurs ordres de grandeur (Figure 23).

<figure><img src="image_273.png" /><figcaption>FIGURE 23 – Représentation schématique de l'évolution de la viscosité d'un manteau partiellement fondu en fonction de la fraction solide.</figcaption></figure>

Ainsi, selon le taux de fusion $\phi_{melt}$, plusieurs expressions pour le calcul de la viscosité du mélange magma/cristaux peuvent être considérées :

1. $\phi_{melt} = 1$ : Au dessus du liquidus, le matériel mantellique se comporte comme un liquide. Dans ce cas, la viscosité de la phase liquide $\eta_l$ dépend essentiellement de sa composition, de la température $T$ et de la pression $P$. Karki and Stixrud (2010) ont développé une approche calculatoire et expérimentale afin de contraindre la viscosité de la phase liquide $MgSiO_3$ à des pressions et températures compatibles avec celles du manteau terrestre. Ils obtiennent ainsi la relation suivante :
$$ \eta = \eta_t = \exp\left(-7.75 + 0.005P(\text{GPa}) - 0.00015P(\text{GPa})^2 + \frac{5000+135P(\text{GPa})+0.23P(\text{GPa})^2}{T-1000}\right) \quad (4) $$

Selon les estimations de Karki and Stixrude (2010), la viscosité de la phase silicatée liquide pourrait varier entre 0.001 et 1 Pa.s. Cette viscosité augmente en fonction de la pression et diminue avec la température.

2. $\phi_{melt} = 0$ : La viscosité de la fraction solide dans le manteau profond est un paramètre difficile à contraindre que ce soit pour le manteau profond actuel (Cížková et al., 2012) ou pour le manteau primitif même si la viscosité d'un manteau profond « à bridgmanite » est de plus en plus documentée (Boioli et al., 2017; Reali et al., 2019). Les deux mécanismes principaux de déformation associés à de la convection mantellique sont le fluage par diffusion et le fluage par dislocation. Dans le cas du fluage par diffusion, la déformation est liée à de la diffusion atomique de défauts cristallins et cette diffusion est activée thermiquement. Dans le cas du fluage par dislocation, la déformation est accommodée par des glissements le long de dislocations qui sont des lignes d'imperfections dans la matrice cristalline. Les modèles expérimentaux et théoriques pour calculer la viscosité s'accordent sur l'expression suivante et commune aux deux mécanismes :

$$ \eta = \eta_s = \frac{1}{A} \exp\left(\frac{E + PV}{RT}\right), \quad (5) $$

avec $P$ la pression, $R$ la constante des gaz (= 8.314 J K$^{-1}$ mol$^{-1}$) et $T$ la température absolue. $A$ est le pré-facteur de viscosité, $E$ est l'énergie d'activation, et $V$ est le volume d'activation pour le mécanisme de fluage approprié (diffusion ou dislocation). Ces 3 paramètres ($A$, $E$ et $V$) peuvent être contraints dans le manteau profond actuel par les vitesses de subduction de la lithosphère (Cížková et al., 2012). En raison de l'absence d'anisotropie sismique dans le manteau inférieur actuel, le fluage par diffusion est généralement considéré comme le mécanisme de déformation dominant dans le manteau profond (Boioli et al., 2017; Reali et al., 2019). Au début de l'histoire de la Terre, la température du manteau était plus élevée, et la diffusion ionique (donc le fluage de diffusion) devrait avoir été encore plus importante dans la fraction du manteau solide (Frost and Ashby, 1982) pour le manteau supérieur et inférieur. Par conséquent, nous pouvons supposer que la déformation de la phase solide interconnectée se produit uniquement par fluage par diffusion. De manière intéressante, la taille des grains intervient dans la valeur du pré-facteur de viscosité $A$ (Monteux et al., 2020).

3. $0 < \phi_{melt} < 1$ : Pour calculer la viscosité du manteau partiellement fondu, il convient de développer un formalisme tenant compte à la fois des propriétés de la phase solide et de la phase liquide. Nous détaillons ici le formalisme proposé par Lebrun et al. (2013). Lors des dernières étapes de solidification de l'océan magmatique, la fraction de matière solide augmente jusqu'à atteindre un seuil ($\phi_{melt} = \phi_{crit}$), qui sépare le régime "liquide" du régime "solide" (Solomatov, 2007). Cette fraction critique de matériel fondu peut être estimée par des expériences de migration de matériel fondu dans
des agrégats de matériel solide (Riley and Kohlstedt, 1991) ou par des experiences de compression tri-axiales à hautes températures (Scott and Kohlstedt, 2006). $\phi_{crit}$ dépend de nombreux paramètres (pression, température, composition chimique, taille des grains...). Scott and Kohlstedt (2006) ont estimé que $\phi_{crit}$ était compris entre 0.25 et 0.3 pour une péridotite à une pression de 300 MPa et une température d'environ 1500 K. Pour un océan magmatique profond, $\phi_{crit}$ pourrait atteindre 0.4 (Solomatov, 2000) (Figure 23).

Pour $\phi_{crit} < \phi_{melt} < 1$, la fraction en cristaux reste relativement faible et le mélange cristaux-magma se comporte comme un liquide. Ainsi la viscosité du matériau partiellement fondu est fonction de la viscosité du magma fondu $\eta_l$ (Roscoe, 1952) :

$$ \eta = \frac{\eta_l}{\left(1 - \left(\frac{1-\phi_{melt}}{1-\phi_{crit}}\right)\right)^{2.5}} \quad (6) $$

Pour $0 < \phi_{melt} < \phi_{crit}$, le mélange cristaux-magma a plutôt tendance à se comporter comme un solide du fait de la forte interaction entre les cristaux. Dans ce cas, la viscosité du matériel dépend cette fois-ci de la viscosité de la phase solide $\eta_s$ et peut s'exprimer de la façon suivante :

$$ \eta = \eta_s \exp(-26\phi_{melt}) \quad (7) $$

Au cours du refroidissement, dès que le seuil de fraction fondue ($\phi_{melt} = \phi_{crit}$) est atteint dans la majeure partie du manteau, l'efficacité de refroidissement du manteau primitif diminue considérablement, même si le manteau peut rester partiellement fondu à certaines profondeurs (Monteux et al., 2016a, 2020).

### 6.2.2 Dynamique de refroidissement du manteau au stade "mushy"

Nous avons vu précédemment que le refroidissement d'un manteau complètement fondu était un processus rapide qui ralentissait brutalement lorsque qu'une partie significative du manteau passait en dessous de 40-50% de taux de fusion. Ce point critique est atteint dans des échelles de temps comprises entre $10^3$ et $10^5$ ans (Abe, 1997; Monteux et al., 2016a). Lorsque ce taux de fusion est atteint dans l'essentiel du manteau, celui-ci passe du stade "océan magmatique" au stade "mush" (mélange de solide et de liquide) (Figure 24).

Malgré une augmentation importante de la viscosité à la suite d'une solidification significative, les nombres de Rayleigh atteints dans le manteau restent très élevés et la convection reste vigoureuse. Les hauts nombres de Rayleigh mis en jeu lors du refroidissement rendent difficile la caractérisation de la dynamique du manteau dans son ensemble en 2 ou 3 dimensions. Ainsi les modèles 1D peuvent encore une fois être pertinents pour ce genre d'approche (Abe, 1997; Lebrun et al., 2013; Monteux et al., 2016a). Les profils de température vont avoir tendance à suivre des profils adiabatiques. Comme illustré dans la Figure 24, l'essentiel du manteau
<figure><img src="image_274.png" /></figure>

FIGURE 24 – Évolution temporelle de la fraction de matériel dans le manteau fondu en fonction du rayon depuis le centre de la Terre au cours du stade océan magmatique (gauche) et du stade mush (à droite). À partir de Monteux et al. (2016a) et Monteux et al. (2020).

et notamment du manteau inférieur va se solidifier en moins de 90 Ma même si la zone en contact avec la frontière manteau noyau (i.e. la couche limite thermique basale) peut atteindre une température supérieure au solidus du manteau (~ 4000 K). Dans le manteau supérieur, les adiabats peuvent être localement supérieurs au solidus et des zones partiellement fondues peuvent subsister.

Les modèles récents montrent que la durée de vie et les profondeurs des dernières zones fondues dans le manteau supérieur sont fortement dépendantes des paramètres considérés pour calculer la viscosité du manteau supérieur (voir Eq.5). Selon le profil de viscosité dans le manteau, il est ainsi possible de maintenir deux couches partiellement fondues en profondeur : (1) une couche située à des profondeurs comprises entre 800 et 1200 km pendant 40 à 400 Ma et (2) une couche plus superficielle (30-50 km de profondeur) pendant plus de 500 millions d'années (Monteux et al., 2020). En associant ces échelles de temps obtenues par des modèles numériques à des observations géologiques, il est ainsi possible d'apporter des contraintes sur la viscosité du manteau primitif profond, sur les paramètres utilisés pour calculer cette viscosité (voir Eq.5) et sur certains paramètres tels que la taille des grains (Monteux et al., 2020).

## 6.3 Stabilité d'une couche de cristaux dans un océan magmatique

La flottabilité de la phase fondue est de fait gouvernée par la propension qu'auront les éléments plutôt lourds tels que le fer à intégrer la phase solide ou liquide (Funamori and Sato, 2010). Cette capacité peut être estimée en laboratoire dans des conditions de P/T caractéristiques du manteau profond en calculant le coefficient de partage ($D_i$) d'un élément $i$ entre
une phase solide silicatée (bridgmanite ou les polymorphes de l'olivine) et la phase liquide. Si les coefficients de partage sont aujourd'hui de mieux en mieux contraints pour une large gamme de pressions/températures, les implications géodynamiques et géochimiques de la différenciation du manteau primitif restent encore discutées. Ainsi, les travaux sur le partage du fer à hautes pressions peuvent conduire à des interprétations différentes : Nomura et al. (2011) suggèrent un liquide très riche en fer et plus dense que le manteau, alors que Andrault et al. (2012) suggèrent la remontée de la phase liquide vers la surface, tout en formant un résidu solide appauvri en éléments incompatibles en profondeur.

La différence de densité entre le solide et le liquide peut entrainer une séparation entre les 2 phases (Caracas et al., 2019) via des processus de sédimentation de cristaux dans l'océan magmatique. Cela peut produire une ségrégation chimique majeure entre les réservoirs géochimiques plus ou moins profonds que la différence de flottabilité soit positive ou négative. Les implications géodynamiques seront très différentes selon les cas. La vigueur de la convection s'opposera potentiellement à la sédimentation de la phase solide d'où la nécessité de modéliser la dynamique globale de cette ségrégation. L'ascension ou la descente de la phase liquide dans un manteau partiellement fondu pourra contribuer respectivement à du volcanisme en surface ou à la mise en place d'un océan magmatique en base de manteau.

Au début de la solidification de l'océan magmatique, les premiers cristaux formés ont pu se déposer en profondeur ou être efficacement mélangés par les mouvements de convection induits par le refroidissement de l'océan magmatique. Nous avons récemment développé des modèles numériques pour suivre l'évolution d'une fine couche de cristaux dans un océan magmatique (Monteux et al., 2023). Ces modèles ont caractérisé l'influence de la taille des cristaux, de la viscosité du magma et de la différence de densité entre les cristaux et le magma. Ces modèles ont mis en évidence l'importance du nombre de Shields (*Sh*) qui représente le rapport entre la contrainte tangentielle due à l'écoulement convectif et la contrainte de flottabilité liée à la différence de densité entre les cristaux et l'océan magmatique (Solomatov et al., 1993). Nous avons notamment contraint le nombre de Shields critique (*Sh*<sub>c</sub>) séparant le régime où les cristaux sédimentent au fond du réservoir et le régime où les cristaux restent en suspension. Nos résultats montrent qu'au cours de la solidification d'un océan magmatique profond, les cristaux sont plus facilement mis en suspension que déposés au fond de l'océan magmatique car les mouvements de convection entrainent efficacement les cristaux malgré leur poids (Figure 25). Cette découverte tend à renforcer l'idée que la cristallisation de l'océan magmatique s'est faite à l'équilibre plutôt que de manière fractionnée.

La solidification de l'océan magmatique a ainsi joué un rôle déterminant dans la composition chimique du manteau primitif et contraint ainsi quels types de minéraux ont pu se former en fonction de la profondeur dans le manteau primitif. Pour mieux contraindre ce phénomène, nous avons comparé la différence de densité entre l'océan magmatique et les cristaux solides pour différentes compositions chimiques et la différence de densité critique à partir de laquelle la ségrégation des cristaux est susceptible de se produire. Nos résultats montrent que la ségrégation de la bridgmanite, première phase minérale à se former dans l'océan magmatique, est plus susceptible de se produire au fond d'un océan magmatique si celui-ci est enrichi en SiO<sub>2</sub>.
<figure><img src="image_275.png" /></figure>

FIGURE 25 – Valeur du nombre de Shields (Sh) en fonction de la viscosité de la phase continue fondue ($\eta_c$) (ou Ra pour l'axe x supérieur). Les symboles carrés correspondent aux modèles utilisés dans notre étude. Les symboles noirs vides représentent les cas où le régime de suspension est atteint. Les symboles verts remplis représentent les cas où le régime de sédimentation est atteint. Les deux régimes sont séparés par des lignes pointillées (vertes pour $\eta_c \le 1$, bleues pour $\eta_c \ge 1$) illustrant le nombre de Shields critique (Sh$_c$). À partir de Monteux et al. (2023)

## 6.4 Perspectives : Couplage entre le manteau primitif terrestre et la Lune

Plusieurs dizaines de millions d'années après le début de son accrétion, la Terre a subi un impact géant qui a conduit à la formation du système Terre-Lune (Hartmann and Davis, 1975) (voir Section 2.4). La modélisation numérique de cet impact géant a montré que cet événement a entraîné une forte augmentation de la température et la formation d'un océan de magma profond (Nakajima and Stevenson, 2015). À la suite de cet impact catastrophique, un énorme transfert de chaleur et une ségrégation chimique se sont produits pendant la cristallisation de ce manteau primitif (Monteux and Arkani-Hamed, 2016). Cependant, les hétérogénéités isotopiques identifiées dans plusieurs éléments radiogéniques formés par des systèmes à courte durée de vie suggèrent que le manteau terrestre n'a jamais été mélangé au point d'être homogénéisé (Boyet and Carlson, 2005). La préservation de telles hétérogénéités est en contradiction
avec les événements de fusion et de mélange de grande ampleur attendus d'un impact géant (Nakajima and Stevenson, 2015) et avec la capacité de la convection mantellique à préserver des réservoirs non traités à grande échelle pendant plusieurs années Gyr. S'il existe un consensus sur l'importance des épisodes d'océan magmatique sur l'évolution thermochimique précoce de la Terre, des questions majeures demeurent sur l'évolution de la fraction fondue lors du refroidissement et sur les modalités de différenciation chimique au sein d'un manteau partiellement fondu, qui se sont probablement poursuivies des centaines de millions d'années après l'impact géant.

Plusieurs années après l'impact, la Lune a été entièrement accrétée par les débris de la collision (Kokubo and Ida, 2000). À cette époque, la Lune était beaucoup plus proche de la Terre qu'aujourd'hui (Chen and Nimmo, 2016), ce qui a généré d'énormes forces de marée sur les deux corps. La Lune a potentiellement commencé à se recristalliser à une distance allant de 5 à 15 rayons terrestres. Le forçage de marée induit par la présence de la Lune génère un cisaillement oscillatoire caractérisé par un taux de cisaillement maximal et une fréquence de marée. Un tel cisaillement oscillatoire peut influencer le flux convectif avec l'océan magmatique et le processus de cristallisation (Ogilvie and Lesur, 2012).

Le but de ce projet est de répondre aux questions suivantes :

1. Quelle est l'influence de la jeune Lune sur la dynamique de refroidissement et de solidification du manteau terrestre en fusion ?

2. Le chauffage hétérogène du manteau terrestre dû aux fortes interactions de marée avec la jeune Lune pourrait-il influencer la différenciation chimique et expliquer la signature géochimique que nous mesurons aujourd'hui ?

3. Comment la dissipation dans les intérieurs partiellement fondus de la Terre et de la Lune primitives a-t-elle influencé leur évolution spin-orbite à long terme ?

Pour répondre à ces questions fondamentales, nous proposons de rassembler des expertises complémentaires en physique minérale, géochimie, modélisation numérique en mécanique des fluides, mécanique des milieux continus et mécanique céleste, à l'interface entre les sciences de la Terre et l'astrophysique, afin de développer une approche numérique innovante et collaborative qui contraindra l'influence de la jeune Lune sur l'évolution thermochimique précoce de la Terre.

La dynamique de refroidissement d'un océan de magma profond est actuellement un sujet très débattu (Ballmer et al., 2017; Maurice et al., 2017; Laneuville et al., 2018; Agrusta et al., 2020). Cependant, aucune des études récentes n'a considéré le couplage mécanique entre la ségrégation des cristaux, les effets de marées et la convection thermique turbulente en 2D/3D tout en considérant des contraintes expérimentales réalistes et récentes relatives aux courbes de fusion, au fractionnement chimique, aux diffusivités thermiques et aux paramètres thermo-dynamiques. Nous proposons pour la première fois de résoudre de manière cohérente le chauffage par les marées, l'évolution de la solidification et la différenciation chimique dans le manteau hadéen depuis le stade de cristallisation de l'océan magmatique jusqu'à la transition archéenne. Nous combinerons les principaux résultats des deux tâches suivantes pour
<figure><img src="image_276.png" /></figure>

FIGURE 26 – Principe général de l'approche numérique prévue pour la modélisation thermo-orbitale couplée du système Terre-Lune. La dissipation de marée pour un intérieur partiellement fondu sera calculée à l'aide du code LPG-Tide (TÂCHE 1) à partir des profils thermiques, de densité et de fusion moyens prédits par les modèles de transfert de chaleur et de matière fondue (TÂCHE 2). La puissance dissipée calculée affectera les bilans thermique et orbital. Le changement d'orbite et de rotation du système Terre-Lune sera calculé à l'aide des modèles d'évolution orbitale développés à l'IMCCE à partir de la puissance dissipée totale de la Terre et de la Lune, qui à son tour modifiera le potentiel de marée.

proposer les premiers modèles d'évolution couplés complets du système Terre-Lune pendant l'ère Hadéenne.

Le projet DEMOTIC (Figure 26) vise à développer de nouveaux modèles numériques pour suivre l'évolution de la température et de la solidification de la Terre. En particulier l'objectif est de contraindre : (1) les effets réciproques des marées liés à la distance Lune-Terre et (2) l'influence des marées sur l'évolution thermochimique de la Terre hadéenne, y compris le fractionnement des éléments chimiques pendant la cristallisation de la MO. Pour la première fois, nous allons résoudre de manière cohérente le chauffage par les marées, l'évolution de la matière en fusion et la différenciation chimique dans le manteau Hadéen depuis le stade de cristallisation de la MO jusqu'à la transition archéenne. Ainsi, nous comprendrons si les marées et le chauffage hétérogène du manteau terrestre dû aux fortes interactions de marée avec la jeune Lune (1) peuvent influencer la différenciation chimique et expliquer la signature géochimique que nous mesurons aujourd'hui à partir des hétérogénéités isotopiques observées dans
la Terre profonde et (2) peuvent avoir une influence sur l'évolution orbitale e à long terme du système Terre-Lune. Dans ce but, le projet DEMOTIC rassemblera l'expertise multidisciplinaire complémentaire du LMV, du LPG Nantes et de l'IMCCE, à l'interface entre les sciences de la Terre et l'astrophysique, afin de développer une approche numérique innovante, collaborative et unique. Ce projet dont je suis le porteur principal vient d'obtenir un financement INSU (Origines) et fera l'objet d'une demande de financement ANR PRC en 2024.

Références personnelles

Monteux, J., Andrault D., Samuel H., 2016, On the cooling of a deep terrestrial magma ocean. Earth and Planetary Science Letters vol.448, p.140-149. LIEN HAL

Monteux, J., Andrault D., Guitreau M., Samuel H., Demouchy S., 2020, A mushy Earth's mantle for more than 500 Myr after the magma ocean solidification, Geophysical Journal International vol.221, p.1165-1181, 2. LIEN HAL

Monteux, J., Qaddah B., Andrault D., 2023, Conditions for segregation of a crystal-rich layer within a convective magma ocean, Journal of Geophysical Research, Volume 128, Issue 5 LIEN HAL
<figure><img src="image_277.png" /></figure>

<figure><img src="image_278.png" /></figure>

# On the cooling of a deep terrestrial magma ocean

J. Monteux<sup>a,*</sup>, D. Andrault<sup>a</sup>, H. Samuel<sup>b</sup>

<figure><img src="image_279.png" /></figure>

CrossMark

## ARTICLE INFO

**Article history:**

Received 18 November 2015

Received in revised form 21 April 2016

Accepted 7 May 2016

Available online 24 May 2016

Editor: C. Sotin

**Keywords:**

early Earth

thermal evolution

magma ocean

numerical modeling

## ABSTRACT

Several episodes of complete melting have probably occurred during the first stages of the Earth's evolution. We have developed a numerical model to monitor the thermal and melt fraction evolutions of a cooling and crystallizing magma ocean from an initially fully molten mantle. For this purpose, we numerically solve the heat equation in 1D spherical geometry, accounting for turbulent heat transfer, and integrating recent and strong experimental constraints from mineral physics. We have explored different initial magma ocean viscosities, compositions, thermal boundary layer thicknesses and initial core temperatures.

We show that the cooling of a thick terrestrial magma ocean is a fast process, with the entire mantle becoming significantly more viscous within 20 kyr. Due to the slope difference between the adiabats and the melting curves, the solidification of the molten mantle occurs from the bottom up. In the meantime, a crust forms due to the high surface radiative heat flow, the last drop of fully molten silicate is restricted to the upper mantle. Among the studied parameters, the magma ocean lifetime is primarily governed by its viscosity. Depending on the thermal boundary layer thickness at the core-mantle boundary, the thermal coupling between the core and magma ocean can either insulate the core during the magma ocean solidification and favor a hot core or drain the heat out of the core simultaneously with the cooling of the magma ocean. Reasonable thickness for the thermal boundary layer, however, suggests rapid core cooling until the core-mantle boundary temperature results in a sluggish lowermost mantle. Once the crystallization of the lowermost mantle becomes significant, the efficiency of the core heat loss decreases. Since a hotter liquidus favors crystallization at hotter temperatures, a hotter deep mantle liquidus favors heat retention within the core. In the context of an initially fully molten mantle, it is difficult to envision the formation of a basal magma ocean or to prevent a major heat depletion of the core. As a consequence, an Earth's geodynamo sustained only by core cooling during 4 Gyr seems unlikely and other sources of motion need to be invoked.

© 2016 Elsevier B.V. All rights reserved.

## 1. Introduction

Geochemical evidence (Touboul et al., 2012; Rizo et al., 2013) suggests that the Earth's mantle has experienced several episodes of global melting during its early evolution, leading to the formation of the early continental crust and facilitating the core formation (Kleine et al., 2009). These episodes were probably enhanced by giant impacts occurring during the late stages of planetary formation (Agnor et al., 1999). Although not yet clearly established, it is likely that these giant impacts, such as the one that is thought to have formed the Earth-Moon system, could have melted 30 to 100% of the Earth's mantle depending on the impactor/target mass ratio and on the pre-impact ther-

mal state of the target (Canup, 2012; Ćuk and Stewart, 2012; Nakajima and Stevenson, 2015). During the cooling and the subsequent crystallization of a magma ocean (MO), compatible elements (e.g. Mg, Cr) were preferentially collected in the solid phase while the incompatible elements (e.g. Al, Na, Fe) selectively partitioned into melts. In addition to temperature, the degree of solid-melt fractionation is highly sensitive to a variety of physical parameters, including pressure (Nomura et al., 2011; Andrault et al., 2012). Hence, characterizing the cooling of a deep terrestrial magma ocean and in particular the timescale and depth at which the last drop of melt solidifies are of first importance to understand the current chemical composition of the Earth's mantle and the dating of its major differentiation events (Boyet and Carlson, 2005).

The composition and the rheology of such a magma ocean directly affect its lifetime, but remain poorly constrained (Solomatov, 2007). The magma ocean is composed of low viscosity molten sili-

<footnote><sup>a</sup> Laboratoire Magmas et Volcans, Université Blaise Pascal, CNRS, IRD, Clermont-Ferrand, France</footnote>

<footnote><sup>b</sup> Institut de Recherche en Astrophysique et Planétologie, CNRS, Toulouse, France</footnote>

<footnote>* Corresponding author.
E-mail address: j.monteux@opgc.univ-bpclermont.fr (J. Monteux).</footnote>
cate material but its chemical composition remains uncertain, with a MgO/SiO₂ ratio around those of chondritic or peridotitic compositions (Ringwood, 1966; Allègre et al., 1995; Javoy et al., 2010). Recent high-pressure laboratory measurements report the solidus and liquidus of both a chondritic and peridotitic mantle compositions up to pressures that are compatible with the Earth's lowermost mantle conditions (Fiquet et al., 2010; Andrault et al., 2011). Moreover, recent shock experiments now provide important constraints on the thermodynamic parameters used to determine the adiabatic profiles in the magma ocean up to 140 GPa (Mosenfelder et al., 2009; Thomas et al., 2012; Thomas and Asimow, 2013). Since the difference between their slopes governs the depth at which crystallization is initiated, both the liquidus and the adiabats play a key role in the cooling of the magma ocean. If the adiabat had a steeper slope than the liquidus in the mid-mantle (Mosenfelder et al., 2007; Stixrude et al., 2009), solidification would start at mid-mantle depth. In this case, a lowermost magma ocean would cool and solidify much more slowly because of the thermal blanketing of the overlying solid mantle (Labrosse et al., 2007). However, if the mantle liquidus had a steeper slope than the adiabat through the whole mantle (Thomas et al., 2012), solidification would start from the CMB thus reducing the likeliness of a basal magma ocean, unless invoking an enrichment in dense incompatible elements in the residual liquid. In any case, the important dynamical change does not occur when the adiabat crosses the liquidus, because the mantle keeps its liquid behavior, but rather when the degree of partial melting decreases below a critical value from which the mantle behaves as a solid. Therefore, the recent determination of melting curves and elastic parameters of silicate melts up to core-mantle boundary (CMB) conditions offers a great opportunity to improve our knowledge of the cooling dynamics of a deep terrestrial magma ocean.

The magma oceans such as the one generated by the Moon-forming impact participated to the core-formation process. The early thermal state of the core remains poorly constrained. It results from the contribution of the accretionary processes (Safronov, 1978; Kaula, 1979), including giant impact (Tonks and Melosh, 1992) and radiogenic heating (Yoshino et al., 2003) as well as the conversion of potential energy into heat via viscous dissipation during the metal/silicate separation (Ke and Solomatov, 2009; Monteux et al., 2009; Ricard et al., 2009; Samuel et al., 2010). The combined processes leading to core formation can yield a wide range of possible early thermal states, depending on the nature and timescale of core formation processes. The core could initially have had a temperature close to the deep mantle temperature if thermal equilibration was efficient. Alternatively, it could have been hotter than the mantle if the gravitational potential energy released during core formation was largely retained within the core itself, a situation which would be followed by a strong heating of the lowermost mantle from this superheated core (Samuel et al., 2010). In turn, the thermo-mechanical properties of the magma ocean can have a strong influence on the early evolution of the heat repartition between the core and the mantle. A key question is to determine how much a deep magma ocean can enhance core cooling. This can have important consequences on the duration and the generation of the Earth's dynamo (Monteux et al., 2011).

The low magma ocean viscosities resulting from the hot early temperatures imply that the cooling of such a deep molten mantle was highly turbulent (Solomatov, 2007). Studies of the early mantle have either characterized the cooling of a magma ocean restricted to the first 1000 km (Abe, 1997) or did not consider the presence of a molten layer just above the core-mantle boundary, and its effect of the CMB heat flow (Nakagawa and Tackley, 2014). However, the hypothesis of an early largely molten mantle combined with the determination of solidus/liquidus and thermodynamical properties of silicate melts up to 140 GPa now allow a more accurate

characterization of the cooling of a deep terrestrial magma ocean and the thermal coupling with its underlying core. The aim of this work is to constrain the lifetime of a deep magma ocean and to determine the pressure at which the magma ocean crystallization finished. To achieve these goals, we have developed a numerical model to characterize the early evolution of (i) the temperature and melt fraction of an initially fully molten isochemical mantle and (ii) the temperature of the core. We incorporate in our models recent and strong experimental constraints on the solidus and liquidus profiles and on the thermodynamical properties of silicate melts up to ~140 GPa. We explore different core temperatures, magma ocean compositions and viscosities.

## 2. Convective cooling of the magma ocean

Miller et al. (1991) characterized the cooling and the subsequent crystallization of a magma ocean with a chondritic composition as a sequence of isentropes with decreasing potential temperature. Later on, Abe (1997) investigated the thermal evolution of magma ocean using a one-dimensional heat transfer model. However, these studies were restricted to the first 1000 km and did not integrate the mutual influence of the magma ocean and its underlying material on the cooling. Labrosse et al. (2007) studied the cooling of a stable dense molten layer above the CMB overlaid by a solid mantle. In their model they consider the crystallization of a single-component (forsterite) magma ocean assuming a solidification proceeding from the top to the bottom according to Mosenfelder et al. (2007). More recently, Nakagawa and Tackley (2014) characterized the coupled thermal evolution of Earth's early mantle and core considering a 2900 km thick viscous mantle but ignoring the potential presence of a molten layer just above the core-mantle boundary, and its effect of the CMB heat flow. Here, we model the secular cooling of an initially fully molten magma ocean by convective transport of heat in a 1-D spherically symmetric geometry. We assume a multicomponent chemically homogeneous magma ocean made of a combination of forsterite, enstatite, fayalite, anorthite and diopside. In the following sections, we describe the model setup and equations.

### 2.1. Physical model for planetary thermal evolution

We model the thermal evolution of a 2900 km thick isochemical silicate mantle overlying an iron core by solving the conservation of energy in a one-dimensional, spherically symmetric domain (with a radius ranging from 3500 to 6400 km):

$$ \rho C_p \frac{\partial T}{\partial t} = \nabla \cdot (k \nabla T), \quad (1) $$

with $\rho$ the density, $C_p$ the mantle heat capacity, $T$ the temperature, $t$ the time and $k$ the thermal conductivity. Among the heat sources that have potentially delivered the energy required for significant melting in the early Earth, the decay of short-lived radioactive isotopes such as $^{26}$Al and $^{60}$Fe have probably played a major role especially for 10 to 100 km size objects (Yoshino et al., 2003). However, their half-life times (0.73 My and 1.5 My respectively) (Carlson and Lugmair, 2000) are much shorter than the time at which the Moon forming impact is supposed to have occurred (between 30 and 100 Myrs after the formation of the first solids of the Solar System) (Kleine and Rudge, 2011). Concerning the long-lived radioactive elements such as $^{40}$K, Th or U, their concentrations were certainly significant at the time of the Moon-forming impact, but their heat production rates are much smaller. Hence the contribution form the long-lived radio-active elements during the magma ocean lifetime is negligible. Thus, we can reasonably neglect radiogenic heating in our models.
Thermal energy is transferred by convection in the region where the temperature gradient is steeper than the adiabatic temperature gradient, or by conduction elsewhere. To account for the heat transfer within a vigorously convecting magma ocean, in Eq. (1), the thermal conductivity $k$ is the sum of the intrinsic thermal conductivity $k_c$ and an effective conductivity due to thermal convection $k_v$. Following Neumann et al. (2014), the latter is estimated as follows:

$$k_v = F_{conv}L/\Delta T, \qquad (2)$$

where $L$ is the thickness of the magma ocean at time $t$, $F_{conv}$ is the convective heat flux at radius $r$ and time $t$ and $\Delta T$ is the difference between the temperature profile and the adiabatic profile $T_{ad}$ with a potential temperature corresponding to the temperature of the surface of the magma ocean.

At radius $r$, the convective heat flux $F_{conv}$ depends on the local Rayleigh number $Ra$:

$$Ra = \frac{\alpha g C_p \rho^2 \Delta T L^3}{k_c \eta}, \qquad (3)$$

where $\alpha$ is the thermal expansion coefficient of the magma ocean, $g$ is the gravitational acceleration assumed to be constant through the whole mantle and $\eta$ is the local dynamic viscosity.

In the context of an initially fully molten mantle, the convective velocities are estimated to $u_0 \sim 10$ m/s (Solomatov, 2007) leading to Reynolds number values ($Re = \rho u_0 L/\eta$) in the order of $10^9$ (with $\eta = 100$ Pa s and $\rho = 4000$ kg m$^{-3}$) and, hence, to highly turbulent convective cooling. The low magma ocean viscosities induce Rayleigh numbers as large as $10^{30}$ (Solomatov, 2007). In such a context, two regimes arise depending on $Ra$ (Solomatov, 2007; Neumann et al., 2014): the “soft-turbulent” regime (if $Ra < 10^{19}$) where

$$F_{conv} = 0.089 \frac{k_c \Delta T}{L} Ra^{1/3}, \qquad (4)$$

and the “hard-turbulent” regime (if $Ra \ge 10^{19}$) where

$$F_{conv} = 0.22 \frac{k_c \Delta T}{L} Ra^{2/7} Pr^{-1/7} \lambda^{-3/7}, \qquad (5)$$

where $Pr$ is the local Prandtl number ($= C_p \eta / k_c$) and $\lambda$ is the aspect ratio for the mean flow. We assume that $\lambda = 1$ and that rotation does not have any significant effect on the heat flow (Solomatov, 2007).

Depending on their ability to migrate either towards the solid phase (compatible) or towards the liquid phase (incompatible), the relative abundance of some chemical elements can potentially modify the buoyancy of both the liquid and solid phases during the crystallization of a magma ocean. This chemical fractionation process could be enhanced by the fractional crystallization of the magma ocean at least in the upper mantle (Solomatov, 2007) and would lead to a liquid residual phase that is more or less buoyant than the solid phase. A dense liquid phase could favor a basal magma ocean (Labrosse et al., 2007) or large scale cumulate overturns (Elkins-Tanton et al., 2003, 2005) while a dense solid phase would enhance volcanic activity (Moyen and Martin, 2012; Martin et al., 2014). However, the values of the partition coefficients at high pressures of the elements that have a key contribution in the density of each phase (e.g. Fe) are still debated (Andrault et al., 2012). In addition, crystals may also gravitationally separate with the dense cumulates sinking and the lighter ones floating toward the surface leading to a supplementary segregation process. In a highly turbulent environment such as a thick magma ocean, the vigor of the convection probably prevents from any chemical segregation especially during the early stages of the magma ocean solidification. Indeed, Tonks and Melosh (1990) have

shown that the Rayleigh number for a planetary scale magma ocean was so high that crystals would remain entrained in the magma which would effectively preclude crystal separation. Therefore, we do not consider any spatial variation in the chemical composition and we currently consider neither the mass flux of melt owing to gravitational separation nor the mass flux of melt owing to convective mixing.

Viscosity is strongly dependent on the local melt fraction $\phi$, which is calculated as follows:

$$\phi = \frac{T - T_{sol}}{T_{liq} - T_{sol}}, \qquad (6)$$

where $T_{liq}$ and $T_{sol}$ are the liquidus and solidus temperatures, respectively. Following Abe (1997) the viscosity of partially molten silicates $\eta$ varies between a fully molten end-member $\eta_m$ (assumed to be constant) and a solid end-member $\eta_s$ that is temperature dependent:

$$\eta = \text{MIN} \left[ \frac{(1-\phi)\rho_m\eta_s + \phi\rho_s\eta_m}{(1-\phi)\rho_m + \phi\rho_s}, \right. \\ \left. \eta_m \left\{ \frac{(1-\phi)\rho_m + \phi\rho_s}{(1-A)(1-\phi)\rho_m + \phi\rho_s} \right\}^{2.5}, 10^{21} \right], \qquad (7)$$

where $\rho_m$ is the density of the molten material, $\rho_s$ is the density of the solid material and $A = 1.67$ (Abe, 1997). When the temperature, $T$, is lower than the solidus temperature of mantle, the viscosity of the solid mantle $\eta_s$ is estimated as follows (Abe, 1997):

$$\eta_s = \eta_{s,0} \exp \left( B \frac{T_{liq}}{T} \right). \qquad (8)$$

We used $\eta_{s,0} = 256$ Pa s, and $B = 25.17$ based on the olivine rheology (Karato and Wu, 1993; Abe, 1997). In the above equation, as for an Arrhenius relation, the viscosity of the solid mantle increases with the pressure. Assuming an adiabatic temperature profile with a potential temperature of 1600 K (Tackley, 2012) leads to a viscosity value of $\sim 10^{23}$ Pa s in the lowermost mantle compatible with estimates of the present-day mantle viscosity profiles (Čižková et al., 2012). Considering a relationship that involves the solidus rather than the liquidus in Eq. (8) would not affect our results since both the liquidus and the solidus used in our models have a similar trend. Since we consider here the evolution of a fully molten to partially molten magma ocean, the cooling timescale of the magma ocean is mainly governed by $\eta_m$ and is weakly dependent on $\eta_s$ as we shall see later.

A strong increase in viscosity occurs when the melt fraction equals 40% (Abe, 1997). Hence, when the melt fraction approaches this critical value, the presence of crystals drastically reduces the efficiency of the magma ocean cooling. In our models, the mantle is considered as part of the magma ocean as long as its melt fraction is larger than 50% (Neumann et al., 2014) and we stop our simulations when the thickness of the magma ocean drop below 100 km.

The pressure profile $P$ is obtained fitting the PREM model (Dziewonski and Anderson, 1981) with a quadratic function of the radius $r$ and is assumed to remain constant with time:

$$P = 4.0074 \times 10^{11} - 91862r + 0.0045483r^2. \qquad (9)$$

## 2.2. Melting curves and adiabats

Due to the uncertainties related to the chemical composition of the magma ocean, we consider in our study two models for melting curves and elastic parameters: the F-peridotic model and the A-chondritic model. We describe these two models in the following sections.
### 2.2.1. Melting curves

The solidus and liquidus play a major role in the early thermal evolution of the magma ocean. Recent laboratory experiments now constrain the liquidus and solidus of mantle-like material up to pressures compatible with the CMB conditions (Fiquet et al., 2010; Andrault et al., 2011). However, the differences between these two studies, in particular the liquidus temperatures, appear too large to be solely due to the difference in chemical compositions between the two types of mantle materials. Regardless of the controversy, we performed calculations using the melting curves from both reports, leading to a F-peridotitic model (Fiquet et al., 2010) and a A-chondritic model (Andrault et al., 2011). The solidus and liquidus profiles are obtained fitting experimental results with a modified Simon and Glatzel equation (Simon and Glatzel, 1929). For pressures below $P = 20$ GPa, we use experimentally determined solidus and liquidus temperatures of chondritic mantle from Herzberg and Zhang (1996):

$$T_{sol} = 1661.2 \left( \frac{P}{1.336 \times 10^9} + 1 \right)^{(1/7.437)}, \quad (10)$$

$$T_{liq} = 1982.1 \left( \frac{P}{6.594 \times 10^9} + 1 \right)^{(1/5.374)}, \quad (11)$$

with $T_{liq}$ is the liquidus temperature and $T_{sol}$ the solidus temperature.

Since the difference between the F-peridotitic and A-chondritic solidus for pressures larger than $P = 20$ GPa is not large, we use the experimentally determined A-chondritic solidus from Andrault et al. (2011):

$$T_{sol} = 2081.8 \left( \frac{P}{101.69 \times 10^9} + 1 \right)^{(1/1.226)} \quad (12)$$

For pressures above 20 GPa, the difference between the F-peridotitic and the A-chondritic liquidus is more important. We use the following expression

$$T_{liq} = c_1 \left( \frac{P}{c_2} + 1 \right)^{(1/c_3)}, \quad (13)$$

with ($c_1 = 78.74$, $c_2 = 4.054 \times 10^6$, $c_3 = 2.44$) for F-peridotitic liquidus (Fiquet et al., 2010) and ($c_1 = 2006.8$, $c_2 = 34.65 \times 10^9$, $c_3 = 1.844$) for A-chondritic liquidus (Andrault et al., 2011).

### 2.2.2. Thermodynamical parameters

The thermodynamical parameters for the molten magma ocean are closely related to its chemical composition. Volume and elastic parameters of silicate liquids has been recently characterized up to 140 GPa using shock compression experiments (Mosenfelder et al., 2007, 2009; Thomas et al., 2012; Thomas and Asimow, 2013). We assume here two multicomponent systems for (i) a A-chondritic composition (62% enstatite + 24% forsterite + 8% fayalite + 4% anorthite + 2% diopside) and (ii) a F-peridotitic composition (33% enstatite + 56% forsterite + 7% fayalite + 3% anorthite + 0.7% diopside). Using fourth-order Birch–Murnaghan/Mie–Grüneisen equation of state fits for molten silicate liquids from Thomas and Asimow (2013), we obtain the melt density $\rho_m$, the volumetric thermal expansion $\alpha$ as a function of pressure as well as the specific heat $C_p$ of the molten material for these two multicomponent assemblages. The density of the solid phase is then calculated as:

$$\rho_s = \rho_m + \Delta\rho, \quad (14)$$

with $\Delta\rho$ the density difference between solid and liquid (see Tables 1 and 2 for values).

<table><caption>**Table 1**

Constant and fixed parameter values for numerical models.</caption><tbody><tr><td>Earth radius</td><td>R</td><td>6370 km</td></tr><tr><td>Core radius</td><td>R<sub>core</sub></td><td>3470 km</td></tr><tr><td>Mantle thickness</td><td>R - R<sub>core</sub></td><td>2900 km</td></tr><tr><td colspan="3"><strong>Mantle properties</strong></td></tr><tr><td>Solid density</td><td>ρ<sub>s</sub></td><td>= ρ<sub>m</sub> + Δρ</td></tr><tr><td>Density contrast</td><td>Δρ/ρ</td><td>1.5% (Tosi et al., 2013)</td></tr><tr><td>Specific enthalpy change</td><td>ΔH</td><td>4 × 10<sup>5</sup> J/kg (Ghosh and McSween, 1998)</td></tr><tr><td>Viscosity of melt phase</td><td>η<sub>m</sub></td><td>1–10<sup>4</sup> Pa s</td></tr><tr><td>Bottom TBL thickness</td><td>e<sub>TBL</sub></td><td>10<sup>-3</sup>–10<sup>-3</sup> m</td></tr><tr><td colspan="3"><strong>Core properties</strong></td></tr><tr><td>Density</td><td>ρ<sub>Fe</sub></td><td>10 000 kg m<sup>-3</sup></td></tr><tr><td>Heat capacity</td><td>C<sub>p,Fe</sub></td><td>800 J kg<sup>-1</sup> K<sup>-1</sup></td></tr></tbody></table>

For a zone of partial melting, the density $\rho'$, the coefficient of volumetric thermal expansion $\alpha'$ and the specific heat $C'_p$ are given as follows (Solomatov, 2007):

$$\frac{1}{\rho'} = \frac{1 - \phi}{\rho_s} + \frac{\phi}{\rho_m}, \quad (15)$$

$$\alpha' = \alpha + \frac{\Delta\rho}{\rho(T_{liq} - T_{sol})}, \quad (16)$$

$$C'_{p} = C_{p} + \frac{\Delta H}{T_{liq} - T_{sol}}, \quad (17)$$

where $\Delta H$ is the latent heat released during solidification.

### 2.2.3. Adiabats

In vigorously convecting systems such as magma oceans, the temperature distribution is nearly adiabatic and isentropic (Solomatov, 2007). In one-phase systems, such as a completely molten or a completely solid layer, the equation for an adiabat is

$$\left(\frac{\partial T}{\partial r}\right)_s = -\frac{\alpha g}{C_p} T. \quad (18)$$

In two-phase systems, the effects of phase changes need to be considered (Solomatov, 2007). The equation for such adiabat is given by:

$$\left(\frac{\partial T}{\partial r}\right)_s = -\frac{\alpha' g}{C'_p} T. \quad (19)$$

This leads to a two-phase adiabat that is steeper than the purely liquid or solid one-phase adiabats (Solomatov, 2007). The adiabatic temperature profiles are calculated by numerical integration of Eq. (18) and Eq. (19) using a fourth-order Runge–Kutta method (Press et al., 1993). These adiabatic temperature profiles are used to calculate at each depth and, when super-adiabatic, the temperature difference $\Delta T$ from Eq. (3). The liquidus and solidus profiles as well as the adiabatic profiles obtained from Eq. (18) for temperatures ranging between 1400 and 4000 K are shown in Fig. 1. We start our models assuming an adiabatic temperature profile with a potential surface temperature $T_p = 3200$ K.

### 2.3. Boundary conditions

Large impacts can generate a rock vapor atmosphere that can last for some years until its energy is radiated to space (Svetsov, 2005). The presence of an atmosphere is expected to slow down the radiation of heat to space (Hamano et al., 2013; Lebrun et al., 2013). However, most of the pre-impact atmosphere is likely to be eroded after a giant impact (Shuvalov, 2009). Thus, in our models,
<table><caption>**Table 2**
Variable and non-dimensional parameter values for numerical models.</caption><tr><td>Melt density</td><td>ρ<sub>m</sub></td><td>A-model: 2684–5274 kg m<sup>-3</sup><br/>F-model: 2679–5378 kg m<sup>-3</sup></td><td>Computed from Thomas and Asimow (2013)</td></tr><tr><td>Heat capacity</td><td>C<sub>p</sub></td><td>A-model: 1742 J kg<sup>-1</sup> K<sup>-1</sup><br/>F-model: 1800 J kg<sup>-1</sup> K<sup>-1</sup></td><td>Computed from Thomas and Asimow (2013)</td></tr><tr><td>Thermal expansion coefficient</td><td>α</td><td>A-model: 1.3 × 10<sup>-5</sup>–7.9 × 10<sup>-5</sup> K<sup>-1</sup><br/>F-model: 2 × 10<sup>-5</sup>–9.6 × 10<sup>-5</sup> K<sup>-1</sup></td><td>Computed from Thomas and Asimow (2013)</td></tr><tr><td>Viscosity of solid phase</td><td>η<sub>s</sub></td><td>1–10<sup>21</sup> Pa s</td><td>From Eq. (8) with η<sub>s,0</sub> = 256 Pa s and B = 25.17</td></tr><tr><td>Viscosity of the magma ocean</td><td>η</td><td>5–10<sup>7</sup> W m<sup>-1</sup> K<sup>-1</sup></td><td>From Eq. (7)<br/>= k<sub>c</sub> + k<sub>v</sub></td></tr><tr><td>Total conductivity</td><td>k</td><td></td><td></td></tr><tr><td>Rayleigh number</td><td>Ra</td><td>at t = 0: 1 × 10<sup>27</sup>–3 × 10<sup>27</sup></td><td>Computed from Eq. (3)</td></tr><tr><td>Prandtl number</td><td>Pr</td><td>350–3.6 × 10<sup>24</sup></td><td>= C<sub>p</sub>η/k<sub>c</sub></td></tr><tr><td>Reynolds number</td><td>Re</td><td>at t = 0: Re ~ 10<sup>9</sup></td><td>From Solomatov (2007)</td></tr></table>

<figure><img src="image_280.png" /><figcaption>**Fig. 1.** Adiabats (with $T_p$ ranging between 1400 and 4000 K) computed for the A-chondritic model (left) and the F-peridotitic model (right). The corresponding solidus and liquidus are represented in green and red respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

we neglect the effect of thermal blanketing, and impose a radiative heat flux boundary condition at the surface:

$$F_{surf} = \sigma (T_{surf}^{4} - T_{eq}^{4}), \qquad (20)$$

with $T_{surf}$ the temperature at the surface of the MO, $\sigma$ the Stefan-Boltzman constant and $T_{eq} = 273$ K the expected equilibrium surface temperature.

At the base of the silicate mantle, core-mantle thermal coupling is accounted for via a conductive heat flux imposed at the core-mantle boundary (CMB):

$$F_{core} = \frac{k_c(T^{core} - T^{mantle}_{CMB})}{e_{TBL}}, \qquad (21)$$

where $T^{core}$ is the average core temperature at the CMB (i.e. we neglect the thermal boundary layer within the core) and $T^{mantle}_{CMB}$ is the mantle temperature right above the CMB. $e_{TBL}$ is the thickness of the thermal boundary layer at the bottom of the mantle where the heat is extracted from the core by conduction. $T^{mantle}_{CMB}$ is obtained solving Eq. (1) while $T^{core}$ is obtained from:

$$V_{core} \rho_{Fe} C_{p,Fe} \frac{dT^{core}}{dt} = S_{core} F_{core}, \qquad (22)$$

where $V_{core}$ is the core volume, $S_{core}$ is the core surface, $\rho_{Fe}$ is the core density, $C_{p,Fe}$ is the core heat capacity and $F_{core}$ is the heat flux through the CMB. This formulation allows to follow the evolution of core temperature as a function of time, based on the CMB heat flux.

## 2.4. Numerical model

Eq. (1) is discretized using a semi-implicit predictor-corrector Finite Difference scheme, of second-order in both space and time

<figure><img src="image_281.png" /><figcaption>**Fig. 2.** Temperature evolution from an initially adiabatic temperature profile with $T_p = 3200$ K and $T_0^{core} = 5000$ K. The liquidus and solidus used in our models are those obtained for the A-chondritic model and are represented respectively with red and green curves. In this model $e_{TBL} = 1$ m and $\eta_m = 100$ Pa s. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

(Press et al., 1993). Our scheme was successfully benchmarked against steady and unsteady analytical solutions for diffusion problems (Crank, 1975). The mantle is discretized using $n = 2900$ equally spaced grid points resulting in a constant spatial resolution $\Delta r = 1$ km. Non-linear effects are handled via a fixed-point/picard iteration procedure. The variable time step is determined as $\Delta t = \min(\Delta r^2/\kappa)$, where $\kappa(r) = k/(\rho C_p)$ is the effective diffusivity.
<figure><img src="image_282.png" /><figcaption>**Fig. 3.** Left panel: Melt fraction evolution from an initially completely molten magma ocean and corresponding to the case illustrated in Fig. 2. A melt fraction of 0.4 is a major discontinuity for the magma ocean viscosity (see text). Right panel: same with a F-peridotitic model and $T_0^{core} = 7000$ K.</figcaption></figure>

<figure><img src="image_283.png" /><figcaption>**Fig. 4.** Left panel: Time evolution of the magma ocean thickness (where the melt fraction is larger than 50%) for different initial core temperatures $T_0^{core}$ and different initial compositions (with $e_{TBL} = 1$ m and $\eta_m = 100$ Pa s). Right panel: Time evolution of the core temperature for different initial core temperatures and different initial compositions (with $e_{TBL} = 1$ m, $\eta_m = 100$ Pa s). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

## 3. Results

### 3.1. Thermal evolution of a deep magma ocean

We follow the thermal evolution of a deep magma ocean with an initially adiabatic temperature profile with $T_p = 3200$ K. Before we study the effect of each magma ocean parameter, we consider the following model as a reference case: a A-chondritic model, a 1 m thick bottom thermal boundary layer ($e_{TBL}$), a melt viscosity value of $\eta_m = 100$ Pa s and an initial core temperature $T_0^{core} = 5000$ K. Fig. 2 shows that the temperature rapidly decreases from the surface where heat is efficiently removed by radiative cooling even if a thin solid crust is formed within this upper thermal boundary layer. In the deepest part of the mantle, the temperature profile remains adiabatic but the cooling is slower. After 5000 yr, solidification occurs from the CMB where the liquidus is steeper than the adiabatic profile. As cooling proceeds, the melt fraction decreases and the last parcel with 100% melt starts to solidify in the upper mantle (see Fig. 3, left). Finally, the whole magma ocean drops below a 50% melt value in $t_{MO} \sim 150$ kyr with $t_{MO}$ being the magma ocean lifetime. This time ranges between the two characteristic timescales mentioned in Solomatov (2000): $10^3$ yr when crystallization starts from the bottom and $10^8$ yr when crystallization of the last drop of melt occurs in the shallow magma ocean. Clearly, our magma ocean lifetime is much shorter than the cooling timescale of ~4 Gyr proposed by Labrosse et al. (2007). This is due to the fact that, in our model, the solidification occurs from

the bottom-up which prevents our magma ocean from being overlaid by a thick insulating solid mantle. In Abe (1997), the magma ocean was restricted to a 1000 km-deep domain and the melting curves were steeper than the one used in our model. Hence, within 150 kyr, most of the mantle temperature profile computed from Abe (1997) is well below the solidus.

We monitored the thickness of the magma ocean (i.e. the thickness of the material having a melt fraction larger than 50%) as a function of time. Fig. 4 (left, black line) shows that after a short period (~20 kyr) where the mantle remains mostly molten, the thickness of the magma ocean rapidly decreases from 2900 km to 200 km, with a change of slope for a thickness of ~2000 km for the A-chondritic model. At the change of slope, the melt fraction reaches the critical value of 40% in the lowermost mantle, which induces an abrupt increase of its viscosity. Then the melt fraction progressively decreases at all mantle depths until we stop our simulations when the magma ocean thickness drops below 100 km. Fig. 4 (right, black line) shows the core temperature as a function of time for our reference case. In this model, the core temperature decreases monotonically from its initial value of 5000 K to a value of 4430 K, which corresponds to a melt fraction $\phi \sim 40\%$ at $P = 140$ GPa (i.e. when the abrupt change in viscosity occurs).

We monitored the time evolution of the heat flow coming at the CMB from the core $F_{core}$ and the heat loss at the surface $F_{surf}$ (Fig. 5, a). During the first 20 kyr, the heat flow at the CMB rapidly increases from $\sim 10^{16}$ to a nearly constant value of $\sim 10^{17}$ W. In the meantime the surface heat flow decreases from
<figure><img src="image_284.png" /><figcaption>**Fig. 5.** (a): Time evolution of the CMB (black) and surface (red) heat flows for a A-chondritic model (assuming $e_{TBL} = 1$ m, $\eta_m = 100$ Pa s and $T_0^{core} = 5000$ K except for the black dashed line where $T_0^{core} = 7000$ K). (b): Same with $e_{TBL} = 1000$ m. (c): Same with $e_{TBL} = 100$ mm. (d): Same with $e_{TBL} = 1$ mm. In all these models, the magma ocean lifetimes (time at which the plots ends up) are very close. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

~10¹⁹ (in agreement with the value proposed by Solomatov, 2000) to ~10¹⁷ W. Then, both the surface and the core heat flows decrease down to a value of ~10¹⁶ W until the end of the magma ocean stage (within $t_{MO} = 150$ kyr).

## 3.2. Influence initial of core temperature

We monitored the influence of the initial core temperature on the cooling of the deep magma ocean considering two values: (1) $T_0^{core} = 5000$ K, which is equal to the temperature at the bottom of the magma ocean for $T_p = 3200$ K in the A-chondritic model and (2) $T_0^{core} = 7000$ K, which corresponds to a core significantly hotter than the mantle. The comparison of the black and green curves in Fig. 4 (left) shows that the initial temperature has a negligible influence on the evolution of the magma ocean thickness as well as on its cooling timescale. When $T_0^{core} = 7000$ K, the initial core heat flow is large (~10¹⁸ W) (Fig. 5, a) because of the initial temperature contrast (= 2000 K) with the lowermost mantle. Then the heat flow decreases as a consequence of both the progressive solidification of the overlying mantle and the core heat depletion. When $T_0^{core} = 5000$ K, the initial CMB heat flow is much smaller than when $T_0^{core} = 7000$ K, but it increases rapidly due to rapid cooling of the lowermost mantle. In both cases, the heat flow decreases when the lowermost mantle has cooled sufficiently to reach the critical melt fraction of 40%. A significantly hotter initial core leads to an increase of the core heat flow, by a factor of 3 between 5000 K and 7000 K (Fig. 5, a). However, this increase of the initial core temperature, which should lead to an increase of the

surface heat flow, is not visible on the surface heat flow evolution because of the stronger efficiency of the surface cooling.

Also, a core initially 2000 K hotter than the lowermost mantle ends up 170 K hotter at the end of the MO stage (Fig. 4, right). When $T_0^{core} = 7000$ K the final core temperature is $T^{core} = 4600$ K which is slightly larger than the temperature at which the melt fraction of the lowermost mantle reaches the 40% critical value in the A-chondritic model.

## 3.3. A-chondritic vs. F-peridotitic model

Here we compare the evolutions of temperature and melt fraction between the A-chondritic and F-peridotitic models (using the corresponding liquidus and thermodynamical parameters). We focus on cases where $T_0^{core} = 7000$ K (with $T_0^{core} = 5000$ K, the core temperature would be lower than the liquidus for the F-peridotitic model at the CMB pressure). Fig. 4 (left) shows that for an initial core temperature $T_0^{core} = 7000$ K, the magma ocean thickness decreases more rapidly in the F-peridotitic model (red curve) than in the A-chondritic model (green curve). This is the direct consequence of a liquidus being significantly higher for a F-peridotitic model than for a A-chondritic model (Fig. 1). Hence, during the cooling of an initial fully molten magma ocean, the onset of mantle crystallization occurs earlier and the melt fraction decreases more rapidly in the F-peridotitic case (Fig. 3, right). The comparison of the A-chondritic and the F-peridotitic models shows a peak of the melt faction in the latter case occurring at a depth of ~600 km. This corresponds to the important slope change in the F-peridotitic liquidus that occurs at 20 GPa (see Fig. 1, right). Since the melt