2022

[1]Franck Selsis

[1]Laboratoire d’astrophysique de Bordeaux, Univ. Bordeaux, CNRS, Pessac,
France

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

3]Observatoire Astronomique de l’Université de Genève, Chemin Pegasi 51,
Sauverny 1290, Switzerland 4]Center for Life in the Universe, Faculty of
Science, University of Geneva, Geneva, Switzerland

# A cool runaway greenhouse without surface magma ocean

<EMAIL_ADDRESS>Jérémy Leconte<EMAIL_ADDRESS>Martin Turbet<EMAIL_ADDRESS>Guillaume Chaverot
<EMAIL_ADDRESS>Émeline Bolmont<EMAIL_ADDRESS>*  [  [
[

###### Abstract

Water vapour atmospheres with content equivalent to the Earth’s oceans,
resulting from impacts2014RSPTA.37230172S  or a high
insolation1969JAtS…26.1191I ; 1988Icar…74..472K , were found to yield a
surface magma ocean2013Natur.497..607H ; 2013JGRE..118.1155L . This was,
however, a consequence of assuming a fully convective
structure1969JAtS…26.1191I ; 1988Icar…74..472K ; 1993Icar..101..108K ;
2007A&A…476.1373S ; 2013Natur.497..607H ; 2013JGRE..118.1155L ;
2013NatGe…6..661G ; 2013ApJ…765..131K ; 2016ApJ…822…24D ; 2021PSJ…..2..207G .
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
works2013Natur.497..607H ; 2013JGRE..118.1155L . Around the reddest stars
($T_{eff}<3000$ K), surface magma oceans cannot form by stellar forcing alone,
whatever the water content. These findings 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 H2O 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 2019ApJ…886..140K
and 3D simulations have been achieved up to 10-30 bars2021Natur.598..276T ,
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 $T_{\textrm{surf}}$ 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 limit1992JAtS…49.2256N  of $\sim 280-290$ Wm-2
(depending on the H2O mixing ratio2022A&A…658A..40C ,
clouds2021Natur.598..276T  and gravity2014ApJ…787L..29K ) to a value of $\sim
274$ Wm-2 over a broad range of surface temperatures: from about 300 K to 1700
K. Above 1700 K the OTR increases again with T${}_{\textrm{surf}}$. 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-2
and a T${}_{\textrm{surf}}$ 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
$p_{\textrm{surf}}$ between 0.1 and 270 bars, stellar effective temperature
$T_{eff}$ between 2600 and 7000 K, OTR between 274 and $10^{5}$ Wm-2 and
internal heat flux between 0 and 1000 Wm-2. It is therefore necessary to self-
consistently calculate the atmospheric thermal profile, without imposing
convection, following the tracks of pioneering works on runaway atmospheres
1984E&PSL..68….1W  and formation/condensation of post-accretion atmospheres
1988JAtS…45.3081A . This is required to reassess results obtained assuming
fully convective profiles on the divergent fate of Venus and the
Earth1988Icar…74..472K ; 2013Natur.497..607H , the circumstellar region where
surface liquid water can be stable1993Icar..101..108K ; 2013ApJ…765..131K  and
the occurrence of surface magma oceans2013Natur.497..607H ;
2013JGRE..118.1155L ; 2015ApJ…806..216H ; 2017JGRE..122.1458S ;
2021JGRE..12606711L .

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 2021A&A…645A..20L  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_{\mathrm{eff}}=2600$ K) to a
F1-type star ($T_{\mathrm{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_{\textrm{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 boukrouche_beyond_2021 . 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 d$B$/d$T$ $\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 2019ApJ…884L…6T . 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_{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.

Figure 1:  Steam atmosphere thermal profiles. Each group of five profiles
shares the same bolometric emission (274 and 10,000 Wm-2, 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${}_{\textrm{eff}}=1$ corresponds
to Earth’s present insolation: 341.5 Wm-2). Dashed lines indicate dry
convection.

Compared with computed equilibrium profiles, convective structures
dramatically overestimate $T_{\textrm{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).

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_{\textrm{int}}$. Dashed curves are obtained
with Earth geothermal flux. For T-1 and $\phi_{\textrm{int}}<0.01$ 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_{\textrm{surf}}$ is below the critical temperature
of H2O and $p_{\textrm{surf}}$ is given by the saturation vapor pressure of
H2O. These curves are based on the profiles shown in Extended Data Figure 2

When using convective profiles, runaway extends to surface-melting
temperatures 2013ApJ…765..131K ; 2014ApJ…787L..29K ; 2013NatGe…6..661G ;
2015AsBio..15..362G ; 2015ApJ…806..216H ; 2019ApJ…875…11N ; 2019A&A…628A..12T
; boukrouche_beyond_2021 . However, in all converged cases and as seen in
Figure 2, planetary cooling becomes sensitive to $T_{\textrm{surf}}$ at lower
temperature: 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
2018NatGe..11..139A . At 1620 K the melt fraction reaches 40% and dominates
the rheology Salvador_Samuel22 . This melting temperature requires
S${}_{\textrm{eff}}>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 internal heat flows below 2, 10 and 100 Wm-2 for
p${}_{\textrm{surf}}$ = 270, 100 and 10 bars of steam respectively. The post-
accretion atmosphere of Venus must have consisted in other species than H2O so
we 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 findings2013Natur.497..607H ;
2013JGRE..118.1155L .

We find a stable steam atmosphere for a solar flux as low as
S${}_{\textrm{eff}}=0.98$, a consequence of a known bistability
2013NatGe…6..661G ; 2021Natur.598..276T : a higher ISR is required to vaporize
an ocean (S${}_{\textrm{eff}}=1.11$ 2013Natur.504..268L ) than to condense a
steam atmosphere (S${}_{\textrm{eff}}=0.92$ 2021Natur.598..276T ). 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_{\textrm{surf}}$ by convective
profiles is considerable: 1200 K for S${}_{\textrm{eff}}=0.8$, just above the
runaway threshold, and 2700 K for S${}_{\textrm{eff}}=58$. No 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 2013Natur.497..607H .

This work shows that ultracool dwarfs ($T_{eff}\leq 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 H2O (see
Figure 2 and Extended Data Figure 6).

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_{\textrm{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_{\textrm{int}}=0$ or 5 Wm-2 (a range of estimated tidal dissipations).

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 $\sim 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 $\mu$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 2019A&A…628A..12T  by the steam atmosphere yielded an underestimation
of the maximum water content of T-1b, c and d 2020A&A…638A..41T ;
2021PSJ…..2….1A ; 2021A&A…647A..53A ; 2021ApJ…922L…4D . With inverse climate
modeling, the internal heat flux only affects the value of $T_{\textrm{surf}}$
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-22018A&A…612A..86T ; 2018ApJ…857..142M ;
2018A&A…613A..37B ; 2018MNRAS.476.5032P ; 2020A&A…644A.165B  adds $\sim 100$
km (1.4%) to the radius. The current 1-$\sigma$ relative precision on the
radius of T-1b is 1%2021PSJ…..2….1A .

Our calculations thus have important consequences for the interpretation of
the exoplanet radius gap 2017AJ….154..109F  and of the true nature of sub-
Neptune size planets 2021JGRE..12606639B ; 2022Sci…377.1211L . Indeed, we show
that water-rich atmospheres are significantly less extended than previously
calculated2020A&A…638A..41T ; 2020ApJ…896L..22M ; 2021ApJ…914…84A , 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 H2O 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 H2O 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 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
runaway2013Natur.504..268L , a simple test demonstrates that convection cannot
be maintained by stellar radiation beyond the runaway threshold. At
equilibrium, net ($\uparrow-\downarrow$) energy fluxes must satisfy at each
level:

$F_{*}(p)+F_{\textrm{E}}(p)+F_{\textrm{C}}(p)=\phi_{\textrm{Int}}$

where $p$ is the pressure level, $F_{*}$ is the net stellar flux,
$F_{\textrm{E}}$ is the net thermal emission flux and $\phi_{\textrm{int}}$ is
the internal heat flux. As convection can only transport heat upward, this
convention implies that the convective flux, $F_{\textit{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_{\textrm{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-2 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_{\textrm{surf}}>700K$ for
$P_{\textrm{surf}}$=270 bars, a solar spectrum, and Earth gravity.

Extended Data Figure 1: Converged versus fully convective profiles. The steam
atmosphere is modeled here for a solar flux of 445.7 Wm-2
($S_{\textrm{eff}}=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 equilibrum. Panel c) Evolution from the initial dashed
“convective” profile towards the converged state, computed with the Exo_k
suite2021A&A…645A..20L  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.

###  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_{\textrm{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 library2021A&A…645A..20L  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 approximationTMA89 . 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
atmospheres2022A&A…658A..40C .

#### 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 atmospheres2013Natur.504..268L . 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 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
climates2021Natur.599…74S .

#### Time integration

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

$\displaystyle c_{p}\frac{\partial\hat{T}}{\partial t}=\hat{H}(\hat{T}),$ (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}_{\mathbb{K}}$, the radiative heating rates can be
linearized through

$\displaystyle\hat{H}^{\mathrm{rad}}(\hat{T}_{\mathbb{K}}+\delta\hat{T})=\hat{H}^{\mathrm{rad}}(\hat{T}_{\mathbb{K}})+\mathbb{K}\cdot\delta\hat{T},$
(2)

where the kernel $\mathbb{K}$ is the Jacobian matrix of the heating rates (of
dimension $N_{\mathrm{lay}}\times N_{\mathrm{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}^{\mathrm{rad}}(\hat{T}_{\mathbb{K}})=\hat{0}$.
Combining the above equations yields the following evolution for the
temperature perturbation:

$\displaystyle c_{p}\frac{\partial}{\partial
t}\delta\hat{T}=\mathbb{K}\cdot\delta\hat{T}.$ (3)

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

$\displaystyle\delta\hat{T}=\delta\hat{T}_{0}\cdot
e^{-\hat{\mathrm{Diag}}(\mathbb{K})t/c_{p}},$ (4)

where $\hat{\mathrm{Diag}}(\mathbb{K})$ is the vector formed by all the
diagonal terms in $\mathbb{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{\tau}=c_{p}/\hat{\mathrm{Diag}}(\mathbb{K})$ that differs from layer
to layer. In the baseline evolution, the smallest radiative timescale in our
atmosphere, $\tau_{\mathrm{rad}}=\min(\hat{\tau})$, 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_{\mathrm{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 $\mathbb{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.

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
evolution, 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 $\Delta 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 reevaporation of rains.

Acceleration in radiative zones. We start by defining a base timescale,
$\tau_{\mathrm{b}}$, as the smallest radiative timescale of a radiative layer
in the atmosphere. Then, for any radiative layer $\mathrm{n}$, the heating
rate is multiplied by $\tau_{\mathrm{n}}/\tau_{\mathrm{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
redistributed 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.

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.

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_{\mathrm{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_{\mathrm{av}}/\tau_{\mathrm{rad}}$.

###  2.3 Simulations

Unless indicated otherwise, simulations presented in this study were run for a
pure H2O atmosphere, a surface gravity of 1g and a surface pressure of 270
bars corresponding to vaporized Earth oceans. Our radiative transfer uses the
H2O linelists from HITRAN20162017JQSRT.203….3G  and
HITEMP2010rothman_hitemp_2010 . 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 H2O-H2O continuum parametrization (version
4.0.12012RSPTA.370.2520M ) 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 linelists2013NatGe…6..661G : we use HITRAN at wavelengths
above 3$\mu$m and HITEMP below 3$\mu$m. Opacities are treated with a
correlated-k approach in 69 spectral bands from 0.28 to 1000
$\mu$m2021Natur.598..276T . 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
spectra2016A&A…596A..33R  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
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
GCM2011ApJ…733L..48W ; 2021Natur.598..276T . In addition to an existing
simulation done for the Sun2021Natur.598..276T , 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 states2021Natur.598..276T . Both 1D and 3D simulations were done for
$p_{\textrm{surf}}=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.

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 H2O and 5% of N2
and the instellation is 500 Wm-2 ($S_{\textrm{eff}}=1.42$) in all cases. The
opacities used in both models are the same, differences come from circulation
and cloud radiative effects.

Exo_k and Generic PCM are in close agreement for the stellar heating rates.
This match is allowed by the absence of clouds on the
dayside2021Natur.598..276T . Thermal profiles show an excellent match at
$p<10^{5}$ 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 $\sim 20$ K at the surface. This small departure
is attributed to the greenhouse warming of nightside clouds2021Natur.598..276T
. 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
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.

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.

###  2.6 Sensitivity to opacity

Radiative-convective profiles computed in this study depend on H2O 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 H2O
self continuum. Compared to version 3.5, the temperature dependency in the
version 4.0.1 was modified to account for new measurements and
parametrizations2018AMT….11.2159L . 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_{\textrm{surf}}$ do not exceed 50 K. Second, it however reminds us that
new data will inevitably alter the quantitative results.

Extended Data Figure 5: Sensitivity to a significant change in the H2O-H2O
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.

###  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 H2O), but not in layers where dry
convection keeps the thermal gradient at the dry adiabatic lapse rate
$\left.\frac{\partial\\!\\!\ \ln T}{\partial\\!\\!\ \ln
p}\right|_{\mathrm{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_{\textrm{surf}}$ that
does not exceed 50 K unless large internal heat fluxes ($>10$ Wm-2) 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 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
reaching 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 Figure 2.
The surface pressure determines 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 diffusivity 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-2 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.

Extended Data Figure 6: Influence of the surface pressure and internal heat
flux. a) Profiles for Venus gravity and insolation 4.5 Gyrs ago. b) Profiles
for an Earth gravity with an ISR of 378 Wm-2 and T-1 spectrum. At
$\phi_{\textrm{int}}=0$ both cases have an OTR of 378 Wm-2. Dashed lines
indicate convective layers. Dots indicate the surface of each individual
profile.

###  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
N2 and 0.0001 % of CO2 1988Icar…74..472K ), the atmosphere released by the
accretion of a terrestrial planet or by a giant impact would be a mixture of
outgassed volatiles. The composition of the resulting atmosphere depends on
the composition of the (proto)planet, accreted bodies or large impactor, and
on the distribution of the volatile compounds between the interior and the
atmosphere. Although most studies of such atmospheres assume a mixture of H2O
and CO21986Natur.322..526M ; 1986Natur.319..303M ; 1988JAtS…45.3081A ;
2013JGRE..118.1155L ; 2015ApJ…806..216H , the ratio between the two species is
however poorly constrained, as well as the presence of other constituents like
N2, CO, H2 2023E&PSL.60117894S  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 convection 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 Salvador_Samuel22 . LABEL:tab:ED_Table_ISR_melt 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 to have its surface
melted, including during the pre-main-sequence phase.

Star |  $S_{\textrm{eff}}$ at solidus |  $S_{\textrm{eff}}$ for melting(∗) |  $S_{\textrm{eff}}$ at liquidus
---|---|---|---
Trappist 1 | 128 | 339 | -
Proxima | 24 | 100 | 519
M3 | 8.3 | 24 | 167
K5 | 3.5 | 11 | 52
Sun | 2.1 | 7.4 | 32
+$\phi_{\textrm{int}}=0.09$ Wm-2 | 1.8 | 7.2 | 31
F1 | 1.6 | 6.2 | 28
(∗) 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.astrophy.u-bordeaux.fr/~jleconte/exo_k-doc/index.html. The
Generic PCM (Generic Global Climate Model, formerly known as LMDZ.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/LMDZ.GENERIC/. 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 interaction 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 H2O-CO2 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 H2O and CO2 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: H2, H2O, CO2, CH4, CO, O2, and N2 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 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).