# A cool runaway greenhouse without surface magma ocean

Franck Selsis, Jérémy Leconte, Martin Turbet, Guillaume Chaverot, Émeline Bolmont

► To cite this version:

Franck Selsis, Jérémy Leconte, Martin Turbet, Guillaume Chaverot, Émeline Bolmont. A cool runaway greenhouse without surface magma ocean. Nature, 2023, 620 (7973), pp.287-291. 10.1038/s41586-023-06258-3. hal-04288510

HAL Id: hal-04288510

https://hal.science/hal-04288510v1

Submitted on 24 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A cool runaway greenhouse without surface magma ocean

Franck Selsis¹*, Jérémy Leconte¹, Martin Turbet², Guillaume Chaverot³ and Émeline Bolmont³,⁴

¹*Laboratoire d'astrophysique de Bordeaux, Univ. Bordeaux, CNRS, Pessac, France.

²Laboratoire de Météorologie Dynamique/IPSL, CNRS, Sorbonne Université, École Normale Supérieure, PSL Research University, École Polytechnique, 75005 Paris, France.

³Observatoire Astronomique de l'Université de Genève, Chemin Pegasi 51, Sauverny 1290, Switzerland.

⁴Center for Life in the Universe, Faculty of Science, University of Geneva, Geneva, Switzerland.

*Corresponding author(s). E-mail(s): franck.selsis@u-bordeaux.fr; Contributing authors: jeremy.leconte@u-bordeaux.fr; martin.turbet@lmd.ipsl.fr; guillaume.chaverot@unige.ch; emeline.bolmont@unige.ch;

## Abstract

Water vapour atmospheres with content equivalent to the Earth's oceans, resulting from impacts[1] or a high insolation[2, 3], were found to yield a surface magma ocean[4, 5]. This was, however, a consequence of assuming a fully convective structure[2, 3, 6, 7, 4, 5, 8, 9, 10, 11]. Here we report, using a consistent climate model, that pure steam atmospheres are commonly shaped by radiative layers, making their thermal structure strongly dependent on the stellar spectrum and internal heat flow. The surface is cooler when an adiabatic profile is not imposed: melting Earth's crust requires an insolation several times higher than today, which will not happen during the main-sequence of the Sun. Venus' surface can solidify before the steam atmosphere escapes, which is opposite to previous works[4, 5]. Around the reddest stars ($T_{eff} < 3000$ K), surface magma oceans cannot form by stellar forcing alone, whatever the water content. These findings
2 *Cool runaway greenhouse*

affect observable signatures of steam atmospheres and exoplanet mass-
radius relationships, drastically changing current constraints on the water
content of Trappist-1 planets. Unlike adiabatic structures, radiative-
convective profiles are sensitive to opacities. New measurements of
poorly constrained high-pressure opacities, in particular far from the
H₂O absorption bands, are thus necessary to refine models of steam
atmospheres, which are important stages in terrestrial planet evolution.

**Keywords:** Runaway greenhouse, magma ocean, Earth, Venus

# 1 Main article

Modeling the steam atmosphere resulting from the vaporization of the Earth oceans (270 bars) is challenging. If a few consistent 1D [12] and 3D simulations have been achieved up to 10-30 bars[13], most models relied on the hypothesis that steam atmospheres are fully convective below their photosphere. This assumption allowed to 1) building the thermal structure from a surface temperature Tsurf using dry and saturated adiabats and 2) determining the corresponding incoming stellar radiation (ISR) by *inverse climate modeling*, i.e. by insuring a radiative balance at the top of the atmosphere based on the computed outgoing thermal radiation (OTR) and Bond albedo A. In this framework, the OTR decreases slowly from the Nakajima-Simpson limit[14] of ~280 – 290 Wm⁻² (depending on the H₂O mixing ratio[15], clouds[13] and gravity[16]) to a value of ~274 Wm⁻² over a broad range of surface temperatures: from about 300 K to 1700 K. Above 1700 K the OTR increases again with Tsurf. Therefore, once the absorbed flux exceeds the runaway threshold, the surface rapidly evolves from *habitable* conditions into a magma ocean. Conversely, when a hot and cooling steam atmosphere above a magma ocean reaches an OTR of 274 Wm⁻² and a Tsurf of 1700 K, it enters a fast cooling phase, leading to solidification of the crust and condensation into an ocean.

However, as shown in the section *Most steam atmospheres are not fully convective* (2.1), balancing the radiative fluxes at the top of a convective atmosphere, as done in *inverse climate modeling*, does not insure that all layers are at thermal equilibrium as they should. Additionnaly, we found that they are not in most of the parameter space we explored: surface pressure psurf between 0.1 and 270 bars, stellar effective temperature Teff between 2600 and 7000 K, OTR between 274 and 10⁵ Wm⁻² and internal heat flux between 0 and 1000 Wm⁻². It is therefore necessary to self-consistently calculate the atmospheric thermal profile, without imposing convection, following the tracks of pioneering works on runaway atmospheres [17] and formation/condensation of post-accretion atmospheres [18]. This is required to reassess results obtained assuming fully convective profiles on the divergent fate of Venus and the Earth[3, 4], the circumstellar region where surface liquid water can be stable[6, 9] and the occurrence of surface magma oceans[4, 5, 19, 20, 21].
To calculate the vertical structure consistently with the incoming stellar radiation and underneath heat flow, we developed an evolution code within the Exo-k library [22] specifically to handle the coexistence of very large evolution timescales at high depth with very short ones in the upper atmosphere. With Exo-k (described in details in the Methods), we computed atmospheric profiles in equilibrium with different stellar fluxes/spectra and compared them in Figure 1 with convective profiles that share the same bolometric emission. We considered a sample of stars from the ultracool dwarf Trappist-1 (hereafter T-1, $T_{\text{eff}}$ = 2600 K) to a F1-type star ($T_{\text{eff}}$ = 7000 K). Profiles with an OTR of 274 Wm$^{-2}$ are at the transition between the hot regime where the OTR is sensitive to $T_{\text{surf}}$ and the unstable plateau of the runaway or condensation (depending whether the evolution proceeds toward increasing or decreasing surface temperature). In the full-convection scenario, the planet exits runaway mostly when the dry troposphere reaches and warms the photosphere, increasing the OTR [23]. For a given OTR, all profiles cross within the photospheric region, where the temperature controls the emission, and which lies within saturated layers for an OTR of 274 Wm$^{-2}$.

Here, we identify two processes keeping steam atmospheres mostly stable against convection and allowing the atmosphere to equilibrate at lower surface temperatures. First, a radiative zone can form above the surface depending on the surface pressure and the internal heat flux and can potentially extend to the whole atmosphere, depending on the stellar type and incoming flux (see Extended Data Figure 2). This results from the vanishing stellar flux at depth, coupled with an upward radiative transport whose efficiency increases with temperature. Indeed, the thermal gradient needed to radiate a given thermal flux decreases with increasing temperature at fixed opacity, as the radiative diffusivity scales as $dB/dT \sim T^3$, where B is Planck's spectral radiance. This effect is well known in hot Jupiter atmospheres where radiative layers extend to pressures as high as 100-1000 bars depending on instellation and heat flow [24]. Second, the absorption of stellar light directly warms the middle atmosphere so that the temperature of the photosphere is not set only by either the saturation pressure $p_{\text{sat}}(T)$ or an adiabat. This reduces further the contrast between the photospheric and surface temperature for a given OTR. At higher ISR, condensation is totally inhibited making some profiles entirely radiative as can be seen in Extended Data Figure 2. The departure from purely convective profiles is not restricted to cases with high pressures blocking the penetration of the stellar radiation toward the surface. Indeed, the intensity and spectral distribution of the stellar flux influence the thermal structure of steam atmospheres whatever the surface pressure as illustrated in Extended Data Figure 6.

Compared with computed equilibrium profiles, convective structures dramatically overestimate $T_{\text{surf}}$ by 400 K (Sun) to 1100 K (T-1) at the exit of the runaway greenhouse (OTR=274 Wm$^{-2}$), and by 1500 K (Sun) to 2300 K (T-1) for an OTR of 10,000 Wm$^{-2}$ for the Earth's surface water reservoir (270 bars).
4 Cool runaway greenhouse

<figure><img src="image_2.png" /><figcaption>**Figure 1** Steam atmosphere thermal profiles. Each group of five profiles shares the same bolometric emission (274 and 10,000 Wm⁻², thin and thick lines respectively). The profiles computed for different stellar spectra are compared with convective profiles of same emission. For each computed profile, the Bond albedo and instellation are given (S<sub>eff</sub> = 1 corresponds to Earth's present insolation: 341.5 Wm⁻²). Dashed lines indicate dry convection.</figcaption></figure>

When using convective profiles, runaway extends to surface-melting tem-
peratures [9, 16, 8, 25, 19, 26, 27, 23]. However, in all converged cases and as
seen in Figure 2, planetary cooling becomes sensitive to T<sub>surf</sub> at lower temper-
ature: under solar forcing, the runaway phase ends at a surface temperature
around 1250 K. For a typical mantle composition, the low-pressure solidus
and liquidus are found at 1400 and 1950 K, respectively [28]. At 1620 K
the melt fraction reaches 40% and dominates the rheology [29]. This melting
temperature requires S<sub>eff</sub> > 7 under solar forcing.

Neither the Earth or Venus have been or will be exposed to such high ISR
during the main sequence of the Sun. As a result, a surface magma ocean
can be sustained only as long as the interior is able to output a sufficient
geothermal flux, either from radiogenic heating or leftover accretion energy. We
thus expect that surface magma oceans solidify more rapidly than predicted
by previous calculations. With the insolation of Venus 4.5 Gyrs ago, the crust
would solidify for internal heat flows below 2, 10 and 100 Wm⁻²
for p<sub>surf</sub> = 270, 100 and 10 bars of steam respectively. The post-accretion
atmosphere of Venus must have consisted in other species than H₂O so we
<figure><img src="image_3.png" /><figcaption>**Figure 2** Outgoing thermal emission as a function of surface temperature. The dashed-dotted line is the classic runaway curve obtained with fully convective profiles, independently from the stellar radiation spectrum and internal flux. With self-consistent profiles, this function depends on the type of the host star and the internal flux $\phi_{\text{int}}$. Dashed curves are obtained with Earth geothermal flux. For T-1 and $\phi_{\text{int}} < 0.01 \text{ Wm}^{-2}$, the runaway greenhouse can end without full vaporization of the ocean. Low-pressure liquidus, solidus and melting (melt fraction = 40%) temperatures are indicated for an Earth mantle composition. Surface pressure is 270 bars except in the blue area where $T_{\text{surf}}$ is below the critical temperature of H$_2$O and $p_{\text{surf}}$ is given by the saturation vapor pressure of H$_2$O. These curves are based on the profiles shown in Extended Data Figure 2</figcaption></figure>

should not draw conclusion from these numbers, but the standard scenarios for Venus early history will have to be revisited with consistent profiles. In particular, surface solidification may not require the loss of water through atmospheric escape in contrast with previous findings[4, 5].

We find a stable steam atmosphere for a solar flux as low as $S_{\text{eff}} = 0.98$, a consequence of a known bistability [8, 13]: a higher ISR is required to vaporize an ocean ($S_{\text{eff}} = 1.11$ [30]) than to condense a steam atmosphere ($S_{\text{eff}} = 0.92$ [13]). The difference between the 3D simulation including clouds and our cloud-free 1D model is attributed to the greenhouse warming by nightside clouds.

Profiles computed for the star T-1 are fully radiative except below the saturated layers, where dry convection is locally triggered by rain evaporation cooling. The steam atmosphere is more opaque to its red radiation and only 0.001% of the incident flux reaches 100 bars, resulting in low surface temperatures. The overestimation of $T_{\text{surf}}$ by convective profiles is considerable: 1200 K for $S_{\text{eff}} = 0.8$, just above the runaway threshold, and 2700 K for $S_{\text{eff}} = 58$. No
6 Cool runaway greenhouse

<figure><img src="image_4.png" /></figure>

**Figure 3** Spectra of Trappist-1b with a steam atmosphere. Left: Top-of-the-atmosphere emission and reflection spectra of T-1b with a 270-bars steam atmosphere. Spectra are computed with a converged profile (blue) and with a $T_{\text{surf}} = 2410$ K convective profile satisfying OLR = ISR(1-A) (red). Right: transmission spectra given as apparent radius variations. Converged profiles are in blue and given for either 1 or 3 Earth oceans and $\phi_{\text{int}} = 0$ or $5 \text{ Wm}^{-2}$ (a range of estimated tidal dissipations).

planet is close enough to T-1 to be exposed to an instellation large enough to melt the surface under a steam atmosphere, even during the pre-main-sequence (PMS) phase. This has important consequences as it weakens the coupling between reservoirs in the atmosphere and the mantle. Losses to space during PMS can exhaust the steam atmosphere while keeping intact the internal reservoir. In previous models assuming a steam atmosphere in equilibrium with a surface magma ocean, atmospheric escape draws water from the mantle and can desiccate the whole planet [4].

This work shows that ultracool dwarfs ($T_{\text{eff}} \le 2600$ K) can theoretically host a hot liquid ocean under a dense steam atmosphere despite an ISR beyond the runaway limit. However, this case remains somehow academic because this requires a high surface pressure (>100 bars), which produces a thick blanket for the internal heat. A fraction of Earth's geothermal flux is then sufficient to heat the surface above the critical temperature of H$_2$O (see Figure 2 and Extended Data Figure 6).

The interpretation of planetary spectra and transit lightcurves differs whether one assumes convective structures or calculates actual profiles. This applies to all exoplanets beyond the runaway limit that could possess a significant water reservoir but let us illustrate this point with the three innermost planets of Trappist-1, which offer incomparable observation conditions for planets with this size and instellation. Occultations of planets b and c and transits of all T-1 planets are observed during the first cycle of the James Webb Space Telescope.

At a given ISR, the higher temperatures found with convective profiles at pressures higher than ~100 mbar are responsible for the two main differences, shown in Figure 3 for the conditions of T-1b. First, it produces a strong emission between the water bands below 2 µm and in particular, in the visible, where it exceeds the flux reflected by the planet. With converged profiles this shortwave emission is basically absent while the stellar flux absorbed in the
upper atmosphere is re-emitted at low temperature, enhancing the longwave emission. Due to larger scale heights at higher temperature, the altitude of the highest opaque layers is increased in convective profiles.

As a consequence, we can see in Figure 3 that convective profiles shift the transmission spectra by several hundreds kilometers for the same surface pressure. This overestimation of the runaway greenhouse radius inflation effect [27] by the steam atmosphere yielded an underestimation of the maximum water content of T-1b, c and d [31, 32, 33, 34]. With inverse climate modeling, the internal heat flux only affects the value of Tsurf corresponding to top-of-the-atmosphere radiative balance. Outside the runaway OTR plateau, only internal heat flows reaching a significant fraction of the instellation affect atmospheric and spectral properties. With actual profiles, the deep atmosphere that is deprived of stellar light is very sensitive to the internal heat flux, which becomes a parameter controlling the planet radius. For T-1b, with one and three Earth oceans, the convective hypothesis causes an overestimation of the radius, in all observation bands, by 300 km (4.2%) and 400 km (5.6%), respectively. Accounting for an internal heat flow due to tidal dissipation and estimated at 5 Wm⁻²[35, 36, 37, 38, 39] adds ~100 km (1.4%) to the radius. The current 1-σ relative precision on the radius of T-1b is 1%[32].

Our calculations thus have important consequences for the interpretation of the exoplanet radius gap [40] and of the true nature of sub-Neptune size planets [41, 42]. Indeed, we show that water-rich atmospheres are significantly less extended than previously calculated[31, 43, 44], which challenges the interpretation of the bulk population of sub-Neptune size planets as being water worlds.

*Inverse climate modeling* based on the full convection assumption is independent of the H₂O opacity at high pressure, sensitive only to the Bond albedo and OTR controled by low-pressure layers. Profiles consistently at equilibrium with the incoming stellar radiation and internal heat flow prove dependent on opacities at all pressures. These opacities control conditions in the deep atmosphere and at the surface (see 2.6 in Methods). The current lack of experimental constraints, in particular for H₂O high-pressure far-wing opacities, thus becomes a critical issue. Our conclusions that radiative layers must be considered and that magma oceans do not result from stellar forcing alone around solar and subsolar stars are robust against orders of magnitude changes in these high-pressures opacities. This is mainly due to the strong temperature feedback that makes radiative transport more efficient at high temperature. Nevertheless, computed thermal profiles remain quantitatively uncertain and we hope this study will motivate measurements of these data, needed to better model the properties and evolution of a broad range of planets, including our own.
# 2 Methods

## 2.1 Most steam atmospheres are not fully convective

Although a convective structure is well established at the onset of the runaway[30], a simple test demonstrates that convection cannot be maintained by stellar radiation beyond the runaway threshold. At equilibrium, net (↑ − ↓) energy fluxes must satisfy at each level:

$$F_*(p) + F_E(p) + F_C(p) = \phi_{\text{Int}}$$

where *p* is the pressure level, $F_*$ is the net stellar flux, $F_E$ is the net thermal emission flux and $\phi_{\text{int}}$ is the internal heat flux. As convection can only transport heat upward, this convention implies that the convective flux, $F_C(p)$, must be positive in a physically consistent model. Extended Data Figure 1.a shows the net radiative thermal and stellar fluxes computed for a convective steam atmosphere with $T_{\text{surf}} = 2000$ K and for an ISR insuring a null net radiative flux at the top of the atmosphere. The profile constructed with adiabats is not at local equilibrium and balancing the net radiative flux would require an unphysical downward convective flux. The solar flux reaching the lower atmosphere is insufficient to compensate for the strong thermal cooling resulting from the assumed convective gradient: only 12, 0.3 and 0.025 Wm⁻² of solar radiation reach a depth of 10 and 100 and 270 bars, respectively. This net flux inconsistency, illustrated here on a specific profile, is systematically observed for $T_{\text{surf}} > 700$K for $P_{\text{surf}}=270$ bars, a solar spectrum, and Earth gravity.

## 2.2 The model

To address this issue we use the tools from the Exo-k evolution package. Panels b and c of Extended Data Figure 1 show the evolution from an initial convective profile at top-of-the-atmosphere radiative equilibrium towards a structure in complete equilibrium (at all pressure levels). In the converged case, $T_{\text{surf}}$ is no less than 700 K cooler than in its assumed convective counterpart. The advantage of Exo-k is that it uses several physical-acceleration techniques able to accurately reach a radiative-convective equilibrium within short CPU times.

The Exo-k library[22] includes a class dedicated to model self consistently the evolution of planetary atmospheres. It accounts for radiation, dry convection, turbulent diffusion, and moist processes (condensation, moist convection, and precipitations) for any number of species. We provide a brief description of the salient features of the model below. More details can be found in the online documentation at http://perso.astrophy.u-bordeaux.fr/~jleconte/exo-k-doc/index.html.

### Radiative transfer

The radiative transfer is computed using a 2-stream approximation[45]. The opacities of the radiative layers are computed using the Exo-k opacity library
and can use all the types of opacity sources available through it (correlated-k, cross sections, collision-induced absorptions, Rayleigh scattering, aerosol Mie scattering, etc.). The radiative transfer module has already been tested extensively for runaway greenhouse atmospheres[15].

## Dry convective adjustment

The dry convective adjustment scheme looks for convectively unstable regions in the atmosphere where the virtual potential temperature decreases upward. The potential temperature and composition of unstable layers are fully mixed over a single timestep while conserving total mass, tracer mass, and enthalpy.

## Moist convective adjustment and condensation

Our scheme to handle moist convection is similar to the one used in previous models for runaway greenhouse atmospheres[30]. For water vapor dominated atmospheres as the ones treated here, the pseudo-adiabat reduces to the saturation vapor-pressure curve of pure steam. Condensation of vapor can also happen in a layer in absence of convective processes, when there is diabatic cooling, for example. To account for this, at each timestep, this parametrization brings any supersaturated layer back to vapor equilibrium. Subsaturated layers can also evaporate condensates if they are present. This saturation adjustment is performed iteratively until equilibrium conditions are found at

<figure><img src="image_5.png" /></figure>

**Extended Data Figure 1** Converged versus fully convective profiles. The steam atmosphere is modeled here for a solar flux of 445.7 Wm⁻² (S<sub>eff</sub> = 1.3), no internal heat flux, a 1g gravity and a 270 bars surface pressure (1 vaporized Earth ocean). **Panel a:** Net fluxes computed for a “convective” P-T profile used as the initial structure in panels c and d. **Panel b:** Net fluxes at equilibrium. Blue (resp. red) area indicate a downward (resp. upward) net radiative flux. Convection transporting energy upward can only balance a downward radiative flux. Red area therefore indicate a departure from thermal equilibrium. **Panel c)** Evolution from the initial dashed “convective” profile towards the converged state, computed with the Exo_k suite[22] evolution package. It takes more than 40,000 yrs to reach the converged state and more than 280 millions steps and 40 hrs of CPU time. **Panel d)** Evolution computed with an acceleration mode (see Methods) of Exo_k, in 165,000 steps and less than 1 min of CPU time. In the acceleration mode, iterations and intermediate profiles do not correspond to physical times and structures.
constant moist enthalpy.

## Rains

In the current version, all condensates precipitate instantaneously upon condensation, always leaving a cloud-free atmosphere. However, as we are interested in deep atmospheres where precipitations are unlikely to reach the surface (if it exists), we implemented a simple scheme to reevaporate falling precipitations. We start from the model top layer and collect all precipitations downward. Whenever an unsaturated layer is met, a fraction f of the condensate that would need to be evaporated to saturate the layer is effectively evaporated. Precipitation are fully evaporated when they reach a layer where the temperature is above the boiling temperature of the falling species. All remaining precipitations when the surface layer is reached, if any, are added to the surface reservoir.

In our study, precipitations are always reevaporated just below the saturated convection region and force a small, dry convective region there. This region experiences some level of time variability that is reminiscent of the results of 3D convection resolving models of hothouse Earth-like climates[46].

## Time integration

Once we are able to compute heating rates in each layer ($\hat{H} \equiv \{H_n\}_{n=1,N_{\text{lay}}}$; hatted symbols refer to vectors that denote the values of a given quantity in all the $N_{\text{lay}}$ atmospheric layers at the same time), the temperature in each layer ($\hat{T}$) is integrated in time until thermal equilibrium is reached following

$$c_p \frac{\partial \hat{T}}{\partial t} = \hat{H}(\hat{T}), \qquad (1)$$

where $c_p$ is the specific heat capacity of the gas. In our case, only the equilibrium state of the atmosphere, where $\hat{H}(\hat{T}) = \hat{0}$, is of interest. Therefore, a particular attention has been devoted to the development of schemes to accelerate the convergence toward such steady-states.

**Adaptive timestep.** As radiative transfer is usually the most expensive part of a 1D model, extra care has been taken to compute it as seldom as possible. As a first step in that direction, Exo_k has an adaptive timestep that is based on the radiative timescale of the atmosphere. This radiative timescale is computed as follows. We start by saying that, around an arbitrary thermal state, say $\hat{T}_K$, the radiative heating rates can be linearized through

$$\hat{H}^{\text{rad}}(\hat{T}_K + \delta\hat{T}) = \hat{H}^{\text{rad}}(\hat{T}_K) + \mathbb{K} \cdot \delta\hat{T}, \qquad (2)$$

where the kernel $\mathbb{K}$ is the Jacobian matrix of the heating rates (of dimension $N_{\text{lay}} \times N_{\text{lay}}$) and $\delta\hat{T}$ is a small temperature perturbation.
This kernel matrix tells us how heating in a layer is related to a temperature change anywhere in our atmosphere. The biggest terms are the diagonal ones, which are usually negative. This shows that any layer whose temperature is increased will tend to emit more and cool in response. The terms directly above and below the diagonal – the coupling terms between adjacent layers – are usually positive but smaller in magnitude. Other off-diagonal terms, which are long range couplings, are much smaller, especially in optically thick atmospheres.

Now, let us consider the radiative evolution of the system around radiative equilibrium — meaning $\hat{H}^{\text{rad}}(\hat{T}_{\mathbb{K}}) = \hat{0}$. Combining the above equations yields the following evolution for the temperature perturbation:

$$c_p \frac{\partial}{\partial t} \delta \hat{T} = \mathbb{K} \cdot \delta \hat{T}. \qquad (3)$$

Keeping only the most important terms (the diagonal ones), the layers decouple and the solution reads

$$\delta \hat{T} = \delta \hat{T}_0 \cdot e^{-\text{Diag}(\mathbb{K})t/c_p}, \qquad (4)$$

where Diag(K) is the vector formed by all the diagonal terms in K. So we see that in each layer, an initial perturbation, $\delta \hat{T}_0$, will be radiated away on a timescale equal to $\hat{t} = c_p/\text{Diag}(K)$ that differs from layer to layer. In the baseline evolution, the smallest radiative timescale in our atmosphere, $\tau_{\text{rad}} = \min(\hat{t})$, is used as timestep to ensure that the radiative evolution is well sampled throughout the atmosphere. We find that this condition is sufficient to ensure a stable evolution in most cases, but the user can always use smaller timesteps by specifying a factor $\alpha$ so that $\Delta t = \alpha \tau_{\text{rad}}$.

**Using the Jacobian to compute fluxes.** Another advantage of having computed the Jacobian is that, as long as the current state of the atmosphere is sufficiently close to the last state for which we computed K, Eq. 2 can be used to compute the radiative heating rates extremely rapidly through a simple matrix multiplication. In practice, we use this method as long as the maximum difference between the current temperature and $\hat{T}_{\mathbb{K}}$ is below some user defined threshold. Otherwise, the full radiative transfer is computed. This is especially efficient when the atmosphere is approaching equilibrium and requires a large number of small temperature increments to equilibrate its optically thick parts while the upper atmosphere requires small timesteps to remain stable.

### Convergence acceleration

Despite the approach described above, the convergence time of an atmosphere remains prohibitive in many cases. This is usually due to the fact that the upper atmosphere has a short radiative timescale that requires small timesteps for stability while the deep atmosphere is very opaque and has a large thermal inertia.
12

To circumvent these issues, the library offers several ways to accelerate
the convergence toward an equilibrium state, as illustrated in Extended Data
Figure 1. It should be stressed, however, that the evolution trajectory followed
by the atmosphere toward this state cannot be regarded as a temporal evolu-
tion, even tough we will refer to every evolution step as a *timestep* for sake of
brievity.

For each timestep, we first compute the heating rates for all the processes
in the regular way described above over a duration Δt. We then identify
radiative zones as groups of adjacent layers where the heating rates of all the
processes are zero, except for the radiation. A radiative zone can be composed
of a single layer. The remaining layers of the atmosphere are grouped in
stacks of adjacent layers that we will call convective zones for convenience,
even though the energy exchange in these zones might be due to formation or
reeevaporation of rains.

**Acceleration in radiative zones.** We start by defining a base timescale, $\tau_b$, as the smallest radiative timescale of a radiative layer in the atmosphere. Then, for any radiative layer n, the heating rate is multiplied by $\tau_n/\tau_b$ before the timestepping is performed.

In a purely radiative atmosphere this would be equivalent to advancing
every layer independently using its own radiative timescale. In practice, this
still allows deep radiative zones with very long timescales to converge in as
many timesteps as their counterparts in the upper atmosphere.

**Acceleration in convective zones.** In convective zones, the layers in the
zone are not independent as they directly exchange energy in a conservative
way: except for radiation, energy that is taken in a part of the zone is redis-
tributed in some other part. Using the scheme above in such zones completely
alters the balance.

To understand why, let us consider a simple 2-layer convective zone where
radiation heats the base layer and cools the upper one. When the layer becomes
unstable, dry convection carries the surplus of heat from the bottom to the top
layer. Equilibrium is reached when the convective energy flux equals both the
net radiative heating of the base layers and the net radiative cooling of the top
layer. Now, one can see that if the radiative heating rates are multiplied by two
independent factors, our equilibrium solution is not an equilibrium anymore
and the system will reach a new unphysical equilibrium.

To avoid this problem, we treat an entire convective zone as one layer
during the acceleration part. We compute the average radiative heating rate
in the zone by computing the net radiative fluxes at the top and bottom of the
convective zone and dividing it by its total mass. Then we compute a radiative
timescale for the whole zone, $\tau_{av}$. For the moment, we use the smallest radiative
timescale of any layer inside the zone to be conservative. Finally, in addition
to the heating rate computed in the initial step, each layer receives the same
average radiative heating rate multiplied by $\tau_{av}/\tau_{rad}$.
<figure><img src="image_6.png" /><figcaption>Extended Data Figure 2 Steam atmosphere converged P-T profiles as a function of instellation and stellar type. Dashed lines indicate dry convection. Adiabatic profiles satisfying top-of-the-atmosphere radiative balance for the minimum and maximum fluxes are in gray.</figcaption></figure>

## 2.3 Simulations

Unless indicated otherwise, simulations presented in this study were run for a pure H₂O atmosphere, a surface gravity of 1g and a surface pressure of 270 bars corresponding to vaporized Earth oceans. Our radiative transfer uses the H₂O linelists from HITRAN2016[47] and HITEMP2010[48]. Above 700 K, we use HITEMP over the whole wavelength range. Below 700 K, we use HITRAN as much as possible because of considerable CPU time required to compute high-resolution spectra with HITEMP but mainly because the MT_CKD H₂O-H₂O continuum parametrization (version 4.0.1[49]) has been determined consistently with HITRAN (a continuum being basically the difference between the opacity measured and the opacity calculated with a given linelist when truncating the line profiles). Wavelength domains have been selected based on a comparison of opacities computed with both linelists[8]: we use HITRAN at wavelengths above 3µm and HITEMP below 3µm. Opacities are treated with a correlated-k approach in 69 spectral bands from 0.28 to 1000 µm[13]. The model includes condensation and transport/evaporation of condensates but not their radiative effects. For the Sun, Trappist-1 (hereafter, T-1) and Proxima Centauri we use synthetic spectra[50] while we approximate other stars (M3, K5 and F1) with a blackbody.

Extended Data Figure 2 shows converged profiles at different ISR for the Sun
14 Cool runaway greenhouse

and other stars and Extended Data Figure 6 presents converged profiles for
various surface pressures and internal heat fluxes.

2.4 Comparison with a GCM

To test the results of Exo_k, we compared them with those of the Generic PCM (Planetary Climate Model), previously known as the LMD Generic GCM[51, 13]. In addition to an existing simulation done for the Sun[13], we performed new runs with the same parameters but colder stellar spectra and synchronous rotation. Convergence of 3D simulations was verified by starting from both hot and cold initial states[13]. Both 1D and 3D simulations were done for pSurf = 10 bars as thicker atmospheres require prohibitively large CPU times in 3D. The comparison is presented in Extended Data Figure 3 and shows a good agreement. First of all, it confirms our main findings: the thermal structure - does depend on the spectral distribution of the stellar radiation,

- is radiative in the lower atmosphere,

- and more sub-adiabatic for cooler stars.

Exo_k and Generic PCM are in close agreement for the stellar heating rates.

<figure><img src="image_7.png" /><figcaption>**Extended Data Figure 3 1D/3D comparison.** Panel a, b and c) Thermal profiles obtained with Exo_k (solid blue lines) and the 3D Generic PCM (dashed red lines) for Trappist-1, Proxima and a M3 star. For the Generic PCM, the spatial and temporal average is shown as well as the range of variations (red area). Panel d) Stellar heating rates with Exo_k (solid lines) and with the Generic PCM (dashed lines). The 10 bars atmosphere consists of 95% of H₂O and 5% of N₂ and the instellation is 500 Wm⁻² (Sₐff = 1.42) in all cases. The opacities used in both models are the same, differences come from circulation and cloud radiative effects.</figcaption></figure>

This match is allowed by the absence of clouds on the dayside[13]. Thermal profiles show an excellent match at *p* < 10⁵ Pa, thanks to large-scale circulation maintaining a moist adiabat on the dayside at 10 mbar despite the absence of condensation. The Generic PCM is systematically hotter by ~ 20 K at the surface. This small departure is attributed to the greenhouse warming of nightside clouds[13]. Neither model finds convective layers at *p* > 10 mbar.

Note that these 1D and 3D simulations were performed with the MT_CKD
3.5 continuum. The 3D runs were done before the release of the 3.6 and 4.0.1
<figure><img src="image_8.png" /><figcaption>**Extended Data Figure 4** Influence of the $c_p$ value on the converged T profiles.
The black curves are the nominal P-T profiles, obtained with a $c_p$ set by iterations to its value at the mean temperature in the dry convection layers. To show the sensitivity to the $c_p$ value we used the lowest and highest temperatures ($T_{min}$ and $T_{max}$) found in these dry convective layers and computed the blue profile with $c_p(T_{min})$ and the red profile with $c_p(T_{max})$. Dashed lines indicate dry convection.</figcaption></figure>

versions. While they are based on MT_CKD 4.0.1 everywhere else in the study, the 1D runs shown in Extended Data Figure 3 are deliberately run with the same opacities as in the 3D model so that differences only result from 3D dynamics and clouds. By enabling the use of a 1D model, this comparison initiated the present study.

## 2.5 Spectra

Exo-k can produce emission/reflection and transmission spectra like the synthetic spectra of the exoplanet Trappist-1b shown in figure Figure 3. We used the same spectral bands to compute the atmospheric structure and the spectra, as it is sufficient to display the effects described in the study. Spectral and profile calculations both account for the variation of gravity with altitude.

## 2.6 Sensitivity to opacity

Radiative-convective profiles computed in this study depend on H₂O opacities that remain uncertain. In particular the radiative transport in optically thick layers is sensitive to far-wing collisional absorption in atmospheric windows. The continuum that controls this absorption is extrapolated from measurements made at temperatures and pressures far below those prevailing in the depths of a vaporized ocean. In Extended Data Figure 5 we show how the calculated profiles are changed by the most recent update of the MT_CKD H₂O self continuum. Compared to version 3.5, the temperature dependency in the version 4.0.1 was modified to account for new measurements and
<figure><img src="image_9.png" /><figcaption>Extended Data Figure 5 Sensitivity to a significant change in the H₂O-H₂O continuum. Profiles are computed for a solar spectrum and three different insolations with versions 3.5 and 4.0.1 of the MT_CKD continuum. Dashed lines indicate dry convection.</figcaption></figure>

parametrizations[52]. This comparison does not constitute in any way an estimate of the uncertainty on these profiles but they illustrate two points. First, despite orders of magnitude change in absorption in some windows, the strong temperature feedback on radiative diffusivity, variations of $T_{surf}$ do not exceed 50 K. Second, it however reminds us that new data will inevitably alter the quantitative results.

## 2.7 Sensitivity to non-ideal behaviors of the equation of state

In our simulations (either 1D or 3D), the gas specific heat capacity $c_p$ has the same value at all levels while it should depend on the temperature. The thermal profile is independent of $c_p$ in radiative regions, as well as in condensation layers (the gas being pure H₂O), but not in layers where dry convection keeps the thermal gradient at the dry adiabatic lapse rate $\frac{\partial \ln T}{\partial \ln p} \bigg|_{\text{ad}} = R/\mu c_p$. The value of $c_p$ must therefore be chosen consistently with the modeled case. When a deep dry-convection zone is present, we thus iteratively adjust $c_p$ during the convergence process to the value of the heat capacity at the mean temperature of the convective layers.

Although these convective regions are often limited in size, we did evaluate the uncertainty caused by using a constant $c_p$ by bracketing our nominal profile with profiles obtained with the lowest and highest values of $c_p$ found in the nominal dry-convective region. Results are shown in Extended Data Figure 4 for the solar case where convective zones are rather extended, so that it can be
seen as a worst case scenario. The resulting changes on the equilibrium profiles
are limited with a variation on T<sub>surf</sub> that does not exceed 50 K unless large
internal heat fluxes (> 10 Wm<sup>-2</sup>) yield extended convective regions spanning
a broad range of temperatures, which increases this sensitivity.

The library also uses the ideal gas law to compute the density of water vapor from pressure and temperature. This assumption is rather valid for high temperatures, but can underestimate the density by a significant factor near the critical point, which affects the computed altitudes and, hence, the transit spectra. Although the library does not handle non-ideal equations of state in general, the altitudes used for the transmission spectra shown in figure 3 have been corrected for the non-ideal behavior of water vapor at high pressure. The largest effect was on the coldest Trappist-1b case that passes very close to the critical point and reduced the height of the atmosphere (transit radius) by less than 8% (0.3%).

## 2.8 Effect of the internal heat flow and surface pressure

When the geothermal flux is comparable to or larger than the stellar flux reach-
ing the surface, it contributes in heating the lower atmosphere. For a given
heat flow and instellation, its effect thus strongly depends on the spectral type.
Around the reddest dwarfs, the stellar flux at the surface can fall below the
value of Earth present geothermal flux. As a consequence, the influence of an
even low heat flux must be accounted for, while it can be neglected in the case
of Sun-like irradiation, as shown in figure 2. The surface pressure deter-
mines both the penetration of the stellar flux and the opacity of the deepest
atmospheric layers, which control, with the temperature, the thermal gradient
required to transport the surface heat flow upward. Extended Data Figure 6
shows how stellar spectrum, surface pressure and internal heat flux influence
together surface temperature and the nature of the thermal profile (convective
vs radiative). In some cases, increasing the geothermal flux can stabilize the
lower atmosphere against convection thanks to the increase of radiative dif-
fusivity with temperature. An early Venus with a steam atmosphere of more
than 100 bars would not, for instance, develop a fully convective structure
even for internal heat fluxes as high as 1000 Wm<sup>-2</sup> as its deepest layers would
remain radiative. Under a pure steam atmosphere, the solidification of Venus
surface would have occurred as soon as its internal heat flow decreases below
a value that depends on the amount of water: 1000, 10 and 2 W/m2, for 10,
100 and 270 bars respectively.

## 2.9 Composition of the atmosphere

While the vaporisation of the Earth ocean by the increase of solar luminosity
would produce, at least at first, a nearly pure steam atmosphere (with 0.3%
of N₂ and 0.0001 % of CO₂ [3]), the atmosphere released by the accretion
of a terrestrial planet or by a giant impact would be a mixture of outgassed
18 Cool runaway greenhouse

volatiles. The composition of the resulting atmosphere depends on the com-
position of the (proto)planet, accreted bodies or large impactor, and on the
distribution of the volatile compounds between the interior and the atmo-
sphere. Although most studies of such atmospheres assume a mixture of H₂O
and CO₂[53, 54, 18, 5, 19], the ratio between the two species is however poorly
constrained, as well as the presence of other constituents like N₂, CO, H₂ [55]
or sulfur-bearing compounds. The structure of atmospheres certainly varies a
lot within this variety of compositions and conclusions derived for pure steam
atmospheres cannot be generalized to post-accretion atmospheres. Exploring
the thermal structure of these atmospheres goes far beyond the goal of the
present study. However, the lesson learned from this work is that the convec-
tion hypothesis should at least be tested through the procedure described in
2.1. Whenever inverse climate modeling cannot be justified, Exo_k can be used
with appropriate opacities to compute a consistent structure.

2.10 Stellar flux required to melt the surface

When the melt fraction exceeds the critical value of 0.4, a silicate mixture with
a composition similar to that of the Earth mantle starts acquire the rheology
of a fluid. In the low-pressure limit this fraction is reached at 1620 K [29]. ??
gives the instellation that raises the surface to the *melting* temperature as a
function of the host-star type. Around Trappist-1, no planet is close enough

<figure><img src="image_10.png" /><figcaption>Extended Data Figure 6 Influence of the surface pressure and internal heat flux.</figcaption></figure>

a) Profiles for Venus gravity and insolation 4.5 Gyrs ago. b) Profiles for an Earth gravity with an ISR of 378 Wm⁻² and T-1 spectrum. At φint = 0 both cases have an OTR of 378 Wm⁻².
Dashed lines indicate convective layers. Dots indicate the surface of each individual profile.
to have its surface melted, including during the pre-main-sequence phase.

<table><thead><tr><th>Star</th><th>S<sub>eff</sub> at solidus</th><th>S<sub>eff</sub> for melting<sup>(*)</sup></th><th>S<sub>eff</sub> at liquidus</th></tr></thead><tbody><tr><td>Trappist 1</td><td>128</td><td><b>339</b></td><td>-</td></tr><tr><td>Proxima</td><td>24</td><td><b>100</b></td><td>519</td></tr><tr><td>M3</td><td>8.3</td><td><b>24</b></td><td>167</td></tr><tr><td>K5</td><td>3.5</td><td><b>11</b></td><td>52</td></tr><tr><td>Sun</td><td>2.1</td><td><b>7.4</b></td><td>32</td></tr><tr><td>+&phi;<sub>int</sub> = 0.09 Wm<sup>-2</sup></td><td>1.8</td><td><b>7.2</b></td><td>31</td></tr><tr><td>F1</td><td>1.6</td><td><b>6.2</b></td><td>28</td></tr></tbody></table>

(*) Melt fraction > 0.4, which implies a fluid rheology.

**ACKNOWLEDGEMENTS.** The authors sincerely thank Kevin Zahnle, Raymond T. Pierrehumbert, Robin Wordsworth and an anonymous referee for their particularly thorough and constructive review. JL and FS acknowledge funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 679030/WHIPLASH), and from the french state: CNES, Programme National de Planétologie (PNP), the ANR (ANR-20-CE49-0009: SOUND), and in the framework of the Investments for the Future programme IdEx, Université de Bordeaux/RRI ORIGINS.

FS and MT acknowledge support from BELSPO BRAIN (project B2/212/PI/PORTAL).

MT acknowledges support from the Tremplin 2022 program of the Faculty of Science and Engineering of Sorbonne University.

The authors thank the Generic PCM team for the teamwork development and improvement of the model.

MT acknowledges the use of the High-Performance Computing (HPC) resources of Centre Informatique National de l'Enseignement Superieur (CINES) under the allocation A0080110391 made by Grand Equipement National de Calcul Intensif (GENCI), which was essential to compute the 3-D GCM simulations presented in this work.

EB and GC acknowledge support from the Swiss National Science Foundation grant 200021_197176. Their work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606.

This research has made use of NASA's Astrophysics Data System.

**DATA AVAILABILITY.** Data generated by the atmospheric codes Exo_k and Generic PCM and used in this study are available at https://doi.org/10.5281/zenodo.6877001.

**CODE AVAILABILITY.** Exo_k is an open-source software. A complete documentation on how to install and use it can be found at http://perso.
20 REFERENCES

astrophy.u-bordeaux.fr/~jleconte/exo_k-doc/index.html.

The Generic PCM (Generic Global Climate Model, formerly known as LMDZwat generic) used in this work is the version 2528 that can be downloaded with documentation from the SVN repository at https://svn.lmd.jussieu.fr/Planeto / trunk / LMDZGENERIC/. More information and documentation are available at http://www-planets.lmd.jussieu.fr.

**AUTHOR CONTRIBUTIONS.** F.S. and J.L. laid the foundations for the study: F.S. wrote most of the paper, performed the 1D runs and made the figures – with guidance from J.L, M.T, E.B. and G.C. – while J.L. developed the code exo-k that made the study possible and wrote a significant part of the text, including exo-k description and manual. M.T. performed the 3D simulations, which were important to validate exo-k 1D profiles and pointed to radiative deep atmospheres. M.T., J.L. and G.C. worked on the spectroscopic data and formatted them for the study. All authors contributed to the response to the four reviewers and revisions of the article.

**COMPETING INTERESTS.** The authors declare that they have no competing financial interests.

**CORRESPONDENCE.** Correspondence and requests for materials should be addressed to F.S. (franck.selsis@u-bordeaux.fr).

**References**

[1] Sleep, N. H., Zahnle, K. J. & Lupu, R. E. Terrestrial aftermath of the Moon-forming impact. *Philosophical Transactions of the Royal Society of London Series A* **372**, 20130172–20130172 (2014).

[2] Ingersoll, A. P. The Runaway Greenhouse: A History of Water on Venus. *Journal of Atmospheric Sciences* **26**, 1191–1198 (1969).

[3] Kasting, J. F. Runaway and moist greenhouse atmospheres and the evolution of Earth and Venus. Icarus **74**, 472–494 (1988).

[4] Hamano, K., Abe, Y. & Genda, H. Emergence of two types of terrestrial planet on solidification of magma ocean. Nature **497**, 607–610 (2013).

[5] Lebrun, T. et al. Thermal evolution of an early magma ocean in inter-action with the atmosphere. *Journal of Geophysical Research (Planets)* **118**, 1155–1176 (2013).

[6] Kasting, J. F., Whitmire, D. P. & Reynolds, R. T. Habitable Zones around Main Sequence Stars. Icarus **101**, 108–128 (1993).

[7] Selsis, F. et al. Habitable planets around the star Gliese 581? A&A **476**, 1373–1387 (2007).

[8] Goldblatt, C., Robinson, T. D., Zahnle, K. J. & Crisp, D. Low simulated radiation limit for runaway greenhouse climates. *Nature Geoscience* **6**, 661–667 (2013).

[9] Kopparapu, R. K. et al. Habitable Zones around Main-sequence Stars: New Estimates. ApJ **765**, 131 (2013).
[10] Ding, F. & Pierrehumbert, R. T. Convection in Condensible-rich Atmospheres. ApJ 822, 24 (2016).

[11] Graham, R. J., Lichtenberg, T., Boukrouche, R. & Pierrehumbert, R. T. A Multispecies Pseudoadiabat for Simulating Condensable-rich Exoplanet Atmospheres. Planet. Sci. J. 2, 207 (2021).

[12] Koll, D. D. B. et al. Identifying Candidate Atmospheres on Rocky M Dwarf Planets via Eclipse Photometry. ApJ 886, 140 (2019).

[13] Turbet, M. et al. Day-night cloud asymmetry prevents early oceans on Venus but not on Earth. Nature 598, 276-280 (2021).

[14] Nakajima, S., Hayashi, Y.-Y. & Abe, Y. A study on the 'runaway greenhouse effect' with a one-dimensional radiative-convective equilibrium model. Journal of Atmospheric Sciences 49, 2256-2266 (1992).

[15] Chaverot, G., Turbet, M., Bolmont, E. & Leconte, J. How does the background atmosphere affect the onset of the runaway greenhouse? A&A 658, A40 (2022).

[16] Kopparapu, R. K. et al. Habitable Zones around Main-sequence Stars: Dependence on Planetary Mass. ApJ 787, L29 (2014).

[17] Watson, A. J., Donahue, T. M. & Kuhn, W. R. Temperatures in a runaway greenhouse on the evolving Venus: implications for water loss. Earth and Planetary Science Letters 68, 1-6 (1984).

[18] Abe, Y. & Matsui, T. Evolution of an Impact-Generated H₂O-CO₂ Atmosphere and Formation of a Hot Proto-Ocean on Earth. Journal of Atmospheric Sciences 45, 3081-3101 (1988).

[19] Hamano, K., Kawahara, H., Abe, Y., Onishi, M. & Hashimoto, G. L. Lifetime and Spectral Evolution of a Magma Ocean with a Steam Atmosphere: Its Detectability by Future Direct Imaging. ApJ 806, 216 (2015).

[20] Salvador, A. et al. The relative influence of H₂O and CO₂ on the primitive surface conditions and evolution of rocky planets. Journal of Geophysical Research (Planets) 122, 1458-1486 (2017).

[21] Lichtenberg, T. et al. Vertically Resolved Magma Ocean-Protoatmosphere Evolution: H₂, H₂O, CO₂, CH₄, CO, O₂, and N₂ as Primary Absorbers. Journal of Geophysical Research (Planets) 126, e06711 (2021).

[22] Leconte, J. Spectral binning of precomputed correlated-k coefficients. A&A 645, A20 (2021).

[23] Boukrouche, R., Lichtenberg, T. & Pierrehumbert, R. T. Beyond runaway: initiation of the post-runaway greenhouse state on rocky exoplanets. arXiv e-prints arXiv:2107.14150 (2021). URL https://ui.adsabs.harvard.edu/abs/2021arXiv210714150B.

[24] Thorngren, D., Gao, P. & Fortney, J. J. The Intrinsic Temperature and Radiative-Convective Boundary Depth in the Atmospheres of Hot Jupiters. ApJ 884, L6 (2019).

[25] Goldblatt, C. Habitability of Waterworlds: Runaway Greenhouses, Atmospheric Expansion, and Multiple Climate States of Pure Water Atmospheres. Astrobiology 15, 362-370 (2015).

[26] Nikolaou, A. et al. What Factors Affect the Duration and Outgassing of
22 REFERENCES

the Terrestrial Magma Ocean? ApJ **875**, 11 (2019).

[27] Turbet, M., Ehrenreich, D., Lovis, C., Bolmont, E. & Fauchez, T. The runaway greenhouse radius inflation effect. An observational diagnostic to probe water on Earth-sized planets and test the habitable zone concept. A&A **628**, A12 (2019).

[28] Andrault, D. et al. Deep and persistent melt layer in the Archaean mantle. *Nature Geoscience* **11**, 139–143 (2018).

[29] Salvador, A. & Samuel, H. Convective outgassing efficiency in planetary magma oceans: insights from computational fluid dynamics. *arXiv e-prints* arXiv:2209.06199 (2022).

[30] Leconte, J., Forget, F., Charnay, B., Wordsworth, R. & Pottier, A. Increased insolation threshold for runaway greenhouse processes on Earth-like planets. Nature **504**, 268–271 (2013).

[31] Turbet, M. et al. Revised mass-radius relationships for water-rich rocky planets more irradiated than the runaway greenhouse limit. A&A **638**, A41 (2020).

[32] Agol, E. et al. Refining the Transit-timing and Photometric Analysis of TRAPPIST-1: Masses, Radii, Densities, Dynamics, and Ephemerides. *Planet. Sci. J.* **2**, 1 (2021).

[33] Acuña, L. et al. Characterisation of the hydrospheres of TRAPPIST-1 planets. A&A **647**, A53 (2021).

[34] Dorn, C. & Lichtenberg, T. Hidden Water in Magma Ocean Exoplanets. ApJ **922**, L4 (2021).

[35] Turbet, M. et al. Modeling climate diversity, tidal dynamics and the fate of volatiles on TRAPPIST-1 planets. A&A **612**, A86 (2018).

[36] Makarov, V. V., Berghea, C. T. & Efroimsky, M. Spin-orbital Tidal Dynamics and Tidal Heating in the TRAPPIST-1 Multiplanet System. ApJ **857**, 142 (2018).

[37] Barr, A. C., Dobos, V. & Kiss, L. L. Interior structures and tidal heating in the TRAPPIST-1 planets. A&A **613**, A37 (2018).

[38] Papaloizou, J. C. B., Szuszkiewicz, E. & Terquem, C. The TRAPPIST-1 system: orbital evolution, tidal dissipation, formation and habitability. MNRAS **476**, 5032–5056 (2018).

[39] Bolmont, E. et al. Solid tidal friction in multi-layer planets: Application to Earth, Venus, a Super Earth and the TRAPPIST-1 planets. Potential approximation of a multi-layer planet as a homogeneous body. A&A **644**, A165 (2020).

[40] Fulton, B. J. et al. The California-Kepler Survey. III. A Gap in the Radius Distribution of Small Planets. AJ **154**, 109 (2017).

[41] Bean, J. L., Raymond, S. N. & Owen, J. E. The Nature and Origins of Sub-Neptune Size Planets. *Journal of Geophysical Research (Planets)* **126**, e06639 (2021).

[42] Luque, R. & Pallé, E. Density, not radius, separates rocky and water-rich small planets orbiting M dwarf stars. *Science* **377**, 1211–1214 (2022).
[43] Mousis, O. et al. Irradiated Ocean Planets Bridge Super-Earth and Sub-Neptune Populations. ApJ **896**, L22 (2020).

[44] Aguichine, A., Mousis, O., Deleuil, M. & Marcq, E. Mass-Radius Relationships for Irradiated Ocean Planets. ApJ **914**, 84 (2021).

[45] Toon, O. B., McKay, C. P., Ackerman, T. P. & Santhanam, K. Rapid calculation of radiative heating rates and photodissociation rates in inhomogeneous multiple scattering atmospheres. J. Geophys. Res. **94**, 16287–16301 (1989).

[46] Seeley, J. T. & Wordsworth, R. D. Episodic deluges in simulated hothouse climates. Nature **599**, 74–79 (2021).

[47] Gordon, I. E. et al. The HITRAN2016 molecular spectroscopic database. J. Quant. Spec. Radiat. Transf. **203**, 3–69 (2017).

[48] Rothman, L. S. et al. HITEMP, the high-temperature molecular spectroscopic database. *Journal of Quantitative Spectroscopy and Radiative Transfer* **111**, 2139–2150 (2010). URL https://ui.adsabs.harvard.edu/abs/2010JQSRT.111.2139R. ADS Bibcode: 2010JQSRT.111.2139R.

[49] Mlawer, E. J. et al. Development and recent evaluation of the MT_CKD model of continuum absorption. *Philosophical Transactions of the Royal Society of London Series A* **370**, 2520–2556 (2012).

[50] Rajpurohit, A. S. et al. Spectral energy distribution of M-subdwarfs: A study of their atmospheric properties. A&A **596**, A33 (2016).

[51] Wordsworth, R. D. et al. Gliese 581d is the First Discovered Terrestrial-mass Exoplanet in the Habitable Zone. ApJ **733**, L48 (2011).

[52] Lechevallier, L. et al. The water vapour self-continuum absorption in the infrared atmospheric windows: new laser measurements near 3.3 and 2.0 µm. *Atmospheric Measurement Techniques* **11**, 2159–2171 (2018).

[53] Matsui, T. & Abe, Y. Impact-induced atmospheres and oceans on Earth and Venus. Nature **322**, 526–528 (1986).

[54] Matsui, T. & Abe, Y. Evolution of an impact-induced atmosphere and magma ocean on the accreting Earth. Nature **319**, 303–305 (1986).

[55] Sossi, P. A., Tollan, P. M. E., Badro, J. & Bower, D. J. Solubility of water in peridotite liquids and the prevalence of steam atmospheres on rocky planets. *Earth and Planetary Science Letters* **601**, 117894 (2023).