{'fraction_non_alphanumeric': 0.030463393309097603, 'fraction_numerical': 0.023330725691408587, 'mean_word_length': 4.477660848194593, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Magma ocean evolution of the TRAPPIST-1 planets

Recent observations of the potentially habitable planets TRAPPIST-1 e, f, and g suggest that they possess large water mass fractions of possibly several tens of wt% of water, even though the host star's activity should drive rapid atmospheric escape. These processes can photolyze water, generating free oxygen and possibly desiccating the planet. After the planets formed, their mantles were likely completely molten with volatiles dissolving and exsolving from the melt. In order to understand these planets and prepare for future observations, the magma ocean phase of these worlds must be understood. To simulate these planets, we have combined existing models of stellar evolution, atmospheric escape, tidal heating, radiogenic heating, magma ocean cooling, planetary radiation, and water-oxygen-iron geochemistry. We present MagmOc, a versatile magma ocean evolution model, validated against the rocky Super-Earth GJ 1132b and early Earth. We simulate the coupled magma ocean-atmospheric evolution of TRAPPIST-1 e, f, and g for a range of tidal and radiogenic heating rates, as well as initial water contents between 1 and 100 Earth oceans. We also reanalyze the structures of these planets and find they have water mass fractions of 0-0.23, 0.01-0.21, and 0.11-0.24 for planets e, f, and g, respectively. Our model does not make a strong prediction about the water and oxygen content of the atmosphere of TRAPPIST-1 e at the time of mantle solidification. In contrast, the model predicts that TRAPPIST-1 f and g would have a thick steam atmosphere with a small amount of oxygen at that stage. For all planets that we investigated, we find that only 3-5% of the initial water will be locked in the mantle after the magma ocean solidified.


Introduction
As of March 2021, there are 4692 confirmed exoplanets 1 . Of those exoplanets, 60 are optimistically classified as potentially habitable 2 , meaning they orbit in the star's habitable zone and are rocky. The habitable zone is defined as the region around a star where liquid water can be present on the surface of a rocky planet, assuming a greenhouse atmosphere (Huang 1959;Kasting et al. 1993a;Kopparapu et al. 2013). According to population statistics from the Kepler exoplanetary survey, a planet should likely be less massive than 5 M C (Otegi, J. F. et al. 2020) and smaller than 1.7R C (e.g. Ginzburg et al. 2018) to be rocky. Only 14 of these 60 planets are in the conservative habitable zone and have a mass smaller than 5 M C and are, therefore, currently the most likely candidates for being habitable rocky worlds. Here, we use the definition of the conservative habitable zone by Kopparapu et al. (2014) with the inner edge being the mass-dependent runaway greenhouse limit and the outer edge the maximum greenhouse limit. However, it is important to note that a planet's position in the habitable zone is only one of many factors that determine the habitability of said planet.
Most of these potentially habitable planets are orbiting M dwarfs, including the planets closest to Earth, e.g. Proxima Centauri b (Anglada-Escudé et al. 2016), TRAPPIST-1 e, f, & g , and Teegarden's Star c (Zechmeister, M. et al. 2019), which are prime candidates for observations with upcoming space-and ground-based instruments like the James Webb Space Telescope (JWST) (Barstow & Irwin 2016;Snellen et al. 2017;Lincowski et al. 2018) or the Extremely Large Telescope (ELT) (Snellen et al. 2015;Meadows et al. 2017). To define the goals for these measurements and identify detectable tracers for life, a firm understanding of the evolution of terrestrial planets around M dwarfs is required (Meadows 2017;Catling et al. 2018).
Planets orbiting M dwarfs are in a very different environment than Earth. These stars are active and a large number of flares reach these rocky planets, which orbit their stars on a much tighter orbit than Earth orbits the Sun (Vida et al. 2017;Günther et al. 2020). The intense X-ray and UV irradiation (XUV) emitted by active M dwarfs leads to enhanced atmospheric erosion (Watson et al. 1981;Lammer et al. 2003Lammer et al. , 2009Erkaev et al. 2007;Owen & Jackson 2012) that could harm organisms living on the surface of the planet (Lammer et al. 2007). During the pre-main sequence phase, M dwarfs are much more luminous, with their potentially habitable planets receiving XUV fluxes that are up to 100 times more powerful than on the modern Earth (Lammer et al. 2007).
Studies on the composition of the TRAPPIST-1 planets by Dorn et al. (2018) and Unterborn et al. (2018a) suggest that some of these planets today have water mass fractions of possibly several tens of wt%, but with large uncertainties. These abundances are orders of magnitude higher than the estimated current water mass fraction of Earth (0.07 -0.23 wt% water or three to ten "terrestrial oceans 3 ," e.g. Elkins-Tanton (2008); Pearson et al. (2014); Schmandt et al. (2014)). This discrepancy suggests a very different formation and evolution history for the TRAPPIST-1 planets compared to the rocky Solar System planets, which is still the subject of ongoing research. Therefore, we need new models to simulate the evolution of terrestrial planets around M dwarfs to be able to explain how much water the TRAPPIST-1 planets could have held on to during their evolution and how much water they must have initially received during their formation in the first place.
Due to the heat produced by accretion, core formation, as well as radiogenic and tidal heating, most planetary-sized bodies in the Solar System appear to have been in a fully molten state right after formation, during which the mantle vigorously convected while gases passed between the atmosphere and liquid rock with relative ease (Wood et al. 1970;Solomon 1979;Wetherill 1990;Lammer et al. 2018). This planetary state, called a magma ocean, drives rapid geochemical evolution that dramatically alters the redox state of the mantle and atmosphere. A thick greenhouse atmosphere heated by the surface and stellar irradiation moderates the cooling of the magma ocean. In many cases, an active nearby star releases high energy photons that dissociate water and remove the atmosphere. Hydrogen escapes more easily than oxygen, so the latter will either enter the melt and bond with iron, or accumulate in the atmosphere (e.g. ). This conceptual model of a magma ocean, which is greatly simplified, is shown schematically in Fig. 1. Note that this is assuming an oxidized magma ocean which is outgassing H 2 O rather than H 2 .
Such a scenario was tackled with simulations for rocky Solar System planets initially for an early atmosphere dominated by water vapour by Matsui & Abe (1986); Zahnle et al. (1988); Abe (1993Abe ( , 1997; Solomatov (2000); Zahnle et al. (2015) and by Elkins-Tanton (2008) for an early H 2 O and CO 2 atmosphere. Schaefer et al. (2016) were the first to investigate magma ocean outgassing and atmosphere erosion for a rocky exoplanet around an M dwarf star, the hot, un-inhabitable Super-Earth GJ 1132b. They found that around a pre-main sequence M dwarf, atmosphere erosion can be so effective that abiotic free O 2 builds up in the early atmosphere, as much as several thousands of bars.
Abiotic oxygen build-up on rocky planets around M dwarfs was already proposed by , but without taking into account oxygen solubility back into the magma ocean. Schaefer et al. (2016) considered a pure steam atmosphere and did not take into account the contribution of CO 2 . They also did not take into account that additional heating sources like tidal heating (e.g. Driscoll & Barnes 2015) and magnetic induction (Kislyakova et al. 2017) could potentially delay solidification of the magma ocean and thus reduce abiotic oxygen build-up.
In this study, we consider magma ocean volatile outgassing and atmosphere erosion for the potentially habitable TRAPPIST-1 planets e, f, g considering an early H 2 O dominated atmosphere. We consider initial water contents from 1 -100 terrestrial oceans and compare the final water content to the results obtained by different interior structure models (Noack et al. 2016; Barr et al. 2018;Dorn et al. 2018;Unterborn et al. 2018b). We also investigate how different levels of interior heating, from radiogenic and tidal sources, affect the final atmospheric composition.
To perform our simulations, we have augmented the open source software package VPLanet (Barnes et al. 2020) to include the physics and chemistry necessary to model the TRAPPIST-1 planets during the magma ocean phase. VPLanet already includes "modules" for radiogenic heating (RadHeat), tidal heating (EqTide), stellar evolution (STELLAR) and water photolysis and hydrogen escape (AtmEsc). Here we introduce a new module MagmOc that calculates the thermal evolution of a magma ocean, volatile fluxes through the surface, and radiative cooling through an H 2 O atmosphere.
In Section 2, we present the magma ocean model, which builds from Schaefer et al. (2016) and Elkins-Tanton (2008). We further describe different calculations of the atmospheric net flux that leads to the cooling of the planet. We discuss our model in comparison to previous work by Schaefer et al. (2016) for GJ 1132b and Elkins-Tanton (2008) and Hamano et al. (2013) for Earth in Section 3. In Section 4, we present the results obtained for TRAPPIST-1 e, f, and g. We discuss the ramification of our results for the TRAPPIST-1 planets in Section 5. We summarize our results in Section 6 and conclude with an outlook, where we offer our magma ocean model to the community as the MagmOc module of VPLanet (Section 7).

The Model
In this section, we describe the set-up of the MagmOc module. We stress that we describe here a relatively simple initial set-up to establish an open-source magma ocean model that is useful for application to extrasolar and Solar System rocky planets and for a broad community of researchers. In the near future, this model will be expanded to include more complex descriptions of the magma ocean.  taken from Schaefer et al. (2016). A concise list of the most important parameters in the thermal model is provided in Table 1. Their values, if appropriate, used for our simulations can be found in the appendix in Table 7. If not otherwise noted, these values are adapted from Schaefer et al. (2016).

Thermal Model
We start our simulations always with a completely molten mantle, which solidifies from the base to the top. In addition to the primordial heat from accretion, the magma ocean is heated by the decay of radioactive isotopes and tidal interactions. Tidal heating was not considered by Schaefer et al. (2016) even though it can be a very large source of energy (Barnes et al. 2010(Barnes et al. , 2013Driscoll & Barnes 2015). The magma ocean cools by turbulently convecting internal energy to the surface, where it is ultimately radiated into space. The atmosphere, in our case a thick atmosphere composed of water, is convective in itself. It cools by emitting infrared radiation (IR) into space, while it is heated by stellar radiation (SR). Part of SR is reflected back to space due to the albedo of a water vapour atmosphere. A schematic of such a planet and the processes in our model is shown in Fig. 1 For the age of the star at the beginning of our simulations, we use the estimated lifetime of the protoplanetary disk, since we are not taking into account further impacts during the magma ocean evolution. Ribas et al. (2014) give an estimate for the lifetime of 4.2´5.8 Myr for disks of low-mass stars, based on statistics of the solar neighborhood. Therefore, we use an age of the star of 5 Myr at the beginning of our simulations.
The heating and cooling processes of the planet are controlled by the thermal model which tracks the potential temperature T p of the mantle and the solidification radius r s for every time step. T p is the temperature of the uppermost layer of the magma ocean and defines the temperature profile of the mantle (see Eq. 4). As convection in the molten part of the mantle is very efficient, the magma ocean and surface are in thermal equilibrium (see Sec. 2.1.1). We therefore assume the surface temperature to be T surf " T p . In addition, the thermal model keeps track of the melt fraction ψ and the time derivative of r s which are needed for the volatile model. The initial parameters are T p " 4000 K and r s " r c , where r c is the core radius. This temperature was chosen to ensure that the mantle is completely molten at the start of the simulation. (See Section 2.1.2 for more information on the treatment of the solidification radius.) More precisely, the magma ocean temperature is described by the following equation (adapted from Schaefer et al. 2016): with the first term on the right hand side describing the heat produced by the fusion of silicates (latent heat), the second term the radiogenic and tidal heating and the third term is the atmospheric net flux cooling the planet off to space. This equation is based on the model by Hamano et al. (2013, Eq. 2) but includes additional heat sources on the right hand side that are also partly included in Schaefer et al. (2016). For simplification of the model we assume the additional heat by radiogenic and tidal heating is distributed over the whole mantle and not just the magma ocean. Therefore, we apply this equation on the whole mantle as well, i.e. the average mantle temperature T m . This temperature is related to the potential temperature by T m " m T p with m " 1.19 for planets of approximately one Earth mass (Schaefer & Sasselov 2015). Following previous models (Elkins-Tanton 2008;Hamano et al. 2013;Lebrun et al. 2013;Schaefer et al. 2016;Nikolaou et al. 2019) we do not include the heat flux from the core into the mantle. In the following, we describe the terms on the right hand side of Eq. (1) that describes the thermal evolution of the planet.

Parameterization of convection
The cooling and solidification of the mantle is driven by very strong convection due to the low viscosity of the liquid melt. Only when the surface melt fraction drops below a critical value ψ c does the viscosity increase. The viscosity thus depends on the melt fraction ψ. The melt fraction at the surface is expressed for T solidus ă T p ă T liquidus by ψ surf " pT p´Tsolidus,surf q pT liquidus,surf´Tsolidus,surf q , with T solidus,surf " 1420 K and T liquidus,surf " 2020 K. Lebrun et al. (2013) showed that ψ " 0.4 marks the phase transition towards a thicker, viscous thermal boundary layer and the effective end of the magma ocean stage 5 . For the viscosity profile, therefore, two different regimes of the surface melt fraction are considered: Above the critical melt fraction ψ c " 0.4, the viscosity parameterization from Lebrun et al. (2013) is used, while for ψ surf ă 0.4, the solid-like viscosity needs to be applied and we follow Lebrun et al. (2013) and Schaefer et al. (2016). The dynamic viscosity η is thus defined by: with η l " 0.000 24 Pa s, η s " 3.8ˆ10 9 Pa s, the activation energy E a " 350 kJ mol´1, and the ideal gas constant R. This returns a viscosity of η ă 0.1 Pa s for a fully molten surface and only for a surface melt fraction of just above 0.4 the viscosity increases to η ą 10 3 Pa s. We therefore conclude that our approximation of efficient convection is valid for most of the solidification process 6 . Please note that we did not include the effect of water in the melt on the viscosity of the magma ocean. Studies by Nikolaou et al. (2019) show that the effect of water dependent changes of the viscosity on the solidification time of the mantle is negligible.

Mantle solidification
To calculate changes in the solidification radius r s , we define the adiabatic temperature structure T prq of the convective mantle with respect to the potential temperature T p : T prq " T p r1`α g c p pr p´r qs.
Here, g is the planetary surface gravity, c p is the specific heat capacity of the magma ocean and α is the thermal expansion coefficient.
To parameterize the solidification, we use the solidus and liquidus profiles from Hirschmann (2000): where a " 104.42 K GPa´1 and b " 1420 K for p ă 5.2 GPa, and a " 26.53 K GPa´1 and b " 1825 K for p ą 5.2 GPa. The depth z dependent pressure is given by ppzq " ρ m gpz´z 2 {2r p q. The liquidus temperature is assumed to be larger than the solidus by 600 K. The solidification radius r s is defined as the radius where the mantle temperature T equals T solidus where the melt fraction is ψ " 0. Equations (4) and (5) are used to calculate the change of the solidification radius r s : We caution the reader that r s only changes when the temperature at the core-mantle boundary r c drops below T solidus . When choosing initial conditions such that the whole mantle is molten, i.e. ψpr c q " 1, we set r s " r c until the base of the mantle starts to solidify.

Radiogenic and tidal heating
Radiogenic heat in the Earth is generated by the decay of 238 U, 235 U, 232 Th, and 40 K in the core, mantle, and crust. The radiogenic power produced by species i in reservoir j is where λ i,1{2 " ln 2{τ i,1{2 , τ i,1{2 is the halflife, t is time, and Q i, j p0q is the initial heat production at t " 0. This heat source is included 6 See Sec. 5 for a discussion of the final 2% of mantle solidification. in our model by connecting it to the RadHeat module of VPLanet (Barnes et al. 2020). Although the planets may form quickly enough for 26 Al heating to be important (Barnes et al. 2016;Lichtenberg et al. 2019), we ignore this effect here 7 . Similarly, the tidal heating rate Q t is calculated by the module EqTide of VPLanet. It is based on the equilibrium tide theory (Darwin 1880;Ferraz-Mello et al. 2008;Leconte et al. 2010), assuming the gravitational potential of the tide raiser on an unperturbed spherical surface can be expressed as the sum of surface waves and that the elongated equilibrium shape of the perturbed body is slightly misaligned with respect to the line that connects the two centers of mass (Barnes et al. 2020). The advantage of this model is that it is semi-analytic and reduces the tidal effects to a single parameter (the tidal quaility factor Q), which is especially valuable for exoplanets where little is known about their interior composition. Self-consistent models, however, would need threedimensional models including the rheology of the interior. Ocean worlds would even need three-dimensional models of the ocean currents (Carone 2012). In our simulations we use the "constantphase-lag" (CPL) model in EqTide where, regardless of orbital and rotational frequencies, the angle between the perturber and the tidal bulge remains constant. When the body is assumed to behave like a harmonic oscillator, the damping must be independent of the frequency for the tidal waves to be linearly summed. The lag between the line connecting the two centers of mass and the direction of the tidal bulge is inversely proportional to the "tidal quality factor" Q.
In this study, we only include the tidal effect on the planet by the star. However, in multi-planet systems tidal interactions between individual planets might be important (Wright 2018; Hay & Matsuyama 2019). Further, we did not include tidal effects raised by a moon because close in planets are unlikely to have moons due to tidal effects between the star, planet, and moon (Barnes & O'Brien 2002;Sasaki & Barnes 2014).

Atmospheric Flux
The simplest way to calculate the outgoing long-wave radiation is by implementing a grey atmosphere model. We base this part of our model on the the grey atmosphere model of Elkins-Tanton (2008), which was applied to Earth, and on the model of Carone et al. (2014), developed for different atmosphere compositions on rocky exoplanets. In the grey radiative transfer approach, the opacity is not calculated for every wavelength bin. Instead, the average infrared opacities of the individual species are combined, weighted by the partial pressure of the specific gas.
The atmospheric net flux is then calculated by multiplying the difference between the Stefan-Boltzmann emission from the surface (σT 4 surf ) and the top of the atmosphere (σT 4 eq ) with the emissivity : F " σpT 4 surf´T 4 eq q.
T eq is the equilibrium temperature of the atmosphere due to the stellar irradiation. The emissivity is calculated with the optical depth τ˚which is the sum of the optical depths of all atmospheric species. The emissivity is given by 7 For the purpose of this study we assume that the start point of our simulation is late enough for 26 Al heating to be negligible with the optical depth of the atmospheric species i M atm i is the atmospheric mass of species i. To calculate the optical depth, we use the atmospheric absorption coefficient k 0,H 2 O " 0.01 m 2 kg´1 for water (Elkins-Tanton 2008) at the reference pressure p 0 " 101 325 Pa (Yamamoto et al. 1952).
One aspect of atmospheric behavior that is not included in the grey atmosphere model is the runaway greenhouse or blanketing effect. This effect prevents the planet from emitting more than 280 W m´2 when greenhouse gasses in the atmosphere are present, unless the surface temperature is larger than 1800 K (Goldblatt et al. 2013;Kopparapu et al. 2013). Once the surface temperature of a rocky planet with a greenhouse atmosphere enters the runaway greenhouse range (600´1800 K), the planet will ultimately heat up because it cannot emit enough radiation to cool. For surface temperatures of T surf 1600 K the temperature of the atmospheric layer at which the optical depth is unity (τ " 1) is between 250 and 300 K, which limits the emitted radiation to the Planck functions of those temperatures. This limit of " 280 W m´2 is called the Simpson-Nakajima radiation limit (Simpson 1928;Nakajima et al. 1992). When the surface temperature exceeds 1600´1800 K, the temperature of the emitting layer (τ " 1) reaches 400 K in the 4 µm water vapor window. The new emerging emission peak allows the planet to emit more than 280 W m´2 and the planet cools efficiently again.
We treat this effect by adapting the grey atmosphere model to set the outgoing flux to 280 W m´2 if the surface temperature is in the range of 600´1800 K and the water pressure is greater than 10 mbar. If the absorbed stellar radiation is larger than this flux, the surface temperature will stay constant at " 1800 K until the atmosphere is desiccated due to the atmospheric escape of hydrogen or the absorbed stellar radiation drops below the runaway greenhouse limit.
To benchmark our adapted grey atmosphere model, we used the petitCODE (Mollière et al. 2015(Mollière et al. , 2017 to calculate the net outgoing flux from a pure steam atmosphere. The petitCODE is a onedimensional planetary atmosphere code that calculates the vertical structure of an optically thick (at the bottom) planetary atmosphere using the correlated-k method. It includes condensation and cloud formation as well as an equilibrium chemistry model for a wide range of atmospheric compositions.
We ran simulations for a range of surface temperatures and atmosphere pressures and calculated the Outgoing Longwave Radiation (OLR) for a 100% water vapor atmosphere on the extrasolar Super-Earth GJ 1132b (Berta- Thompson et al. 2015). Based on this, we created an interpolated look-up table for a given surface temperature and atmospheric pressure to calculate the OLR during runtime and implemented it into our model.
In Section 3.1, we show the magma ocean evolution on GJ 1132b for both radiative transfer models (grey & petitCODE) and compare the results to those from Schaefer et al. (2016) who used a line-by-line radiative transfer model.

Volatile Model
The volatile model controls the exchange of water and oxygen between the different layers of the system. Unlike for a cool planet with a solid surface, the partitioning of volatiles between the interior and atmosphere are in equilibrium during a magma ocean phase. Most importantly, the mass fraction of water in the melt is set by the partial pressure of water in the atmosphere. Therefore, the system is characterized by, and very sensitive to, sources and sinks.
At the base of the magma ocean, melt solidifies due to the cooling of the planet, locking water in the solid rock. However, the solid rock can only keep a small part of the water that was dissolved in the melt. The majority of the water stays in the melt increasing the mass fraction of water and, therefore, water is outgassed into the atmosphere to increase the pressure and restore the equilibrium.
At the top of the atmosphere water is lost into space because the star emits XUV radiation which dissociates water in the atmosphere into hydrogen and oxygen. While the hydrogen and part of the oxygen escape to space, the rest of the oxygen stays in the atmosphere and dissolves into the magma ocean where it oxidizes FeO to Fe 2 O 3 . When all of the FeO in the melt is oxidized or the mantle has solidified completely, oxygen starts to build up in the atmosphere.
The volatile model uses the thermal model to calculate the amount of water and oxygen that is stored in the atmosphere + magma ocean system (MOA) and the solid mantle (SOL). We assume that in the beginning the MOA contains all the water and no free oxygen is present. In equilibrium, the partial pressure of water vapour (in Pascals) in the atmosphere p H 2 O and the water mass fraction in the magma ocean F H 2 O are related by (Papale 1997): An overview of the most important parameters of the volatile model is provided in Table 2. The values used in our simulations are provided in the appendix in Table 7. The mass balance of the water contained in the magma ocean and atmosphere system is given by the following equation: As already mentioned, there are two sinks of water in the system: the solidification of the magma ocean which traps water in the solid mantle and the photolyzation of water followed by atmospheric es-cape of hydrogen and oxygen at the top of the atmosphere. Therefore, the water mass in the different reservoirs is governed by two differential equations, taking into account both sinks: dM moa We assume a constant, mantle-averaged partition coefficient k H 2 O for water between the melt and solid. Even though the value of this parameter should vary for different materials, for the purpose of this simplified model a constant coefficient is sufficient. Due to the XUV flux from the star, the water is dissociated and the hydrogen escapes at the rate φ 1 . The remaining oxygen goes into the magma ocean and oxidizes FeO to Fe 2 O 3 . The oxygen can be dragged along with the hydrogen and can escape at a rate φ 2 in the case of very high hydrogen loss rate 8 , which is likely to happen for planets in close orbits around early M dwarfs due to the particularly high XUV flux emitted by these stars. The escape rates of hydrogen and oxygen, φ 1 and φ 2 , respectively, are calculated by the VPLanet module AtmEsc, which is based on the atmospheric escape model of .
The amount of oxygen in the different reservoirs is governed by a similar set of differential equations: However, unlike water which dissolves into the magma ocean following Eq. (11), the oxygen reacts with the FeO in the melt, oxidizing it to Fe 2 O 3 . The composition used by Schaefer et al. (2016) and in MagmOc is the Bulk Silicate Earth (O'Neill & Palme 1998). To simplify our model, we assume that no oxygen builds up in the atmosphere as long as the reservoir of FeO in the magma ocean is not fully oxidized. When there is no FeO left in the magma ocean, no more oxygen dissolves into the melt but it builds up in the atmosphere.

Validation
In the following, we present MagmOc simulations for two scenarios to validate our model: one set-up for the Super-Earth GJ 1132b and one for Earth and compare to previous work by Schaefer et al.

GJ 1132b
GJ 1132b is a Super-Earth, with a size of 1.2 Earth radii orbiting an M dwarf at a distance to Earth of about 12 pc. With an equilibrium temperature of 409 K, assuming an albedo of 0.75, it is too hot to sustain liquid water on the surface and be habitable. Still, it is an interesting target for observing atmospheric composition and Albedo (steam atm.) 0.75 a assumption, r c,C " 3400 km dynamics on a rocky exoplanet (Berta- Thompson et al. 2015). The physical parameters for GJ 1132b used for this work can be found in the Appendix in Table 3.
In order to test MagmOc, we compare the results of VPLanet for GJ 1132b to those from Schaefer et al. (2016), in the following called Schaefer model. We calculated the atmosphere and magma ocean evolution for 0.1 -200 terrestrial oceans and compared the results using petitCODE (Mollière et al. 2015(Mollière et al. , 2017 and the grey model (Elkins-Tanton 2008; Carone et al. 2014) for the atmosphere cooling part with those obtained by Schaefer et al. (2016). Figure 2 shows as an example of the resulting cooling curves of the mantle temperature for 100 TO initial water content from the Schaefer model and VPLanet with both atmospheric models: the grey atmosphere and the petitCODE. For most of the evolution, the mantle cools more slowly in the Schaefer model but it solidifies earlier compared to our model when the atmosphere is desiccated. Only for the first Myr does the Schaefer model cools faster. Note, however, the logarithmic time axis which skews the perception towards the very early evolution.
The reason for the different cooling behavior is different implementations of the climate models: the Schaefer model uses a line-by-line radiative climate model with a simplified temperatureatmosphere profile at the base that follows a dry adiabat from the surface (set by the magma ocean temperature) to the tropopause. We use a similar temperature structure but a simplified radiative treatment in our grey atmosphere model. The petitCODE uses the correlated-k method (Mollière et al. 2015(Mollière et al. , 2017 for radiative transfer and iterates the atmospheric pressure temperature profile selfconsistently until emission and absorption are in balance for each vertical atmosphere layer in the radiative region. Below that, in the convective region, the atmosphere follows a moist adiabatic profile. The petitCODE produces a cooling curve that matches more closely the results from the Schaefer model than the grey atmosphere does, especially after the first million years. But it leads to a later solidification of the magma ocean like the grey atmosphere model. In other words, even though the cooling curve is much steeper with the petitCODE and the grey atmosphere model compared to the Schaefer model for large parts of the simulation, i.e. the first 10 Myr, the magma ocean needs more time to solidify. That is because the time of solidification mainly depends on the atmospheric desiccation, i.e. the loss of all water from the atmosphere towards the end of the magma ocean state. Apparently, the atmosphere in the Schaefer model desiccates more rapidly than both our grey and petitCODE models. A possible reason for this deviation might be a different age of the star at the beginning of the magma ocean simulation. We use a stellar age of 5 Myr while Schaefer et al. (2016) seem to use no offset. Because we start the magma ocean evolution later, we receive a lower XUV flux during the evolution and   Fig. 7). The necessary input files to reproduce the results in this paper can be found on our GitHub repository linked in the figure captions: Õ GitHub thus need more time to desiccate the atmosphere. When starting our simulation with a stellar age of 0, the atmosphere is desiccated after 43 Myr which is only " 15% larger than the result of Schaefer et al. (2016). We argue that the estimated lifetime of the protoplanetary disk ("5 Myr, Ribas et al. 2014) is a more realistic starting time for the simulation of the magma ocean evolution. We note that we did not not take the effect of large impactors into account and how these could re-initiate the magma ocean stage.
We also examined the oxygen pressure in the atmosphere at the end of the atmospheric escape for different initial water masses and for the different models (Fig. 2, lower right). Here, we find in our models that the amount of oxygen that builds up in the atmosphere depends crucially on the slope of the cooling curve towards the end of the magma ocean evolution as it fully solidifies. Due to the faster cooling with the grey atmosphere and the petitCODE compared to the Schaefer model, the FeO in the magma ocean becomes fully oxidized much faster for the last few percent of the mantle solidifi-cation. Therefore, there is still water escaping from the atmosphere when the oxygen buffer in the mantle is filled and free oxygen starts to build up in the atmosphere.
In the grey atmosphere model, the planet cools down the fastest (see Fig. 2). Therefore, oxygen starts to build up already for initial water masses of 70 TO. With the petitCODE, only for initially more than 100 TO of water in the system, oxygen builds up in the atmosphere, while the Schaefer model does not start to build up oxygen for initial water masses smaller than 200 TO.
We did not run simulations for initial water masses larger than 200 TO with the petitCODE and 100 TO with the grey atmosphere. For one of the most extreme cases we considered (100 TO), we find that the steep cooling curve of the mantle also drives catastrophic outgassing that leads to a mass fraction of water in the melt of nearly 20%. Since even 20% is a significant fraction of the mantle and our model does not take into account changes in the density or the planet's radius, we did not run simulations with more than 100 TO with VPLanet/grey and 200 TO with VPLanet/petitCODE. For 1000 TO initial water mass, the mass fraction of water in the melt would reach even 100%, which we consider as unphysical in this set-up. To consider such high values, we would need to account for changes in density in the equations of state, which are currently not implemented. When cooling is slower, as is the case in the Schaefer model, the atmospheric escape can balance the outgassing better and the water fraction in the melt stays at a smaller level, even for very high water contents. Therefore, it was possible for Schaefer et al. (2016) to run simulations up to initial water masses of 1000 TO.
We conclude that the grey atmosphere model based on Elkins-Tanton (2008) and Carone et al. (2014) was successfully adapted and benchmarked with the correlated-k petitCODE (Mollière et al. 2015(Mollière et al. , 2017 to be applicable to a primordial steam atmosphere on a magma ocean world. Generally, we reproduce the results of Schaefer et al. (2016) with MagmOc using both the petitCODE and the grey atmosphere model (Fig. 2): • Larger initial water masses lead to longer solidification times, where our solidification times for both atmosphere models are consistently larger by a little less than one order of magnitude compared to Schaefer et al. (2016), • Abiotic oxygen builds up for about 100 TO initial water mass. The exact initial water mass limit, for which onset of oxygen build-up happens, varies within half an order of magnitude (between 80 TO for the grey atmosphere model to 200 TO for Schaefer et al. (2016)).
We will show henceforth MagmOc results for the TRAPPIST-1 planets and Earth using the grey atmosphere model. We also validated our magma ocean-atmosphere model for the case of the young Earth. Similar to other terrestrial planets, Earth is believed to have had a deep magma ocean due to all the energy deposited during accretion (e.g. Elkins-Tanton 2008;Lammer et al. 2018). During the solidification of this magma ocean, some of the interior water outgassed into the atmosphere. Unlike GJ 1132b, which is orbiting an M dwarf, Earth is orbiting a G star at a much larger distance and the Sun is not as active as an M dwarf. Therefore, the atmospheric escape on early Earth was much weaker than on the young GJ 1132b.

Earth
Earth's initial water content is weakly constrained because even the Earth's current total water content is uncertain. In addition to the one Earth ocean of water located on the Earth's surface, much more water might be locked for example in the Earth's mantle Hamano et al. (2013) present simulations of terrestrial planets, including Earth, with a water atmosphere containing 5 TO and in a magma ocean phase: Their results show that Earth solidifies rapidly, after " 4 Myr for an initial water mass of 5 TO (T ini " 3000 K, A " 0.3). Figure 3 shows the results of simulations with MagmOc for Earth with the same initial water content (T ini " 3000 K, A " 0.75). The red arrows indicate how the solidification of the magma ocean, i.e. the decreasing size of the melt reservoir, leads to the very strong or catastrophic outgassing of water into the atmosphere (Elkins-Tanton 2008;Lammer et al. 2018).
Even though we use a different model for the outgoing flux, which leads to a different shape of the cooling profile, the solidification time differ only by a factor of 4 (MagmOc: t sol " 0.9 Myr, Hamano et al. (2013): t sol " 4 Myr). Also, the resulting water masses in the different reservoirs are the same for both simulations (See Fig. 3 and Hamano et al. (2013, Fig. 1)).
For GJ 1132b, the magma ocean is initially prevented from effective cooling by the extreme greenhouse effect of the water vapor atmosphere, since the planet is so close to its star. It can, however, cool off after atmospheric escape of hydrogen has reduced the atmospheric mass. Thus, atmospheric escape is in this case the driving factor for magma ocean solidification. In our GJ 1132b simulations, we tend to observe similar but slightly later (within one order of magnitude) solidification times compared to Schaefer et al. (2016). Since Earth is farther away from its host star than  (2008) presents simulations of the evolution of a magma ocean on Earth for different amounts of water and CO 2 in the atmosphere. For our simulations of the early Earth, we used the values as listed in Table 4 that are based on Elkins-Tanton (2008). We also chose a magma ocean depth of 2000 km to be consistent with Elkins-Tanton (2008). However, we did not include the effect of CO 2 in our model. We compare our results with a pure water atmosphere to her results in Table 5.
Especially for the water rich case (20 TO water and 4 TO CO 2 ), the results are very similar and differ by less than a factor of two compared to our Earth simulation with a pure water atmosphere simulation. We derive, however, consistently larger final atmosphere pressures and solidification time. We attribute these differences mainly to a different prescription of equations of states for H 2 O, because we adopted the density profiles used by Schaefer et al. (2016), which differ from those used by Elkins-Tanton (2008).
The results presented in this section are qualitatively -within one order of magnitude -similar to Schaefer et al. (2016), Elkins-Tanton (2008), andHamano et al. (2013). This shows that our model is able to simulate the magma ocean evolution of different kinds of rocky planets: in-and outside of their star's habitable zone. The quantitative differences are within observational uncertainties. We therefore conclude that our model successfully reproduces past results and is validated.
We further note here, that for both, GJ 1132b and also for Earth, we find that only 3´5% of the initial water will be locked in the mantle after the magma ocean solidified. Thus, it appears that for Earth, the intermediately wet formation scenario regime (10 -100 TO) is favored to derive the current low water content of 3-10 TO. This conclusion, however, comes with the caveat that we did not consider the influence of the Moon impactor and other smaller impactors. E.g. Zahnle et al. (2020) recently showed that smaller impactors can substantially change mantle and atmosphere chemistry and remove water from the system. Also, the moon forming impact was only the last of multiple events leading to large-scale mantle melting (Jacobson et al. 2014;Bottke et al. 2015). Two episodes of magma oceans are needed to explain the isotopic difference of 3 He/ 22 Ne between the depleted mantle and a primitive deep reservoir. The D/H ratio of sea water indicates that the silicate Earth during the last global magma ocean was oxidized enough to outgas H 2 O and CO 2 and not H 2 and CO (Pahlevan et al. 2019). The latter, reducing atmosphere mixture would heave facilitated atmospheric escape of hydrogen leading to a higher D/H ratio. Similarly, the D/H ratio implies a short-lived steam atmosphere, i.e. a short solidification time of the magma ocean (Stüeken et al. 2020). Both circumstances, the small effect of atmospheric escape of hydrogen and the short solidification times, are represented in our model. We now turn to applying MagmOc to the potentially habitable planets of the TRAPPIST-1 system.

Results for the TRAPPIST-1 planetary system
Since TRAPPIST-1 e, f, and g are orbiting TRAPPIST-1 in its habitable zone ( (Table 6). With these constraints, it is possible to give estimates on the density and composition of the planets (Dorn et al. 2018;Unterborn et al. 2018b; Barr et al. 2018).
To have a better understanding of the range of water fractions in TRAPPIST-1 e, f, and g, depending on several random parameters such as the interior and surface temperatures, we used for comparison our interior structure model (Noack et al. 2016). With this model, we investigated the range of compositions for TRAP-PIST-1 e, f, and g based on the masses and radii as listed in Agol et al. (2021) for a 1σ error range. Our interior structure model solves composition-dependent equations of state for iron, different rock assemblages, high-pressure ices and liquid water. The tran-  sition from liquid water to ice as well as the correct high-pressure ice phase (ice IV, VI, VII, VIII or X) is determined self-consistently for an adiabatic temperature profile. We assume arbitrary compositions for the three planets, where the iron weight fraction and the water fraction are randomly selected between 0 and 70 wt%. We also account for different iron distributions between mantle and core, where we allow between 0 and 20% of the mantle minerals to be iron-based. Other parameters, such as mass (within the known error range), surface pressure (between 1 and 100 bar), interior temperature profile and the surface temperature (effective temperature plus up to 50 degrees to account for greenhouse gases) are also selected randomly. Figure 4 shows the entire range of radii that are predicted for TRAPPIST-1 e, f and g by our model plotted over the input mass. The colour coding indicates the water content, and the size of the dots show the planet iron fraction. Grey circles indicate the range of masses and radii from the estimates by Agol et al. (2021) reported in Table 6. The average predicted water concentration together with a 1σ standard deviation and the total interval of observed water fractions is given in Table 6.
This evaluation as well as the comparison to published data reported in Table 6 show, that the uncertainty in the current water content remains large for the TRAPPIST-1 planets. We thus simulate the magma ocean and volatile evolution of the planets TRAP-PIST-1 e, f, and g for initial water masses from 1´100 TO until the planets become desiccated or until atmospheric escape stops.
In this work, we simulate the dominant process that leads to significant atmospheric erosion during the pre-main sequence of the host star, which is contained in the VPLanet module AtmEsc, described in more detail in Barnes et al. (2020). More precisely, AtmEsc uses energy-and diffusion limited hydrogen and oxygen escape following ; . Hydrogen is produced when the absorbed XUV flux of the host star photodissociates water, where we use an XUV absorption efficiency of XUV " 0.3 9 (see also Table 6). For the TRAPPIST-1 system, we use the stellar evolution model contained in the VPLanet module STELLAR. It assumes the constant-followed-by-power-law model of Ribas et al. (2005) for XUV-flux and calculates the total bolometric luminosity from grids based on Baraffe et al. (2015). The stellar luminosity of TRAPPIST-1 decreases with time. At some point during the evolution, typically after the magma ocean solidification, the incoming flux at TRAPPIST-1 e, f, and g reaches the threshold S eff , where we use equation (2)  To set up the model for the TRAPPIST-1 simulations, we choose the following specific assumptions as input for the equations described in Sections 2.1. For the abundance of radioactive isotopes we use Earth abundances, scaled by the mass of the planet. In addition, we consider tidal heating due to eccentric orbits with fixed eccentricities, as shown in Table 6. We did run simulations with different initial rotation periods and found that rotational tidal heating does not effect the long term evolution of the magma ocean as the rotation decays over a very short time scale. To decrease computational requirements, we therefore decided to start our simulations with the planets tidally locked and use a tide model with a constant phase lag (CPL) and set Q to 100. A concise compilation of the input parameters for the TRAPPIST-1 planets and the settings for VPLanet can be found in Table 6.

Different Evolution Scenarios
We identify in this work three evolutionary scenarios for the TRAP-PIST-1 planets, depending on the initial water content and their distance to TRAPPIST-1: These are scenario 1 (dry, ă 2 TO for planet g), 2 (intermediate, 2´100 TO), 3 (extremely wet, ą 100 TO). The left panel of Fig. 5 shows the solidification time of the magma ocean for TRAPPIST-1 g for initial water masses of 1´100 TO. Colored areas show how long the atmospheric escape continues (yellow) and whether the evolution results in a dry (red) or wet (blue) atmosphere. For each of the three scenarios, representative cases, described in more detail in the next subsections and in Fig. 6, are indicated by vertical lines. This figure further shows that a TRAPPIST-1 g scenario with low water content and short magma ocean solidification times runs the risk of complete desiccation (red region). Higher water content and solidification times larger than 14 Myr, on the other hand, would keep water vapour in the atmosphere after the end of the magma ocean state.
In addition, we show the solidification times and input water ranges for the three different scenarios for all three planets in the right panel of Fig. 5. Note that the border between the scenarios moves to lower water masses with increasing distance to the star. Figure 6 shows the evolution of the temperature and the atmospheric pressures for three cases identified in Fig. 5 (left panel) with numbers indicating important phases of the evolution. These phases are described in detail in the following for the three representative cases indicated in Fig. 5

(left):
Example for scenario 1 (dry): Initial water content 1 TO for TRAPPIST-1 g An example of the dry evolution (scenario 1) for TRAPPIST-1 g is shown in Fig. 6, top panel. The evolution of the temperatures and volatile reservoirs passes through four stages: I ă 0.1 Myr: Rapid onset of solidification and outgassing The solidification of the magma ocean leads to a rapid increase in the atmospheric pressure as the solid rock cannot store as much water as the liquid melt.
II ą 0.1 Myr: Runaway Greenhouse When the surface temperature drops to 1800 K, the ab-  6 for more information). Right: Solidification and desiccation times for planets e, f, and g (solid lines in red, yellow and blue, respectively). Range of the different scenarios is indicated with colored bars. The grey arrow indicates the distance to the star which increases from planet e to g. The range of the different scenarios depends on the orbital distance of the planet to the star which determines the incoming stellar radiation and therefore the cooling flux. Õ GitHub sorbed stellar radiation is larger than the greenhouse limit of 280 W m´2 (Goldblatt et al. 2013). Thus, the planet is not able to cool any further.
III Atmospheric escape drives outgassing Since water is escaping at the top of the atmosphere and the magma ocean has not yet completely solidified (the magma melt fraction ψ is still larger than 0.4), more water must outgas into to the atmosphere to maintain the equilibrium between the pressure of water vapour in the atmosphere and the water mass fraction in the magma ocean. However, since atmospheric escape is more efficient than the outgassing from the interior, the pressure in the atmosphere decreases rapidly.
IV,V " 6 Myr: Atmospheric desiccation & solidification Once the atmosphere loses all its water, the planet's cooling rate is no longer restricted by the runaway greenhouse limit and the magma ocean cools rapidly. We stop the evolution when atmospheric water pressure becomes zero, which prevents abiotic oxygen accumulation.
In scenario 1, the planet ends its magma ocean evolution with neither water vapour nor abiotically created O 2 in the atmosphere. We did not consider here the effect of N 2 or CO 2 gas left in the system, which could form the main constituents of the atmosphere in this scenario. II ą 0.5 Myr: Runaway greenhouse Further, cooling of the magma ocean surface to temperatures below 1800 K is again prohibited by the runaway greenhouse limit.
II Atmospheric escape drives outgassing Escape at the top of the atmosphere drives water out of the magma ocean. Since the atmosphere is very thick compared to scenario 1, atmospheric escape is in this scenario not able to decrease the pressure efficiently as long as the mantle is not yet solidified.
V " 14 Myr: Mantle solidification The stellar flux has decreased such that the planetary thermal emission (which is in energy balance with the input stellar irradiation) drops below the greenhouse limit. This allows the magma ocean surface to cool, the temperature to drop below 1800 K and the mantle to solidify. The magma melt fraction ψ becomes smaller than 0.4, marking the end of the magma IV " 50 Myr: Atmospheric desiccation The solidification of the mantle stops the efficient outgassing of water into the atmosphere. Atmospheric escape now removes water from the atmosphere efficiently. At the same time abiotic oxygen builds up.
In scenario 2, the planet ends its magma ocean evolution with a partially or completely eroded steam atmosphere and as a direct consequence with a significant amount of abiotically created O 2 (ą 100 bar). Since the magma ocean is already completely solidified at this point, O 2 can not be efficiently removed from the atmosphere in our model anymore. In the specific case depicted in the middle panels of Fig. 6, with a relatively low initial water mass (5 TO), the steam atmosphere is completely eroded and only oxygen remains.
Example for scenario 3 (extremely wet): Initial water content 100 TO for TRAPPIST-1 g Here, an example of an evolution track for the very wet evolution (scenario 3) is shown for TRAPPIST-1 g, where we depict the magma ocean and surface temperature, water vapour and O 2 atmosphere content in the bottom panel of Fig. 6.

I ą 1 Myr: Onset of solidification and outgassing
The extremely thick steam atmosphere slows down the solidification rate such that the runaway greenhouse limit (II) is not reached before the end of the simulation in this scenario. Efficient outgassing that increases the atmospheric pressure is also delayed.
V " 20 Myr: Mantle solidification The mantle eventually solidifies, slightly later than in scenario 2 because of the slow cooling rate through the thick steam atmosphere.
IV Atmospheric escape continues By now, the stellar XUV flux is very low and thus the atmospheric escape is inefficient. The atmospheric pressure decreases only sightly. At the same time, abiotically created O 2 builds up. However, since this scenario started with a much larger water reservoir than scenario 2, the final oxygen pressure is much lower than the remaining water pressure.
VI " 76 Myr: H 2 O atmospheric escape stops The stellar incoming flux at the location of the planet drops below the moist greenhouse threshold. H 2 O condenses out, the stratosphere dessicates and the atmospheric escape of H 2 O stops. A thin magma ocean is still present when the planet enters this regime. However, the viscosity of the remnant magma layer is so high that it cannot efficiently heat the atmosphere anymore. Water condensation at the surface will likely lead to rapid solidification of the remnant magma ocean at this point.
In scenario 3, the planet ends its magma ocean evolution, when an extremely thick steam atmosphere (ą 10000 bar) with about 1% free O 2 starts to condense out to form an ocean, leaving potentially a 100 bar O 2 atmosphere behind. However, we did not implement processes that could draw down O 2 into the mantle or into a water ocean after solidification.

Influence of Additional Heat Sources
We also tested the influence of additional internal heat sources on the magma ocean evolution of the TRAPPIST-1 planets by increasing the radiogenic and tidal heating. We found that we need to increase the internal heating enormously to generate substantial changes in the magma ocean evolution.
We implement this extreme heating by fixing the eccentricity to 0.1, thus mimicking the perturbation so of other planets, and by setting the abundance of 40 K to 1000 times the abundance of Earth. These adjustments are likely unrealistic, but serve as a useful bound on the plausible magma evolution of these worlds. The last two rows of Table 6 show the heating powers for both heating scenarios (reference and extreme).
Since TRAPPIST-1 e is closest to the star, the effect of the increased tidal heating is largest for this planet; the heating power in the high, undamped eccentricity case is larger by nearly three orders of magnitude compared to the nominal case with low, dampening eccentricity. The proximity to the host star also means that TRAPPIST-1 e will be most affected by additional heat sources via star-planet interactions such as induction heating (Kislyakova et al. 2017), which are not modeled here. Figure 7 shows the influence of the extreme heating on the magma ocean evolution and the volatile budgets of TRAPPIST-1 e. After an initial cooling phase, the internal heat flux is able to balance the cooling flux through the thick steam atmosphere. The mantle does not continue to cool but a magma ocean remains on the surface at the end of the simulation (after " 250 Myr). In other words, a prolonged magma ocean stage is an extreme version of scenario 3, but where the surface remains hotter than 2000 K. The prolonged magma ocean stage also affects the evolution of the volatile budgets: most notably, no oxygen builds up in the atmosphere. Even though atmospheric escape produces large amounts of oxygen, the liquid magma ocean can store all of it by oxidizing FeO to Fe 2 O 3 . Furthermore, the high surface temperature prevents the formation of surface oceans. Thus, the planet does not enter evolution phase VI even after 250 Myr. Phase VI (see Fig. 6, bottom left panel) is determined by low incoming stellar irradiation that allows for condensation of H 2 O, if the surface is cool enough, thus removing H 2 O from the upper atmosphere and stopping atmospheric escape (moist greenhouse limit).
Due to the hot surface that prevents condensation, about 10% of the total water content that is not dissolved in the magma ocean remains as water vapour in the atmosphere, where it undergoes continuous atmospheric escape. This atmospheric escape reenables outgassing from the magma ocean (evolution phase III). In other words, the prolongation of a thick magma ocean keeps water vapour in the mantle longer, but it also leads to a prolongation of the phase with atmospheric escape driven by outgassing and, thus, counter-intuitively, to further H 2 O loss from the system. At the same time, the prolonged presence of a thick magma ocean efficiently prevents the build-up of O 2 in the atmosphere.

Overview of TRAPPIST-1 e, f, g results
A condensed overview of the results obtained from all our simulations are presented in Fig. 8: The upper left panel shows the mantle solidification and atmospheric desiccation time of the planets. We define the desiccation time to be when no water is left in the atmosphere. However, oxygen might still be present in the atmosphere.   (Table 6). For larger initial water masses, water will remain in the atmosphere at the end of the simulation. Upper right: Total water remaining in the system (solid) and water locked in the solid mantle (dashed). The grey line indicates the initial water mass as comparison and upper limit. Lower left: Partial water pressure in the atmosphere at the end of the simulation. Lower right: Partial oxygen pressure at the end of the simulation. ÕGitHub The time when the atmospheric escape of H 2 O stops is indicated with a dot. For small initial water masses (e.g. the dry scenario 1), the planet solidifies quickly once the atmosphere is desiccated. In this scenario, the more water is in the system, the longer it takes to dessicate the atmosphere and thus to solidify the magma ocean. If the initial water mass is large enough (the intermediate wet scenario 2 for more than 5´10 TO), there is still a thick H 2 O atmosphere present when the stellar flux has decreased such that the absorbed radiation at the location of the TRAPPIST-1 planets drops below the runaway greenhouse limit. In this case, the outgoing thermal radiation of the planets is limited to 280 W m´2 until the atmosphere desiccates (Section 2.1.4, Goldblatt et al. 2013). From that point on, the planets are able to cool even with a thick atmosphere. The O 2 that is building up in the atmosphere has no greenhouse effect. This process prevents the solidification time from changing for a large range of initial water masses in scenario 2. It only depends on the distance of the planet to the star. For larger initial water masses (ě 100 TO for TRAPPIST-1 g), the thick steam atmosphere prolongs the lifetime of the magma oceans until they reach the limit of the runaway greenhouse, which can occur after 80 Myr for TRAPPIST-1 g and after 240 Myr for TRAPPIST-1 e, if the surface is cool enough for H 2 O condensation to occur (the extremely wet scenario 3). The upper right panel in Fig. 8 shows the remaining water in the system at the end of our simulations, either when the atmosphere is desiccated or the water condenses out to form an ocean. In addition to the total water amount, the amount of water locked into the mantle is shown. Furthermore, we show, as an upper limit, the initial water mass. Up to initial water masses of " 10 TO (dry scenario 1), all the water in the atmosphere is lost. Only a few percent of the inital water mass remains locked in the solid mantle, where it could be outgassed via volcanism in its further evolution. When starting with more water (intermediate scenario 2 and very wet scenario 3 for TRAPPIST-1 g), some will remain in the atmosphere at the end of the magma ocean stage. We caution that, due to the logarithmic scale, it may appear that almost all water is retained for 100 TO initial water mass after the evolution, but actually 20% of the initial water mass (20 TO) is also lost in this case. In absolute terms, the planets in the very wet scenarios (upper range of scenario 2 for TRAPPIST-1 e, f and scenario 3 for TRAPPIST-1 g) lose approximately the same amount of water or more than in the dryer cases. Still a much larger fraction of the initial water content remains in the planets (atmosphere and mantle) compared to the very dry scenarios 1.
The lower left panel in Fig. 8 is very closely related to the upper right one: it shows the partial water pressure in the atmosphere which is left in the atmosphere when atmospheric escape of H 2 O stops because water vapour does not reach the upper atmosphere anymore.
The lower right panel shows the oxygen pressure that builds up in the atmosphere due to photolysis of H 2 O and the subsequent escape of hydrogen. For initial water masses smaller than 5´10 TO (dry scenario 1), the produced oxygen can be stored completely in the mantle, bound in Fe 2 O 3 . For larger initial water masses (scenarios 2 and 3) up to 1000 bar of oxygen can build up in the atmosphere. We note that with increasing distance to the star the atmospheric escape of hydrogen is less efficient and less O 2 builds up in the atmosphere. Furthermore, for TRAPPIST-1 g we notice that the final oxygen pressure drops in the very wet scenario 3. The prolonged lifetime of the magma ocean allows more oxygen to dis-solve and oxidize the melt than in the intermediate scenario 2.

Discussion
In this work, we presented simulations of the magma ocean and volatile evolution for the planets TRAPPIST-1 e, f, and g. Studies by Dorn et al. (2018) suggest a very large water content of the TRAPPIST-1 planets with water mass fractions of up to 15%, calculations with the interior structure model by Noack et al. (2016) even higher mass fractions of up to 18%, while Agol et al. (2021) suggest a much drier composition with a water mass fraction of less than 5%. In Fig. 9, we show an overview of the final water content of the planets TRAPPIST-1 e, f, and g from all our simulations: for initial water masses from 1´100 TO and for two different heating scenarios -reference: Earth's abundance of radioactive iostopes and small, damped eccentricities; extreme: 1000 times Earth abundance of 40 K and fixed eccentricities (e " 0.1). The atmospheric composition is 100% H 2 O.
We compare our results to estimates of the current water mass by Dorn et al. (2018) and with results from the model by Noack et al. (2016) (Fig. 9, right). This comparison suggests a very wet formation, i.e. with an initial water mass of more than 100 TO, especially for the planets f and g. As described in Sec. 3, we did not run simulations for initial water contents of more than 100 TO. For such large water masses, the water fraction in the melt becomes very large and eventually 1, which is unphysical. To simulate these scenarios correctly, an adaption of the physical equations would be necessary. To get an impression of how the magma ocean evolution would continue for larger water contents, we extrapolated the remaining water budgets up to an initial water mass of 1000 TO. We based these extrapolations on the assumption that an increase in the water content of the mantle and the atmosphere would not increase the lifetime of the magma ocean. However, we can see already for TRAPPIST-1 g that 100 TO of water slightly prolongs the magma ocean lifetime and damps the build-up of oxygen. Therefore, we can conclude that the extrapolations presented in Fig. 9 show only upper limits.
We rather expect the extremely wet scenarios with ą100 TO to behave similar to the cases including extreme heating. For planet g, the solidification time scale rises with more than 70 TO initial water content, whereas the magma ocean solidification time remains relatively constant at about 14 Myr for 2´70 TO (see Fig. 8, upper left panel). The extreme heating cases were found to significantly extend the magma ocean life time (to at least 250 Myr for planet e) and to prevent abiotic oxygen build-up (Fig. 7). Further, the planet remains relatively long in a stage, where atmosphere erosion drives interior outgassing, since the surface is too hot to allow for condensation of water. This leads to an increase in the total water loss for a given planet compared to a case where the surface can solidify. With a solid surface, a water ocean can form on the surface, thus saving H 2 O from being photolysed in the upper atmosphere via condensation. Testing such a hypothesis with an improvement of our model, where we can extend the simulation beyond 250 Myr, and with a more physically plausible heating mechanism could provide further insight into habitability of the TRAPPIST-1 planets.
Even though we most likely overestimate the remaining water content in the system (especially for very wet cases), we can show that only 3´5% of the initial water will be locked in the mantle due to the small partition coefficient between melt and solid. Therefore, Simulations for all planets with 100% steam atmospheres (red, yellow, blue). Red and purple bars indicate range of initial water fractions that lead to abiotic O 2 build-up. Values for initial water fractions ą 20´30 wt% (grey area) are extrapolations on the assumption that solidification times do not increase with more water. Probability distributions show the current estimates for water content of the TRAPPIST-1 planets by Dorn et al. (2018) and those calculated with the interior structure model by Noack et al. (2016, this study) with 1σ error range for mass and radius (Table 6). ÕGitHub it is not possible in our simulations to sequester large amounts of water in the mantle without several thousand bars of water in the atmosphere or at the surface. Consequently, with a water content today of several hundred terrestrial oceans and only a few percent of it being locked in the mantle, planets f and g will most likely be worlds with very deep oceans. It should be noted that it has been postulated that instead of the magma ocean freezing from the bottom of the mantle to the top, due to the high temperature contrast between lower mantle and core, a basal magma ocean could remain during the magma ocean crystallization process (Labrosse et al. 2007). In addition, Nomura et al. (2011) suggest that such a basal magma ocean could have extended to up to 1000 km thickness due to enrichment in iron, making the magma denser than the solidifying rock above. Such a basal magma ocean could then in theory also be a reservoir for volatiles, which, on the other hand, leads to a decrease of the melt density and might therefore destabilize such a basal magma layer.
A 300 TO water body on the surface of TRAPPIST-1 g would result in a " 670 km deep ocean, most likely forming high pressure ice layers in the ocean (Dorn et al. 2018). It is still an objective of current research to determine whether such ocean worlds are potentially habitable (e.g. Noack et al. 2016). Most notably for the further evolution of the planetary atmosphere, deep oceans may suppress further outgassing due to the high temperatures and pressures at the ocean floor even if volcanic activity is present (Noack et al. 2016), but transport through the ice layer might still be possible (Kalousová & Sotin 2018). Any further evolution can therefore only be driven by chemical interactions between the oxygen-rich atmosphere and the deep water ocean. Glaser et al. (2020) suggest that, in the long-term evolution of a water world, the biotic and abiotic (by photolysis) rates of oxygen production might be very similar, further devaluing oxygen's role as a potential biosignature. However, the masses of the TRAPPIST-1 planets are still debated and our results may have to be revised in the near future.
We also show that the TRAPPIST-1 planets start to build up several hundred bars of oxygen for initial water masses larger than 5 TO. First estimates by Lingam (2020) indicate that once the magma ocean solidifies it would be very difficult to reduce abiotic oxygen build-up assuming present Earth sinks and sources. Though Harman et al. (2018) present that lightning might be an effective way of reducing the oxygen abundance in such an atmosphere, as long as water vapour is present in the atmosphere. We find, however, that in scenario 2 (5 TO for TRAPPIST-1 g, Fig. 6) it is in principle possible to build up significant amounts of free oxygen and completely erode water vapour via atmospheric escape. While such a dry oxygen-rich atmosphere is unlikely to occur for the very wet TRAPPIST-1 g, it could occur for the drier TRAP-PIST-1 e. If these high oxygen pressures are a physical reality, it is possible to detect them with JWST (Lincowski et al. 2018).
However, we also show that large heating rates can prevent the abiotic build-up of oxygen in the atmosphere and therefore could facilitate the habitability of these planets. Heating rates of 1000 times Earth abundance of 40 K seem extremely high. However, 40 K presents a powerful and continuous heat source for the long-term evolution (several 100 Myr) of a rocky planet and can thus help us to understand how a planet can sustain a magma ocean over a long time or even permanently. Also, there are further heat sources like additional tidal (e.g. via planet-planet tides that were neglected here) or magnetic heating (e.g. Kislyakova et al. 2017), which we did not consider in this work.
Our results favor a very wet formation of the planets TRAP-PIST-1 f and g, which provides necessary context to discuss how the delivery of such a high water mass fraction can be achieved during formation. E.g. Lichtenberg et al. (2019) show that radiogenic heating of 26 Al leads to desiccation in rocky planetesimals and postulates that only planet building blocks poor in 26 Al can lead to water-rich rocky planets. It will be needed to investigate, however, if water delivery via cometary material from the outer disk could explain a water-rich planet like TRAPPIST-1 g, even if the rocky planetesimals were water-poor. Miguel et al. (2020) provide dedicated formation scenarios for TRAPPIST-1, where they also propose that TRAPPIST-1 g likely formed very water-rich, TRAP-PIST-1 e could have formed dry to water-rich, and TRAPPIST-1 f lies in between with respect to water mass fraction. All these scenarios compare well with our comparison between current water content estimates after magma ocean evolution and initial water content before magma ocean evolution. However, we note that the uncertainties and systematic errors that underpin current water estimates are large. Thus, more work is needed to tighten the possible ranges of water content scenarios.
One assumption that can affect the magma ocean evolution is the starting time of the simulation, i.e. the age of the star at t " 0. As mentioned in Section 2.1, we assume an age of the star of 5 Myr. As a late impact such as the moon forming impact on Earth can reset the magma ocean evolution, this is an important parameter to discuss. However, the compact architecture of the TRAPPIST-1 system suggests that the disk was cleared of larger rocks very early in the system's evolution as a late impactor would have severely disturbed the stability of the system (Gillon et al. 2017;Tamayo et al. 2017). So far our model assumes efficient cooling due to the strong convection of the magma ocean. This is only true until the melt fraction at the surface drops below 0.4 and the viscosity of the melt becomes solid-like. Lebrun et al. (2013) proposed that for viscosi-ties below 0.4 the magma ocean could enter into a 'mush stage' where a thermal boundary layer (TBL) develops on the surface. A similar idea is presented by Schaefer et al. (2016). Other models either do not take this effect into account (e.g. Hamano et al. 2013) or stop the simulation of the solidification process at that point (Elkins-Tanton 2008;Nikolaou et al. 2019). Debaille et al. (2009) propose that a thick thermal boundary layer, including possibly even a lithosphere, may be neglected towards the end of the magma ocean since iron-rich minerals that are the last components to crystallize from the magma ocean likely lead to overturning near the surface and a re-setting of the thermal boundary layer.
Because we aimed in this work for a versatile and numerically efficient magma ocean description, we neglected the 'mush stage' towards the end of the magma ocean stage, following the argumentation of Debaille et al. (2009) andHamano et al. (2013). As Lebrun et al. (2013) note, the introduction of a TBL requires an additional iterative loop to balance the heat flux out of the TBL and the atmospheric heat flux. Nevertheless, it would still be interesting to investigate in a future update of the code how the volatile budget is modified towards the end of the magma ocean evolution by a TBL.
Even when taking into account the TBL, our model assumes that the convection in the mantle is described by a single Rayleigh number, which is a very simplifying assumption. As the viscosity strongly increases with decreasing melt fraction, a layered treatment of the convection would be more accurate. As a result, the thermal evolution and outgassing rates would slow down. However, since the viscosity increases dramatically only for the lowest layers of the magma ocean where the melt fraction is below 0.4, we would not expect our results to change for most of the mantle solidification. Only when the magma ocean becomes very shallow, most of the magma ocean will be at a melt fraction of around and below 0.4. This period is represented by the 'mush stage' in Lebrun et al. (2013).

Conclusion
We investigated the magma ocean evolution of the planets TRAP-PIST-1 e, f, and g and identified three general scenarios (Fig. 5): • Scenario 1, dry: The steam atmosphere prevents the mantle from solidifying because the absorbed stellar radiation is higher than the runaway greenhouse flux (Sec. 2.1.4 and Goldblatt et al. 2013). Only when the planet has lost all its atmospheric water via atmospheric escape within a few Myr does the magma ocean solidify. No oxygen builds up in the atmosphere because the magma ocean efficiently removes it from the atmosphere.
• Scenario 2, moderately wet: The magma ocean solidifies while there is still water left in the atmosphere if the absorbed flux has dropped below the runaway greenhouse limit. At this point, the magma ocean is no longer able to remove the oxygen from the atmosphere. As the atmospheric escape of hydrogen continues, several hundred bars of oxygen can build up in the atmosphere. At the time of mantle solidification, the atmosphere will consist of free oxygen and potentially water vapour. There is also the possibility that the magma ocean ends with an oxygen-rich and dry atmosphere.
• Scenario 3, extremely wet: The thick steam atmosphere prolongs the lifetime of the magma ocean, which absorbs the photolytically produced oxygen. At the end of the magma ocean stage, the planet ends up with a thick steam atmosphere and a smaller relative contribution of oxygen.
In Scenario 3, the magma ocean phase can, in principle, be extended up to 250 Myr for TRAPPIST-1 e, if the planet receives continuous interior heating of 2.8ˆ10 4 TW. The thick magma ocean prevents oxygen build-up and 90% of H 2 O remains dissolved in the magma ocean. The whole planet is, however, after 250 Myr in a state where water vapour is kept in the atmosphere due to the very hot surface (2000 K) and is eroded via photolysis. The water loss in the atmosphere is instantly replenished with efficient outgassing from the magma ocean that is still very thick (interior driven outgassing). That is, the planet still suffers continuous water loss even after 250 Myr due to the presence of a thick magma ocean. In contrast to that, without strong interior heating, the stellar irradiation would be so low and the surface so cool that water can condense at the surface, thus removing water from the atmosphere and keeping it from getting eroded. We found consistently for all investigated planets in this work (including Earth and GJ 1132b) that only 3´5% of the initial water will be locked in the mantle after the magma ocean solidified.
Comparison with the current estimates of the water content of the TRAPPIST-1 planets by recent studies (Noack et al. 2016;Barr et al. 2018;Dorn et al. 2018;Unterborn et al. 2018b) suggest that the potentially habitable TRAPPIST-1 planets have water mass fractions of 0-0.23, 0.01-0.21, and 0.11-0.24 for planets e, f, and g, respectively. This suggests that planet g followed scenario 3 and planet f scenario 2 or 3, meaning that both planets are most likely covered by a thick ocean and a potentially oxygen rich atmosphere. However, significant evolution after solidification, as well as processes not included in our model, can also significantly affect atmospheric abundances. Planet e, however, might have followed scenario 1 or 2 and could therefore be the most similar to Earth. The water-richness required to explain the high water content of TRAPPIST-1 f even after the magma ocean evolution, requires a very wet formation scenario. That may mean that the planetesimals from which this planet formed were poor in 26 Al (Lichtenberg et al. 2019). The possible initial water content that is necessary to start a magma ocean evolution and still end up with the current observed water mass fraction in the TRAPPIST-1 e, f, and g planets are also well in line with the formation model of Miguel et al. (2020). Those authors predict a high likelihood of water-richness for TRAPPIST-1 f after formation and that TRAPPIST-1 e could be formed dry -in line with our results for this planet when we try to reproduce the current water mass estimates by Dorn et al. (2018).
Since only 3´5% of the initial water will be locked in the mantle after the magma ocean solidified, we also conclude that for Earth the intermediate wet formation scenario regime (10 -100 TO) appears to be favored to derive the current low water content of 1-10 TO. However, we did not consider influence of the Moon impactor and other smaller impactors as Zahnle et al. (2020).
Our work also has major implications for the habitability of the TRAPPIST-1 planets, in particular for TRAPPIST-1 g. 300 TO water body on the surface of TRAPPIST-1 g, as appears likely from the combination of formation (Miguel et al. 2020), magma-ocean evolution (this work) and current water estimates (Dorn et al. 2018;Unterborn et al. 2018a) would result in a " 670 km deep ocean. Such deep ocean worlds are potentially uninhabitable (e.g. Noack et al. 2016). Furthermore, it could build up several hundred bars of oxygen at the end of its magma ocean stage, which could be very difficult to remove (Lingam 2020). Thick oxygen atmospheres, however, could be detectable with JWST (Lincowski et al. 2018), constraining evolution models further.
To conclude, MagmOc is a powerful tool inside the VPLanet code to simulate the thermal and volatile evolution of terrestrial (exo)planets during their magma ocean phase. It is especially applicable to simulate a large parameter space for various heating rates, orbital configurations, initial compositions, and stellar types.

Outlook
Even though our model successfully reproduces previous results for magma ocean evolution of Earth and the Super-Earth GJ1132b (Elkins-Tanton 2008;Hamano et al. 2013;Schaefer et al. 2016) and produces valuable results for the TRAPPIST-1 planets, there are numerous directions for future research.
The next step will be to add the powerful greenhouse gas CO 2 into the model. As the example of Venus shows, CO 2 can be an important constituent of atmospheres of rocky planets. CO 2 was also a major constituent of early Earth's atmosphere (e.g. Kasting 1993; Deng et al. 2020). The amount of CO 2 present in the atmosphere determines whether condensation of water to an ocean can occur. The higher the CO 2 partial pressure the more difficult it is for the water to condense (Lebrun et al. 2013;Massol et al. 2016;Salvador et al. 2017;Stüeken et al. 2020). To fully understand the early evolution of rocky planets, it is necessary to include CO 2 into this model. Our model so far can only simulate water composition up to 100 TO. Since TRAPPIST-1g has likely a higher water content than 100 TO, future work must reformulate the equation of state in a future version of MagmOc to coherently tackle extremely wet scenarios.
Furthermore, in this work, we assumed bulk silicate Earth composition, i.e. an oxidizing mantle composition. However, even for early Earth this assumption has come under scrutiny (Kasting et al. 1993b;Gaillard & Scaillet 2014;Ehlmann et al. 2016;Fegley et al. 2016;Schaefer et al. 2016;Del Genio et al. 2020;Armstrong et al. 2019). Further work may be warranted to understand if and under which conditions a reducing mantle composition could occur, which would lead to outgassing of CO and H 2 instead of CO 2 and H 2 O as assumed here (Katyal et al. 2020;Ortenzi et al. 2020). In addition to its effect on the redox state of the outgassed atmosphere, the mantle, and in particular the crust, composition have an important influence on the stability of liquid water on the surface after the magma ocean solidified (Herbort et al. 2020).
We invoked an extreme interior heating scenario for TRAP-PIST-1 e via unphysically high radiogenic heating rates to explore if a magma ocean can last more than 100 Myr. There are, potentially, more physical additional heat sources for the TRAPPIST-1 system, which we did not consider in the model, and which could achieve the necessary heating rate of 3ˆ10 4 TW. Due to the close proximity of planets orbiting M dwarfs to their host star, e.g. electromagnetic induction heating due to the inclination of the magnetic dipole of the star (Kislyakova et al. 2017) and enhanced tidal heating due to the decay of the obliquity (Heller et al. 2011) could be relevant. Also, Millholland & Laughlin (2019) note that obliquity tides in a tight planetary system like TRAPPIST-1 could create up to 3ˆ10 3 TW heating. Future work will show if these additional mechanisms can indeed provide sufficient, sustained heating to lead to a long magma ocean phase, preventing the abiotic build-up of oxygen as long as the magma ocean is thick enough. In this case, we will also need to extend our model to be able to perform simulations to tackle coupled magma ocean-atmosphere evolution over Gyrs time scale. In this context, we will also need to investigate how much water will be removed from the planet, because the prolonged magma ocean presence also prolongs water erosion.
The comparison between our simulations and the current water content of the TRAPPIST-1 planets are based on composition estimates from Dorn et al. (2018), which are in turn based on planetary radii and masses constraints of Delrez et al. (2018) and Grimm et al. (2018), respectively, as well as updated measurements by Agol et al. (2021). If new mass constraints occur that lead to new water content estimates, then these new results can be easily compared to our simulation results because we agnostically covered a large range of water abundance.
In this work, we only simulated the planet's evolution until the magma ocean stage stops. Another obvious further line of research is to extend our model to consider outgassing from the solidified mantle or the thermal evolution of the planet's interior after the end of the magma ocean phase (Driscoll & Bercovici 2013, Garcia et al., in prep.). Then we might be able to confirm that outgassing from the mantle and crust continues even when the atmosphere was completely eroded (our scenario 1) to build up a new, potentially habitable atmosphere (Godolt et al. 2019). This approach would also allow predictions for the atmospheric content today that could be tested by JWST (Wunderlich et al. 2019). However, this would require adding the treatment of the solidified part of the mantle to the thermal evolution model. Furthermore, taking into account the specific location of radiogenic and tidal heat production rather than using a mantle averaged approach would advance our model, as well as including the heat flux from the core into the mantle. This is also true for the partition coefficient of water between solid and melt which we assume to be constant throughout the mantle. However, this parameter depends on the material of the mantle and should therefore vary with depth.
We also did not take into account yet the influence of external objects like the Moon forming impactor and possible other smaller impactors that could further modify atmosphere and mantle composition as recently proposed by Zahnle et al. (2020) to explain the early Earth evolution. This may be another possible avenue to expand the model.
For exoplanets, the Super-Earth 55 Cancri e may be currently the most observationally accessible planet that has even today a magma ocean on its surface (Demory et al. 2016).The atmospheric composition of this planet is currently still debated based on tentative atmospheric HCN observations with the Hubble Space Telescope. HCN would imply a reducing, N 2 dominated atmosphere instead of an oxidizing atmosphere as assumed in this study (Tsiaras et al. 2016;Hammond & T. Pierrehumbert 2017;Zilinskas et al. 2020). Future observations with JWST or the ELT will target 55 Cancri e as an example to study the evolution of magma ocean planets and the influence of these magma oceans on the atmosphere. Future research should reproduce the atmosphere composition of this exoplanet with an adapted model, which also takes a reducing mantle composition into account.
Similarly, new results on Io may provide additional insight and constraints on models of magma oceans. While our work suggests TRAPPIST-1 e-g are not likely to possess magma oceans today, Io may still possess one (Khurana et al. 2011;Tyler et al. 2015). As new research on Io becomes available, it may also inform future interior modeling of the TRAPPIST-1 planets. In particular, NASA's Io Volcanic Observatory (McEwen et al. 2014), should it launch, will provide invaluable data in our understanding of extremely heated bodies. Future work could also apply MagmOc to Io for calibration and to predict features that may be detected by flybys.
The still debated composition of the TRAPPIST-1 planets as well as the lava world 55 Cancri e demonstrate that there is a real need to reformulate and expand geophysical models that were mainly developed for the Solar System planets, especially Earth. We have demonstrated that MagmOc not only has merit to simulate the magma ocean stage of Earth, shedding light on how our world became habitable, it also is versatile enough to be applicable to extrasolar planets like the TRAPPIST-1-planets. For rocky exoplanets, MagmOc provides the ideal link between formation models such as that of Miguel et al. (2020), interior structure studies like Noack et al. (2016); Dorn et al. (2018);Unterborn et al. (2018a), to habitability studies like that of Lingam (2020), as well as to studies which explore extreme geophysical scenarios such as Kislyakova et al. (2017) and Millholland & Laughlin (2019). MagmOc is thus a valuable tool to explore the diversity of rocky words in the Solar System and beyond.