RETRACTED: Redox Evolution of the Crystallizing Terrestrial Magma Ocean and Its Influence on the Outgassed Atmosphere

Magma oceans (MOs) are episodes of large-scale melting of the mantle of terrestrial planets. The energy delivered by the Moon-forming impact induced a deep MO on the young Earth, corresponding to the last episode of core-mantle equilibration. The crystallization of this MO led to the outgassing of volatiles initially present in the Earth’s mantle, resulting in the formation of a secondary atmosphere. During outgassing, the MO acts as a chemical buffer for the atmosphere via the oxygen fugacity, set by the equilibrium between ferrous- and ferric-iron oxides in the silicate melts. By tracking the evolution of the oxygen fugacity during MO solidification, we model the evolving composition of a C-O-H atmosphere. We use the atmospheric composition to calculate its thermal structure and radiative flux. This allows us to calculate the lifetime of the terrestrial MO. We find that, upon crystallizing, the MO evolves from a mildly reducing to a highly oxidized redox state, thereby transiting from a CO- and H2-dominated atmosphere to a CO2- and H2O-dominated one. We find the overall duration of the MO crystallization to depend mostly on the bulk H content of the mantle, and to remain below 1.5 millions yr for up to nine Earth’s water oceans’ worth of H. Our model also suggests that reduced atmospheres emit lower infrared radiation than oxidized ones, despite of the lower greenhouse effect of reduced species, resulting in a longer MO lifetime in the former case. Although developed for a deep MO on Earth, the framework applies to all terrestrial planet and exoplanet MOs, depending on their volatile budgets.


Introduction
During their formation, terrestrial planets are thought to go through phases of large-scale melting of their silicate mantles referred to as "magma oceans" (MOs; Elkins-Tanton 2012).Such MOs would represent the most important episodes of rapid chemical equilibration between the mantle and the atmosphere in a planet's history (Abe & Matsui 1985).Conversely, the solid mantle is chemically isolated, and chemical exchanges with the atmosphere proceed through extremely sluggish processes (for instance, volcanism or subduction of volatile-rich oceanic crust).Therefore, outgassing taking place during MO solidification largely sets the initial stage for the planetary atmospheric evolution.Characterizing it properly is thus key to understanding the differences between terrestrial planets in the solar system.Understanding MO-equilibrated atmosphere is also crucial to interpreting upcoming observations of exoplanetary atmospheres.
Terrestrial MO atmospheres, over the temperature range considered in this study (1500-3000 K), are thought to be largely dominated by the most abundant volatile elements (e.g., H and C).In addition to being of prime interest as life-essential elements, volatile elements also form gaseous species with strong absorption features in the infrared, affecting the radiative state of the atmosphere.They are, therefore, the focus of many MO-atmosphere fractionation studies (Kuramoto & Matsui 1996;Dasgupta 2013;Hirschmann 2016;Dasgupta & Grewal 2019).
Early work by Abe & Matsui (1985, 1988) laid the ground work for research on the link between the terrestrial MO and an outgassed H 2 O-CO 2 atmosphere, delivered by volatiles from accreting planetesimals.They pointed out that the thermal blanketing effect of a thick steam and CO 2 atmosphere prevents efficient cooling of the MO.Indeed, the presence of greenhouse gases in the atmosphere reduces the thermal radiative heat flux from the MO surface, while leaving the high wavelength instellation largely unaltered, thereby reducing the net cooling of the planet and slowing down its solidification.Later studies by Elkins-Tanton (2008) and Lebrun et al. (2013), considering the same gaseous species, emphasized the difference in the outgassing patterns of the two gases: CO 2 , having a low solubility in silicate melts, is readily outgassed from the MO since the beginning, while water, being much more soluble, remains largely dissolved until the late stages of crystallization.Addressing the thermal evolution of the MO, these studies showed that the initial volatile budget, especially that of water, has a first-order influence on the MO lifetime, whereby more volatiles induce a longer crystallization.The MO lifetime can even become infinite when the incoming stellar flux exceeds a threshold, associated with the buffering effect of water condensation on the atmospheric infrared emission (Hamano et al. 2013;Salvador et al. 2017), above which the MO no longer cools.This scenario might have been applicable to Venus, in which case only the dessication of its steam atmosphere by H 2 escape could have ended the MO phase.More recently, Bower et al. (2019) pointed out a mechanism by which the late outgassing of water decreases the partial pressure of CO 2 , casting a more complex picture of the outgassing pattern than observed thus far.In two related studies, Nikolaou et al. (2019) and Katyal et al. (2019) also investigated the sensitivity of the terrestrial MO duration to a series of parameters including the initial volatile budget, the convective regime in the MO, and the parameterization of the radiative flux through the atmosphere, confirming the prime influence of H budget.
All of these studies focused on oxidized atmospheres (composed of H 2 O and CO 2 ).Lately, a significant effort has been undertaken to broaden the scope to more general chemistry.Estimating the MO's oxygen fugacity ( f O 2 ) based on novel experiments on molten peridotite, Sossi et al. (2020) suggested that the early terrestrial atmosphere inherited from MO outgassing might have been akin to the atmosphere of present-day Venus, i.e., dominated by CO 2 and N 2 .Bringing S into the picture while treating f O 2 as a free parameter, Gaillard et al. (2022b) showed that reduced atmospheres are dominated by H-C species while oxidized ones are dominated by C-N-S species.They showed in particular that the present-day surficial reservoirs of C and N on the Earth can be matched in the atmosphere at equilibrium with a deep MO, inferring that later processes might have played only minor roles in the surface inventory of these elements.Other studies focused on quantifying the radiative state of the outgassed atmosphere, beyond CO 2 and H 2 O.For example, Lichtenberg et al. (2021) analyzed the effect of different gaseous species individually (in single-gas model atmospheres).They distinguished between three groups of volatiles: O 2 , N 2 , and CO have both a low solubility and poor greenhouse effect, yielding immediately outgassed atmospheres and short-lived MOs.CO 2 , CH 4 , and H 2 O have a stronger greenhouse effect and various solubilities (generally higher), inducing longer MO and evolving atmospheric pressure.Finally H 2 , was found to have the strongest effect on the lifetime of the MO, due to its strong collisioninduced absorption at high pressure.Focusing on the H-and C-bearing systems (i.e., H 2 O, H 2 , CO 2 , and CO), Katyal et al. (2020) calculated the redox speciation in the atmosphere overlying different stages of the MO evolution computed by Nikolaou et al. (2019).By varying the f O 2 of the MO at equilibrium with the atmosphere as a free parameter, they found that the transition from reduced to oxidized conditions translates to a decrease of the radiative cooling flux.Taking it one step further, Bower et al. (2022) computed the chemical evolution of an atmosphere outgassed throughout the crystallization of an MO with an oxygen fugacity at 0.5 log-unit below the iron-wüstite buffer (following Sossi et al. 2020).By investigating the bulk C/H ratio as well as the volatile budget, they found a general evolution toward a more oxidized atmosphere as crystallization proceeds.Furthermore, the combined outgassing of H-and C-bearing species generally results in a transition from a CO-dominated early atmosphere to an H 2 O-dominated late one, although a major part of the H can remain trapped in interstitial melts in the mantle.
These previous studies, therefore, established that MO oxidation state is a crucial parameter governing the chemistry of atmospheres.As long as the MO remains in equilibrium with the atmosphere and dominates it in terms of mass, one can assume that the MO chemistry buffers the oxidation state of the atmosphere.The redox chemistry in the MO is initially set by the equilibrium between metallic iron and iron oxide during the formation of the core.For an Earth-like composition, this event buffered the f O 2 to two log-units below the iron-wüstite buffer (setting ΔIW = −2, where ΔIW is the deviation between the actual f O 2 of the MO and that corresponding to the iron-wüstite buffer).Subsequently, after metallic iron segregates to the core, the redox state of the MO is governed by the equilibrium between ferrous-(Fe 2+ ) and ferric-(Fe 3+ ) iron oxides (i.e., FeO and FeO 1.5 ).Hirschmann (2012) proposed that the pressure-dependence of this equilibrium should induce a negative gradient in the f O 2 versus depth profile in an MO where the ratio between the FeO and FeO 1.5 concentrations is homogenized by convective mixing.In this case, depending on the effective pressure of core-mantle equilibrium, the initial f O 2 at the surface (where the MO equilibrates with the atmosphere) may deviate from ΔIW = −2.Recent studies using experimental (Zhang et al. 2017;Armstrong et al. 2019) or first principle molecular dynamics (Deng et al. 2020) methods have confirmed it.For a constant ferric-to-total iron oxide ratio (Fe 3+ /ΣFe) throughout the MO, the oxygen fugacity increases (relative to the iron-wüstite buffer) with decreasing pressure (except in the top ∼10 GPa).On Earth, the last major event of core-mantle equilibration is thought to be associated with the Moon-forming impact, which produced a deep MO.It is thus possible that the surface f O 2 was higher than IW-2 (Wood et al. 2006).
Previous studies on MO-atmosphere interactions have considered these uncertainties on the redox state by varying the f O 2 as a free parameter (Katyal et al. 2020;Sossi et al. 2020;Gaillard et al. 2022b;Bower et al. 2022).However, the quantity Fe 3+ /ΣFe, which governs the MO f O 2 , is expected to evolve as the MO crystallizes, because ferric iron is more incompatible than ferrous iron in silicate minerals (Boujibar et al. 2016;Davis & Cottrell 2021;Rudra & Hirschmann 2022).Therefore, Fe 3+ /ΣFe is expected to increase, making the MO more oxidized as it crystallizes.In other words, the redox evolution is expected to be coupled to the thermal evolution of the MO.In the present study, we self-consistently compute the evolution of the oxygen fugacity in the MO by tracking Fe 3+ /ΣFe.To do this, we use the parameterization linking f O 2 and Fe 3+ /ΣFe obtained by Hirschmann (2022), fitted on a P-Tcomposition space encompassing high pressures (from Deng et al. 2020), high (MO-like) temperatures, and peridotitic composition relevant for the MO (from Sossi et al. 2020).Based on the surface f O 2 of the MO, we compute the redox speciation in an H-C-O atmosphere.From the composition of the outgassed atmosphere, we calculate the cooling flux of the MO using a convective-radiative model allowing us to capture influential effects on the thermal evolution of the coupled atmosphere and MO system, which were not accessible to simpler atmospheric models.We vary the crystallization regime (i.e., the efficiency of crystals settling versus entrainment), the bulk volatile budget, and the bulk C/H ratio, to investigate their influences on the composition of the outgassed atmosphere, the thermal evolution and lifetime of the MO crystallization, and the retention of volatile in the mantle.We apply this approach to the case of the post-Moon-forming impact terrestrial MO.

Model
This study has three main model components.First, there is the crystallization model, which, depending on the dynamical regime of settling of the crystallizing minerals, allows us to relate the potential temperature of the MO to its depth, and to calculate elemental fractionation between solids and liquids.
Second, there is the redox chemistry model, which, based on the evolution of the Fe 3+ /ΣFe ratio in the MO, allows us to calculate the f O 2 and, in turn, the chemical speciation of the atmosphere in equilibrium with the MO.Considering mass conservation of volatile elements as well as solubility equilibria, we track the outgassing of the MO and the formation of the atmosphere.The outputs of the chemical model are the partial pressures of each species present in the atmosphere and their concentration in the MO (as well as f O 2 and Fe 3+ /ΣFe).Third, there is the thermal evolution model, which allows us to calculate the cooling flux of the MO radiating through its outgassed atmosphere, and thus constrain its lifetime.
In the following, Section 2.1 presents the technical definition of the MO depending on the dynamical regime of crystals settling.Section 2.2 presents the redox chemical model for both iron oxides and volatile species.Finally, Section 2.3 introduces the thermal evolution model.

Magma Ocean Crystallization
The crystallization of the terrestrial MO is generally thought to proceed from the bottom up to the surface as a consequence of the stronger pressure gradient of the melting curves compared to the temperature profile in the MO (Solomatov 2015), and/or of the higher density of solid crystals.Limits to the validity of both of these assumptions (Stixrude et al. 2009;Caracas et al. 2019) have been pointed out, which could lead to a middle-out crystallization of the MO, and possibly to the formation of a basal MO (Labrosse et al. 2007).In this work, we restrict the study to a terrestrial MO 55 GPa-equivalent in depth (Deng et al. 2020), as it represents the effective pressure of core-mantle equilibration.It is likely too shallow for a middle-out crystallization pathway to apply.
Under the assumption that crystals efficiently settle at the bottom of the MO and form a compact cumulates pile, the distinction between the MO and its solid cumulates is straightforward.However, when crystals remain in suspension and a partially molten domain exists, this distinction becomes blurred.Here, we envisage both cases, and refer to them as "fractional" and "equilibrium" crystallization, respectively.In the latter case, the "MO" refers to the layer in a liquid-like dynamical regime (undergoing turbulent convection), while "solid cumulates" refers to the layer in a solid-like dynamical regime (undergoing creep deformation).The bottom of the MO is then defined as the depth at which this transition between two dynamical regimes occurs, as described in Section 2.1.1.
The fate of crystals is thus essential to the distinction between fractional and equilibrium crystallization.The settling or entrainment of the crystallizing minerals in an MO is the result of a competition between convective currents and buoyancy forces (Solomatov & Stevenson 1993a;Patocka et al. 2020Patocka et al. , 2022)).It is poorly constrained, mainly due to the difficulties in reproducing MO-like dynamical regimes either experimentally or numerically.
Notations in use in Section 2.1 are listed in Table 1 2.1.1.Equilibrium Crystallization In this scenario (represented in Figure 1(a)), we assume that convective currents in the MO are stronger than the buoyancy force acting on the crystallizing solids.Thus, crystals remain entrained, and the solidifying medium has a melt fraction parameterized as a function of pressure, temperature, and the melting curves: where P is the ambient (lithostatic) pressure, T is the corresponding temperature on the MO adiabat (see Section 2.1.3),and T sol and T liq are the corresponding solidus and liquidus temperatures, respectively (Fiquet et al. 2010).
The medium behaves dynamically as a liquid until the ambient melt fraction reaches a threshold value, referred to as the rheologically critical melt fraction (RCMF; Solomatov & Stevenson 1993b).When the melt fraction crosses the RCMF, solid crystals start to form an interconnected lattice, which supports the deformation, resulting in a drastic increase of the effective viscosity (Costa et al. 2009).The value of the RCMF varies depending on multiple parameters like the shape and size of the crystals.In line with previous studies, we use a value of 0.4 (Solomatov 2015;Monteux et al. 2016;Nikolaou et al. 2019) for the main part of this study, but we also vary it in Section 4.1.The MO extends from the surface down to a melt fraction of 0.4, and always contains a significant mass of freefloating crystals.Similarly, the solidifying cumulates retain a significant quantity of trapped melt, which is assumed to freeze in place as the solid mantle further cools down, eventually contributing to the mantle's volatile budget as discussed in Section 4.1.The crystals in suspension in the MO remain in chemical equilibrium with the melt, while the melt trapped in the solid cumulates is chemically isolated from the MO, and its composition is no longer affected by MO fractionation.The initial temperature profile of the equilibrium crystallization scenario is represented in Figure A1(b).

Fractional Crystallization
In this scenario (represented in Figure 1(b)), we assume that the buoyancy force acting on the crystallizing phases is stronger than the convective currents in the MO.Thus, the crystals settle at the bottom of the MO and form a compact cumulates pile from which melt is entirely extracted.The mass of solid phases formed over one crystallization step is used to compute an equivalent-volume layer thickness at the top of the Note.Notice that in Tables 1, 2, and 3, "variable" means spatially changing, while "evolving" means temporally changing (when both apply, only the first one is indicated).

R e t r a c t e d
current cumulate pile and update the depth of the MO accordingly.
We assume that, at the onset of the MO's crystallization, its bottom is at the liquidus temperature at a pressure of 55 GPa.Thus the fractional crystallization cases exhibit a higher initial potential temperature (2939 K) than the equilibrium crystallization ones (2571 K), where the bottom of the MO is at the RCMF at 55 GPa.The initial temperature profile of the fractional crystallization scenario is represented in Figure A1(b).The minerals are in chemical equilibrium with the MO only during the crystallization step they are formed in, and then become chemically isolated from the next step on.

MO Adiabat
The temperature profile in the MO is used to calculate T-dependent quantities, in particular the melt fraction (Equation (1)), but also the equations of state (see below) or the MO viscosity (Equation (B1)).The temperature follows an adiabat (Solomatov 2015):

Chemical Model
The chemical model relies on three considerations: (1) mass conservation of Fe 2+ , Fe 3+ , H, and C, (2) redox equilibria between FeO (i.e., Fe 2+ ) and FeO 1.5 (i.e., Fe 3+ ), H 2 and H 2 O, and CO and CO 2 , respectively, and (3) dissolution equilibria: crystal-melt partitioning of Fe 2+ and Fe 3+ , and gas-melt dissolution of H 2 and H 2 O, and CO and CO 2 .Notice that we consider both H and C perfectly incompatible.Redox equilibrium of the FeO-FeO 1.5 system are calculated using the parameterization from Hirschmann (2022; Equation (A1)), and redox equilibria between the H-bearing and C-bearing species are calculated using temperature-dependent reaction constants from JANAF tables (Chase 1998;Appendix A.4) as well as the calculated f O 2 for.Crystal-melt partitioning of FeO and FeO 1.5 is calculated using a set of partition coefficients (Appendix A.3), and gas-melt equilibrium is calculated using solubility laws (Equation (11), Table 5).
The pathway of the chemical model is as follows: from Fe 2+ and Fe 3+ mass conservation, we track the relative abundances (assumed homogeneous in the MO) of both iron oxides, expressed as the molar ferric-to-total iron ratio, Fe 3+ /ΣFe.From Fe 3+ /ΣFe, we calculate the f O 2 using the redox equilibrium parameterization.We then use f O 2 to calculate the relative abundances of the gaseous species in the atmosphere.Finally, the system is closed by considering dissolution equilibria of volatile species and mass conservation of C and H.
Notations in use in Section 2.2 are listed in Table 2 2.2.1.Mass Conservation Mass conservation is enforced on a system consisting of the crystallizing MO (as defined in either Section 2.1.1 or 2.1.2) and the atmosphere.Notice that, for the sake of generality, here the mass conservation model is described for a hypothetical element that can enter solid crystals, melt, and atmosphere, while in our case Fe 2+ and Fe 2+ are distributed between crystals and melt while C and H are distributed between melt and atmosphere only.At each crystallization step, the inventory of element (or ion for Fe) e is distributed between three reservoirs: entrained crystals in the MO, melt in the MO, and R e t r a c t e d atmosphere: where M e syst is the total mass of e in the system (MO + atmosphere), M crys , M liq , and M atm are the masses of crystals in the MO, liquid in the MO, and gas in the atmosphere, respectively, and X e crys , X e liq , and X e atm are the mass fractions of e in the corresponding reservoirs.Notice that in the fractional crystallization scenario, M liq corresponds to the mass of the MO and M crys = 0.
The MO, considered as an open system, loses a mass dM syst to the cumulates pile during one crystallization step (see Appendix A.6).We neglect atmospheric loss, which would be another mass sink for the system in the case of volatiles; however, we discuss its potential effects in Section 4.3.dM syst is composed of a fraction f RCMF of liquid and a fraction (1 − f RCMF ) of crystals (notice that this also applies to the fractional crystallization scenario by considering f RCMF = 0 in this case).Accordingly, the mass of each element e in the system decreases by: . 4 e e e syst RCMF liq RCMF crys syst Notice that in the fractional crystallization case, X e crys represents the content e in the crystals that form and settle during the current crystallization step, and can be determined in the same way as for equilibrium crystallization (namely, Equation (7) below).

f O2 in the Magma Ocean
During the core formation event associated with mantle melting and MO creation, equilibrium between Fe-rich metallic melt and Fe-bearing silicate melt (Reaction 5) sets the f O 2 of the metal-bearing MO to ΔIW = −2 (Wood et al. 2006;Rubie et al. 2011), where ΔIW stands for the difference between the actual f O 2 and that set by the iron-wüstite buffer, in log-units.The f O 2 value corresponding to ΔIW = −2 is thus used to compute the initial Fe 3+ /ΣFe in the MO silicate liquid, from the equilibrium between ferrous FeO and Fe 2 O 3 (Reaction 6).
The equilibrium constant of Reaction (6) (given by Equation (A1) of Hirschmann 2022) depends on pressure and temperature (and composition) (Hirschmann 2012;Armstrong et al. 2019;Deng et al. 2020;Hirschmann 2022), so the effective P-T conditions of core-mantle equilibration need to be assumed to calculate the Fe 3+ /ΣFe value corresponding to ΔIW = −2.

R e t r a c t e d
We consider two end-member cases: (1) the "reduced" case, in which metal-silicate equilibration occurs at the surface of the MO, simulating the scenario of an MO previously well drained of its metallic iron, with a small impactor emulsifying and supplying metallic melt at near-surface conditions.In this case, ΔIW = −2 at the surface correspond to Fe 3+ /Σ Fe = 0.05 for fractional crystallization, and Fe 3+ /Σ Fe = 0.04 for equilibrium crystallization.
(2) the "oxidized" case, in which we set ΔIW = −2 at the base of the MO at 55 GPa (Deng et al. 2020), i.e., where final ponding of metallic melt could have taken place.This yields Fe 3+ /Σ Fe = 0.11 for fractional crystallization, and Fe 3+ /Σ Fe = 0.08 for equilibrium crystallization.These various end members for core-mantle equilibration conditions and associated redox state are illustrated in Figure A1.
Using the same parameterization (Equation ( A1)), we then reverse the reasoning and compute the subsequent evolution of the f O 2 profile in the MO from the evolution of Fe 3+ /ΣFe (assumed homogeneous in the MO), as the MO crystallizes.To do so, we track the independent fractionation of Fe 3+ and Fe 2+ between solid and melt, using Equations (3) and (4), and their partition coefficients definition: where e is either FeO or FeO 1.5 , and D e is their respective crystal-silicate melt partition coefficient, specified below.The calculation of initial mass fraction of each oxides from initial molar Fe 3+ /ΣFe and the bulk silicate Earth (BSE) total iron oxides mass fraction is given in Appendix A.3.
Ferric iron oxide has a partition coefficient dependent on P, calculated based on those for the individual mineral-melt partition coefficients of crystallizing minerals (based on the crystallization sequence of the 2000 km deep MO from Elkins-Tanton 2008, truncated to 55 GPa; see Appendix A.3).
We consider a partition coefficient value of 1 for ferrous iron oxide.While ferrous iron oxide is incompatible in most minerals forming the crystallization sequence, the increase of the ferric iron mass fraction only yields very high total iron oxide mass fractions (>20 wt%) in the late stages of fractional crystallization, which would be yet amplified by considering ferrous iron oxide incompatible.This is a natural feature of the fractional crystallization case, which leads to maximum partitioning and is meant to be an end-member case rather than a realistic description.In the equilibrium crystallization case, a value of + D Fe 2 close to 1 (>0.9)leads to a final iron oxides mass fraction of 12 wt%, similar to that of evolved basaltic magmas on Earth.Therefore, for the sake of simplicity, we adopt a partition coefficient of unity, thus amplifying the effect of differential fractionation of Fe 2+ and Fe 3+ .We do not expect, however, the use of a more realistic value for + D Fe 2 to significantly affect our results.
In the fractional crystallization case, where fractionation is much more important, doing so yields an unrealistically high final mass fraction of total iron oxide.We thus set the partition coefficient of FeO to 1, and let the total iron oxide mass fraction increase up to a maximum of 20% due to the sole effect of FeO 1.5 fractionation.While this might artificially increase the effect of iron oxides differential partitionning on MO redox evolution, this value is not far from that used in the equilibrium crystallization case, suggesting that it still captures a realistic behavior of FeO.
Knowing the evolution of the mass fraction of both oxides (and thus of the molar Fe 3+ /ΣFe, following the inverse pathway of Appendix A.3), we use Hirschmann's (2022) parameterization to calculate the corresponding f O 2 at the MO potential temperature and surface pressure, i.e., the conditions at which further equilibria considering f O 2 are assumed to take place.

Redox Speciation in the Atmosphere and Volatiles Solubility
At the high temperatures met at the surface of an MO (>1500 K), methane is only present in negligible amounts.Thus the C-O-H system is governed by the two following redox equilibria in the atmosphere: In this section, e refers to one of the volatile elements (H or C) while s refers to one of the molecular species (H 2 O,H 2 ,CO 2 , or CO; O 2 being excluded from this section's considerations unless explicitly stated).The fugacities of the different gaseous species in the atmosphere are related through the value of the reaction constants of equilibria in Equations ( 8) and (9), and the oxygen fugacity: where f s is the fugacity of species s involved in the reaction (including O 2 , whose fugacity calculation was presented in the previous section, unlike other gases), c s is its corresponding stoichiometric coefficient (signed: positive for products and negative for reactants), K is the equilibrium constant, and T is the effective temperature of the equilibrium.We consider that the effective temperature at which equilibrium is achieved (i.e., the temperature at which K and f O 2 are evaluated) is the potential temperature of the MO.This has the advantage of making atmospheric speciation a state function of T pot , making the comparison between different models easier.This choice is discussed in Section 3. We consider the fugacities to be equal to the partial pressures (expressed in bar) and replace f s with p s in the following.
From the partial pressures of the gaseous species, we calculate their content dissolved in the melt X s liq (expressed as mass fraction) using species-specific solubility laws taking the general form: where α s and β s are the solubility law coefficients of species s (values listed in Table 5).We consider volatiles to be species to be perfectly incompatible (i.e., having a partition coefficient equal to zero), and as a result, = X 0 s sol holds at all time.The total atmosphere mass is: where R p and g are the planetary radius and gravity, respectively, and p atm is the total atmospheric pressure.The mass fraction of species X in the atmosphere is expressed as a 6 The Planetary Science Journal, 4:31 (22pp), 2023 February Maurice, Dasgupta, & Hassanzadeh R e t r a c t e d function of its partial pressure as: where μ s is the molecular weight of species s, and μ atm is the average molecular weight of the atmosphere.Finally, as an element e can be carried by several species s (e.g., H 2 O and H 2 for H), its mass inventory (Equation ( 3)) can be written as a function of the partial pressures of gaseous species carrying this element, using Equations (11), ( 12), and (13) and summing over all concerned species: , 14 where μ e is the atomic mass of element e and l s e is the number of atoms of e in a molecule of s.
We solve Equation (10) two times (once for Equilibrium (8) and once for Equilibrium (9)) and Equation (14) two times (once for conservation of H and once for conservation of C), resulting in four equations to solve for the partial pressures of the four gaseous species (O 2 not included).These equations are solved using the "hybrid" root finding method implemented in the "root" function of scipy (Virtanen et al. 2020).Contrarily to H and C, O is not considered conserved, but the partial pressure of O 2 (i.e the f O 2 ) is buffered at the value imposed by Equilibrium (6).This assumption only breaks down when the masses of the MO and of the atmosphere become comparable, which is largely past the 99 wt% MO solidification that we use as our final state.The remaining MO's mass is then 2.32 × 10 22 kg, i.e., the equivalent of a 4457 bar, pure-CO 2 atmosphere, and higher pressure for a mixed atmosphere.This is a factor 2 larger than the final pressure of our most volatile-rich case.

Thermal Evolution Model
The thermal evolution model is used to compute the time evolution of the MO crystallization and associated atmospheric outgassing.It consists of two main parts: heat conservation in the planet (including core, solid mantle and MO), and radiative heat flux through the greenhouse atmosphere, the latter providing a boundary condition to the former.Notations in use in Section 2.3 are listed in Table 3.

Planetary Heat Conservation
The general form of the heat conservation equation is as follows: where E th is the thermal energy of the planet (the volumeintegral of ρc p T), L is the latent heat of crystallization of silicate, M sol is the mass of solids in the whole mantle, H int is the internal heating, and F p is the heat flow through the planetary surface (defined positive outward).The lefthand-side term of Equation (15) corresponds to the secular cooling of the planet.The first term on the left-hand side corresponds to the release of latent heat associated with MO crystallization, the second term corresponds to the radioactive decay of heat-producing elements, and the third term corresponds to the heat flow out through the planetary surface.
An accurate parameterization of the secular cooling term involves resolving various complex heat transfer mechanisms within the young Earth, which are, for some, not well understood (e.g., the heat flux at the melting interface between the MO and solid mantle; Labrosse et al. 2018;Agrusta et al. 2020).For the sake of simplicity, here we consider that the thermal energy of the planet decreases linearly with T pot between the initial and final thermal states (see Appendix B.1).The secular cooling term can thus be written: , where dE th /dT pot is a constant (Table 3).This approximation amounts to considering that the whole interior of the Earth convects, initially along a single adiabat, and a boundary layer progressively builds up at the CMB, offsetting the adiabat in the core, preventing it from cooling below the present-day estimated temperature (Andrault et al. 2016).The final thermal state is similar to (if not colder than) the present-day geotherm and hence likely underestimates the temperature in the Earth at the end of MO crystallization.We thus overestimate the secular cooling term, and in turn, predict an upper bound for the MO lifetime.
M sol depends only on T pot , so we can apply the chain rule and write dM sol /dt = dM sol /dT pot × dT pot /dt.Contrarily to the secular cooling term, dM sol /dT pot is not constant but varies as an expression of the shape of the melting curves.We do not account for latent heat release associated with iron crystallization in the core, as inner core crystallization most likely did not start before the end of the MO crystallization (Bono et al. 2019; an agreement with our final thermal state).

R e t r a c t e d
The internal heating term corresponds to the heat produced by the decay of the long-lived radioactive elements 235 U, 238 U, 232 Th, and 40 K, distributed in the whole mantle, with abundances indicated in Table 4.2 of Turcotte & Schubert (2014), calculated at 100 Myr after CAIs formation (our estimate for the Moon-forming impact; Maurice et al. 2020).Considering the short lifetime of the MO compared to the halflives of these isotopes, we neglect their decrease with time due to radioactive decay and consider H int as a constant (Table 3).The heat flux through the MO is the focus of Section 2.3.2.
Equation (15) can thus be recast as the evolution of the MO's potential temperature: The denominator of Equation ( 16) is always positive (since dE th /dT pot > 0 and dM sol /dT pot < 0), so the potential temperature decreases with time as long as the heat flux in the MO dominates internal heating.

Convective Heat Flux through the MO
Following Lebrun et al. (2013), the heat flux through the planetary surface (F p ) is obtained by equating the convective heat flux from the MO (F conv ) with the radiative heat flux through the atmosphere (F rad ).Using results from boundary layer theory (Siggia 1994), the former is parameterized as a function of the temperature drop across the thermal boundary layer at the top of the MO (Lebrun et al. 2013): where T sfc is the temperature at the surface of the MO, k and d are the thermal conducitivity in the MO and its depth, respectively, and Ra is the Rayleigh number of the MO, defined as: where η is the viscosity of the MO, and κ = k/(ρc p ) its thermal diffusivity.Following Nikolaou et al. (2019), η depends on T following a Vogel-Fulcher-Tammann equation, further affected by the melt fraction in the equilibrium crystallization case (see Appendix B.2 for the expression of η).As ρ, η, α, and c p are not homogeneous in the MO, Ra is calculated using average values for those quantities (arithmetically averaged for ρ, α, and c p , and harmonically averaged for η, following Nikolaou et al. 2019).At high Rayleigh numbers, another convection regime is theorized (referred to as the "hard" turbulence regime, as opposed to the "soft" turbulence regime described by Equation (17)), which follows a different scaling: the study, we discuss the effect of changing to the hard turbulence scaling in Section 4.4.The radiative flux F rad is also a function of T sfc (see Section 2.3.3),so we solve F conv (T sfc ) = F rad (T sfc ) − F Sun (where F Sun is the shortwave flux from the Sun) for T sfc using the brentq method implemented in scipy (Virtanen et al. 2020; function root_scalar), and use the value at which these fluxes coincide for F p .

Convective-radiative Atmosphere
We compute the radiative heat flux through an atmosphere composed of H 2 , H 2 O, CO, and CO 2 .Similar to previous studies (Marcq et al. 2017;Katyal et al. 2019), the atmosphere consists (from the surface to the top of the atmosphere) of a dry troposphere, a moist troposphere (where water condenses), and a stratosphere.In the dry troposphere, the temperature profile is given by the dry adiabat: where p is the gas (hydrostatic) pressure, and cp is the average molar heat capacity of the gas.When the saturation pressure of water at the ambient temperature reaches the partial pressure of water, condensation occurs and the temperature gradient is given by the moist pseudo-adiabat (assuming complete rain-out of the condensate; Pierrehumbert 2010): ( ) where p a is the sum of the partial pressures of all gases other than water (dry background air), c p,a is the molar heat capacity of the dry background air, and L H O 2 and c p,H O 2 are the latent heat of vaporization and molar heat capacity of water, respectively.Equation (21) gives the slope of the (pseudo) adiabat (lapse rate) as a derivative in p a rather than total pressure.The system is closed by the Clausius-Clapeyron relation, which gives the partial pressure of water as a function of temperature (assuming that water is always at saturation in the moist troposphere; see Appendix B.3 for the expression of the saturation pressure).Equations (20) and (21) are solved using Runge-Kutta fourth-order total-and dry air pressurestepping, respectively, with T sfc as boundary condition.
Finally, following previous studies (Marcq et al. 2017;Katyal et al. 2019), we consider an isothermal stratosphere at 200 K extending from the pressure level where the tropospheric temperature reaches 200 K up to the top of the atmosphere (defined at 1 Pa).While we assume that water is the only condensing species, we verified a posteriori this assumption by comparing the partial pressures of all other gases present with their respective saturation pressures at ambient temperature, and found that in all cases water is indeed the only saturated gas.
The resulting temperature profile and mass mixing ratios are then used to compute the bolometric flux (outgoing longwave radiation, OLR) using the radiative transfer code petitRAD-TRANS (Mollière et al. 2019).The heat flux is computed as the wavelength-integrated flux (between 0.1 and 251 μm) using the correlated-k method (Lacis & Oinas 1991) with a coarse wavelength binning (10 bins).The coarse binning allows for rapid calculation, as we do not need the wavelength-resolved spectrum.We also account for the collisional opacity of H 2 , as it has been shown to be significant for high partial pressures of H 2 in MO atmospheres (Lichtenberg et al. 2021).

Results
We calculate the evolution of the f O 2 in the MO, the corresponding redox speciation of the outgassed C-O-H atmosphere, and the thermal evolution of the MO cooling through this greenhouse atmosphere.We investigate the two crystallization scenarios presented in Sections 2.1.1 and 2.1.2and the two initial conditions on f O 2 discussed in Section 2.2.2.The initial volatile budget of the terrestrial MO is poorly constrained.Marty (2012) and Hirschmann (2018) provided estimates for the BSE volatile content (11.6 and 2 EOequivalent of H, respectively, and both find C/H ∼ 1.7), but it is still debated whether these volatiles were delivered during the late veneer (Albarède 2009) or were already present in the building blocks of the Earth (Piani et al. 2019;Grewal et al. 2021).In the former case, estimates of the BSE volatile content do not reflect the content in the terrestrial MO, while in the latter case they provide a lower bound.Catling & Kasting (2017) suggested that, accounting for the time evolution of the Sun's extreme-UV (XUV) flux, if energy-limited escape had been operating throughout Earth's life (which is not realistic but provides an upper bound of the escape flux), it would have depleted about 10 Earth oceans (EOs)-worth of H.As this provides a generous overestimate of the amount of H lost, here we vary the H budget of the BSE budget between 1 and 9 EOs of water.In order to cover a large parameter space, we vary the mass C/H ratio between 0.5 and 2. Notice that the volatile budgets are indicated as the BSE contents, i.e., contents in the combined (initial) MO and unmolten lower mantle, assuming a homogeneous distribution (as can be expected if volatiles where delivered during Earth's accretion).Thus, since the MO represents initially ∼58 wt% of the mantle, the actual MO volatile content is in the same proportion (e.g., 0.58 EO-worth of H content for the minimum H end member).The volatile content in the lower mantle is not relevant to the present study.

Evolution of f O 2 during Magma Ocean Crystallization
We calculate the evolution of Fe 3+ /ΣFe in the liquid throughout MO crystallization for both crystallization scenarios, and both initial f O 2 conditions (Figure 2(a)).The different initial conditions on f O 2 translate to different initial conditions for Fe 3+ /ΣFe (see Section 2.2.2 and Appendix A.2).The difference in initial Fe 3+ /ΣFe between the two crystallization scenarios is due to the difference in initial potential temperature, and thus the P-T conditions at which ΔIW = −2 (see Figure A1(b)).The higher incompatibility of Fe 3+ relative to Fe 2+ leads to an increase of Fe 3+ /ΣFe in the liquid as the MO crystallizes for all cases.Fractionation is stronger for fractional crystallization because no melt is retained in the solid cumulates, while they retain a fraction f RCMF in equilibrium crystallization.Thus, a significantly higher Fe 3+ /ΣFe is reached at the end of fractional crystallization cases.
Using liquid's Fe 3+ /ΣFe, we compute the surface f O 2 (evaluated at T pot ) during MO crystallization (Figure 2(b)).Sossi et al. (2020) posited that the temperature dependence of f O 2 as set by Reactions (5) and (6) are similar, which would dismiss the influence of temperature on ΔIW.However, Hirschmann (2022) argued that taking into consideration the difference in heat capacity between FeO and FeO 1.5 , Reaction (6) has an increased temperature dependence past 2000 K.As a matter of fact, under a pressure of 1 bar and for a fixed Fe 3+ /ΣFe, ΔIW as calculated from Hirschmann's (2022) parameterization increases by ∼3 log-units upon decreasing the temperature from 3000 to 1500 K (roughly the range addressed in the present study).Thus, the trend observed in Figure 2  On both panels, solid curves correspond to fractional crystallization and dashed curves to equilibrium crystallization, while orange curves correspond to an initial ΔIW = −2 at the surface (reduced cases) and blue curves to an initial ΔIW = −2 at the bottom of the MO (resulting in a higher surface f O2 , hence oxidized cases).Notice that the fraction of solidified MO (x-axis) is different from the crystal fraction in the case of equilibrium crystallization because the former includes the melt trapped in the solid cumulates while the latter does not.9 The Planetary Science Journal, 4:31 (22pp), 2023 February Maurice, Dasgupta, & Hassanzadeh R e t r a c t e d due to both increase in Fe 3+ /ΣFe as well as the decrease in potential temperature as the MO crystallizes.
The effects of the temperature and of Fe 3+ /ΣFe on f O 2 act in the same direction and yield a final increase of ΔIW of several log-units in all cases.The equilibrium crystallization cases (dashed curves in Figure 2) exhibit the lowest increase in Fe 3+ /ΣFe, and most of their increase in ΔIW is accounted for by the decrease in T pot .Conversely, in the fractional crystallization scenario (solid curves in Figure 2), Fe 3+ /ΣFe reaches high values (FeO 1.5 becoming the major iron oxide) and induces a highly oxidized final state.

Atmospheric Speciation and Magma Ocean Volatile Content
The oxidation state of the MO obtained in the previous section is used to calculate the composition of the outgassed atmosphere.Three main effects govern the evolution of atmospheric pressure and composition during MO crystallization.
First, as both H and C are incompatible elements, they are partitioned into silicate melts during crystallization.Thus, MO solidification increases their concentrations in the residual melt.This increase is accommodated by outgassing until the partial pressures of the volatile-bearing species in the atmosphere and their concentration in the silicate melt reach an equilibrium given by solubility laws (Equation ( 11), assumed to be an instantaneous process).Second, when an element exists in different species in the gaseous phase (e.g., oxidized or reduced), the evolution of its outgassing is further complicated by the different solubilities of the gaseous species.For instance, in the H-bearing redox system, H 2 O has a significantly higher solubility in silicate melts than H 2 , allowing for more H to dissolve when the f O 2 in the MO increases, all other things being equal.Third, when several species are present in the atmosphere, variations in the partial pressure of one of them influence the others: the average molecular weight in the atmosphere is altered, which modifies the relationship between the mass of one species in the atmosphere and its partial pressure (as per Equation ( 14)).If the average molecular mass of the atmosphere decreases (for instance if H 2 O is outgassed in a CO 2 -dominated atmosphere), the partial pressures of other species will decrease to accommodate the same mass in the atmosphere (e.g., Bower et al. 2019).Actually, this effect is associated with a mass redistribution between the atmosphere and the MO, since a modification of the partial pressure induces a modification of the solubility and hence of the dissolved mass.This leads to the surprising configuration where a reduction of the partial pressure of a gas is correlated with its net outgassing, because its solubility is accordingly decreased.
How these three effects operate is controlled by specific features of the H and C systems, in particular the difference in solubility between oxidized and reduced species and the temperature dependence of the equilibrium constant of their respective redox reactions.For a BSE content of 3 EOs-worth of H and a C/H ratio of 1, undergoing fractional crystallization, with an initial f O 2 buffered at IW-2 at the surface (reduced case), Figure 3(a) represents the evolution of the partial pressures of the different volatile-bearing species in the atmosphere, and Figure 3(b) shows the solubilities of H and C. We chose this case to illustrate the evolution of the volatile speciation because its f O 2 crosses the stability domains of most species, contrary to the more oxidized cases.
In both the C and H redox systems, the oxidized species has a higher solubility than the reduced one.Thus the increase in f O 2 with MO crystallization reported in the previous section competes with the partitioning of volatiles in the remaining melt: the former reduces outgassing (first effect mentioned above) while the latter promotes it (second effect mentioned above).However, in the case of C, both species have a fairly low solubility, and most C has been readily outgassed since the beginning of MO solidification.The effect of redox speciation on C solubility-as accounted for by our model-is thus marginal.The possibility of graphite saturation (that would reduce C solubility and sequester a part of the C budget) is discussed in Section 4.2.The solubility difference between oxidized and reduced species is much larger in the case of the H system, where water has a solubility in silicate melt higher than molecular hydrogen by several orders of magnitude (Hirschmann 2016).Even at the beginning of the crystallization, when the MO is the most reduced and H 2 is 10 times more abundant than water in the atmosphere, it represents only ∼10% of the total dissolved H (Figure 3(b)).Thus, the high solubility of water compensates for the unfavorable redox conditions and the effective solubility of H is set by water at all times.
Contrarily to other species whose solubility follow Henry laws, H 2 O's partial pressure increases as its solubility to the power . Thus, partitioning of water in the remaining MO leads to a sharper increase in p H O 2 compared to p CO 2 (the other species favored by redox evolution) in the late part of the crystallization.In particular, following the third effect mentioned above, the delayed outgassing of water decreases the partial pressure of CO 2 in the last 20% (see flattening and even final decrease of the solid orange curve in Figure 3(a)).
The initial BSE abundances of C and H influence the evolution of the MO outgassing.Figure 4 represents the evolution of the molar mixing ratios of the atmospheric species throughout the MO crystallization over all of the combinations of BSE H budgets and C/H ratio, for a "reduced" MO undergoing fractional crystallization.Increasing the C and H mass budgets by the same amount (i.e., keeping C/H constant) induces a significant increase in H 2 ʼs molar mixing ratio.This is due to the small molecular mass of this species: the same increase in mass in the atmosphere corresponds to a larger increase in partial pressure for light species than for heavy ones.Perfect fractional crystallization being assumed in this case, all volatiles are eventually outgassed.Since all cases represented follow the same f O 2 evolution where the final redox state is oxidized enough for reduced species to be negligible, they all converge toward similar final atmospheric compositions, where the relative molar fractions only depend on C/H.
As already pointed out in the case of iron oxides partitioning, enrichment in incompatible elements in the MO is much lower in the case of equilibrium crystallization compared to fractional crystallization (Figure 2(a)).This effect is even stronger for volatiles, which we assume to be perfectly incompatible.It follows that volatiles are significantly less outgassed in the former case, and a large amount is retained in melts trapped in the solid cumulates.This applies particularly to H, due to the high solubility of H 2 O. Furthermore, the lower final f O 2 reached in the case of equilibrium crystallization induces a more reduced final atmosphere (although still dominated by H 2 O and CO 2 for all cases).

R e t r a c t e d
Figure 5 represents the evolution of the atmospheric composition for both crystallization scenarios and both initial redox states, for a case with the same volatile inventories as that presented in Figure 3.The late outgassing of water induces the final decrease in the average molecular mass in the atmosphere (black curves), observed in all cases.However, the final mixing fraction of water is much smaller in the equilibrium crystallization scenario, because a large quantity remains sequestered in the trapped melts.This translates into a milder final decrease of the molecular mass in the atmosphere, as well as a dichotomy in the net change in atmospheric molecular mass over MO solidification: negative in the case of fractional crystallization (because H 2 O, the dominant gas in the final atmosphere, is lighter than CO, the dominant gas in the initial atmosphere), and positive in equilibrium crystallization cases, whose final atmospheres are dominated by heavy CO 2 .
An oxidized initial state leads to a drastic reduction of H 2 and an increase of CO 2 (H 2 O staying largely dissolved in the MO), and thus an increased solubility of H compared to C. Because C-bearing gases are heavier than H-bearing ones, this translates into an increased average molecular mas of the atmosphere (compare black curves in each column).

Timescale of Magma Ocean Crystallization
Previous studies found that MO outgassing of greenhouse species (in particular H 2 O and CO 2 ) increases the opacity in the atmosphere, thereby decreasing the cooling flux of the MO and prolonging its lifetime (Abe & Matsui 1985, 1988;Elkins-Tanton 2008;Lebrun et al. 2013;Salvador et al. 2017;Bower et al. 2019;Nikolaou et al. 2019).As shown in the previous section, redox speciation in the atmosphere influences volatiles solubilities, and thus atmospheric outgassing.Furthermore, the different species of the redox couples have different radiative as well as thermodynamic properties, affecting the infrared opacity and lapse rate of the atmosphere.Therefore, redox speciation in the atmosphere further influences the time evolution of MO crystallization.
Figure 6 represents the time evolution of the total atmospheric pressure for all investigated values of the C/H ratio, keeping a BSE H budget of 3 EOs (panel (a)), as well as for all of the BSE H budgets envisaged, keeping C/H = 1 (panel (b)).Among the oxidized species, H 2 O is a stronger infrared absorber than CO 2 .Among the reduced ones, H 2 and CO are only mild greenhouse gases, but Lichtenberg et al. (2021) found that H 2 collision-induced absorption might be significant.It may thus be anticipated that H will be more influential than C in setting the MO lifetime, and hence varying the H budget will be crucial.This is verified in Figure 6(a), where increasing the H budget from 1 to 9 Earth ocean masses of water prolongs the MO crystallization by ∼1 order of magnitude.We notice that accounting for the H 2 -H 2 collision-induced absorption does not significantly prolong the MO lifetime, even in cases where the partial pressure of H 2 is the highest (∼145 bar).Lichtenberg et al. (2021) pointed out the important thermal blanketing effect of a thick H 2 atmosphere, but for a higher H 2 pressure (∼200 bar).The influence of C/H (Figure 6(b)) also follows an intuitive trend, whereby increasing the C/H ratio while keeping the H inventory constant (i.e., increasing the C inventory) translates into a prolongation of the MO lifetime.This prolongation is only marginal compared to that induced by increasing the H budget. Furthermore, as was already observed, fractional crystallization of the MO allows for a much more extensive outgassing than equilibrium crystallization and thus, in turn, for a slower cooling.Indeed, we consistently observe that, all other

R e t r a c t e d
things being equal, the former (solid curves) yields a significantly longer MO than the latter (dashed curves) in Figure 6.Over the entire parameter space covered, the shortest MO (initially reduced MO undergoing equilibrium crystallization, with a BSE H budget of 1 EO and C/H = 0.5) solidifies in 67 kyr, and the longest MO (initially oxidized MO undergoing fractional crystallization, with a BSE H budget of 9 EOs and C/H = 2) solidifies in 1.530 Myr.One noticeable feature in Figure 6 is the consistently longer lifetime of initially reduced MOs (thick lines) compared to initially oxidized ones (thin lines).This effect seems to contradict the stronger infrared absorption of oxidized gases compared to reduced ones.Indeed, Katyal et al. (2019) found that, for a fixed snapshot in the MO evolution, reduced cases exhibit larger OLR than oxidized ones, which should lead to their shorter lifetimes.The surface temperatures for which they calculated the OLR correspond to the beginning (3000 K) and the end (1500 K) of MO crystallization (these temperatures are roughly similar to those in the present study).As described below, we do observe that the OLR of oxidized cases becomes larger than that of reduced ones for intermediate to high temperatures, but it is similar for such end-member temperatures.Bower et al. (2022) also calculated the OLR for both reduced and oxidized MO atmospheres, and found the former to crystallize faster than the latter (in agreement with Katyal et al.'s 2020 results).However, the simple atmospheric model they used (from Abe & Matsui 1985) cannot capture the influence of atmospheric chemistry on the temperature profile, which we show is key to the effect responsible for the stronger radiation of oxidized atmospheres.
The radiative heat flux through the atmosphere is the result of species-specific (wavelength-dependent) absorption features and temperature-controlled emission of all pressure levels.For a fractional crystallization case (BSE H budget of 3 EOs, C/H = 1), Figure 7 represents the temperature and OLR contribution profiles in the atmosphere at three successive snapshots of MO crystallization: (1) beginning of crystallization, (2) after 75 wt% of the MO has crystallized (corresponding to a surface temperature of 2000 K, represented by white circles in Figure 6(b)), and (3) end of MO crystallization.The OLR contribution is the normalized, wavelength-integrated contribution of each pressure level to the OLR.Generally, at low (subsolidus) temperature, emission from the deepest layers of thick (>10 bar) atmospheres is fully absorbed, and only shallower layers' emission effectively contributes to the OLR.However here, the extremely hot surface of the MO does contribute significantly to the initial OLR, despite the relatively thick atmosphere.Figure 7(a) shows that the surface of the MO accounts for 75% of the OLR at the beginning of MO crystallization.
As crystallization proceeds, the surface contribution diminishes, and the bulk OLR contribution is offset toward higher levels (Figure 7(b)).While the height of the radiating levels is controlled by the absorption of overlying layers (and thus by the composition of the atmosphere), the intensity of their radiation is governed by their temperature, which in turn is set by the lapse rate in the troposphere (the stratospheric contribution always being negligible).Among the gaseous  .Although the bulk of the OLR-contributing pressure levels lies at slightly lower p in the oxidized case, their temperature is higher than in the reduced case due to the difference in the lapse rate, resulting in lower OLR of the latter case.
Eventually, as the MO crystallizes, both cases evolve toward a strongly oxidized atmosphere, and they become almost indistinguishable.At the same time, the surface temperature further cools down and the atmosphere becomes thicker, so the contribution from the surface vanishes, and the main contributing layers spread over a larger pressure range, mostly in the moist troposphere (Figure 7(c)).The stronger OLR exhibited by oxidized atmospheres is consistent throughout the parameter space explored, resulting in a general shorter lifetime of initially oxidized MOs.

Volatiles Retention in the MO Cumulates
In the case of fractional crystallization, volatiles are virtually entirely outgassed (since we neglect partitionning of volatiles in solids).In the equilibrium crystallization cases, however, the large amount of interstitial melt trapped within the solid mantle allows for sequestration of a significant mass of volatiles.We expect chemical speciation to play a role in the retention of volatiles in the mantle, by indirectly setting the effective elemental solubilities.H in particular can be retained due to the high solubility of water.Figure 3 showed that even in the reducing case, H solubility is set high by water.As a result, a large quantity of H is retained in the interstitial melt trapped in the mantle.In contrast, both CO and CO 2 have a low solubility, and C is mostly outgassed since the beginning of the MO crystallization.However, CO 2 has a solubility approximately twice that of CO, so the redox speciation must affect whatever little C is retained in the mantle.
Figures 8(a) and b present the fraction of the initial MO volatiles mass (i.e., 58% of the indicated BSE budget) that is retained in the solid cumulates of the MO over the whole parameter space (for equilibrium crystallization).We observe that for all cases, more than 70% (and up to ∼93%) of the H, and from ∼1.8% to ∼5.2% of C, are retained in the cumulates.Increasing the bulk quantity of H has an opposite effect to increasing the C budget on the retention of both elements in the MO cumulates.Increasing the bulk H budget decreases the final relative amount of both H and C retained.This is due to the corresponding increase in H 2 in the atmosphere reported in Section 3.2 (see Figure 4), which promotes the outgassing of other species by lowering the average molecular mass of the atmosphere.On the contrary, increasing the bulk inventory of C leads to an increased volatiles retention in the MO cumulates because both CO and CO 2 are heavier than H 2 O (and H 2 ), so increasing the bulk C (which is largely partitioned in the atmosphere) increases the average molecular mass of the atmosphere, deferring outgassing.Figures 8(c) and (d) present the final H and C contents of the MO cumulates.The small decrease in the fraction of retained volatiles with increasing BSE volatiles budget is largely offset by the absolute increase in BSE volatile budget.The H content in the cumulates increases from ∼32 to ∼270 ppm as the BSE H budget increases from 1-9 EO.Similarly, the retained fraction of C increases from less than 1 ppm up to ∼33 ppm.The oxidation state also plays an important role for C. As already pointed out, water sets the effective solubility of H, even in the most reduced phase of MO crystallization.For C, however, CO 2 solubility is about twice that of CO, so oxidized conditions promote C retention.Furthermore, the transition between a CO-dominated and a CO 2 -dominated atmosphere occurs after approximately 80% of the MO crystallized in reduced cases, and after only 40% in the oxidized case (see Figure 2(b), when IW ∼ 1).Hence, the evolution of the MO toward oxidized conditions accentuates the favorable effect of f O 2 on C retention.Overall, the final quantity of C in the mantle is multiplied by 2 upon switching from the reduced to the oxidized case, while the retention of H increases by at most 7%.
While we consider an RCMF value of 0.4, the volume fraction of trapped melt in the cumulates may be lower, either because of an intrinsically lower value of the melt fraction at which the dynamical regime transition occurs, or because of melt extraction from the cumulates pile by percolation (Hier-Majumder & Hirschmann 2017).The amount of volatiles retained in the solid mantle scales with the the amount of trapped melt, i.e., the RCMF.Adopting a lower value of the RCMF leads to enhanced fractionation (tending toward the fractional crystallization case as the RCMF becomes vanishingly small), and thus strengthens the increase in f O 2 as the MO crystallizes.Although the increase of f O 2 at cumulates crystallization associated with a lower RCMF can affect the solubility of volatiles, based on our results, we expect the MO cumulates retention to be roughly linear with RCMF value.We find that decreasing the RCMF to 0.2 consistently halves the amount of C in the mantle throughout the parameter space, while the amount of H retained is multiplied by ∼0.7, without affecting the trends described above.The nonlinear behavior of H retention is attributable to the nonlinear solubility of water, which is largely setting H solubility.
These results confirm the low efficiency of H outgassing associated with MO crystallization (Miyazaki & Korenaga 2022a;Bower et al. 2022), and suggest that, beyond Earth, terrestrial planets and exoplanets that accreted from volatile-rich material start from a wet interior post-MO stage.This can have drastic consequences on their long-term mantle dynamics, as water significantly decreases both the melting temperatures and the viscosity of silicate rocks.If a plate tectonic regime starts directly after MO crystallization, this could promote convection and thus result in a high plate velocity, accelerating the drawdown of atmospheric CO 2 and bringing the planet out of the runaway greenhouse state and toward habitable conditions (Miyazaki & Korenaga2022a, 2022b).If the planet is in the stagnant lid regime, a wet mantle Figure 6.Time-series of the total pressure for all H BSE budgets keeping C/H constant at 1 (a), and for all C/H values keeping the BSE H budget constant to 3 EOs (b).Panel (c) is a zoomed inset of the cases having the a BSE H content of 1 EO.Thick and thin lines represent reduced and oxidized cases, respectively, while solid and dashed lines correspond to fractional and equilibrium crystallization cases, respectively.Orange curves coincide on both panels (notice that time axes differ between the panels).White circles correspond to the 75 wt% MO solidification, reported in the middle panel (b) of Figure 7 (panels (a) and (c) correspond to the origin and the end of the corresponding lines, respectively).would lead to abundant melting.However, the latter may not necessarily translate into efficient outgassing if the atmosphere is thick enough to increase the solubility of water in basatic melt (Gaillard & Scaillet 2014;Tosi et al. 2017;Godolt et al. 2019).Thus, the MO-outgassed atmosphere could prevent further outgassing of stagnant lid planets.

Graphite Saturation
The volatile content of the system can be further altered by processes that have been neglected in the model.The low f O 2 reached in the depth of the MO in the reduced cases (down to IW − 4; see Figure A1(a)) decreases the CO and CO 2 concentrations at graphite saturation (Eguchi & Dasgupta 2018;Yoshioka et al. 2019), making precipitation of C possible.If the concentration of either (or both) species, as predicted by our model, increases above saturation, graphite (or diamond) would be produced to buffer the MO concentration at this value.Depending on the buoyancy of the precipitate, it could either settle at the bottom and leave the system (similar to the silicate crystals in the fractional crystallization case), remain entrained in the MO (as the silicate crystals in the equilibrium crystallization case), or accumulate at the surface and form a floating graphite layer (in which case the atmosphere redox state might be controlled by the CCO buffer; Keppler & Golabeck 2019).However, we notice that all of our cases also contain H, which could form CH species (particularly at high pressure; Chi et al. 2014;Li et al. 2015;Gaillard et al. 2022a) whose solubility could offset C saturation.Gaillard et al. (2022a) suggested that above 35 GPa (corresponding to the most reducing domain of the MO), C-content at C saturation should be above 100 ppm in H-bearing silicate melts, which is largely above the C concentration in our most C-rich case.However, it is important to keep in mind that this result relies on parameterizations calibrated for much lower temperatures and pressures than those met in a deep MO.

Atmospheric Escape
Another process that has not been accounted for and could yet affect the H budget is atmospheric erosion.In the atmosphere, H 2 , the lightest species, is particularly vulnerable to thermal escape, triggered by the enhanced XUV flux of the young Sun.If H 2 is abundant in the exosphere (the region of the atmosphere from where it can readily escape), escape is limited by the energy provided by the XUV flux (energylimited escape) while if diffusion of H 2 into the atmosphere is slow (and/or the XUV flux is low), H 2 availability in the exosphere becomes the bottleneck setting the escape rate (diffusion-limited escape).The energy-limited case provides a theoretical maximum for the escape flux, when all of the incoming solar XUV is used for H 2 escape.Following Watson et al. (1981), the H 2 mass-loss rate in the energy-limited case is: where R XUV is the radius at which XUV energy is effectively deposited, ò is an XUV heating efficiency factor, F XUV (t) is the (time-dependent) received XUV flux,  is the gravitational constant, and M p is the mass of the Earth.The time-dependence in F XUV accounts for the increased XUV activity of the young Sun.Following Katyal et al. (2020), we calculate M H 2  taking R XUV = R p , ò = 0.5 and F XUV (t) = 5 × (4.5/t) 1.24 × 10 −3 where t is the time elapsed after CAIs formation (i.e., the age of the Sun) in gigayears, and F XUV is expressed in W/m 2 .Evaluating Equation (22) for a 100 Myr-old Sun, we find an energy-limited H 2 mass flux of 6.32 × 10 5 kg s −1 .Using this value, we can calculate a posteriori the maximum mass of H lost over the lifetime of the MO for all cases in our parameter space.We find that it represents at most 4% of the initial MO content, which is too small to have a significant impact on the chemical and/or thermal evolution of the atmosphere and the MO.We note, however, that this effect might be important for exoplanets orbiting M-dwarfs, which have a much larger XUV flux during their early stages (Luger & Barnes 2015;Barth et al. 2021).In such cases, escape of H results in a net oxidation of the MO+atmosphere system (Wordsworth et al. 2018), which would alter the equilibrium between ferrous and ferric iron oxides.Another case where atmospheric escape could be relevant is when the incoming stellar radiation is higher than the OLR limit imposed by water condensation.In this case, the planetary heat balance cancels out, and the MO no longer crystallizes.H 2 escape is then the main mechanism to alter the atmosphere, and increase the OLR until mantle eventually solidifies.This mechanism is sensitive to the stellar flux (hence the orbital distance), and for Sun-like stars, a threshold in the semimajor axis has been suggested to be located between Venus's and Earth's orbits (Hamano et al. 2013;Salvador et al. 2017).A similar mechanism operating around other stars could allow for observable, longlived MO exoplanets.

Turbulent Convection Scaling
The MO lifetimes obtained in Section 3 are too short the for atmospheric loss to be significant.However, Nikolaou et al. (2019) found that adopting the hard turbulence (Equation ( 19)) rather than the soft turbulence (Equation ( 17)) scaling for the convective heat flux results in an increased MO lifetime.The two essential differences between these two scalings are the dependence of the heat flux on the Rayleigh number and the Prandtl number.In particular, for a given Ra, the hard turbulence scaling yields a lower heat flux than the soft one as long as Ra 1.7 × 10 8 × Pr 3 , which is always the case in our fractional crystallization simulations, and holds until 97% crystallization in the equilibrium cases.Conversely, sustaining the same heat flux requires a higher Ra in the hard turbulence scaling.
For the same set of parameters (volatiles budget, initial redox state, and crystallization scenario), at a given stage of MO crystallization (i.e., a given T pot ), all quantities entering the definition of Ra are fixed, except T sfc , while the Prandtl number , for all equilibrium crystallization cases.We note again that in our model, the volatile retention in solid cumulates is in the form of trapped melts and not via incorporation in the crystal lattice.Colors distinguish C/H while linestyles distinguish initial redox state (solid for reduced, dashed for oxidized).
The Planetary Science Journal, 4:31 (22pp), 2023 February Maurice, Dasgupta, & Hassanzadeh R e t r a c t e d is completely determined.Thus, the lower heat flux obtained in the hard turbulence scaling can be offset by decreasing T sfc (i.e., increasing ΔT, the temperature drop through the thermal boundary layer).In turn, decreasing T sfc decreases the temperature in the atmosphere (where the composition is also fully characterized by T pot ), which decreases the radiative flux.This latter effect induces a negative feedback, and the convective and radiative fluxes are balanced at a slightly lower value than in the soft turbulence scaling, although the Rayleigh number is higher.This effect is illustrated for an MO undergoing equilibrium crystallization with 1 EO H in the BSE and C/H = 0.5, in Figures 9(a)-(c).A given position on the x-axis corresponds a given T pot .For the convective flux to balance the radiative flux, a higher Ra is needed when the hard turbulence scaling is adopted (dashed curves, Figure 9 We anticipate this effect to lose importance as the heat flux (and accordingly the Ra) diminishes, until vanishing when Ra ∼ 1.7 × 10 8 Pr 3 .We confirm it by sweeping through our parameter space: we show that the prolongation of the MO lifetime due to adopting the hard turbulence scaling is important when the volatile budget is low (i.e., the radiative flux is generally high), and lessens as the volatile budget increases, yielding a thicker atmosphere.Figures 9(d Overall, for equilibrium crystallization cases (where the atmosphere is in general thinner), adopting the hard turbulence scaling yields an increase in MO lifetime from 73.5% for the most volatile-depleted case (represented in Figures 9(a)-(c)), to 10.8% for the most volatile-rich case (represented in Figures 9(d)-(f)).For fractional crystallization, where outgassing is much more extensive, this increase is limited between 8.3% (most volatile-depleted case) and 1.3% (most volatile-rich case).All of these results are for initially oxidized cases; the lower heat flux resulting from initially reduced cases described in Section 3.3 yields slightly lower MO lifetime prolongations.While the effect of using the hard turbulence scaling is more important for thinner atmospheres, which are more vulnerable to atmospheric escape, the prolongations are still not large enough for the energy-limited H-escape to alter significantly the planetary volatile budget during the MO lifetime.
Finally, it should be added that the dependence of the heat flux on the planetary rotation has also been omitted.According to Solomatov (2015), it becomes significant when Ra drops below 10 19 , which only occurs during the few last wt% crystallization of the MO (see solid curves on Figures 9(a) and (d)).Thus, this effect will likely not affect the bulk of the MO solidification, although we notice that the last stages of MO As the rotation of the Earth after the Moon-forming impact is unknown, we defer the investigation of the effect of the rotation to future studies.

Concluding Remarks
By tracking the evolution of f O 2 in the MO during its crystallization and accounting for its influence on atmospheric speciation, we investigated the couplings between redox chemistry and volatile solubility, and their impact on MO outgassing and greenhouse effect.We reached the following conclusions: 1.The incompatibility of Fe 3+ in the minerals in the crystallization sequence of the terrestrial MO results in its progressive enrichment, and in turn, in the increase of the f O 2 throughout the crystallization.2. Depending on its initial redox state (set by the effective conditions of core-mantle equilibration), the MO can either start from a mildly reduced state (ΔIW = −2 at the surface) and reach a mildly to highly oxidized one (ΔIW = 3 to ΔIW = 8), or start slightly more oxidized (ΔIW = 0) and reach higher values.3. The atmosphere in equilibrium with the MO at the final stages of crystallization is consistently oxidized, dominated by H 2 O and CO 2 .4. The bulk H budget has a first-order influence on the lifetime of the MO, yielding a crystallization duration between 67 kyr and 1.53 Myr over the parameter space investigated here.The C budget and crystallization scenario have a second-order impact.The enhanced fractionation resulting from fractional crystallization yields a more extensive outgassing, as well as a more oxidized final state. 5. Contrarily to previous studies, we find that initially reduced MOs have a longer lifetime than initially oxidized ones, in spite of the stronger greenhouse effect of oxidized species.This is due to the steeper lapse rate in reduced atmosphere, which yields a lower OLR and thus slower cooling.6.The retention of volatiles in melts trapped in the solid cumulates can amount up to 93% of the bulk volatile budget.This mass of volatiles constitutes a reservoir than will possibly undergo further redox speciation influenced by solid-mantle buffers, and can alter the chemistry of the post-MO atmosphere if released at the surface.7. The most volatile-depleted cases are most sensitive to the scaling used for the turbulent heat flux through the MO: adopting the hard turbulence scaling rather than the soft one yields an increase in the MO lifetime of up to 75% for these cases.
While these conclusions have been drawn in the context of the terrestrial MO formed after the Moon-forming giant impact, the processes investigated apply to early evolution of all terrestrial planets.Some parameters, such as initial volatile contents can vary significantly depending on the accretion history, but the trend of increasing oxidation state with ongoing MO crystallization and its associated consequences on MOoutgassed atmosphere is a general feature of rocky planets' and exoplanets' MOs.
The authors thank the two anonymous reviewers for providing constructive comments, which helped improve the quality of this study.M.M. thanks D. Bower for insightful discussions at an early stage of this project, J. Deng for sharing the Python routines used to calculate the pressure-dependence of the f O 2 , and P. Mollière for his help with using petitRADTRANS.
where X i mol is the molar concentrations of oxide i, f O 2 is the oxygen fugacity, R is the gas constant, P 0 = 10 5 Pa is the reference state pressure, and a, b, c, ΔC p , T 0 , and γ 1,K,9 are the fitted parameters given in Table 2 of Hirschmann (2022).The composition in all oxides is taken from the BSE composition of Palme & O'Neill (2003).Notice that for iron oxide, only the initial value of FeO+FeO 1.5 corresponds to Palme & O'Neill (2003).
In order to determine the initial ferric-to-total iron ratio, we evaluate Equation (A1) at the P-T-conditions met at the surface and at the bottom of the magma ocean (MO) in the "reduced" and "oxidized" cases, respectively (see Figure A1), and at ΔIW = −2 at the corresponding P-T conditions, where IW(P, T) is given by Hirschmann (2021).
Subsequently, to calculate the evolution of the f O 2 in the atmosphere, we use Equation (A1) to compute f O 2 at the P-T conditions met at the surface of the magma ocean, with X X given by the differential partitioning of each iron oxide (see Appendix A.3).

A.2. Initial Redox State
The initial f O 2 profile of the oxidized and reduced fractional crystallization cases are represented in Figure A1.The respective values of the initial Fe 3+ /ΣFe are calculated at the effective conditions of core-mantle equilibration (represented by the markers in both panels), to match ΔIW = −2.The difference in f O 2 profiles between the fractional and equilibrium

R e t r a c t e d
crystallization cases sharing the same core-mantle equilibration depth (i.e., between curves of the same color) is due to the colder initial temperature profile in the latter case, which intercepts the liquidus at 55 GPa, while the former intercepts the RCMF.

A.3. Ferric Iron Oxide Partition Coefficients
We calculate the partition coefficient of ferric iron oxide at each crystallization step as the mass-weighted average of the mineral-melt partition coefficient of ferric iron oxide in each mineral present in the composition of the crystallized cumulates layer.We first calculate a pressure profile of ferric iron oxide partition coefficient for the crystallization sequence of Elkins-Tanton (2008; for a 2000 km deep MO, truncated at 55 GPa).In the equilibrium crystallization scenario, the effective partition coefficient of ferric iron oxide is equal to this profile averaged over all pressures between the bottom of the MO and the pressure at which the adiabat intercepts the liquidus.In the fractional crystallization scenario, the partition coefficient of ferric iron oxide is equal to the value on this profile at the pressure of the bottom of the MO.Individual mineral-melt partition coefficients for ferric iron oxide are listed in Table 4.

A.4. Volatiles Redox Equilibria
Equilibirum constants for the volatile redox reactions are calculated form the Gibbs-free energy of reaction (Δ f G 0 ), itself calculated from the Gibbs-free energy of each of the reactants and products (Δ r G X ): where T is the temperature, and R is the gas constant.The values of the species formation Gibbs-free energy are interpolated at the desired temperature from the NIST tables (https://webbook.nist.gov/).

A.5. Volatiles Solubilities
The volatile solubility laws are given in Equation (11).The species-specific α and β coefficients are given in Table 5.

A.6. Mass Conservation
In both crystallization scenarios, the potential temperature entirely defines the advancement in crystallization and thus the mass of the chemical system M syst .Hence, we can express the mass change dM syst as a function of dT pot : dM syst = dM syst /dT pot × dT pot .We obtain dT pot from the thermal evolution model (Equation ( 15)).How dM syst /dT pot is obtained in either crystallization scenario is described below.
In the equilibrium crystallization scenario, the bottom boundary of the chemical system is defined by the intercept between the adiabat and the RCMF temperature; we thus have: , where R RCMF (T pot ) is the radius at which the adiabat (having potential temperature value T pot ) intercepts the RCMF.Taking the derivative of this expression yields: which can be calculated from the shape of the melting curves and the equations of state, by mapping R RCMF onto T pot and taking the derivative of the function thus obtained.
In the fractional crystallization scenario, the decrease in the chemical system's mass is exactly equal to the decrease in melt mass (i.e., the increment in crystals mass) due to crystallization, which is calculated by integrating the melt fraction , and thus: where R MO (T pot ) is the bottom of the MO when the potential temperature is T pot .When evaluating Equation (A5), f is obtained from Equation (1), even though it does not represent the actual profile of the melt fraction in the fractional crystallization scenario, since crystals are assumed to settle (the melt fraction profile is just a step function).Actually, Equation (1) does represent the melt fraction profile in the MO for the vanishingly small (virtual) instant between cooling of the adiabat and settling of the crystals, as represented in Figure 1.

B.1. Secular Cooling
The secular cooling term in Equation ( 16) is calculated as ΔE th /ΔT pot , where ΔE th is the difference between the thermal energy of the initial (hot) and final (cold) states, as represented in Figure B1.Here we present the calculation of ΔE th based on the fractional crystallization case, which differs from ΔE th in the equilibrium crystallization case, but the ratio ΔE th /ΔT pot is assumed to be the same.In the initial state, the whole planet follows a single adiabat whose slope depends on the equation K, corresponding to a 55 GPa-deep MO).In final state, the temperature in the mantle is given by the adiabat from a potential temperature is equal to the surface temperature of the solidus ( = T 1362 f pot K; Fiquet et al. 2010); hence, ΔT pot = 1577 K), and the temperature in the core follows an adiabat from the estimated present-day CMB temperature (4400 K; Andrault et al. 2016).This yields a total thermal energy difference of 1.31 × 10 31 J, and in turn, a secular cooling term of 8.34 × 10 27 J/K.This value, calculated with initial and final states corresponding to the fractional crystallization, applies to the equilibrium crystallization as well, since dE th /dT pot is independent from the potential temperature span covered.

B.2. MO Viscosity
The viscosity of silicate melts depends on T following a Vogel-Fulcher-Tammann law (as in Nikolaou et al. 2019): with A K = 2.4 × 10 −4 Pa s, B K = 4600 K and C K = 1000 K.We do not consider the dependence of η on the water content in the melt.While we find that this content can increase by 1 order of magnitude during MO crystallization, reaching the order of ∼2 wt% (thus possibly decreasing η by about 1 order of magnitude), this effect is small compared to the general increase of η by several orders of magnitude due to temperature decrease as the MO crystallizes (see Figure 2(a) of Nikolaou et al. 2019).Furthermore, no model for the influence of water on silicate melt viscosity has been calibrated on a peridotitic composition, to our knowledge.Nikolaou et al. (2019) used a parameterization calibrated on basanitic composition (with a modified pre-factor), which differs significantly with the BSE composition that we consider in this study.
In the equilibrium crystallization case, the viscosity in the MO is altered by the presence of crystals.The effective viscosity (used to compute the convective heat flux in the MO) depends on f and reads (following Roscoe 1952): where η(T) is given Equation (B1).

B.3. Water Saturation Pressure
In the moist troposphere, the partial pressure of water (p H O 2 ) is equal to its saturation pressure at ambient temperature, given by Pierrehumbert (2010): where p ref and T ref correspond to a reference point on the saturation curve of water, here taken to be the triple point (see Table 3).Equation (B3) is derived from the Clausius-Clapeyron equation with the assumption that the latent heat is constant.

Appendix C Sensitivity of the MO Lifetime
Many factors affect the lifetime of the terrestrial MO (Monteux et al. 2016;Nikolaou et al. 2019).In this study, simplifications were made for the sake of interpretability of the model results.We discuss the anticipated consequences of accounting for the neglected processes in this section.

C.1. Approximations in the Secular Cooling
Our thermal evolution is simple, and the secular cooling term in particular lacks complexity.It aims to provide an upper bound for the MO lifetime, which can then be used as a conservative estimate, e.g., to address the effect of atmospheric escape while the MO is extent.The secular cooling term, as it is defined (Section B.1) and used (Section 2.3.1), has three main sources of inaccuracies or simplifications.First, dE th /dT pot is probably not constant, as intrinsically time-dependent processes will be at play during MO crystallization, which have timescales comparable with the MO lifetime, for instance the onset of solid-state convection in the solidified mantle (Ballmer et al. 2017;Maurice et al. 2017).As a result, separating the time and temperature dependence of the secular cooling term as implied in going from Equations (15) to (16) might be inaccurate.However, resolving these processes is largely beyond the scope of this study, and taking the dE th /dT pot constant is a first-order approach that we expect to yield correct order-of-magnitude results.Second, the initial thermal state of the post-giant-impact Earth is largely unconstrained.This thermal state is the result of the primordial accretion heating of the planet (subsequently dissipated by up to 150 Myr of thermal evolution prior to the Moon-forming impact) as well as internal heating provided by radioactive elements (operating over the same time span), and the heating from the impact itself.Here, the initial thermal state is entirely defined by the choice of a 55 GPa-deep MO as well as a temperature profile following one single adiabat.We make the assumption that such a thermal state is compatible with the anterior thermal evolution of the proto-Earth.Third, the final thermal state, at the end of the MO crystallization, is also the result of hypotheses, namely that the mantle and core follow each a single adiabat.The mantle might instead have a hotter temperature corresponding to the transient regime of solidstate convection onset (although solid-state convection might have time to be fully developed, in which case the temperature profile is likely adiabatic indeed).This assumption is in line with previous studies on terrestrial MO solidification (Lebrun et al. 2013;Nikolaou et al. 2019), easing the comparison between their results and ours.It is more likely that the core follows an adiabat, as inner core growth will not have initiated yet.However, the potential temperature of the core's adiabat (defined at the CMB) is uncertain.We evaluated the secular cooling term for a post-MO core temperature of 4400 as suggested by Andrault et al. (2016).While a lower temperature is not excluded, we remark that the cooling of the core by 1000 K between the initial and final thermal states represents only about 10% of the secular cooling, so an error on the core temperature of a few 100 K is unlikely to dramatically affect our results.

C.2. Atmospheric Evolution
In line with our assumption of total rain-out of the moist troposphere (following Pierrehumbert 2010), we disregard cloud formation.Clouds are arguably the most complicated process influencing heat transfer through the atmosphere, as they can absorb infrared radiation (thus decrease the planet's cooling rate) as well as increase the albedo (thus increase the planet's cooling rate).The former effect can be accounted for in 1D column atmospheric models like the one we use, but requires cloud microphysics parameterizations that add several new parameters (e.g., turbulence diffusion in the clouds and size distribution of droplets).The latter effect requires 3D circulation models that can resolve the cloud distribution.Pluriel et al. (2019) ran such simulations for CO 2 -H 2 O atmospheres and proposed a parameterization of the albedo as a function of the composition of the atmosphere and the surface temperature.However, here we also consider reduced species, which can alter the structure of the troposphere (in particular, the extension of the moist troposphere).Hence, the parameterization of Pluriel et al. (2019) may not be suitable for such atmospheres.

Figure 1 .
Figure 1.Two different crystallization scenarios: (a) equilibrium crystallization where crystals remained entrained by the flow and the bottom of the MO is where the melt fraction reaches the RCMF, and (b) fractional crystallization where crystals settle at the bottom and form a compacted cumulates pile, whose top corresponds to the bottom of the MO.Curves are simply illustrative and do not correspond to the actual melting curves and MO adiabat.
(b) is

Figure 2 .
Figure2.Evolution of Fe 3+ /ΣFe (a) and the surface f O2 (b) throughout MO crystallization.On both panels, solid curves correspond to fractional crystallization and dashed curves to equilibrium crystallization, while orange curves correspond to an initial ΔIW = −2 at the surface (reduced cases) and blue curves to an initial ΔIW = −2 at the bottom of the MO (resulting in a higher surface f O2 , hence oxidized cases).Notice that the fraction of solidified MO (x-axis) is different from the crystal fraction in the case of equilibrium crystallization because the former includes the melt trapped in the solid cumulates while the latter does not.

Figure 3 .
Figure 3. Evolution: (a) of the partial and total pressures and (b) the solubilities throughout MO solidification.In panel (b), thick curves represent the total solubility of an element, while the thin curves represent the species-specific solubilities (corresponding to panel (a)'s legend).The surface ΔIW at the corresponding fraction of solidified MO is reported on the top x-axis.

Figure 4 .
Figure 4. Evolution of the molar mixing ratios for BSE H contents between 1 and 9 EOs-worth of H and C/H ratios from 0.5-2 for an MO undergoing equilibrium crystallization (initial ΔIW = −2 at the surface).

Figure 5 .
Figure 5. Molar fraction evolutions (colored areas, read off left y-axis) and average molecular mass in the atmosphere (black curve, read off right y-axis) for a case with 3 EOs-worth of H in the BSE and C/H = 1, for fractional (upper row) and equilibrium (lower row) crystallization and initially reduced (left column) and oxidized (right column) MO.

Figure 7 .
Figure7.Profiles of temperature (solid lines, read off bottom x-axis) and OLR contributions (dashed lines, read off top x-axis) for atmospheres overlying an initially oxidized (blue lines) and an initially reduced (orange lines) MO, after 0 wt% (panel (a)), 75 wt% (panel (b)), and 98 wt% (panel (c)) of the MO solidified.The f O2 of the atmosphere at each corresponding stage is indicated in the legend.Notice that the y-axes are the same for all three panels, while their corresponding x-axes (both top and bottom) vary from one panel to the other.

Figure 8 .
Figure 8. Fraction of volatiles initially in the MO+atmosphere retained in the final cumulates (H (a) and C (b)), and their corresponding contents in the solid cumulates (H (c) and C (d)), for all equilibrium crystallization cases.We note again that in our model, the volatile retention in solid cumulates is in the form of trapped melts and not via incorporation in the crystal lattice.Colors distinguish C/H while linestyles distinguish initial redox state (solid for reduced, dashed for oxidized).
(a)), which is provided by a larger ΔT (Figure 9(b)).The resulting cooling flux of the MO is lower in the hard turbulence case during the whole MO crystallization (Figure 9(c)).
)-(f) are similar to Figures 9(a)-(c), but for a case having the maximum volatiles budget.While there still is a clear difference in the Ra (Figure 9(d)) and associated ΔT (Figure 9(e)) between hard and soft turbulence scalings, the resulting MO cooling fluxes are much more similar (Figure 9(f)).

Figure 9 .
Figure 9.For the "soft" (blue curves) and "hard" (orange curves) turbulence scalings, evolution of (panels (a) and (d)) the Rayleigh number, (panels (b) and (e)) the temperature drop across the thermal boundary layer, and (panels (c) and (f)) the OLR.Panels (a)-(c) correspond to the most volatile-depleted case, and panels (d)-(f)correspond to the most volatile-rich case of the parameter space.In all panels, the x-axis is the mass fraction fraction of solidified MO (as in Figures2 and 3).

Figure A1 .
Figure A1.Initial profiles of the f O2 (expressed as ΔIW) (a) and the temperature (b) for all distinct cases.The linestyles in both panels and the colors in panel (a) correspond to those used in Figure2: orange for reduced cases and blue for oxidized ones, solid for fractional crystallization and dashed for equilibrium crystallization.The markers indicate the P-T-f O2 conditions of core-mantle equilibrium (ΔIW = −2): orange for reduced cases (equilibration at the surface) and blue for oxidized ones (equilibration at 55 GPa), diamonds for fractional crystallization and circles for equilibrium crystallization (notice that markers of the same color overlap in panel (a)).

Table 1
Notation Used in Section 2.1 Notes.Notice that elements are generally designated by the letter e and gaseous species by the letter s.The average molecular mass in the atmosphere is the sum of the molecular mass of each species weighted by its molar fraction.
a Negative as the MO crystallizes.b Fugacities are equal to partial pressures normalized by a reference pressure of 1 bar.c d Molecular masses are generally given in g mol −1 and atomic masses in atomic mass unit (amu), which convert into each other with a factor unity; here we use SI units and convert both in kg mol −1 .The Planetary Science Journal, 4:31 (22pp), 2023 February Maurice, Dasgupta, & Hassanzadeh

Table 3
Notation Used in Section 2.3

Table 4
Mineral-melt Partition Coefficients for Ferric-iron Oxide

Table 5
Solubility Law Coefficients for the Different Gaseous Species α is expressed in ppm/Pa 1/β .Notice that α values of H 2 , CO, and CO 2 are given for element solubility in Hirschmann (2016), and are here corrected for species solubility.The Planetary Science Journal, 4:31 (22pp), 2023 February Maurice, Dasgupta, & Hassanzadeh R e t r a c t e d of state of the ambient material (liquid MgSiO 3 in the mantle; de Koker et al. 2013; liquid iron in the core; Saxena & Eriksson 2015), given by the initial potential temperature a Figure B1.Initial (blue) and final (orange) thermal states from which the secular cooling term is calculated.20