{'fraction_non_alphanumeric': 0.035731096159703185, 'fraction_numerical': 0.020752000443004678, 'mean_word_length': 4.936472715318869, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Convective outgassing efficiency in planetary magma oceans: insights from computational fluid dynamics

Planetary atmospheres are commonly thought to result from the efficient outgassing of cooling magma oceans. During this stage, vigorous convective motions in the molten interior are believed to rapidly transport the dissolved volatiles to shallow depths where they exsolve and burst at the surface. This assumption of efficient degassing and atmosphere formation has important implications for planetary evolution, but has never been tested against fluid dynamics considerations. Yet, during a convective cycle, only a finite fraction of the magma ocean can reach the shallow depths where volatiles exsolution can occur, and a large-scale circulation may prevent a substantial magma ocean volume from rapidly reaching the planetary surface. Therefore, we conducted computational fluid dynamics experiments of vigorous 2D and 3D Rayleigh-B\'enard convection at Prandtl number of unity to characterize the ability of the convecting fluid to reach shallow depths at which volatiles are exsolved and extracted to the atmosphere. Outgassing efficiency is essentially a function of the magnitude of the convective velocities. This allows deriving simple expressions to predict the time evolution of the amount of outgassed volatiles as a function of the magma ocean governing parameters. For plausible cases, the time required to exsolve all oversaturated water can exceed the magma ocean lifetime in a given highly vigorous transient stage, leading to incomplete or even negligible outgassing. Furthermore, the planet size and the initial magma ocean water content, through the convective vigor and the exsolution depth, respectively, strongly affect magma oceans degassing efficiency, possibly leading to divergent planetary evolution paths and resulting surface conditions. Overall, despite vigorous convection, for a significant range of parameters, convective degassing appears not as efficient as previously thought.


Introduction
By controlling both the surface temperature and pressure, atmospheres play a fundamental role in planetary evolution and habitability (defined as the ability of a planet to sustain liquid water at its surface; e.g., Kasting et al., 1993).
The formation of atmospheres is thought to have occurred early in planetary history, through the outgassing of the molten interior, during the magma ocean phase. During this period, accretional heating (e.g., Safronov, 1978;Kaula, 1979;Tonks and Melosh, 1993;Nakajima et al., 2021), radiogenic heating by the decay of short-lived elements (e.g., Merk et al., 2002;Bhatia, 2021), and the conversion of gravitational potential energy into heat via viscous dissipation during core formation (e.g., Sasaki and Nakazawa, 1986;, likely lead rocky bodies to experience at least one episode of global mantle melting (Wood et al., 1970;Stevenson, 1981;Elkins-Tanton, 2012).
In the classical scenario for magma ocean progressive cooling and crystallization (e.g., Solomatov, 2015) the magma ocean undergoes different major evolution stages that are essentially governed by the rheology of the meltcrystal mixture and by the relationships between the evolving geotherm and mantle melting curves (Figure 1a), leading to changes in the dynamics and in the thickness L of the fully (φ " 1, where φ is the melt fraction) or essentially (φ c ă φ ă 1, where φ c is a critical melt fraction at which the rheology of the partially molten mixture becomes dominated by that of the solids) liquid magma ocean layer: (i) a global fully/essentially molten stage (Figure 1b), where the entire mantle thickness is either fully molten (because the temperature is hotter than the liquidus at all magma ocean depths) or is partially molten (crystals form and accumulate at the base of the mantle where the geotherm intersects the liquidus first) but has a melt fraction that exceeds the critical value φ c « 0.4 at all depths (Lejeune and Richet, 1995;Abe, 1997). In this case, the partially molten mixture behaves as a liquid throughout the entire mantle thickness, i.e.,  Figure 1: (a) Viscosity of the silicate magma ocean as a function of its evolving temperature/crystal fraction (modified from Lejeune and Richet, 1995;Salvador et al., 2017). (b-d) Schematic representation of the succession of three main stages during the progressive solidification of the magma ocean, from (a) a fully (φ " 1)/essentially (φc ă φ ă 1) molten stage of global extent, to (b) a mixed essentially molten upper unit coexisting with a lower mushy unit, to (c) a global mushy stage. Red and yellow areas correspond to regions whose rheologies are dominated by the melt (φ ą φc) and by the crystals (φ ă φc), respectively. The rheological transition and associated viscosity jump occur at the critical liquid fraction φc (black dashed lines in (a) and (c)). During the magma ocean bottom-up solidification, from top to bottom in (b-d), the convective vigor decreases and the inertial effects increase as the viscosity increases and the temperature decreases.
become interconnected, leading to a rheology of the partially molten mixture dominated by the solids. This transition yields an abrupt increase in viscosity ( Figure 1a). The appearance of this evolving rheological front (where φ " φ c ) marks the separation between two units that convect on very distinct time scales ( Figure 1c): the lower, deep mushy unit characterized by a relatively high crystal fraction and a considerably higher viscosity, and the upper, essentially liquid, unit with relatively low crystal fraction, whose rheology remains dominated by that of the melt (L ă R planet´Rcore ). The viscosity contrast between these two units imply a considerably higher convective vigor with inertial effects in the upper unit, which is the one of interest here. However, as the rheological front progresses upwards, the thickness L of the upper essentially liquid unit decreases, leading to a progressive decrease of convective vigor. (iii) When the rheological front reaches the surface, the entire mantle thickness enters a global mushy stage (Figure 1d) characterized by a relatively high viscosity (dominated by the rheology of the interconnected solids) with negligible inertial effects and convective vigor compared to that of the earliest magma ocean stage (Abe, 1993;Lejeune and Richet, 1995).
The stages described above and depicted in Figure 1 assume a bottom-up magma ocean crystallization (e.g., Andrault et al., 2011) but the reasoning would remain the same if the solidification initiates at the mid-mantle as proposed for example in Labrosse et al. (2007) and focus on the region where the solidification front propagates upward. The distribution of chemical elements would however be considerably affected in the presence of several isolated melt reservoirs.
These major stages are composed of a succession of smaller "transient stages", where the motions of the convecting fluid and its properties gradually evolve from a melt-to a crystal-dominated rheology. The entire magma ocean solidification sequence is thus made out of the succession of all these transient stages. The latter correspond to the time intervals during which magma ocean degassing efficiency is evaluated (analogous to the computational time steps used to discretize the magma ocean evolution in numerical simulations).
The duration of each major stage varies according to several more or less constrained quantities such as the initial thermo-chemical state, the planet size, the proximity to the star, the type of orbited star, or the initial amount of volatiles and their influence on key quantities, such as temperature, density, and viscosity. The first part of the earliest stage (i.e., the fully molten stage), is the shortest in duration (10-500 years; e.g., Solomatov, 2000;Lebrun et al., 2013;Salvador et al., 2017;Nikolaou et al., 2019;Bower et al., 2019) but the most vigorous one. Despite its short lifetime, this stage remains crucial from the point of view of heat transfer and volatiles distribution because the vigor of convective motions favors efficient mixing throughout the entire molten mantle.
Even though non-negligible outgassing can continue after the magma ocean stage (Gillmann and Tackley, 2014;Gaillard and Scaillet, 2014;Tosi et al., 2017;Gillmann et al., 2020;Ortenzi et al., 2020;Gaillard et al., 2021), the absence of large-scale melting at the surface and the presence of a stiff, impermeable, and insulating lid on top of the mantle (stagnant-lid or tectonic plates), or even large amounts of water on top of the solid surface of water worlds, limit direct exchanges between the interior and the atmosphere and considerably reduce the outgassing rate and the percolation of volatiles towards the surface (Kite et al., 2009;Noack et al., 2017;Dorn et al., 2018;Kite and Ford, 2018;Krissansen-Totton et al., 2021b). Consequently, planetary outgassing during the earliest magma ocean stages essentially sets up the initial conditions from which rocky planets and their atmospheres evolve until the present day. These early stages of planetary formation may entirely determine the habitability of rocky planets (e.g., Kasting, 1988;Zahnle et al., 2007;Hamano et al., 2013;Foley and Driscoll, 2016;Salvador et al., 2017;Driscoll, 2018;Krissansen-Totton et al., 2021a;Korenaga, 2021, 2022).
Volatile materials are essentially delivered to a planet during its accretionary phase (e.g., Morbidelli et al., 2012;Marty, 2012;O'Brien et al., 2018;Venturini et al., 2020) and are first contained within the molten silicate mantle before being outgassed from the interior to progressively build up an atmosphere according to the following sequence. (i) Volatile species initially dissolved within Solidified mantle φ < φ c Figure 2: Schematic representation of the magma ocean degassing process and associated velocities. v sol is the upward solidification velocity associated with the ascent of the rheological front R rheol (black dashed line) where φ " φc (assuming that solidification occurs from bottom-up), vconv is the liquid-state convection velocity, vascent is the ascent velocity of the gas bubbles rising towards the surface. d exsol is the exsolution depth (i.e., the dissolved volatile saturation limit) obtained from solubility laws and above which oversaturated volatiles are exsolved out of the melt and start forming gas bubbles. The rheological front separates the low viscosity molten mantle with zero or relatively small crystal fraction (red) and the underlying mushy mantle with solid-dominated rheology (yellow). The thickness L of the vigorously convective magma ocean extends from the surface of the planet down to the rheological transition front.

Molten mantle
the molten mantle are transported by convective motion (black arrows in Figure 2). (ii) As the magma travels upward towards the planetary surface, the decrease in the lithostatic pressure allows the volatile species to exsolve out of the melt and form gas bubbles ( Figure 2). The depth at which volatiles are exsolved, d exsol (blue dashed line in Figure 2), corresponds to the saturation limit of the melt and is given by the solubility of the volatile species (i.e., the maximum amount of volatiles that can be dissolved into a silicate melt) (Carroll and Holloway, 1994;Papale, 1997;Berlo et al., 2011). (iii) The gas bubbles subsequently rise up through the upper molten mantle (blue arrows in Figure 2) and eventually burst at the surface (e.g., Patočka et al., 2020), thus progressively forming the atmosphere (e.g., Massol et al., 2016;Ikoma et al., 2018). The solubility of the volatile species mainly depends on the pressure and on the volatile concentration but is also a function of melt composition, the type of volatile species mixtures, and temperature (e.g., Pan et al., 1991;Moore et al., 1998;Papale, 1999;Wallace and Anderson Jr., 1999;Keppler and Bolfan-Casanova, 2006;Berlo et al., 2011). The outgassing process continues as the melt becomes more and more enriched in incompatible volatile species (i.e., the chemical species that preferentially partition into the melt phase) with the gradual solidification of the magma ocean (Elkins-Tanton, 2008;Lebrun et al., 2013). The sequence described above implicitly refers to secondary (formed) atmospheres rather than atmospheres resulting from impact degassing (e.g., Lange and Ahrens, 1982;Matsui and Abe, 1986a,b;Zahnle et al., 1988;Hashimoto et al., 2007;Schaefer and Fegley, 2010;Sakuraba et al., 2019;Zahnle et al., 2020). Conversely, primary or primordial atmospheres are composed of the solar nebula gases gravitationally captured by the planet (Hayashi et al., 1979;Sasaki and Nakazawa, 1990;Ikoma and Genda, 2006;Lammer et al., 2014;Olson and Sharp, 2019;Lammer et al., 2020). Primary atmospheres are mainly composed of light elements, i.e., H 2 and He, that are sensitive to escape processes such as impact erosion or solar extreme-ultraviolet-driven planetary wind, and these atmospheres would have dissipated when the solar nebula disappeared (e.g., Zahnle et al., 2007Zahnle et al., , 2010Schlichting et al., 2015;Schlichting 8 and Mukhopadhyay, 2018;Sakuraba et al., 2019;Lammer et al., 2020).
So far, magma ocean degassing, and thus atmosphere formation, has been considered efficient and quasi-instantaneous. This hypothesis is primarily based on one a priori robust -yet debatable-assumption that has not yet been fully tested by means of fluid dynamics experiments. It stems from the fact that a magma ocean with a large thickness and very low viscosity convects very vigorously (Figure 1), leading to large convective velocities (v conv , black arrows in Figure 2) and therefore short transit times. The solidification front velocity of the magma ocean (v sol , yellow arrow in Figure 2) is expected to be significantly smaller than the rapid convective velocities which would then carry volatiles-supersaturated liquids to pressures small enough such that volatiles above melt saturation exsolve into gas bubbles much faster than the magma ocean solidification time (i.e., v sol ! v conv ). Bubbles would then rise up to the surface fast enough that they would burst into the atmosphere before being carried back down by convective motions (i.e., v conv ! v ascent , with v ascent being the bubbles ascent velocity, blue arrows in Figure 2) (Elkins-Tanton, 2008;Ikoma et al., 2018). The assumption of rapid ascent of gas bubbles seems reasonable for magma oceans, as high temperatures, low viscosity magmas allow gas to easily segregate from the melt in the form of small bubbles that can easily merge leading to larger bubbles with a faster ascending speed (Sparks, 1978;Sparks et al., 1994;Lesher and Spera, 2015). Consequently, most of the magma ocean would have "seen" the surface at least once and lost its volatiles in excess of saturation on that occasion. Based on estimates of convective velocities only, complete circulation of the magma ocean up to the surface, and hence exsolution and degassing of all oversaturated volatiles is thought to occur in one to three weeks (e.g., Elkins-Tanton, 2008).
However, for a vigorously convecting fluid in the presence of inertial effects (the case of a magma ocean), the fluid motions can be organized according to large-scale circulations (e.g., Krishnamurti and Howard, 1981;Castaing et al., 1989;Lohse and Grossmann, 1992;Siggia, 1994), which can isolate a significant part of the magma ocean from the surface, thereby limiting the amount 9 of volatile-supersaturated melt reaching the exsolution depth, and restricting outgassing to the magma ocean fraction reaching shallow depths only.
The implicit assumption of the complete and instantaneous degassing of oversaturated volatile species in a convecting magma ocean has strong implications on the study of planetary atmosphere formation, and on the coupled interior-atmosphere cooling, solidification, and chemical differentiation of magma oceans. For example, coupled magma ocean-atmosphere numerical models implicitly assume instantaneous degassing of oversaturated volatile species in excess of saturation at all times. This amounts to computing the transfer of volatile species between the planetary interior and the atmosphere via the equilibrium between their volatile content in the gas phase from either side of the planetary surface at each computational time step (Zahnle et al., 1988;Elkins-Tanton, 2008;Hamano et al., 2013;Lebrun et al., 2013;Hamano et al., 2015;Schaefer et al., 2016;Salvador et al., 2017;Nikolaou et al., 2019;Bower et al., 2019;Lichtenberg et al., 2021;Barth et al., 2021;Krissansen-Totton et al., 2021a), such that degassing is assumed to be "instantaneous" at each modeling time step. Although the assumption of efficient outgassing in vigorously convecting fluids relying on the large value of the convective velocities may seem sound in the frame of magma ocean modeling, the degassing efficiency in such a context has never been tested against fluid dynamics considerations. This is particularly important because the early atmospheres and planetary thermo-chemical state exert a considerable influence on the subsequent evolution of rocky bodies until the present day.
Addressing the fundamental question of magma ocean degassing efficiency therefore requires to constrain the amount of melt (containing dissolved volatile species) that reaches the volatile exsolution depth as a function of time, while considering the convective patterns taking place during magma ocean cooling. To this end, we conducted computational fluid dynamics experiments of Rayleigh-Bénard convection accounting for inertial effects at a finite Prandtl number value of unity. These experiments can characterize the magma ocean outgassing efficiency as a function of the parameters governing its dynamics.
Then, the resulting description of outgassing efficiency allows one to predict the amount of oversaturated volatiles outgassed during a magma ocean stage as a function of time.
The paper is organized as follows: Section 2 provides details of the model and numerical procedure. Section 3 presents the results of the numerical experiments and relates the observed convective dynamics to the amount of fluid reaching a given depth. This allows estimating the time required for the entire fluid volume to reach arbitrary exsolution depths required for the oversaturated volatile species to fully exsolve and outgas. In Section 4 prior to the conclusion we apply our model predictions to rocky planet magma oceans in a vigorously convecting stage to compare the water depletion timescales and follow the time evolution of magma ocean degassing. To this end, we consider different initial water contents and associated realistic exsolution depths, along with two different planet sizes: an Earth-sized planet and a five Earth mass super-Earth, where we account for both an enhanced and a weak degassing limit scenarios in 200 addition to a more conservative magma ocean parameters setup.

Methods
To model the dynamics of the magma ocean, we performed numerical experiments using StreamV, a finite-volume code that models the evolution of a convecting fluid by solving the Navier-Stokes and energy equations (Samuel, 2012a), under the Boussinesq approximation. These equations written below in non-dimensional form, are the conservation of mass: the conservation of momentum: 11 and the conservation of internal energy: where u is the velocity vector, t is time, p is the dynamic pressure, T is the potential temperature, e g is a vertical unit vector pointing upwards, and DX{Dt " BX{Bt`u¨∇X is the Lagrangian derivative for the scalar or vector quantity X.
The characteristic scales used to obtain the above non-dimensional equations are the magma ocean thickness L for distance, the diffusion scale L 2 {κ for time (with κ the thermal diffusivity), κ{L for velocities, the superadiabatic temperature difference between the base and the top of the magma ocean ∆T for temperatures, and the viscous scale ηκ{L 2 for pressure and stresses.
As seen in Equation (2), the convective dynamics is entirely governed by two dimensionless numbers. One is the Prandtl number, P r, which expresses the ratio of momentum to heat diffusion: where ν " η{ρ is the kinematic viscosity, η is the dynamic viscosity, and ρ is the fluid density. The second governing dimensionless parameter is the thermal Rayleigh number, Ra, that compares the thermal buoyancy force that drives convection to the resisting effects of thermal and momentum diffusion, and therefore represents a measure of the convective vigor: where g is the gravitational acceleration, and α is the thermal expansion coefficient. Note that all the physical parameters that enter the P r and Ra numbers are assumed to be homogeneous, i.e., no temperature or pressure dependence is considered.
The Rayleigh number for a terrestrial magma ocean can be as large as 10 31 (Solomatov, 2000). Such high value is out of reach for our numerical experiments, because they would imply resolving extremely thin boundary layers that would require a prohibitive associated amount of computational time. However, our calculations span Rayleigh numbers from 10 8 to 10 12 and therefore reach the highly turbulent Rayleigh-Bénard convection regime (e.g., Siggia, 1994;Grossmann and Lohse, 2000;Chillà and Schumacher, 2012), which allows for reasonable extrapolations of our results to higher Rayleigh number values.
Contrary to the current solid-state planetary mantles (Schubert et al., 2001) where inertial effects are negligible compared to viscous effects (i.e., the infinite Prandtl number approximation), inertia cannot be neglected in the frame of a vigorously convective magma ocean. We therefore considered a constant Prandtl number P r " 1 in our numerical experiments (consistent with experimental constraints on the small viscosity of silicate melts; e.g., Urbain et al., 1982) and varied systematically the Rayleigh number.
Most of our experiments were conducted in a 2D (x, z) Cartesian domain, with a few cases in 3D Cartesian geometry (x, y, z), where the z-axis is parallel to the gravity vector and oriented upward. The top and bottom horizontal boundaries of the domain are isothermal, with dimensionless temperatures T set to 0 and 1, respectively, corresponding to Rayleigh-Bénard convection configuration where the fluid is heated from below and cooled from above, which is relevant for magma oceans. Along the vertical sidewalls and at the surface, we considered free-slip boundary conditions, while the bottom boundary is rigid (i.e., no-slip velocity boundary conditions), which accounts for the location where the silicates are solid (or above the rheological threshold at which the solid phase dominates the rheology; e.g., Lejeune and Richet, 1995).
The domain is discretized using a uniform grid spacing that is fine enough to properly resolve the thermal and viscous boundary layers and small-scale con- In 2D geometry, using a pure stream function formulation defined as pBψ{Bz,´Bψ{Bxq " pu x , u z q, the mass and momentum Equations (1) and (2) are recast into a single equation involving a general biharmonic operator (e.g., Kupferman, 2002). The algebraic system resulting from the discretization of this Equation is solved using a geometric multigrid method (Brandt, 1982). This implementation has been successfully benchmarked against various analytical and numerical solutions (Samuel, 2012a(Samuel, ,b, 2014Tosi et al., 2015;Samuel, 2018), including high Reynolds numbers, turbulent flows (Ghia et al., 1982). For 3D domains, StreamV relies on primitive variable formulation of the Navier-Stokes equations. The latter are solved using a second-order in time projection method (Chorin, 1968) on a Staggered grid. This decomposition yields a Poisson equation for the dynamic pressure, which is solved using a parallel geometric multigrid method (Brandt, 1982) relying on V-cycles with Jacobi relaxations. The gradient of the obtained pressure is then used to correct the velocity to satisfy the incompressibility constraint (Equation (1)). Even though the benchmark of the approach in the limit of infinite Prandtl number (which relies on another numerical strategy) is documented elsewhere , the StreamV implementation in the finite Prandtl number case has not been documented in 3D geometry, which is therefore shown in Appendix C following a standard convection benchmark (Fusegi et al., 1991).
For each case run, we first let the convecting system evolve until statistical steady-state, which ensures that the system has forgotten its initial condition.
We describe below the convective dynamics and detail the characteristics of the dynamical regime and associated scalings. Then, in the next Section, we relate our experimental results to the amount of fluid reaching any given depth to make quantitative predictions for the efficiency of magma ocean outgassing. with the Rayleigh number through the development of numerous small-scale vortices ( Figure 3). Therefore, the downwelling and upwelling structures responsible for the large-scale circulation switch from major extended plumes at low Ra numbers (Figures 3a and 3b) to small-scale vortices detached from the boundary layers at high Ra values (Figures 3c-e). These small-scale vortices are also present all over the domain at high Rayleigh numbers, as a result of the increased turbulence. Such features are consistent with the "hard" turbulent convective regime (Heslot et al., 1987), characterized by the existence of a large-scale coherent circulation within the turbulent fluid. This is further confirmed by the exponential probability density function in the temperature distribution observed in our experiments, which is another characteristic of this regime. Features associated to this convection regime and emphasized here have been extensively described in the literature (e.g., Castaing et al., 1989;Yuen, 1999, 2000;Rogers et al., 2003).
We derived N u´Ra and Re´Ra scaling relationships, linking respectively the dimensionless Nusselt number, N u, which measures the convective efficiency (transport of heat) relative to diffusion (if N u " 1, the heat flux is purely conductive), and the Reynolds number, Re, which is the ratio of inertial forces to viscous forces (associated with momentum transport), to the Rayleigh number.   (Grossmann and Lohse, 2000, referred to as G00) and numerical models (Boussinesq model of Rogers et al., 2003, R03) are plotted in blue and pink, respectively. The N u´Ra scaling usually considered in magma ocean studies is shown in green (e.g., Siggia, 1994;Solomatov, 2015).
Re " a Re Ra γ Re .
The dependence between these parameters relies on the still controversial values of both the exponents γ and pre-factors a (e.g., Malkus, 1954;Grossmann and Lohse, 2000;Chillà and Schumacher, 2012;Stevens et al., 2013Stevens et al., , 2018. For the calculations conducted here, we found a N u " 0.51 and γ N u " 0.25 ( Figure 4a, black line), using the data points shown in Figure 4a (black crosses), which correspond to time-averaged values over the 100 convective transits achieved during the second step of our numerical procedure (i.e., after statistical steady-state has been reached). Because N u fluctuates with time around a value at statistical steady-state, the associated range for the power exponent is γ N u " 0.21´0.32 (data points and corresponding uncertainty bars in Figure 4a). This is in good agreement with the analytic γ N u " 1{4 scaling (regime I l in Grossmann and Lohse, 2000, blue line in Figure 4a), and encompasses the classical experimentally-derived N u " Ra 2{7 scaling (Castaing et al., 1989;Wu and Libchaber, 1992) as well as outcomes of numerical simulations (Boussinesq model of Rogers et al., 2003, pink line in Figure 4a), that are closer to the canonical Ra 1{3 scaling (green line in Figure 4a, also corresponding to regime IV u in Grossmann and Lohse, 2000) often considered in magma ocean parameterized convection studies (e.g., Lebrun et al., 2013;Solomatov, 2015). The upper bound N u " Ra 1{3 scaling would induce significantly larger outgoing heat fluxes at the top of the cell when extrapolating to magma oceans relevant Rayleigh numbers (up to Ra " 10 31 ), thus implying a more efficient heat transport and cooling. Yet, an exponent lower than γ N u " 1{3 has been suggested for largescale circulations (Lebrun et al., 2013) and would thus be in agreement with the convective patterns observed here ( Figure 3). For the Reynolds-Rayleigh scaling (Equation (7)), the power-law relationship obtained here takes a Re " 1.32 and (1997) and with the numerical simulations of Rogers et al. (2003), where γ Re is found to be γ Re " 0.49 for their Boussinesq model ( Figure 4b, pink line).
These studies agree with the regime I l scaling of Grossmann and Lohse (2000) (γ Re " 1{2, blue line in Figure 4b), but the Re " Ra 4{9 scaling of their IV u regime, associated with higher Rayleigh number values, is even closer to our results. Note that both N u´Ra and Re´Ra scalings are sensitive to boundary conditions, aspect ratio, and values of Rayleigh and Prandtl numbers considered (e.g., Ahlers et al., 2009;Chillà and Schumacher, 2012;Roche, 2020). In our experiments, the velocity and temperature boundary conditions are different from the aforementioned previous studies that consider either rigid, free-slip or periodic boundary conditions for the horizontal and vertical surfaces. Similarly, the aspect ratio and the geometry differ among the different studies. Given these differences and taking into account the uncertainty ranges, we can confidently consider that our scalings are compatible with those obtained in previous works.
The obtained scaling laws will be used to extrapolate our findings to magma oceans relevant Ra numbers and corresponding parameters in Section 4.

Characterization of the fluid ability to reach shallow depths
As mentioned above, the outgassing efficiency of the magma ocean first relies on the ability of the silicate magma to reach the exsolution depth, regardless of the ability of exsolved volatiles to reach the surface prior to possible reinjection and re-dissolution of oversaturated volatile species at higher lithostatic pressures via downwelling motions. For this reason, we characterize the ability of vigorously convecting fluid parcels to reach shallow depths as a function of time. This can be seen as the optimistic limiting case where gas bubbles, once formed at a given exsolution depth, instantaneously reach the surface. In this case, the degassing depth and the exsolution depths are equivalent.
where U is the time-dependent, dimensionless, spatially-averaged magnitude of the flow velocity. Therefore, the time required to achieve one transit, N transit " 1, corresponds to the time required for a convection-driven fluid parcel to travel a distance equivalent to that between one horizontal surface to the other. More To quantify the amount of fluid that reaches a given depth as a function of time, we use the minimum depth reached by the tracers averaged over the entire population of the particles used, xd min y: where N tracers is the total number of tracers. Figure 6a shows the time evolution of xd min y for the different Rayleigh number values considered in this study.
xd min y appears to strongly depend on the value of the Rayleigh number. The decrease of xd min y is faster when the Rayleigh number is larger. This is due to the fact that convective velocities, and thus convective efficiency, increase with the Rayleigh number so that the time needed to achieve the same degree of mixing and homogenization is much shorter at higher Rayleigh numbers.
However, using the number of transits (Equation (8) We observed two distinct phases for the evolution of xd min y: first, a strong and rapid decrease occurs, during the first two transit times (required for the system to "forget" the initial position of the tracers), followed by a more slowly decreasing trend of xd min y (Figure 6b). For example, for Ra " 10 12 , within will be used as the characteristic property of the distribution of d min in the next Section (in Equation (12)).
Previous studies reported on the fact that a third spatial degree of freedom may affect both the global (i.e., N u and Re numbers) and local properties of the convective flow, especially at low P r number values (e.g., Schmalzl et al., 2004;van der Poel et al., 2013). In particular, Schmalzl et al. (2002Schmalzl et al. ( , 2004 suggested that differences between 2D and 3D convection may appear as a result of (i) the constrained motion of the large-scale circulation in 2D flows compared to 3D geometry, and (ii) the increasing partitioning of energy in the toroidal component of the velocity in 3D flows at P r ă 1, otherwise absent in 2D geometry, thus strengthening the large-scale component of the flow . However, other parameters such as the boundary conditions, the aspect ratio, or the existence of multiple states and potential multi-stability observed around P r " 1 may also play a significant, yet not sufficiently constrained, role in these discrepancies (e.g., van der Poel et al., 2011van der Poel et al., , 2012van der Poel et al., , 2013. Therefore, to test the aforementioned influences on our results we performed a series of numerical experiments using a 3D geometry. These experiments considered a larger aspect ratio (of two in the horizontal directions) to ensure that the flow motions were not artificially constrained to a unit cell in our 2D experiments. Figure 8 (10) and (11)) are in good agreement with the 3D calculations. This demonstrates that the ability of the convecting fluid to transport material to shallow depths remains comparable in 3D. This confirms that the main quantity that affects the convecting outgassing efficiency is the magnitude of vertical velocities, regardless of the existence of toroidal motion.
The observed good agreement between 2D and 3D cases implies that we can confidently apply the scalings derived from our 2D experiments.

Description of the global distribution of the minimum depth reached by fluid parcels
In addition to the values of xd min y and σ dmin , a more general description of the distribution of the minimum depth reached by the tracers (Figures 5e to 5h) is required to accurately predict the magma ocean outgassing efficiency. As seen in Figures 5e to 5h, when the number of transit times increases, the distribution of d min becomes asymmetric as a result of convective motion. Using σ dmin to constrain its spread, such a distribution can be approximated by a truncated asymmetric Laplace distribution, whose probability density function writes: where the parameters µ " 0 (leading to a truncated asymmetric distribution) and b " b σ 2 dmin {2, with σ dmin given by Equation (11). In principle, Equation (12) needs to be integrated between d min "´8 and d min "`8 for its integral to equal unity, i.e., for encompassing the entire range of values that can be taken, which is nonphysical in our case since d min ranges between 0 and 1. Therefore, we normalize Equation (12) such that the probability density function equals unity when integrated for depths ranging between 0 and 1:

L "
Lapd min q ş dmin"1 dmin"0 Lapd min q dd min . 30 The above normalized probability density function, L, is shown by the blue curve in Figures 5g, 5h, and 9b.

Magma ocean dynamic outgassing model and simple predictions
Equation (13) can be used to estimate the amount of fluid that has reached a certain critical depth, d c , for a given number of transit times, by computing its integral between d min " 0 and d min " d c (yellow shaded area in Figures 9a and   9b). If we assume this critical depth to be the exsolution depth, i.e., d c " d exsol , we can compute the exsolved fraction of the fluid denoted X exsolved : which corresponds to the fraction of the melt that has reached the exsolution depth and thus has been able to exsolve its oversaturated volatiles. This predicted value for X exsolved can be compared to the value obtained in our numerical experiments, X exp exsolved , given by: where N tracers pd min ď d c q is the number of tracers that have reached the critical depth. This comparison allows one to estimate the uncertainty on the predicted magma ocean exsolved fraction computed, ∆X exsolved " X exsolved´X exp exsolved , illustrated by the deviation between the theoretical Laplace distribution (blue curve in Figures 5g, 5h, and 9b) and the observed, empirical distribution of d min (red histogram in Figures 5g, 5h, and 9b). Note that ∆X exsolved depends both on the critical depth, d c , and on the number of convective transits considered.
Out of 200 critical depth values for more than 50 number of transits tested, we found that the standard deviation of the error of the theoretical magma ocean exsolved fraction is ∆X exsolved "˘6%. Even though the error was computed for Ra " 10 12 only, choosing our maximum value for Ra is reasonable given the fact that Rayleigh number values for global magma oceans are high.  (11)). (a) The yellow area represents the magma ocean exsolved fraction corresponding to the tracers that have reached a given critical depth since the initial time t " 0, dc « 8.85ˆ10´3, corresponding to 26.55 km when considering a 3000 km deep magma ocean (horizontal dashed line). (b) The corresponding theoretical value is computed through the integration of the yellow shaded area (Equation (14)).  Table 1.  Figure 10). Details on the relationship between the number of transits and the duration of a given magma ocean transient stage, t MO , are given in the next Section. For the hypothetical 3000 kilometers deep magma ocean considered in our reference case (Table 1)

Application to Outgassing in Vigorously Convecting Silicate Planetary Magma Oceans
In the previous Section our estimates of the magma ocean outgassing efficiency rely on two sets of parameters arbitrarily chosen: (1) the parameters that govern the magma ocean dynamics considered to estimate the time spent by the magma ocean in a given regime, t MO (chosen here to correspond the hard turbulent convection regime, see Table 1), and (2)

Estimation of the exsolution depths
Realistic exsolution depths can be determined out of solubility laws, as a function of the abundance of volatile species present in the melt. These laws give the lithostatic saturation pressure of a given volatile species as a function of its dissolved amount within the melt, which is usually expressed in terms of volatile mass fraction relative to the mass of the melt, X MO vol . Thus, they control the partitioning of a volatile species between the silicate melt and gas phase according to a power-law relationship between saturation pressure and volatile abundance in the magma following (e.g., Carroll and Holloway, 1994;Holloway and Blank, 1994;Papale, 1997;Wallace and Anderson Jr., 1999): where X MO vol " M vol {M MO is the mass fraction of the volatile species expressed as the mass of the volatile species within the melt, M vol , relative to the mass of the melt, M MO , ζ vol is a constant, and β vol is the power law coefficient (see Table 2 for values considered here). For a given species, the saturation pressure, P sat vol , is the overburden pressure threshold below which (above in terms of depth) no more volatile can be dissolved within the melt, i.e., the maximum amount of dissolved volatiles within the liquid is reached, and any extra oversaturated volatile is then extracted out of the solution (i.e., "exsolved") to form gas bubbles. In coupled magma ocean-atmosphere studies, the saturation pressures computed out of solubility laws are directly considered to be equivalent to the atmospheric gas partial pressures, i.e., P sat vol " P atm vol , and thus sometimes referred to as the gas pressures within the melt, and are used to compute the corresponding atmospheric mass of the volatile species, M atm vol , following (e.g., Elkins-Tanton, 2008;Pierrehumbert, 2010;Hamano et al., 2013;Lebrun et al., 2013;Hamano et al., 2015;Schaefer et al., 2016;Salvador et al., 2017; where µ vol {μ is the ratio of the molar mass of the volatile species to the mean molar mass of the atmosphere, which equals one when considering single-component atmospheres and is thus omitted in such a case. As discussed earlier, this direct correspondence between the volatile species mass fraction in the magma ocean and its partial pressure in the atmosphere, and the resulting implicit partial pressures equilibrium between the gas phase in the melt and the atmosphere implies the instantaneous, and thus efficient, outgassing of any oversaturated volatile species at each computational time step (efficient outgassing hypothesis), thus neglecting the effect of magma ocean convective dynamics on the amount of melt actually reaching the exsolution depth and able to form gas bubbles. Again, here we revisit this assumption of efficient and instantaneous outgassing by realistically computing the amount of melt reaching the exsolution depth as a function of time. The depth at which volatile species are exsolved out of the melt and start forming gas bubbles is thus the exsolution depth, d exsol , which is calculated from the corresponding saturation pressure. Assuming hydrostatic equilibrium, the saturation pressure of a given volatile species, P sat vol , is: For the sake of simplicity, the above quantity is assumed to be constant The coexistence of several volatile species dissolved into the melt, such as H 2 O and CO 2 , influence their respective solubility, leading to more complex relationships between the volatile abundances and saturation pressures (e.g., Papale, 1999;Berlo et al., 2011;Massol et al., 2016). For the sake of simplicity, we will restrict our analysis to the case of a single-component atmosphere, made of pure water vapor H 2 O, which is one of the most important greenhouse gas and plays a fundamental role in planetary habitability (e.g., Wallace and Hobbs, 2006;Pierrehumbert, 2010;Westall and Brack, 2018). For the solubility of water within the magma ocean, we use the law classically referred to as Carroll and Holloway (1994) (whose parameters are given in Table 2), which is commonly used in coupled magma ocean-atmosphere studies (e.g., Lebrun et al., 2013;Salvador et al., 2017;Nikolaou et al., 2019;Bower et al., 2019). Implications of other water solubility laws are discussed in Appendix A. Although the uncertainty related to the choice of the solubility law is not negligible, it would not affect the conclusions. We assume a given water budget, expressed as a number of current Earth ocean mass that translates to a given mass fraction of dissolved water relative to the melt mass, X MO H2O , which, together with the solubility law (Equation (16)), constrains the value of the exsolution depth, d exsol (Equation (18)).

Magma ocean governing parameters and convective timescales
The number of transits achieved, N transits , relates to the time spent in a given magma ocean transient stage, t MO , following: where τ transit is the transit timescale corresponding to the time required for a parcel of the convecting fluid to travel from one horizontal surface to the other, which is written: where L is the thickness of the molten layer, and U Ω refers to the convective velocity accounting for the effect of planetary rotation. The latter may influence convection and the magnitude of the associated fluid motion (Solomatov, 2000;Hansen, 2015, 2019). Therefore, following Solomatov (2000, Equation (11)), this effective convective velocity is written: where Ω is the planet's angular velocity, C P is the specific heat at constant pressure, v conv " Re pρL{ηq´1 is the convective flow velocity computed using our Re´Ra scaling (Equation (7)), and F is the convective heat flux at the surface of the magma ocean: where N u is the Nusselt number computed using Equation (6)   * Using the mass-radius relationship given in Kopparapu et al. (2014) and based on the exoplanets.org known exoplanets database (Wright et al., 2011). velocities, we considered the uncertainty on two poorly constrained parameters that affect the magma ocean dynamics: its thermal expansion, α, and the temperature difference across the upper thermal boundary layer, ∆T . To account for the possible range of variation of these three parameters (U Ω , α, and ∆T ) and their influence on the magma ocean dynamics, we considered two end-member cases in addition to our reference case: (i) an "enhanced degassing limit", where convection is the most efficient due to high convective velocities (3 U Ω ), a large temperature difference across the thermal boundary layer (∆T " 1000 K; e.g., Lebrun et al., 2013), and a large thermal expansion (α " 10´5 K´1), and (ii) a "weak degassing limit" associated with less vigorous convection, and with smaller convective velocities (U Ω {3), smaller temperature difference (∆T " 20 K; e.g., Lebrun et al., 2013), and smaller (α " 10´6 K´1) thermal expansion (Table 4).
We compared the time required to fully exsolve all the oversaturated water from the magma ocean for the three different cases to the time spent by the magma ocean in a given transient stage, t MO . To do so, we evaluate the degassing efficiency of the global, fully molten magma ocean (corresponding to our transient stage) and consider the pair of Ra and P r numbers and associated parameters as constant, for a given water content. This corresponds to the earliest and most vigorously convecting magma ocean stage described in Figure 1b. Coupled magma ocean-atmosphere models (e.g., Solomatov, 2000;Lebrun et al., 2013;Nikolaou et al., 2019) show that depending on the ability of the growing atmosphere to efficiently loose heat to space, the duration of such a deep, vigorously convecting, fully molten stage can vary from t MO min « 10 years, when no atmosphere is present, to up to t MO max « 500 years when a thick and thermally opaque atmosphere overlies the magma ocean and buffers the heat loss to space (e.g., Lebrun et al., 2013). Considering these two end-member scenarios, the complete solidification of the terrestrial magma ocean, commonly defined as the time at which the rheological transition front reaches the surface (i.e., when the "mush" stage is reached; panels a and d in Figure 1) ranges from a thousand years to several millions years (e.g., Elkins-Tanton, 2008; Lebrun et al., 2013;Hamano et al., 2013Hamano et al., , 2015Salvador et al., 2017;Nikolaou et al., 2019;Bower et al., 2019;Krissansen-Totton et al., 2021a).

Results: outgassing model predictions
The outgassing efficiency of a magma ocean depends on the duration of the magma ocean stage in the considered transient stage, and on the ability of convective motions to bring molten material to shallow depths at which the oversaturated volatiles can exsolve within such a time interval. We considered these two aspects as described below. In addition, we make the assumption that water rapidly reaches the surface once exsolved, without being recycled/ingassed back deeper in the mantle by downwelling currents, and we neglect the influence of water concentration on magma ocean viscosity.  Tables 1 and 3, respectively, for the reference cases (plain blue and red curves, respectively), and in Table 4 for the enhanced and for the weak degassing limits (dashed blue and red curves, respectively). The magma ocean thickness for the super-Earth case is scaled based on Earth's mantle-to-planet radius ratio. For both planets, the water content is fixed to one Earth ocean mass, 1 M EO " 1.41 0 21 kg, initially dissolved and homogeneously distributed within the magma ocean. For the 3000 km deep Earth-sized magma ocean, this corresponds to a mass fraction of dissolved water relative to the magma ocean mass equal to X MO H2O " 100ˆ1 MEO MMO " 3.6ˆ10´2 wt%, which yields an exsolution depth, d exsol , of 5ˆ10´3 km (Equations (16) and (18)). For the 5 Earth mass rocky planet, X MO H2O " 7.73ˆ10´3 wt%, i.e., the same mass of water diluted within a larger magma ocean volume corresponds to a smaller exsolution depth d exsol " 3.14ˆ10´4 km. The plausible volatiles exsolved fraction range corresponding to the plausible variation range of parameters explored (U Ω , α, and ∆T ) lies within the blue and red areas of Figures 11a and 11b, respectively.
For the reference cases (Tables 1 and 3), 767 and 969 convective transits are required for the entire magma oceans volumes to reach the corresponding exsolution depths (Equation (14)), i.e., for the entire oversaturated water to be exsolved and outgassed, for the Earth-sized planet and for the super-Earth planet, respectively. This corresponds to 61 and 59 years (Equation (19)), respectively. Despite significantly different exsolution depths, the resulting time required to fully exsolve oversaturated H 2 O from the melt is similar. This is due to the fact that the larger magma ocean thickness and gravity on larger planets increase the Ra number, leading to more vigorous convective motions . The exsolution depth is computed considering the mass fraction of 1 Earth ocean of water dissolved within a magma ocean of thickness an Earth-like planet-to-mantle radius ratio. Reference cases (plain lines) parameters for the Earth-sized planet and super-Earth are given in Tables 1 and 3, respectively. Parameters considered for the enhanced and weak degassing limits (blue and red dashed lines) are listed in Table 4. Vertical black dashed lines indicate the minimum and maximum duration of the magma ocean in the fully molten transient stage from Lebrun et al. (2013) coupled magma ocean-atmosphere model. on super-Earths. This offsets the effect of shallower exsolution depths (requiring more convective transits to be achieved to outgas volatiles), and results in a more efficient degassing on super-Earths. Overall, magma ocean exsolution 800 starts slightly later on super-Earths than on Earth-sized planets, but full degassing is reached after similar amounts of time spent in the considered magma ocean transient stage, t MO , leading to a slightly faster magma ocean degassing on super-Earths.
For both types of planets (Figures 11a and 11b, respectively), the timing of degassing for the reference case lies between the shortest and longest magma ocean fully molten stage duration [t MO min; t MO max]. This implies that none of the potentially available oversaturated water would have been exsolved and outgassed in the fast cooling case (t MO min), while it would be fully exsolved and outgassed in the slow cooling scenario (t MO max). Furthermore, this indicates that for a conservative set of magma ocean governing parameters, complete or partial outgassing remain possible. However, if one considers a weak degassing limit, exsolution never starts before the end of the magma ocean transient stage considered, even for its longest duration, t MO max, regardless of the planet size (Figures 11a-b). In contrast, the enhanced degassing limit shows the opposite behavior since the complete outgassing always occurs before the end of the magma ocean transient stage considered, even for its shortest duration, t MO min, associated with the fastest cooling rate (Figures 11a-b). Yet, the cooling rate itself is controlled by the amount of greenhouse gases in the atmosphere and thus by the outgassing rate. Therefore, if the outgassing is not efficient and slow as in the weak degassing limit scenario or for small planets, the cooling rate will be high, allowing less time for the magma ocean to outgas, thus making the degassing efficiency weaker and so on. Conversely, for the fast, enhanced degassing limit scenario and for large planets, the amount of greenhouse gases outgassed will slow down the solidification, making the degassing more efficient by increasing the duration of the magma ocean transient stages and leaving more time for the melt to reach the exsolution depth.
Therefore, between the extreme enhanced and weak degassing limits, various combinations of parameters governing magma ocean dynamics would result in either a fast and complete outgassing, an incomplete degassing, or no degassing at all. These results emphasize the fact that magma ocean degassing strongly depends on its dynamics and associated parameters, and underline once again that complete outgassing during the most vigorously convecting, fully molten magma ocean stage (and by extension, during subsequent, less vigorous transient stages) is far from being certain even for magma oceans lasting for several thousands of years in the fully molten state (Figures 11a-b).
Since we assumed that the same amount of water is dissolved within a smaller magma ocean volume, the Earth-sized magma ocean is more concentrated in water compared to that of the super-Earth case, as it would be in a super-Earth magma ocean at an advanced crystallization stage (because water partitions preferentially into the remaining melt, whose volatile concentration increases with solidification). This, in addition to its smaller magma ocean thickness, Rayleigh number, and associated convective vigor, makes the Earth-sized case the analog of a more advanced crystallization stage than the super-Earth case.
While we focused on the most vigorous fully molten stage, similar outgassing timings between these two different planetary cases, that would be representative of two slightly different convection efficiency/crystallization stages, suggest that our conclusions would apply for later transient crystallization stages and associated less vigorous convective dynamics (characterized by other pairs of Ra´P r values experienced during subsequent cooling stages). This suggests that degassing efficiency would decrease with time, along the solidification of the magma ocean due to the decrease of its convective vigor buffering the increase of the exsolution depth associated with water enrichment of the melt.  (Figures 12b, 12d, and B.4).

Influence of the initial volatile content
For completeness, the influence of the exsolution depth on the exsolution time is discussed in Appendix B. This behavior is even more pronounced on highly water-concentrated Earth-sized magma oceans, where a tiny fraction of their volume almost instantaneously reaches the large exsolution depth associated to 10 M EO (lower left corner in Figure 12b). Conversely, smaller water contents decrease the exsolution depth ( Figure A.1) and slow down degassing, regardless of the planet size (Figures 12a and 12c). For magma oceans with such a low water content, the time required to exsolve 100% of oversaturated water is generally longer than the time spent in the fully molten stage (Figures 12a and   12c), implying that water degassing is more likely to be incomplete or may not even have started before the end of the fully molten stage.
The fact that early degassing of water may be incomplete or even nonexistent in some cases has strong implications regarding the entire cooling history of the magma ocean, the resulting planetary surface conditions, and even the chemical differentiation and subsequent long-term evolution of the planet (Maurice et al., 2017;Ballmer et al., 2017). Indeed, if degassing is not efficient, the weak thermal blanketing effect of the resulting tenuous steam atmosphere will allow for efficient heat escape to space and rapid magma ocean cooling, making degassing less efficient. Such a feedback would favor water storage in the mantle, resulting in early dry surface conditions, after rapid magma ocean solidification, with a water-rich solid-state mantle. This may have important implications regarding its properties and long-term evolution, including dynamics, convective regime, melting and degassing (e.g., Lange, 1994;Regenauer-Lieb et al., 2001;Demouchy and Bolfan-Casanova, 2016;Ni et al., 2016;Ohtani, 2020;Korenaga, 2021, 2022). Conversely, a large initial concentration of water would induce an efficient and rapid exsolution of the magma ocean oversaturated water, which would, in turn, slow down the magma ocean cooling because of the efficient greenhouse effect of a thick steam atmosphere overlying the molten surface, and therefore leave time for the magma ocean to fully outgas its oversaturated volatile content during the slow solidification process. The slow solidification of a magma ocean due to the presence of its thick steam atmosphere and the resulting loss of water through atmospheric escape has been proposed, as a consequence of the orbital distance, to explain the divergence between the Earth and Venus (Hamano et al., 2013). In addition, the atmospheric water content affects the spectral properties of exoplanets (Robinson, 2020;Smith et al., 2020) and high water contents atmospheres even likely experience a runaway greenhouse radius inflation that changes their mass-radius relationship in a possibly future detectable manner (Turbet et al., 2019(Turbet et al., , 2020.
These aspects emphasize the crucial role that magma ocean outgassing efficiency may play in planetary evolution, habitability, and exoplanets detectability and characterization.

Model limitations and future directions
By highlighting the influence of the dynamics of the fully molten magma ocean stage on its early degassing, our results motivate further investigations to characterize magma ocean degassing efficiency during its entire lifetime.
This implies extending the approach presented here to values of Prandtl numbers larger than one, to better cover the different transient stages experienced throughout magma ocean solidification. One should note that due to the possible differences in convective dynamics in 2D and 3D simulations at P r ă 1 Breuer et al., 2004), the validity of our results may be limited to Prandlt numbers of unity and larger. Nevertheless, such values are relevant for terrestrial magma oceans (e.g., Höink et al., 2006;Solomatov, 2015;Hansen, 2015, 2019, andreferences therein). While such a task is out of the scope of the present paper, our numerical experiments and outgassing predictions show that for a significant range of magma ocean governing parameters, the degassing is not as efficient as previously thought (localized and minimal to no degassing scenarios proposed in Ikoma et al., 2018). Therefore, the assumption of instantaneous outgassing of oversaturated volatile species may not always be valid and should only be considered with care.
Since the cooling rate of the magma ocean is controlled by its degassing via the amount of outgassed atmospheric greenhouse gases, our results suggest that magma ocean simulations shall consider the influence of magma ocean dynamics to self-consistently infer the magma ocean lifetime and the associated surface conditions.
The expressions we derived to predict the spectrum of the minimum depth/pressure reached by the convective fluid with time can be used in coupled atmospheremantle models. This may reveal additional complexities during and after the magma ocean evolution, due to the influence of volatile content on several key physical quantities that govern magma ocean and the subsequent solidified man-50 tle convective dynamics, such as viscosity or density (e.g., Lange, 1994;Litasov and Ohtani, 2007;Karki and Stixrude, 2010;Nikolaou et al., 2019). Ultimately, this coupled approach will provide a considerably more accurate picture of the outgassing history of terrestrial magma oceans together with the formation of planetary atmospheres.

Conclusions
We investigated the influence of magma ocean convective dynamics on their degassing efficiency using numerical fluid dynamics modeling of thermal con- We used our experimentally-derived expressions to estimate the magma ocean volume fraction that reaches any given exsolution depth as a function of time. We found that a substantial number of convective transits need to be achieved for the entire vigorously convecting magma ocean volume to reach the exsolution depth, and thus, for all oversaturated volatile species to be exsolved out of the melt, form gas bubbles and degas at the surface. This is particularly important for volatile species that are highly soluble in the melt, such as water. Depending on the value of the governing parameters considered, the time required to fully exsolve oversaturated volatile species can be comparable or even longer than the time spent in a fully molten transient stage. Therefore, in some cases, a fully molten mantle can start to crystallize before its entire volume has reached the exsolution depth. For example, despite vigorous convection, we found that complete outgassing of water can require more than 10 4 years in a fully molten global terrestrial magma ocean, a duration that may be

Appendix A. Effect of solubility laws
Solubility laws provide the saturation pressure, P sat vol , of any given volatile species as a function of its concentration in the melt. Thus, they give the overburden lithostatic pressure where a given volatile species has reached its maximum concentration in the melt and below which no more can be dissolved into the melt and start forming gas bubbles. Assuming hydrostatic equilibrium, the corresponding exsolution depth, d exsol , is obtained from Equation (18). Figure A.1 shows different solubility laws for water in silicate melts (see Table A.1 for the corresponding solubility laws coefficients). In this study, we considered the Carroll and Holloway (1994) solubility law for water (Equation (16), Table 2, black line in Figure A.1), which has been used in several magma ocean modeling studies (e.g., Lebrun et al., 2013;Salvador et al., 2017;Nikolaou et al., 2019;Bower et al., 2019) and gives intermediate saturation pressures compared to others. Solubility laws obtained for basaltic systems are generally more appropriate for magma oceans as basalts, i.e., relatively low silica content melts, have high temperatures, low viscosities and high gas diffusivities (e.g., Sparks et al., 1994;Lesher and Spera, 2015). However, they tend to have more modest volatile contents than more viscous, silicic melts such as rhyolitic or dacitic melts. However, note that the viscosity of mafic and ultramafic melts (such as basalt and komatiite, respectively), is comparatively much less affected by the dissolved water content than that of polymerized, low temperatures, rhyolitic or dacitic melts (e.g., Lesher and Spera, 2015, their Holloway and Blank (1994); Papale (1997), respectively), and derived from rhyolitic melts (green line in Figure A.1, from Table 1 in Lichtenberg et al. (2021) and references therein). For each solubility law considered, symbols mark the saturation pressure for three different water contents spanning three orders of magnitude, from 0.1 (circles) to 10 Earth ocean mass (squares), with one Earth ocean mass (1 M EO ) equal to 1.4ˆ10 21 kg. Note that for the same mass of water, the mass fraction of dissolved water, X MO H2O , which is the absolute mass of water relative to the mass of the melt, is larger for Earth-sized planets than for super-Earths as the same amount of water is dissolved and diluted within a larger, more massive, molten mantle in the latter. Considering such similar water contents for planets of different sizes implicitly makes the assumption that water endowment does not scale with planetary mass during the accretion sequence. Consequently this does not account for the possible effect of planet size on volatile abundance and assumes that water delivery is relatively less efficient on larger planets. This accretion scenario could be in agreement with the stochastic nature of volatile delivery, which may then be weakly dependent on planetary mass (e.g., Raymond et al., 2006Raymond et al., , 2007. However, the uncertainties associated to volatile delivery sources and timing do not allow for a conclusive statement yet (e.g., Marty, 2012;Morbidelli et al., 2012;O'Brien et al., 2018;Venturini et al., 2020). As seen in Figure A.1, for a given water content, the solubility law chosen strongly influences the saturation pressure and thus the exsolution depth at which gas bubbles start to form. Conversely, the amount of water that can be dissolved into and outgassed out of the melt strongly varies depending on the solubility law chosen. These differences become more and more important when decreasing the mass fraction of dissolved water and are thus (i) larger for low water contents, and (ii) always larger for super-Earths than for Earth-sized planets for a given dissolved water mass. As a result, on super-Earths, the saturation pressure of 0.1 Earth ocean mass can vary over more than four orders of magnitude solely depending on the arbitrary solubility law choice (red circles). As shown in this study, the exsolution depth, and thus the choice of the solubility law, play a fundamental role regarding magma ocean degassing efficiency as the latter is determined by comparing the time required for the entire molten volume to reach this very depth and the time spent in a given transient stage (of given convective dynamics). The solubility law considered is thus of fundamental importance when constraining the timing of atmospheric formation and magma ocean solidification, as well as the surface  Earth-sized planets exsolution depth, d exsol [km] 10 -4 10 -3 10 -2 10 -1 10 0 10 1  (1994) Papale (1997) Lichtenberg et al. (2021) Carroll  x 2.4 . At a given cooling transient stage, an ascending fluid undergoes decompression, i.e., a decrease of ambient lithostatic pressure, and goes from the right hand side to the left hand side of the diagram. As solidification proceeds, the silicate melt volume decreases and water partitions readily into the melt such that evolving magma ocean liquids are enriched in water, and the mass fraction of dissolved water X MO H 2 O increases along the saturation curve, from the bottom to the top of the diagram. In this study, we used the Carroll and Holloway (1994) solubility law (black, parameters given in Table 2). Parameters of the other solubility laws are given in   Tables 2 and A.1 for solubility laws parameters). Reference cases (plain lines) parameters for the Earth-sized planet and super-Earth are listed in Tables 1 and 3, respectively. Parameters considered for the enhanced and weak degassing limits (blue and red dashed lines) are listed in Table 4. Vertical black dashed lines indicate the minimum and maximum duration of the magma ocean in the fully molten transient stage from Lebrun et al. (2013) coupled magma ocean-atmosphere model.  Tables 2 and A.1 for solubility laws parameters). Reference cases (plain lines) parameters for the Earth-sized planet and super-Earth are given in Tables 1 and 3, respectively. Parameters considered for the enhanced and weak degassing limits (blue and red dashed lines) are given in Table 4. Vertical black dashed lines indicate the minimum and maximum duration of the magma ocean in the totally molten transient stage from Lebrun et al. (2013) coupled magma ocean-atmosphere model.

Appendix B. Influence of the exsolution depth on the magma ocean exsolution time
One can estimate the efficiency of magma ocean outgassing using the exsolution time, t exsol , defined as the time required for 99.99% of the magma ocean to reach a given exsolution depth at which its oversaturated water is fully exsolved and outgassed. This quantity is displayed as a function of the exsolution depth in Figure B.4 for the Earth-like and super-Earth reference cases (blue and red plain lines, whose parameters are given in Tables 1 and 3, respectively), with the corresponding plausible ranges included between the weak and efficient degassing limits (blue and red areas in Figures B.4a and B.4b, respectively, whose parameters are listed in Table 4). As discussed above, an increase of the exsolution depth can result either from an increase of the initial concentration of dissolved water within a magma ocean of given size, or from the water incompatibility and resulting enrichment of the residual melt during magma ocean crystallization. Regardless of planetary size, the timing for complete magma ocean exsolution decreases when increasing the exsolution depth ( Figures B.4a and B.4b), since deep exsolution pressures are more rapidly reached, implying that magma ocean degassing is more efficient for larger water contents, as discussed above. Yet, one cannot state that magma ocean degassing becomes more and more efficient as crystallization proceeds only by looking at the exsolution depths and based on the fact that melt becomes more water-concentrated. Indeed, the decrease of the convection vigor with solidification would buffer the effect of larger exsolution depths. This underlines the importance of considering these two competing effects to accurately characterize magma ocean degassing efficiency. Finally, depending on the duration of the fully molten stage, almost all scenarios in the parameters space explored predict that the magma ocean complete exsolution time is larger than t MO min " 10 years, thus generally preventing complete magma ocean degassing for this transient stage. Even when considering slow magma ocean cooling associated with the longest duration of the fully molten stage (t MO max " 500 years), about half of the scenarios, ex-  Tables 1 and 3, respectively. These parameters have been used to compute the reference cases exsolution times (plain lines), and parameters for the weak (upper dashed lines), and enhanced (lower dashed lines) degassing limits are listed in Table 4.
cluding the reference case, lead to incomplete or no degassing (t exsol ą t MO ).
In such cases, the amount of water that remained dissolved within the melt (because part or none of the latter have reached the exsolution depth) must be accounted for when considering the entire magma ocean cooling sequence.
Appendix C. Benchmark results of the StreamV code in 3D geometry and in the finite Prandtl number context We considered a convection benchmark proposed in Fusegi et al. (1991), which is an extension of the 2D test proposed in Davis and Jones (1983) to 3D geometry. It consists in a cubic domain in which two opposite vertical sidewalls (chosen here to be located perpendicular to the x-coordinate, while the z-axis is vertical and oriented opposite to the gravity vector) have homogeneous distinct temperatures. The four other walls are insulating, and no-slip conditions are imposed on all boundaries. The material properties are homogeneous, and in this case the system is governed by the Prandtl number and the Rayleigh num-61 ber already defined (with Ra being defined here in terms of the temperature difference between the two isothermal vertical sidewalls). The Prandtl number is set to 0.71, but the Rayleigh number is varied between 10 3 and 10 6 . Convective currents progressively develop and the system is evolved until reaching steady state at which several quantities are compared: (i) the average Nusselt number taken at the isothermal boundaries. In the conservative finite-volume formulation, the latter is the same on both isothermal boundaries at the steady sate; (ii) the maximum values for the x and z components of the velocity field: v MP x max " maxrv x py " 0.5qs, v MP z max " maxrv z py " 0.5qs in the symmetry plane (y=0.5); (iii) the maximum values for each component of the velocity field in the entire domain: v xmax , v ymax , v zmax . We performed a comparison with the study by Labrosse et al. (2000) that documents all the diagnostics mentioned above. Cases for which Ra is lower than 10 5 are computed on a grid using 64ˆ64ˆ64 cubic cells, while the other cases use 128ˆ128ˆ128 cubic cells. The results displayed in Table C.2 show that the agreement is good with a relative error smaller than one percent for each quantities for each case considered. • Volatile outgassing in a vigorously convecting magma ocean may be considerably weaker than previously thought