# The effect of out-of-equilibrium outgassing on the cooling of planetary magma oceans

Arthur Walbecq, Henri Samuel, Angela Limare

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

## To cite this version:

Arthur Walbecq, Henri Samuel, Angela Limare. The effect of out-of-equilibrium outgassing on the cooling of planetary magma oceans. Icarus, 2025, 434, 10.1016/j.icarus.2025.116513. insu-05080884

HAL Id: insu-05080884

https://insu.hal.science/insu-05080884v1

Submitted on 23 May 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Research Paper

# The effect of out-of-equilibrium outgassing on the cooling of planetary magma oceans

Arthur Walbecq* , Henri Samuel, Angela Limare

Université Paris Cité, Institut de Physique du Globe de Paris, CNRS, 1 rue Jussieu, Paris, 75005, France

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

## ARTICLE INFO

Keywords:

Atmospheres

Evolution

Experimental techniques

Planetary dynamics

Terrestrial planets

Thermal histories

## ABSTRACT

Rocky planet mantles likely experienced at least one global magma ocean stage. These magma oceans cooled and solidified rapidly, setting the conditions under which rocky planet mantles evolved in the long-term. Upon magma ocean solidification, vigorous convective motions are commonly thought to efficiently outgas dissolved volatiles, progressively forming a secondary atmosphere. Through the atmospheric blanketing effect, exsolved volatiles can then significantly slow down the solidification of magma oceans, which can alter the final thermo-chemical state of the mantles and their long-term evolution. While the efficient outgassing is a common hypothesis in magma ocean evolution models, the outgassing efficiency may be limited by the fact that fluid parcels containing dissolved volatiles need to reach small pressures corresponding to shallow exsolution depths, $d_{exsol}$, allowing bubbles to form and volatiles to be outgassed. A recent numerical study (Salvador and Samuel, 2023) revealed that small exsolution depths can delay exsolution of volatiles beyond the solidification time of the mantle, in the case of a magma ocean of constant extent and for a constant exsolution depth.

Here, we extended this work by performing computational fluid dynamics and analog experiments at various Rayleigh and Prandtl numbers, which govern convective motions. We derived for the first time a general law that describes the flux of exsolved volatiles for a magma ocean of evolving thickness, $D$, and exsolution depth, $d_{exsol}$. For sufficiently high Reynolds numbers, in a magma ocean, convecting with velocities $v$, the time evolution of the fraction of exsolved volatiles $\chi$ obeys a first order differential equation $\dot{\chi} = v C_b d_{exsol}(1 - \chi)/D^2$, with $C_b = 6.6$.

We implemented this parameterized flux into a coupled magma ocean-atmosphere evolution model to more realistically assess the effect of convective motions on the formation of secondary atmospheres, along with their mantle evolution.

We found that for a broad range of parameters (e.g., high rotation rates, large planetary masses, or relatively low initial volatile contents) inefficient outgassing can lead to a mantle solidification more than one order of magnitude faster than in the case with an atmosphere in equilibrium with the magma ocean. This may impact not only the thermal state but also the composition of the solidified mantle in major and trace elements, and the subsequent billions of years of evolution of terrestrial planets.

## 1. Introduction

During the early stages of planetary evolution, heat sources provided by giant impacts (Tonks and Melosh, 1993; Canup, 2008), and the conversion of potential energy to heat during metal-silicate differentiation and core formation (Senshu et al., 2002; Monteux et al., 2009; Samuel et al., 2010; Rubie et al., 2015) are large enough to bring planetary mantles to temperatures up to thousands of Kelvins, implying that rocky planets likely experienced one or several partial or total mantle melting episodes during their history, in the form of magma

pools or magma oceans. The present-day mantles of rocky planets then result from the cooling of these magma oceans.

Magma ocean cooling dynamics can be reasonably approximated by a Rayleigh-Bénard convection problem, in which a volume of fluid is both cooled from above and heated from below. The dimensionless numbers that govern this problem are the Rayleigh number, that compares buoyancy forces to dissipating effects:

$$Ra = \frac{\alpha g \Delta T D^3}{\kappa v}, \qquad (1)$$

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: walbecq@ipgp.fr (A. Walbecq).</footnote>
and the Prandtl number, that compares viscous to thermal dissipation:

$$ Pr = \frac{\nu}{\kappa}, \qquad (2) $$

where $\alpha$ is the thermal expansivity, $g$ is the acceleration of gravity, $\Delta T$ is the superadiabatic temperature difference between the top and the bottom surfaces, $D = R_p - R_b$ is the height of convecting fluid, with $R_p$ the surface radius of the magma ocean and $R_b$ its bottom radius, $\kappa$ is the thermal diffusivity, and $\nu$ is the kinematic viscosity. In this context, the output quantities that characterize respectively the heat transfer and the flow regime are the dimensionless averaged heat flux at the boundaries and the dimensionless average velocity. These are represented by the Nusselt number $Nu = Q/(\lambda\Delta T/D)$ and by the Reynolds number $Re = UD/\nu$, where $Q$ is the mean heat flux, measured at the top or at the bottom surface of the domain, $\lambda$ is the thermal conductivity and $U$ is a representative velocity for the fluid. Depending on the definition, $U$ can refer to a large-scale speed, to a vertical or horizontal averaged speed, or to the Root Mean Squared (RMS) speed. For a global terrestrial magma ocean, the Rayleigh number can reach values as high as $10^{31}$ for a temperature difference of $10^3$ K (Nikolaou et al., 2019). As a result of these high $Ra$ values, heat transport is efficient, leading to large surface heat flux and rapid early planetary cooling. Upon solidification, a rheological front that marks the transition between a liquid-like to a solid-like behavior progresses upwards, and volatile content in the melt increases due to the preferential partitioning of (incompatible) volatiles in the residual melt (Fig. 1). The exsolution of supersaturated volatiles leads to the formation of a secondary atmosphere. Due to greenhouse gases such as $CO_2$ and $H_2O$, the atmosphere plays an important role in controlling the surface temperature and the outgoing heat flux of the magma ocean. A dense atmosphere strongly reduces the outgoing heat flux and increases the solidification time of a terrestrial magma ocean from a few thousands years to up to a million years (Lebrun et al., 2013; Nikolaou et al., 2019). The interplay between solidification dynamics and atmospheric formation can lead to a variety of surface conditions on rocky planets (Hamano et al., 2013; Salvador et al., 2017; Miyazaki and Korenaga, 2022). Moreover, depending on the duration of magma ocean solidification, the mantle can result in various compositional structures once it is fully solidified (Maurice et al., 2017; Boukaré et al., 2018).

Parameterized models have been developed over decades to assess the time evolution of magma oceans in contact with an evolving atmosphere (Abe, 1997; Elkins-Tanton, 2008; Lebrun et al., 2013; Salvador et al., 2017; Nikolaou et al., 2019; Bower et al., 2019; Lichtenberg et al., 2021; Miyazaki and Korenaga, 2022; Massol et al., 2023). In these models, atmosphere and magma ocean evolution are coupled through the outgassing of dissolved volatiles and the blanketing effect of the atmosphere. To compute volatile outgassing, parameterized magma ocean-atmosphere models systematically assume that due to the high magnitude of the convective velocities, the magma ocean is always in chemical equilibrium with its atmosphere (Abe, 1997; Elkins-Tanton, 2008; Lebrun et al., 2013; Salvador et al., 2017; Nikolaou et al., 2019; Miyazaki and Korenaga, 2022; Massol et al., 2023). However, magma ocean outgassing first requires that oversaturated fluid parcels reach shallow depths (typically a few hundred meters), where pressure is sufficiently low to allow for bubble formation via exsolution (Fig. 1). Additionally, fluid parcels that have crossed the saturation depth can remain saturated because of the energy barrier needed to form bubbles. Depending on the oversaturation pressure, oversaturated fluid parcels could travel upward hundreds of meter before the nucleation rate is sufficient for bubbles to be effectively formed (Sparks et al., 1994). Once bubbles are formed, they need to travel sufficiently fast to reach the surface and to be effectively outgassed. These three mechanisms can potentially significantly reduce magma ocean outgassing. So far, the validity of the equilibrium outgassing hypothesis has only been tested once in numerical simulations (Salvador and Samuel, 2023). By

comparing the saturation (or equilibrium) pressure of a given volatile species and the hydrostatic pressure, one can estimate the minimum depth at which the saturated volatiles can be exsolved. For a terrestrial magma ocean, the saturation depth corresponds to a shallow layer of tens of meters. Using passive tracers, the authors estimated the minimum depth reached over time by fluid parcels in convection experiments at a fixed finite Prandtl number, and varying systematically the convective vigor ($Ra$). By extrapolating to the terrestrial magma ocean case, they concluded that the outgassing time can significantly exceed the time spent in a fully liquid magma ocean.

While these results show that the ability of convection to transport dissolved volatiles to low pressures is a strong prerequisite to magma ocean outgassing, out-of-equilibrium outgassing has so far not been considered in coupled magma ocean-atmosphere models. Massol et al. (2023) recently discussed these results, making however a misinterpretation that led them to suggest that out-of-equilibrium outgassing was negligible because low-density volatiles bubbles would rapidly reach the surface of the magma ocean. Indeed, they misconstrued that Salvador and Samuel (2023) focused on the ability of convective dynamics to bring magma ocean parcels to shallow pressures where dissolved volatiles can exsolve and form bubbles, not on the ascent of volatiles bubbles to the surface of the magma ocean, which is a subsequent and radically different process.

Even though Salvador and Samuel (2023) showed that the equilibrium outgassing hypothesis is not systematically a reliable assumption, the approach using the minimum depth reached only allows to draw conclusions on a magma ocean of constant thickness and fixed exsolution depth. Upon volatile exsolution, the saturation pressure is bound to decrease and with it, the exsolution depth. As a result, the outgassing efficiency should decrease over time, leading to longer outgassing times. However, as an atmosphere forms above the magma ocean, its blanketing effect extends the lifetime of the magma ocean in its liquid phase, leading to more complete outgassing. Therefore, a coupled model is needed to study the effect of out-of-equilibrium outgassing on magma ocean evolution. In addition, the study of Salvador and Samuel (2023) considered only a fixed value of the Prandtl number of one, indicating that the thermal and momentum diffusivities are comparable, while Prandtl number values both larger or lower than unity are also plausible for terrestrial magma oceans.

In this work, we follow an approach similar to the one described above by measuring the ability of convective motions to bring material to shallow depths. In contrast with Salvador and Samuel (2023) that relied exclusively on numerical modeling of convective motion at a fixed Prandtl number of one, we also performed laboratory experiments that validate our numerical simulations, and we systematically vary the Prandtl number in addition to the Rayleigh number. We derived a general law to describe the flux of outgassed volatiles as a function of the governing parameters (the Rayleigh and Prandtl numbers), and the exsolution depth. We implemented this scaling law in a parameterized magma ocean cooling model, coupled to its atmosphere, allowing for the first time to account for out-of-equilibrium outgassing in this context.

The paper is structured as follows. The next Section 2 describes laboratory and numerical experiments with rigid boundaries and numerical experiments with free-slip boundaries. Section 3 presents the results of these experiments focusing on the quantification of the outgassing efficiency. The following Section 4 describes our parameterized magma ocean model, along with the new outgassing scheme we developed. Section 5 discusses the results of the parameterized magma ocean model where we investigated the effects of out-of-equilibrium outgassing on magma ocean solidification style and solidification times. Section 6 discusses the applicability and the validity of our model prior to concluding remarks in Section 7.

## 2. Convection experiments

To model magma oceans dynamics, we carried out convection experiments, with both analog and computational fluid dynamics
<figure><img src="image_5.png" /><figcaption>Fig. 1. Schematic representation of convective outgassing and solidification in a magma ocean of evolving thickness $D = R_p - R_{b,t}$, in the case where solidification starts at the Core-Mantle Boundary. The magma ocean convects with a characteristic velocity $V_{RMS}$. When volatiles dissolved into fluid parcels cross the exsolution depth, bubbles can form, ascend, and outgas into the atmosphere. The magma ocean cools down with heat flux $F_{conv}$, radiated to space by the atmosphere $F_{atm}$. As the magma cools down, the melt fraction $\phi$ diminishes. When the melt fraction becomes smaller than the critical fraction $\phi_c \approx 0.4$, the magma behaves essentially as a solid.</figcaption></figure>

approaches. The first set of analog and numerical experiments have rigid conditions on all boundaries and serves as a mutual benchmark of both approaches, allowing us to understand fluid motions in our convective experiments. The second set of numerical experiments considers free-slip conditions at vertical walls and at the top boundary, which are more suitable for reproducing the convective dynamics of cooling magma oceans.

## 2.1. Benchmark of numerical simulations against laboratory experiments

Here, we will briefly describe how we successfully benchmarked numerical experiments carried out with a previously described numerical code by using controlled laboratory experiments and particle velocimetry techniques (Schanz et al., 2013, 2016). These results are described in greater detail in Walbecq et al. (2024).

Our numerical experiments are performed using StreamV, a 3D Finite-Volume numerical code that solves the Navier-Stokes equations in a Cartesian geometry, coupled to the energy equation, under the Boussinesq approximation (Samuel and Evonuk, 2010; Samuel, 2012). The code was benchmarked against various analytical and numerical solutions, including turbulent flows (Salvador and Samuel, 2023 and references therein). The set of governing equations is written below in dimensionless form:

$$ \nabla \cdot \mathbf{v} = 0, \qquad (3) $$

$$ \frac{1}{Pr_0} \frac{D\mathbf{v}}{Dt} = -\nabla p + \nabla \cdot [\eta (\nabla \mathbf{v} + \nabla \mathbf{v}^T)] + Ra_0 \alpha T \mathbf{e}_z, \qquad (4) $$

$$ \frac{DT}{Dt} = \nabla^2 T, \qquad (5) $$

where $\mathbf{v}$ is the velocity vector, $t$ is the time, $p$ is the dynamic pressure, $\eta$ is the dynamic viscosity, $\alpha$ is the thermal expansion coefficient, $T$ is the temperature. $D/Dt$ expresses the Lagrangian derivative. $\mathbf{e}_z$ is the vertical unit vector pointing upwards, and $Ra_0$ and $Pr_0$ are the Rayleigh and Prandtl number, computed at the top boundary. The aspect ratio

reproduces the physical dimensions of the tank used for the laboratory experiments. Following Walbecq et al. (2024) the viscosity and thermal expansivity have the same dependence with the temperature as the ones measured for the fluid used to conduct the analog experiments, and can be expressed as:

$$ \eta(T) = \exp[-\ln(C_{\eta})T], \qquad (6) $$

$$ \alpha(T) = 1 + C_{\alpha}T, \qquad (7) $$

where $C_{\eta}$ and $C_{\alpha}$ are determined using fits in the temperature range of each analog experiment. Reproducing the same aspect ratio and temperature dependencies of the governing parameters facilitates the comparison between analog and numerical experiments. Table 1 lists the values of the governing parameters for each numerical and laboratory experiment.

Laboratory experiments were performed in a plexiglas tank of dimensions 30 cm × 10 cm × 15 cm. The tank is heated from below and cooled from above by thermostated water baths. We measure the top and bottom temperature and wait until the thermal steady state is reached before starting the velocity measurements. The latter are performed by seeding the fluid homogeneously with nearly neutrally buoyant polyamid particles. Because of their small size (60 µm) and their nearly neutral buoyancy, they can be considered as passive tracers of the flow. We then follow the fluid motions by particle tracking using a fast blinking LED and Particle Tracking Velocimetry (Schanz et al., 2016). Once we obtain the particle trajectories, we interpolate the velocities onto a fine grid, with a regular spacing of 2.27 mm. We correct the boundary velocities to match rigid boundary conditions while constraining the divergence of the velocity field to zero following the approach detailed in Walbecq et al. (2024). This procedure yields long time series of the complete three-dimensional measured velocity field, well adapted to the study of advection problems.

Fig. 2 displays the velocity fields for analog and numerical experiments conducted at the same $Ra$ and $Pr$ values. Panel (a) corresponds
<caption>**Table 1**</caption>

Set of parameters used for each experiment and for its numerical counterpart. The laboratory and numerical experiments are described in Section 2.1. Each experiment is defined by the values of the Rayleigh and Prandtl numbers, Ra<sub>exp</sub> and Pr<sub>exp</sub>. T<sub>top</sub> and T<sub>bot</sub> are the temperature at the top and the bottom of the tank respectively. The Rayleigh and Prandtl numbers are set by the average temperature and the properties of the fluid used. Ra<sub>0</sub> and Pr<sub>0</sub> are the Rayleigh and Prandtl numbers at the top boundary, which are used for the numerical experiments (see Eq. (4)). Ra<sub>sim</sub> and Pr<sub>sim</sub> are the Rayleigh and Prandtl numbers computed at T = (T<sub>top</sub> + T<sub>bot</sub>)/2.

<table><thead><tr><th colspan="5">Laboratory Experiments</th><th colspan="6">Numerical Experiments</th></tr><tr><th>N*</th><th>T<sub>top</sub>(°C)</th><th>T<sub>bot</sub>(°C)</th><th>Ra<sub>exp</sub></th><th>Pr<sub>exp</sub></th><th>Ra<sub>0</sub></th><th>Pr<sub>0</sub></th><th>C<sub>η</sub></th><th>C<sub>α</sub></th><th>Ra<sub>sim</sub></th><th>Pr<sub>sim</sub></th></tr></thead><tbody><tr><td>1</td><td>16.9</td><td>38.2</td><td>1.43 · 10<sup>9</sup></td><td>11.6</td><td>9.04 · 10<sup>8</sup></td><td>15.89</td><td>1.78</td><td>0.32</td><td>1.40 · 10<sup>9</sup></td><td>11.9</td></tr><tr><td>2</td><td>16.2</td><td>29.4</td><td>0.73 · 10<sup>9</sup></td><td>13.3</td><td>5.42 · 10<sup>8</sup></td><td>16.23</td><td>1.46</td><td>0.21</td><td>0.72 · 10<sup>9</sup></td><td>13.4</td></tr><tr><td>3</td><td>43.7</td><td>49.7</td><td>0.78 · 10<sup>9</sup></td><td>7.6</td><td>7.13 · 10<sup>8</sup></td><td>8.00</td><td>1.13</td><td>0.16</td><td>0.82 · 10<sup>9</sup></td><td>7.5</td></tr><tr><td>4</td><td>21.4</td><td>47.5</td><td>2.29 · 10<sup>9</sup></td><td>9.8</td><td>1.36 · 10<sup>9</sup></td><td>13.86</td><td>1.89</td><td>0.40</td><td>2.24 · 10<sup>9</sup></td><td>10.1</td></tr><tr><td>5</td><td>30.8</td><td>56.8</td><td>3.10 · 10<sup>9</sup></td><td>8.0</td><td>1.99 · 10<sup>9</sup></td><td>10.69</td><td>1.75</td><td>0.42</td><td>3.18 · 10<sup>9</sup></td><td>8.1</td></tr><tr><td>6</td><td>10.7</td><td>15.5</td><td>0.17 · 10<sup>9</sup></td><td>17.9</td><td>1.52 · 10<sup>8</sup></td><td>19.32</td><td>1.15</td><td>0.10</td><td>0.17 · 10<sup>9</sup></td><td>18.0</td></tr><tr><td>7</td><td>16.5</td><td>59.6</td><td>4.29 · 10<sup>9</sup></td><td>9.0</td><td>1.79 · 10<sup>9</sup></td><td>16.10</td><td>3.02</td><td>0.69</td><td>4.19 · 10<sup>9</sup></td><td>9.3</td></tr><tr><td>8</td><td>16.4</td><td>38.9</td><td>2.94 · 10<sup>8</sup></td><td>50.8</td><td>1.66 · 10<sup>8</sup></td><td>73.6</td><td>2.02</td><td>0.50</td><td>2.95 · 10<sup>8</sup></td><td>51.8</td></tr><tr><td>9</td><td>16.9</td><td>47.8</td><td>5.04 · 10<sup>8</sup></td><td>44.0</td><td>2.35 · 10<sup>8</sup></td><td>72.3</td><td>2.53</td><td>0.63</td><td>4.94 · 10<sup>8</sup></td><td>45.5</td></tr><tr><td>10</td><td>15.9</td><td>29.7</td><td>1.43 · 10<sup>8</sup></td><td>59.2</td><td>0.99 · 10<sup>8</sup></td><td>74.7</td><td>1.57</td><td>0.32</td><td>1.44 · 10<sup>8</sup></td><td>59.6</td></tr><tr><td>11</td><td>43.5</td><td>49.7</td><td>1.77 · 10<sup>8</sup></td><td>30.6</td><td>1.59 · 10<sup>8</sup></td><td>32.8</td><td>1.16</td><td>0.07</td><td>1.77 · 10<sup>8</sup></td><td>30.6</td></tr><tr><td>12</td><td>10.9</td><td>56.3</td><td>7.80 · 10<sup>8</sup></td><td>42.5</td><td>2.47 · 10<sup>8</sup></td><td>89.0</td><td>3.88</td><td>1.06</td><td>7.43 · 10<sup>8</sup></td><td>45.2</td></tr><tr><td>13</td><td>8.7</td><td>15.5</td><td>0.39 · 10<sup>8</sup></td><td>85.3</td><td>0.32 · 10<sup>8</sup></td><td>95.9</td><td>1.23</td><td>0.21</td><td>0.39 · 10<sup>8</sup></td><td>86.5</td></tr><tr><td>14</td><td>11.3</td><td>60.9</td><td>9.46 · 10<sup>8</sup></td><td>39.7</td><td>2.76 · 10<sup>8</sup></td><td>87.7</td><td>4.31</td><td>0.14</td><td>9.01 · 10<sup>8</sup></td><td>42.2</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 2.** Velocity amplitude for (a) a laboratory and (b) a numerical experiment, with rigid boundary conditions. The Rayleigh number is equal to 1.7 $10^8$ for both experiments. The Prandtl number is equal to 17.9 and 18.0 for panels a and b respectively. For each panel, the slices are located along the planes X = 0 cm, X = 10 cm, X = 20 cm, Y = 5 cm, Y = 10 cm and Z = 0 cm. The vectors represent the velocity field in each plane, whose norms are proportional to the velocity magnitudes.</figcaption></figure>

to a laboratory experiment and panel (b) displays the numerical counterpart. The conclusions shown below stand throughout all the Ra-Pr parameter space we explored.

The laboratory experiment (Fig. 2a), displays a large-scale structure characterized by two convecting cells. The velocity amplitude is therefore highest in the middle of the tank and in the vicinity of the edges of the domain. In between, in convection cell centers, the velocity is considerably lower. The numerical counterpart (Fig. 2b) displays the same flow structure with two convecting cells and a similar velocity amplitude pattern. The average velocity amplitude is close for each case, the experimental one being slightly smaller than the numerical one (2.9 $10^{-3}$ vs. 3.0 $10^{-3}$ m/s, respectively). Overall, we found that all fourteen experiments (laboratory and numerical) presented similar velocity structures, and comparable average amplitude velocities. The laboratory experiments always yielded slightly lower velocity amplitude, but we did not find more than a 5% difference. We interpret this difference in amplitude as an effect of the recording frequencies for the analog experiments. Particles moving too fast relative to the recording frequency cannot be properly tracked and the velocimetry technique underestimates their velocity. The good agreement observed between numerical and analog experiments cross-validates both approaches.

Because an atmosphere lies on top of the magma ocean, the above experiments with rigid boundaries may not entirely reproduce the dynamics of the magma ocean whose surface is closer to free-slip

boundary conditions. In particular, one can expect velocities near the surface to be higher in the free-slip case than in the rigid case. Similarly, the velocity gradients near the surface can also be quite different between rigid and free-slip boundaries. These might play an important role in the advection of particles, and eventually in the exsolution dynamics. As in Salvador and Samuel (2023) we therefore conducted numerical experiments with free-slip boundaries except at the bottom as described in the following Section.

## 2.2. Free-slip numerical experiments

We performed experiments for Rayleigh numbers ranging from $10^6$ to $10^9$ and for Prandtl numbers ranging from 0.1 to 100. The entire dataset consisting of 14 experiments is listed in Table 2. The grid size for each case is uniform along each direction and was set as to obtain at least four points in the thermal boundary layers, resulting in grids ranging from (256 × 64 × 128) to (384 × 96 × 192) cells, depending on the Rayleigh and Prandtl numbers. We kept the same aspect ratio as in the previous experiments with rigid boundary conditions. Due to the presence of convective rolls, the chosen geometry can have some influence on the dynamics of thermal convection. The heat transport can vary slightly with increasing aspect ratios and even present a maximum (Wagner and Shishkina, 2013; Stevens et al., 2024). However,
<caption>**Table 2**</caption>

Set of parameters for each carried out numerical experiment. Ra, Pr, Nu and Re are respectively the Rayleigh, Prandtl, Nusselt, and Reynolds numbers and have been defined in Section 1. C_b is the exsolutio coefficient, determined for each experiment and defined in Section 3.1.

<table><thead><tr><th>Experiment</th><th>Ra</th><th>Pr</th><th>Re</th><th>Nu</th><th>C<sub>b</sub></th></tr></thead><tbody><tr><td>1</td><td>10<sup>6</sup></td><td>1</td><td>319.9</td><td>12.47</td><td>2.723</td></tr><tr><td>2</td><td>3 10<sup>6</sup></td><td>1</td><td>520.0</td><td>16.65</td><td>3.147</td></tr><tr><td>3</td><td>10<sup>7</sup></td><td>1</td><td>820.0</td><td>22.82</td><td>3.886</td></tr><tr><td>4</td><td>3 10<sup>7</sup></td><td>1</td><td>1315</td><td>30.35</td><td>4.472</td></tr><tr><td>5</td><td>10<sup>8</sup></td><td>1</td><td>2124</td><td>41.67</td><td>5.375</td></tr><tr><td>6</td><td>3 10<sup>8</sup></td><td>1</td><td>3530</td><td>56.05</td><td>5.770</td></tr><tr><td>7</td><td>10<sup>9</sup></td><td>1</td><td>5558</td><td>79.98</td><td>6.461</td></tr><tr><td>8</td><td>10<sup>6</sup></td><td>10</td><td>45.5</td><td>12.77</td><td>1.75</td></tr><tr><td>9</td><td>10<sup>6</sup></td><td>100</td><td>4.33</td><td>11.84</td><td>1.752</td></tr><tr><td>10</td><td>10<sup>8</sup></td><td>10</td><td>442.3</td><td>47.33</td><td>3.434</td></tr><tr><td>11</td><td>10<sup>8</sup></td><td>100</td><td>70.31</td><td>48.67</td><td>2.250</td></tr><tr><td>12</td><td>10<sup>6</sup></td><td>0.1</td><td>1563</td><td>9.26</td><td>3.543</td></tr><tr><td>13</td><td>10<sup>8</sup></td><td>0.1</td><td>9381</td><td>30.71</td><td>5.866</td></tr><tr><td>14</td><td>10<sup>9</sup></td><td>0.1</td><td>26322</td><td>56.58</td><td>6.612</td></tr></tbody></table>

because the aspect ratio is the same for all of our experiments, it unlikely affects the trends we observe and our results in general.

Fig. 3 displays the velocity and the temperature fields for experiments at the same Rayleigh number Ra = 10⁸ and Prandtl numbers Pr = 0.1, 1, and 10 (for the top, middle and bottom row respectively). Figs. 3a–b show that the flow fields are turbulent and display a chaotic pattern. There is no clear large-scale structure emerging, neither in the velocity nor in the temperature field. Because Pr is relatively low, the Reynolds number Re is high (Re = 9381). The velocity field is thus turbulent and the bulk is efficiently mixed, resulting in a quasi-homogeneous temperature field.

For a higher Prandtl number Pr = 1 (Figs. 3c–d), the flow tends to be more structured. The velocity amplitude is high near the X = 0 cm and X = 30 cm walls, as well as in the middle of the tank. Between these high velocity zones, the velocity amplitude drops. The flow is going downwards at the left and right boundaries, and upwards in the middle of the tank. This correlates with the temperature field (Fig. 3d), which is less homogeneous than in Fig. 3b and presents a large hot patch in the middle as well as smaller cold patches on the sides. Because the Prandtl number is higher, the flow is less turbulent, allowing the hot patch to be maintained in time. It rises, due to high positive buoyancy and the cold patches on the sides fall, due to negative buoyancy. Therefore, the interplay between the temperature field and a less turbulent velocity field leads to this flow structure, with two large convecting cells.

With an even higher Prandtl number Pr = 10, the structure described above is more pronounced. The two-cell structure is also present (Fig. 3e), with the same circulation pattern (flow going upwards in the middle, downwards on the sides). Namely, the temperature field presents a well-marked rising thermal plume in the middle, which is thinner than in the previous case. This transition from a convection regime with turbulent flow and diffuse thermal plumes to a regime dominated by sharp plumes with increasing Pr numbers has been previously described in the literature and is therefore expected (Verzicco and Camussi, 1999; Schmalzl et al., 2002). For a fixed Pr number, increasing the Rayleigh number leads to a transition from extended plumes to small-scale structures. The flow structures in our numerical experiments are thus compatible with previous results from the literature (Castaing et al., 1989; Ahlers et al., 2009; van der Poel et al., 2011; Chillà and Schumacher, 2012; Schumacher and Bodenschatz, 2012).

We computed the Nusselt and Reynolds numbers. Their dependence on the governing parameters Ra and Pr is shown in Fig. 4.

A double fit on Ra and Pr yields a power law Nu ∝ Ra<sup>0.27±0.024</sup> Pr<sup>0.05±0.03</sup>, close to Nu ∝ Ra<sup>2/7</sup> (Shraiman and Siggia, 1990; Grossmann and Lohse, 2000). While our data exhibit a slight dependence on Pr, it is expected that this dependence vanishes in the magma ocean, with Pr > 1 (Verzicco and Camussi, 1999). We

therefore performed a second fit, on Ra alone, that yields a power law Nu ∝ Ra<sup>0.26±0.04</sup>, displayed in Fig. 4. We also represent the scaling law frequently used in magma ocean evolution models, Nu = 0.089 Ra<sup>1/3</sup> (Lebrun et al., 2013; Solomatov, 2015; Nikolaou et al., 2019). All three curves (red, yellow and dashed black lines) present a fair agreement in the Ra range we study. Our data also present a fair agreement with the results from Shraiman and Siggia (1990). The slight difference in amplitude between our fit and scalings from Shraiman and Siggia (1990) can be explained by differences in boundary conditions and aspect ratio. While the differences in heat flux are small over the range of Rayleigh numbers displayed here, they could be important when extrapolated to the magma ocean case. The influence of the chosen power law on magma ocean cooling and solidification has been previously investigated in Nikolaou et al. (2019). They found that when using the 2/7 exponent instead of the 1/3 exponent, the total solidification time only increased by 25%. Indeed, in these parameterized cooling models, the atmospheric temperature adjusts to a higher heat flux (and conversely for a lower heat flux) with an increase in temperature (resp. with a decrease in temperature), which mitigates the effect of the choice of the power law. For simplicity and comparison purposes, we will therefore only use the 1/3 scaling law, corresponding to the so called 'soft turbulence' regime (Solomatov, 2015).

Next, we analyzed the dependence of the Reynolds number on our governing parameters Ra and Pr. Because we expect a transition between laminar and turbulent regimes, we limited the following analysis to relatively large Reynolds number (Re > 200). Because of this restriction, the Prandtl number ranges from 0.1 to 10. The limited range then makes it difficult to perform a fit on Pr. Following the scaling proposed by Shraiman and Siggia (1990), we therefore prescribed the Re dependency on Pr. We performed a fit on RePr<sup>5/7</sup> as a function of the Rayleigh number. This yields a scaling law RePr<sup>5/7</sup> = (1.09<sup>+0.45</sup><sub>-0.32</sub>)Ra<sup>0.41<sup>+0.02</sup><sub>-0.02</sub></sup>. This law is close to the one from Shraiman and Siggia (1990), that predicted an exponent equal to 3/7 ≈ 0.43. The scalings from Shraiman and Siggia (1990) include however a logarithmic correction that we cannot test due to our limited range on ln(Re). Instead, we compare our results with Solomatov (2015), who predicts Re = 0.6Ra<sup>3/7</sup>Pr<sup>-5/7</sup>. Fig. 4b shows the comparison of the two laws, as well as our experimental data. The two curves display a similar behavior. The 3/7 law stands below our best fit. This difference in amplitude can be explained by differences in aspect ratio and in boundary conditions. Differences in the velocity used in the definition of Re can also cause additional discrepancies. Fig. 4b also indicates that the correction Pr<sup>5/7</sup> on Re is a good choice, as it collapses the data points on a single curve. The only exception to this are data points at low Reynolds numbers (i.e., either low Ra or high Pr), which are not considered for this fit. The fact that these data points do not collapse on the curve RePr<sup>5/7</sup> ∝ Ra<sup>0.41</sup> indicates that their exclusion from the fit is a valid hypothesis.

The reasonable agreement between our results and those reported in the literature indicates that our findings are compatible with scalings commonly used in magma ocean cooling studies. In the following, we therefore resort to classical scalings that are appropriate for extrapolation to Ra and Pr corresponding to terrestrial magma oceans.

### 3. Estimation of the outgassing rate

Outgassing efficiency in the magma ocean does not depend on convective velocities only, but rather on the ability of the flow to transport material above a given shallow exsolutio depth. In the following, we consider the most optimistic scenario for outgassing i.e., once fluid particles reach the exsolutio depth, oversaturated volatiles are exsolved in the form of bubbles that are able to reach the surface and be outgassed instantly.

Our approach follows the one used in Salvador and Samuel (2023) that provides the typical outgassing time for a given exsolutio depth,
<figure><img src="image_7.png" /><figcaption>Fig. 3. Velocity amplitudes (left column) and temperature fields (right column) for three numerical experiments at statistical steady state, at $Ra = 10^8$ and $Pr = 0.1, 1$, and 10 (panels a–b, c–d and e–f respectively). The boundary conditions are set to free-slip at all boundaries except the bottom one, which is set to rigid. For each panel, the slices are located along the planes $X = 0$ cm, $X = 10$ cm, $X = 20$ cm, $Y = 5$ cm, $Y = 10$ cm and $Z = 0$ cm. The arrows represent the velocity field in each plane, the arrow length being proportional to the magnitude of the velocity field.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 4. Comparison of our results with scaling laws from the literature. (a) Nusselt number as a function of the Rayleigh number, in the form $Nu \propto Ra^\gamma$. Multiple y-values correspond to simulations with different $Pr$. The best fit to our data (blue dots) indicates an exponent $\gamma = 0.263$ (red line), close to the value $\gamma = 2/7$ from Shraiman and Siggia (1990) (yellow line). The dashed lines indicate scaling proposed in Eq. (A.11) (Solomatov, 2015). (b) Reynolds number as a function of the Rayleigh number. Multiple y-values correspond to simulations with different $Pr$. The red line is the best fit from our data, performed by setting the $Pr$ exponent to $-5/7$ and the yellow line correspond to $RePr^{5/7} = 0.6 Ra^{3/7}$ (Eq. (19), Solomatov, 2015). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

but has a limited applicability for the description of the time evolution of a solidifying magma ocean. An alternative to overcome this limitation is a parameterized description of the exsolution rate, rather than the exsolution time scale. Therefore, instead of tracking the minimum depth ever reached by fluid parcels, we adopted a more direct approach

where we considered the exsolution depth as a variable and estimate the amount of fluid parcels having reached this depth over time.

Each numerical experiment is seeded with $N_{tot}$ passive tracers that are then advected. Each tracer $i$ is associated with a counter $c_i$, initialized to 0. If over time, a tracer travels above a given exsolution
<figure><img src="image_9.png" /><figcaption>Fig. 5. Evolution of the exsolved fraction $\chi$ as a function of the number of transits $N_t$. (a) Evolution for a single experiment at $Ra = 3 \cdot 10^8$ and $Pr = 1$ and various exsolution depths. (b) Evolution for a fixed exsolution depth and various numerical experiments at $Pr = 1$.</figcaption></figure>

depth, the counter is set to 1. We considered 20 exsolution depths for each experiment, ranging from $d_{exsol,min}$ to $d_{exsol,max}$. The minimum exsolution depth $d_{exsol,min}$ is set independently for each experiment depending on the grid resolution. We performed tests for different tracer densities, ranging from 1 to 128 tracers per cell and found no significant difference in the results. For safety, we used a density of 8 tracers per cell for all experiments. By dividing the sum of the counters by the total number of tracers, one can approximate the number of particles having reached the surface over time. The fraction $\chi$ of exsolved particles associated with the exsolution depth $d_{exsol}$ at time $t$ can then be estimated as:

$$ \chi(d_{exsol}, t) = \frac{1}{N_{tot}} \sum_{i=1}^{N=N_{tot}} c_i(d_{exsol}, t). \qquad (8) $$

We expect outgassing efficiency to scale with convective velocities. Following Salvador and Samuel (2023), we introduce the number of transits:

$$ N_t = \int_0^t \frac{V_{RMS}}{D} dt, \qquad (9) $$

where $V_{RMS}$ is the RMS of the velocities and $D$ is the distance between the top and bottom boundaries. This dimensionless time allows one to compare the time spent in a given thermal state to the time taken for a fluid particle traveling at velocity $V_{RMS}$ to move from one horizontal surface to the other. Using an integral rather than the time ratio allows one to estimate the number of transits in situations where the velocity or the depth change over time.

Fig. 5a shows the time evolution of the exsolved fraction for a single experiment, at $Ra = 3 \cdot 10^8$ and $Pr = 1$ and compares various exsolution depths. Fig. 5b shows the time evolution of the exsolved fraction for various Rayleigh numbers, ranging from $Ra = 10^6$ to $Ra = 10^9$ and for a fixed exsolution depth $d_{exsol} = 7 \cdot 10^{-3}$ and a fixed Prandtl number $Pr = 1$. We observe that the curves computed with the smallest exsolution depths increase more slowly and the curves computed with the largest exsolution depths increase more rapidly. The exsolution time therefore appears to be inversely proportional to the exsolution depth. The yellow and orange curves (for $d_{exsol} = 4.3 \cdot 10^{-2}$ and $d_{exsol} = 1.9 \cdot 10^{-2}$) converge towards the asymptotic value $\chi = 1$ after 30 and 60 transits, respectively, consistent with the results reported in Salvador and Samuel (2023). For smaller exsolution depths, the exsolution time is too large to be directly observed in absence of very long time series of the velocity field which require prohibitively long computing time and large data storage.

The experiments considered in Fig. 5b display a wide range of velocity amplitudes, with Reynolds numbers ranging from $Re = 319.9$ (for $Ra = 10^6$, $Pr = 1$) to $Re = 5558$ (for $Ra = 10^9$, $Pr = 1$). The

exsolution curves display however very similar dynamics, with close characteristics exsolution times. This observation shows that while the number of transits is the adequate variable to describe exsolution dynamics, additional work is needed to address small differences that may remain between distinct experiments. These differences are addressed in the following section.

### 3.1. Characterization of the volatile exsolution efficiency

To extrapolate the observed exsolution dynamics to the magma ocean case, we need to parameterize each exsolution curve as a function of the Rayleigh and Prandtl numbers and as a function of the exsolution depth. By plotting the derivative of the exsolved fraction with respect to the number of transits against the exsolved fraction (Fig. 6a), we have observed that regardless of $Ra$, $Pr$ and $d_{exsol}$, the form of the exsolution remains the same and can be described by a first-order differential equation:

$$ \frac{d\chi}{dN_t} = A\chi + B, \qquad (10) $$

where the coefficients $A$ and $B$ depend on $Ra$, $Pr$ and on $d_{exsol}$. Because all experiments are expected to reach $\chi = 1$ when $N_t$ tends to infinity, all exsolution curves obey $A = -B$, indicating that the parameter $B$ is the inverse of the exsolution time. Our experiments show that for $d_{exsol} < 0.01$, the parameter $B$ depends linearly on the exsolution depth (Fig. 6b). For global planetary magma oceans, the exsolution depths are small compared to the magma ocean depth. We therefore limit our analysis to small exsolution depths and we express parameter $B$ as $B = C_b d_{exsol}$. We subsequently refer to $C_b$ as the exsolution coefficient. By construction, the exsolution coefficient depends on $Ra$ and $Pr$ only. Fig. 7a shows exsolution curves reconstructed by solving Eq. (10), using the fitted exsolution coefficient from our experiment and Fig. 7b shows the exsolution coefficient as a function of $Re$.

All curves in Fig. 7 exhibit good agreement between the raw data from the advection of tracers and reconstructed exsolution curves. The agreement is better for small exsolution depths and the reconstruction tends to slightly overestimate the exsolved fraction for larger exsolution depths ($d_{exsol} = 8.5 \cdot 10^{-3}$). Because we limit the fit to small values of $d_{exsol}$, relevant to the magma ocean case, the exsolution coefficient does not accurately describe exsolution curves for large $d_{exsol}$. The misfit for $d_{exsol} = 8.5 \cdot 10^{-3}$ is hence expected. For $d_{exsol} < 0.01$, the exsolution dynamics is therefore well described by a first order differential equation where the exsolution coefficient is the only governing parameter. Fig. 7b shows that the exsolution coefficient varies between experiments, for varying Reynolds numbers. The exsolution coefficient
<figure><img src="image_10.png" /><figcaption>**Fig. 6.** (a) $\frac{d\chi}{dN}$ plotted against $\chi$ for various exsolution depths, for experiment 6 in Table 2. All curves exhibit a linear behavior. For the lowest exsolution depth (red curve), extrapolation of the fit indicates that the curve tends to $(\frac{d\chi}{dN}, \chi) = (0, 1)$. (b) Evolution of the fitting parameter *B* with the exsolution depth. The linear fit is computed for $d_{\text{exsol}} < 0.01$ D. Above this value, the parameter *B* is lower than the predicted parameter.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 7.** Left panel shows the exsolved fraction of experiment 6 in Table 2 ($Ra = 3 \times 10^8$ and $Pr = 1$) for various exsolution depths using an exponential fit; blue stars represent the data and red line represents the fit. The right panel shows the evolution of the exsolution coefficient $C_b$ with the Reynolds number $Re$. Each dot is color-coded according to the corresponding Prandtl number. Error bars on $Re$ correspond to time-averaging uncertainties and error bars on $C_b$ correspond to uncertainties on the fitted model. Both correspond to a 95% confidence interval. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

$C_b$ first increases as a function of the Reynolds number and seems to saturate at high Reynolds numbers, above $Re = 5 \times 10^3$, possibly less. The existence of such plateau may be explained by the fact that for sufficiently high Reynolds numbers, the flow is fully turbulent and the interior is well mixed, resulting in efficient transport of fluid towards shallow depths. Moreover, one can expect the upward velocity to play an important role in transporting fluid particles to shallow depths. At low Reynolds numbers, in our configuration, the flow is composed of a large single cell, making the vertical velocity small, relative to the $V_{RMS}$. We could therefore expect more efficient outgassing for high Reynolds numbers, when the flow is more isotropic. This is further supported by the fact that $dC_b/dRe$ decreases for values greater than $Re \approx 5 \times 10^3$. Therefore in the following, we use the maximum exsolution coefficient observed $C_b = 6.6$ to estimate outgassing efficiency for magma ocean, where we expect the Reynolds number to be sufficiently high with turbulent fluid motions.

Our convective outgassing model, through a parameterized first order differential equation, allows for estimations of volatile fluxes

between an evolving magma ocean and its atmosphere. In the following Section we apply this description to parameterized magma ocean cooling model coupled to an evolving atmosphere.

## 4. Application to planetary magma ocean-atmosphere coupled evolution

This section describes our parameterized magma ocean cooling and solidification model. It consists of a 1D spherically symmetric magma ocean, coupled with its growing atmosphere. The main parameters defining the magma ocean are described in Fig. 1.

Our approach follows that of previous works (Elkins-Tanton, 2008; Lebrun et al., 2013; Salvador et al., 2017; Nikolaou et al., 2019), with however noticeable differences in the treatment of non-equilibrium outgassing that is absent in all previous studies, and with distinct assumptions regarding the differential cooling of the magma ocean above and below the rheological front. For completeness, we provide a full description of the governing equations for this model in Appendix A. By comparing the temperature profile $T(r)$ with experimentally measured melting curves, one can compute the melt fraction $\phi$ in the mantle along
the planetary radius *r*. The melt fraction $\phi(r)$ defines the structure of the mantle between essentially (or fully) molten and essentially solidified zones. This step allows to compute the partitioning of volatiles elements between the solid mantle, the melt and the growing atmosphere. Once atmospheric pressures are known, the outgoing heat flux, $F_{\text{conv}}$, can be computed (see appendix for the detailed expressions) and the conservation of energy that describe the thermal evolution of the silicate mantle writes:

$$ \frac{d}{dt} \int_{R_0(t)}^{R_p} \rho T \left( C_p + \Delta H \frac{d\phi}{dT} \right) r^2 dr = - R_p^2 F_{\text{conv}}, \quad (11) $$

where $T(r)$ is the temperature profile along the planetary radius, $t$ is the time, $R_p$ is the radius of surface of the magma ocean, $\rho$ is the magma density, $C_p$ is the heat capacity, $\Delta H$ is the specific enthalpy difference for fusion/solidification. The radius $R_0(t)$ defines the time-dependent volume over which the conservation of heat is integrated. In previous studies, this radius was set to $R_0(t) = R_c$, with $R_c$ the radius of the core-mantle boundary, implying the heat surface loss is distributed over the entire mantle, regardless of it being liquid or not (Lebrun et al., 2013; Salvador et al., 2017; Nikolaou et al., 2019). However, due to the high viscosity of the mantle at relatively low melt fractions, convective motions in the mantle are slow and thus inefficient at transporting heat. Therefore, we instead integrate Eq. (11) only in the essentially liquid part, with $R_0(t) = R_{b,t}$. This represents a major difference compared with previous models, which considerably influences the solidification time. The detailed justification for our choice can be found in Appendix A and Appendix C. The above energy balance can be expressed as a function of the potential temperature, $T_p$ (see Appendix A). Once $T_p$ is updated, we can compute the new temperature profile. This loop is then repeated until the rheological front reaches the surface.

## 4.1. Volatile outgassing

We describe atmospheric formation coupled to magma ocean cooling. We consider two styles of atmospheric formation. The first one is equilibrium outgassing commonly assumed in Elkins-Tanton (2008), Lebrun et al. (2013) or in Nikolaou et al. (2019). The second one stems from our exsolution model and is based on computation of the exsolution depths at each time step. We consider two volatiles species for the atmosphere: CO₂ and H₂O but other species could be considered as well (Lichtenberg et al., 2021). To primarily facilitate the comparison with previous studies (Elkins-Tanton, 2008; Lebrun et al., 2013; Salvador et al., 2017; Nikolaou et al., 2019) we do not consider the effect of volatile coupling reported in Bower et al. (2019) or in Massol et al. (2023), which represents better the atmospheric composition and the effect of partial pressures in volatiles. However, we evaluated the influence of such a simplification in the discussion Section, which we found to be small for the cases we considered.

### 4.1.1. Equilibrium outgassing

Because volatiles are more soluble in the melt than in the solid phase, the volatile content in the melt increases over time upon magma ocean cooling and solidification. In the frame of equilibrium outgassing, the outgassing of oversaturated volatiles is assumed to be instantaneous. The atmospheric pressure of a given volatile species is then given by its saturation pressure. We compute the saturation pressure through solubility laws from Holloway et al. (1994) for CO₂:

$$ P_{\text{sat,CO}_2} = \frac{X_{\text{CO}_2}}{4.4 \ 10^{-12}}, \quad (12) $$

and from Pan et al. (1991) for H₂O:

$$ P_{\text{sat,H}_2\text{O}} = \left( \frac{X_{\text{H}_2\text{O}}}{6.8 \ 10^{-8}} \right)^{1/0.7}, \quad (13) $$

where X refers to the mass concentration of a volatile species (CO₂ or H₂O). While more recent solubility measurements exist (Sossi et al.,

2023), we chose to use the same solubility laws as in previous parameterized cooling models (Lebrun et al., 2013; Salvador et al., 2017; Nikolaou et al., 2019), which facilitates comparison. The effect of using different solubility laws is discussed in Section 6.

By introducing the partition coefficient $\kappa$ between melt and solid phases for each volatile species, we can express the mass balance for each volatile:

$$ X_0 M_0 = M_{\text{atm,t}} + M_{\text{s,pv,t}} \kappa_{\text{pv}} X_t + M_{\text{s,lhz,t}} \kappa_{\text{lhz}} X_t + M_{\text{l,t}} X_t, \quad (14) $$

where $X_0$ is the initial volatile content, $M_0$ is the initial mass of the magma ocean, $M_{\text{atm,t}}$ is the atmospheric mass, $M_{\text{s,pv,t}}$ and $M_{\text{s,lhz,t}}$ are respectively perovskite and lherzolite masses, $\kappa_{\text{pv}}$ and $\kappa_{\text{lhz}}$ their partition coefficients, $M_{\text{l,t}}$ the liquid mass and $X_t$ the volatile content in the magma ocean. The subscripts $t$ indicate that the quantity is evaluated at time step $t$. The atmospheric mass is given by:

$$ M_{\text{atm}} = P_{\text{sat}}(X) \frac{4\pi R_p^2}{g}. \quad (15) $$

For a given thermal state of the magma ocean, Eq. (14) depends only on the equilibrium volatile content $X_t$ and can be solved analytically or numerically, depending on the expression of $P_{\text{sat}}(X)$. We chose to solve Eqs. (14) and (15) together with Eq. (12) or (13) for each volatile using a Newton solver. In the case of equilibrium outgassing, the volatile content in each reservoir is thus entirely determined by the thermal state of the mantle.

### 4.1.2. Out-of-equilibrium outgassing

For an out-of-equilibrium outgassing case, magma ocean exsolution is limited by the large number of overturns needed to transport fluid parcels to shallow outgassing depths. This results in a magma ocean that is not in equilibrium with its atmosphere. To account for this, we still assume that the liquid and solid mantle are in equilibrium with each other but we consider that the atmosphere is not in equilibrium with the liquid mantle. Because of the high viscosity contrast between the solid and liquid mantle, the assumption that these phases are chemically equilibrated is a strong hypothesis. However, due to the small values of the partition coefficient (see Table 3), partitioning of volatiles between the solid and liquid mantle phases only plays a marginal role on outgassing and on the thermal evolution of the magma ocean. Dropping this hypothesis and setting the partition coefficients to zero in Eq. (14) would lead to slightly higher final atmospheric pressures, but do not change the picture we describe below. An alternative scenario for mantle-magma chemical equilibration would be to consider that due to inefficient compaction at the rheological front, a non-negligible fraction of volatiles dissolved in the (partial) melt can remain trapped in the mantle (Hier-Majumder and Hirschmann, 2017; Sim et al., 2024). This scenario represents another end-member, where not all the initial volatile content is available for outgassing. In that sense, it is similar to the basal magma ocean scenario that we discuss in Section 5.3.

For each volatile species, the atmospheric mass is then smaller than the equilibrium mass computed with Eq. (14). We define the difference between these two masses as the potentially exsolvable mass $M_{\text{exsol,t+\delta t}} = M_{\text{atm,eq,t+\delta t}} - M_{\text{atm,t}}$. We consider that out-of-equilibrium volatiles species can only be outgassed at shallow depths, where the volatile saturation pressure is higher than the hydrostatic pressure. This leads to the following expression for the saturation depth:

$$ d_{\text{sat}} = \frac{P_{\text{sat}}(X) - P_{\text{atm,t}}}{\rho g}, \quad (16) $$

where $P_{\text{sat}}(X)$ is the saturation pressure, $P_{\text{atm,t}}$ is the atmospheric pressure and $\rho$ is the melt density. Here, we consider that upon crossing the saturation depth, all saturated volatiles form bubbles instantaneously. Under this conservative hypothesis, the exsolution depth is equal to the saturation depth for each volatile. Contrary to what was considered in the previous sections, the exsolution depth is therefore no longer a
<caption>Table 3</caption>

Parameters values used in the parameterized cooling model. The last column refers to the first equation in which these parameters are defined and used. Several values are used for the absorption coefficient of CO₂, which are defined in Section 5.3. The Earth's radius, core-mantle boundary and surface gravity are used as references for a planet of one Earth mass and are used in the set of Eqs. (22).

<table><thead><tr><th>Parameter</th><th>Description</th><th>Value</th><th>Definition</th></tr></thead><tbody><tr><td>α<sub>0</sub></td><td>Surface thermal expansivity</td><td>10<sup>-4</sup> K<sup>-1</sup></td><td>Eq. (A.2)</td></tr><tr><td>K<sub>0</sub></td><td>Surface bulk modulus</td><td>200 10<sup>9</sup> GPa</td><td>Eq. (A.2)</td></tr><tr><td>K'</td><td>Pressure derivative of bulk modulus</td><td>4</td><td>Eq. (A.2)</td></tr><tr><td>A<sub>k</sub></td><td>Viscosity parameterization coefficient</td><td>2.4 10<sup>-4</sup> Pa s</td><td>Eq. (A.9)</td></tr><tr><td>B<sub>k</sub></td><td>Viscosity parameterization coefficient</td><td>4600 K</td><td>Eq. (A.9)</td></tr><tr><td>C<sub>k</sub></td><td>Viscosity parameterization coefficient</td><td>1000 K</td><td>Eq. (A.9)</td></tr><tr><td>φ<sub>c</sub></td><td>Critical melt fraction</td><td>0.4</td><td>Eq. (A.10)</td></tr><tr><td>κ<sub>pu,CO<sub>2</sub></sub></td><td>Partition coefficient</td><td>5 10<sup>-4</sup></td><td>Eq. (14)</td></tr><tr><td>κ<sub>ihz,CO<sub>2</sub></sub></td><td>Partition coefficient</td><td>2.1 10<sup>-3</sup></td><td>Eq. (14)</td></tr><tr><td>κ<sub>pu,H<sub>2</sub>O</sub></td><td>Partition coefficient</td><td>1 10<sup>-4</sup></td><td>Eq. (14)</td></tr><tr><td>κ<sub>ihz,H<sub>2</sub>O</sub></td><td>Partition coefficient</td><td>1.1 10<sup>-2</sup></td><td>Eq. (14)</td></tr><tr><td>λ</td><td>Thermal conductivity</td><td>4 W m<sup>-1</sup> K<sup>-1</sup></td><td>Eq. (A.11)</td></tr><tr><td>A</td><td>Planetary albedo</td><td>0.3</td><td>Eq. (A.13)</td></tr><tr><td>S<sub>0</sub></td><td>Incident stellar flux</td><td>1361 W m<sup>-2</sup></td><td>Eq. (A.13)</td></tr><tr><td>σ</td><td>Stefan-Boltzmann constant</td><td>5.67 10<sup>-8</sup> W m<sup>-2</sup> K<sup>-4</sup></td><td>Eq. (A.13)</td></tr><tr><td>k<sub>0,CO<sub>2</sub></sub></td><td>Absorption coefficient for CO<sub>2</sub></td><td>[1 · 10<sup>-4</sup> – 5 · 10<sup>-2</sup>] m<sup>2</sup> kg<sup>-1</sup></td><td>Eq. (A.15)</td></tr><tr><td>k<sub>0,H<sub>2</sub>O</sub></td><td>Absorption coefficient for H<sub>2</sub>O</td><td>1 · 10<sup>-2</sup> m<sup>2</sup> kg<sup>-1</sup></td><td>Eq. (A.15)</td></tr><tr><td>C<sub>p</sub></td><td>Heat capacity</td><td>1000 J kg<sup>-1</sup> K<sup>-1</sup></td><td>Eq. (A.22)</td></tr><tr><td>ΔH</td><td>Specific enthalpy difference</td><td>4 · 10<sup>5</sup> J kg<sup>-1</sup></td><td>Eq. (A.19)</td></tr><tr><td>Ω</td><td>Rotation rate</td><td>10<sup>-3</sup> rad s<sup>-1</sup></td><td>Eq. (20)</td></tr><tr><td>R<sub>p,0</sub></td><td>Earth radius</td><td>6371 km</td><td>–</td></tr><tr><td>R<sub>c,0</sub> = R<sub>b,t=0</sub></td><td>Core-Mantle boundary</td><td>3481 km</td><td>–</td></tr><tr><td>g<sub>0</sub></td><td>Surface gravity</td><td>9.81 m s<sup>-2</sup></td><td>–</td></tr></tbody></table>

free parameter and is imposed by the volatile content and the pressure gradient.

Solving Eq. (10) over an infinitesimal transit $δN$ then indicates the exsolved fraction:

$$ \chi(d_{\text{exsol}}, \delta N) = 1 - \exp(-B\delta N) \quad (17) $$

Rewriting Eq. (17) in its dimensional form, assuming the saturation depth equals the exsolution depth then yields the exsolved mass $δM$ between consecutive time steps $t$ and $t + δt$:

$$ \delta M = M_{\text{exsol},t} \chi(d_{\text{exsol}}, \delta N) = M_{\text{exsol},t} \left[ 1 - \exp \left( - \frac{V_{\text{RMS}} C_b d_{\text{sat}} \delta t}{D^2} \right) \right], \quad (18) $$

where $V_{\text{RMS}}$ is the mean flow velocity, $C_b$ is the exsolution coefficient described in Section 3, $D = R_p - R_{b,t}$ is the evolving magma ocean depth, with $R_{b,t}$ the bottom of the magma ocean at time $t$ and $δt$ is the time step size.

# 5. Results

## 5.1. Temporal evolution of the equilibrium outgassing case

In this Section, we compute the time evolution of the reference, equilibrium outgassing case that consists of a terrestrial magma ocean, starting with a potential temperature of $T_p = 4000$ K, with no rotation. We adopt two pairs of values for the initial volatile content of CO₂ and H₂O. The first pair, later designated as a high volatile content, corresponds to initial volatile content of 410 and 130 ppm for H₂O and CO₂ respectively. These values correspond to the one adopted in Elkins-Tanton (2008), Lebrun et al. (2013) and Nikolaou et al. (2019). The second pair, later designated as a low volatile content, corresponds to initial volatile content of 41 and 13 ppm for H₂O and CO₂ respectively (Nikolaou et al., 2019). For these two cases, the emissivity is computed using Eqs. (A.14) to (A.17). The absorption coefficient for H₂O is set to $k_{0,H_2O} = 10^{-2}$ m²/kg and the absorption coefficient for CO₂ is set to $k_{0,CO_2} = 10^{-3}$ m²/kg as in Nikolaou et al. (2019). For these reference cases, we considered equilibrium outgassing of the magma ocean. We also performed computations for a black body radiation scenario. In this latter scenario, the emissivity is set to $\epsilon = 1$

throughout the entire mantle evolution. Fig. 8 shows the time evolution of atmospheric pressures, emissivities, temperatures, heat fluxes, and magma ocean depth. The blue curves refer to the high volatile end-member and the red curves refer to the low volatile end-member. The black body case is shown in black, when relevant.

Fig. 8a displays the evolution of the atmospheric pressures. The high and low volatile cases show similar trends, except shifted in time. In both cases the pressures are constant for the first thousand years (7303 and 922 years for high and low volatile content respectively. While some crystals form during this first phase, the melt fraction remains above $\phi_c = 0.4$, therefore the entire mantle remains essentially liquid and all terms in Eq. (14) remain comparable. After this first essentially liquid phase, the melt fraction passes below $\phi_c = 0.4$. Because of the small values of the melt-solid partition coefficients for volatile species, the creation of a solid phase shifts the equilibrium towards higher atmospheric masses. The atmospheric pressures then gradually increase, following the evolution of the mantle. Despite a higher initial content for H₂O, both cases show higher atmospheric pressures in CO₂ in the first phase. Because H₂O has a higher melt solubility than CO₂, it remains dissolved while the mantle is mostly liquid. This results in an atmosphere first dominated by CO₂. However, the atmosphere changes to being dominated by H₂O when the mantle becomes mostly solid. At this moment, the right hand side of Eq. (14) is dominated by the atmospheric term. The atmosphere then reflects the initial volatile content.

The atmospheric evolution directly affects the emissivity (Fig. 8b). It is constant at first and drops once the atmospheric pressures start to increase. For lower volatile contents, higher emissivities are expected (see Eq. (A.17) in the appendix), consistent with what we observed. The high volatile content end-member contains ten times more volatiles than in the low volatile end-member. This results in a more opaque atmosphere, and in emissivities that are approximately ten times lower.

Fig. 8c shows that despite being shifted in time, the potential temperatures follow the same evolution for all three cases (high volatile content, low volatile content and black body). The potential temperature gradually decreases until the rheological front reaches the surface. The three cases exhibit considerably different surface temperatures. The lowest surface temperature is that of the black body case, the next colder one is that of the low volatile content end-member and the hottest is the one of the high volatile content end-member. In
<figure><img src="image_12.png" /><figcaption>**Fig. 8.** Time evolution for the equilibrium outgassing case (a) the atmospheric pressures, (b) the emissivity, (c) the temperature, (d) the heat flux and (e) the magma ocean depth. The blue lines refer to the high volatile content case, for which initial content of CO₂ and H₂O are 130 ppm and 410 ppm, respectively. The red lines refer to the low volatile content case, for which initial content of CO₂ and H₂O are 13 ppm and 41 ppm respectively. In panels (c), (d) and (e), dark lines refer to a black-body case, where the emissivity has been set to 1 and the volatile content therefore plays no role. In panel (a), the solid lines show the atmospheric pressure of CO₂ and the dashed lines show that of H₂O. All curves were computed with $k_0 c_{O_2} = 10^{-3}$ m²/kg. In panel (c), the solid lines represent the potential temperature $T_p$ and the dashed lines represent the surface temperature $T_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>

our gray (or black) body model, the surface temperature can range anywhere from $T_{eq}$ to $T_p$ and the atmosphere adjusts its temperature immediately. If we consider an atmosphere with a high emissivity, the surface temperature will be close to $T_{eq}$ (see appendix Eq. (A.18)) and the atmosphere will still output a high heat flux. If the emissivity is smaller, the heat flux is limited by Eq. (A.18). The surface temperature will then increase and will adjust to temperatures closer to $T_p$. In turn, the temperature difference $T_p - T_s$ will be smaller and the convective heat flux will also adapt. Therefore, in our coupled model, high emissivities are associated with large temperature differences and thus, with low surface temperatures. On the contrary, small emissivities lead to small temperature differences, and therefore result in high surface temperatures.

The observed temperature differences result in larger fluxes for cases with a large emissivity, as displayed in Fig. 8d. The black body case presents the highest heat flux, followed by the low volatile content case. The high volatile content case presents the smaller heat flux.

**Fig. 8e** shows how the magma ocean depth evolves over time for each case. Here again, the three cases display the same behavior, albeit shifted in time. The depth remains constant during the first years, because the melt fraction remains above the rheological transition in the entire mantle. Once the mantle starts to solidify, the depth diminishes gradually. We observe a plateau in melt fraction when the depth reaches 573 km, which corresponds to the rheological front propagating in the transition zone. As shown by the melting curves displayed in Appendix Fig. A.13, this discontinuity is included in our model and therefore expected. Once the rheological front has reached the surface, we stop modeling the magma ocean and we consider the solidification time as the time taken for the rheological front to reach the surface. The solidification time is thus a result of the formation of the atmosphere, temperature adjustment and heat fluxes. As expected, the black body model is the faster to solidify, with a solidification time of 342 years. Next comes the low volatile content case with a
solidification of 2252 years. Finally, the high volatile content case displays the largest solidification time, with 21 723 years about 60 and 10 times longer than for the two previous cases.

Overall, the cases illustrated in Fig. 8 show how the atmosphere forms in a cooling magma ocean, at equilibrium with its atmosphere. The atmosphere is first dominated by the less soluble volatile (here, CO₂) before being dominated by the most abundant volatile (here, H₂O) towards the last stages. A magma ocean with low volatile content displays small opacities and high emissivities, associated with large temperature differences. This results in a high heat flux and a fast cooling. We thus illustrate that the coupling between a magma ocean and its growing atmosphere sets its evolution and its solidification time, consistent with previous studies (Lebrun et al., 2013; Salvador et al., 2017; Nikolaou et al., 2019; Lichtenberg et al., 2021).

## 5.2. Temporal evolution in the out-of-equilibrium case

In this Section, we show an illustration for an out-of-equilibrium outgassing case. Except of the outgassing scheme, all parameters correspond to the low volatile content case described in the previous Section. In the out-of-equilibrium outgassing case, the initial quantities in each reservoir need to be specified. Therefore, we initialize our computation with an atmospheric pressure of zero bar for every volatile, i.e., the volatile content in the magma ocean is equal to the initial volatile content $X_{0,l}$. During magma ocean evolution, CO₂ and H₂O will progressively outgas from the liquid mantle to the atmosphere. Because the partition coefficients between the liquid and solid mantle are small (see Table 3), knowing the atmospheric pressure in a given volatile allows one to estimate the volatile content in the liquid mantle. The convective velocities needed to account for volatile outgassing efficiency are computed using the following expression (Solomatov, 2015):

$$v_{\text{conv},0} = 0.6 \left( \frac{\alpha g D F_{\text{conv}}}{\rho C_p} \right)^{1/3} . \qquad (19)$$

Due to the high range of Ra over which we extrapolate, and following Solomatov (2015), we associate this expression with uncertainties of a factor three, $v_{\text{conv}} = v_{\text{conv},0} \times [1/3; 3]$, that we take into account in our computations. The results for the case computed with $C_p = 6.6$ are shown with solid blue lines in Fig. 9, while cases that take into account uncertainties on $v_{\text{conv}}$ are shown with fading blue or red shades.

For reference, we include an equilibrium outgassing case, with the same parameters. This case corresponds to the low initial volatile content, previously discussed in Section 5.1. Figs. 9a and 9b show that for both CO₂ and H₂O, the atmospheric pressure lags behind the pressure in the equilibrium case. The pressures start from a zero bar atmosphere and are slowly converging towards the equilibrium pressure. For CO₂, the equilibrium pressure is almost reached when the rheological front starts to progress upwards. For H₂O however, atmospheric pressures in the out-of-equilibrium case are still a hundred times lower than in the equilibrium case when the rheological front starts to progress upwards. Because the solubility of water is higher than that of CO₂, it has a lower saturation pressure and therefore, shallower exsolution depths (Eq. (16)).

Once the rheological front starts to progress upwards (after 573 years), both atmospheric pressures present a discontinuity in growing rate and start to increase faster. In the equilibrium case, the atmospheric pressure can only increase as fast as the rheological front progresses. For the out-of-equilibrium case, once the rheological front progresses upwards, the atmospheric pressures increase faster than the rheological front progresses, before slowing down at the end of the cooling stage. As the magma ocean solidifies, the Rayleigh number and the velocity also decrease. Additionally, the relative exsolution depth increases with decreasing the magma ocean thickness. The volatile concentration also increases, due to decreasing magma ocean extent. Small changes to velocity and viscosity combined with the increase of

the exsolution depth therefore explain the enhanced outgassing rate at the start of the solidification. As a result of this increased outgassing rate, the final atmospheric pressures of both CO₂ and H₂O are close for both equilibrium and out-of-equilibrium outgassing. Therefore, the main difference between the two outgassing schemes is not the final volatile amount in each reservoir but rather the evolutionary paths that quantities follow.

Fig. 9c shows that the surface temperature is initially substantially lower in the out-of-equilibrium case. While it decreases monotonously in the equilibrium case, it first increases in the out-of-equilibrium case, going from 3000K to 3130K in a few hundred years. Once the atmospheric pressure builds up, the atmosphere gradually cools down and follows the same evolution as in the equilibrium case. The overall temperature difference is always higher in the out-of-equilibrium case than in the equilibrium case.

As a result of inefficient outgassing during the first phase of the magma ocean, the atmosphere of the out-of-equilibrium magma ocean presents a higher emissivity (Fig. 9), leading to higher heat fluxes and faster cooling. The rheological front starts to progress upwards sooner for the out-of-equilibrium case (573 years vs. 922 years). The emissivity then gradually diminishes and reaches values similar to the values of the equilibrium case. The magma ocean evolution is then similar in both cases. The rheological front reaches the surface in 1688 and 2252 years for the out-of-equilibrium and equilibrium case respectively. The time difference between the two cases thus mainly comes from the first phase of the magma ocean, where the outgassing efficiency is low and cooling is fast.

These results show that for a magma ocean out of equilibrium with its atmosphere, the outgassing efficiency is limited in the first phase of the magma ocean where the entire magma ocean is essentially liquid. The most soluble volatile are less exsolved, independently of their respective abundances, and the atmosphere is first dominated by the least soluble ones. Tenuous atmospheres are associated with large temperature differences and the surface temperature can evolve in a non-monotonous way. The delayed atmospheric formation results in higher emissivities and a faster cooling. Because of a smaller depth, outgassing then increases and the atmospheric formation catches on. The difference in solidification time is then mostly set by the time spent in the first totally liquid phase. Delayed atmospheric formation results in faster solidification times, compared to the equilibrium case. While being for a particular case, the remarks made here are general and we find that out-of-equilibrium outgassing and delayed atmospheric formation always shorten the solidification time. However, depending on the choice of parameters, the difference in solidification time can vary considerably and a more systematic study is needed. This is addressed in the following Section.

## 5.3. Parameter study

To fully appreciate the potential influence of out-of-equilibrium outgassing, we apply the previous approach to a broader range of cases. We first include the presence of rotation in our model, with a rotation rate $\Omega = 10^{-3}$ rad s⁻¹. This value corresponds to a fast rotating Earth (Cuk and Stewart, 2012). The heat flux is then parameterized using Eq. (A.12). We also introduce the following length scale to Eq. (19).

$$l \approx \frac{v_{\text{conv}}}{\Omega}. \qquad (20)$$

Solving Eqs. (19) and (20) leads to the following expression for the convective velocities:

$$v_{\text{conv},0} \propto \left( \frac{\alpha g F_{\text{rot}}}{\rho C_p \Omega} \right)^{1/2} . \qquad (21)$$

As in the previous section, following Solomatov (2015), this expression is associated with uncertainties of a factor three, $v_{\text{conv}} = v_{\text{conv},0} \times [1/3; 3]$, that we take into account in our computations.
<figure><img src="image_13.png" /><figcaption>Fig. 9 Time evolution of the atmospheric pressure of (a) CO₂ and (b) H₂O, (c) of the temperatures, (d) of the emissivity and (e) of the magma ocean depth. The dashed lines show a case where the magma ocean is at equilibrium with its atmosphere and correspond to the low volatile content case in Fig. 8. The blue and red fading shades show cases where the outgassing is limited by convective motions and the magma ocean is not in equilibrium with its atmosphere. Blue and red fading shades were computed by taking into account uncertainties on v conv,0 and C s. In panel (c), the blue curves represent the potential temperature T p while red curves show the surface temperature T s. Apart from the outgassing scheme, all curves were computed with the same parameters. All curves were computed with k 0,CO 2 = 10 − 3 m $^{2}$/kg. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

To gain a more general insight, we considered terrestrial plan-
ets of different sizes (Ida and Lin, 2004). In our computations, the
key parameters of a planet are the radius R<sub>p</sub>, the gravity g and the
magma ocean depth. The latter is defined through the radius of the
core mantle boundary R<sub>c</sub>. Following scalings proposed in Valencia
et al. (2006), we set the planetary mass M<sub>p</sub> and compute the defining
parameters (Stamenković et al., 2012):

$$
R_p = R_{p,\text{Earth}} \left( \frac{M_p}{M_{\text{Earth}}} \right)^{0.27} \qquad (22)
$$

$$
R_c = R_{c,\text{Earth}} \left( \frac{M_p}{M_{\text{Earth}}} \right)^{0.247} \qquad (23)
$$

$$
g = g_{\text{Earth}} \left( \frac{M_p}{M_{\text{Earth}}} \right)^{0.46} . \qquad (24)
$$

The melting curves displayed in Appendix A (Fig. A.13) are only
valid for terrestrial mantle pressures, up to 140 GPa. For planets with
larger masses and mantle pressures, we extend the melting curves
with an expression for pressures higher than 140 GPa. We adopt the
parameterization of the solidus curve proposed in Stixrude (2014), that
we adjust to connect to the previous solidus curves at 140 GPa.

$$
T_{\text{sol}} = 5200 \frac{(P/140)^{0.48}}{1 - \ln(0.79)}, \qquad (25)
$$

where *P* is expressed in GPa. For the liquidus curve, we consider the
same variations with pressure, but shifted in temperature:

$$
T_{\text{liq}} = 5200 \frac{(P/140)^{0.48}}{1 - \ln(0.79)} + 1100\text{K}, \qquad (26)
$$

where *P* is expressed in GPa. The shift in temperature corresponds to the difference between liquidus and solidus temperature at 140 GPa, computed from Eq. (A.6). As shown in Appendix Figs. A.13 and F.18, temperature profiles intersect the melting curves at *P* ≈ 90 GPa, which might lead to a basal magma ocean when the resulting density profile allows for it (Labrosse et al., 2007). Until now, for simplicity, we considered that heavier crystals sank instantaneously to the bottom of the magma ocean and that no basal magma ocean formed. Here, we also consider the possibility for the development of a basal magma ocean (BMO). When a basal magma ocean is created, we consider that there is no mass exchange between the basal and surface magma ocean and that volatiles in the basal melt can therefore no longer be outgassed. We therefore update Eq. (14) by subtracting the volatile mass in the BMO to the total mass. The mass of volatile *i* segregated in the BMO is then expressed as:

$$
M_{\text{loss},i} = M_{\text{BMO}} X_{i,i=\text{BMO}}, \qquad (27)
$$
<figure><img src="image_14.png" /><figcaption>**Fig. 10.** Ratio of solidification times (equilibrium vs. non-equilibrium) for various planetary masses in magma oceans influenced by rotation. Red curves (resp. purple) correspond to cases where volatiles were segregated in a basal magma ocean (resp. where there was no BMO formation). The solid black lines show cases computed with $v_{\text{conv}} = v_{\text{conv},0}$ and $C_b = 6.6$. The fading lines show cases computed by accounting for uncertainties on both $v_{\text{conv}}$ and $C_b$. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

where $M_{\text{BMO}}$ is the mass of the BMO, determined by the melting curves, and $X_{\text{I,t}=\text{BMO}}$ is the volatile content at the time of the formation of the BMO. In other words, the formation of a basal magma ocean traps volatiles located at depths.

For each set of parameters, we compared equilibrium outgassing and out-of-equilibrium outgassing. We computed the solidification time – defined as the time when the rheological front reaches the surface – for equilibrium and out-of-equilibrium cases $t_{\text{sol,eq}}$ and $t_{\text{sol,neq}}$. We first discuss the effect of rotation rate, BMO trapping and planetary mass prior to showing the effect of volatile content and absorption coefficient of CO₂.

Fig. 10 shows the ratio of the solidification times for planetary masses ranging from one to ten Earth masses. We considered a low initial volatile content, and an absorption coefficient of CO₂, $k_{\text{O,CO}_2} = 5 \times 10^{-2} \text{ m}^2 \text{ kg}^{-1}$. We considered rotation, and the presence or the absence of BMO trapping. The median case ($v_{\text{conv}}$, $C_b$) = ($v_{\text{conv},0}$, 6.6) is shown in solid line while cases that account for uncertainties are shown with fading colors.

Fig. 10 shows that the ratio of solidification time increases with increasing planetary mass. Because of larger masses and heat reservoir, the equilibrium solidification time increases with increasing planetary mass. However, the out-of-equilibrium solidification time does not increase as much with increasing planetary mass, leading to less efficient outgassing for larger planets. Because of large Rayleigh numbers and deep magma ocean, one could expect the convective velocities to be sufficiently high to guarantee efficient outgassing. However, as shown above, the outgassing time is set by the magnitude of the velocities, but also by the exsolution depth. With increasing planetary mass and gravity, the pressure gradient also increases. For a given volatile content, the exsolution depth therefore decreases with planetary mass (Eq. (16)).

It is therefore the competition between high velocities and shallow exsolution depths that leads to the observed faster solidification and less efficient outgassing for large planetary masses.

By comparing the red and purple curves in Fig. 10, one can assess the influence of BMO trapping on thermal evolution of planetary bodies. The ratio of solidification time increases in the presence of BMO trapping. Following the same previous reasoning above, this means that outgassing is less efficient in the presence of BMO trapping. Because BMO trapping isolates by definition a part of the magma ocean from the atmosphere, it decreases the volatile masses available for outgassing, resulting in a thinner atmosphere than in the case without BMO trapping. The difference between the equilibrium and out-of-equilibrium cases stems from the fact that upon trapping, volatile contents are different. As mentioned above, in the out-of-equilibrium case, the atmospheric pressures are lagging, compared to the equilibrium case. This means that at a given thermal state, the volatile contents in the magma ocean are higher for the out-of-equilibrium case. BMO trapping in the out-of-equilibrium case thus traps melt that is richer in volatiles, thereby amplifying the difference in solidification times.

Fig. 10 also illustrates the effect of convective velocities on outgassing. The overall trends are similar for all parameters (BMO trapping or not, planetary mass) and the comments previously made stand. There is however a strong difference in the solidification time ratios for different velocity estimations. When using the median velocity estimation, the ratios of solidification times range between 5 and 19 for all cases. When using the lower bound of velocity and exsolution coefficient estimations, the ratios are higher in each case, ranging from 6 to 30 (Fig. 10). Because in the equilibrium case, the convective velocities play no role, comparing the central line and faded red or purple shades is equivalent to comparing the effect of velocity amplitude on out-of-equilibrium outgassing only. All parameters equal otherwise, a smaller
<figure><img src="image_15.png" /><figcaption>**Fig. 11.** Parameter study of solidification times for equilibrium and non-equilibrium outgassing. The parameters made to vary were the volatile content in both CO₂ and H₂O, and the CO₂ absorption coefficient. Red and blue curves (resp. yellow and purple) correspond to cases where volatiles were sequestered in a basal magma ocean (resp. where there was no BMO formation). Red and purple curves (resp. blue and yellow) correspond to a rotating case (resp. a non-rotating case). Contrary to Fig. 10, these plots do not take into account uncertainties in the velocity amplitude and only consider the exsolution coefficient C_b = 6.6. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

velocity amplitude thus leads to slower outgassing. This illustrates the importance of velocity amplitudes on outgassing efficiency.

As demonstrated in Fig. 8, the volatile content can have tremendous effects on the solidification time of magma oceans. We therefore varied the initial volatile content, as done above (Fig. 8). For H₂O, we use initial contents of 41 ppm (low volatile content) and 410 ppm (high volatile content). For CO₂, we use initial contents of 13 ppm (low volatile content) and 130 ppm (high volatile content). Another key parameter is the absorption coefficient of CO₂. In parameterized coupled magma ocean-atmosphere models, different values have been adopted, ranging from k₀,CO₂ = 10⁻⁴ m²/kg to k₀,CO₂ = 5 10⁻² m²/kg (Elkins-Tanton, 2008; Lebrun et al., 2013; Nikolaou et al., 2019). The absorption coefficient remains unconstrained and can play a large role on atmospheric thermal equilibrium and the associated outgoing heat flux. We therefore study its role on thermal evolution of the magma ocean.

We adopt the same approach as above: for each set of parameters, we compute the time evolution of the magma ocean in the equilibrium and the out-of-equilibrium case and compute the ratio of the solidification times. Fig. 11 shows the evolution of this ratio for planetary masses ranging from 1 to 10 Earth mass. We also compare rotating and no-rotating cases and we include the possibility for basal magma ocean trapping.

The left column of Fig. 11 shows the results for a low volatile content, while the right column shows the results for a high volatile content. The first row was obtained for the lowest value of the absorption coefficient, k₀,CO₂ = 10⁻⁴ m²/kg, and the bottom row for k₀,CO₂ = 5 10⁻² m²/kg. For all cases, increasing the planetary mass increases the ratio of solidification time. The outgassing efficiency drops for deep magma oceans and therefore, for large planets. The solidification time is always shorter for the out-of-equilibrium cases. Note that contrary to Fig. 10 these plots do not take into account uncertainties on the velocity amplitude and on the exsolution coefficient. Consequently, the influence of non-equilibrium outgassing would be even more pronounced than what is displayed in Fig. 11 for cases computed with lower velocities and exsolution coefficient.

For high volatile contents (right column in Fig. 11), the median solidification time ratios range between 1 and 4. The non-rotating cases are the most efficient at outgassing volatiles (yellow and blue curves), and the rotating cases are the least efficient at outgassing volatiles (purple and red curves). For high volatile contents, larger values of k₀,CO₂ are associated with solidification time ratios closer to 1, indicative of a more efficient outgassing in the non-equilibrium case. Non-rotating cases are less affected by the increase of k₀,CO₂ than fast rotating cases. When the volatile content is high, even for inefficient outgassing, the atmosphere reaches relatively high pressures rapidly (a few bars of CO₂ in a few tens of years). Because of this non-zero atmospheric pressure, increasing the absorption coefficient effectively reduces the emissivity and slows down cooling. This decrease in solidification ratio with increasing absorption coefficient is not surprising.

For low volatile content however (left column in Fig. 11), the solidification ratios increase with increasing the absorption coefficient. They range between 2 and 5 for low values of k₀,CO₂ and between 4 and 20 for high values of k₀,CO₂. Because the volatile content is low, the atmospheric pressure is much lower than in the high volatile content case for out-of-equilibrium outgassing. At low atmospheric pressures, the opacity is small regardless of the value of the absorption coefficient. Changing the absorption coefficient therefore only plays a minor role in out-of-equilibrium thermal evolution of planetary bodies with low volatile content. The increase in solidification time ratio is then driven by the increase in solidification time for the equilibrium cases, which is well explained by increasing the absorption coefficient.

For the non-rotating cases (blue and yellow curves), the effect of out-of-equilibrium outgassing is not quite visible when only comparing the solidification time ratios. For these cases, high velocities lead to early outgassing of CO₂ which lowers the emissivity rapidly and extends the magma ocean lifetime. The prolonged lifetime then enables more efficient water outgassing. The evolutionary path of volatile quantities in non-rotating cases are then close to the ones in the equilibrium case. Therefore, upon formation of the basal magma ocean, the differences between the equilibrium and non-equilibrium cases are small, which
<figure><img src="image_16.png" /><figcaption>Fig. 12. Time evolution of the ³He mass normalized by its initial abundance in each reservoir (yellow, red and blue curves for magma ocean, basal magma ocean and atmosphere respectively) for an out-of-equilibrium evolution of a terrestrial magma ocean with low volatile content. The mass is normalized with respect to the initial mass of ³He. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

explains the closeness of the blue and yellow curves. When looking only at the solidification time, we are however blind to the outgassing path of rare elements and other later-outgassed volatiles, such as ³He. The effect of non-equilibrium outgassing on rare volatile species and their distribution is discussed below.

## 5.4. Non-equilibrium outgassing of noble gases: The case of ³He

To illustrate the effect of out-of-equilibrium outgassing on the outgassing path of ³He, we perform an additional computation. We chose ³He because it is a key component to the interpretation of geochemical data that suggest the presence of hidden long-lived reservoirs rich in rare gases sampled by modern mantle plumes (Weis et al., 2023). Here, we focus on the Earth's mantle case, without accounting for rotation, which is the least favorable case for volatile exsolution. We set the initial volatile content in CO₂ to 13 ppm and the initial volatile content in H₂O to 41 ppm. The absorption coefficient for CO₂ is set to k₀,CO₂ = 5 10⁻² m² kg⁻¹. We consider the case where a basal magma ocean is created. We compute the saturation pressure of ³He as Heber et al. (2007):

$$P_{\text{sat},^3\text{He}} = \frac{X_{^3\text{He}}}{10^{-15}} \qquad (28)$$

Following Özgür el and Caracas (2023), we set the initial ³He concentration to 4 10⁻⁵ ppm. Fig. 12 shows the time evolution of ³He in each reservoir (magma ocean, basal magma ocean and atmosphere) for the out-of-equilibrium case. The equilibrium case is described but not shown, as it results in a state that does not evolve with time.

Helium is initially entirely dissolved in the magma ocean (yellow patch). As the magma ocean convects, volatiles are outgassed at a slowly varying rate (Eq. (18)), and the mass of ³He in the atmosphere grows (blue patch). This growth occurs at the expense of the volatile content in the magma ocean that steadily decreases over time. The temperature profile crosses the curve $\phi = \phi_c$ for the first time at a depth

2328 km, after 3700 years. This results in the creation of a basal magma ocean, that extends from this depth to the Core-Mantle Boundary. At this stage, we consider that there is no mass exchange between the two oceans (surface and basal). This results in a fraction of the dissolved ³He being trapped in the basal magma ocean (red patch). Then, the trapped part does not evolve over time and the content evolution of the surface magma ocean and atmosphere are linked through Eqs. (14) and (18), with a modified mass budget.

The amount of trapped ³He depends on the volatile content in the magma ocean at the time of creation of the basal magma ocean and on the extent of the basal magma ocean. In this case, the volume fraction of the BMO is equal to 11%. In the case shown in Fig. 12, the relative volatile content in the magma ocean is high, because small volatile content and low saturation pressures reduce outgassing efficiency. This leads to nearly 10% of the initial mass of ³He being trapped in the basal magma ocean, considerably affecting the mass budget. Conversely, in the equilibrium case, because ³He presents a low solubility (Eq. (28)), it is almost entirely outgassed instantaneously from the entire magma ocean. The volatile content in the magma ocean is therefore extremely low when the basal magma ocean forms and only a small fraction of ³He is trapped. We computed the trapped fraction to be equal to 7.7 10⁻⁴ % of the initial mass in ³He. This is orders of magnitudes lower than in the out-of-equilibrium case. Our model therefore illustrates the effect of out-of-equilibrium outgassing on the creation of chemical reservoirs in a basal magma ocean. Despite the extreme simplicity of the BMO trapping mechanism we considered, this illustrates how out-of-equilibrium outgassing can significantly contribute to the creation of rare gases reservoirs that may persist to the present day, as suggested by the analyses of Ocean Island Basalt samples (e.g., Weis et al., 2023 and references therein). Note that the outcome described above for Helium would be similar if we had considered a high initial volatile content in CO₂ and H₂O, as well as smaller absorption coefficient for CO₂. Therefore, regardless of the various possible conditions experienced by
terrestrial magma oceans, accounting for out-of-equilibrium outgassing likely significantly influences the thermal and/or chemical evolution of the magma ocean.

# 6. Discussion

Our modeling approach relies on several assumptions that we review below, prior to discussing the influence inefficient outgassing could have on planetary evolution beyond solidification times.

## 6.1. Model limitations

**Validity of the linear law.** The outgassing time is inversely proportional to the exsolution depth (Fig. 6b). This proportionality is valid for small exsolution depths, and remains valid for the major part of the life of the magma ocean. However, when the rheological front approaches the surface, its rapid ascent combined with inefficient outgassing leads to a super-saturated and extremely shallow magma ocean. The exsolution depth therefore rapidly increases above $d_{exsol} = 0.01$ D and the proportionality is no longer valid. Because of the way the exsolution coefficient evolves with the exsolution depth, we therefore tend to overestimate the outgassing efficiency of the magma ocean during its last evolutionary stages. A more accurate description of the exsolution coefficient would lead to smaller outgassing efficiency and would strengthen difference with the equilibrium case. Yet, since most of the volatiles have already outgassed when the rheological front approaches the surface, the emissivity is already low and the overestimated outgassing would only slightly decrease the already low emissivities. We therefore expect that we would not see significant changes with a higher order description of the exsolution coefficient.

**Magma ocean regime.** We found that the exsolution coefficient depends only on the Reynolds number. It first increases with the Reynolds number, before seemingly reaching a plateau at $C_b = 6.6$, for Reynolds number larger than $Re \sim 5 \times 10^3$. We interpreted this plateau as a result of turbulent flow. When the flow is not turbulent, mixing is inefficient and some regions may be separated from the surface. Moreover, the vertical velocity is relatively low, leading to small exsolution coefficients. For higher Reynolds number and more turbulent flows, mixing is more efficient and these regions break down, leading to more efficient exsolution. A more accurate description and characterization of separated regions of the flow as a function of the Reynolds number and experiments at higher $Re$ should be conducted to better explain the observed plateau. These will be the focus of a future study. However, because the flow patterns of convection have been largely described in the literature and because of the threshold value that separates increasing exsolution coefficient value from the plateau ($Re_c \sim 5 \times 10^3$), we are confident that turbulence is developed enough and that transport will be the main limiting phenomenon in the magma ocean.

**Interaction between volatile species.** Our parameterized cooling models do not account for the interactions between different volatiles (Bower et al., 2019; Massol et al., 2023), where the atmospheric mass of a given volatile species depends on the mean molecular mass of the atmosphere. While accounting for volatile interactions through partial pressures is more physically sound, results from Massol et al. (2023) tend to show that the influence on solidification times is relatively small in terms of magma ocean solidification time. Indeed, when modeling volatile interactions, Massol et al. (2023) found a 16.4% increase in solidification time, due to larger atmospheric pressures and larger opacities. They also found that volatile interactions modify the final volatile repartition between the atmosphere and the mantle, leading to higher water pressures. We therefore expect the type of modeling of the atmosphere to play only a marginal role, in particular on the cooling timescales of the magma ocean. We verified that this was indeed the case by performing the same calculations, but accounting this time for the effect of volatile interactions, following Bower et al. (2019)

and Massol et al. (2023). The results of these calculations are shown in Appendix B. We found no significant differences in magma ocean solidification times for cases with equilibrium or out-of-equilibrium outgassing. This is because during the first stages of magma ocean solidification, only $CO_2$ is outgassed and the atmosphere is therefore essentially composed of a single species. Later, when water begins to outgas, multi-species populate the atmosphere. However, because the emissivity is already low at this later stage, changes in composition due to multi-species interactions play little to no role on the outgoing heat flux given our assumption of a gray atmosphere. While the role of volatile interaction may be more pronounced for prolonged stages occurring after the solidification front reaches the magma ocean's surface and including atmospheric escape, our present study did not investigate this time period and such processes.

**Importance of the solubility laws.** To compare our results with previous studies (Lebrun et al., 2013; Salvador et al., 2017; Nikolaou et al., 2019), we described water solubility using the same solubility law (Eq. (13)). Sossi et al. (2023) however showed that water solubility is higher for compositions closer to that of the terrestrial magma ocean. For completeness, we therefore performed the same calculations using these more accurate solubility laws. As a result of the higher solubility, we found that water remains dissolved in the magma ocean longer for both cases (equilibrium and non-equilibrium). This resulted in overall slightly higher emissivities and slightly decreased solidification times for all cases (from 2257 years to 1986 years for the equilibrium case and from 1712 years to 1567 years for the non-equilibrium case). Because a higher water solubility will lead to smaller exsolution depths, we therefore expect water outgassing to be less efficient when using the solubility law from Sossi et al. (2023). This would therefore amplify the differences between the equilibrium and non-equilibrium cases.

**Magma ocean mixing.** Our outgassing model consists of an initially homogeneous magma ocean that vigorously convects and transports fluid parcels to progressively shallower exsolution depths. This constitutes a reasonable approximation for most of the magma ocean lifetime, because we expect mixing to be efficient. Upon magma ocean solidification, oversaturated volatiles in fluid parcels are created close to the rheological front, to the bottom of the magma ocean, and are exsolved at the surface of the magma ocean, where the pressure is low. For a magma ocean that is not perfectly well mixed, one should expect a saturation gradient along with lateral variations in the magma ocean. This would correspond to a situation different than the one we modeled. To more accurately describe outgassing efficiency, a thermochemical model including mixing and chemical evolution would be needed.

**Basal magma ocean and chemical equilibration.** In our model, we considered that while the atmosphere is out of equilibrium with the magma ocean, the solidified mantle fully equilibrates chemically with the melt through partition coefficient $\kappa$ (Eq. (14)), which represents an end-member for mantle equilibration. However, during mantle solidification, a non-negligible fraction of the melt can remain trapped in the solid due to inefficient compaction (Hier-Majumder and Hirschmann, 2017; Sim et al., 2024). In this scenario, large fractions of the initial volatile content can therefore remain sequestered in the solid mantle after solidification. In that sense, our basal magma ocean scenario is similar to the melt trapping one, because it prevents a fraction of the volatiles from being outgassed, and the lower part of the mantle is not chemically equilibrated with the shallow magma ocean and atmosphere. Still, the two scenarios present differences that may be important. In the BMO scenario, we used a simplified model, where upon creation of a mid-ocean solid layer, the bottom part remains separated from the atmosphere. The trapped fraction therefore depends only on the relative steepness of the melting curves and of the adiabat and on the extent of the magma ocean. The partial melt trapping scenario, on the other hand, is based on the physical understanding of compaction dynamics and may offer a more consistent volatile segregation model and coupling the out-of-equilibrium outgassing with melt trapping could be an interesting future prospect.
**Geometry and rotation.** In our model, we estimate the outgassing time scales from convection experiments in a Cartesian geometry. Planetary magma oceans present a spherical geometry quite different from that of our Cartesian experiments. Moreover, while our parameterized cooling models account for rotation, our fluid dynamics experiments do not take into account the effect of rotation on the dynamics to derive our scalings. We believe that because the value of the exsolution coefficient is determined by turbulent mixing efficiency, a change in geometry would only moderately affect the value of the exsolution coefficient as long as the aspect ratio of the magma ocean remains on the order of unity. Such changes would be negligible compared to uncertainties on the velocity, and should therefore not alter our conclusions. While we are confident that the overall outgassing dynamics would be similar in the presence of rotation (first order differential equation and outgassing time scales inversely proportional to the exsolution depth), the value of the exsolution coefficient $C_b$ could be different.

Moreover, rotational dynamics introduces different regimes (rotation-dominated vs. convection-dominated), and the behavior of the exsolution coefficient could depend on these regimes. The magma ocean may cool down through different regimes and present various solidification and outgassing rates. We plan to extend this work by explicitly including the effects of rotation on outgassing dynamics, including spherical geometry.

**Other contributions to inefficient volatile outgassing.** We characterized convective outgassing efficiency by focusing on the ability of magma ocean fluid parcels to reach exsolution depths. As discussed in the introduction section, we did not consider the inertia associated with the formation of bubbles above the saturation depth, which can be substantial, nor the possibility of volatile dissolution due to vigorous downwelling motions. For these reasons, the results and the influence of out-of-equilibrium volatile outgassing presented here should be considered as lower bound, and one could expect that out-of-equilibrium volatile outgassing in magma oceans to be even more pronounced. These aspects are the focus of ongoing investigations.

## 6.2. Implications of our results

Our results show that contrary to the assumption of equilibrium outgassing systematically made in previous studies, magma ocean outgassing can be quite inefficient, with strong effects on the solidification time scales of terrestrial bodies, but also on their cooling paths.

When the mantle solidifies, iron partitioning between melt and solidified silicates can lead to the formation of a gravitationally unstable density stratification, that might lead to mantle overturning (Elkins-Tanton et al., 2003). Depending on the timing of the onset of convection in the solid mantle (e.g. before or after mantle solidification is complete) the mantle can exhibit vastly different structures (Maurice et al., 2017; Ballmer et al., 2017; Boukaré et al., 2018). If the solid mantle starts to convect before solidification is complete, early mantle heterogeneities will be progressively homogenized during the course of the solidification, resulting in a relatively well mixed mantle at the end of solidification. On the other hand, if the solid mantle does not convect while it solidifies, this might result in a single overturn event that would then lead to a stable layered mantle with a composition distinct from the overlying mantle in both major (i.e., iron), minor, and trace elements. Such a differentiation mechanism has been advocated to explain the presence of different degrees of compositional heterogeneity within the Earth's mantle sampled by mantle plumes (Weis et al., 2023), or more recently detected on Mars by seismic waves (Samuel et al., 2023).

Because large planets with low volatile content are less efficient at outgassing their atmospheres, they should solidify rapidly and they could therefore exhibit different mantle structures and long-term evolution from smaller planets with higher volatile content and slow cooling. In turn, the structure of the mantle would have tremendous effects on the long-term evolution of the mantle.

Additionally, due to the decrease of viscosity with water content, the final water content in the mantle can also play an important role in the dynamics and the long-term evolution of the mantle (Regenauer-Lieb et al., 2001; Foley and Driscoll, 2016).

The partitioning of volatiles at the end of magma ocean solidification is also essential to assess the habitability of exoplanets. The relative abundances of CO₂ and H₂O in planetary atmospheres after mantle solidification set the surface conditions of rocky planets in terms of pressure and temperature and the possibility of condensation of a liquid water ocean (Salvador et al., 2017; Massol et al., 2023). Mantle water storage due to inefficient outgassing then makes it difficult to condense a water ocean once the planet has cooled. Taking into account updated volatile budgets due to inefficient outgassing could then shift the conditions for which a water ocean can condense and therefore influence the potential habitability of a planet (Hart, 1979; Kasting et al., 1993).

The evolutionary path of exoplanets is also responsible for their observed diversity. Hamano et al. (2013) theorized the existence of two types of terrestrial exoplanets, depending on the distance to their star. Type I planets correspond to planets such as the Earth, relatively far from their star, that cooled rapidly. For type II planets however, closer to their star, the radiation limit of the atmosphere balances the incoming radiation from the star, resulting in a small outgoing heat flux and the solidification time of the magma ocean is therefore considerably longer. The large solidification time therefore leads to massive water loss through hydrodynamic escape, resulting in a large decrease in the total water inventory. Recently, using a more realistic climate model, Selsis et al. (2023) showed that atmospheres of exoplanets should be colder than what was proposed in Hamano et al. (2013), and that planets would therefore need to receive higher stellar radiation to be labeled type II. Even if the phenomena responsible for the water loss (hydrodynamic escape and radiation limit) are not taken into account in the present study, we can still hypothesize that inefficient outgassing would change the outcome described in Hamano et al. (2013). In the case of inefficient outgassing, solidification is expected to be fast, thereby reducing the water loss by hydrodynamic escape. Moreover, because water remains dissolved in the magma ocean for most of its lifetime, the radiation limit is only expected to be reached at the end of magma ocean solidification. A type II planet could therefore solidify more rapidly than estimated in Hamano et al. (2013), resulting in a smaller water loss. Lastly, water storage in the mantle for longer time scales due to inefficient outgassing could prevent water from hydrodynamic escape, leading to atmospheres that are more water-rich at the end of solidification. Because we expect outgassing to be efficient for large water contents, the evolution of type II planets is likely to be similar to that described in Hamano et al. (2013). On the other hand, type II planets with lower water contents would be more likely to follow the evolution path of type I planets.

## 7. Conclusions

Using laboratory and numerical experiments, we reproduced the convective dynamics of planetary magma oceans, and we characterized the evolution of the exsolved fraction of volatiles species as a function of the governing parameters. We introduced a parameterization of the flux of outgassed volatiles depending on the exsolution depth. For a magma ocean of thickness $D$ convecting with velocities $v$, the evolution of the fraction of exsolved volatiles $\chi$ obeys a first degree differential equation $\dot{\chi} = v B (1-\chi)/D^2$. We found that the inverse of the exsolution time was proportional to the exsolution depth through $B = C_b d_{\text{exsol}}$. The exsolution coefficient $C_b$ increases as a function of the Reynolds number $Re$ and reaches a plateau $C_b = 6.6$ for Reynolds number values higher than $Re \sim 5 \cdot 10^3$. Our model describes both the outgassing time scales and the outgassing flux, which represents an improvement and a generalization of the model in Salvador and Samuel (2023) that estimated outgassing time scales, which were only applicable for a constant magma ocean depth and a fixed Prandtl number value of unity.
We applied our exsolution model to a magma ocean cooling model, coupled with its atmosphere. This coupled model accounts for the first time for the influence of out-of-equilibrium outgassing on the thermal evolution of the magma ocean and in particular on the solidification time of planetary bodies. We showed that out-of-equilibrium outgassing can reduce solidification times by more than one order of magnitude depending on the context, due to delayed atmosphere formation and high subsequent emissivities.

The effect of out-of-equilibrium is stronger for large planets that present a basal magma ocean, high rotation rates and low volatile content. For low volatile content, we find that higher absorption coefficients of CO₂ are associated with larger effects, with respect to the equilibrium case. For higher volatile contents, relatively fast outgassing of CO₂ into the atmosphere inverses the effect of higher absorption coefficients.

Overall, in line with recent findings (Salvador and Samuel, 2023) our results confirmed that in various contexts the common assumption of systematic equilibrium outgassing in terrestrial magma oceans is not a valid hypothesis.

### CRediT authorship contribution statement

**Arthur Walbecq:** Writing – review & editing, Writing – original draft, Visualization, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Henri Samuel:** Writing – review & editing, Validation, Supervision, Software, Methodology, Funding acquisition, Conceptualization. **Angela Limare:** Writing – review & editing, Validation, Supervision, Methodology, Funding acquisition, Conceptualization.

### 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 thank the editor, Paolo Sossi, and two anonymous referees, whose comments improved our manuscript. This study was funded by the IdEx Université Paris Cité, ANR-18-IDEX-0001, by the Programme National de Planétologie (PNP) of CNRS/INSU, and by the CNRS/INSU Programme et Équipements Prioritaires de Recherche (PEPR-Origins). We acknowledge partial support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement no. 101019965-ERC advanced grant SEPtiM). Numerical computations were performed on the S-CAPAD/DANTE platform, IPGP, France.

### Appendix A. Magma ocean cooling model

We describe below the set of equations that we used for our 1D parameterized model. We consider a planet of radius $R_p$ composed of a metallic core of radius $R_c$ and a convecting mantle, that is initially completely molten and perfectly mixed. The bottom of the magma ocean is then initially described by $R_{b,0} = R_c$. The temperature profile in the mantle is then adiabatic and is described by:

$$ \frac{dT}{dP} = \frac{\alpha T}{\rho C_p}, \qquad (A.1) $$

where T is the temperature profile, which varies with the radius r, P is the pressure, α is the thermal expansivity, ρ is the density and C_p is the heat capacity. The thermal expansivity varies with pressure (Abe, 1997):

$$ \alpha = \alpha_0 \left[ \frac{P K'}{K_0} \right]^{\frac{1-K'}{K'}}, \qquad (A.2) $$

where $α_0$ is the surface thermal expansivity, $K_0$ the surface bulk modulus and $K'$ is the pressure derivative of $K_0$. We numerically solve Eq. (A.1) assuming a hydrostatic pressure profile $P = \rho g(r - R_p)$ with constant gravity g and obtain the temperature profile associated with the potential temperature $T_p$. At all times we determine the melt fraction profile using $T(r)$ and the melting curves. We adopt a pressure-dependent parameterization of the melting curves referred to as synthetic parameterization in Nikolaou et al. (2019). These melting curves are a combination of different experimental data for the range of pressure of the terrestrial mantle. The melting curves are expressed as follows:

$$ T_{\text{sol}} = 1120.661 + 273.15 + 132.899P - 5.904P^2, \qquad (A.3) $$

for $P \le 2.7$ GPa (Hirschmann, 2000).

$$ T_{\text{sol}} = 1086 + 273.15 - 5.7P + \ln(P), \qquad (A.4) $$

for $2.7 \le P \le 22.5$ GPa (Herzberg et al., 2000).

$$ T_{\text{liq}} = 2014.497 + 37.743P - 0.472P^2, \qquad (A.5) $$

for $P \le 22.5$ GPa (Zhang and Herzberg, 1994).

$$ T_{\text{sol}} = 1762.722 + 31.595P - 0.102P^2 \qquad (A.6) $$

$$ T_{\text{liq}} = 1803.547 + 50.810P - 0.185P^2, \qquad (A.7) $$

for $22.5 \le P \le 140$ GPa (Fiquet et al., 2010).

The melt fraction φ is estimated as (e.g., Abe, 1997):

$$ \phi = \frac{T - T_{\text{sol}}}{T_{\text{liq}} - T_{\text{sol}}}, \qquad (A.8) $$

where $T_{sol}$ and $T_{liq}$ are the solidus and liquidus temperatures at pressure P. For $T > T_{liq}$, φ is equal to 1 and for $T < T_{sol}$, φ is equal to 0. For any melt fraction between 0 and 1, crystals are present in the melt. When the melt fraction is close to 1, the presence of crystals only slightly influences the rheology and the mantle behaves like a liquid (Lejeune and Richet, 1995). For melt fractions close to 0, a larger crystal fraction is present and the behavior changes to a solid-like behavior. The critical melt fraction that separates these two regimes ranges between 0.3 (Hier-Majumder and Hirschmann, 2017; Maurice et al., 2017) and 0.5 (Monteux et al., 2016; Ballmer et al., 2017) in the literature. Here, we adopt the fixed value $φ_c = 0.4$, often used in previous parameterized cooling models Lebrun et al. (2013), Bower et al. (2018) and Nikolaou et al. (2019). Fig. A.13 shows the melting curves as a function of pressure, as well as three adiabats, at $T_p = 4000$ K, $T_{p} = 3300$ K and $T_p = 2500$ K (red, orange and yellow curves, respectively).

Upon mantle cooling, the temperature profile initially above the liquidus shifts towards colder values and eventually crosses the melting curves. As a result, the melt fraction at a given pressure decreases with time, and crosses the critical value $φ = φ_c$ (shown in gray in Fig. A.13). Depending on the relative steepness of the adiabat and of the melting curves, they might intersect at the bottom or in the middle of the mantle. If the curves intersect at the bottom, the rheological front will progress upwards and the evolution of the magma ocean depth will be continuous. If the curves intersect in the middle of the mantle, the way the rheological front propagates will depend on the density structure. If crystals that form in the middle of the magma ocean are denser than the surrounding liquid, they may sink to the bottom of the mantle and the rheological front will progress upwards. On the other hand, if the density profile is stable, crystals can remain close to the location where they are formed, leading to the formation of a basal magma ocean and a rheological front that progresses both upwards and downwards from the middle of the mantle.

#### A.1. Magma ocean viscosity

One of the key governing parameters of magma ocean dynamics is the viscosity, η. It affects both the Rayleigh and Prandtl numbers,
<figure><img src="image_17.png" /><figcaption>Fig. A.13. Melting curves and temperature profiles for three isotherms $T_p = 4000$ K, $T_p = 3300$ K and $T_p = 2500$ K. The black dotted line shows the liquidus, the black dashed line shows the solidus and the gray dashed line shows the rheological transition for a critical liquid content $\phi_c = 0.4$. The expressions used are described in Appendix A. All the curves were computed by solving simultaneously Eqs. (A.1) and (A.2) along with a lithostatic pressure profile, with $g_{\text{planet}} = g_{\text{Earth}} = 9.81 \text{ m s}^{-2}$.</figcaption></figure>

the dimensionless parameters that govern the convective motion. It varies with temperature and increases abruptly when the melt fraction drops below the critical melt fraction. For melt fractions above $\phi_c$, we model the viscosity profile using the following Vogel-Fulcher-Tammann (VFT) equation (Karki and Stixrude, 2010):

$$ \eta_l = A_k \exp\left(\frac{B_k}{T - C_k}\right), \qquad (A.9) $$

where the values of $A_k$, $B_k$ and $C_k$ are listed in Table 3. We include the effect of crystal content using the following expression (Roscoe, 1952):

$$ \eta = \frac{\eta_l}{\left(1 - \frac{1-\phi}{1-\phi_c}\right)^{2.5}}. \qquad (A.10) $$

In Nikolaou et al. (2019), the viscosity was computed through a volumetric harmonic mean to account for viscosity variations in the magma ocean. Because the convective motions in this context are essentially caused by destabilization of the upper thermal boundary layer, we preferred as in Lebrun et al. (2013) to use the viscosity and associated Rayleigh numbers computed at $T = T_p$. The difference between the two expressions are estimated in Appendix E, and were however found to be relatively minor. The above expressions are only used to describe the liquid-like behavior of the mantle for melt fraction $\phi > \phi_c$. For melt fraction $\phi < \phi_c$, the mantle is considered as solid and we neglect thermal evolution compared to that of the magma ocean. This assumption is justified by the fact that the solid mantle has a convective vigor that is orders of magnitude smaller than that of the molten magma ocean and will be further discussed below.

## A.2. Heat fluxes and thermal evolution

The convective heat flux is parameterized using two distinct scaling laws. The first one corresponds to the so-called ‘soft’ turbulence regime and yields a power law (Solomatov, 2015):

$$ F_{\text{soft}} = 0.089 \frac{\lambda(T_p - T_s)Ra^{1/3}}{D}. \qquad (A.11) $$

In its early history, the Earth might have rotated faster, with rotation rate up to ten times higher than today’s rotation rate (Agnor et al., 1999; Kokubo and Ida, 2007; Kokubo and Genda, 2010; Ćuk and Stewart, 2012). The presence of rotation considerably affect the dynamics of Rayleigh-Bénard convection, both in terms of heat flux and convective velocities, and might have played an important role during mantle solidification (Moeller and Hansen, 2013; Maas and Hansen, 2015; Gastine et al., 2016; Maas and Hansen, 2019; Thomas et al., 2024). For a rapidly rotating terrestrial planet, we therefore introduce the scaling proposed in Gastine et al. (2016):

$$ F_{\text{rot}} = 0.15 \frac{\lambda(T_p - T_s)}{D} Ra^{3/2} Ek^2 Pr^{-1/2}, \qquad (A.12) $$

where $Ek = \nu/(\Omega D^2)$ is the Ekman number, with $\Omega$ the rotation rate in rad/s. The smaller the Ekman number, the larger are the effects of rotation on the dynamic of the flow.

Following Elkins-Tanton (2008), we used a gray atmosphere model and considered two gases, H₂O and CO₂ that absorb and radiate strongly. For H₂O, we considered a single absorption coefficient at normal atmospheric conditions ($T_0, P_0$) $k_{0,H_2O} = 10^{-2} \text{ m}^2/\text{kg}$. For CO₂, we adopt different absorption coefficients at normal atmospheric conditions ($T_0, P_0$), ranging from $k_{0,CO_2} = 10^{-4} \text{ m}^2/\text{kg}$ to $k_{0,CO_2} = 5 \ 10^{-2} \text{ m}^2/\text{kg}$. The effect of the absorption coefficient is discussed more extensively in Section 5.3.

At the top of the atmosphere, the incoming flux from the star is set by the incident radiation $S_0$.

$$ F_{\text{Star}} = (1 - A) \frac{S_0}{4} = \sigma T_{\text{eq}}^4, \qquad (A.13) $$

where A is the albedo, $\sigma$ is the Stefan-Boltzmann constant and $T_{\text{eq}}$ is by definition the black body equilibrium temperature, i.e., the temperature at which the planet would be in equilibrium if it radiated a flux equal to the one received from the star. Because the planet is cooling down, its surface temperature is higher than the equilibrium temperature $T_p > T_{\text{eq}}$, and it radiates a flux $F_{\text{gray}}$. This flux depends on the volatile content
<figure><img src="image_18.png" /><figcaption>Fig. A.14. Time evolution of the heat flux computed at the surface of an Earth-sized planet (blue) and the heat flux estimated at the interface between the mushy/essentially solid mantle and the fully/essentially molten magma ocean (i.e., located at $r = R_f$ (in red). Gray shades correspond to a temperature difference across the interface $r = R_f$ ranging from 1 K to $10^3$ K. This plot corresponds to one of the cases shown in Fig. 8: a low volatile content (13 ppm and 41 ppm in CO₂ and H₂O respectively), with non-equilibrium outgassing and a gray atmosphere. The CO₂ absorption coefficient was set to $k_{0,CO_2} = 10^{-3}$ m² kg⁻¹. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

in the atmosphere. This dependency is set by the emissivity $\epsilon$, using the following set of equations. We compute the opacity of the atmosphere $\tau_i$ for each volatile:

$$ \tau_i = \frac{3k_i P_i}{2g}, \qquad (A.14) $$

where $k_i$ is the absorption coefficient of species i under pressure $P_i$. It is linked to the absorption coefficient at normal atmospheric conditions ($T_0, P_0$) $k_{0,i}$ through:

$$ k_i = \left( \frac{k_{0,i} g}{3 P_0} \right)^{1/2} . \qquad (A.15) $$

Eq. (A.14) can be described in terms of the atmospheric mass in volatile species i as:

$$ \tau_i = \frac{3k_i M_{\text{atm,i}}}{8\pi R_p^2}. \qquad (A.16) $$

The total opacity $\tau$ is expressed as the sum of opacities and the emissivity is then:

$$ \epsilon = \frac{2}{2 + \tau} = \frac{2}{2 + \sum_j \tau_j}. \qquad (A.17) $$

This leads to the following expression for the atmospheric flux:

$$ F_{\text{gray}} = \sigma \epsilon (T_s^4 - T_{\text{eq}}^4). \qquad (A.18) $$

In this model, we considered that because of the relative low density and viscosity of the atmosphere, convective motions in the atmosphere are efficient and the atmosphere is always at thermal equilibrium. This means that heat radiated at the top of the magma ocean is instantaneously radiated at the top of the atmosphere and the atmosphere adjusts its temperature. As a consequence, we consider that the convective heat flux (computed from Eqs. (A.11) or (A.12)) is equal to the atmospheric heat flux (Eq. (A.18)). For a given thermal state of the mantle and of the atmosphere, we numerically solve this set of equations for the unknown $T_s$, the surface temperature, using a Newton method.

Once the thermal structure of the mantle and the heat fluxes are known, the time evolution of the mantle is obtained by integration of the energy conservation equation:

$$ \frac{d}{dt} \int_{V_{\text{MOO}}} \rho T \left( C_p + \Delta H \frac{d\phi}{dT} \right) dV = -4\pi R_p^2 F_{\text{conv}}, \qquad (A.19) $$

where $V_{\text{MO}}$ is the volume of the essentially liquid magma ocean extending between the time-dependent rheological front (located at the radius $r = R_{b,i}$), and the radius of the planetary surface ($r = R_p$). $T(r)$ is the temperature, $\rho$ is the magma density, $C_p$ is the heat capacity, $\Delta H$ is the specific enthalpy difference for fusion/solidification. Since $T$ is a function of $r$, $\frac{d\phi}{dT}$ may also vary along the planetary radius. The values of each aforementioned parameter are listed in Table 3. The above equation can be expanded and simplified by replacing $dV$ with $4\pi r^2 dr$:

$$ \frac{d(\epsilon T_p)}{dt} = -R_p^2 F_{\text{conv}}, \qquad (A.20) $$

where $\epsilon$ contains the radial dependence of various quantities, and is defined as:

$$ \epsilon = \frac{1}{T_p} \int_{R_{b,i}}^{R_p} \rho T(r) \left( C_p + \Delta H \frac{d\phi(r)}{dT} \right) r^2 dr, \qquad (A.21) $$

Then, similar to previous works (e.g., Stevenson et al., 1983; Stamenković et al., 2012; Lebrun et al., 2013; Nikolaou et al., 2019) Eq. (A.20) is further simplified to:

$$ \epsilon \frac{dT_p}{dt} = -R_p^2 F_{\text{conv}}. \qquad (A.22) $$

We insist on the fact that the volume integration in Eqs. (A.19), and (A.21) above is performed only in the essentially liquid part of the mantle in contact with the surface. These expressions therefore contrast with previous parameterized models of magma ocean cooling (Lebrun
<figure><img src="image_19.png" /><figcaption>**Fig. B.15.** Time evolution of the (a) CO₂ pressure, (b), H₂O pressure, (c) emissivity, (d) surface temperature, (e) potential temperature and (f) rheological front in the case of out-of-equilibrium outgassing for an Earth-sized planet. The black dashed lines represent a case where no volatile interaction is taken into account, while the solid red lines represent cases computed with the model from Bower et al. (2019). All curves were computed for one earth mass, without rotation and with no BMO trapping. The initial content in CO₂ and H₂O were 13 ppm and 41 ppm respectively. All curves were computed with $k_{\text{O,CO}_2} = 10^{-3} \text{ m}^2/\text{kg}$.</figcaption></figure>

et al., 2013; Salvador et al., 2017; Nikolaou et al., 2019; Massol et al., 2023). These previous models assumed that the heat loss was distributed over the entire mantle, regardless of the local rheology. This hypothesis relies on the assumption that heat can be distributed evenly throughout the mantle, through porous flow within the mush, or through convective overturn in the solidified mantle (Massol et al., 2023). However, in the solidified mantle, the viscosity increase associated with the rheological transition is expected to largely hinder fluid motions. Moreover, because the Rayleigh number decreases with increasing viscosity, the solid mantle convects at considerably larger time scales than the liquid mantle. This can be seen in Fig. A.14 that compares the time evolution of the heat flux estimated at the location of the rheological transition with the heat flux computed at the surface of the magma ocean.

The bottom heat flux is always orders of magnitudes smaller than the surface heat flux, confirming that the time scale for the cooling of the magma ocean is considerably shorter than that of the cooling of the mushy mantle underneath it. Such distinct time scales make it difficult to support the claim that the heat is distributed evenly to the solid mantle. We therefore consider instead that the solid mantle cooling is considerably less efficient than that of the liquid mantle and that the heat loss only affects the liquid convecting part of the mantle. This assumption is more suitable as long as the magma ocean below

the rheological fraction exists. The influence of such assumption on the solidification time scale is significant and is discussed in Appendix C. Note also that in the above expression (Eq. (A.22)), we do not take into account heating by the decay of radioactive elements or incoming heat flux from the core-mantle boundary, which are negligible compared to the convective heat flux when a magma ocean is present (Salvador et al., 2017; Nikolaou et al., 2019).

## Appendix B. Effect of volatile interaction on the cooling time scales

Our parameterized cooling models do not take into account gas mixing effects. While such assumption was proved to be incorrect (Bower et al., 2019; Massol et al., 2023), we primarily made this simplification to more easily discuss the influence of each volatile (H₂O and CO₂) separately and/or to compare our results with previous studies (Lebrun et al., 2013; Salvador et al., 2017; Nikolaou et al., 2019). To estimate the impact of this simplification on our model results, Fig. B.15 illustrates the effect of gas mixing on the cooling history of the magma ocean in an out-of-equilibrium case, computed with the equations from Bower et al. (2019).

The results show that despite distinct evolution towards the end of the magma ocean solidification, the differences in magma ocean
<figure><img src="image_20.png" /><figcaption>Fig. C.16. Time evolution of (a) the magma ocean depth, (b) the outgoing heat flux, (c) the Rayleigh number and (d) the potential and surface temperatures. All curves represent the evolution of a terrestrial magma ocean, with black body radiation starting from $T_p = 4000K$. The black and blue curves show the evolution of a terrestrial magma ocean where Eq. (14) is integrated over the entire mantle. For the red curve, this integration is restricted to the liquid part of the magma ocean ($\phi > \phi_c = 0.4$). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

cooling history are small. Therefore, we do not expect that accounting for gas mixing effects will significantly affect our model results and conclusions.

## Appendix C. Integration of the conservation of heat

In our parameterized model, we obtain the time evolution of the mantle by integration of the energy conservation equation (Eq. (A.22)). As discussed above, throughout the existence of a fully/essentially molten magma ocean, we consider that due to the solid-like viscosity of the lower mantle, the integration should only be performed on the liquid part of the mantle. Eq. (A.22) should therefore include a second term on the left-hand side which relates to the time evolution of the considered volume. Additionally, as the magma ocean solidifies, forming crystals are being transferred from the magma ocean to the solid mantle, resulting in a heat loss from the magma ocean. Therefore, similar to Abe (1997), Eq. (A.22) should read:

$$ \epsilon \frac{dT_p}{dt} + \frac{dE_{\text{MO}}}{dt} = -R_p^2 F_{\text{conv}} + F_{\text{crystals}}, \qquad (C.1) $$

where $E_{\text{MO}} \propto C_p T_p M_{\text{MO}}$ is the heat associated with the evolving (essentially) liquid magma ocean mass, $M_{\text{MO}}$, and $F_{\text{crystals}}$ is the heat flux associated with the flux of mass out of the system. These two terms exactly compensate each other, leading to the expression for heat conservation previously written (Eq. (A.22)).

To assess the importance of the volume considered upon solving Eq. (A.22), we compared outputs of our calculations with results from Nikolaou et al. (2019) for the same conditions. Time evolutions are computed using a black-body model. Fig. C.16 shows the time evolution of the important characteristics of magma ocean for three cases. The black curves correspond to data from Nikolaou et al. (2019) and the blue curves correspond to our model with the same parameters. The red curves correspond to the same parameters as the blue curves,

except that the integration of Eq. (A.22) is performed only in the essentially liquid part of the mantle (i.e., for $r \ge R_{b,t}$).

Our results show a good agreement with results from Nikolaou et al. (2019). In particular, both models (black and blue curves) yield similar solidification times for the magma ocean, of about 2000 years. The slight differences observed are likely due to the different choices made for computing the viscosity that enters the governing dimensionless parameters (see Appendix E). The red curves however display considerably shorter solidification time scales, with a solidification time of 341 years. All three curves display the same time evolution before the depth of the magma ocean starts to decrease. Then, the red curves all evolve more rapidly than the blue and black curves, which follow each other closely. When integrating Eq. (A.22) over the volume of the essentially liquid magma ocean, as soon as the magma ocean size decreases, the heat loss is distributed over a smaller volume, resulting in faster cooling and shorter solidification time scales. This shows that the model in Nikolaou et al. (2019) (as well as other previous models Lebrun et al., 2013; Salvador et al., 2017; Massol et al., 2023) overestimate solidification time.

## Appendix D. Solidification times

See Tables D.4 and D.5.

## Appendix E. Choice of the reference viscosity

The viscosity is one of the key parameters controlling the dynamics of the magma ocean. Not only does it appear in both the governing parameters of Rayleigh-Bénard convection, the Rayleigh and Prandtl numbers, its sharp increase due to low melt content controls the rheological transition between a liquid magma ocean and a solidified mantle. The adopted definition of the viscosity therefore strongly influences the time evolution of the magma ocean. In Nikolaou et al. (2019),
<caption>Table D.4</caption>

Solidification times for cases displayed in Fig. 11. The solidification times are shown in years. The gray cells represent equilibrium cases (labeled 'eq') while white cells represent out-of-equilibrium cases (labeled 'neq'). All temporal evolution were computed with an absorption coefficient for CO₂ equal to k₀CO₂ = 10⁻⁴m²kg⁻¹.

<table><thead><tr><th rowspan="2">Equilibrium state</th><th rowspan="2"></th><th colspan="2">1</th><th colspan="2">2</th><th colspan="2">3</th><th colspan="2">4</th><th colspan="2">5</th><th colspan="2">6</th><th colspan="2">7</th><th colspan="2">8</th><th colspan="2">9</th><th colspan="2">10</th></tr><tr><th>eq</th><th>neq</th><th>eq</th><th>neq</th><th>eq</th><th>neq</th><th>eq</th><th>neq</th><th>eq</th><th>neq</th><th>eq</th><th>neq</th><th>eq</th><th>neq</th><th>eq</th><th>neq</th><th>eq</th><th>neq</th><th>eq</th><th>neq</th></tr></thead><tbody><tr><td rowspan="4">Low volatile content</td><td rowspan="2">No BMO trapping</td><td>No rotation</td><td>1002</td><td>780</td><td>1226</td><td>901</td><td>1411</td><td>998</td><td>1573</td><td>1079</td><td>1717</td><td>1152</td><td>1849</td><td>1218</td><td>1974</td><td>1279</td><td>2093</td><td>1338</td><td>2208</td><td>1395</td><td>2319</td><td>1450</td></tr><tr><td>With rotation</td><td>667</td><td>234</td><td>865</td><td>292</td><td>1022</td><td>339</td><td>1156</td><td>378</td><td>1274</td><td>413</td><td>1382</td><td>444</td><td>1483</td><td>472</td><td>1578</td><td>499</td><td>1669</td><td>524</td><td>1757</td><td>548</td></tr><tr><td rowspan="2">BMO trapping</td><td>No rotation</td><td>984</td><td>765</td><td>1130</td><td>808</td><td>1262</td><td>849</td><td>1383</td><td>890</td><td>1495</td><td>930</td><td>1600</td><td>968</td><td>1703</td><td>1008</td><td>1802</td><td>1046</td><td>1899</td><td>1086</td><td>1994</td><td>1126</td></tr><tr><td>With rotation</td><td>651</td><td>222</td><td>777</td><td>221</td><td>882</td><td>225</td><td>977</td><td>231</td><td>1064</td><td>239</td><td>1145</td><td>247</td><td>1224</td><td>257</td><td>1300</td><td>267</td><td>1374</td><td>277</td><td>1446</td><td>287</td></tr><tr><td rowspan="4">High volatile content</td><td rowspan="2">No BMO trapping</td><td>No rotation</td><td>7317</td><td>6748</td><td>9456</td><td>8787</td><td>11101</td><td>10341</td><td>12481</td><td>11633</td><td>13681</td><td>12750</td><td>14766</td><td>13754</td><td>15771</td><td>14680</td><td>16717</td><td>15548</td><td>17614</td><td>16369</td><td>18474</td><td>17155</td></tr><tr><td>With rotation</td><td>6725</td><td>3989</td><td>8821</td><td>5014</td><td>10418</td><td>5699</td><td>11752</td><td>6219</td><td>12908</td><td>6639</td><td>13952</td><td>6997</td><td>14915</td><td>7312</td><td>15820</td><td>7596</td><td>16677</td><td>7857</td><td>17498</td><td>8100</td></tr><tr><td rowspan="2">BMO trapping</td><td>No rotation</td><td>7016</td><td>6467</td><td>7990</td><td>7376</td><td>8852</td><td>8155</td><td>9654</td><td>8871</td><td>10400</td><td>9533</td><td>11106</td><td>10157</td><td>11801</td><td>10771</td><td>12464</td><td>11355</td><td>13121</td><td>11936</td><td>13769</td><td>12507</td></tr><tr><td>With rotation</td><td>6428</td><td>3704</td><td>7371</td><td>3566</td><td>8191</td><td>3442</td><td>8949</td><td>3363</td><td>9653</td><td>3312</td><td>10319</td><td>3281</td><td>10972</td><td>3277</td><td>11596</td><td>3275</td><td>12214</td><td>3293</td><td>12822</td><td>3322</td></tr></tbody></table>
<caption>Table D.5</caption>

Solidification times for cases displayed in Fig. 11. The solidification times are shown in years. The gray cells represent equilibrium cases (labeled 'eq') while white cells represent out-of-equilibrium cases (labeled 'neq'). All temporal evolution were computed with an absorption coefficient for CO₂ equal to k₀CO₂ = 5 · 10⁻² m² kg⁻¹.

<table><thead><tr><th colspan="2">Planetary mass ratio</th><th colspan="2">1</th><th colspan="2">2</th><th colspan="2">3</th><th colspan="2">4</th><th colspan="2">5</th><th colspan="2">6</th><th colspan="2">7</th><th colspan="2">8</th><th colspan="2">9</th><th colspan="2">10</th></tr><tr><th colspan="2">Equilibrium state</th><th>eq</th><th>neq</th><th>eq</th><th>neq</th><th>eq</th><th>neq</th><th>eq</th><th>neq</th><th>eq</th><th>neq</th><th>eq</th><th>neq</th><th>eq</th><th>neq</th><th>eq</th><th>neq</th><th>eq</th><th>neq</th><th>eq</th><th>neq</th></tr></thead><tbody><tr><td rowspan="4">Low volatile content</td><td>No BMO trapping</td><td>10234</td><td>8866</td><td>13930</td><td>11585</td><td>16850</td><td>13639</td><td>19340</td><td>15331</td><td>21527</td><td>16773</td><td>23518</td><td>18059</td><td>25378</td><td>19244</td><td>27136</td><td>20353</td><td>28815</td><td>21403</td><td>30426</td><td>22405</td></tr><tr><td>BMO trapping</td><td>9578</td><td>2753</td><td>13213</td><td>3581</td><td>16071</td><td>4105</td><td>18503</td><td>4486</td><td>20635</td><td>4780</td><td>22574</td><td>5024</td><td>24382</td><td>5228</td><td>26091</td><td>5405</td><td>27721</td><td>5563</td><td>29258</td><td>5702</td></tr><tr><td>No rotation</td><td>10094</td><td>8714</td><td>13173</td><td>10769</td><td>15573</td><td>12270</td><td>17638</td><td>13515</td><td>19470</td><td>14587</td><td>21157</td><td>15559</td><td>22762</td><td>16481</td><td>24288</td><td>17351</td><td>25767</td><td>18198</td><td>27202</td><td>19022</td></tr><tr><td>With rotation</td><td>9440</td><td>2409</td><td>12464</td><td>2100</td><td>14806</td><td>1874</td><td>16814</td><td>1720</td><td>18593</td><td>1611</td><td>20229</td><td>1529</td><td>21784</td><td>1476</td><td>23260</td><td>1425</td><td>24691</td><td>1394</td><td>26079</td><td>1373</td></tr><tr><td rowspan="4">High volatile content</td><td>No BMO trapping</td><td>96715</td><td>95774</td><td>133250</td><td>131900</td><td>161940</td><td>160220</td><td>186320</td><td>184260</td><td>207680</td><td>205300</td><td>227110</td><td>224400</td><td>245220</td><td>242210</td><td>262330</td><td>259010</td><td>278640</td><td>275040</td><td>294290</td><td>290390</td></tr><tr><td>BMO trapping</td><td>95633</td><td>79192</td><td>132070</td><td>103840</td><td>160660</td><td>121860</td><td>184940</td><td>136390</td><td>206220</td><td>148570</td><td>225560</td><td>159290</td><td>243590</td><td>169080</td><td>260610</td><td>178170</td><td>276850</td><td>186740</td><td>292410</td><td>194870</td></tr><tr><td>No rotation</td><td>95193</td><td>94281</td><td>125180</td><td>123880</td><td>148430</td><td>146770</td><td>168400</td><td>166380</td><td>186100</td><td>183750</td><td>202380</td><td>199710</td><td>217860</td><td>214880</td><td>232560</td><td>229290</td><td>246810</td><td>243250</td><td>260640</td><td>256790</td></tr><tr><td>With rotation</td><td>94115</td><td>77456</td><td>124010</td><td>94727</td><td>147170</td><td>106740</td><td>167040</td><td>116470</td><td>184660</td><td>124720</td><td>200860</td><td>132130</td><td>216250</td><td>139160</td><td>230870</td><td>145770</td><td>245050</td><td>152230</td><td>258800</td><td>158510</td></tr></tbody></table>
<figure><img src="image_21.png" /><figcaption>**Fig. E.17.** Time evolution of the viscosity for two different definitions. The black curve corresponds to the case where the viscosity is averaged over the liquid magma ocean depth using an harmonic mean and the blue curve corresponds to the case where the viscosity is estimated at $T = T_p$. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. F.18.** Melting curves and temperature profiles of a planet with 10 Earth masses, for three isotherms $T_p = 4000$ K, $T_p = 3300$ K and $T_p = 2500$ K. The black dotted line shows the liquidus, the black dashed line shows the solidus and the gray dashed line shows the rheological transition for a critical liquid content $\phi_c = 0.4$. The expressions used are described in Appendix A. All the curves were computed by solving simultaneously Eqs. (A.1) and (A.2) along with a lithostatic pressure profile, with $g_{\text{planet}} = g_{\text{Earth}} \left(\frac{M_*}{M_{\text{Earth}}}\right)^{0.46} = 28.29 \text{ m s}^{-2}$.</figcaption></figure>
the viscosity was computed using a volumetric harmonic mean of the viscosity profile. In this study, because the convective motions are determined by destabilization of a thin layer at the top of the magma ocean, we considered that the Rayleigh and Prandtl numbers should be computed using the viscosity at the top of the magma ocean. We show the time evolution of the magma ocean viscosity for both definitions in Fig. E.17.

While our definition yields slightly higher values of viscosity, it shows a good agreement with the definition in Nikolaou et al. (2019), with a slow and steady increase as the magma ocean cools down, before a sharp increase when the melt fraction decreases towards the critical melt fraction and the rheological front gets close to the surface. The viscosity computed using the definition in Nikolaou et al. (2019) suffers however from several discontinuities. These discontinuities stem from the fact the definition of the viscosity depends strongly on the domain considered, due to the use of the volumetric harmonic mean, yielding a non-physical evolution of the viscosity. Because these discontinuities in viscosity are only artifacts of the calculation method, we chose to use our definition of the viscosity to describe the magma ocean.

## Appendix F. Melting curves for super-earths

When applying our results to planets larger than Earth, we resort to melting curves from Stixrude (2014), which we adapt to the synthetic melting curves from Nikolaou et al. (2019), for continuity purposes. The melting curves' evolution with pressure is shown in Fig. F.18, along with temperature profiles.

Regardless of the planetary mass, the adiabat intersects the $\phi = \phi_c$ curve at a pressure of about 90 GPa. If the density profile allows it, this can result in the formation of a basal magma ocean, whose extent is therefore larger for larger planetary masses.

## Data availability

Data will be made available on request.

## References

* Abe, Y., 1997. Thermal and chemical evolution of the terrestrial magma ocean. Phys. Earth Planet. Inter. 100 (1-4), 27-39.

* Agnor, C.B., Canup, R.M., Levison, H.F., 1999. On the character and consequences of large impacts in the late stage of terrestrial planet formation. Icarus 142 (1), 219-237.

* Ahlers, G., Grossmann, S., Lohse, D., 2009. Heat transfer and large scale dynamics in turbulent Rayleigh-Bénard convection. Rev. Modern Phys. 81 (2), 503-537.

* Ballmer, M.D., Lourenço, 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. Geosystems 18 (7), 2785-2806.

* Boukaré, C.-E., Parmentier, E., Parman, S., 2018. Timing of mantle overturn during magma ocean solidification. Earth Planet. Sci. Lett. 491, 216-225.

* Bower, D.J., Kitzmann, D., Wolf, A.S., Sanan, P., Dorn, C., Oza, A.V., 2019. Linking the evolution of terrestrial interiors and an early outgassed atmosphere to astrophysical observations. Astron. Astrophys. 631, A103.

* Bower, D.J., Sanan, P., Wolf, A.S., 2018. Numerical solution of a non-linear conservation law applicable to the interior dynamics of partially molten planets. Phys. Earth Planet. Inter. 274, 49-62.

* Canup, R.M., 2008. Accretion of the Earth. Philos. Trans. R. Soc. A: Math. Phys. Eng. Sci. 366 (1883), 4061-4075.

* Castaing, B., Gunaratne, G., Heslot, F., Kadanoff, L., Libchaber, A., Thomae, S., Wu, X.-Z., Zaleski, S., Zanetti, G., 1989. Scaling of hard thermal turbulence in Rayleigh-Bénard convection. J. Fluid Mech. 204, 1-30.

* Chillà, F., Schumacher, J., 2012. New perspectives in turbulent Rayleigh-Bénard convection. Eur. Phys. J. E 35, 1-25.

* Čuk, M., Stewart, S.T., 2012. Making the Moon from a fast-spinning Earth: A giant impact followed by resonant despinning. Science 338 (6110), 1047-1052.

* Elkins-Tanton, L.T., 2008. Linked magma ocean solidification and atmospheric growth for Earth and Mars. Earth Planet. Sci. Lett. 271 (1-4), 181-191.

* Elkins-Tanton, L.T., Parmentier, E., Hess, P., 2003. Magma ocean fractional crystallization and cumulate overturn in terrestrial planets: Implications for Mars. Meteorit. Planet. Sci. 38 (12), 1753-1771.

* Fiquet, G., Auzende, A., Siebert, J., Corgne, A., Bureau, H., Ozawa, H., Garbarino, G., 2010. Melting of peridotite to 140 gigapascals. Science 329 (5998), 1516-1518.

* Foley, B.J., Driscoll, P.E., 2016. Whole planet coupling between climate, mantle, and core: Implications for rocky planet evolution. Geochem. Geophys. Geosystems 17 (5), 1885-1914.

* Gastine, T., Wicht, J., Aubert, J., 2016. Scaling regimes in spherical shell rotating convection. J. Fluid Mech. 808, 690-732.

* Grossmann, S., Lohse, D., 2000. Scaling in thermal convection: A unifying theory. J. Fluid Mech. 407, 27-56.

* Hamano, K., Abe, Y., Genda, H., 2013. Emergence of two types of terrestrial planet on solidification of magma ocean. Nature 497 (7451), 607-610.

* Hart, M.H., 1979. Habitable zones about main sequence stars. Icarus 37 (1), 351-357.

* Heber, V.S., Brooker, R.A., Kelley, S.P., Wood, B.J., 2007. Crystal-melt partitioning of noble gases (helium, neon, argon, krypton, and xenon) for olivine and clinopyroxene. Geochim. Cosmochim. Acta 71 (4), 1041-1061.

* Herzberg, C., Raterron, P., Zhang, J., 2000. New experimental observations on the anhydrous solidus for peridotite KLB-1. Geochem. Geophys. Geosystems 1 (11).

* Hier-Majumder, S., Hirschmann, M.M., 2017. The origin of volatiles in the Earth's mantle. Geochem. Geophys. Geosystems 18 (8), 3078-3092.

* Hirschmann, M.M., 2000. Mantle solidus: Experimental constraints and the effects of peridotite composition. Geochem. Geophys. Geosystems 1 (10).

* Holloway, J., Blank, J., Carroll, M., 1994. Volatiles in magmas. Rev. Miner. 30.

* Ida, S., Lin, D.N., 2004. Toward a deterministic model of planetary formation. I. A desert in the mass and semimajor axis distributions of extrasolar planets. Astrophys. J. 604 (1), 388.

* Karki, B.B., Stixrude, L.P., 2010. Viscosity of MgSiO₃ liquid at Earth's mantle conditions: Implications for an early magma ocean. Science 328 (5979), 740-742.

* Kasting, J.F., Whitmire, D.P., Reynolds, R.T., 1993. Habitable zones around main sequence stars. Icarus 101 (1), 108-128.

* Kokubo, E., Genda, H., 2010. Formation of terrestrial planets from protoplanets under a realistic accretion condition. Astrophys. J. Lett. 714 (1), L21.

* Kokubo, E., Ida, S., 2007. Formation of terrestrial planets from protoplanets. II. Statistics of planetary spin. Astrophys. J. 671 (2), 2082.

* Labrosse, S., Hernlund, J., Coltice, N., 2007. A crystallizing dense magma ocean at the base of the Earth's mantle. Nature 450 (7171), 866-869.

* Lebrun, T., Massol, H., Chassefière, E., Davaille, A., Marcq, E., Sarda, P., Leblanc, F., Brandeis, G., 2013. Thermal evolution of an early magma ocean in interaction with the atmosphere. J. Geophys. Res. Planets 118 (6), 1155-1176.

* Lejeune, A.-M., Richet, P., 1995. Rheology of crystal-bearing silicate melts: An experimental study at high viscosities. J. Geophys. Res. Solid Earth 100 (B3), 4215-4229.

* Lichtenberg, T., Bower, D.J., Hammond, M., Boukrouche, R., Sanan, P., Tsai, S.-M., Pierrehumbert, R.T., 2021. Vertically resolved magma ocean-protoatmosphere evolution: H₂, H₂O, CO₂, CH₄, CO, O₂, and N₂ as primary absorbers. J. Geophys. Res. Planets 126 (2), e2020JE006711.

* Maas, C., Hansen, U., 2015. Effects of Earth's rotation on the early differentiation of a terrestrial magma ocean. J. Geophys. Res. Solid Earth 120 (11), 7508-7525.

* Maas, C., Hansen, U., 2019. Dynamics of a terrestrial magma ocean under planetary rotation: A study in spherical geometry. Earth Planet. Sci. Lett. 513, 81-94.

* Massol, H., Davaille, A., Sarda, P., 2023. Early formation of a water ocean as a function of initial CO₂ and H₂O contents in a solidifying rocky planet. J. Geophys. Res. Planets 128 (8), e2023JE007848.

* Maurice, M., Tosi, N., Samuel, H., Plesa, A.-C., Hüttig, C., Breuer, D., 2017. Onset of solid-state mantle convection and mixing during magma ocean solidification. J. Geophys. Res. Planets 122 (3), 577-598.

* Miyazaki, Y., Korenaga, J., 2022. Inefficient water degassing inhibits ocean formation on rocky planets: An insight from self-consistent mantle degassing models. Astrobiology 22 (6), 713-734.

* Moeller, A., Hansen, U., 2013. Influence of rotation on the metal rain in a hadean magma ocean. Geochem. Geophys. Geosystems 14 (4), 1226-1244.

* Monteux, J., Andrault, D., Samuel, H., 2016. On the cooling of a deep terrestrial magma ocean. Earth Planet. Sci. Lett. 448, 140-149.

* Monteux, J., Ricard, Y., Coltice, N., Dubuffet, F., Ulvrova, M., 2009. A model of metal-silicate separation on growing planets. Earth Planet. Sci. Lett. 287 (3-4), 353-362.

* Nikolaou, A., Katyal, N., Tosi, N., Godolt, M., Grenfell, J.L., Rauer, H., 2019. What factors affect the duration and outgassing of the terrestrial magma ocean? Astrophys. J. 875 (1), 11.

* Ozgurel, O., Caracas, R., 2023. The magma ocean was a huge helium reservoir in the early Earth. Geochem. Perspect. Lett. 25, 46-50.

* Pan, V., Holloway, J.R., Hervig, R.L., 1991. The pressure and temperature dependence of carbon dioxide solubility in tholeiitic basalt melts. Geochim. Cosmochim. Acta 55 (6), 1587-1595. http://dx.doi.org/10.1016/0016-7037(91)90130-W.

* van der Poel, E.P., Stevens, R.J., Lohse, D., 2011. Connecting flow structures and heat flux in turbulent Rayleigh-Bénard convection. Phys. Rev. E-Stat. Nonlinear, Soft Matter Phys. 84 (4), 045303.

* Regenauer-Lieb, K., Yuen, D.A., Branlund, J., 2001. The initiation of subduction: Criticality by addition of water? Science 294 (5542), 578-580.

* Roscoe, R., 1952. The viscosity of suspensions of rigid spheres. Br. J. Appl. Phys. 3 (8), 267.
Rubie, D.C., Jacobson, S.A., Morbidelli, A., O'Brien, D.P., Young, E.D., de Vries, J., Nimmo, F., Palme, H., Frost, D.J., 2015. Accretion and differentiation of the terrestrial planets with implications for the compositions of early-formed solar system bodies and accretion of water. *Icarus* 248, 89–108.

Salvador, A., Massol, H., Davaille, A., Marcq, E., Sarda, P., Chassefière, E., 2017. The relative influence of H₂O and CO₂ on the primitive surface conditions and evolution of rocky planets. *J. Geophys. Res. Planets* 122 (7), 1458–1486.

Salvador, A., Samuel, H., 2023. Convective outgassing efficiency in planetary magma oceans: Insights from computational fluid dynamics. *Icarus* 390, 115265.

Samuel, H., 2012. A re-evaluation of metal diapir breakup and equilibration in terrestrial magma oceans. *Earth Planet. Sci. Lett.* 313, 105–114.

Samuel, H., Drilleau, M., Rivoldini, A., Xu, Z., Huang, Q., Garcia, R.F., Lekić, V., Irving, J.C., Badro, J., Lognonné, P.H., et al., 2023. Geophysical evidence for an enriched molten silicate layer above Mars's core. *Nature* 622 (7984), 712–717.

Samuel, H., Evonuk, M., 2010. Modeling advection in geophysical flows with particle level sets. *Geochem. Geophys. Geosystems* 11 (8).

Samuel, H., Tackley, P., Evonuk, M., 2010. Heat partitioning in terrestrial planets during core formation by negative diapirism. *Earth Planet. Sci. Lett.* 290 (1–2), 13–19.

Schanz, D., Gesemann, S., Schröder, A., 2016. Shake-the-box: Lagrangian particle tracking at high particle image densities. *Exp. Fluids* 57, 1–27.

Schanz, D., Schröder, A., Gesemann, S., Michaelis, D., Wieneke, B., 2013. 'Shake the box': A highly efficient and accurate tomographic particle tracking velocimetry (TOMO-PTV) method using prediction of particle positions. In: PIV13; 10th International Symposium on Particle Image Velocimetry, Delft, the Netherlands, July 1–3, 2013.

Schmalzl, J., Breuer, M., Hansen, U., 2002. The influence of the Prandtl number on the style of vigorous thermal convection. *Geophys. Astrophys. Fluid Dyn.* 96 (5), 381–403.

Schumacher, J., Bodenschatz, E., 2012. Focus on new perspectives in high-Rayleigh-number turbulent convection. *New J. Phys.* 14 (9), 095004.

Selsis, F., Leconte, J., Turbet, M., Chaverot, G., Bolmont, É., 2023. A cool runaway greenhouse without surface magma ocean. *Nature* 620 (7973), 287–291.

Senshu, H., Kuramoto, K., Matsui, T., 2002. Thermal evolution of a growing Mars. *J. Geophys. Res. Planets* 107 (E12), 1.

Shraiman, B.I., Siggia, E.D., 1990. Heat transport in high-Rayleigh-number convection. *Phys. Rev. A* 42 (6), 3650.

Sim, S.J., Hirschmann, M.M., Hier-Majumder, S., 2024. Volatile and trace element storage in a crystallizing martian magma ocean. *J. Geophys. Res. Planets* 129 (8), e2024JE008346.

Solomatov, V., 2015. Magma oceans and primordial mantle differentiation. In: Treatise on Geophysics, second ed. pp. 81–104.

Sossi, P.A., Tollan, P.M., Badro, J., Bower, D.J., 2023. Solubility of water in peridotite liquids and the prevalence of steam atmospheres on rocky planets. *Earth Planet. Sci. Lett.* 601, 117894.

Sparks, R., Barclay, J., Jaupart, C., Mader, H., Phillips, J., Sturtevant, B., 1994. Physical aspects of magma degassing. *Miner. Soc. Am. Rev. Miner.* 30, 413–445.

Stamenković, V., Noack, L., Breuer, D., Spohn, T., 2012. The influence of pressure-dependent viscosity on the thermal evolution of super-earths. *Astrophys. J.* 748 (1), 22–41. http://dx.doi.org/10.1088/0004-637X/748/1/41.

Stevens, R.J., Hartmann, R., Verzicco, R., Lohse, D., 2024. How wide must Rayleigh-bénard cells be to prevent finite aspect ratio effects in turbulent flow? *J. Fluid Mech.* 1000, A58.

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

Stixrude, L., 2014. Melting in super-earths. *Philos. Trans. R. Soc. A: Math. Phys. Eng. Sci.* 372 (2014), 20130076.

Thomas, B., Samuel, H., Farnetani, C., Aubert, J., Chauvel, C., 2024. Mixing time of heterogeneities in a buoyancy-dominated magma ocean. *Geophys. J. Int.* 236 (2), 764–777.

Tonks, W.B., Melosh, H.J., 1993. Magma ocean formation due to giant impacts. *J. Geophys. Res. Planets* 98 (E3), 5319–5333.

Valencia, D., O'Connell, R.J., Sasselov, D., 2006. Internal structure of massive terrestrial planets. *Icarus* 181 (2), 545–554.

Verzicco, R., Camussi, R., 1999. Prandtl number effects in convective turbulence. *J. Fluid Mech.* 383, 55–73.

Wagner, S., Shishkina, O., 2013. Aspect-ratio dependency of Rayleigh-Bénard convection in box-shaped containers. *Phys. Fluids* 25 (8).

Walbecq, A., Samuel, H., Limare, A., 2024. Fully determined three-dimensional velocity field in a divergence-free convection experiment with rigid boundary conditions. *Exp. Fluids* 65 (5), 1–17.

Weis, D., Harpp, K.S., Harrison, L.N., Boyet, M., Chauvel, C., Farnetani, C.G., Finlayson, V.A., Lee, K.K.M., Parai, R., Shahar, A., Williamson, N.M.B., 2023. Earth's mantle composition revealed by mantle plumes. *Nat. Rev. Earth Env.* 4, 604–625. http://dx.doi.org/10.1038/s43017-023-00467-0.

Zhang, J., Herzberg, C., 1994. Melting experiments on anhydrous peridotite KLB-1 from 5.0 to 22.5 GPa. *J. Geophys. Res. Solid Earth* 99 (B9), 17729–17742.