Magma Ocean Evolution of the TRAPPIST-1 Planets

1. Introduction

As of March 2021, there are 4692 confirmed exoplanets¹. Of those exoplanets, 60 are optimistically classified as potentially habitable², 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_⊕ (Otegi et al.,2020) and smaller than 1.7 R_⊕ (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_⊕ 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, for example, Proxima Centauri b (Anglada-Escudé et al.,2016), TRAPPIST-1 e, f, and g (Gillon et al.,2016), and Teegarden's star c (Zechmeister et al.,2019), which are prime candidates for observations with upcoming space- and ground-based instruments like the James Webb Space Telescope (JWST) (Barstow and 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 ultraviolet irradiation (XUV) emitted by active M dwarfs leads to enhanced atmospheric erosion (Watson et al.,1981; Lammer et al.,2003, 2009; Erkaev et al.,2007; Owen and 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 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 weight percent, 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³,” for example, 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., Luger and Barnes, 2015). 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₂O rather than H₂.

Such a scenario was tackled with simulations for rocky Solar System planets initially for an early atmosphere dominated by water vapor by Matsui and Abe (1986), Zahnle et al. (1988), Abe (1993, 1997), Solomatov (2000), Zahnle et al. (2015), and by Elkins-Tanton (2008) for an early H₂O and CO₂ 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, uninhabitable super-Earth GJ 1132b. They found that around a pre-main sequence M dwarf, atmosphere erosion can be so effective that abiotic free O₂ builds up in the early atmosphere, as much as several thousands of bar.

Abiotic oxygen buildup on rocky planets around M dwarfs was already proposed by Luger and Barnes (2015), 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₂. They also did not take into account that additional heating sources like tidal heating (e.g., Driscoll and Barnes, 2015) and magnetic induction (Kislyakova et al.,2017) could potentially delay solidification of the magma ocean and thus reduce abiotic oxygen buildup.

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₂O-dominated atmosphere. We consider initial water contents from 1 to 100 terrestrial oceans (TO) 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 (Barnes et al.,2020) to include the physics and chemistry necessary to model the TRAPPIST-1 planets during the magma ocean phase. already includes “modules” for radiogenic heating (), tidal heating (), stellar evolution (), and water photolysis and hydrogen escape (). Here we introduce a new module that calculates the thermal evolution of a magma ocean, volatile fluxes through the surface, and radiative cooling through an H₂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 module of (Section 7).

The necessary input files to reproduce the results in this paper can be found in our GitHub repository linked in the figure captions. Colored versions of the figures can be found in the same repository and the version of this paper on arXiv.

2. The Model

In this section, we describe the setup of the module. We stress that we describe here a relatively simple initial setup 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.

2.1. Thermal model

Our thermal model is based in large part on the coupled stellar-atmosphere-interior model for terrestrial exoplanets by Schaefer et al. (2016) and for rocky Solar System planets by Elkins-Tanton (2008). If not otherwise noted, the equations in this section are 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 Appendix Table A1. If not otherwise noted, these values are adapted from Schaefer et al. (2016).

Table 1. Parameters for the Thermal Model

Symbol	Parameter
Tp	Potential temperature of the mantle
Tsurf	Surface temperature
Tsolidus	Solidus temperature
Tliquidus	Liquidus temperature
rs	Solidification radius
	Mantle bulk density
cp	Silicate heat capacity
	Heat of silicate fusion
Qr	Heat by radioactive decay [W kg−1]
Qt	Heat by tidal interactionsa [W kg−1]
F	Heat flux leaving atmosphere (OLR – ASRb)
	Thermal expansion coefficient
	Dynamic viscosity [Pa s]
	Magma ocean averaged melt fraction
	Surface melt fraction

^aNot included in Schaefer et al. (2016).

^bOLR = outgoing long-wave radiation; ASR = absorbed stellar radiation.

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, 2013; Driscoll and 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 stellar radiation is reflected back to space due to the albedo of a water vapor 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) gave 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 Section 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 = 400 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):

(1)

with the first term on the right-hand side describing the heat produced by the fusion of silicates (latent heat) and the second term describing the radiogenic and tidal heating; 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 the work of 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, that is, the average mantle temperature T_m. This temperature is related to the potential temperature by T_m = T_p with for planets of approximately 1 Earth mass (Schaefer and 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.

2.1.1. 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 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

(2)

with T_solidus,surf = 1420 K and T_{liquidus,surf} = 2020 K.

Lebrun et al. (2013) showed that  = 0.4 marks the phase transition toward a thicker, viscous thermal boundary layer and the effective end of the magma ocean stage⁵. For the viscosity profile, therefore, two different regimes of the surface melt fraction are considered: above the critical melt fraction  = 0.4, the viscosity parameterization from Lebrun et al. (2013) is used, while for 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

(3)

with , , the activation energy , and the ideal gas constant R.

This returns a viscosity of for a fully molten surface, and only for a surface melt fraction of just above 0.4 the viscosity increases to . We therefore conclude that our approximation of efficient convection is valid for most of the solidification process⁶. 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) showed that the effect of water-dependent changes of the viscosity on the solidification time of the mantle is negligible.

2.1.2. Mantle solidification

To calculate changes in the solidification radius r_s, we define the adiabatic temperature structure T(r) of the convective mantle with respect to the potential temperature T_p:

(4)

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):

(5)

where and for , and and for . The depth z dependent pressure is given by . The liquidus temperature is assumed to be larger than the solidus by 600. 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:

(6)

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, that is, , we set until the base of the mantle starts to solidify.

2.1.3. Radiogenic and tidal heating

Radiogenic heat in Earth is generated by the decay of ²³⁸U, ²³⁵U, ²³²Th, and ⁴⁰K in the core, mantle, and crust. The radiogenic power produced by species i in reservoir j is

(7)

where , is the half-life, t is time, and is the initial heat production at t = 0. This heat source is included in our model by connecting it to the module of (Barnes et al.,2020). Although the planets may form quickly enough for ²⁶Al heating to be important (Barnes et al.,2016; Lichtenberg et al.,2019), we ignore this effect here⁷.

Similarly, the tidal heating rate Q_t is calculated by the module of . 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 three-dimensional 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 “constant-phase-lag” (CPL) model in 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 multiplanet systems tidal interactions between individual planets might be important (Wright, 2018; Hay and 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 and O'Brien, 2002; Sasaki and Barnes, 2014).

2.1.4. 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 () and the top of the atmosphere () with the emissivity :

(8)

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

(9)

with the optical depth of the atmospheric species i

(10)

is the atmospheric mass of species i. To calculate the optical depth, we use the atmospheric absorption coefficient for water (Elkins-Tanton, 2008) at the reference pressure (Yamamoto and Onishi, 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 when greenhouse gases 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 ≲ 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⁻² 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⁻², and the planet cools efficiently again.

We treat this effect by adapting the grey atmosphere model to set the outgoing flux to 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, 2017) to calculate the net outgoing flux from a pure steam atmosphere. The petitCODE is a one-dimensional planetary atmosphere code that calculates the vertical structure of an optically thick (at the bottom) planetary atmosphere by 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 long-wave radiation 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 outgoing long-wave radiation 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 and petitCODE) and compare the results to those from Schaefer et al. (2016), who used a line-by-line radiative transfer model.

2.2. 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 is 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; 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 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₂O₃. When all 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 vapor (in pascals) in the atmosphere and the water mass fraction in the magma ocean are related by (Papale, 1997)

(11)

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 Appendix Table A1.

Table 2. Parameters for the Volatile Model

Symbol	Parameter
	Water mass fraction in liquid melt
	Water mass magma ocean + atmosphere
	Water mass crystals in magma ocean
	Water mass in the liquid melt
	Water mass in the atmosphere
	Water mass in the solidified mantle
	Water part. coeff. melt-solid
	XUV-driven atm. mass-loss rate of H [kg/m2/s]
	XUV-driven atm. mass-loss rate of O [kg/m2/s]

The mass balance of the water contained in the magma ocean and atmosphere system is given by the following equation:

(12)

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 escape 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:

(13)

(14)

We assume a constant, mantle-averaged partition coefficient 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 . The remaining oxygen goes into the magma ocean and oxidizes FeO to Fe₂O₃. The oxygen can be dragged along with the hydrogen and can escape at a rate in the case of very high hydrogen loss rate⁸, 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, and , respectively, are calculated by the module , which is based on the atmospheric escape model of Luger and Barnes (2015).

The amount of oxygen in the different reservoirs is governed by a similar set of differential equations:

(15)

(16)

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₂O₃.

The composition used by Schaefer et al. (2016) and in is the Bulk Silicate Earth (O'Neill and 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.

3. Validation

In the following, we present simulations for two scenarios to validate our model—one setup for the super-Earth GJ 1132b and one for Earth—and compare to previous work by Schaefer et al. (2016), Elkins-Tanton (2008), and Hamano et al. (2013), respectively.

3.1. 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 dynamics on a rocky exoplanet (Berta-Thompson et al.,2015). The physical parameters for GJ 1132b used for this work can be found in Table 3.

To test , we compare the results of for GJ 1132b to those from Schaefer et al. (2016), referred to in the following as the Schaefer model. We calculated the atmosphere and magma ocean evolution for 0.1–200 TO and compared the results using petitCODE (Mollière et al.,2015, 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 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 million years does the Schaefer model cool faster. Note, however, the logarithmic time axis which skews the perception toward 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 temperature-atmosphere 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, 2017) for radiative transfer and iterates the atmospheric pressure temperature profile self-consistently 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, that is, 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, that is, the loss of all water from the atmosphere toward 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) seemed to use no offset. Because we start the magma ocean evolution later, we receive a lower XUV flux during the evolution and 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 toward 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 solidification. 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 /grey and 200 TO with /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 setup. 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 the works of Elkins-Tanton (2008) and Carone et al. (2014) was successfully adapted and benchmarked with the correlated-k petitCODE (Mollière et al.,2015, 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 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 that of 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 buildup happens, varies within half an order of magnitude (between 80 TO for the grey atmosphere model and 200 TO for Schaefer et al.,2016).

We will show henceforth results for the TRAPPIST-1 planets and Earth using the grey atmosphere model.

3.2. Earth

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's initial water content is weakly constrained because even Earth's current total water content is uncertain. In addition to the one Earth ocean of water located on Earth's surface, much more water might be locked for example in Earth's mantle transition zone at a depth of 410–660 km (Pearson et al.,2014; Schmandt et al.,2014).

Hamano (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 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, that is, 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 differs only by a factor of 4 (: 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 those of Schaefer et al. (2016). Since Earth is farther away from its host star than GJ 1132b, the solidification time is mainly determined by the implementation of the cooling flux. In our Earth simulation, this leads to a similar but slightly earlier (by a factor of 3) solidification time compared to that of Hamano et al. (2013).

Elkins-Tanton (2008) presents simulations of the evolution of a magma ocean on Earth for different amounts of water and CO₂ 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₂ in our model. We compare our results with a pure water atmosphere to her results in Table 5.

Table 4. Input Parameters for Earth

Parameter	Initial value
	2–20 TO
	4000 K / 3000 K
	0.3
Atmosphere model	grey
modules used

Table 5. Results for Earth: Comparison between Elkins-Tanton (2008, Tab. 3, Earth [2000 km]) and (without CO₂)

Initial volatile mass	Dry case: 2 TO H2O / 0.4 TO CO2		Wet case: 20 TO H2O / 4 TO CO2
Model	Elkins-Tanton	(w/o CO2)	Elkins-Tanton	(w/o CO2)
Fraction of initial water content degassed into atmosphere [%]
	70	94	91	96
Final atmospheric pressure (sum of partial pressures of H2O and CO2) [bar]
	240	499	3150	5083
Time to reach 98% solidification, for and [Myr]
	0.06	0.3	2.4	3.0
Water content of liquids remaining at 98% solidification [wt%]
	1.5	1.7	5.3	9.5

→ GitHub (https://github.com/pbfeu/Trappist1_MagmOc/tree/master/Tab_Earth_Elkins-Tanton)

Especially for the water-rich case (20 TO water and 4 TO CO₂), the results are very similar and differ by less than a factor of 2 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₂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 those of Schaefer et al. (2016), Elkins-Tanton (2008), and Hamano et al. (2013). This shows that our model is able to simulate the magma ocean evolution of different kinds of rocky planets—inside and outside 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 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. For example, 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 ³He/²²Ne between the depleted mantle and a primitive deep reservoir. The D/H ratio of seawater indicates that the silicate Earth during the last global magma ocean was oxidized enough to outgas H₂O and CO₂ and not H₂ and CO (Pahlevan et al.,2019). The latter, reducing atmosphere mixture would have facilitated atmospheric escape of hydrogen leading to a higher D/H ratio. Similarly, the D/H ratio implies a short-lived steam atmosphere, that is, 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 to the potentially habitable planets of the TRAPPIST-1 system.

4. Results for the TRAPPIST-1 Planetary System

Since TRAPPIST-1 e, f, and g are orbiting TRAPPIST-1 in its habitable zone (e.g., Kasting et al.,1993a; Abe et al.,2011; Kopparapu et al.,2014; Catling et al.,2018; Turbet et al.,2018; Chen et al.,2019), we concentrated our efforts on simulating the magma ocean and volatile evolution for these three planets. Agol et al. (2021) and Grimm et al. (2018) provided new constraints on the radii and masses of the TRAPPIST-1 planets (Table 6). With these constraints, it is possible to give estimates on the density and composition of the planets (Barr et al.,2018; Dorn et al.,2018; Unterborn et al.,2018b).

Table 6. Physical and Run Parameters for TRAPPIST-1 e, f, and g

Planet	TRAPPIST-1 e	TRAPPIST-1 f	TRAPPIST-1 g
rp [R⊕]a	0.920	1.045	1.129
rc [R⊕]b	0.490	0.557	0.602
Mp [M⊕]a	0.692	1.039	1.321
a [au]a	0.0293	0.0385	0.0468
ec	0.005	0.01	0.002
Escape stop time [Myr]	253.2	129.4	76.3
Albedo	0.75
	1–100 TO
	4000 K
	0.3
Stellar age at t = 0	5 Myr
Atmospheric model	grey
modules used
Heating power [TW] for reference/extreme heating at beginning of simulation (t = 0):
Radiogenic heating
Tidal heating
Current water fraction estimates [wt%] (values in brackets show min and max water mass fraction):
Dorn et al. (2018)d
Barr et al. (2018)e	0–98	32–99	32–60
Unterborn et al. (2018b)d	0–1	2–11	9–23
Agol et al. (2021)f
This study, g

^aAgol et al. (2021).

^b

^cGrimm et al. (2018).

^dBased on Grimm et al. (2018) data.

^eBased on Gillon et al. (2017) data.

^fFor a core mass fraction of 25%.

^gSee Fig. 4, model by Noack et al. (2016), based on Agol et al. (2021) data.

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 TRAPPIST-1 e, f, and g based on the masses and radii as listed in the work of 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 transition 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), 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 color coding indicates the water content, and the size of the dots shows 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 shows 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 TRAPPIST-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 module , described in more detail in the work of Barnes et al. (2020). More precisely, uses energy- and diffusion-limited hydrogen and oxygen escape following the works of Luger and Barnes (2015) and Luger et al. (2015). Hydrogen is produced when the absorbed XUV flux of the host star photodissociates water, where we use an XUV absorption efficiency of (see also Table 6).⁹ For the TRAPPIST-1 system, we use the stellar evolution model contained in the module . 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 the work of 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 Eq. 2 of Kopparapu et al. (2013) and choose their values for the runaway greenhouse limit. This limit depends on the mass of the planet (Kopparapu et al.,2014). At this point, we assume that the pure steam atmosphere condenses out. Furthermore, the stratosphere desiccates and atmospheric escape of H₂O stops because water vapor no longer reaches the upper atmosphere to be photolyzed.

To set up the model for the TRAPPIST-1 simulations, we choose the following specific assumptions as input for the equations described in Section 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 can be found in Table 6.

4.1. Different evolution scenarios

We identify in this work three evolutionary scenarios for the TRAPPIST-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 vapor 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:

In scenario 1, the planet ends its magma ocean evolution with neither water vapor nor abiotically created O₂ in the atmosphere. We did not consider here the effect of N₂ or CO₂ gas left in the system, which could form the main constituents of the atmosphere in this scenario.

Example for scenario 2 (intermediate): Initial water content 5 TO for TRAPPIST-1 g

Here, an example of the intermediate water content case (scenario 2) is shown for TRAPPIST-1 g, as shown in the middle panel of Fig. 6.

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₂ (>100 bar). Since the magma ocean is already completely solidified at this point, O₂ cannot 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 vapor and O₂ atmosphere content in the bottom panel of Fig. 6.

In scenario 3, the planet ends its magma ocean evolution, when an extremely thick steam atmosphere (>10,000 bar) with about 1% free O₂ starts to condense out to form an ocean, leaving potentially a 100 bar O₂ atmosphere behind. However, we did not implement processes that could draw down O₂ into the mantle or into a water ocean after solidification.

4.2. 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 perturbations of other planets, and by setting the abundance of ⁴⁰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 3 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₂O₃. 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₂O, if the surface is cool enough, thus removing H₂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 vapor in the atmosphere, where it undergoes continuous atmospheric escape. This atmospheric escape re-enables outgassing from the magma ocean (evolution phase III). In other words, the prolongation of a thick magma ocean keeps water vapor 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₂O loss from the system. At the same time, the prolonged presence of a thick magma ocean efficiently prevents the buildup of O₂ in the atmosphere.

4.3. Overview of TRAPPIST-1 e, f, g results

A condensed overview of the results obtained from all our simulations is 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. The time when the atmospheric escape of H₂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₂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 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₂ 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₂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 scenario 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₂O stops because water vapor 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₂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₂O₃. 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₂ 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 dissolve and oxidize the melt than in the intermediate scenario 2.

5. 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) suggested 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's abundance of ⁴⁰K and fixed eccentricities (e = 0.1). The atmospheric composition is 100% H₂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, that is, with an initial water mass of more than 100 TO, especially for the planets f and g. As described in Section 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 buildup 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 lifetime (to at least 250 Myr for planet e) and to prevent abiotic oxygen buildup (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₂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, it is not possible in our simulations to sequester large amounts of water in the mantle without several thousand bar 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) suggested 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á and 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) suggested 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 bar 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 buildup assuming present Earth sinks and sources, though Harman et al. (2018) presented that lightning might be an effective way of reducing the oxygen abundance in such an atmosphere, as long as water vapor 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 vapor 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 TRAPPIST-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 buildup of oxygen in the atmosphere and therefore could facilitate the habitability of these planets. Heating rates of 1000 times Earth abundance of ⁴⁰K seem extremely high. However, ⁴⁰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 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 TRAPPIST-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. For example, Lichtenberg et al. (2019) showed that radiogenic heating of ²⁶Al leads to desiccation in rocky planetesimals and postulates that only planet building blocks poor in ²⁶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) provided dedicated formation scenarios for TRAPPIST-1, where they also propose that TRAPPIST-1 g likely formed very water-rich, TRAPPIST-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, that is, 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 viscosities 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) proposed that a thick TBL, including possibly even a lithosphere, may be neglected toward 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 TBL.

Because we aimed in this work for a versatile and numerically efficient magma ocean description, we neglected the “mush stage” toward the end of the magma ocean stage, following the argumentation of Debaille et al. (2009) and Hamano et al. (2013). As Lebrun et al. (2013) noted, 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 toward 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).

6. Conclusion

We investigated the magma ocean evolution of the planets TRAPPIST-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 (Section 2.1.4, Goldblatt et al.,2013). Only when the planet has lost all its atmospheric water via atmospheric escape within a few million years 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 bar of oxygen can build up in the atmosphere. At the time of mantle solidification, the atmosphere will consist of free oxygen and potentially water vapor. 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 . The thick magma ocean prevents oxygen buildup, and 90% of H₂O remains dissolved in the magma ocean. The whole planet is, however, after 250 Myr in a state where water vapor 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) suggests 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 ²⁶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 is 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 was the case for 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 (Noack et al.,2016). Furthermore, it could build up several hundred bar 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, is a powerful tool inside the 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.

7. Outlook

Even though our model successfully reproduces previous results for magma ocean evolution of Earth and the super-Earth GJ 1132b (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₂ into the model. As the example of Venus shows, CO₂ can be an important constituent of atmospheres of rocky planets. CO₂ was also a major constituent of early Earth's atmosphere (e.g., Kasting, 1993; Deng et al.,2020). The amount of CO₂ present in the atmosphere determines whether condensation of water to an ocean can occur. The higher the CO₂ 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₂ into this model. Our model so far can only simulate water composition up to 100 TO. Since TRAPPIST-1 g has likely a higher water content than 100 TO, future work must reformulate the equation of state in a future version of to coherently tackle extremely wet scenarios.

Furthermore, in this work, we assumed bulk silicate Earth composition, that is, an oxidizing mantle composition. However, even for early Earth this assumption has come under scrutiny (Kasting et al.,1993b; Gaillard and Scaillet, 2014; Ehlmann et al.,2016; Fegley et al.,2016; Schaefer et al.,2016; Armstrong et al.,2019; Del Genio et al.,2020). 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₂ instead of CO₂ and H₂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 has 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 TRAPPIST-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 . Due to the close proximity of planets orbiting M dwarfs to their host star, for example, 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 and Laughlin (2019) note that obliquity tides in a tight planetary system like TRAPPIST-1 could create up to 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 buildup 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 a billions-of-years 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 is 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 and Bercovici, 2013; Garcia et al., unpublished data). 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₂-dominated atmosphere instead of an oxidizing atmosphere as assumed in this study (Tsiaras et al.,2016; Hammond and 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, f, and 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 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 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, provides the ideal link between formation models such as that of Miguel et al. (2020), interior structure studies like those of 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 those of Kislyakova et al. (2017) and Millholland and Laughlin (2019). is thus a valuable tool to explore the diversity of rocky worlds in the Solar System and beyond.

Author Disclosure Statement

No competing financial interests exist.

Associate Editor: Sara Seager

Appendix

Table A1. Values for Parameters Used in

Symbol	Parameter	Value
a	Mantle bulk density
cp	Silicate heat capacity
	Heat of silicate fusion
	Thermal expansion coefficient
	Water part. coeff. melt-solid	0.01

If not otherwise indicated, values are taken from Schaefer et al. (2016).

^aLebrun et al. (2013).