# Stability of high-temperature salty ice suggests electrolyte permeability in water-rich exoplanet icy mantles

Jean-Alexis Hernandez, Razvan Caracas, Stéphane Labrosse

► To cite this version:

Jean-Alexis Hernandez, Razvan Caracas, Stéphane Labrosse. Stability of high-temperature salty ice suggests electrolyte permeability in water-rich exoplanet icy mantles. Nature Communications, 2022, 13, 10.1038/s41467-022-30796-5. insu-03748521

HAL Id: insu-03748521

https://insu.hal.science/insu-03748521

Submitted on 9 Aug 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
ARTICLE

https://doi.org/10.1038/s41467-022-30796-5

OPEN

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

Check for updates

# Stability of high-temperature salty ice suggests electrolyte permeability in water-rich exoplanet icy mantles

Jean-Alexis Hernandez ID 1,2,3✉️, Razvan Caracas ID 2,3,4 & Stéphane Labrosse ID 2

Electrolytes play an important role in the internal structure and dynamics of water-rich satellites and potentially water-rich exoplanets. However, in planets, the presence of a large high-pressure ice mantle is thought to hinder the exchange and transport of electrolytes between various liquid and solid deep layers. Here we show, using first-principles simulations, that up to 2.5 wt% NaCl can be dissolved in dense water ice at interior conditions of water-rich super-Earths and mini-Neptunes. The salt impurities enhance the diffusion of H atoms, extending the stability field of recently discovered superionic ice, and push towards higher pressures the transition to the stiffer ice X phase. Scaling laws for thermo-compositional convection show that salts entering the high pressure ice layer can be readily transported across. These findings suggest that the high-pressure ice mantle of water-rich exoplanets is permeable to the convective transport of electrolytes between the inner rocky core and the outer liquid layer.

<footnote n="1">¹European Synchrotron Radiation Facility, Grenoble, France. ²CNRS, Ecole Normale Supérieure de Lyon, Université de Lyon, Laboratoire de Géologie de Lyon LGLTPE UMR 5276, Lyon 69364, France. ³Centre for Earth Evolution and Dynamics, University of Oslo, Oslo 0315, Norway. ⁴Université de Paris, Institut de Physique du Globe de Paris, CNRS, 1 rue Jussieu, Paris 75005, France. ✉email: jean-alexis.hernandez@esrf.fr</footnote>
Water-rich exoplanets and icy moons are promising environments for the development of biological reactions. The extent and the composition of a potential ocean depends on the solubility and the availability of the various electrolytes. Both in situ measurements by planetary missions and aqueous alteration processes observed in chondrites suggest that salts based on Na⁺, Mg²⁺, Cl⁻ and SO₄²⁻ are present in extraterrestrial planetary aqueous environments¹⁻⁶. While an interface between liquid water and a silicate mantle favors hydrothermal reactions releasing nutrients in the ocean, the presence of an underlying high-pressure ice mantle might prevent the transport of electrolytes⁷ due to their insolubility in ice VI⁸. Nevertheless molecular ice VII has been shown to include non-negligible amounts of halides within their structures at conditions similar to large icy moon interiors⁸⁻¹⁴. These experimental results have recently been confirmed by the discovery of salty ice VII inclu- sions trapped in diamonds formed in the Earth mantle¹⁵.

The ability of ice VII to retain salts in its structure results from its body-centered cubic (bcc) sub-lattice of O atoms, which is a common characteristic of all ice phases between 2 GPa and 300 K to at least 270 GPa and 2000 K¹⁶. At these conditions, bcc H₂O ice undergoes multiple transitions due to changes in the dis- tribution and the dynamics of the H atoms. Below 30–40 GPa and 800 K, ice VII (space group *Pn-3m*) is the stable phase. It has a static orientational disorder of the H₂O molecules. The H atoms occupy the half-diagonals of the bcc O sub-lattice and the arrangement of H₂O molecules form a strong hydrogen-bond (HB) network with O-H···O bonds, where O-H denotes the intramolecular covalent bond and H···O the intermolecular bond. Increasing pressure in ice VII results in a progressive dis- sociation of the H₂O molecules due to the symmetrization of the O-H···O bonds. The compression first activates a dynamic translational disorder of H atoms along the O-H···O bonds (O- H···O↔O···H-O), this is ice VII'. Above 60–80 GPa, the symmetrization is completed, the H atoms are located at mid- distance between the two O atoms of the bcc half-diagonal, and the H₂O molecules are fully dissociated; this ionic ice is labelled ice X¹⁷,¹⁸. Increasing the temperature in ice VII activates a dynamic orientational disorder of the H₂O molecules¹⁹. Increasing both pressure and temperature leads to the combina- tion of both orientational and translational motions of the H atoms and results in an onset of H diffusion characterizing the superionic regime¹⁶,¹⁹⁻²⁶.

These transitions strongly influence the compressibility, elas- ticity and transport properties of bcc ice¹⁹,²³,²⁴,²⁷⁻³⁰. For instance, the transition between superionic ice VII'' and ice VII' is asso- ciated with a first-order transition, a dip in the elastic constants and a strong decrease of the electrical conductivity¹⁹,²⁴. The completion the hydrogen bond symmetrization makes ice X stiffer than other bcc ices and preliminary investigation of its rheology predict a change of slip-system around 250 GPa²⁸.

The inclusion of salt in ice VII decreases its melting temperature³¹ and breaks the connectivity of its hydrogen bond network¹¹,¹³. In particular, the presence of halides increases the rotational disorder of water molecules and decreases the role of quantum effects in the O-H···O bond symmetrization³², which shifts the transition to ice X by about 40 GPa at ambient temperature¹². Nevertheless, the stability of salt-bearing ice is not known at conditions typical of water-rich super-Earths, i.e., up to a few megabars and few thousand Kelvin. The inclusion mechanism of NaCl in bcc ice remains debated and its influence on the development of the superionicity have not been investigated yet.

Here we compute the solubility and the properties of NaCl- bearing water ice along the 1600 K isotherm at relevant condi- tions for the interiors of water-rich exoplanets presenting a high-

pressure ice mantle. Our approach is based on first-principles molecular dynamics simulations (DFT-MD) and extensive ther- modynamic modeling. Then, based on our results and analytical considerations, we investigate how the presence of salt would influence the convection and the transport of electrolytes in the high-pressure ice layer of large water-rich exoplanets.

## Results and discussion

**Inclusion mechanism for Na⁺ and Cl⁻.** Several experimental studies⁹,¹³,³¹ suggested three different and competing structural mechanisms for the inclusion of the Na⁺ and Cl⁻ ions in the bcc structure of ice VII. Based on experiments on NaCl-bearing ice, Frank⁹ inferred that both Na⁺ and Cl⁻ ions occupy octahedral face-centered sites. Based on experiments on LiCl-bearing ice, another group¹¹,¹²,³³ proposed that Cl⁻ substitutes a water molecule and Na⁺ remains in an octahedral site. Based on energy and size considerations, the same group suggested that both Na⁺ and Cl⁻ substitute H₂O molecules on the bcc sites (double sub- stitution mechanism). X-ray diffraction (XRD) measurements of natural (K,Na)Cl-bearing ice VII at ambient temperature¹⁵ revealed the same substitution sites on the bcc lattice.

Consequently, as a prerequisite to our study, we determine the inclusion mechanism for Na⁺ and Cl⁻ ions in cubic water ice at 1600 K. We perform a direct crystallization of a salty water solution. We start with two equilibrated solutions in the fluid domain at 1600 K and 30 GPa, one NaCl·126H₂O and one NaCl·30H₂O (see Table 1 for notation and transformation between wt% vs molar%), and directly compress them at about 70 GPa, in the solid domain. The high-pressure simulations last for 12 picoseconds and both solutions crystallize into a bcc ice with both Na⁺ and Cl⁻ occupying the sites of the water molecules (Fig. 1). While the stoichiometry of the solutions may have favored the bcc arrangement, the outcome of these simulations could have been also very different, including a glassy state. But the fact that we eventually obtain bcc ice adds a further confirmation of the double substitution mechanism for incorporating salt impurities into dense ices, with both Na and Cl substituting H₂O molecules.

**Properties and phase boundaries of the salty ice.** At each volume, we perform several DFT-MD in the canonical ensemble (volume and temperature fixed) with different configurations of Na and Cl in the bcc O sub-lattice and the resulting properties are averaged over the different configurations.

First, we evaluate the effect of NaCl inclusion on the bcc sub- lattice. The analysis of the radial pair distribution functions reveals that the substitution of water molecules by Na⁺ and Cl⁻ ions distorts the bcc sub-lattice (Fig. S13). At fixed volume (a = 2.525 Å), the shortest O-Na distance (i.e., the half-diagonal of the bcc unit-cell) is ~10% shorter than the O-O distance due to the stronger attraction between the two oppositely charged ions.

<table><caption>Table 1 The equivalence between mole% and wt% for the various NaCl-bearing H₂O terms used in the simulations denoted as NaCl·R·H₂O.</caption><thead><tr><th>Mole fraction NaCl</th><th>Mass fraction NaCl</th><th>R</th><th>Supercell</th><th>N<sub>NaCl</sub></th></tr></thead><tbody><tr><td>0.004</td><td>0.013</td><td>248</td><td>5 × 5 × 5</td><td>1</td></tr><tr><td>0.008</td><td>0.025</td><td>126</td><td>4 × 4 × 4</td><td>1</td></tr><tr><td>0.016</td><td>0.050</td><td>62</td><td>4 × 4 × 4</td><td>2</td></tr><tr><td>0.019</td><td>0.059</td><td>52</td><td>3 × 3 × 3</td><td>1</td></tr><tr><td>0.032</td><td>0.098</td><td>30</td><td>4 × 4 × 4</td><td>4</td></tr><tr><td>0.067</td><td>0.189</td><td>14</td><td>4 × 4 × 4</td><td>8</td></tr><tr><td>0.149</td><td>0.351</td><td>6</td><td>4 × 4 × 4</td><td>16</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Crystallization of a NaCl-126H₂O solution during a DFT-MD simulation at 70 GPa. The simulations are done in the NVT ensemble (constant number of particles, volume and temperature). During crystallization both Na⁺ and Cl⁻ ions occupy body-centered cubic lattice sites and both substitute a water molecule.</figcaption></figure>

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

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

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

<figure><img src="image_7.png" /><figcaption>**Fig. 2** H-bonding regimes in NaCl-bearing bcc ice at 1600 K. **a** Pressure dependence of the shortest O-H bond length. For a given NaCl concentration, OH reaches a maximum at the transition towards ice X. See text for more details. **b** OH-rotation rates (dots) and H-translation rates (squares) under compression in salty ice. **c** H diffusion coefficients ($D_H$). The inset represents a zoom in the 30-100 GPa pressure range showing the increasing pressure of the onset of diffusion with increasing NaCl content. When non-visible error bars are smaller than the size of the markers. **d** Stereographic projection of the density distribution of O-H bond orientations in ice X. In pure and salty, the four maxima corresponds to the (111) directions as expected for O-H bonds aligned with the half-diagonals of the bcc sub-lattice. In salty ice the orientation distribution of the O-H bonds is broader than in pure ice, indicating that H atoms tend to delocalize more out of the (111) direction (OO axis) than in pure ice. The color scale is the same for the two distributions and the white lines contours the same density of O-H bond orientations.</figcaption></figure>

O and Cl, having charge with the same sign, repulse each other resulting in longer O-Cl bonds. At low densities (2.361g.cm⁻³) and for a given NaCl concentration (NaCl·30H₂O), the O-Na and the O-Cl distances respectively differ from the O-O distance by -14% and by +6%. The differences change to -9.5% and +8% at higher densities (4.221g.cm⁻³). These distortions of the bcc sub-lattice are expected to decrease the elastic constants of NaCl-bearing ice compared to pure water ice.

As mentioned above, the various phases of bcc ice differ by the distribution and the dynamics of the protons (H atoms). To identify

the different regimes, we perform a three-body analysis¹⁹,²⁴ of the O-H···O bonds based on the atomic trajectories.

The comparison of the shortest OH distance to the shortest OO distance allows us to monitor the symmetrization of the O-H···O bonds, and the transition to ice X when OH = OO/2. Figure 2a shows that increasing concentration of NaCl shifts the transition to ice X towards higher pressures compared to pure ice. For instance, there is an upwards shift of 50 GPa in the ice X transition pressure when adding 35.1 wt% NaCl.

Counting the number of H-translations (O_a-H···O_b → O_a···H-O_b) and OH-rotations (O_a-H···O_b → O_a-H···O_c)
allows us to quantify dynamical disorders and to distinguish all other bcc phases in our simulations. Ice VII is not dynamically disordered and both OH-rotation and H-translation rates would be equal to zero. Plastic ice (free rotations of H₂O molecules) would have a positive OH-rotation rate and a H-translation rate equal to zero. Non-diffusive ices VII' and X would have H-translations without any OH-rotation. The presence of both H-translations and OH-rotations indicates that H atoms are diffusing: this characterizes the superionic regime. In this case, superionic ice VII'' is distinguished from superionic ices VII' and X by a OH-rotation rate greater than the H-translation rate. Figure 2b shows that in pure H₂O ice both H-translations and OH-rotations are activated up to ~150 GPa, indicating that H atoms diffuse below this pressure at the time-scale of our simulations and that pure ice is in a superionic regime. Below 52 GPa, pure ice is affected by more OH-rotations than H-translations and is in the superionic VII'' phase. Above 150 GPa, only H-translations are observed in pure ice indicating that ice is in its non-diffusive ice X form (OH = OO/2 at these conditions).

At a given pressure, increasing concentrations of NaCl strongly enhances the orientational disorder in bcc ice (increase of the OH-rotation rate and decrease of the H-translation rate). As a consequence, the pressure at which the OH-rotation rate equals the H-translation rate, i.e., the transition pressure between superionic ice VII'' and superionic ice VII' (or superionic ice X depending on the relation between the OH and OO distances), increases with the NaCl concentration. Also, above 2.5 wt% NaCl, we found that the OH-rotation rate never reaches zero, indicating that salty ice is superionic at all sample pressure at 1600 K. The enhanced orientational disorder in salty ice is caused by both the inability of Na to form any bonds with H, and by the absence of covalent Cl-H bonds. Only O-H···Cl bonds exist and present a strong asymmetry compare to O-H···O bonds as both the O-H bond and is H···Cl bonds are ~20-25% shorter than in normal O-H···O bonds (see Fig. S14).

These changes in bonding regime are reflected in the variation of the proton diffusion (Fig. 2c) and electrical conductivity $\sigma_e$ (see Fig. S15). In particular, superionic ice VII'' is characterized by the highest diffusion coefficients with a plateau at 1.5–2 m².s⁻¹ ($\sigma_e$=40–60 S.cm⁻¹) in both pure and salty ices. Under pressure, $D_H$ and $\sigma_e$ decrease down to statistically zero (i.e., $D_H$ < 10⁻¹⁰ m².s⁻¹) in ice with NaCl concentrations lower than 5.0 wt% although both OH-rotational and H-translational disorder still persist (Fig. 2b). For large NaCl concentrations ($R$=14 and 6), $D_H$ presents a second maximum at 100–120 GPa related to the enhanced H delocalization out of the O–H···O bonds (Fig. 2d).

In summary, the inclusion of salt distorts the bcc sub-lattice, increases the orientational disorder and weakens the hydrogen-bond network of the ice. This shifts the superionic VII'' - VII' - X transition sequence towards higher pressure with increasing salt concentrations (see Fig. 3). Therefore, at given pressure-temperature conditions, salty ice is less stiff, has a greater electrical conductivity and is likely to be less viscous than pure ice.

**Solubility of NaCl in high-temperature water ice.** The entire phase diagram of ice at high pressure and temperature is affected by the presence of salt, with the phase boundaries pushed to higher pressures (Fig. 3). The amplitude of the shift depends on the amount of salt that can actually be dissolved into the ice. Its determination needs to pass through the evaluation of the Gibbs free energies for the ice-salt solid solutions. This involves the estimation of the mixing terms, which can considerably lower the energy of the intermediate terms with respect to the mechanical mixture of the pure end-member terms.

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

**Fig. 3 Phase diagram of bcc water ices.** Solid-solid phase boundaries in pure ice are indicated by black lines. Corresponding phase boundaries in NaCl-bearing ice (NaCl·RH₂O) ice are shown at 1600 K and 300 K. Grey areas indicate where pure H₂O is superionic, with darker areas corresponds to regions with higher H-diffusion coefficients. The phase boundaries in pure H₂O ice are drawn based on our DFT-MD simulations¹⁹,²⁴ and on both experimental⁶⁴ and theoretical works¹⁸,⁶⁵ that take into account quantum effects at 300 K. The melting line of pure ice as obtained in our single-phase simulations is shown by the red dashed line. Experimental melting points from Schwager et al.⁶⁶. and Queyroux et al.²⁶. are shown as red circles and red empty squares respectively. Black empty squares correspond to the experimental determination of first-order transition between superionic ice VII'' and ice VII'²⁶, and black dots refer to a solid-solid transition⁶⁶, likely corresponding to the same transition. Blue dashed lines show adiabatic profiles in the H₂O layer of hypothetical ocean exoplanets with habitable pressure-temperature conditions in the surface ocean and no or thin atmosphere (see Supplementary Information). The different profiles are anchored at surface temperatures of 273 K (temperature of crystallization of pure ice Ih), 300 K and 350 K (representative temperature for an optimal growth of thermophile and hyperthermophile organisms⁶⁷) respectively. Circled numbers 1, 2, 3 and 4 indicate approximate pressures at the bottom of the high-pressure ice mantle for 1, 2, 3, and 4 $M_\oplus$ planets (i.e., ~1.3, 1.5, 1.7, 1.9 $R_\oplus$) with 50 wt% H₂O and Earth-like cores as calculated in Sotin et al. (2007) for a surface temperature of 300 K.

In solid solutions, the configurations that present the major contribution are generally separated by high energy barriers although their own energies can be very close or similar. As such, the transition paths between configurations cannot be sampled directly by standard molecular dynamics. Then statistical methods need to be employed to estimate the configuration entropy, as described in the methodology section, which depends on the number of configurations accessible to the system, i.e., the Na, Cl, and O arrangements in the bcc sub-lattice. As we investigate the solubility and properties of salty ice along the 1600 K isotherm, we validate the statistical approach by a second set of calculations exploring all possible arrangements for two concentrations at 100 GPa and 1600 K.

At 1600 K and 100 GPa, the Gibbs free energy of mixing is negative up to 5–10 wt% NaCl depending on the estimation of the configuration entropy and presents a minimum at 2.5 wt% NaCl (NaCl·126H₂O). Assuming that H₂O-bearing B2-type NaCl is the only phase that forms close to the NaCl end-member term, and that its free energy of mixing has only a small negative incursion, the common tangent to the two free energy curves intersects the curve of NaCl-bearing H₂O free energy of mixing close to its minimum. This suggests that the inclusion of NaCl in water ice at 1600 K between 80 and 200 GPa is thermodynamically favorable up to ~2.5 wt% NaCl. At these conditions, salty superionic ice always presents both a larger volume (by +0.41% for 2.5 wt% NaCl and by +0.86% for 5.88 wt% NaCl, both at 1600 K and
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 4** Thermodynamics of NaCl-RH₂O mixture at 1600 K and 100 GPa.</figcaption></figure>

a Gibbs free energy of mixing $\Delta g_{mix}$ as computed from a fixed pressure approach and a complete sampling of the reduced configuration space (blue circles). In this case the configuration entropy is exact. Black (with correction for the GGA volume overestimation in B2 NaCl) and red (without) lines correspond to $\Delta g_{mix}$ calculated from the fixed volume approach (see text and Supplementary Information) and a random sampling of the entire configuration space. Due to the limited sampling, the configuration entropy is less constrained in this approach. Inset: $\Delta g_{mix}$ at 200 GPa and 1600 K. b Relative volume (blue) and density (red) differences between NaCl·RH₂O and H₂O. Crosses corresponds to individual configurations and pluses to Boltzmann averages obtained from the complete sampling of the reduced configuration space. Lines correspond to the model derived from the fixed volume approach.

100 GPa) and a larger density than pure water ice (Fig. 4a). Such amounts are in agreement with recent experimental and theoretical works performed on the inclusion of various halides (NaCl, LiCl, MgCl₂, RbI) in ice VII at lower pressures and temperatures⁸,⁹,¹¹⁻¹⁴. A solubility of 2.5 wt% is still consistent with the observation of temperature-induced exsolution of some NaCl above 500 K in NaCl-doped ice VII whose initial concentration was higher¹⁰. The volume increase due to the addition of salt contrasts with the XRD measurements of Frank et al.⁹,¹⁰, but agrees with all previously mentioned studies⁸,¹¹⁻¹³.

**Transport of electrolytes in a high-pressure ice mantle.** Salts are observed at the surface of small icy bodies of the Solar system and may be produced from hydrous alteration processes taking place before and during the accretion, and/or from the interactions between the water-rich layer and a rocky core. Although our DFT

calculations are done at 1600 K, the applications of our results are not restricted to this exact temperature. Indeed, the combination of our high-pressure and high-temperature results with previous studies done at milder conditions⁸,¹³ shows that bcc ices are able to include non-negligible amounts of salt in their structure over a broad range of conditions relevant for the interiors of water-rich super-Earths³⁴ ($1R_⊕ < R < 1.8R_⊕$), mini-Neptunes ($1.8R_⊕ < R < 4R_⊕$), and Neptune-like exoplanets³⁵,³⁶, according to the classification made in ref. ³⁷,³⁸, based on the Kepler exoplanet catalogue.

As shown above, the inclusion of salt affects the phase diagram and the transport properties of bcc ices (e.g., extension of the superionic field and of the transition pressures between the different phases) and increases the density by about 0.5% compared to pure ice, as the positive change in specific volume due to the addition of salt is compensated by its higher molecular weight. These changes are expected to influence the dynamics of a high-pressure ice mantle, in a variety of different ways, depending on the planet characteristics. A full understanding of salt transport between the rocky core and the ocean by convection in the high-pressure ice layer requires complex fluid dynamical double-diffusive models including complexities such as phase changes, pressure-, temperature- and composition-dependence of physical properties, and possibly partial melting³⁹⁻⁴¹. This falls clearly beyond the scope of the present paper but, the feasibility of the process can be evaluated using the current knowledge on mantle convection, which includes similar ingredients⁴².

Given that salts and other electrolytes are a necessary (but of course not sufficient) condition for habitability in a (sub-)surface ocean, it is interesting to evaluate if salt can be involved in a global cycle, i.e., if salt and other electrolytes could be transported throughout the high-pressure ice mantle by convection. In the following, we develop our reasoning and provide numerical applications to a hypothetical planet of 1 $M_⊕$ with 50 wt% H₂O, an Earth-like core, absent or very thin atmosphere, and a surface temperature of 300 K (case 1 in Fig. 3) allowing for the existence of a surface ocean with mild enough temperature conditions to sustain life as we know it⁴³. Such case (~1.3 $R_⊕$) is chosen to fit within the possible range of super-Earths ($R < 1.6 - 1.8R_⊕$). It is depicted in Fig. 5, and presents a direct ocean/ice VII boundary as ice VI melts below 350 K. Planets with radii above the so-called radius valley ($R > 1.8R_⊕$) are thought to have thick atmospheres which would result in higher temperatures in the ocean³⁷,³⁸. However, it is worth mentioning that the case of larger water-rich super-Earths with similar surface ocean conditions cannot be excluded. Exoplanets with orbital periods longer than 100 days and receiving less than 10 times the Earth's Solar flux are not observed because of biases in observations (see Supplementary Information for a discussion on statistical results based on the Kepler catalogue). It is depicted in Fig. 5, and presents a direct ocean/ice VII boundary as ice VI melts below 350 K.

Considering an adiabatic profile in the H₂O layers (see Fig. 3, and Supplementary Information for calculation details), we find that the temperature at the bottom of the ocean is 367 K, resulting in a liquid/ice VII boundary at the bottom of the ocean and about 800 K at the bottom of the ice mantle in the absence of thermal boundary layer with the underlying rocky mantle. This indicates that the whole high-pressure ice mantle is composed of bcc ice as depicted in Fig. 5. Thermal convection in a high pressure ice layer depends on the value of the dimensionless Rayleigh number which, in the case of an imposed heat flux $q_c$ from the rocky core, takes the form³⁹:

$$Ra = \frac{\alpha \rho g q_c d^4}{k \kappa \eta} = 3 \times 10^{13} \left(\frac{R_c}{5000 \text{km}}\right) \left(\frac{q_c}{10 \text{mW/m}^2}\right) \left(\frac{d}{3000 \text{km}}\right)^4 \left(\frac{\eta}{1.3 \times 10^{17} \text{Pas}}\right)^{-1} \quad (1)$$
<figure><img src="image_11.png" /><figcaption>Fig. 5 Transport of salt through the high-pressure ice mantle of a hypothetical water-rich exoplanet with $1M_{\odot}$, 50 wt% $H_2O$, and a surface temperature of 300 K (corresponds to case 1 in Fig. 3). A thermal plume creates a focused upward flux of salty ice that melts at the boundary with the ocean. The crystallization at the bottom of the ocean over a broad area produces a diffuse return flow of salt in the mantle. Depending on the initial conditions and on the partitioning coefficients of NaCl between the ice and the ocean, the icy mantle may acts either as a well or a source of electrolytes for the ocean. Label C illustrates a possible thermal boundary layer at the bottom of the icy mantle at the contact with the hotter rocky mantle.</figcaption></figure>

with $\alpha$ the coefficient of thermal expansion, $\rho$ the density, $g$ the gravitational attraction, $d$ the thickness of the ice layer, $\eta$ the viscosity, $k$ and $\kappa$ the thermal conductivity and diffusivity, respectively. However, transport properties for bcc ices are quite uncertain$^{28}$. To our knowledge no viscosity measurement in ice VII has been reported yet and we highly encourage any experimental effort in this direction. Here we use values quoted for ice VI from Choblet et al.$^{39}$, as a guideline and keep on the right-hand-side of Eq. 1 the parameters most uncertain or varying among planets. In the present case, geometric factors of the planet are from Sotin et al.$^{44}$, but these values can be changed depending on the planet of application. The gravity acceleration is computed assuming a rocky core of radius $R_c$ and density $\rho_c = 6 \times 10^3 \text{kg.m}^{-3}$. The reference values for the various parameters are quite conservative and the Rayleigh number is expected to be well above the critical value $Ra_c = 352$ for the onset of solid state convection in the presence of a solid-liquid phase change at the surface$^{45}$, validating the choice of an adiabatic temperature profile in the bulk of the ice layer. The lowermost part of the layer would see a steeper temperature gradient associated with the convective boundary layer. For such large values of Ra, the typical flow velocity in the ice layer scales as$^{46}$ $V_{rms} = 0.2Ra^{0.66}\kappa/d \sim 700\text{m/yr}$ for the reference values used above, which can very efficiently transport salty ice from the ice/core interface to the ocean. The most important shortcoming in these estimates is neglecting the effect of salts on density which should impede thermal convection if we consider that it enters at the bottom of the solid. This effect can be quantified by introducing the buoyancy number, $B = \Delta\rho_\chi/\rho\alpha\Delta T$, with $\Delta T$ the temperature difference driving thermal convection and $\Delta\rho_\chi/\rho$ the density increase owing to salt enrichment$^{47}$. A large value of B leads to a strong effect of composition and possible stratification, while efficient entrainment is possible for $B \ll 1$. The exact value for the transition in regimes has not been studied for the problem discussed here but is generally around 1. For convection in the solid ice layer surrounded by an ocean, $\Delta T$ is the temperature drop across the bottom boundary layer$^{46}$ and can be estimated from the scaling of heat transfer as function of the Rayleigh number (see Methods). The value of $\Delta\rho_\chi/\rho$ must be lower than

the solubility limit, $\delta\rho_m = 0.7\%$, but also depends on the mode of transfer from the rocky core. The largest flow of salt is likely to happen if ice melting occurs at the interface with the rocky core. In this case, salts may be transported to the ocean by porous flow$^{39-41}$. We here consider the situation, a priori less favorable for transfer of salts to the ocean, where no melting occurs and salts flow to the ice layer by diffusion, with a diffusion coefficient $D$. The relevant value entering $B$ is the density increase from salts in the thermal boundary layer, with a thickness $\sqrt{\kappa t_*}$ with $t_*$ the time needed to make the boundary layer unstable$^{48}$. In the same time, the concentration of salts increases in a layer of thickness $\sqrt{Dt_*}$ making the average relative density increase from salts in the thermal boundary layer equal to $\delta\rho_m \sqrt{D/\kappa}/2$. Combining all that (see Methods) leads to:

$$B = 3 \times 10^{-2} \left( \frac{D}{3.5 \times 10^{-12} \text{m}^2 \text{s}^{-1}} \right)^{1/2} \left( \frac{R_c}{5000 \text{km}} \right)^{1/3} \left( \frac{q_c}{10 \text{mWm}^{-2}} \right)^{-2/3} \left( \frac{d}{3000 \text{km}} \right)^{1/3} \left( \frac{\eta}{1.3 \times 10^{17} \text{Pas}} \right)^{-1/3} \quad (2)$$

Taking the same reference values as before and a tentative value for the maximum salt diffusivity equal to that of O atoms, $3.5 \times 10^{-12} \text{m}^2 \text{s}^{-1}$ (see Methods), leads to a very small buoyancy number which implies a very limited effect of salts against thermal convection. Moreover, the style of convection expected in an ice layer bounded above by an ocean shows focused hot plumes with a diffuse return flow$^{46}$; the up-welling material melts when it reaches the ocean$^{39,45}$, therefore releasing salts there (Fig. 5). The presence of a high pressure ice mantle in a water planet should therefore not be considered as a chemical insulator for the transport of NaCl and likely other electrolytes, towards a shallow ocean$^{39}$. We suggest future assessments of the habitability potential of large water-rich planets to consider the present results.

Both the efficiency and the duration of the electrolyte transport to the shallow ocean depend on the actual size, composition and surface temperature of the considered planet, which might result in different scenarios at the interface between the ocean and the ice mantle, and between the ice mantle and the rocky core$^{8,44,49,50}$. In the case of an exoplanet with a thick atmosphere, most of the ocean would reach a supercritical state, and the habitability would not be preserved. In such planets, the ocean would also be in contact with bcc ice (VII ou superionic VII" depending on the temperature) given that the total amount of $H_2O$ is large enough to reach the liquid/solid transition. If the surface temperature is very high, the conditions in the water layer may never intersect the melting line. In the case of an exoplanet with a very low surface temperature without atmosphere, a layer of ice VI could cap the ice VII layer, and the phase change VII to VI in up-welling plumes could lead to ex-solving NaCl to make a polycrystalline assemblage$^{8}$. The gross density of the assemblage should be similar to that of salt-doped ice VII and, according to the discussion surrounding Eq. (2), should not prevent convection to proceed.

Finally, depending on the planetary formation scenario, the initial electrolyte concentration in the water-rich envelopes may vary. If water is brought after the crystallization of the rocky mantle, water may remain relatively free from alteration products. If the accretion of the icy and rocky materials is synchronous, the water-based layers would likely contain a certain initial concentration of electrolytes due to the production of salts by hydrous alteration of chondritic material$^{6}$. In this case, the evolution of the ocean salinity would depend on the partitioning of salts between liquid water and ice during crystallization at the bottom of the ocean. The enrichment in salt of the icy mantle
would be influenced by the nature of the water/rock interface. In the case of a solid-solid interface, the enrichment of the ice layer depends on the partitioning coefficients of each chemical species between the silicates and the cubic ice. If partial melting occurs at the bottom of the high-pressure ice mantle due to the temperature increase in a thermal boundary layer, the incorporation of impurities in the ice could be much more efficient than in the case of a solid-solid interface³⁹⁻⁴¹. Future studies need to estimate the partitioning coefficients of essential electrolytes between liquid water, superionic ice, and silicates at conditions comparable to the base of the ice layer in ocean planets to improve our understanding of the internal functioning of these planets and to better assess their habitability.

## Methods

### First-principles molecular dynamics

**Computational details.** We study the behavior of NaCl-bearing dense H₂O ice by performing first-principles molecular dynamics (DFT-MD) simulations with the Vienna Ab initio Simulation Package⁵¹,⁵² (VASP v4.6 and v5.2). To model the atoms, we use the PAW potentials⁵³,⁵⁴ with the following electronic configurations: O:[He]2s²2p⁴, H:1s², Na:[Ne]3s¹ and Na:[He]2s²2p⁶3s¹, and Cl:[He]2s²2p⁵. The electronic structure is calculated within the density functional theory (DFT) framework and the Perdew-Burke-Ernzerhof general gradient approximation⁵⁵ is used to approximate the exchange-correlation term. Most DFT-MD simulations are performed with 4 × 4 × 4 bcc supercells (see Table 1) with a cutoff energy of 550 eV for the plane-waves that ensures that internal energy is converged within 5 meV per atoms, and 1.5 GPa precision in pressure with respect to a calculation with a cutoff energy of 1000 eV. We calculate electronic wavefunctions at the Γ-point of the Monkhorst k-point grid⁵⁶. Time steps of 0.5 or 1 fs are employed in all simulations. DFT-MD are either performed in the NVT, i.e., constant Number of particles, Volume, and Temperature, and NPT, i.e., constant Number of particles, Pressure, and Temperature, ensembles. For NVT sampling, temperature is controlled using velocity rescaling (isokinetic) as no difference was observed for averaged properties when using a Nose-Hoover thermostat⁵⁷,⁵⁸. Isobaric-isothermal simulations are performed with Raman-Parrinello molecular dynamics⁵⁹ and a Langevin thermostat.

**Determination of the free energy of mixing.** The estimation of the thermodynamic stability of the binary system requires the calculation of the Gibbs free energy of mixing $Δg_{mix}(w)$ as a function of salt concentration $w$, pressure $P$, and temperature $T$. $Δg_{mix}(w, P, T)$ is related to the Gibbs free energies of the pure phases ($g_{H₂O}(P, T)$, $g_{NaCl}(P, T)$) and to the configuration averaged Gibbs free energy of the mixture ($g_{sol}(w, P, T)$) as:

$$Δg_{mix}(w, P, T) = g_{sol}(w, P, T) - (wg_{NaCl}(P, T) + (1 - w)g_{H₂O}(P, T)) \quad (3)$$

$g_{sol}(w, P, T)$ includes a configuration entropy term, $s_{conf}(w, P, T)$.

First we compute the free energies of the two end-members: H₂O and B₂-type NaCl. We estimate the vibrational contribution to the entropy and the corrections for the nuclear quantum effects by computing the partial vibrational spectra from the velocity autocorrelation functions and using the two phases thermodynamic memory function model (2PT-MF)²³,⁶⁰⁻⁶². We estimate the configuration entropy for the H₂O-rich part of the solid-solution taking advantage of recent developments in the computation of mixture thermodynamics based on reduced configuration space sampling⁶³.

Then we perform a series of NVT simulations with different concentrations of NaCl at different volumes (Table 1); at each given volume we consider several but not all possible starting NaCl configurations that yield an average Helmholtz free energy $F$ from which we deduced the configuration entropy and weighted thermodynamic properties. From the resulting data points, we construct the Helmholtz free energy surface at 1600 K over the entire volume range of interest. By derivation we obtained the Gibbs free energy of the solution as function of pressure and NaCl concentration (see section S2).

Nevertheless as we sample a subset of the configuration space, the configuration entropy is likely underestimated. Second we perform NVT simulations at 100 GPa for all possible configurations at three NaCl concentrations ($R = 250, 126$ and 52) to have a correct value of the configuration entropy; the two approaches yield about 5% difference in the mass fraction of NaCl that can be dissolved in the ice.

The full description of the free-energy estimation and the fitting procedure to construct the free energy surface are available in Supplementary Information. A numerical implementation can be found at https://osf.io/w64fm/.

**Estimation of NaCl diffusivity.** The diffusivity of NaCl (D) in the ice is required to calculate the buoyancy number B of a thermal boundary layer enriched in salt at the bottom of the ice mantle. In order to provide a reasonable order of magnitude for this quantity, we consider that D is similar to the diffusivity of O atoms $D_O$, possible in our calculations only in the presence of vacancies. To estimate $D_O$, we

perform an additional DFT-MD run in a defective ice X structure with 126 H₂O and two H₂O vacancies. The simulation is done in the NVT ensemble (1600 K and a corresponding pressure of 110 GPa) with a time step of 1 fs. During the duration of this long simulation (97.5 ps of statistics after 2.5 ps of equilibration), we observe a total of 17 jumps, enough to obtain a linear regime in the mean square displacement (see Fig. S16) and a diffusion coefficient $D_O = 3.5 \times 10^{-12} m^2 \cdot s^{-1}$. Given that in natural ice the concentration of vacancies may be lower than 2/128, our estimation of $D_O$, and thus D, likely represents an upper bound.

**Geodynamical scalings.** We use scaling laws obtained in the context of the dynamics of the Earth mantle to estimate the feasibility of solute transport by convection in the high-pressure ice layer from the rocky core of icy exoplanets to the overlying ocean. The first important control dimensionless number for this process, the Rayleigh number (Eq. 1), depends on several ill-constrained and planet-dependent parameters. Some physical parameters for bcc ices are obtained directly from the DFT-MD calculations ($ρ = 2.0 \times 10^3$ kg m⁻³, $α = 2 \times 10^{-4}$ K⁻¹), some are educated guesses based on constraints on ice VI ($η = 1.310^{17}$ Pas, $κ = 4.310^{-7}$ m² s⁻¹, $k = 1.6$ Wm⁻¹ K⁻¹). The gravitational acceleration is computed as $g = 4πGρ_c R_c / 3$ with $G = 6.67 \times 10^{-11}$ m³ kg⁻¹ s⁻² the gravitational constant. Conservative reference values for the planet-dependent parameters are used and the most uncertain or variable parameters are emphasized on the right-hand-side of Eq. (1).

The second important dimensionless control number is the buoyancy number which balances the density change from composition and that from temperature. Thermal convection in the ice layer proceeds by diffusive growth of a boundary layer until its thickness $δ_T \sim \sqrt{kt}$ reaches a critical value, set by the value of the boundary layer Rayleigh number $Ra_δ = Ra(δ_T/d)^n$ reaching the critical value for the onset of convection. This provides a critical time, $t_c$. In the same time, the boundary layer gets enriched in salts by diffusion, on a thickness $δ_χ \sim \sqrt{Dt_c}$, with a saturated value $δ_p = 7 \times 10^{-3}ρ$ at the rock-ice boundary. The temperature and composition profiles in the boundary layer take error function shapes with their respective thicknesses as scaling parameter but can be linearized to estimate that the total mass of salts in the composition boundary layer is $δ_p = \sqrt{Dt_c}/2$ per unit surface, which amounts to an average density increase in the thermal boundary layer from salt enrichment $Δρ_χ = δ_p = \sqrt{D/κ}/2$. The value of D was estimated from DFT-MD as detailed in the previous section. To estimate the temperature difference $ΔT$ across the boundary layer, we use the scaling of the dimensionless heat flux as function of the Rayleigh number⁴⁶ as $Nu = 0.37Ra^{1/3}$ which provides a scaling for the thickness of the boundary layer as $d_{BL} = d/(0.37Ra^{1/3})$. The temperature difference is then $ΔT = q_c d_{BL}/k$. Collecting these relationships gives for the buoyancy number

$$B = \frac{\delta \rho_m}{2\rho} \sqrt{\frac{D}{\kappa} \frac{0.37Ra^{1/3}k}{d q_c \alpha}} \quad (4)$$

Combining this equation with the scaling for the Rayleigh number (Eq. 1) leads to Eq. (2). It shows that the density increase from enrichment in salts is not likely to affect convection in the ice layer, as long as the transport at the interface happens by diffusion. If melting of ice occurs, a much more efficient transport process is possible³⁹⁻⁴¹.

## Data availability

The data generated in this study are provided in the article, in Supplementary Information, and in a public repository at https://osf.io/w64fm/ together with a numerical implementation of the free energy surface of NaCl-bearing bcc ice. Raw data from atomistic simulations are available from the corresponding author upon request.

## Code availability

DFT simulations have been performed with the proprietary code VASP (https://www.vasp.at/).

Received: 9 June 2020; Accepted: 6 May 2022;

Published online: 21 June 2022

## References

1. Zolensky, M., Bodnar, R., Schwandt, C. & Yang, S. V. Halide minerals in the monahans (1998) and zag h chondrite regolith breccias. *Lunar Planet. Sci.* **31**, 1181 (2000).

2. Izawa, M., Nesbitt, H., MacRae, N. & Hoffman, E. Composition and evolution of the early oceans: evidence from the tagish lake meteorite. *Earth Planet. Sci. Lett.* **298**, 443–449 (2010).
3. McCord, T. et al. Salts on europa's surface detected by galileo's near infrared mapping spectrometer. *Science* **280**, 1242-1245 (1998).

4. McCord, T. B., Hansen, G. B. & Hibbitts, C. A. Hydrated salt minerals on ganymede's surface: evidence of an ocean below. *Science* **292**, 1523-1525 (2001).

5. Postberg, F., Schmidt, J., Hillier, J., Kempf, S. & Srama, R. A salt-water reservoir as the source of a compositionally stratified plume on encladus. *Nature* **474**, 620-622 (2011).

6. Zolotov, M. Y. Aqueous fluid composition in ci chondritic materials: chemical equilibrium assessments in closed systems. *Icarus* **220**, 713-729 (2012).

7. Noack, L. et al. Water-rich planets: How habitable is a water layer deeper than on earth? *Icarus* **277**, 215-236 (2016).

8. Journaux, B. et al. Salt partitioning between water and high-pressure ices. implication for the dynamics and habitability of icy moons and water-rich planetary bodies. *Earth Planet. Sci. Lett.* **463**, 36-47 (2017).

9. Frank, M. R. et al. Experimental study of the nacl-h₂o system up to 28 gpa: Implications for ice-rich planetary bodies. *Phys. Earth Planet. Interiors* **155**, 152-162 (2006).

10. Frank, M. R., Scott, H. P., Maglio, S. J., Prakapenka, V. B. & Shen, G. Temperature induced immiscibility in the nacl-h₂o system at high pressure. *Phys. Earth Planet. Interiors* **170**, 107-114 (2008).

11. Klotz, S., Bove, L. E., Strässle, T., Hansen, T. C. & Saitta, A. M. The preparation and structure of salty ice vii under pressure. *Nat. Mater.* **8**, 405-409 (2009).

12. Bove, L. E. et al. Effect of salt on the h-bond symmetrization in ice. *Proc. Natl Acad. Sci.* **112**, 8216-8220 (2015).

13. Ludl, A.-A. et al. Probing ice vii crystallization from amorphous nacl-d₂o solutions at gigapascal pressures. *Phys. Chem. Chem. Phys.* **19**, 1875-1883 (2017).

14. Watanabe, M., Komatsu, K., Noritake, F. & Kagi, H. Structural incorporation of mgcl2 into ice vii at room temperature. *Jpn. J. Appl. Phys.* **56**, 05FB03 (2017).

15. Tschauner, O. et al. Ice-VII inclusions in diamonds: evidence for aqueous fluid in Earth's deep mantle. *Science* **359**, 1136-1139 (2018).

16. Millot, M. et al. Nanosecond x-ray diffraction of shock-compressed superionic water ice. *Nature* **569**, 251-255 (2019).

17. Holzapfel, W. On the symmetry of the hydrogen bonds in ice vii. *J. Chem. Phys.* **56**, 712-715 (1972).

18. Benoit, M., Marx, D. & Parrinello, M. Tunnelling and zero-point motion in high-pressure ice. *Nature* **392**, 258-261 (1998).

19. Hernandez, J.-A. & Caracas, R. Proton dynamics and the phase diagram of dense water ice. *J. Chem. Phys.* **148**, 214501 (2018).

20. Cavazzoni, C. et al. Superionic and metallic states of water and ammonia giant planet conditions. *Science* **283**, 44-6 (1999).

21. Goldman, N., Fried, L. E., Kuo, I.-F. W. & Mundy, C. J. Bonding in superionic phase of water. *Phys. Rev. Lett.* **94**, 217801 (2005).

22. Goncharov, A. F. et al. Dynamic ionization of water under extreme conditions. *Phys. Rev. Lett.* **94**, 125508 (2005).

23. French, M., Desjarlais, M. P. & Redmer, R. Ab initio calculation of thermodynamic potentials and entropies for superionic water. *Phys. Rev. E* **93**, 022140 (2016).

24. Hernandez, J.-A. & Caracas, R. Superionic-superionic phase transitions in body-centered cubic h₂o ice. *Phys. Rev. Lett.* **117**, 135503 (2016).

25. Millot, M. et al. Experimental evidence for superionic water ice using shock compression. *Nat. Phys.* **1**, 297-302 (2018).

26. Queyroux, J.-A. et al. Melting curve and isostructural solid transition in superionic ice. *Phys. Rev. Lett.* **125**, 195501 (2020).

27. Sugimura, E., Komabayashi, T., Ohta, K., Hirose, K. & Ohishi, Y. Experimental evidence of superionic conduction in h₂o ice. *J. Chem. Phys.* **137** (2012).

28. Journaux, B. et al. Elasticity and dislocations in ice x under pressure. *Phys. Earth Planet. Inter.* **236**, 10-15 (2014).

29. French, M. & Redmer, R. Construction of a thermodynamic potential for the water ices vii and x. *Phys. Rev. B* **91**, 014308 (2015).

30. Sun, S., He, Y., Kim, D. Y. & Li, H. Anomalous elastic properties of superionic ice. *Phys. Rev. B* **102**, 104108 (2020).

31. Journaux, B., Daniel, I., Caracas, R., Montagnac, G. & Cardon, H. Influence of nacl on ice {VI} and ice {VII} melting curves up to 6 gpa, implications for large icy moons. *Icarus* **226**, 355-363 (2013).

32. Bronstein, Y. et al. Quantum versus classical protons in pure and salty ice under pressure. *Phys. Rev. B* **93**, 024104 (2016).

33. Klotz, S. et al. Ice vii from aqueous salt solutions: From a glass to a crystal with broken h-bonds. *Scientific Reports* **6**, 32040 (2016).

34. Léger, A. et al. A new family of planets? "ocean-planets". *Icarus* **169**, 499-504 (2004).

35. Zeng, L. et al. Growth model interpretation of planet size distribution. *Proc. Natl Acad. Sci.* **116**, 9723-9728 (2019).

36. Mousis, O. et al. Irradiated ocean planets bridge super-earth and sub-neptune populations. *Astrophysical J.* **896**, L22 (2020).

37. Fulton, B. J. et al. The california-kepler survey. iii. a gap in the radius distribution of small planets. *Astronomical J.* **154**, 109 (2017).

38. Lehmer, O. R. & Catling, D. C. Rocky worlds limited to 1.8 earth radii by atmospheric escape during a star's extreme uv saturation. *Astrophysical J.* **845**, 130 (2017).

39. Choblet, G., Tobie, G., Sotin, C., Kalousová, K. & Grasset, O. Heat transport in the high-pressure ice mantle of large icy moons. *Icarus* **285**, 252-262 (2017).

40. Kalousová, K., Sotin, C., Choblet, G., Tobie, G. & Grasset, O. Two-phase convection in Ganymede's high-pressure ice layer — Implications for its geological evolution. *Icarus* **299**, 133-147 (2018).

41. Kalousová, K. & Sotin, C. Melting in high-pressure ice layers of large ocean worlds - Implications for volatiles transport. *Geophys. Res. Lett.* **45**, 8096-8103 (2018).

42. Schubert, G., Turcotte, D. L. & Olson, P. *Mantle Convection in the Earth and Planets* (Cambridge University Press, Cambridge, 2001).

43. Cockell, C. S. et al. Habitability: a review. *Astrobiology* **16**, 89-117 (2016).

44. Sotin, C., Grasset, O. & Mocquet, A. Mass-radius curve for extrasolar earth-like planets and ocean planets. *Icarus* **191**, 337-351 (2007).

45. Labrosse, S., Morison, A., Deguen, R. & Alboussière, T. Rayleigh-bénard convection in a creeping solid with a phase change at either or both horizontal boundaries. *J. Fluid Mech.* **846**, 5-36 (2018).

46. Agrusta, R. et al. Mantle convection interacting with magma oceans. *Geophys. J. Int.* **220**, 1878-1892 (2019).

47. Le Bars, M. & Davaille, A. Stability of thermal convection in two superimposed miscible viscous fluids. *J. Fluid Mech.* **471**, 339-363 (2002).

48. Sotin, C. & Labrosse, S. Three-dimensional thermal convection of an isoviscous, infinite-Prandtl-number fluid heated from within and from below: applications to heat transfer in planetary mantles. *Phys. Earth Planet. Inter.* **112**, 171-190 (1999).

49. Fu, R., O'Connell, R. J. & Sasselov, D. D. The interior dynamics of water planets. *Astrophysical J.* **708**, 1326 (2009).

50. Zeng, L. & Sasselov, D. The effect of temperature evolution on the interior structure of h2o-rich planets. *Astrophysical J.* **784**, 96 (2014).

51. Kresse, G. & Furthmüller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. *Comput. Mat. Sci.* **6**, 15-50 (1996).

52. Kresse, G. & Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. *Phys. Rev. B* **54**, 11169-11186 (1996).

53. Blöchl, P. E. Projector augmented-wave method. *Phys. Rev. B* **50**, 17953 (1994).

54. Kresse, G. & Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. *Phys. Rev. B* **59**, 1758 (1999).

55. Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradient approximation made simple. *Phys. Rev. Lett.* **77**, 3865-3868 (1996).

56. Monkhorst, H. J. & Pack, J. D. Special points for brillouin-zone integrations. *Phys. Rev. B* **13**, 5188 (1976).

57. Nose, S. A unified formulation of the constant temperature molecular dynamics methods. *J. Chem. Phys.* **81**, 511 (1984).

58. Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. *Phys. Rev. A* **31**, 1695 (1985).

59. Parrinello, M. & Rahman, A. Polymorphic transitions in single crystals: a new molecular dynamics method. *J. Appl. Phys.* **52**, 7182-7190 (1981).

60. Lin, S.-T., Blanco, M. & G, W. A. III The two-phase model for calculating thermodynamic properties of liquids from molecular dynamics: validation for the phase diagram of lennard-jones fluids. *J. Chem. Phys.* **119**, 11792-11805 (2003).

61. Lai, P.-K., Hsieh, C.-M. & Lin, S.-T. Rapid determination of entropy and free energy of mixtures from molecular dynamics simulations with the two-phase thermodynamic model. *Phys. Chem. Chem. Phys.* **14**, 15206-15213 (2012).

62. Desjarlais, M. P. First-principles calculation of entropy for liquid metals. *Phys. Rev. E* **88**, 062145 (2013).

63. Grau-Crespo, R., Hamad, S., Catlow, C. & De Leeuw, N. Symmetry-adapted configurational modelling of fractional site occupancy in solids. *J. Phys.: Condens. Matter* **19**, 256201 (2007).

64. Pruzan, P. et al. Phase diagram of ice in the vii-viii-x domain: vibrational and structural data for strongly compressed ice viii. *J. Raman Spectroscopy* **34** (2003).

65. Bronstein, Y., Depondt, P., Finocchi, F. & Saitta, A. M. Quantum-driven phase transition in ice described via an efficient langevin approach. *Phys. Rev. B* **89**, 214101 (2014).

66. Schwager, B. & Boehler, R. H₂o: another ice phase and its melting curve. *High. Press. Res.* **28**, 431-433 (2008).

67. Madigan, M. T. et al. *Brock biology of microorganisms*, vol. 11 (Prentice hall Upper Saddle River, NJ, 1997).
## Acknowledgements

We acknowledge access to the GENCI supercomputing centers (CINES, IDRIS, CCRT) through a series of eDARI grants for performing the simulations (stl2876). Additional computational resources were provided by the Norwegian infrastructure for high-performance computing (NOTUR grants NN9329K, NN2916K, and NN9697K). R.C. was funded by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement n°681818 - IMPACT) and by the Deep Carbon Observatory under a grant from the Extreme Physics and Chemistry Directorate. J.A.H. and R.C. were funded by the Research Council of Norway through its Centres of Excellence funding scheme, project number 223272. SL was funded by a grant from the Agence Nationale de la Recherche (ANR) under project number ANR-15-CE31-0018-01.

## Author contributions

J.-A.H., R.C. and S.L. conceived the study and wrote the manuscript. J.-A.H. performed the atomistic simulations and analysis, J.-A.H and S.L. performed the geodynamic calculations.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-022-30796-5.

**Correspondence** and requests for materials should be addressed to Jean-Alexis Hernandez.

**Peer review information** *Nature Communications* thanks Sylvain Petitgirard and the anonymous reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022