© The Authors 2025

1 Introduction

After planet formation, rocky planets¹ can start their evolution with a global magma ocean and outgassing from its interior that replaces the primary hydrogen-helium envelope inherited from the protoplanetary disk with a secondary atmosphere (e.g., Lammer et al. 2018; Stüeken et al. 2020). The composition of the secondary atmosphere depends on the redox state of the mantle (Deng et al. 2020; Ortenzi et al. 2020).

It is generally recognized that the magma ocean not only represents a direct link between atmospheric properties and the rocky planet’s mantle, but is further crucial for determining the abundances of important volatiles such as water and carbon dioxide in the mantle and the atmosphere (Chao et al. 2021; Barth et al. 2021; Moore et al. 2023; Krissansen-Totton & Fortney 2022). Therefore, despite its relatively short duration (1–100 Myrs), the magma ocean phase “sets the stage” for the long-term planetary evolution over billion of years (e.g Chao et al. 2021; Lammer et al. 2018; Stüeken et al. 2020).

For Earth-sized planets like Venus, Earth, and the potentially habitable exoplanets TRAPPIST-1 e, f, and g, predominantly H₂O-CO₂ outgassing is expected (Deng et al. 2020). However, a water-vapor-dominated atmosphere is subject to extreme ultra violet irradiation (XUV) and thus H₂O photolysis with subsequent escape of H₂ can occur. Therefore, water-dominated atmospheres are prone to significant mass loss, even in the Solar System with its relatively calm host star. Hamano et al. (2013) demonstrated that atmospheric erosion in the inner Solar System could deplete all of Earth’s oceans within 100 Myrs. Venus experienced a relatively long magma ocean stage of 100 Myrs, during which most of its water was lost, whereas the shorter magma ocean stage on Earth, lasting just a few million years, prevented significant water loss. For Earth-sized planets, the magma ocean stage thus plays a particularly critical role in determining their surface water content, and thus their habitability (e.g., Barth et al. 2021; Hamano et al. 2013; Tian & Ida 2015). The interaction between H₂O outgassing and atmospheric escape during the magma ocean stage is even more crucial for rocky planets orbiting active M dwarf stars, which emit intense XUV flux during a prolonged pre-main-sequence stage (e.g., Johnstone et al. 2021; Tian & Ida 2015).

It was shown by Raymond et al. (2022) that the TRAPPIST- 1 planets – in contrast to Earth – were not strongly modified by large impactors during late accretion, which would alter the composition of the outgassed atmosphere after the magma ocean stage (e.g., Zahnle et al. 2020, 2015). Thus, any CO₂ and H₂O observed in the atmospheres today must be primordial (see e.g., Krissansen-Totton & Fortney 2022). Hence, atmospheric constraints on the rocky planets in the TRAPPIST-1 system may shed light on the volatile budget of the magma ocean, planet formation, and also the potential for habitability.

Initial observations with JWST suggest that TRAPPIST-1 b (Greene et al. 2023) and c (Zieba et al. 2023) do not possess dense CO₂ atmospheres and even indicate the absence of a substantial atmosphere with p_surf ≥ 0.1 bar. However, recent analyses suggest that the MIRI data for TRAPPIST-1 c are consistent with a variety of atmosphere compositions: water vapor (Acuña et al. 2023), O₂ with CO₂ or H₂O (Lincowski et al. 2023) with surface pressures up to 40 bar. These results are, however, based on only five secondary eclipses for TRAPPIST-1 b and four for TRAPPIST-1 c, which may also be affected by stellar contamination (Lim et al. 2023).

The combination of ultraprecise density constraints (Agol et al. 2021) and interior models (Noack et al. 2016; Unterborn et al. 2018; Dorn et al. 2018) suggest that at least TRAPPIST- 1 g currently incorporates >10 wt% of water in its interior structure, which would correspond in mass to more than 100 terrestrial oceans (TOs) of water² (Unterborn et al. 2018; Barth et al. 2021; Raymond et al. 2022). The high inferred volatile content of TRAPPIST-1 g is supported by planet formation models (Miguel et al. 2020; Unterborn et al. 2018; Schoonenberg et al. 2019) that suggest a volatile-rich formation scenario for the outer TRAPPIST-1 planets with initial water mass fractions of up to 50 wt%. Determining the magma ocean solidification time for rocky planets and the resulting distribution of volatiles are thus pressing science questions of current and future missions like PLATO (Turbet et al. 2019; Schlecker et al. 2024) and mission concepts like the Large Interferometer for Exoplanets (LIFE, Bonati et al. 2019).

There are, however, still several open questions in the modeling of magma oceans on diverse rocky planets. Many previous magma ocean models, including the open-source MagmOc1.0 (Barth et al. 2021) as part of the VPLANET modeling suite (Barnes et al. 2020), assumed a pure H₂O steam atmosphere (Hamano et al. 2013; Goldblatt et al. 2013; Lichtenberg et al. 2021). Water is not the only critical molecule in terrestrial planet evolution Elkins-Tanton (2008); Nikolaou et al. (2019); Krissansen-Totton & Fortney (2022). Mixed H₂O-CO₂ outgassing can also modify the volatile distribution as the atmosphere composition changes from an initial CO₂ -dominated atmosphere to a H₂O-dominated atmosphere one (Bower et al. 2019). The impact of simultaneous H₂O and CO₂ outgassing during the magma ocean stage on the potentially habitable TRAPPIST-1 planets e, f, and g together with atmospheric escape has not yet been investigated.

We tackle for the first time the H₂O-CO₂ outgassing feedback between H₂O and CO₂ outgassing on TRAPPIST-1 e, f, and g with an upgraded version of the Barth et al. (2021) model, VPLanet/MagmOc2.0. In Sect. 2.2, we first introduce a new solution for coupled outgassing of two volatiles. We next present an update of the atmospheric evolution model, in which we take into account the vertical extent of a mixed CO₂-H₂O atmosphere (Sect. 2.4).

We employed full radiative transfer (RT) calculations with petitRADTRANS (Mollière et al. 2019) for two different atmosphere treatments: in the RT grid model, we compiled the outgoing longwave radiation (OLR) in a 3D emission grid for different surface pressures, surface temperatures, and water volume mixing ratios and interpolated during the simulation time (Sect. 2.5). In the corrected gray atmosphere model, we used the RT grid calculation to formulate an empirical approximation that reproduces the results of the RT grid atmosphere model generally within 10% accuracy (Appendix B).

We also investigated how the runaway greenhouse limit changes with different surface gravities (7.5–22.5 m/s²; Sect. B.4). After validation of various Earth scenarios and outgassing laws (Appendix C), we applied MagmOc2.0 to the potentially habitable TRAPPIST-1 planets e, f, and g, in which we investigated initial volatile contents between 1 and 100 TO H₂O (Sect. 3). For each initial water scenario, we calculated evolutionary tracks with no CO₂ and initial CO₂ mass equal to 0.3× and 1× initial water mass, respectively, to illustrate the impact of additional CO₂ on the thermal and volatile evolution for a relatively coarse H₂O grid. For each evolution scenario, we further implemented two albedo (α) assumptions: α = 0.75 for a better comparison with Barth et al. (2021), and α = 0, the clear-sky assumption.

We first present for TRAPPIST-1 g and the clear-sky scenario three example evolutions with an initial water mass of 1, 5, and 100 TO (Sect. 3.1) as in Barth et al. (2021). Subsequently, a concise overview of solidification times and remaining water fraction in the solid mantle after the magma ocean evolution is given for TRAPPIST-1 e, f, and g, for an initial water mass between 1 and 100 TO and for the three CO₂ ratios (Sect. 3.2). Finer evolution grids with respect to initial water and CO₂ masses for TRAPPIST-1 e, f, and g are found in a detailed grid of end states of H₂O, CO₂, and O₂ partial pressures as well as sequestered H₂O and O₂ in the mantle. We further investigate the strong impact of CO₂ on the desiccation, and thus the magma ocean lifetime for water-poor composition (≤10 TO H₂O).

We tie the explored magma ocean simulations to new constraints on the interiors of all TRAPPIST-1 planets. Here, the inner planets (b, c, and d) yield a constraint on the iron fraction of 27 wt-% that allows us to place TRAPPIST-1 e in the dry water regime (≤10 TO H₂O), for which a partly desiccated CO₂- atmosphere is expected. The amount of abiotically produced O₂ in these scenarios depends strongly on the magma ocean lifetime. TRAPPIST-1 f and g are expected to be in the water-rich regime (>>10 TO H₂O), based on the interior models, for which the magma ocean ends in a “wet” CO₂ atmosphere.

Next, we explore the impact of an extended H₂O and CO₂ atmosphere on the planetary radii and measured bulk density. We also present a compositional model of the refractory elements and volatiles present during planet formation adjusted for the metallicity of TRAPPIST-1, and compare it to interior structures with respect to the iron fractions and volatile ratios of TRAPPIST-1 e, f, g that are compatible with the measured masses and radii of Agol et al. (2021) (Sect. 3.3). The results of the magma ocean, composition, and interior models are presented in Sect. 4 and discussed in Sect. 5. We conclude in Sect. 6 that comparative studies between TRAPPIST-1 e and g and the inner planets are warranted to identify end states of volatile-poor formation (inner planets), desiccated evolution (TRAPPIST-1 e), and volatile-rich formation with little modification by desiccation (TRAPPIST-1 g). We outline future avenues to improve MagmOc2.0 for a better link to planet formation and also to tackle the full diversity of outgassed atmospheres during the crucial magma ocean stage (Sect. 7).

2 Methods

We utilized and expanded the open source VPLANET framework (Barnes et al. 2020) that connects stellar and planetary processes in order to simulate the evolution of stars and planets over time spans of gigayears. For the present paper, we used the stellar module, which tracks the bolometric luminosity of the star according to the Baraffe et al. (2015) stellar evolution model grid, and the XUV evolution according to the Ribas et al. (2005) model originally developed for solar twins, but which appears compatible with lower-mass stars (see e.g., Richey-Yowell et al. 2022). The VPLANET AtmEsc module tracks water photolysis by the stellar radiation, hydrogen escape via energy- and diffusion-limited escape, and oxygen escape via hydrodynamic drag (Watson et al. 1981; Hunten et al. 1987; Luger & Barnes 2015). The VPLANET EqTide module simulates tidal effects, including frictional heating for the both the constant-phase-lag and constant-time-lag models (Ferraz-Mello et al. 2008; Leconte et al. 2010; Barnes et al. 2013). The RadHeat module tracks the radiogenic heating of the unstable isotopes ⁴⁰K, ²³²Th, ²³⁵U, and ²³⁸U.

2.1 MagmOc approach

MagmOc is a geophysical and geochemical model for the coupled mantle and atmosphere when the mantle is at least partially molten (Barth et al. 2021). We assume a bulk silicate Earth composition from O’Neill & Palme (1998), following Schaefer et al. (2016). Table 1 provides the most relevant geophysical parameters for this work.

Our magma ocean model assumes as outlined in Barth et al. (2021) efficient cooling via convection of the magma ocean, which is only true until the melt fraction at the surface drops below 0.4. This condition is reached in our model typically with surface temperatures of 1650 K and a solidification radius r_s that already comprises about 99% of the planetary radius. Below this temperature, we follow the argumentation of Debaille et al. (2009) that a thick thermal boundary layer may be neglected toward the end of the magma ocean because crystallization of iron-rich minerals lead to overturning near the surface and a resetting of the thermal boundary layer. The module thus switches to a solid-like viscosity (Barth et al. 2021, Eq. (3)) and advances solidification and thermal evolution further until rs is equal the planetary radius, which is typically for surface temperatures of 1400 K. At this point, our model does not further advance mantle and surface temperature evolution in contrast to other models (Lebrun et al. 2013; Bower et al. 2019, 2022; Krissansen-Totton & Fortney 2022; Lichtenberg et al. 2021).

We note that it is worthwhile to extend simulations beyond the mantle solidification time with MagmOcV2.0 for cases with significant atmospheric erosion during the magma ocean stage to capture the final stages of complete desiccation via atmospheric escape that is an upper atmosphere process. The caveat is here that a) surface temperatures and heat fluxes are kept “artificially” high after solidification and b) no water condensation can occur, which would “save” the remaining water from being lost. Thus, water loss due to atmospheric erosion may be overestimated.

2.2 H2O-CO2 volatile budget

Here, we focus on advancing the volatile treatment of MagmOcV1.0, which is needed for coherent modeling of complete devolatilization on rocky exoplanets with a large variety of initial compositions and H₂O-CO₂ ratios during the magma ocean stage. Table 2 shows an overview of the relevant parameters.

The volatile budget was set up similarly to Barth et al. (2021), where we assume at all times that a volatile i (i=H₂O,CO₂) is distributed in the fully coupled magma ocean-atmosphere (MOA) system. This system consists of a liquid magma ocean, M^liq , the atmosphere, M^atm, and crystals forming within the liquid magma ocean, M^crystal ; the latter was set to zero initially.

We further assumed initially that the volatiles, i, in the liquid magma ocean and atmosphere are in balance with each other for a given initial volatile mass, $M_i^{initial}$, which results in an initial volatile mass fraction, $F_i^{initial}$, governing both the amount of volatile outgassed to build an atmosphere and the amount of volatile dissolved in the melt. The mantle begins to solidify from bottom to top; that is, we assumed that the crystals fall to the bottom of the magma ocean, M^crystal, incorporating a small amount of the available volatile mass from the melt (F_iM^liq). Sequestering of volatiles in the solidifying mantle, dM^{S ol}, is regulated by the mantle-averaged constant partition coefficient, k_i. The volatile mass in the solidified mantle $\left(M_i^{Sol}\right)$ is inaccessible for the MOA, and thus comprises a sink term. As the magma ocean depth, and thus the total mass of M^liq , decreases with increasing solidification radius, $\frac{\text{d}{r}_s}{\text{\hspace{0.17em}d}t}$ , the mass fraction, F_i , of the volatile in the melt typically increases, leading to more outgassing.

In summary, we assumed a mass balance in the MOA for a volatile, i, of the form (see also Barth et al. 2021; Bower et al. 2019): $$\begin{array}{l}{M}_i^{\text{moa}}=M_i^{\text{crystal\hspace{0.17em}}}+M_i^{\text{liq}}+M_i^{\text{atm}}\\ =k_i{F}_i{M}^{\text{crystal\hspace{0.17em}}}+F_i{M}^{\text{liq}}+\frac{4\pi r_{\text{p}}^2}{ɡ}{p}_{i,\text{\hspace{0.17em}mass\hspace{0.17em}}},\end{array}$$(1)

where F_i is the volatile mass fraction of the liquid part of the magma ocean, $M_i^{liq}$. In addition, we assumed that the volatile, i, is partly sequestered in the crystallized magma ocean, $M_i^{cry}$. The mass of the volatile in the atmosphere was also set by F_i, where the partial pressure of volatile, i, was determined by a Henrian fit to laboratory data (see Fig. A.1 and references in Appendix A): $$p_{i,\text{part}}={\left(\frac{F_i-c_i}{a_i}\right)}^{b_i},$$(2)

where a_i, b_i, and c_i are the fit coefficients. We note that c_i is a term that suppresses outgassing and is only used in the solubility law of Elkins-Tanton (2008). For the TRAPPIST-1 planet simulations, we used the solubility law of Nikolaou et al. (2019), for which c_i was set to zero. Since we implemented and tested both solubility laws (Appendix A), we show here the generalized form of the Henrian fit with the suppression term.

As is pointed out by Bower et al. (2019), the mass weighted pressure of volatile, i – that is, p_i,mass – enters into the mass balance. The relation between the partial and mass weighted pressure is: $$p_{i,\text{mass}}=\frac{{\mu }_i}{{\overline{\mu }}_{\text{atm}}}{p}_{i,\text{part}},$$(3)

where µ_i is the molecular mass of volatile i and ${\overline{\mu }}_{\text{atm}}$ is the mean molecular mass of all atmosphere constituents.

In addition, atmospheric erosion can remove H₂O from the system, and thus acts as another sink term to the volatile budget, $M_i^{\text{moa}}$ (Fig. 1; see Sect. 2.3 for a more detailed description of atmospheric escape). The two processes that remove volatiles from $M_i^{moa}$ are described by the following differential equations: $$\frac{\text{d}{M}_i^{sol}}{\text{d}t}=k_i{F}_i\underset{\text{d}{M}^{\text{cy\hspace{0.17em}stal\hspace{0.17em}}/\text{d}t}}{\underbrace{4\pi {\rho }_m{r}_s^2\frac{\text{\hspace{0.17em}d}{r}_s}{\text{\hspace{0.17em}d}t}}}$$(4) $$\frac{\text{d}{M}_{{\text{H}}_2\text{O}}^{esc}}{\text{\hspace{0.17em}d}t}=4\pi r_p^2{\varphi }_1\frac{{\mu }_{{\text{H}}_2\text{O}}}{{\mu }_H}$$(5) $$\frac{\text{d}{M}_{{\text{CO}}_2}^{esc}}{\text{\hspace{0.17em}d}t}=0.$$(6)

We refer to Tables 1 and 2 for a concise description of all relevant parameters.

2.2.1 Implementation of outgassing during run time

We used the starting conditions for the magma melt fraction of volatile i – that is, F_i(t = 0) (Fig. 1 left) - as input parameters for MagmOc2.0. We used a root-finding algorithm to solve the mass balance equations (Eq. (1)) with respect to F_i for (t = 0) for a given volatile initial mass, $M_i^{initial}$, to find the starting conditions for a given initial magma ocean depth.

In MagmOc1.0, root finding was performed at all time steps to solve for F_i (t) that satisfies the mass balance, which is still relatively efficient for a single outgassed volatile. In MagmOc2.0, we opted instead to operate on the time derivatives of the mass balance equation, $\frac{\text{d}{M}_i^{moa}}{\text{d}t}$ , which is more numerically efficient for multiple volatiles. We also took into account outgassing feedback that occurs due to changes in the mean molecular mass of the atmosphere, $\frac{\text{d}{\overline{\mu }}_{\text{atm}}}{\text{\hspace{0.17em}d}t}$, when a magma ocean evolves from a CO₂ toward an H₂O-dominated atmosphere (Bower et al. 2019).

During run-time, the derivatives of the volatile mass fractions are advanced using a set of coupled ordinary differential equation: $$\frac{\text{d}{F}_{{\text{H}}_2\text{O}}}{\text{\hspace{0.17em}d}t}=\frac{C_{{\text{H}}_2\text{O}}{A}_{{\text{CO}}_2}-C_{{\text{CO}}_2}{A}_{{\text{H}}_2\text{O}}}{A_{{\text{H}}_2\text{O}}{B}_{{\text{CO}}_2}-C_{{\text{CO}}_2}{B}_{{\text{H}}_2\text{O}}},$$(7) $$\frac{\text{\hspace{0.17em}d}{F}_{{\text{CO}}_2}}{\text{\hspace{0.17em}d}t}=\frac{C_{{\text{CO}}_2}{B}_{{\text{H}}_2\text{O}}-C_{{\text{H}}_2\text{O}}{B}_{{\text{CO}}_2}}{A_{{\text{H}}_2\text{O}}{B}_{{\text{CO}}_2}-A_{{\text{CO}}_2}{B}_{{\text{H}}_2\text{O}}},$$(8)

where A_i, B_i, and C_i comprise different components of the complete derivative of the (volatile) mass balance equations. The full derivation is outlined in Appendix A. The sink terms for the solidified mantle (Eq. (4)) and atmospheric erosion (Eq. (5)) are substracted from $\frac{\text{d}{M}_i^{moa}}{\text{d}t}$.

We note that these coupled differential equations for outgassing of multiple volatiles differ from those of Bower et al. (2019), because we choose the volatile melt fractions F_i and not their partial pressures p_i,part(F_i) as primary variables of integration. We have further verified that MagmOc2.0 with a pure H₂O atmosphere yields the same results as MagmOc1.0 (not shown). Furthermore, the benchmarking that we performed for different Earth scenarios with and without mixed atmospheres further confirms the validation of the code (Appendix C).

2.3 Atmospheric escape

Atmospheric escape was calculated by the VPLanet module AtmEsc, which includes XUV photolysis of H₂O into H and O. We assume that the hydrogen produced from photolysis escapes into space based on the hydrodynamic escape mechanism described in Barnes et al. (2020) (see also Watson et al. 1981, Zahnle & Kasting 1986, and Luger & Barnes 2015). Hydrogen may escape in one of two ways depending on the composition of the atmosphere. In water-dominated atmospheres, hydrogen atoms are liberated from water molecules where the optical depth is approximately unity, which is close to the exobase. Hence, these atoms can escape if the incident photons carry enough energy for the hydrogen atoms to achieve escape velocity, which is generally true for XUV photons. This escape mechanisms assumes as the limiting factor the amount of photon energy deposited in the atmosphere. In this energy-limited regime, hydrogen atoms can also carry away oxygen produced by photodissociation if the hydrogen escapes at sufficiently high velocity (Hunten et al. 1987). Our energy-limited escape model is identical to the one described in Barnes et al. (2020), Appendix A.

If carbon dioxide and/or oxygen accumulates in the atmosphere, however, there may be fewer water molecules than photolyzing photons near the exobase. In this case, the limiting factor is the availability of water molecules, which must diffuse through the background gas(es) to reach the photolyzing layer of the atmosphere. In our model, we assume escape transitions from the energy-limited regime to this diffusion-limited regime when the atmosphere consists of more than 60% carbon dioxide or oxygen. We adopt the diffusion-limited escape model from Luger & Barnes (2015), which is a hybrid of the models originally presented in Zahnle & Kasting (1986) and Hunten et al. (1987). In our model, the flux is given by $$F_{\text{diff}}=\frac{\left(m_{\text{bg}}-m_{\text{H}}\right)\left(1-X_{\text{bg}}\right)b_{\text{bg}}ɡm_{\text{H}}}{k{T}_{\text{flow\hspace{0.17em}}}},$$(9)

where m_b𝑔 is the mass of the background gas (either oxygen or carbon dioxide), mH is the mass of a hydrogen atom, X_b𝑔 is the molar mixing ratio of the background gas a the base of the flow, b_b𝑔 is the binary diffusion coefficient for the dominant background gas, k is the Boltzmann constant, and T_flow is the temperature at the base of the flow. Mason & Marrero (1970) provides the values for the diffusion constants as $$b_{\text{O}}=4.8\times {10}^{19}{T}_{\text{flow\hspace{0.17em}}}^{0.75},$$(10)

and $$b_{{\text{CO}}_2}=8.4\times {10}^{17}{T}_{\text{flow\hspace{0.17em}}}^{0.6}.$$(11)

We assume that T_flow = 400 K for all cases (Luger & Barnes 2015). The only difference between this diffusion-limited escape model and the one described in Appendix A of Barnes et al. (2020) is that here we include CO₂.

2.4 Vertically extended mixed atmospheres of magma oceans

Another improvement of MagmOc2.0 is its treatment of the atmosphere as vertically extended rather than as a single layer. Following the methodology established in previous studies (Goldblatt et al. 2013; Lichtenberg et al. 2021; Schaefer et al. 2016), we constructed pressure-temperature profiles for a given surface pressure, p_{sur f} , surface temperature, T_{sur f} , and atmosphere composition.

For atmospheres consisting of a single volatile, we assumed a dry adiabat for the lower atmospheric layers. For the upper atmospheric layers, we incorporated latent heat release via condensation, described by $$\frac{\text{d\hspace{0.17em}}\ln T}{\text{\hspace{0.17em}d\hspace{0.17em}}\ln p}=R/c_{p,i}(T)\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}dry}$$(12) $$=RT/L_i\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}condensation,\hspace{0.17em}}$$(13)

where p is the gas pressure, T is the gas temperature, R is the ideal gas constant, and c_p,i(T) is the specific heat capacity of the volatile, i (either H₂O or CO₂), using the Shomate equation with values obtained from the NIST data base (Cox et al. 1984; Chase 1998). Values for latent heat, L_i, were taken from Lichtenberg et al. (2021) and Pierrehumbert (2010).

A condensate becomes thermally stable when the p-T profile intersects the relevant condensation curves. For the H₂O condensation curve, we adopted the August-Roche-Magnus formulation described in Alduchov & Eskridge (1996). For pure CO₂, we used the Clausius-Clayperon relation. We validated this method by comparing pressure-temperature profiles for H₂O and CO₂ and T_surf = 500-2000 K with previous work (Goldblatt et al. 2013; Lichtenberg et al. 2021).

For mixed H₂O and CO₂ atmospheres, we adopted the multispecies adiabat formulation with condensation of Graham et al. (2021), where we assumed for simplicity again that CO₂ and H₂O are not removed from an atmospheric layer when condensation occurs (that is, a supersaturation ratio of S = 1). Thus, we assume in a specific atmosphere layer equilibrium between condensation and instant re-vaporization of volatiles. This approach is valid for a hot steam atmospheres with vigorous mixing. We note that, for the calculation of the atmosphere’s pressure-temperature profile, we used, like Lichtenberg et al. (2021), the ideal gas law for the volatiles, which needs to be revisited for very volatile-rich compositions³. However, we first established here a basic framework to expand the boundaries of Solar-System magma ocean simulations by starting with simplified assumptions and offering it as an open-source project to the community. We further note that for the majority of the magma ocean simulations, the radiative properties are set by the water vapor opacities. Water, however, makes the atmosphere in the infrared – that is, for the outgoing long wave irradiation optically thick already for low atmospheric pressures – p > 0.1 bar. Furthermore, we find that the magma ocean solidification occurs, when the planet is in the runaway greenhouse limit. Thus, the cooling of the magma ocean is set by the radiative and thermodynamics properties of upper atmosphere layers, the temperature of which is mostly determined by the latent heat release of water vapor that is not impacted by supercritical conditions at the surface. Also Marcq et al. (2017) have verified that deviations from the ideal gas law are of minor importance for calculating the outgoing long wave irradiation of their magma ocean atmospheres. Supercritical surface layers may impact, however, the dry adiabat and thus the very hot initial conditions of the magma ocean stage. Differences in the initial conditions of the magma ocean, however, were found already for Earth to be of minor importance for the magma ocean lifetime and volatile evolution, as the magma ocean spends the majority of the evolution in the runaway greenhouse regime (Appendix C). In all cases, we utilized a fourth-order Runge-Kutta integrator to compute the pressure-temperature (p_gas, T_gas) profile by integrating upward from surface temperature T_surf = 500–4000 K and p_surf = 0.26–26 000 bar⁴.

Figure 2 (left) shows example profiles for $p_{\text{surf\hspace{0.17em}}}=p_{{\text{H}}_2\text{O}}+p_{{\text{CO}}_2}=260$ bar for different water volume mixing ratios $x_{{\text{H}}_2\text{O}}$. We defined the vertically uniform $x_{{\text{H}}_2\text{O}}$ as the ratio of H₂O partial pressure over the total surface pressure: $$x_{{\text{H}}_2\text{O}}=p_{{\text{H}}_2\text{O}}/p_{\text{surf\hspace{0.17em}}},$$(14)

assuming a well-mixed two-component H₂O-CO₂ atmosphere: $$p_{{\text{H}}_2\text{O}}+p_{{\text{CO}}_2}=p_{\text{surf\hspace{0.17em}}}.$$(15)

We stress again that we operate under the assumption of a highly mixed, hot atmosphere, with instant reevaporation of condensates. Thus, we assumed to the first order that almost all water is in the gas phase.

The mixed H₂O-CO₂ pressure-temperature profiles indicate that water vapor releases significantly more latent heat per volume than CO₂, causing the profiles to closely resemble the pure H₂O profile even when 90% of the atmosphere is comprised of CO₂ (Fig. 2 left). The resulting pressure-temperature profiles for a given H₂ O-CO₂ content are used as input for RT calculations with the open source, atmosphere model petitRADTRANS (Mollière et al. 2019) described in the next subsection.

2.5 Radiative transfer in a vertically extended atmosphere

The thermal evolution of the magma ocean planet is determined by the net outgoing emission F_net , which is defined as the difference between OLR and incoming absorbed stellar radiation (ASR) at the top of the atmosphere: $$F_{\text{net\hspace{0.17em}}}=F_{\text{OLR}}-F_{\text{ASR}}.$$(16)

For the ASR flux, we used $$F_{\text{ASR}}=\sigma T_{\text{eff}}^4,$$(17)

where the planet’s effective or black body temperature, T_eff, was calculated as $$T_{\text{eff}}={\left(\frac{1-\alpha }{4\sigma }\frac{L\left(t_{*}\right)}{4\pi a^2}\right)}^{\frac{1}{4}},$$(18)

with the bolometric luminosity of the star, L(t_*) at stellar age, t_*, where t_* = t + t_ini; that is, simulation t plus initial stellar age of 5 Myrs, the semi-major a and the planetary albedo, α, of the planet. In this work, we we assumed two albedos: 0.75 and 0. The albedo of 0.75 is mostly used to facilitate comparison between MagmOc1.0 and MagmOc2.0. However, our assumption of a well mixed atmosphere with immediate reevaporation of condensates is more consistent with a cloud-free atmosphere with low albedo. Moreover, for planets orbiting ultra cool M dwarfs like TRAPPIST-1, very low scattering on top of the atmosphere (albedo= 0-0.1) is predicted (Kopparapu et al. 2013). Thus, we adopted a clear sky albedo (α) of 0 as the new default.

We calculated the OLR with petitRADTRANS (Mollière et al. 2019) using opacities sampled with the correlated k- method for a wavelength resolution of R=1000. We further assume constant H₂O and CO₂ volume mixing ratios throughout the atmosphere together with the pressure-temperature profiles as outlined in the previous section. The OLR is equal to the integrated emission on top of the atmosphere between 0.2 and 35 microns, similar to the wavelength coverage in Goldblatt et al. (2013).

Similar to previous work (Goldblatt et al. 2013; Lichtenberg et al. 2021; Boukrouche et al. 2021), we emphasize the importance of using the continuum opacities for the dominant greenhouse gases, here H₂O and CO₂ are vital ingredients to correctly model the atmospheres of rocky (exo-)planets in the habitable zone. The opacity sources and their references are provided in Table 3. We compute the H₂O and CO₂ k-tables using the HITRAN2020 line list data and broadening coefficients (Gordon et al. 2022), with a constant line-wing cutoff of 25cm⁻¹. These opacities were converted to a format for input into petitRADTRANS as outlined in Chubb et al. (2021). We have verified that for a pure H₂O atmosphere, the OLR yields the canonical runaway greenhouse limit of 282 W/m² for surface temperatures between 500 and 1800 K (Goldblatt et al. 2013). We further note that for CO₂ we use an additional continuum opacity compared to Lichtenberg et al. (2021), which leads to an overall reduction in emission of 20% in a pure CO₂ atmosphere with p_surf = 260 bar compared to Lichtenberg et al. (2021). Already, Marcq et al. (2017) preformed non-gray RT calculations for mixed H₂O-CO₂ atmospheres on magma oceans. They also retrieve the blanketing effect of the runaway greenhouse limit of 280 W/m2 that is indeed the major factor shaping the magma ocean evolution also for mixed H₂O-CO₂ atmospheres in the simulations presented here. In fact, Marcq et al. (2017) also point out that the runaway greenhouse limit can be lowered by tens of W/m² when CO₂ starts to dominate over H₂O (see also Appendices B and B.4).

The emission curves for p_surf = 260 bar and T_surf for different water volume mixing ratios $x_{{\text{H}}_2\text{O}}$ (Fig. 2 right) are very similar to the emission curve of a pure H₂O atmosphere for $x_{{\text{H}}_2\text{O}}\ge {10}^{-2}$. There are, however, two notable differences: a CO₂ absorption band at 4.3 µm appears as soon as CO₂ is added and its amplitude increases with higher CO₂ abundances. In addition, the overall thermal emission decreases with lower $x_{{\text{H}}_2\text{O}}$ because less water per volume is available that can condense out and thus heat the upper atmospheric layers. Consequently, the upper atmosphere layers cool down and emit less flux (see Fig. 2 left).

For $x_{{\text{H}}_2\text{O}}\le {10}^{-3}$, CO₂ emission begins to dominate, causing the overall emission to increase as $x_{{\text{H}}_2\text{O}}$ decreases, eventually exceeding the runaway greenhouse limit. This effect occurs because the hot dry adiabat of a CO₂-dominated atmosphere extends up to p ≤ 50 mbar into the emitting atmosphere layers.

A closer inspection reveals that the runaway greenhouse limit for a mixed H₂O-CO₂ atmosphere with $x_{{\text{H}}_2\text{O}}<1$ is restricted to a lower surface temperature regime compared to an atmosphere composed of 100% H₂O (see e.g., Goldblatt et al. 2013; Lichtenberg et al. 2021) and further becomes smaller than the canonical value of 282 W/m². Further details on this phenomenon are provided in Appendices B and B.4.

During the runtime of MagmOc2.0, we do not compute full RT; instead, we construct for a given planet a thermal emission grid (see Table 4 for the setup). As the magma ocean evolves, the model interpolates surface temperature linearly, while surface pressure and water volume mixing ratio are interpolated logarithmically within the range ${10}^{-6}\ge x_{{\text{H}}_2\text{O}}\le 1$. For $x_{{\text{H}}_2\text{O}}<{10}^{-6}$, thermal emission for $x_{{\text{H}}_2\text{O}}=0$ is used, representing a pure CO₂ atmosphere. This is done because we interpolate in log-space and thus cannot interpolate to log(0). We also found that for $x_{{\text{H}}_2\text{O}}={10}^{-6}$, the integrated thermal emission deviates by less than 10% from the that of a pure CO₂ atmosphere.

To speed-up computation time and as a “sanity check” for the thermal evolution based on full RT calculations, we have also developed an analytic approximation of the thermal emission of a mixed H₂ O-CO₂ steam atmosphere, including modifications to the runaway greenhouse atmosphere limit with increasing CO₂ content. The equations describing the approximation are outlined in the Appendix B.

3 Magma ocean evolution of TRAPPIST-1 e, f, and g

We tested MagmOc2.0 for different Earth scenarios (Appendix C), confirming previous research findings on the substantial H₂O- CO₂ outgassing feedback (Bower et al. 2019). Furthermore, we found that the addition of CO₂ has a minor impact on the solidification timescale compared to a pure H₂O atmosphere when atmospheric escape is negligible. This similarity between mixed H₂O-CO₂ and pure H₂O simulations is primarily due to the dominance of the runway greenhouse radiation limit for the majority of the magma ocean lifetime. In contrast to the Earth’s magma ocean stage (e.g., Hamano et al. 2013; Barth et al. 2021), atmospheric erosion is not negligible for the potentially habitable TRAPPIST-1 planets e, f, and g.

In this work, we systematically revisit evolution trajectories investigated by Barth et al. (2021) with pure H₂O atmospheres for 1–100 TO initial water mass. We focus on simulations with no CO₂, an initial CO₂ mass equal to 0.3× the initial H₂O mass, and an extreme scenario, where the initial CO₂ mass is assumed to be equal to the initial H₂O mass (see Table 5 for the relevant parameters). We further note that we neglect tidal interactions in accordance with the results by Barth et al. (2021).

The wide range in initial volatile composition encompasses the large uncertainties in volatile content and the H₂O-CO₂ ratio acquired during planet formation (e.g., Bitsch et al. 2019). These diverse scenarios allow us to assess to what extent CO₂ impacts the magma ocean lifetime and volatile distribution on magma oceans with oxidized outgassing in close proximity to an M dwarf star.

We adopted the H₂O and CO₂ outgassing laws of Nikolaou et al. (2019), which better represent the current understanding that CO₂ is already outgassed during the initial magma ocean stage (see Sect. 2.2), whereas H₂O only builds up when the mantle begins to fully solidify. Unless otherwise specified, we adopted the RT atmosphere grid model.

3.1 TRAPPIST-1 g evolution

Following Barth et al. (2021), we selected TRAPPIST-1 g simulations for initial water masses of 1, 5, and 100 TO, respectively, to discuss critical stages in the magma ocean evolution for planets in the habitable zone of TRAPPIST-1. Barth et al. (2021) identified roughly three scenarios for the magma ocean evolution of TRAPPIST-1 g that represent significant stages for understanding the impact of atmospheric erosion on an oxidized magma ocean. We thus explore here the same initial water scenarios for better comparison. To assess the impact of CO₂ on the three outlined scenarios (scenario 1–3), we compare the surface temperature and evolution tracks without CO₂ to those with CO₂ (Fig. 4). We note, however, that Barth et al. (2021) assumed an albedo = 0.75. Here, we focus mostly on simulations with a clear-sky albedo of 0, which is more consistent with our clearsky RT as well as with the low scattering efficiency calculated for planets around M dwarfs (Kopparapu et al. 2013).

We first report the general trends for all scenarios: there is a significant impact on the general volatile distribution (Appendix C.7) in simulations with additional CO₂; the magma ocean starts with a CO₂-dominated atmosphere, then H₂O tends to become the most dominant species as the mantle solidifies. The planet can end its magma ocean stage again with a CO₂ dominated atmosphere, when the planet desiccates. The feedback on the volatile distribution between H₂O and CO₂ during outgassing is mostly evident in changes in CO₂ partial pressures. When H₂O starts to become dominant during late stage mantle outgassing, CO₂ partial pressures drop. If the planet desiccates, CO₂ partial pressures rise. These changes generally align with our results for the Earth magma ocean simulations (Appendix C.3), as well as with those derived by Bower et al. (2019).

Bower et al. (2019) also report that the presence of a thick CO₂ atmosphere delays the outgassing of H₂O during the magma ocean evolution for Earth. We cannot reproduce this claim in our simulations. Instead, we find delayed H₂O outgassing only in the volatile-poor “scenario 1” with albedo 0.75, specifically in the extreme case of adding 1 TO CO₂ (not shown). Closer inspection reveals that the delay in H₂O outgassing arises solely from the slower cooling of the magma ocean and is not attributable to the feedback effect. While we do not find significant changes in time of H₂O outgassing, the thermal evolution and volatile distribution on TRAPPIST-1 g can be impacted with CO₂, depending on initial water content.

In the “dry” scenario with 1 TO H₂O, the atmosphere becomes completely desiccated, which leads to rapid mantle solidification within 7 Myrs with no CO₂ and moderate CO₂ content. If an extreme amount of CO₂ is added, then diffusion limited escape delays desiccation such that complete solidification occurs after 17 Myrs. Furthermore, we find that even in the case of moderate CO₂ content, even a slight delay of solidification by 1 Myrs due to diffusion limited escape facilitates moderate abiotic O₂ build-up of a few bar. Without any CO₂, all abiotically created O₂ is deposited in the mantle before complete solidification.

In the “intermediate” scenario with 5 TO, complete desiccation occurs at a much later date compared to the dry scenario; that is, after 70 Myrs simulation time. Moderate CO₂ content leads again to a slight delay in desiccation such that 10 bar abiot- ically created O₂ can accumulate in the atmosphere compared to the CO₂-free simulation for which only 2 bar O₂ remains. With extreme CO₂ content, complete desiccation is delayed even further by CO₂ diffusion limited escape to 120 Myrs. In this simulation, the planet enters the habitable zone of its host star at 76 Myrs with about 1 TO of H₂O in the atmosphere (200 bar). Beyond 200 Myrs simulation time, the model is no longer suitable to follow the further evolution, as is evidenced by the plateauing of the surface temperature evolution (see Sect. 2).

For the “wet” scenario with 100 TO, all simulations enter the habitable zone at 76 Myrs with more than 2000 bars of H₂O in the atmosphere and the mantle solidifies very late after 150 Myrs. Additional CO₂ extends the solidification time by only few % compared to the H₂O-free simulation time. Generally, additional CO₂ has little effect on the final magma ocean state. The final state is characterized by an accumulation of more than 10000 bars of water in the atmosphere as the mantle solidifies. All abiotically created O₂ enters the mantle already during the magma ocean evolution and thus no O₂ build-up occurs even after several 100 Myrs. The lack of final O₂ build-up beyond solidification may be, however, due to the limitation of the model that is strictly valid “only” until the magma ocean solidification state at 150 Myrs is reached.

A closer look on the evolution of the net fluxes for the three scenarios (Fig. 5) reveals the key processes determining the magma ocean lifetime. For the first, dry scenario, the final thermal evolution is primarily determined by the fate of the water, even in the presence of 200 bar CO₂. The mixed H₂O- CO₂ atmosphere enters the runaway greenhouse radiation limit already after 100 000 years, which reduces the outgoing thermal radiation to less than 300 W/m² as long as water is present. Once all water has been eroded, the remaining 60–200 bar thick CO₂ atmosphere with a 1400 K surface temperature yields net outgoing flux of more than 1000 W/m² that quickly leads to mantle solidification (Sect. 2.5, Lichtenberg et al. 2021). All initial differences with thermal evolution due to differences in the CO₂ content are removed, once the simulations enter the runaway greenhouse limit.

For the second, intermediate scenario, both the evolution of the runaway greenhouse radiation limit and atmospheric erosion determine the mantle solidification time. Here, again the runaway greenhouse radiation limit with less than 300 W/m² is reached after 100 000 years. Because there is more water in the system compared to the dry scenario, the radiation limit is maintained for at least 70 Myrs. The net flux continuously drops during that time mainly because the stellar luminosity continuously decreases as the host star evolves to the main sequence. We further note that simulations with additional CO₂ consistently show a lower runaway greenhouse limit (in particular around 1 Myrs evolution time) by 10s of W/m² as long as the atmosphere remains CO₂ dominated (see Appendix B.4). When the majority of H₂O is outgassed during the final stages of mantle solidification (at 70 Myrs), the thermal evolution curves converge for the moderate CO₂ and CO₂-free simulation until complete desiccation occurs that triggers quasi-instant solidification. For extreme CO₂ content, however, diffusion limited escape can prevent desiccation and thus delays solidification to 120 Myrs.

The third, wet scenario with 100 TO H₂O is entirely determined by the runaway greenhouse limit. Additional CO₂ leads only initially to differences in evolution. Once, the simulations enter the runaway greenhouse limit (after a few 100 000 Myrs), these differences cease to matter. We note here again that the simulations with additional CO₂ have a lower runaway greenhouse radiation limit, as long as CO₂ remains the dominant atmosphere constituent. All simulations converge, however, once the majority of the water is outgassed after 100 Myrs. Even with extreme CO₂ content, solidification occurs only slightly later compared to the CO₂-free simulation.

In Barth et al. (2021), the radiation limit was assumed to apply only for surface temperatures cooler than 1800 K, regardless of surface pressure. Within that framework, thick H₂O-dominated atmospheres entered the runaway greenhouse radiation limit very late leading to an extension of the magma ocean lifetime. In this work, however, we find that magma oceans with thick H₂O dominated atmospheres enter the greenhouse limit for hotter surface temperatures (Appendix B). Thus, we have not reproduced the magma ocean lifetime extension proposed in Barth et al. (2021) for very high volatile content.

In summary, the fate of water effectively dominates magma ocean solidification lifetimes for all cases. CO₂ generally only influences the timescale of water erosion. Because the atmosphere keeps more water vapor, in particular for extreme CO₂ content, the magma ocean lifetime of TRAPPIST-1 g can be significantly delayed for the dry and intermediate scenarios with 1 TO and 5 TO H₂O, respectively.

We caution, however, that we cannot predict how long the 1400 K surface temperatures can be maintained, because the geophysical model terminates the evolution of the mantle at this point. It thus remains to be confirmed if the extreme jump in outgoing flux as the planet desiccates that triggers immediate solidification can be re-captured with other models that follow the mantle evolution beyond the magma ocean solidification state.

3.2 Overview of magma ocean evolution in TRAPPIST-1 e, f, and g

We find that additional CO₂ has a particular strong impact on the timescale of atmospheric water loss for simulations with dry to intermediate initial water content (≤10 TO H₂O). Because desiccation is in our model followed by complete mantle solidification, the impact of CO₂ on atmospheric erosion also impacts the magma ocean solidification timescale. In addition, the magma ocean duration changes with different scattering (albedo) assumptions, which may also impact atmospheric water loss. Thus, we first inspect the change in atmospheric desiccation time for TRAPPIST-1 e, f, and g and for our nominal clearsky (albedo=0) and the high albedo simulations (albedo =0.75) (Fig. 6).

In our simulations that assume clear-sky albedo, F_net is low and thus the magma ocean cools and evolves over longer timescales compared to the high albedo simulations. This difference in cooling timescale by itself does not appear, however, to strongly impact atmospheric mass loss rate as is evident by comparing pure H₂O simulations for albedo 0, and 0.75, respectively (Fig. 6, solid lines). The addition of CO₂, however, leads to significant changes in atmospheric desiccation.

If a clear sky is assumed, scenarios with extreme CO₂ content only lead to complete atmospheric desiccation with very dry initial water content (≤2 TO H₂O). For the high albedo simulations, the shorter magma ocean duration compared to the clear-sky simulations appears to limit the impact of CO₂. Desiccation is only significantly delayed for TRAPPIST-1 e, even with extreme CO₂ content.

The combination of delayed atmospheric water mass loss and low magma ocean cooling rate in the clear-sky case also leads to a significant extension of the magma ocean solidification timescale with intermediate initial water content (5-10 TO H₂O)⁵ (Fig. 7). This is true for TRAPPIST-1 e, f, and g with extreme CO₂ content (Fig. 7 dotted lines, bottom left panel) compared to the pure H₂O and high albedo simulations (Fig. 7 solid lines, upper left panel). The magma ocean solidification times are so long that the planets can enter their habitable zone before atmospheric desiccation is complete.

For the volatile-rich scenarios (≥ 10 TO H₂O) and with clear skies (albedo = 0), TRAPPIST-1 e, f, and g enter the habitable zone with partly molten surfaces and a mixed H₂O-CO₂ dominated atmosphere, because atmospheric escape is of lesser importance. But even there, more CO₂ leads to an extension of mantle solidification time in particular for TRAPPIST-1 e. In these scenarios, the addition of CO₂ reduces the runaway greenhouse radiation limit by a few tens of W/m² (see also Marcq et al. 2017). For simulations with albedo=0.75, the extension of the magma ocean lifetime for the volatile-rich scenarios is less evident and is at most 10% even for TRAPPIST-1 e.

The duration of the magma ocean also shapes the amount of remaining water in the solidified mantle due to atmospheric erosion (Fig. 7, right panels). The impact of atmospheric erosion is already evident by comparing the remaining water fraction for the individual planets. TRAPPIST-1 g can retain the most water because it is farther away from its host star, and thus experiences less erosion compared to TRAPPIST-1 f. TRAPPIST-1 f experiences less atmospheric erosion than TRAPPIST-1 e, and thus retains more water than the latter.

The impact of water mass loss over different timescales explains the differences for remaining water in the solidified mantle when comparing high albedo (short magma ocean phase <100 Myrs) with clear sky simulations (long magma ocean phase, up to 350 Myrs). As has already been pointed out by Barth et al. (2021), the volatiles in a magma ocean planet are embedded in a strongly coupled system. If this system is subject to continuous water erosion for a longer time, then also more water can be removed from the planet, because more water is outgassed from the melt to compensate for the mass loss from the atmosphere. Consequently, for the same planet and the same initial water input, there is more remaining water in the melt in the high albedo simulations compared to clear sky simulations with longer magma ocean lifetimes.

However, only 2–6%⁶ of the initial water can be retained in the melt in the pure H₂O outgassing scenario for albedo 0.75 and 1–5% for albedo 0. Apparently, the impact of extended water erosion over timescales of 100 Myrs and longer does not lead to a drastic reduction in remaining water. This relatively minor impact on the water mantle content with extended magma ocean lifetime indicates that most water is sequestered in the solidified part of the mantle during the very first million years of simulation time, before atmospheric erosion has had time to remove significant amounts of water from the atmosphere and thus from the remaining molten part of the mantle. It takes typically more than a few million years to significantly erode the atmosphere, except for the very dry scenarios with 1 TO H₂O.

Adding CO₂ tends to increase the percentage of remaining water fraction by up to 2%. For TRAPPIST-1 e, the total amount of water stored in the solid mantle can thus be doubled with extreme initial CO₂ content, assuming a clear-sky albedo. The impact of CO₂ on the remaining water is so strong, because more CO₂ in the system leads to higher water melt fractions for the whole duration of the magma ocean evolution. Consequently, more H₂O can be sequestered in the solid part of the mantle throughout the evolution.

In any case, our work shows that it is important to capture the feedback of CO₂ on the atmospheric water loss and the magma solidification timescale. Large amounts of CO₂ (>1000 bar) can prevent or at least delay complete desiccation for the TRAPPIST- 1 e, f and g planets. We further find that among all investigated planets, TRAPPIST-1 e, orbiting closest to its host star, is more strongly impacted by various processes that shape the magma ocean stage: These are the assumed albedo, which represents here cases with low scattering (albedo=0) versus high scattering of incoming stellar light (albedo =0.75), the reduction in the runaway greenhouse radiation limit with high CO₂ abundances, as well as the reduction of atmospheric water loss because atmospheric escape of H₂O is diffusion limited with high amounts of CO₂ in the atmosphere.

3.3 Atmosphere extent and interior modeling constraints

To be able to understand the impact of an extended atmosphere on the observed radii of the TRAPPIST-1 planets, we first made use of an interior structure model (Noack et al. 2017) employing look-up tables created with Perple_x (Connolly 2009) for thermodynamic properties of the silicate mantles. For a first- order estimate on the potential chemical composition of the TRAPPIST-1 planets, we use an adapted version of the condensation model presented in Bitsch & Battistini (2020), and employed the measured metallicity of TRAPPIST-1 to derive the likely chemical composition of the star and planets (depending on their distance of the star; see Appendix E for more information on the numerical models). For the silicate mantle of the TRAPPIST-1 planets, and under the assumption of an Earthlike mantle iron number of 0.1 (i.e., a magnesium number of 0.9), our model predicts the following molar composition: 5.35% FeO, 48.15% MgO, 39.09% SiO₂, 2.85% CaO, 1.83% Al₂O₃, and 1.0% Na₂O. Residual iron as well as any condensed FeS contributes to the metal core, with a predicted mass fraction of 27%. However, in our interior structure model, we vary the coremass fraction further to obtain realistic interior structure profiles for the TRAPPIST-1 planets that can match their observed radii. For simplicity, in our model we consider only pure iron (see Appendix E). For the innermost planets, b, c and d, a coremass fraction of 27% leads to a model radius of 1.114 Earth radii (compared to an observed radius of ${1.116}_{-0.014}^{+0.012}$, Ago et al. 2021), 1.098 Earth radii (compared to ${1.097}_{-0.014}^{+0.012}$) and 0.774 Earth radii (compared to ${0.788}_{-0.011}^{+0.01}$), and therefore our model radii match the measured radii within the observational error. This suggests that our compositional model is able to correctly predict the planetary composition of the TRAPPIST-1 planets. For the outer planets of the system, the appearance of volatiles adds a degeneracy to our interior structure.

For the building blocks of the outer planets TRAPPIST-1 e, f, and g (and h with almost exactly the same composition as g), our model predicts that the fraction of silicates and metals in the planetary building blocks would decrease to 64 wt-% (e, f) and 57 wt-% (g, h) (see Fig. 8). For e and f, the remainder of the mass is water in the form of ice layers (including high-pressure ice) as well as liquid water. For g and h, volatiles are separated into 32 wt-% H₂O as well as 11 wt-% NH₄. However, during planet accretion, due to high-energy impacts and melting processes in the interior, we expect the final planetary compositions to be drier. In addition, planetary orbits may change during the planet formation process, leading to an additional uncertainty in accreted planetary building material, and hence the fraction of volatile materials.

Table 6 lists the range of iron fractions and water/ice fractions that would match the observed masses and radii of TRAPPIST-1 e, f, and g while taking into account the star- derived mantle chemical setup. Figure F.1 shows the possible H₂O ranges for each planet within error bars. While for all three planets, the measured radii could be explained by a dry composition (with 25.2, 20.4 and 15.9 wt-%, respectively), increasing iron fractions are possible when adding a water/ice layer of increasing extend with increasing distance to the star. It should be noted that here we only considered water ice and did not take into account any contribution of NH₄ , which may be abundant in addition to H₂O in TRAPPIST-1 g following our compositional model. Given the apparent lower iron fraction in the system compared to the solar system (as suggested also in Unterborn et al. 2018), we only investigate core-mass fractions of up to 50 wt-%. Within this range, the maximum water fraction is 17.1 wt-% and therefore below 20 wt-%, in accordance with Dorn et al. (2018), Unterborn et al. (2018), Agol et al. (2021), Barth et al. (2021) and Acuña et al. (2021), even though here we apply our TRAPPIST-1 adapted compositional model instead of an Earth-like mineralogy.

In Fig. 9, we display how the planetary radius and density of the planet underneath the atmosphere would change when considering extended atmospheres of different compositions and pressures, calculated as in Ortenzi et al. (2020). We first calculate the predicted atmospheric thickness for the measured mass and radii of each planet and calculate the related radius change. To give examples on how the density of the planets would be affected, we subtract the predicted radius change from the planet radius and calculated the increased density of the sub- atmospheric planet. We note that for today’s measured planet masses and radii, the extended steam atmospheres depicted here are not realistic and the density variations should therefore only be seen as exemplary values.

The atmospheric pressure has a logarithmic effect on the thickness of the atmosphere, leading to increasing atmospheric pressures showing only a weak effect for pressures above a few tens of bars. Water steam atmospheres have approximately twice as strong an influence on the atmosphere thickness as the heavier CO₂ atmospheres. In addition, the atmospheric temperature strongly impacts the atmospheric thickness, with an almost linear effect (a doubled absolute temperature leads to an approximately doubled atmosphere thickness).

4 Results summary

In this study, we introduce MagmOc2.0, a magma ocean model with versatile multi-volatile outgassing and thermal cooling informed by RT calculations in a vertically extended atmosphere. Comparison of different atmosphere models for the Earth with mixed H₂O-CO₂ content shows that using RT calculation in a vertically extended atmosphere generally yields longer solidification times compared to a gray atmosphere model. The solidification times still remain within one order of magnitude when atmospheric escape is neglected (Appendices C.1 and C.3).

For a pure CO₂ atmosphere, however, the magma ocean reaches the solidification surface temperature of 1400 K significantly faster (by more than one order of magnitude) compared to a gray atmosphere model (Appendix C.2). We find that a vertically extended CO₂ dominated atmosphere with surface temperatures larger than 2000 K emits more thermal flux on top of the atmosphere than is captured with the original gray model.

Consequently, adding CO₂ to a water steam atmosphere leads to faster cooling of the magma ocean during the initial magma ocean stage.

As soon as the runaway greenhouse limit is reached (T_surf ≤ 1800 K), however, a second thermal evolution stage begins, lasting much longer than the initial stage. Here, the addition of CO₂ leads to slower cooling because its presence cools the upper atmosphere and this reduces overall thermal emission compared to a pure H₂O atmosphere.

Application to the potentially habitable TRAPPIST-1 planets e, f, and g shows that the albedo is of major importance in determining magma ocean lifetimes. If an albedo of 0.75 is assumed, as for the Earth-like simulations (Appendix C), than the magma ocean lifetime is always shorter than 100 Myrs for TRAPPIST-1 e, f, and g. An albedo of zero, which is probably more appropriate for hot H₂O-CO₂ dominated atmospheres around M dwarfs (Kopparapu et al. 2013) extends the magma ocean lifetime for all investigated TRAPPIST-1 planets such that they can even enter their respective habitable zone with intermediate to water-rich compositions (T > 10 TO initial water mass). The presence of CO₂, on the other hand, results for water-rich scenarios overall only in a small (1–10%) net extension of the magma ocean lifetime. For dry compositions (≤10 TO initial H₂O) and large amounts of CO₂, CO₂ diffusion limited escape can also significantly extend the magma ocean lifetime for TRAPPIST-1 e (albedo 0.75) and TRAPPIST-1 e, f, and g (albedo 0).

The presence of CO₂, however, modifies the volatile distribution in all Earth and TRAPPIST-1 e, f, and g simulations. The magma ocean evolves from an initially CO₂-dominated atmosphere into a H₂O-dominated atmosphere and for water-poor scenarios (≤10 TO H₂O) the end state is again a CO₂-dominated atmosphere. This shift in atmosphere composition is accompanied by a feedback in outgassing as outlined by Bower et al. (2019). Generally, the addition of CO₂ increases the water melt fraction. Larger melt fractions lead further to larger partial pressures of H₂O and may also shift the on-set of water outgassing. The impact of O₂ is less clear, because in our simulations O₂ only accumulates as the magma ocean solidifies. Here, the magma ocean lifetime that can be extended with low albedos and/or diffusion limited escape is important. We thus find that an extension of the magma ocean lifetime with a clear-sky albedo reduces the amount of abiotically created O₂ strongly, as already outlined by Barth et al. (2021).

The increase in the water melt fraction due to the H₂O-CO₂ feedback is particularly evident for TRAPPIST-1 e (Fig. 7, right panels). Larger water melt fractions also increase the amount of water sequestered in the solid mantle. Here, the amount of remaining water can be increased by 100% and more for the clear-sky albedo with equal amounts of H₂O and CO₂ compared to a pure H₂O scenario. However, the fraction of initial water retained in the solid mantle still remains low: only between 1% and 6% of initial water mass can be retained in the mantle.

New interior composition constraints for all TRAPPIST-1 planets reveals that at least the inner planets (b,c, d) are apparently dry and exhibit an iron fraction of 25 wt-%. If the sane iron fraction is assumed, than TRAPPIST-1 e has currently a low water content, consistent with a dry evolution scenario (≤10 initial H₂O). We find that for such a dry composition, the TRAPPIST-1 e magma ocean lifetime can be significantly extended due to CO₂ diffusion limited escape. The planetary atmosphere may not completely desiccate even with initial water masses as low as 5 TO H₂O within the first 280 Myrs of its evolution; that is, until it enters its habitable zone. TRAPPIST- 1 e is thus the ideal planet to test formation theories and the impact of atmospheric erosion on the water budget during the early evolution of rocky planets orbiting M dwarfs.

In this work, we further developed an analytical corrected gray atmosphere model. This model incorporates information from a thermal emission grid obtained via full RT calculations. Compared to the traditional gray atmosphere model, this corrected gray atmosphere represents a significant improvement. It yields sufficiently similar results for all investigated scenarios in this work compared to the RT atmosphere model, including the special scenario for TRAPPIST-1 e. Both methods are computationally efficient, however, the corrected gray model is 10–26 times faster compared to the RT atmosphere model (Appendix D).

5 Discussion

The magma ocean model MagmOc2.0 that is introduced in this work is designed to facilitate efficient testing of various scenarios and assumptions. This capability is particularly valuable in the era of detailed characterization of rocky planets beyond our Solar System. MagmOc2.0 provides a fast approach to study the evolution of rocky planets in a large parameter space. On the other hand, MagmOc2.0 can benefit from more complex magma ocean and atmosphere models, leading to continuous improvement of this model for better overall understanding of planetary processes.

In this work, specifically, MagmOc2.0 contributes to a better understanding of the potential final atmospheric state and its impact on inferred interior structure for the potentially habitable TRAPPIST-1 e, f, and g planets. This research is crucial for unraveling the complex interplay between planet formation, planetary structure, atmospheric composition and ultimately their potential habitability. TRAPPIST-1 e may be very volatilepoor compared to TRAPPIST-1 f and g if we assume an iron mass fraction of 27 wt% as for the inner planets. TRAPPIST- 1 f and g, on the other hand, must have formed volatile-rich to explain large water fractions that we estimate for the same iron mass fraction, resulting potentially in thick CO₂ atmospheres that may prevent the existence of surface liquid water.

If TRAPPIST-1 e is indeed (partly) desiccated, then it may have formed close to the water-ice line, during which the planet migrated rapidly inward, while f and g formed exterior to the ice line and migrated later. Alternatively, the planets did not undergo significant migration, instead the water ice line shifted during the disk evolution outward starting from a location near TRAPPIST- 1 e, leading to the same outcome (Bitsch et al. 2019).

5.1 Mutual benefit between versatile MagmOc2.0 and other magma ocean models

The investigation of mixed outgassing on magma oceans in the habitable zone of TRAPPIST-1 underscores the potential synergies between magma oceans models with a more comprehensive treatment of geophsical processes, more complex models that capture the geochemistry more completely and simpler frameworks such as VPLanet/MagmOc2.0 that can explore a broad range of initial volatile content and the impact of various treatments of atmospheric and outgassing processes. For instance, the incorporation of feedback between H₂O and CO₂ outgassing was motivated by the magma ocean model results of Bower et al. (2019). Despite employing different differential equation and a simpler geophysical model, we generally reproduce the findings of Bower et al. (2019) for Earth (Appendix C.7).

Comparison between results from the Earth magma ocean model of Elkins-Tanton (2008) and MagmOc2.0 further yields agreement in volatile distribution when the same outgassing laws are used that suppress CO₂ outgassing (Table C.1). However, when modern outgassing laws are applied, there is less agreement in volatile distribution, in particular for CO₂. Suppressed CO₂ outgassing leads to 10–100 times higher CO₂ melt fractions at the end of the magma ocean stage compared to a magma ocean evolution, where the majority of CO₂ is outgassed already at the beginning (Table C.3). Thus, MagmOc2.0 can shed light on how initial outgassing conditions can alter the volatile distribution in diverse magma ocean models.

Application of MagmOc2.0 to TRAPPIST-1 e, f, and g reveals that additional CO₂ also significantly influences these planets’ volatile distributions, resulting in increased abiotic O₂ build-up and a higher remaining water fraction in the solid mantle at the end of the magma ocean stage. Furthermore, we find that the inclusion of CO₂ does not consistently lead to delayed H₂O outgassing, as inferred by Bower et al. (2019) for an Earth-like scenario. Instead, the onset of H₂O outgassing may be delayed, start even earlier, or remain unaffected compared to simulations without CO₂. In our model, it is the radiative properties of the mixed H₂O-CO₂ atmosphere at the time of mantle solidification that primarily determine when H₂O outgassing occurs. This insight is facilitated by exploring more volatile scenarios with different atmosphere models (Appendices C.4 and C.5).

We also demonstrate that MagmOc2.0 can identify unexpected feedbacks between atmospheric erosion and mantle solidification due to diffusion limited escape in mixed H₂O-CO₂ outgassing. For TRAPPIST-1 e with 10 TO H₂O initially, the reduction of the runaway greenhouse radiation limit by 10s of W/m² as water erosion from the atmosphere substantially prolongs the magma ocean stage and also diminishes water sequestration in the mantle.

We note that this special case was not identified in a similar magma ocean study with mixed H₂O-CO₂ outgassing for the TRAPPIST-1 planets conducted by Krissansen-Totton & Fortney (2022). Their model did not consider the reduction of the runaway greenhouse limit with CO₂, which is a prerequisite for this scenario to occur. Additionally, it does not account for changes in the volatile distribution due to the evolution in atmospheric composition. Apart from these differences, the model of Krissansen-Totton & Fortney (2022) yields similar magma ocean solidification lifetimes for TRAPPIST-1 e, f, and g compared to our study.

5.2 The concept of solidification time

When comparing different magma oceans simulations, it is important to discuss the concept of ‘solidification time’ with care. In MagmOc, the magma ocean effectively ends with the convectively dominated viscosity regime. This point is reached when the critical the surface melt fraction ψ = 0.4, following Lebrun et al. (2013) with surface temperatures of 1650 K. At this point 99% of the planetary radius is solidified. In principle, the remaining melt is no longer well mixed and in full equilibrium with the atmosphere. For the last % of radius solidification, thus, MagmOc switches to a high viscosity treatment of the mantle with no further outgassing but with further surface cooling and ongoing atmosphere erosion. Complete solidification then occurs when the solidification radius is equal to the planetary radius, which is in our model typically for surface temperatures of 1400 K (see Sect. 2). Thus, solidification temperatures in our model effectively comprises between 1650 and 1400 K.

Since our magma ocean effectively ends at ψ = 0.4 and complete solidification occurs shortly thereafter, MagmOc2.0 stops earlier and at higher surface temperatures of 1400 K compared to the model of Lichtenberg et al. (e.g., 2021). In this model solidification is defined when the rheological front reaches the surface with ψ = 0.22. Since these authors present the full thermal evolution for pure H₂O and CO₂ atmospheres, respectively, we can still compare MagmOc2.0 to their model state at ψ = 0.4 which is reached in their models with surface temperatures between 1700 K and 1500 K, depending on atmosphere composition, with a similar strong temperature decrease in this regime as also observed in our models (Lichtenberg et al. 2021, Fig. 5). Thus, we identify general agreement for the thermal evolution of our Earth magma ocean model and that of Lichtenberg et al. (2021) (Appendix C.2).

Our effective solidification temperature of 1400 K is similar to that of the model of Krissansen-Totton & Fortney (2022) and only 100–150 K cooler compared to the solidification temperatures in the models of Nikolaou et al. (2019) and Hamano et al. (2013), respectively. In the magma ocean model of Elkins-Tanton (2008) solidification is assumed when 98% of the planet’s radius is solid, which aligns well with our choice of ψ = 0.4 with 99% mantle solidification, allowing comparisons between MagmOc2.0 and their Earth simulations (Appendix C).

5.3 The impact of clouds

Efficient cooling of a magma ocean requires a large net flux, F_Net, leaving the top of the atmosphere, which is the difference between OLR and absorbed stellar atmosphere (ASR). Our simulation results show that mainly two effects determine the magma ocean lifetime: The presence of water vapor and scattering on top of the atmosphere via Rayleigh scattering and clouds, encapsulated in our work via the bond albedo α. The presence of water “sets” the OLR, the albedo α sets the ASR.

When comparing simulations with α = 0.75 (Venus-like scattering) to α = 0 (cloud-free scattering), we find that the cloud-free assumption yields a much larger magma ocean lifetime for all planets even without adding CO₂. The addition of CO₂ is in this respect of minor importance, it may prolong, however, the timescale of water loss as outlined in the next section. We further note that Marcq et al. (2017) propose in their model that the inclusion of clouds lowers the runaway greenhouse radiation limit to 197 W/m², whereas for the majority of our simulations the limit of 280 W/m holds unless the atmosphere is highly water-depleted. Kopparapu et al. (2013), on the other hand, argue that water clouds can be neglected for the calculation of the OLR.

In the clear-sky simulations, the magma ocean lifetime is so prolonged that all planets may enter the habitable zone⁷ before complete mantle solidification for water-rich composition (>10 TO H₂O). This result is of outermost importance because many geophysical models assume by default that water can condense out if the incoming stellar flux is low enough. While this is a safe assumption for Earth, this is apparently no longer true for planets that orbit an M dwarf star.

Instead, volatile-rich planets can be ‘nominally’ in the habitable zone, but condensation on the surface is not yet possible. In addition, atmospheric water loss continues for several hundred years after the habitable zone entrance for TRAPPIST-1 e. In principle, this result is similar to that of Turbet et al. (2021) who showed that even in the Solar system, Venus probably never had a chance to form a surface ocean because thick water clouds prevented efficient cooling of the surface. On the other hand, Yang et al. (2014) showed that for tidally locked planets around M dwarf stars, the high albedo from water clouds forming predominantly over the dayside of an Earth-like planet may allow the inner edge of the habitable zone to shift closer to its host star.

In this work, we simplified the problem by adopting two extreme values to outline the scope of the impact. In our model-setup, however, the cloud-free assumption is the most physically consistent, because we assume efficient reevapora-tion after condensation (that is, with a supersaturation ratio of S=1). Furthermore, Kopparapu et al. (2013) point out that water clouds in a water dominated atmosphere have only a minor impact because of the high intrinsic optical thickness of water vapor in the infrared. Conversely, Marcq et al. (2017) find that the runaway greenhouse radiation limit is lowered to 200 W/m², when water clouds are included. CO₂ clouds, on the other hand, may allow for strong scattering. This effect may be relevant in our framework for desiccated CO₂ atmospheres. However, even for relatively low water abundances of x_H2O = 10⁻³, the upper atmosphere is warmed by the latent heat release of water condensation such that CO₂ clouds are not expected to form (see Fig. B.2). In any case, clearly more work is warranted to assess the impact of clouds in particular to investigate if and when the TRAPPIST-1 e, f, and g planets can form surface oceans after the magma ocean stage. This will require further work on the details of precipitation and evaporation of water on rocky planets around M dwarfs. Already for ‘Earth-like’ conditions with temperate temperatures, different cloud models can yield very different results (Sergeev et al. 2022). Furthermore, Sergeev et al. (2024) point out that resolving moist convection on the dayside of potentially habitable planets requires a higher resolution than is adopted so far. In any case, clearly more work is needed to explore the role of clouds on the outgoing long wave irradiation and albedo for the magma ocean stage on rocky planets around red dwarfs.

5.4 Volatile evolution

A magma ocean evolution with a mixed H₂O-CO₂ atmosphere can exhibit feedback between H₂O and CO₂ outgassing that shapes the overall volatile distribution; that is, the distribution of H₂O, CO₂, and O₂. In addition, in a CO₂-dominated atmosphere that is prone to emerge at the beginning and the end of the magma ocean stage, CO₂-diffusion limited escape can limit water loss, and thus prevent complete desiccation on the investigated planets for a longer time. In the following subsections, we discuss the feedback impact for each volatile species in more detail.

5.4.1 H₂O

The fate of water during the magma ocean evolution on rocky planets around M dwarf stars is of utmost importance for constraining their habitability. If there is too little water (<10 TO H₂O), the planets completely desiccate during the magma ocean stage with only a few percent initial water mass remaining in the mantle.

We find that the mantle is driest for our TRAPPIST-1 e simulations compared to simulations of TRAPPIST-1 f and g because planet e is more strongly affected by UV-photolysis of H₂O. Additionally, the mantle tends to become even drier when the magma ocean stage is prolonged, which is the case for the clear-sky (albedo=0) simulations. The desiccation of the mantle of TRAPPIST-1 e can be, however, compensated by large amounts of CO₂ as is evident by comparing the remaining mantle water inventory for the albedo 0 and 0.75 evolution scenario with assumed initial CO₂ mass being equal to the initial H₂O content (Fig. 7). Despite much longer desiccation times, the remaining water content remains of similar order.

A desiccated evolution scenario is further of particular interest for TRAPPIST-1 e because its composition is apparently consistent with a dry Earth-like H₂O content (<<1 wt-%) (see Sect. 3.3). Thus, here our evolution scenarios with an initial water mass of 1–10 TO, estimated for Earth (see Appendix C), are particularly relevant. They show that the atmosphere of TRAPPIST-1 e may retain 5 TO water until the star’s luminosity dropped such that it enters the habitable zone (Fig. 6). Similarly, even if the atmosphere is desiccated, up to 4% of the initial water mass could be stored in the mantle, in principle enabling a secondary surface water reservoir to build up (Godolt et al. 2019). Since such a relatively wet mantle would require large amounts of CO₂, the surface may be Venus-like and thus too warm for surface liquid water. On the other hand, there may be several other ways to retain more water in the melt like delayed H₂O outgassing (Ikoma et al. 2018; Bottinga & Javoy 1990; Lensky et al. 2006) and a basal magma ocean that doesn’t require huge amounts of CO₂ (see also Sect. 5.4). Conversely, the outer planets TRAPPIST-1 f and g may have the opposite problem of having too much water for a habitable surface ocean. A very thick water ocean not only dilutes nutrients, the further enrichment from the crust may also be suppressed due to the high pressures at the bottom of the ocean (Noack et al. 2016).

The result that an extension of the magma ocean stage generally leads to a drier planet is in marked contrasts with the findings from the 0D magma ocean model by Moore et al. (2023). Figure 1 illustrates the reason for this discrepancy: the amount of dissolved water in the magma ocean, our primary variable of integration, $F_{{\text{H}}_2\text{O}}$, is in balance with the outgassed atmosphere. As has previously been pointed out by Barth et al. (2021), as long as such a balance exists, atmospheric escape will drive additional outgassing from the magma ocean to establish a new equilibrium. It is thus mainly the reduction of water loss in the presence of a thick CO₂ atmosphere that can prevent atmosphere desiccation, in particular on TRAPPIST-1 e. It is then the retention of significant amount of water vapor in the atmosphere (>100 bar) that extends the magma ocean lifetime, and not the other way around.

We note that there are alternative mechanisms to ‘save’ water from atmospheric escape during the magma ocean stage apart from CO₂ diffusion limited escape. In MagmOc2.0, water is sequestered by partitioning 1% of available water mass in the solid mantle during solidification $\left(k_{{\text{H}}_2\text{O}}\right)$, consistent with values for relevant minerals such as lherzolite and peridotite (Elkins-Tanton 2008; Johnson & Dick 1992). This water is saved from atmospheric escape because it is no longer part of the fully coupled volatile system. Nevertheless, the assumption of a constant $k_{{\text{H}}_2\text{O}}$ at the bottom of a highly convective magma ocean may be overly simplistic (Ikoma et al. 2018). Additionally, whether crystals settle at the bottom of the magma ocean (fractional solidification), as is assumed in this model, or remain suspended in the melt (batch solidification) may further impact how much water is incorporated into the solid mantle. Moreover, the role of minerals capable of incorporating higher water mass fractions like phyllosilicates has not been accounted for yet (Herbort et al. 2020). Lastly, the efficiency of outgassing in a highly convective magma ocean may be reduced compared to ‘static’ laboratory experiments (Salvador & Samuel 2023).

In this study, we identify an H₂O-CO₂ feedback as one mechanism that increases the amount of water in the mantle. Adding CO₂-outgassing leads to changes in the mean molecular weight of the atmosphere during the magma evolution (Nikolaou et al. 2019; Bower et al. 2019; Lichtenberg et al. 2021): Initially, the magma ocean begins its evolution with a CO₂-dominated rather than a H₂O-dominated one because CO₂ is less soluble in the melt compared to H₂O. One major consequence is that the water melt fraction $F_{{\text{H}}_2\text{O}}$ increases in a magma ocean with a CO₂-dominated atmosphere, leading to increased H₂O out-gassing with higher H₂O partial pressures of H₂O at later times. Simultaneously, because there is more water in the melt, a greater amount of water can be partitioned and sequestered into the solidifying mantle. Generally it holds that the more CO₂ is in the MOA system initially, the more H₂O remains in the mantle. Consequently, the addition of CO₂ results in 15–50% more remaining water in the mantle at the end of the magma ocean stage for all investigated cases.

Other possibilities to retain more water in the melt include the delay or suppression of H₂O outgassing (Ikoma et al. 2018) through mechanisms such a bubble nucleation (Bottinga & Javoy 1990; Lensky et al. 2006). Furthermore, Bower et al. (2022); Maurice et al. (2017) that use more complex models pointed out that the mode of crystallization can significantly increase water retention. In addition, Maurice et al. (2017) suggest that Rayleigh-Taylor instabilities may replenish the mantle with volatiles from the deeper interior. Furthermore, Hier-Majumder & Hirschmann (2017) suggest that volatile-rich may be trapped at the solidification front. Another mechanism to retain water in the melt may also be the early stratification of the magma ocean (O’Rourke 2020; Samuel et al. 2023) such that part of the magma ocean solidifies downward, resulting in a basal magma ocean. If stratification occurs before the majority of H₂O is out-gassed, then a basal magma ocean could, in principle, retain large amounts of water over the course of billion of years (Moore et al. 2023). Therefore, it is very likely that we underestimate the amount of water that can be retained in the planet. Still, we emphasize that increasing the amount of water with the addition of CO₂ will in all these cases also increase the amount of volatiles that can be stored in the planet. It is thus highly important to diagnose in the near-future the CO₂ atmosphere budget on the potentially habitable TRAPPIST-1 planets in the future.

Here, our magma ocean simulations can help as they also incorporate atmospheric escape. Most importantly we find that the inclusion of CO₂ can significantly reduce atmospheric escape compared to a pure H₂O atmosphere, as already pointed out by Moore et al. (2023). On the other hand, Joule heating, which is not considered in this work, may instead lead to increased atmospheric erosion compared to our model, at least for TRAPPIST-1 e (Cohen et al. 2024).

5.4.2 CO₂

The addition of CO₂ has a profound impact on the overall volatile budget, as it is much less soluble in magma and less prone to be partitioned into the solid mantle than H₂O (Table 1). Furthermore, CO₂ can create a ‘diffusion barrier’ for atmospheric water loss as outlined above for water-poor compositions (<10 TO H₂O) and in particular for TRAPPIST-1 e at the inner edge of the habitable zone of its host star. Still, even in this case, generally large amounts of CO₂ are required for a significant impact of the magma ocean; that is, CO₂ mass ≥0.5× initial water mass (see detailed grid). Such extreme scenarios can lead to very thick (≥10³ bar) CO₂-dominated atmospheres. In very volatile-rich cases (>>10 TO H₂O and CO₂), the magma ocean may not even completely solidify, retaining surface temperatures above 1400 K (Fig. 3).

In contrast, Krissansen-Totton & Fortney (2022) assume less extreme values of additional CO₂, with partial pressures of up to 10³ bars, allowing at least the crust to fully solidify. Another difference is that the runaway greenhouse radiation limit depends in our model on the water volume mixing ratio and can be lowered by 10s of W/m², when CO₂ is added (Appendix B.4). Instead, Krissansen-Totton & Fortney (2022) assume that all water condenses out of the atmosphere onto the solidified crust of TRAPPIST-1 e, f, and g, once the ASR drops below the constant runaway greenhouse radiation limit of 282 W/m². Due to the presence of liquid surface water in their model, surface weathering and the carbon cycle remove hundreds of bars of atmospheric CO₂. In our model, the presence of liquid surface water is questionable for the thick remaining CO₂ atmospheres, which appears to be a reasonable outcome for the outer planets TRAPPIST-1 f and g. Furthermore, in contrast to the model of Krissansen-Totton & Fortney (2022), the magma ocean stage in our cloud-free (albedo=0) model can be extended well into the habitable zone. Thus, we caution to automatically assume rainout of water from the atmosphere, only due to low stellar luminosity.

We further note again that the role of clouds (Turbet et al. 2021; Yang et al. 2014) is of outermost importance (see Sect. 5.3). If water cannot condense out on the surface after the magma ocean stage, then the remaining H₂O-CO₂-O₂ atmospheres would only be modified by further atmospheric loss. Constraints of CO₂-H₂O dominated atmospheres on TRAPPIST-1 e, f, and g may thus provide insights into the abundances of volatiles that were present during the magma ocean stage and on conditions to form surface liquid water, as also pointed out by Krissansen-Totton (2023).

5.4.3 O₂

In our work, atmospheric O₂ begins to accumulate as soon as more than 5 TO H₂O initial water are in the system for albedo=0.75, reaching several hundred bars of partial pressure for 100 TO initial H₂O, as has already been outlined by Luger & Barnes (2015) and Barth et al. (2021). O₂ can only be stored in the FeO-buffer in the mantle, which remains accessible only as long as the magma ocean is not solidified. For the high albedo simulations, the magma ocean solidifies within 100 Myrs and thus tends to yield large partial pressures of O₂. Conversely, for magma ocean lifetimes larger than 100 Myrs, fewer than 50 bars of O₂ build up abiotically. Instead, more O₂ is stored in the mantle (see detailed grid). This grid also shows the tendency of more O₂ being stored in the mantle for extended magma ocean lifetimes even for simulations with a high albedo of 0.75 and initial water masses between 5 and 10 TO.

Krissansen-Totton & Fortney (2022) find in their model that even if there is significant abiotically build-up of oxygen, up to 100 bar O₂ can be removed due to reduced outgassing later in the evolution and crustal oxidation. However, the authors also find that in some scenarios 100 bar of atmospheric O₂ may remain on TRAPPIST-1 e, f, and g.

Both models neglect, however, possible O₂ sequestration in the iron-core (Wordsworth et al. 2018). Moreover, abiotic O₂ build-up occurs when the majority of H₂O is outgassed - that is, when the mantle is mostly solidified – already limiting the exchange of volatiles between the atmosphere and the mantle. Therefore, it is questionable if removal of atmospheric O₂ via the even deeper iron-core can play a significant role this late in the magma ocean evolution.

5.5 The need for tighter link to planet formation models

The prospects of obtaining reliable atmospheric constraints for rocky exoplanets in the near future holds great promise. Establishing a coherent link between planet formation, rocky planet evolution (including the crucial magma ocean stage), interior models and atmospheric composition constraints has the potential to unveil the complete history of rocky exoplanets. Studies of the TRAPPIST-1 system are particularly promising because it is postulated to be only weakly modified by late accretion (Raymond et al. 2022).

The water content of the Trappist planets needs to have been delivered during the protoplanetary disk phase. Planet formation during this initial stage can be driven either by pebble accretion (e.g., Ormel et al. 2017; Schoonenberg et al. 2019; Bitsch et al. 2019) or planetesimal accretion (e.g., Coleman et al. 2019), where both scenarios are capable of roughly reproducing the observed structure of the system. Both scenarios require planet migration as a key ingredient due to the observed resonance configuration of the system (Pichierri et al. 2024). This begs the question of whether the planets formed interior or exterior to the water and CO₂ ice lines in the disk, and how their water content is influenced by their formation location and disk cooling. Consequently, planet formation theories can give a range of initial water contents for the TRAPPIST-1 planets, depending on the exact stellar composition of Trappist-1 and on how efficient atmospheric recycling of incoming water rich pebbles is (Müller et al. 2024a).

In the Schoonenberg et al. (2019) scenario, the planets are formed originally at the water ice line and the migrate inward, which sets their water content. The migration direction and speed in comparison with the growth rate then sets the final water content of planets forming exterior to the water ice line (e.g., Bitsch et al. 2019). Planets that start their formation exterior to the water ice line, but then migrate inward very quickly will have a low water content despite their water-rich formation location. Planets that finish their accretion in the outer regions of the disk and that then only migrate across the water ice line once their accretion is complete will have the maximal water content that the system allows, which should be around 35–50% for solar composition, depending on the exact chemical model and stellar composition (Bond et al. 2010; Bitsch & Battistini 2020; Cabral et al. 2023). The exact volatile ratio could also be affected by atmospheric recycling of incoming water ice rich pebbles, where water fractions above 15% are hard to achieve. A high volatile mass fraction appears to be the currently favored scenario for TRAPPIST-1 g inferred from interior models (Unterborn et al. 2018; Barth et al. 2021; Raymond et al. 2022).

To assess the initial CO₂ content of the TRAPPIST-1 planets, the CO₂ evaporation line needs to be considered, which is located further away from the host star due to the lower evaporation temperature of CO₂ compared to water. This difference inherently implies that a planet forming in the outer regions of the disk that accretes CO₂ ice, will also accrete water ice. The ratio of the water to CO₂ content is then set by the initial CO₂ to water ratio of the system and also by the growth and migration behavior of the planet exterior to the CO₂ line, similar to that described above for the water content.

In addition, the effect of pebble evaporation and condensation can dramatically change the water and CO₂ content of growing planets. As pebbles drift inward and cross evaporation fronts, they release their volatile content into the gas phase (e.g., Piso et al. 2015; Aguichine et al. 2020; Schneider & Bitsch 2021), enriching the disk with vapor. This vapor can then diffuse outward and re-condense, resulting in very large fractions of the corresponding molecule in the solids around the evaporation front (Ros & Johansen 2013). This effect can potentially explain the low water to CO ratio of the Comet C/2016 R2 (Mousis et al. 2021). This effect occurs for all volatiles, indicating that planets that form close to the CO₂ evaporation front could have very large CO₂ to water ratios.

Here, we build upon previous work (Bitsch & Battistini 2020) to constrain the composition of the building blocks of the TRAPPIST-1 planets thus to derive constraints on the iron fractions 25 wt% and thus on the volatile composition of the inner planets. Based on our results, TRAPPIST-1 b, c, and d already formed volatile-poor and are thus not expected to currently exhibit any atmosphere.

In contrast to that TRAPPIST-1 e, f, and g started with volatile-rich building blocks, of which they may have lost a significant fraction already before the magma ocean phase. However, of these planets TRAPPIST-1 e is particularly prone to complete desiccation during the magma ocean evolution - in particular, when the already long magma ocean stage of 50 Myrs is extended even further. Interior models confirm that currently measured radii and masses are compatible with a dry interior and (moderately) thick CO₂ atmosphere of 100 bar. TRAPPIST-1 g, on the other hand, may exhibit up to 17 wt% of water, which indicates an evolution that has allowed it to retain significant fractions of volatiles until today.

This work thus shows that a strong link between CO₂ and H₂O content during formation that is modified during the rocky planet evolution justifies the need for a coherent chain between long term planet evolution models, current atmosphere characterization efforts of rocky exoplanets, and planet formation theories.

The TRAPPIST-1 planetary system may be key here to compare the outcome of three different scenarios in volatile accretion during formation and their subsequent loss or retain-ment via JWST observations and future missions like LIFE, ARIEL. PLATO may be highly useful to constrain extended steam atmospheres (up to 10%) in young rocky exoplanets that are at the end of their magma ocean stage (Turbet et al. 2019; Schlecker et al. 2024).

6 Conclusion

MagmOcV2.0 can provide novel insights into feedback between H₂O and CO₂ for outgassing in magma oceans on Earth-mass rocky planets in the habitable zone around TRAPPIST-1. Together with modeling of the refractory and volatile abundances during formation and modeling of the current interior composition, we derived the following key insights for the TRAPPIST-1 system:

• A compositional model adjusted by the measured metallicity of TRAPPIST-1 yields a dry composition with an iron fraction of 27 wt% for the inner TRAPPIST-1 planets (b, c, and d), which is compatible with the measured radii and masses from Agol et al. (2021). This result is also compatible with recent JWST measurements that may indicate the absence of a substantial atmosphere on TRAPPIST-1 b and c (Zieba et al. 2023; Greene et al. 2023).

• Assuming an iron fraction of 27 wt%, interior structure modeling for TRAPPIST-1 f and g yields a substantial water fraction of 2.3 wt%, and 3.8 wt%, respectively. This result aligns well with a volatile-rich formation environment that is inferred for these planets in this work and also aligns with formation modeling (see also Barth et al. 2021; Raymond et al. 2022; Miguel et al. 2020; Schoonenberg et al. 2019).

• Without additional sink terms, thick (>1000 bar), mixed CO₂-H₂O atmospheres remain at the end of the magma ocean stage for volatile-rich scenarios ( >10 TO initial water mass) that are favored for TRAPPIST-1 f and g.

• More than 1000 bar of CO₂ retained after the magma ocean stage may contribute up to 2% to the total planetary radius if water vapor is removed after the magma ocean stage. An H₂O-dominated atmosphere would comprise even up to 10% of the planetary radius, as has already been pointed out by Turbet et al. (2019); Schlecker et al. (2024). While such extended steam atmospheres are incompatible with expectations for the evolved TRAPPIST-1 system, they may be identifiable for younger systems with future space missions like PLATO.

• In all cases, we find that adding CO₂ increases the H₂O mass fraction in the melt. Thus, between 15 and 100% more water can be retained in the solidified crust, which results in the sequestration of between 2 and 6% of the initial water mass in the solidified mantle for TRAPPIST-1 e, f, and g.

• TRAPPIST-1 e may have formed, like TRAPPIST-1 f and g, in a volatile-rich formation environment. If, however, an iron fraction of 27 wt% is assumed, as was derived for the inner planets, then its current radius and mass is consistent with a more desiccated evolution scenario (water content of ${0.6}_{-0.6}^{+2.1}$ wt%). The values also include Earth-like water content (≤ 10 TO H₂O).

• TRAPPIST-1 e magma ocean simulations with Earth-like water content are very sensitive to assumptions about scattering (albedo) as well as to changes in atmospheric water loss brought about via diffusion-limited escape in a CO₂ background atmosphere. Here, CO₂ can prevent, or at least delay, complete atmospheric water loss, even for a long magma ocean stage of 350 Myrs. The atmospheric end result for TRAPPIST-1 e after the magma ocean stage could thus be a either a dry CO₂ or a mixed H₂O-CO₂ atmosphere with surface pressure of several 100 to 1000 bar. The mixed atmosphere could, in principle, further evolve into a more Earth-like state if water condenses out of the atmosphere and if up to 1000 bar of atmospheric CO₂ can be sequestered in the mantle, as is proposed, for example, by Krissansen-Totton & Fortney (2022).

• TRAPPIST-1 e thus emerges as a “Rosetta stone” for deciphering formation as well as evolution processes for rocky planets in the habitable zone of active M dwarfs.

7 Outlook

In this work, we have illustrated how combined H₂O-CO₂ outgassing can alter the magma ocean evolution and volatile distribution on rocky planets in the habitable zone of their Mdwarf host stars, using TRAPPIST-1 e, f, and g as examples. However, future investigation is necessary to test the impact of other factors that we either neglected or assumed to be constant. For example, we kept the albedo fixed at α = 0.75 and α = 0, respectively), neglecting albedo variations as the atmosphere evolves (Pluriel et al. 2019). We also assumed in the vertically extended atmosphere immediate reevaporation after condensation, which is a strong simplification, in particular toward the end of the magma ocean stage. Additionally, exploring the evolution of the mantle’s redox state (Katyal et al. 2020) will be important because it may result in outgassing of reduced volatiles, such as H₂ and CO (Ortenzi et al. 2020; Deng et al. 2020).

H₂ is proposed to be dissolved into the mantle and core during the planetary accretion state, establishing a more direct connection between magma ocean and formation models (Johansen et al. 2023; Young et al. 2023, e.g.,). A significant atmospheric H₂ content during the magma ocean phase may lead to even more substantial outgassing feedback effects, as is tackled in this study, because it would significantly reduce the mean molecular weight of the atmosphere. Moreover, H₂ is a potent greenhouse gas that can significantly prolong the magma ocean stage (Lichtenberg et al. 2021) even more. At the same time, it can allow for surface liquid water on rocky planets at the outer edge of the habitable zone, such as Mars (Pahlevan et al. 2022). Our versatile outgassing formalism, coupled with the RT atmosphere model, allows for the adoption of other outgassing scenarios.

While our current setup lacks a thermal boundary layer (as in Barth et al. 2021), future investigations are planned that will incorporate now more advanced geophysical VPLanet modules, like ThermInt, for the exploration of an Earth-like planet with a stagnant lid configuration (Driscoll & Bercovici 2014). Additionally, plans include integration of the “mush stage” described by Lebrun et al. (2013) after ψ = 0.4 is reached. Similarly, we recognize the necessity for a modification of the atmospheric equation of state that is in this work based on the ideal gas law at the surface for particularly volatile-rich evolution scenarios; for example, those with an initial volatile mass much larger than 10 TO, which is p_{sur f} > 1000 bar.

Still, our work strongly indicates that TRAPPIST-1 e is a particularly interesting planet, the atmospheric composition of which may shed light on key processes of formation and early rocky planet evolution. Furthermore, it will be worthwhile to compare TRAPPIST-1 e with TRAPPIST-1 g and f, which are expected to be more volatile-rich than e. Here, the CO₂ atmospheric content will be key, to investigate if it represents the primordial volatile content or if it is modified after the magma ocean phase; for example, by a carbonate-silicate cycle.

In summary, MagmOc2.0 presents a valuable tool to gain insights into the early stages of rocky planet evolution and to understand the impact of various processes on the outgassed secondary atmosphere. These insights may facilitate rocky planet atmosphere characterization with the James Webb Space Telescope and future missions like the Habitable World Observatory and LIFE, and improve our understanding of the formation and evolution history of rocky planets (Bonati et al. 2019; Turbet et al. 2019; Way et al. 2022; Krissansen-Totton & Fortney 2022; Schlecker et al. 2024).

Data availability

Appendix D is available at https://doi.org/10.5281/zenodo.14442985.

Acknowledgements

We thank the anonymous referee for a highly constructive discussion that improved this work substantially. We thank Tim Lichtenberg and his group for excellent discussions about the opacities and atmospheres on rocky planets. We also thank Paul Mollière and Thomas Henning for support during the master thesis of P.B. that was the starting point for the whole MagmOc development. L.C. acknowledges support by the DFG priority programme SP1833 “Building a habitable Earth” Grant CA 1795/3 and the Royal Astronomical Society University Fellowship URF R1 211718 hosted by the University of St Andrews. R.B. acknowledges support from NASA grants numbered 80NSSC20K0229, 80NSSC18K0829 and 80NSSC18K0261. P.B. acknowledges financial support from the Austrian Academy of Science. Ch.H. and L.C. acknowledge funding from the European Union H2020-MSCA-ITN-2019 under Grant Agreement no. 860470 (CHAMELEON). L.N. acknowledges funding from the European Union (ERC, DIVERSE, 101087755).

Appendix A Derivation of the coupled H₂O-CO₂ outgassing differential equations

We present here an improvement of the volatile model used by Barth et al. (2021) for MagmOcl.O, which originally adopted the magma ocean model including outgassing of Schaefer et al. (2016). To derive the coupled differential equations that drive multi-component outgassing, we use the mass balance in the magma ocean and atmosphere as a starting point as by Barth et al. (2021, Eq. 12) (see also Eq 1).

For clarity, we outline again the setup of the volatile reservoirs. Each volatile i (here H₂O and CO₂ ) is stored across the magma ocean and atmosphere system (MOA) as follows (Fig. 1): a) in the crystallized portion of the magma ocean $\left(M_i^{\text{crystal\hspace{0.17em}}}\right)$, b)in the liquid magma phase $\left(M_i^{\text{liq}}\right)$, c) and in the atmosphere $\left(M_i^{\text{atm\hspace{0.17em}}}\right)$. Additionally, we define, F_i, as the mass fraction of the volatile i solved in the melt. We also assume a constant, mantle-averaged partition coefficient k_i between the melt and the crystal phase. The latter is assumed to sink to the bottom and form the solidified mantle. For the partition coefficient, we adopt $k_{{\text{H}}_2\text{O}}$ from Schaefer et al. (2016) and feco₂ from Lebrun et al. (2013) (Table 1).

Lastly, the mass of volatile i in the atmosphere is calculated via $$M_i^{\text{atm}}=\frac{4\pi r_{\text{p}}^2}{ɡ}{p}_{i,mass},$$(A.1)

where p_i,mass is mass weighted pressure with respect to the partial pressure p_i,part of volatile i (see e.g., Bower et al. (2019)): $$p_{i,\text{mass}}=\frac{{\mu }_i}{{\overline{\mu }}_{\text{atm}}}{p}_{i,\text{part}},$$(A.2)

where µ_i represents the molecular mass of the specific volatile and ${\overline{\mu }}_{atm}$ stands for the mean atmospheric molecular mass, considering the combination of all outgassed volatiles. It is worth noting that in the water-dominated steam atmosphere setup of MagmOc1.0 (Barth et al. 2021), the authors assumed $p_{{\text{H}}_2\text{O},mass}=p_{{\text{H}}_2\text{O}}$ because it held true for most of the simulation time that ${\mu }_{{\text{H}}_2\text{O}}\approx {\overline{\mu }}_{atm}$. This assumption no longer holds in cases where a H₂O dominated atmosphere evolves into CO₂ dominated atmosphere as simulated in this work. We also neglect here O₂ because atmospheric O₂ is in this model typically produced after the end of the magma ocean and therefore does not enter the mass balance equations describing balance between outgassing and dissolved volatile mass in the magma ocean phase.

Thus, the volatile mass balance can be summarized as (see Table 2 for an overview of all the components of the volatile model.): $$\begin{array}{c}{M}_i^{\text{moa}}=M_i^{\text{crystal}}+M_i^{\text{liq}}+M_i^{\text{atm}}\\ =k_i{F}_i{M}^{\text{crystal}}+F_i{M}^{\text{liq}}+\frac{4\pi r_{\text{p}}^2}{ɡ}{p}_{i,\text{\hspace{0.17em}mass,\hspace{0.17em}}}\end{array}$$(A.3)

where M^crystal denotes the mass of the crystallized part of the magma ocean and M^liq is the mass of the liquid phase of the magma ocean. We further assume that the partial pressure of the volatile at the surface is determined by the fraction of the volatile in the melt F_i (Elkins-Tanton 2008; Schaefer et al. 2016; Lichtenberg et al. 2021).

In Barth et al. (2021), the water outgassing law of Schaefer et al. (2016), based on the laboratory data of Papale (1997), was used. Therefore, we employ for the pure water outgassing the same relation to compare with MagmOc1.0: $$p_{{\text{H}}_2\text{O}}={\left(\frac{F_{{\text{H}}_2\text{O}}}{3.44\times {10}^{-8}}\right)}^{1/0.74}[\text{Pa}].$$(A.4)

For a mixed atmosphere model, Elkins-Tanton (2008) assumed the following solubility laws for H₂O and CO₂ that were also based Papale (1997) but have a form that suppresses outgassing below a certain volatile melt fraction: $$p_{{\text{H}}_2\text{O}}={\left(\frac{F_{{\text{H}}_2\text{O}}-3\times {10}^{-3}}{2.08\times {10}^{-6}}\right)}^{1/0.52}[\text{Pa}]$$(A.5) $$p_{{\text{CO}}_2}={\left(\frac{F_{{\text{CO}}_2}-5\times {10}^{-4}}{2.08\times {10}^{-6}}\right)}^{1/0.45}[\text{Pa}].$$(A.6)

Figure A.1 provides an overview of solubility laws used in the literature. The H₂O solubility laws agree with each other within one order of magnitude for 260 bar≤ p < lO⁴ bar. Consequently, it is not surprising that magma oceans with H₂O outgassing yield generally similar results (see Appendix C). However, significant deviations are apparent between the CO₂ solubility law used by Elkins-Tanton (2008) and those used by Lichtenberg et al. (2021) and Nikolaou et al. (2019). Conversely, the CO₂laws used by Lichtenberg et al. (2021) and Nikolaou et al. (2019) generally agree with each other.

In MagmOc1,0 the mass balance equations (Eq. A.3) were solved for $F_{{\text{H}}_2\text{O}}$ at every time step using a root finding method. Using this approach for two or more volatiles, however, is cumbersome and numerically costly.

In the new model MagmOc2.0, we opt instead to solve by numerical integration: $$F_{{\text{H}}_2\text{O}}(t)=F_{{\text{H}}_2\text{O}}(0)+\frac{d{F}_{{\text{H}}_2\text{O}}(t)}{dt}\text{Δ}t$$(A.7) $$F_{{\text{CO}}_2}(t)=F_{{\text{CO}}_2}(0)+\frac{d{F}_{{\text{CO}}_2}(t)}{dt}\mathrm{\Delta }t$$(A.8)

where the starting points of the integration, $F_{{\text{H}}_2\text{O}}(0)$ and $F_{{\text{CO}}_2}(0)$, are input parameters for MagmOc2.0 and are derived via root finding only once when preparing the simulations.

The derivatives of the melt fractions $F_i{}^{\prime }$ can be found by differentiating Eqs. A.3, assuming that the total volatile mass in the system, $M_i^{moa}$, has to remain constant when a new equilibrium in mass for H₂O and CO₂ is established between all reservoirs.

For mixed H₂O-CO₂ outgassing, we also account for changes in the mean molecular mass of the atmosphere ${\overline{\mu }}_{atm}$, given by $${\overline{\mu }}_{atm}=\frac{p_{{\text{H}}_2\text{O}}\cdot {\mu }_{{\text{H}}_2\text{O}}+p_{{\text{CO}}_2}\cdot {\mu }_{{\text{CO}}_2}}{\underset{p_{surf}}{\underbrace{p_{{\text{H}}_2\text{O}}+p_{{\text{CO}}_2}}}},$$(A.9)

where $p_{{\text{H}}_2\text{O}}+p_{{\text{CO}}_2}$ is the surface pressure of the complete atmosphere p_surf.

To make the derivation more manageable, we introduce placeholders for several terms of the mass balance equation, which are listed in Table A.1. The derivatives of the mass balance equations for H₂O and CO₂ are then expressed by $$\frac{d}{dt}\left(a_1+b_1\cdot ɡ+h\cdot ɡ+\frac{d\cdot m}{w_1{\left(\frac{l+u}{v}\right)}^{ex{p}_2}\cdot {\left(\frac{ɡ+e}{f}\right)}^{-ex{p}_1}+1}\right)=0$$(A.10) $$\frac{d}{dt}\left(a_2+b_2\cdot l+h\cdot l+\frac{d\cdot m}{w_2{\left(\frac{ɡ+e}{f}\right)}^{ex{p}_1}\cdot {\left(\frac{l+u}{v}\right)}^{-ex{p}_2}+1}\right)=0$$(A.11)

We point out that m′ is a place holder for the derivative of the total surface pressure p_surf of the H₂O- CO₂ atmosphere and thus depends on both, $F_{{\text{H}}_2\text{O}}^{\prime }={ɡ}^{\prime }$ and $F_{{\text{CO}}_2}^{\prime }=l^{\prime }$. Its derivative is (according to the substitution we chose in Table A.1): $$\begin{array}{l}\frac{dm}{dt}=\frac{d}{dt}{\left(\frac{ɡ+e}{f}\right)}^{ex{p}_1}+\frac{d}{dt}{\left(\frac{l+u}{v}\right)}^{ex{p}_2}\\ =ex{p}_1\left(\frac{{ɡ}^{\prime }}{f}\right){\left(\frac{ɡ+e}{f}\right)}^{ex{p}_1-1}+ex{p}_2\left(\frac{l^{\prime }}{v}\right){\left(\frac{l+u}{v}\right)}^{ex{p}_2-1}\end{array}$$(A.12)

We now replace m′ in the derivation of the mass balance equations with this resolved expression and rearrange the two mass balance equations to gather all terms that contain $F_{{\text{H}}_2\text{O}}^{\prime }={ɡ}^{\prime }$ and $F_{{\text{CO}}_2}^{\prime }=l^{\prime }$ on the right hand side and all other terms on the left hand side. We thus derive for H₂O: $$\begin{array}{l}{ɡ}^{\prime }[b_1+h+d\cdot m\cdot \frac{{\exp }_1{w}_1\frac{1}{f}\cdot {\left(\frac{l+u}{v}\right)}^{{\exp }_2}}{{\left(\frac{ɡ+e}{f}\right)}^{{\exp }_1+1}}\cdot {\left(\frac{w_1{\left(\frac{l+u}{v}\right)}^{{\exp }_2}}{{\left(\frac{ɡ+e}{f}\right)}^{{\exp }_1}}+1\right)}^{-2}\\ +\frac{d\cdot {\exp }_1\frac{1}{f}\cdot {\left(\frac{ɡ+e}{f}\right)}^{{\exp }_1-1}}{w_1{\left(\frac{l+u}{v}\right)}^{{\exp }_2}{\left(\frac{ɡ+e}{f}\right)}^{-{\exp }_1}+1}]\\ +l^{\prime }[-d\cdot m\cdot \frac{{\exp }_2{w}_1\frac{1}{v}\cdot {\left(\frac{l+u}{v}\right)}^{{\exp }_2-1}}{{\left(\frac{ɡ+e}{f}\right)}^{{\exp }_1}}\cdot {\left(\frac{w_1{\left(\frac{l+u}{v}\right)}^{{\exp }_2}}{{\left(\frac{ɡ+e}{f}\right)}^{{\exp }_1}}+1\right)}^{-2}\\ +\frac{d\cdot {\exp }_2\frac{1}{v}{\left(\frac{l+u}{v}\right)}^{{\exp }_2-1}}{w_1{\left(\frac{l+u}{v}\right)}^{{\exp }_2}{\left(\frac{ɡ+e}{f}\right)}^{-{\exp }_1}+1}]\\ =-\left[a_1^{\prime }+ɡ\cdot b_1^{\prime }+ɡ\cdot h^{\prime }\right].\end{array}$$(A.13)

and accordingly for CO₂: $$\begin{array}{l}{l}^{\prime }[b_2+h+d\cdot m\cdot \frac{{\exp }_2{w}_2\left(\frac{1}{v}\right){\left(\frac{ɡ+e}{f}\right)}^{{\exp }_1}}{{\left(\frac{l+u}{v}\right)}^{{\exp }_2+1}}\cdot {\left(\frac{w_2{\left(\frac{ɡ+e}{f}\right)}^{{\exp }_1}}{{\left(\frac{l+u}{v}\right)}^{{\exp }_2}}+1\right)}^{-2}\\ +\frac{d\cdot {\exp }_2\left(\frac{1}{v}\right){\left(\frac{l+u}{v}\right)}^{{\exp }_2-1}}{w_2{\left(\frac{ɡ+e}{f}\right)}^{{\exp }_1}{\left(\frac{l+u}{v}\right)}^{-{\exp }_2}+1}]\\ +{ɡ}^{\prime }[-d\cdot m\cdot \frac{{\exp }_1{w}_2\frac{1}{f}{\left(\frac{ɡ+e}{f}\right)}^{{\exp }_1-1}}{{\left(\frac{l+u}{v}\right)}^{{\exp }_2}}\cdot {\left(\frac{w_2{\left(\frac{ɡ+e}{f}\right)}^{{\exp }_1}}{{\left(\frac{l+u}{v}\right)}^{{\exp }_2}}+1\right)}^{-2}\\ +\frac{d\cdot {\exp }_1\frac{1}{f}{\left(\frac{ɡ+e}{f}\right)}^{{\exp }_1-1}}{w_2{\left(\frac{ɡ+e}{f}\right)}^{{\exp }_1}{\left(\frac{l+u}{v}\right)}^{-{\exp }_2}+1}]\\ =-\left[a_2^{\prime }+l\cdot b_2^{\prime }+l\cdot h^{\prime }\right].\end{array}$$(A.14)

These two equations establish a set of linear equations with respect to ${ɡ}^{\prime }=F_{{\text{H}}_2\text{O}}^{\prime }$ and $l^{\prime }=F_{{\text{CO}}_2}^{\prime }$ of the form: $$A_{{\text{H}}_2\text{O}}{l}^{\prime }+B_{{\text{H}}_2\text{O}}{ɡ}^{\prime }=-C_{{\text{H}}_2\text{O}}$$(A.15) $$A_{{\text{CO}}_2}{l}^{\prime }+B_{{\text{CO}}_2}{ɡ}^{\prime }=-C_{{\text{CO}}_2}.$$(A.16)

Therefore, we can use Cramer’s rule to derive unique solutions for 𝑔′ and l′ as long as $A_{{\text{H}}_2\text{O}}{B}_{{\text{CO}}_2}\ne A_{{\text{CO}}_2}{B}_{{\text{H}}_2\text{O}}$: $${ɡ}^{\prime }=\frac{d{F}_{{\text{H}}_2\text{O}}}{dt}=\frac{C_{{\text{H}}_2\text{O}}{A}_{{\text{CO}}_2}-C_{{\text{CO}}_2}{A}_{{\text{H}}_2\text{O}}}{A_{{\text{H}}_2\text{O}}{B}_{{\text{CO}}_2}-A_{{\text{CO}}_2}{B}_{{\text{H}}_2\text{O}}}$$(A.17) $$l^{\prime }=\frac{d{F}_{{\text{CO}}_2}}{dt}=\frac{C_{{\text{CO}}_2}{B}_{{\text{H}}_2\text{O}}-C_{{\text{H}}_2\text{O}}{B}_{{\text{CO}}_2}}{A_{{\text{H}}_2\text{O}}{B}_{{\text{CO}}_2}-A_{{\text{CO}}_2}{B}_{{\text{H}}_2\text{O}}}.$$(A.18)

Thus, we demonstrate that while the mass balance equation with the parametric equation for volatile outgassing cannot be solved for analytically with respect to $F_{{\text{H}}_2\text{O}}$ and $F_{{\text{CO}}_2}$, their derivatives can. Therefore in MagmOc2.0, we add the two differential equations, Eqs. (A.7) and (A.8), with Eqs. A.17 and A.18. We note that the last term of the sum comprising $A_{{\text{H}}_2\text{O}},\dots ,B_{{\text{H}}_2\text{O}}$ encapsulates changes in surface pressure and also in mean molecular weight (see also Eq. A.2).

We further emphasize, similar to Bower et al. (2019), that the derivation outlined in this section can in principle be generalized to more than two volatiles as long as it is possible to derive a set of n independent linear equations for n volatiles that can be solved for time derivatives of all relevant volatile mass fractions F_{1… n}.

Appendix B A corrected gray atmosphere model for a two component H2O-CO2 atmosphere

As is outlined in Sect. 2.5, the thermal evolution of the planet is determined by the net outgoing flux F_net; that is, the difference between ASR and the thermal flux emitted by the planet, also known as OLR, on top of the atmosphere. In MagmOc1.0 (Barth et al. 2021), the OLR was calculated with a gray atmosphere model: $$F_{OLR}={ϵ}_{{\text{H}}_2\text{O}}\sigma T_{\text{surf\hspace{0.17em}}}^4,$$(B.1)

where σ is the Stefan-Boltzmann constant, and ${ϵ}_{{\text{H}}_2\text{O}}$ is the emissivity of the H₂O atmosphere.

The emissivity of a volatile i can be derived from the optical depth of an atmosphere composed of volatile i with the Rosse- land mean opacity over the infrared wavelength range, κ₀(i), at the reference pressure, p₀. The optical depth, τ_i, further depends on surface gravity, 𝑔, and partial pressure, p_i, at surface as (see e.g., Carone et al. 2014; Elkins-Tanton 2008; Catling & Kasting 2017; Pierrehumbert 2010): $${\tau }_i=p_i\sqrt{0.75\cdot {\kappa }_0(i)/\left(ɡ\cdot p_0\right)}.$$(B.2)

From this relation, the emissivity, є_i, follows as $${ϵ}_i=\frac{2}{{\tau }_i+2}.$$(B.3)

We used this gray formalism as a basis to compare to calculations with full RT in a vertically extended atmosphere and to derive a parametric fit.

B.1 Pure H2O atmosphere for Earth gravity

For a pure H₂O atmosphere, the gray atmosphere analytical model with κ₀(H₂O) = 0.25 m²/kg at p₀ = 1.013 bar and the runaway greenhouse radiation limit of OLR_lim(H₂O) = 282 W/m² by prescribing (Fig. B.1, left) $$F_{OLR}\left({\text{H}}_2\text{O},c_P=\text{const}\right)=\max \left({\text{OLR}}_{lim}\left({\text{H}}_2\text{O}\right),{ϵ}_{{\text{H}}_2\text{O}}\sigma T_{\text{surf\hspace{0.17em}}}^4\right)$$(B.4)

yields an agreement for T_surf = 500 − 4000 K and p_surf = 1 − 26000 bar within one order of magnitude when compared to thermal emission from full RT (Sect. 2.5) in a vertically extended atmosphere (Sect. 2.4) with a constant dry adiabatic lapse rate of $c_{p,{\text{H}}_2\text{O}}=37\text{\hspace{0.17em}J}$ J mole⁻¹ K⁻¹.⁸ The heat capacity of water is, however, not constant, but increases strongly for T > 1000 K, which results in a steeper gradient for temperature profiles for the early stages of the magma ocean evolution (Fig. B.2). The steeper temperature-pressure gradient results in hotter upper atmospheric layers and consequently higher outgoing thermal fluxes (Fig. B.1, left). We can account to first order for increased thermal emission with larger surface temperatures and pressures (Fig. B.1, right) by adopting an emission correction term: $$E^{*}\left({\text{H}}_2\text{O}\right)={\left(\frac{T_{surf}}{1500K}\right)}^{lo{ɡ}_{10}\left(\frac{p_{{\text{H}}_2\text{O}}}{1bar}\right)}.$$(B.5)

The corrected gray atmosphere model for a pure H₂O is then described with $$F_{OLR}\left({\text{H}}_2\text{O}\right)=\max \left({\text{OLR}}_{lim}\left({\text{H}}_2\text{O}\right),{ϵ}_{{\text{H}}_2\text{O}}\sigma T_{surf}^4{E}^{*}\left({\text{H}}_2\text{O}\right)\right).$$(B.6)

B.2 Pure CO2 atmosphere for Earth gravity

For a pure CO₂ atmosphere, we find even larger disagreement between thermal emission based on the gray atmosphere model and the full RT calculations (Fig. B.3, left). The gray atmosphere model (κ₀(CO₂) = 0.001 m²/kg at p₀ = 1.013 bar) yields only good agreement for p_surf = 1 bar and underestimates the OLR by orders of magnitudes for higher surface pressures. Assuming κ₀(CO₂) = 0.005 m²/kg at p₀ = 1.013 bar, as used in Lebrun et al. (2013); Elkins-Tanton (2008), yields even worse results and underestimates the OLR by at least one order of magnitude also for p_surf = 1 bar.

The following analytical prescription is, however, able to mimic the OLR from the full RT calculations in a vertically extended pure CO₂ atmosphere to good accuracy for surface pressures between 1 − 26000 bar and T_surf = 500 − 4000 K. Here, we set the CO₂ radiation limit OLR_lim(CO₂) = 64 W/m² (Fig. B.3, right): $$F_{OLR}\left({\text{CO}}_2\right)=\max \left(OL{R}_{lim}\left({\text{CO}}_2\right),{ϵ}_{{\text{CO}}_2}^{*}\sigma T_{surf}^4\cdot E^{*}\right),$$(B.7)

and used an adapted optical depth for CO₂ and the partial pressure, $P_{{\text{CO}}_2}$: $${\tau }_{{\text{CO}}_2}^{*}={\left(\frac{p_{{\text{CO}}_2}}{p_{ref}\left({\text{CO}}_2\right)}\right)}^{0.75}\cdot p_{ref}\left({\text{CO}}_2\right)\sqrt{0.75\cdot {\kappa }_0\left({\text{CO}}_2\right)/\left(ɡ\cdot p_0\right)},$$(B.8)

where p_ref (CO₂)=1 bar is the reference pressure for the fit. We then derived ${ϵ}_{{\text{CO}}_2}^{\ast }$via Eq. (B.2) by substituting ${\tau }_{{\text{CO}}_2}$ with ${\tau }_{{\text{CO}}_2}^{\ast }$.

As in the pure H₂O-atmosphere case, the steepness of the pressure temperature profile for large surface temperatures in a CO₂ atmosphere requires a correction term for the emission, which we fit as $$E^{*}\left({\text{CO}}_2\right)={\left(\frac{T_{surf}}{1200K}\right)}^{lo{ɡ}_{10}\left(\frac{p_{{\text{CO}}_2}}{1bar}\right)}.$$(B.9)

B.3 Mixed H2O-CO2 atmosphere for Earth gravity

For mixed H₂O-CO₂ atmospheres, the following parametric approximation yields thermal emission within one order of magnitude compared to the full RT model (Fig. B.4) for T_surf = 500 K - 4000 K and p_surf = 1 bar - 2600 bar using: $$F_{\text{OLR}}\left(x_{{\text{H}}_2\text{O}}\right)=\max \left({\text{OLR}}_{\lim }\left({\text{x}}_{{\text{H}}_2\text{O}}\right),{ϵ}_{\text{Mix}}^{*}\sigma {\text{T}}_{\text{surf}}^4\cdot {\text{E}}_{\text{Mix}}^{*}\right),$$(B.10)

with $${ϵ}_{Mix}^{*}=\frac{2}{{\tau }_{Mix}+2},$$(B.11)

where τ_Mix is the sum of the modified optical depths derived in the previous two subsections: $${\tau }_{Mix}={\tau }_{{\text{H}}_2\text{O}}+{\tau }_{{\text{CO}}_2}^{*}.$$(B.12)

The runaway greenhouse radiation radiation, or H₂O OLR limit, decreases when CO₂ is added to a H₂O-dominated atmosphere, as already pointed out by Goldblatt et al. (2013). We systemti- cally investigate this effect for a vertically extended H₂O-CO₂ atmosphere. We find that the following fit to the H₂O OLR limit $OL{R}_{Lim}\left(x_{{\text{H}}_2\text{O}}\right)$ with respect to volume mixing ratio $\left(x_{{\text{H}}_2\text{O}}\right)$ and for Earth surface gravity (𝑔 = 9.81 m/s) yields agreement to first order (Fig. B.5): $${\text{OLR}}_{\text{Lim}}\left({\text{x}}_{{\text{H}}_2\text{O}}\right)=4.9\cdot {\log }_{10}{\left({\text{x}}_{{\text{H}}_2\text{O}}\right)}^2+66.9\cdot {\log }_{10}\left({\text{x}}_{{\text{H}}_2\text{O}}\right)+282,$$(B.13)

where $x_{{\text{H}}_2\text{O}}\le {10}^{-6}$ was set to be identical to a pure CO₂ atmosphere. That is, we assume for such low volume mixing ratios of water in a CO₂-dominated atmosphere that the water contribution can be neglected (within 10% accuracy) when simulating the thermal evolution of the magma ocean (see also Fig. 2, which shows that both the pressure-temperature profile and the emission are very close to that of a 100% CO₂ atmosphere for $x_{{\text{H}}_2\text{O}}={10}^{-5}$).

We adopt the following correction in emission, using the volume fraction of H₂O, $x_{{\text{H}}_2\text{O}}=p_{{\text{H}}_2\text{O}}/p_{surf}$, in a two-component atmosphere setup with $p_{{\text{H}}_2\text{O}}+p_{{\text{CO}}_2}=p_{surf}$: $$\begin{array}{l}{E}_{Mix}^{*}={\left(\frac{T_{surf}}{1200K}\right)}^{lo{ɡ}_{10}\left(\left(\frac{p_{surf}}{1bar}\right)\cdot \left(1-x_{{\text{H}}_2\text{O}}\right)+x_{{\text{H}}_2\text{O}}\right)}\cdot \left(1-x_{{\text{H}}_2\text{O}}\right)+\\ {\left(\frac{T_{surf}}{1500K}\right)}^{lo{ɡ}_{10}(\left(\frac{p_{surf}}{1bar}\right)x_{{\text{H}}_2\text{O}}+\left(1-x_{{\text{H}}_2\text{O}})\right)}\cdot x_{{\text{H}}_2\text{O}}\end{array}$$(B.14)

Two examples of the resulting grid emission and corrected gray emission are shown in (Fig. B.4), using κ₀(CO₂) = 0.0008 m²/kg and κ₀(H₂O) = 0.25 m²/kg at p₀ = 1.013 bar. These show that the corrected gray atmosphere model yields agreement to first order to the full RT calculations even for very different atmosphere mixing. We note that this two-component H₂O-CO₂ atmosphere model also encapsulates 100% H₂O and CO₂ compositions by setting $x_{{\text{H}}_2\text{O}}$ equal to 1 or 0, respectively.

B.4 Generalization of gray atmosphere for rocky exoplanets

So far, we have applied the RT and corrected gray atmosphere model to Earth. In principle, the latter can also be applied to any rocky planet with a different surface gravity 𝑔.

The pressure-temperature profiles will not change with different surface gravities for a given atmospheric composition and surface temperature and pressure set by the hot magma ocean, because the atmospheric lapse rate is determined either by R/c_p,i when the vertical coordinate is chosen, where R is the ideal gas constant and c_p (T) is the combined specific heat capacity of the volatiles in the atmosphere or by latent heat release L_i. Thus, pressure-temperature profiles for rocky planets with the same atmospheric composition, surface temperature, and surface pressure are in this setup identical for different surface gravities. The emission arising from the top of the atmosphere, however, will still change because the planet’s opacity depends on surface gravity (Eq. B.2). Consequently, also the runaway greenhouse limit changes with surface gravity (Fig. B.6, left).

We find that the runaway greenhouse radiation limit for different surface gravities can be parameterized with a simple adjustment to Eq. (B.13) (Fig. B.6, right): $${\mathrm{OLR}}_{\text{Lim}}\left(x_{{\text{H}}_2\text{O}},\text{\hspace{0.17em}g}\right)={\mathrm{OLR}}_{\text{Lim}}\left({\text{x}}_{{\text{H}}_2\text{O}}\right)\cdot \left(1+\ln \left(\frac{\text{g}}{{\text{g}}_{\text{Earth\hspace{0.17em}}}}\right)\cdot 0.1\right).$$(B.15)

Appendix C Validation of MagmOc2.0 for an oxidized Earth

We applied our magma ocean to the Earth regime to make a comparison with the previous results of Barth et al. (2021), who also benchmarked MagmOc1.0 for Earth-like conditions for comparison with Hamano et al. (2013) and Elkins-Tanton (2008). We further exploited the versatility of MagmOc2.0 to test three different atmosphere models and their impact on the thermal evolution of the mag, a ocean: the gray formalism as used in MagmOc1.0, an improved gray formalism, and a model based on full RT calculations in a vertically extended atmosphere. Henceforth, we refer to the first model as the “gray” model, to the second as the “corrected gray” model, and to the third as the “RT” model. Moreover, we investigate the impact of various H₂O and CO₂ outgassing laws. Specifically, we compare the laws used by Elkins-Tanton (2008) and Nikolaou et al. (2019), respectively. Finally, we assess the impact of adding CO₂ to a water-dominated atmosphere for our simulations of the oxidized Earth magma ocean.

In the following, we begin by evaluation the role of the different atmosphere models for pure H₂O and CO₂ atmospheres, respectively. Subsequently, we compare the outcome of simulations using the different atmosphere models for mixed H₂O-CO₂ atmospheres. For these mixed scenarios, we adopt as in Barth et al. (2021), an initial magma ocean depth of 2000 km, consistent with Elkins-Tanton (2008), one of the first models to investigate magma oceans with outgassing of both, H₂O and CO₂. Next, we assess the feedback between H₂O and CO₂ outgassing and compare to the results of Bower et al. (2019).

C.1 Pure H2O atmosphere

For the pure H₂O atmosphere, we follow Barth et al. (2021) by simulating the evolution of the Earth magma ocean with an initial water content of 5 TO for comparison with Hamano et al. (2013). We use the H₂O solubility law of Schaefer et al. (2016) and begin with fully molten mantle down to the core (R_Core = 3400 km), which gives an initial magma ocean depth of 2978 km. We test MagmOc2.0 with the gray atmosphere model and ĸ₀(H₂O) = 0.01 m²/kg, the corrected gray atmosphere model (Appendix B and ĸ₀(H₂O) = 0.25 m²/kg), and the RT atmosphere model (Sect. 2.5).

Figure C.1 (top) illustrates that the gray atmosphere model yields solidification after 0.9 Myrs, as was expected⁹ from Barth et al. (2021). The solidification time of 0.9 Myrs is notably shorter than the Earth magma ocean duration of 4 Myrs reported by Hamano et al. (2013) for the same initial water mass. The simulations with the corrected gray and full RT atmosphere model yield longer solidification times of 1.5 Myrs, aligning more closely to Hamano et al. (2013). The magma ocean simulations reach across all three atmosphere models the runaway greenhouse radiation limit, OLR_lim(1) = 282 W/m², toward the end of the magma ocean evolution. As the simulations progress further, the net flux increases because the Sun’s bolometric flux, and thus the ASR, has substantially decreased after 1 Myrs, whereas the planet’s OLR does not change significantly. We note that the gray atmosphere model consistently assumes OLR_lim(1) = 282 W/m² for T_surf ≤ 1800 K, irrespective of surface pressure and prior OLR evolution. This assumption can result in an abrupt drop in net flux, as is evident in Fig. C.1 (bottom panel).

We further compare our results to the Earth magma ocean simulations of Nikolaou et al. (2019), who also investigated differences between a gray atmosphere model and a full RT model for the pure H₂O atmosphere. However, these authors adopted an initial water content of about 1.0 TO H₂O. For better comparison with this work and that of Elkins-Tanton (2008), we adopt an initial H₂O melt fraction of 0.05 in a 2000 km deep magma ocean, which is equivalent to 1.02 TO initial H₂O mass and aligns with the magma ocean depth of the model by Elkins-Tanton (2008). These simulations yield shorter solidification times ranging between 0.8 and 1 Myrs compared to the 5 TO H₂O simulation, where once again the gray model results again in the shortest solidification time (Fig. C.1, bottom).

In both the 5 TO and 1.02 TO H₂O cases, we consistently find that the OLR of the gray model is higher during the earliest evolution stages (within the range of 10³ – 10⁴ years) compared to the corrected gray and RT atmosphere model. The higher ini- tal OLR of the gray atmosphere model leads to faster overall cooling, thereby explaining the shorter solidification times. A similar deviation in OLR evolution between a gray and RT atmosphere model was reported by Nikolaou et al. (2019, Fig. 3). Consequently, they found that in their simulation with a gray atmosphere model that the magma ocean solidified after several 0.1 Myrs, whereas the simulation using the line-by-line RT model of Katyal et al. (2019) yielded a solidification time of 1 Myrs, in good agreement with our 1.02 TO H2O simulations. These results confirm that the magma ocean lifetime for a pure H₂O atmosphere can differ by several 0.1 million years across different models and atmosphere assumptions (Nikolaou et al. 2019).

However, our simulations show here that the corrected gray and full RT yield evolution tracks that are in good agreement with each other in our magma ocean model. We further find agreement between the OLR of our RT model and the RT calculations of Lichtenberg et al. (2021) for a pure H₂O atmosphere with T_surf = 500 - 3000 K and p_surf = 1 bar and 260 bar. These authors similarly reproduce a magma ocean lifetime of about 1 Myrs with a pure H₂O atmosphere, where we note different definitions for solidification times (see Sect. 3). Furthermore, Lichtenberg et al. (2021) have a cooler initial temperature of 3000 K compared to our initial temperature of 4000 K (Hamano et al. 2013; Schaefer et al. 2016; Nikolaou et al. 2019; Barth et al. 2021). The general agreement between our pure H₂O simulations and that of Lichtenberg et al. (2021) is not surprising, given the usage of very similar H₂O opacities (Table 3) and vertically extended atmospheres for the RT calculations.

C.2 Pure CO2 atmosphere

We simulate the pure CO₂ outgassing scenario for Earth following Elkins-Tanton (2008). Consequently, we assume an initial magma ocean depth of 2000 km and a CO₂ outgassing law with initial saturation in the melt (Sect. 2.2). The simulation assumes an initial CO₂ melt fraction of 0.6, which corresponds to 14.7 TO CO₂.

Figure C.2 demonstrates the thermal evolution of our Earth magma ocean simulations with a dense CO₂ atmosphere. We find that solidification occurs after 0.5 Myrs when the gray model is used, in qualitative agreement with the solidification time of 0.8 Myrs reported by Elkins-Tanton (2008). Table C.1 also illustrates general agreement between the final volatile budget and atmospheric pressure in MagmOc2.0 compared to the results of Elkins-Tanton (2008). The agreement between both models in volatile reservoirs and solidification times with the gray atmosphere model, provides once again confidence that our model, despite its simplification, captures the magma ocean accurately enough for comparison with more complex magma ocean models.

However, significant differences in solidification time compared to the gray atmosphere model arise in simulations that use the corrected gray and RT atmosphere model. More precisely, the solidification time decreases by one order of magnitude to approximately 20 000 - 30 000 years. This discrepancy stems from the fact that the nominal gray model neglects the vertical extent of the atmosphere. The large extension of a hot, dense CO₂ atmosphere is evident from Fig. B.2 and also supported by Lichtenberg et al. (2021, Fig. 3). A pure CO₂ atmosphere has a much steeper temperature gradient in the troposphere compared to H₂O, resulting in hotter atmosphere layers that contribute to the thermal emission on top of the atmosphere. Therefore, a particularly high thermal flux occurs on top of a CO₂-dominated atmosphere for the initial magma ocean stage with high surface temperatures (T_surf ≥ 2000 K) and pressures (p_surf > 1 bar), which are not captured with a gray model (Appendix B).

Figure C.3 provides a detailed illustration of the strong deviation in emission between the gray model and the full RT grid. The RT grid model yields OLR > 10⁵ W/m² for T_surf > 2500 K, whereas the gray model yields for the same temperature and pressure range OLR < 10⁵ W/m². Figure C.2 also demonstrates how the magma ocean simulation with 14.7 TO initial CO₂ mass evolves with the RT grid model (black line) evolves within the emission grid, build with full RT calculations. The comparison reveals that the OLR differs by more than one order of magnitude already at the start of the simulation (Fig. C.2, panel right).

Similarly, Lichtenberg et al. (2021) demonstrate that a magma ocean with a pure CO₂ atmosphere and full RT yields emits large amounts of flux (OLR > 10⁵ W/m²) for T_surf > 2500 K and p_surf = 260 bar. When we further compare our CO₂ magma ocean simulation to that of Lichtenberg et al. (2021), we likewise find that their magma ocean with a 200 - 300 bar dense CO₂ atmosphere reaches surface temperatures between 1500-1600 K, which corresponds to our solidification criterion, after 10 000 years in the same order of magnitude than our CO₂ magma ocean simulation.

While the gray atmosphere model fails to accurately capture the thermal evolution of a magma ocean with a dense CO₂ atmosphere, the agreement between the RT grid simulation and the simulation with the corrected gray atmosphere model is much better (Fig. C.3). The latter yields a solidification time of the same order of magnitude (40 000 years). The results of the corrected gray model shows that our analytical approximation generally reproduce thermal emission even for a hot, vertically extended CO₂ atmosphere (Appendix B).

C.3 Mixed CO2-H2O atmospheres - Different atmosphere models

To validate MagmOc2.0 for a mixed H₂O-CO₂ atmosphere, we follow again Barth et al. (2021) and compare our simulations first to results of the Earth simulations of Elkins-Tanton (2008).

In Barth et al. (2021), the pure H₂O atmosphere already showed promising agreement with Elkins-Tanton (2008).

We test here two scenarios: one with melt fractions 0.05 and 0.01 for H₂O and CO₂, respectively, which we call henceforth the dry scenario, the other with melt fractions 0.5 and 0.1 for H₂O and CO₂, respectively, which we call henceforth the wet scenario. With an initial magma ocean depth of 2000 km, these melt fractions correspond for the dry case to an initial water mass 1.02 TO of H₂O and for the wet case to an initial mass of 10.2 TO H₂O. The initial CO₂ mass is 0.21 and 2.12 TO, respectively. These melt fractions (Table C.1) also correspond to a 5:1 ratio between H₂O and CO₂ that is assumed by Elkins-Tanton (2008). We further adopt the same solubility laws used by Elkins-Tanton (2008) that differ from the H₂O outgassing law used in the previous section, where instead the law by Schaefer et al. (2016) was used (see Sect. 2.2 for an overview of outgassing laws used in this work).

First, we investigate the impact of the gray atmosphere and the full RT atmosphere model for mixed H₂O-CO₂ outgassing. We emphasize, however, that we “switched off” the runaway greenhouse radiation limit in the gray atmosphere model for these specific simulations. This choice is made because Elkins-Tanton (2008) did not include this limiting factor for the OLR. In all other gray simulations in this work, we consider the radiation limit. Figure C.4 (solid lines) and Table C.1 illustrate that our gray simulations without the runaway greenhouse radiation limit are in general agreement with the solidification times reported by Elkins-Tanton (2008): The magma ocean in the wet case solidifies after 1.7 Myrs, and in the dry case after only 0.1 Myrs.

In the wet case, the differences between simulations with the gray atmosphere model and the RT model do not significantly impact the magma ocean solidification times. However, the RT model simulation yields a hotter surface during the early evolution due to a lower initial net flux, resulting in less outgassing of volatiles initially. The hotter magma ocean persists until solidification at 1.7 Myrs, after which the gray and RT evolution tracks converge. This convergence occurs because the net flux in the RT atmosphere simulation is larger compared to the gray simulation in the end, thereby offsetting initial differences.

In the dry case, the magma ocean begins its evolution essentially without a significant atmosphere, when the outgassing laws of Elkins-Tanton (2008) are used. This dry, hot initial state once again leads to a hotter magma ocean in the simulation with the RT atmosphere model compared to the gray simulation. Unlike in the wet case, the evolution tracks do not converge at the end. This is because the net flux of the RT simulation mostly remains below that of the gray simulation until solidification. Consequently, the magma ocean stage is prolonged with the RT atmosphere model compared to simulations using the gray model.

While the lifetime of the Earth magma ocean lifetime may vary with different atmosphere models, this variability does not appear to affect the final volatile content. Simulations using the RT and gray atmosphere model show virtually identical quantities of volatiles in the end (Table C.1).

We further compare the RT grid simulations for H₂O-CO₂ atmospheres with simulations using the corrected gray atmosphere model. Our analysis reveals that in the mixed H₂O-CO₂ setup, the corrected gray atmosphere model underestimates the OLR in the initial evolution stage, during which surface temperatures are very high and surface pressures are particularly low ( < 10 bar). This tendency is illustrated in Fig. B.4 (left): The corrected gray atmosphere model consistently underestimates the OLR compared to the RT grid for high surface temperatures (T_surf ≥ 2000 K) and low pressures (p_surf ≤ 26 bar).

However, these discrepancies in flux are limited to the first thousand years of simulation. More importantly, they do not lead to significant differences in later stages of the magma ocean evolution. Throughout the majority of the magma ocean evolution, in both the dry and wet cases, the surface temperature and atmosphere surface pressures agree well between the RT and corrected gray simulation.

Finally, we present RT atmosphere simulations for the dry and wet cases that exclude CO₂ (Fig. C.5 dash-dotted line). We find that the absence of CO₂ has a minimal influence on the magma ocean evolution, when the outgassing laws of Elkins-Tanton (2008) are used. However, this setup severely suppresses outgassing of CO₂, such that it remains a minor constituent throughout the entire magma ocean evolution. In the following, we thus tackle the impact of alternative solubility laws.

In summary, we generally demonstrate with MagmOc2.O agreement in H₂O and CO₂ distribution with Elkins-Tanton (2008), and also in solidification times, when the gray atmosphere model without runaway greenhouse limit is used. Different atmosphere models may yield variations in solidification time, confirming Nikolaou et al. (2019). Generally, we find that the differences between atmosphere models are largest during the initial stage and we confirm again that the runaway greenhouse radiation limit has a strong impact on the later magma oecean evolution and thus the solidification times.

C.4 Mixed CO2-H2O atmospheres - Different outgassing laws

More recent solubility laws (e.g., Nikolaou et al. 2019; Lichtenberg et al. 2021) suggest that CO₂ should dominate the atmosphere during the initial magma ocean stage due its significantly lower solubility in the magma compared to H₂O (Sect. 2.2). To test the impact of such a scenario, we implement in MagmOc2.0 the outgassing laws of Nikolaou et al. (2019). In the following, we denote simulations using the outgassing laws of Elkins-Tanton (2008) as “ET” and simulations using outgassing laws of Nikolaou et al. (2019) as “Ni.”

To facilitate a comparison of the impact of different outgassing laws with results from previous sections, we once again perform dry and wet case simulations (Table C.2), maintaining an initial magma ocean depth of 2000 km. In the dry-case simulation, volatile inventories are broadly consistent with the setup of the simulations outlined in Nikolaou et al. (2019).

Figure C.6 illustrates the impact of different outgassing laws on the dry case, where major differences are evident for CO₂ surface pressures. In the ET-simulation, CO₂ is outgassed relatively late and never becomes a dominant component of the atmosphere, consistent with results from the previous section. Conversely, in the Niko-simulations, a significant portion of CO₂ is already outgassed at the onset of the magma ocean stage, in agreement with magma ocean simulations of Nikolaou et al. (2019) and Lichtenberg et al. (2021). Consequently, CO₂ initially dominates the atmosphere, with substantial H₂O outgassing occurring only after 10 000 years of simulation time. Eventually H₂O becomes the dominant atmospheric species.

Because the magma ocean evolution in the Ni-simulation begins in the dry case with a thick atmosphere rather than a “bare rock scenario,” the magma ocean cooling is less efficient (Fig. C.6 dash-dotted red line). Consequently, the Ni-simulation reaches solidification later, at 0.8 Myrs, compared to the solidification time of 0.5 Myrs in the ET-simulation (Fig. C.6 solid red line).

We conduct a similar comparison between different outgassing laws and atmosphere models for a volatile-rich wet-case Earth scenario. Here, all initial masses are scaled from the dry case by a factor of ten (Fig. C.7, Table C.2). In this scenario, the thermal evolution of the magma ocean across diverse scenarios is quite similar, converging at around 2 Myrs. Significant differences are again obtained for the H₂O and CO₂ partial pressures. In the ET-simulations, CO₂ consistently remains a minor constituent in the atmosphere. Conversely, in the Ni-simulations, the atmospheric composition initially CO₂-dominated and eventually becomes H₂O-dominated.

Evaluating again the performance of the corrected gray atmosphere model against the RT model in the Ni-simulations, we find that the former overestimates the OLR of a hot, CO₂- dominated atmosphere by a factor of two to three (Figs. C.6 and C.7, dash-dotted purple line). This discrepancy has been already previously noted (Fig. B.4, right panel and Fig. C.2), where it was found that the corrected gray model can deviate from the RT model by up to half an order of magnitude for thick CO₂ atmospheres with high surface temperatures (T_surf > 2500 K). In any case, also here we find that the surface temperature evolution tracks with the corrected gray atmosphere model eventually converge with those with the RT atmosphere. In the dry case, convergence occurs after 10,000 years of simulation time, while in the wet case, convergence occurs after 0.1 Myrs. Despite these discrepancies in surface temperatures, the outgassed H₂O and CO₂ volatile content generally agrees throughout the entire magma ocean evolution between the Ni-simulations with the corrected gray and RT atmosphere model.

Correctly capturing the amount of outgassed volatiles is crucial for assessing the impact of atmospheric erosion on magma oceans around M dwarf stars. For the Earth simulations, the corrected gray atmosphere model, benchmarked against the numerically more costly RT atmosphere model (Sect. D), appears to adequately capture the volatile evolution in Earth magma oceans. This good performance also holds true with outgassing laws that result in atmospheric composition changing from CO₂-dominated to H₂O-dominated.

C.5 Mixed CO2-H2O atmospheres - Impact of CO2

In the ET-simulations, CO₂ consistently remains a minor atmospheric constituent with a negligible influence on the magma ocean evolution. However, in simulations using the outgassing laws used of Nikolaou et al. (2019), CO₂ is initially the dominant atmospheric species that can affect the initial magma ocean stage, as is evidenced by the prolongation of the magma ocean stage in the Ni-simulations compared to the ET-simulations in the dry case. To assess the impact of CO₂ in the Ni-simulations more coherently, we conducted several simulations without CO₂ (Figs. C.6 and C.7, dotted red curves).

A detailed inspection of the thermal evolution of the dry case shows that a mixed CO₂ -H₂ O atmosphere exhibits higher thermal emission during the first 0.3 Myrs of the evolution compared to a pure H₂O atmosphere. This difference arises because a CO₂-rich atmosphere emits more flux than an equivalent H₂O- atmosophere for high surface temperatures (T_surf > 2000 K, see Sections C.2 and B). However, beyond 0.3 Myrs, the system reaches the runaway greenhouse limit that is set primarly determined by the water content in the atmosphere.

As is discussed in Sect. 2.4, the addition of CO₂ tends to cool the emitting atmosphere layers, thereby reducing the OLR limit and diminishing thermal cooling in a mixed atmosphere compared to a pure H₂O atmosphere. In other words, the addition of CO₂ for a sufficiently oxidized mantle (Ortenzi et al. 2020) tends to shorten the magma ocean lifetime outside of the runaway greenhouse limit and tends to extend the solidification time within it. In the dry case, both effects counterbalance each other, resulting in a slightly longer solidification time for the mixed atmosphere (by a few tens of thousands of years or a few percent).

The higher thermal emission of a mixed atmosphere outside of the runaway greenhouse limit is more pronounced when more volatiles are in the system as in the wet case (Fig. C.7, dotted red line). However, even in this scenario, the runaway greenhouse limit reached toward the end of the evolution results after 2 Myrs in a convergence of solidification time compared to a magma ocean model with a pure H₂O atmosphere. Upon close examination, the mixed atmosphere exhibits a slightly shorter solidification time, albeit by only a few percent.

Our assessment of a minimal impact of CO₂ on the solidification times of magma oceans disagrees with the results of Nikolaou et al. (2019) who report an extension of magma ocean lifetime with additional CO₂. These authors used, however, a gray atmosphere model. As is discussed in Appendix B, a gray atmosphere model consistently overestimates thermal emission for thick H₂O atmospheres (Fig. B.1) and underestimates thermal emission for thick CO₂ atmospheres (Fig. B.3) because it does not adequately capture the thermal emission in a vertically extended atmosphere, as already outlined by Lichtenberg et al. (2021).

C.6 Mixed CO2-H2O atmospheres - Importance of initial conditions

We compare, similarly to Table C.1, the volatile budget at the end of the magma ocean between the ET and Ni-simulations (Table C.3). The most notable difference is the very low remaining CO₂ content in the latter cases. This discrepancy can be explained by the different initial conditions, driven by the different CO₂ outgassing laws.

In the model of Elkins-Tanton (2008), CO₂ is mostly dissolved in the melt initially, resulting in a relatively large initial melt fraction, FCO , compared to the outgassing laws used by Nikolaou et al. (2019). As the magma ocean solidifies and the thickness of the magma ocean decreases, as a result significant enrichment of volatile mass fraction occurs, assuming no substantial sink terms. Only a very small fraction of CO₂ is partitioned from the melt into the solid mantle (Table 1) and atmospheric erosion for Earth is negligible during the magma ocean evolution. Thus, the volatile mass fraction at the end of the magma ocean stage can only be larger than the initial value.

Because CO₂ outgassing is suppressed in the ET-simulations, initial values of $F_{{\text{H}}_2\text{O}}$ are relatively high already. Consequently, the magma ocean ends with large CO₂ mass fractions of 0.2 - 0.7 wt%. Conversely, the Ni-simulations begin with a substantially lower mass fraction in the magma ocean and consequently end the magma ocean stage with CO₂ mass fractions of 0.003 - 0.3 wt% , which is up to two orders of magnitudes lower compared to the ET-simulations. The differences are much smaller for H₂O, because the ET and Ni-outgassing laws assume both a high solubility of H₂O, resulting in similar initial values of $F_{{\text{H}}_2\text{O}}$.

We thus conclude that the initial conditions, which determine the initial volatile mass fraction in the magma ocean, are critical for determining the overall volatile budget at the end of the magma ocean phase and how much of the volatiles can be retained in the mantle.

C.7 Mixed CO2-H2O atmospheres - Outgassing feedback

We confirm here that the evolution of the atmospheric CO₂ and H₂O exhibit a feedback effect, when the mean molecular weight of the atmosphere changes during the magma ocean evolution (Bower et al. 2019). A change in atmospheric composition occurs during the magma ocean evolution because CO₂ is less soluble than H₂O and thus dominates outgassing at the beginning of the magma ocean stage, as captured by the outgassing laws of Nikolaou et al. (2019) (Sect. C.4). The model of Elkins-Tanton (2008) avoids this effect by suppressing CO₂ outgassing and hence these simulations show no such feedback.

Figure C.8 provides more detailed insights about the outgassing feedback in the Nicolaou (short: Niko) simulations. In the mixed volatile setup (solid lines), the atmospheric composition undergoes drastic changes from CO₂-dominated (molecular weight of CO₂: µ = 44 g/mol) to H₂O-dominated (molecular weight of H₂O: µ = 18 g/mol), as also reported by Nikolaou et al. (2019); Bower et al. (2019). Consequently, the mean molecular weight of the atmosphere decreases as the magma ocean solidifies. Additionally, the melt fraction of CO₂ decreases as H₂O begins to accumulate in the atmosphere (Fig. C.8, left panel, dotted red lines). In the absence of H₂O, CO₂ outgassing would continue to increase steadily as the magma ocean solidifies. In contrast, the melt fraction of H₂O, and thus the partial pressure of H₂O, increases in the mixed atmosphere case. In the absence of CO₂ , the melt fraction still increases as the mantle solidifies, albeit not as strongly (Fig. C.8, left panel, dotted blue line).

The feedback between atmosphere composition and outgassing affects the volatile distribution at the end of the magma ocean. For the dry case, the final CO₂ partial pressure decreases from $p_{{\text{CO}}_2}=100$ bar in a pure CO₂ atmosphere to $p_{{\text{CO}}_2}=67$ bar in a mixed atmosphere as H₂O becomes the dominant species. The interplay between atmospheric composition and outgassing results in higher H₂O outgassing ($p_{{\text{H}}_2\text{O}}=206\text{bar}$) compared to the scenario with no CO₂ in the system, where ’only’ $p_{{\text{H}}_2\text{O}}=206$ bar is outgassed. The changes in H₂O and CO₂ pressures for the dry case are in first order agreement with Bower et al. (see 2019, their Fig. 7a) for a similar volatile content. Nikolaou et al. (2019) did not account for the outgassing feedback effect and thus report an increase in both, CO₂ and H₂O partial pressures as the magma ocean solidifies.

We further find, for the wet-case magma ocean, that volatile outgassing is modified similarly to the dry case. We also emphasize again that a drastic change in atmosphere composition, and thus in outgassing, does not occur when the outgassing laws of Elkins-Tanton (2008) are used, where H₂O is always the dominant volatile in the atmosphere.

We conclude that MagmOc2.0 adequately accounts for outgassing feedback between H₂O and CO₂, and that the corrected gray and RT atmosphere model capture thermal emission in vertically extended thick atmosphere. Even more, MagmOc2.0 is well-equipped to also tackle planets in the habitable zone around M dwarf stars, where atmospheric erosion of H₂O mainly by XUV photolysis (see e.g., Luger & Barnes 2015) has the potential to further change the composition of a mixed H₂O-CO₂ atmosphere during the magma ocean stage.

Appendix D Numerical stability and runtime

The code MagmOc2.0 is designed to efficiently explore diverse planetary and stellar scenarios. To achieve this goal, the code solves various sets of ordinary differential equations (ODEs), including equations that describe outgassing (Sect. 2.2). In this appendix, we provide details on numerical stability, convergence tests, and the typical wall clock times of the simulations discussed earlier. We focus here on the albedo=0.75 simulations.

One parameter in these specific ODEs is the decrease in liquid magma ocean $\frac{d{M}^{liq}}{dt}$ as the solidification radius increases $\frac{d{r}_s}{dt}$, which undergoes, however, phase state transitions for the critical mantle melt fraction ψ_c = 0.4 (Barth et al. 2021, Sections 2.1.1 and 2.1.2). To tackle numerical instabilities introduced by non-continuous changes of melt fractions $F_{{\text{CO}}_2}(t)$ and $F_{{\text{H}}_2\text{O}}(t)$ in time t, we carefully monitor the volatile masses in the different reservoirs of the coupled MOA system ($M_i^{\text{cystal\hspace{0.17em}}}+M_i^{liq}+M_i^{\text{atm\hspace{0.17em}}}$) to ensure that their sum is equal to $M_i^{\text{moa}}$ within reasonable limits (differences of less than 5%). $M_i^{\text{moa}}$ also has to be initially equal to the prescribed initial volatile mass, $M_i^{\text{ini}}$. To avoid numerical “leakage of mass,” we added two terms, $\mathrm{\Delta }{M}_{{\text{H}}_2\text{O}}^{\text{Corr\hspace{0.17em}}}$ and $\mathrm{\Delta }{M}_{{\text{CO}}_2}^{\text{Corr\hspace{0.17em}}}$, that provide source or sink terms to the partial derivative of $M_{{\text{H}}_2\text{O}}^{\text{moa}}$ to ensure mass conservation. These two terms take the form $$\begin{array}{l}\mathrm{\Delta }{M}_{{\text{H}}_2\text{O}}^{\text{Corr\hspace{0.17em}}}=\frac{M_{{\text{H}}_2\text{O}}^{\text{ini}}-M_{{\text{H}}_2\text{O}}^{\text{moa\hspace{0.17em}}}}{\mathrm{\Delta }{t}_{\text{curr\hspace{0.17em}}}}\\ \mathrm{\Delta }{M}_{{\text{CO}}_2}^{\text{Corr\hspace{0.17em}}}=\frac{M_{{\text{CO}}_2}^{\text{ini}}-M_{{\text{CO}}_2}^{\text{moa\hspace{0.17em}}}}{\mathrm{\Delta }{t}_{\text{curr\hspace{0.17em}}}}\end{array}$$(D.1)

where ∆t_curr is the current time step during the magma ocean evolution that ensures that the volatile mass budget during runtime does not deviate more than 5% from the initial mass budget. We note that $\mathrm{\Delta }{M}_{{\text{H}}_2\text{O}}^{\text{Corr\hspace{0.17em}}}$ and $\mathrm{\Delta }{M}_{{\text{CO}}_2}$ are rates in units of [kg/s] to counterbalance numerical mass loss.

That is, the differentiation of the mass balance equations for volatile i is modified during simulation run time such that: $$\frac{d{M}_i^{\text{moa\hspace{0.17em}}}}{dt}+\mathrm{\Delta }{M}_i^{\text{Corr\hspace{0.17em}}}=\frac{d{M}_i^{\text{crystal\hspace{0.17em}}}}{dt}+\frac{d{M}_i^{\text{liq}}}{dt}+\frac{d{M}_i^{\text{atm}}}{dt}.$$(D.2)

In terms of the substitution framework (Table A.1), we add $\mathrm{\Delta }{M}_i^{\text{Corr\hspace{0.17em}}}$ to $a_1^{\prime }$ (for H₂O) and $a_1^{\prime }$ (for CO₂), respectively. To avoid overcompensation, we limit the correction term to values smaller than 1% of the total volatile budget. Despite these measures, sometimes noticeable fluctuations in the outgassed volatiles can occur due to fluctuations in F_i.

We demonstrate such fluctuations for the TRAPPIST-1 g magma ocean evolution tracks for 1 TO, 5 TO, and 100 TO initial mass H₂O with CO₂ initial mass scaled by H₂O mass with 0x, 0.3x and 1x the inital H₂O mass simulated with different relative accuracies ϵ of the Runge-Kutta integrator (Figs. D.1, D.2, D.3).

It is evident that atmospheric pressures and net flux can strongly fluctuate for the lowest integrator accuracy of ϵ = 10^–1. In all investigated cases, however, these fluctuations have a negligible impact on the surface temperature and overall volatile evolution. This is confirmed by comparison with simulations of higher numerical accuracy. More precisely, a relative accuracy of ϵ = 10^–3 in the Runge-Kutta integrator results in stable simulations that can be used to benchmark a specific evolution track in case of numerical issues.

We further show the run times of the 18 TRAPPIST-1 g simulations with 1 TO, 5 TO, and 100 TO initial mass of H₂O and an equivalent quantity of CO₂ mass with different numerical accuracies (Table D.1). The simulations were carried out on a single CPU on a AMD Ryzen Threadripper PRO 5955WX with 16 cores. The magma ocean simulations using the full RT atmosphere model may require up to 28 mins to complete. Simulations with the corrected gray atmosphere model, benchmarked with full RT calculations, are much more efficient, with runtimes of 1 minute and less, depending on numerical accuracy. The difference in runtime between ϵ = 10^–3 and 10^–1 is just 15 to 20 seconds. Thus, with the corrected gray atmosphere model, the simulations can be carried out with highest numerical stability without excessive numerical costs. MagmOc2.0 as part of the VPLanet-software framework is thus ideally suited to explore the magma ocean evolution for rocky exoplanets in the habitable zone of their host stars, spanning a wide range in fundamental parameters and exploring different processes to account for the large diversities in exoplanet evolution.

Appendix E Chemical composition

Since stars and their planetary accretion disks are formed by the collapse of the same interstellar dust cloud, the composition of a star can be used as a first estimate for the upper limit of the composition of the accretion disk. However, the stellar composition of TRAPPIST-1 has not yet been determined. We therefore derive elemental abundances for the main planet-forming elements from the stellar metallicity. For this, large-scale astronomical surveys play a crucial role in providing the necessary data for understanding the chemical composition of stars. One such survey is the Galactic Archaeology with HERMES (GALAH) project. The GALAH survey focuses on determining the detailed chemical composition of Milky Way stars, contributing to a better understanding of Galactic chemical evolution. In this study, the third release of the GALAH survey by Buder et al. (2021) serves as the main resource, offering a rich dataset that facilitates the calculation of metallicity-dependent compositional variations. However, the GALAH survey has limitations in providing data for all chemical elements of interest to our study. To address this gap, the Hypatia catalog was incorporated specifically for elements such as nitrogen and sulfur. The Hypatia Catalog, compiled from 84 literature sources, presents spectroscopic abundance data for 50 elements across 3058 stars in the solar neighborhood (within 150 pc of the Sun). Employing a binning strategy with intervals set at [Fe/H] = 0.05, stars were systematically organized to facilitate thorough data analysis. The mean metallicity and chemical composition for each bin was computed, providing a representative value within that specific metallicity range. For TRAPPIST-1, we use a metallic- ity value of 0.04±0.08 (Gillon et al. 2017), which leads to the predicted stellar composition reported in Sect. 3.3.

Bitsch & Battistini (2020) suggested a stoichiometric model to obtain a first-order estimate on the compositional variation in planetary building material depending on the local temperature within an accretion disk. In this approach, the gas within the accretion disk is assumed to have achieved the state of chemical equilibrium before condensation, with the complete set of molecules preexisting in the gas. Consequently, the relative abundance of molecules can be calculated stochiometrically and based on their condensation temperature Lodders (2003).

For a first estimate on the temperature profile within the accretion disk, we employ a simple power law Williams & Cieza (2011): $$T_{disk}(r)=T_S{\left(\frac{r}{x{R}_i}\right)}^{-3/4},$$(E.1)

where T_S represents the sublimation temperature and is set to 1500 K at the inner edge of the disk Ri (following observations by Dullemond & Monnier 2010). x is a free scaling parameter which we set here to x = 2 for a good agreement of our profile compared to Jorge et al. (2022). Using this estimate for the disk temperature profile allows us to determine the composition of the planetary building blocks at different distances of the star. For this, we extended the model of Bitsch & Battistini (2020) to include nitrogen, aluminum and calcium species. The full set of molecules, their condensation temperatures and stoichiometric calculations are listed in table E.1. Due to the overabundance of hydrogen, the elemental abundances are expressed as X/H.

Appendix F Interior structure models

We integrate the compositional information in our interiorstructure models for all three planets. It is important to note here that while we directly employed the refractory elements condensing out of the accretion disk in our interior-structure model, we allowed for the core-mass fraction (considering here for simplicity a pure iron core) as well as the volatile fraction (considering here for simplicity only water, since it is the most abundant volatile molecule) to vary, to explore the range of potential interior structures in TRAPPIST-1 e, f, and g. For a first justification of our model, we applied our interior structure model without any volatile material to the three innermost planets, and obtain modeled planet radii consistent with observations (see Sect. 3.3). Figure F.1 displays the full range of water content within error bars of planet radii and mass measurements.

The interior structure model follows Noack et al. (2016), where we dissect the planet into 1000 shells, which we fill with core, mantle or ice/water material for given planet mass as well as core and water mass fractions. We assume an adiabatic temperature profile from surface to the center of the planet, while adopting a temperature increase at the core-mantle boundary following Stixrude (2014) and Noack & Lasbleis (2020). The surface temperature is set at 300 K, which means that we consider here for simplicity at the surface a liquid water layer. High- pressure ice forms for all three planets for water fractions starting at a few wt-%.

For the core and water thermodynamic properties (such as density or heat capacity), we employ equations of state as outlined in Noack et al. (2016). For the silicate mantle, we improved the approach of that study by developing lookup-tables for all relevant thermodynamic properties for silicate mantles of varying composition calculated with Perple_X (Connolly 2009) and interpolated for our predicted mantle composition as well as local temperature and pressure conditions in each shell. We use a different resolution for our look-up tables for low, intermediate, and high pressures, where we vary temperature in steps of 25 K between 200 and 5000 K and pressure in steps of 1250 bar for pressures below 25 GPa, 1.9 GPa for pressures between 25 and 400 GPa, and 5.5 GPa for pressures above 400 GPa going up to 1500 GPa (which ensures coverage of the entire pressure range encountered in the planetary interior in this study).

References

---

All Tables

All Figures