# The Runaway Greenhouse Effect on Hycean Worlds

Hamish Innes Atmospheric, Oceanic and Planetary Physics, Department of
Physics, University of Oxford, UK  Shang-Min Tsai Atmospheric, Oceanic and
Planetary Physics, Department of Physics, University of Oxford, UK  Department
of Earth Sciences, University of California, Riverside, California, USA
Raymond T. Pierrehumbert Atmospheric, Oceanic and Planetary Physics,
Department of Physics, University of Oxford, UK Hamish Innes
<EMAIL_ADDRESS>
###### Abstract

Hycean worlds are a proposed subset of sub-Neptune exoplanets with substantial
water inventories, liquid surface oceans and extended hydrogen-dominated
atmospheres that could be favourable for habitability. In this work, we aim to
quantitatively define the inner edge of the Hycean habitable zone using a 1D
radiative-convective model. As a limiting case, we model a dry hydrogen-helium
envelope above a surface ocean. We find that 10 to 20 bars of atmosphere
produces enough greenhouse effect to drive a liquid surface ocean
supercritical when forced with current Earth-like instellation. Introducing
water vapour into the atmosphere, we show the runaway greenhouse instellation
limit is greatly reduced due to the presence of superadiabatic layers where
convection is inhibited. This moves the inner edge of the habitable zone from
$\approx 1$ AU for a G-star to 1.6 AU (3.85 AU) for a Hycean world with a
H2-He inventory of 1 bar (10 bar). For an M-star, the inner edge is
equivalently moved from 0.17 AU to 0.28 AU (0.54 AU). Our results suggest that
most of the current Hycean world observational targets are not likely to
sustain a liquid water ocean. We present an analytical framework for
interpreting our results, finding that the maximum possible OLR scales
approximately inversely with the dry mass inventory of the atmosphere. We
discuss the possible limitations of our 1D modelling and recommend the use of
3D convection-resolving models to explore the robustness of superadiabatic
layers.

††journal: ApJ

##  1 Introduction

###  1.1 Sub-Neptunes and Hycean Worlds

Theorists have long predicted the existence of small, water-rich planets
forming outside the ice line with no solar system analogue (Kuchner, 2003;
Léger et al., 2004). With the launch of the Kepler space telescope, the
population of detected sub-Neptunes increased by orders of magnitude (Batalha,
2014). They are considered one of the most abundant types of planet in the
galaxy (alongside super-Earths) (Petigura et al., 2013; Marcy et al., 2014;
Winn & Fabrycky, 2015). Limited data leaves their inner compositions largely
unconstrained, with structures ranging from almost pure water to rocky/iron
cores producing the same mass and radius when given a hydrogen-dominated
envelope (Rogers & Seager, 2010; Rogers et al., 2011; Marcy et al., 2014;
Lopez et al., 2012; Rogers, 2015; Madhusudhan et al., 2020). Population
demographics provide insights into the composition of smaller planets. The
discovery of the radius valley (Fulton et al., 2017) separating the population
of smaller super-Earths and larger sub-Neptunes provides a way to distinguish
these two types of planet. By studying how the position of the valley varies
with instellation, there is consensus that it is shaped by atmospheric loss,
either via core-powered mass loss (Gupta & Schlichting, 2019) or
photoevaporative loss (Owen & Wu, 2013).

“Hycean Worlds” (Madhusudhan et al., 2021) are a proposed subset of water-rich
sub-Neptunes with hydrogen-dominated atmospheres. The sub-Neptune K2-18 b
(Cloutier et al., 2017, 2019) has been the subject of multiple studies owing
to the claimed detection of water vapour in its atmosphere (Benneke et al.,
2019; Tsiaras et al., 2019), with Madhusudhan et al. (2020) showing that its
mass and radius were consistent with being a Hycean world. Further work on the
habitability (Madhusudhan et al., 2021) and interior structure (Nixon &
Madhusudhan, 2021) of Hycean worlds confirm there is a wide range of mass-
radius parameter space potentially occupied by this type of planet. Their
larger radii and masses compared to terrestrial-type habitable planets makes
them more amenable to observational characterisation, with several temperate
candidates having already been identified (see Table 1 of Madhusudhan et al.,
2021).

###  1.2 Planetary Habitability and the Runaway Greenhouse Effect

####  1.2.1 Background

Most models of planetary habitability focus on atmospheres where the mean
molecular weight (MMW) of the background atmosphere is greater than the
condensible component (and for good reason – such is the case on Earth).
Pioneering work related to Venus’s atmosphere (Simpson, 1929; Komabayasi,
1967; Ingersoll, 1969; Nakajima et al., 1992) established a maximum
instellation above which a planet can no longer increase its infrared cooling
into space. This is known as the “runaway greenhouse” limit (Ingersoll, 1969).
As high surface temperatures make water vapour a non-negligible component of
the atmosphere, there comes a point where the atmosphere becomes optically
thick to infrared radiation, decoupling the surface temperature from a fixed
photospheric temperature. If a planet’s instellation exceeds this value, its
water inventory evaporates into the atmosphere, with photodissociation and
atmospheric escape eventually driving it into space. The inner edge of the
habitable zone (HZ) is typically calculated assuming the runaway greenhouse
effect is the limiting factor influencing habitability (Kasting et al., 1993;
Kopparapu, 2013).

####  1.2.2 The Impact of an H2-He Background and Convective Inhibition

As interest around hydrogen-dominated habitable worlds has grown, it is
important to assess the impact of hydrogen as a background gas as opposed to
nitrogen or other high MMW gas. Hydrogen gas’s main source of opacity comes
from its collision-induced absorption (CIA) spectrum. Since the strength of
absorption scales as the density of the atmosphere squared, CIA becomes an
important factor in H2 atmospheres of around 1 bar or thicker, with 40 bars of
pure H2 allowing habitable surface temperatures out to 1.5 AU for an M star
and 10 AU for a G star (Pierrehumbert & Gaidos, 2011). Accounting for
absorption and water vapour feedbacks is necessary to further constrain the
inner edge of the habitable zone. In Madhusudhan et al. (2021), inner edge
calculations are, for the most part, made assuming a constant water mixing
ratio of 10% (or less, if condensation occurs). Calculations with a water-
saturated atmosphere have also been performed (e.g., Figure 9 of Piette &
Madhusudhan, 2020), however their model neglects compositional effects (see
below). When hydrogen is used as the background gas in classical runaway
greenhouse calculations (Koll & Cronin, 2019), the presence of H2 introduces
novel behaviour on account of its low MMW. Unlike for high MMW background
gases which tend to raise the runaway greenhouse limit and lead to non-
monotonic changes in outgoing longwave radiation (OLR) with surface
temperature, H2 atmospheres monotonically approach the pure steam limit with
increasing surface temperature.

Our work aims to build on Koll & Cronin (2019) by introducing the effect of
MMW-induced convective inhibition on the temperature profiles of hydrogen-
dominated planets. For atmospheres with a condensing component heavier than
the background gas, decreasing temperatures with altitude leads to a sharp
decrease in MMW between the lower atmosphere and the upper atmosphere. If the
concentration of the condensible is high enough, compositional gradients
stabilise the atmosphere to convection (Guillot, 1995; Li & Ingersoll, 2015;
Leconte et al., 2017) and double-diffusive instabilities (Leconte et al.,
2017) even if the lapse rate is super-adiabatic. Molecular gradient-induced
convective inhibition has been invoked to explain periodic storms in Saturn’s
atmosphere (Li & Ingersoll, 2015), stable wave ducts for gravity waves in
Jupiter’s atmosphere (Ingersoll & Kanamori, 1995) and explain the step-like
behaviour of methane’s abundance in Uranus’s condensation layer (Irwin et al.,
2022). Neglecting this effect results in an underestimation of the deep
atmospheric temperature in gas giant planets (Leconte et al., 2017). Recent
studies have investigated convective inhibition due to condensing silicates in
a H2-He envelope, which alters the radius and cooling times of young super-
Earths and sub-Neptunes (Misener & Schlichting, 2022; Markham et al., 2022)

The atmospheric structure of runaway atmospheres is often assumed to be on a
moist adiabat integrated upwards from the surface. However, if the surface
temperature is high enough, high moisture contents will inhibit convection and
lead to radiative layers in the lower atmosphere. We will discuss this effect
in further detail in Section 3. Our aim is to calculate new limiting
instellations for sub-Neptunes hosting a liquid water ocean with a hydrogen-
dominated atmosphere.

####  1.2.3 Super-Runaway States

Even a liquid water ocean can be too hot to be habitable for life as we know
it, but from the standpoint of planetary structure a second important
transition occurs when the surface temperature of the water layer reaches the
critical point. At this point, the liquid-gas phase transition of water
disappears. More importantly for our purposes, H2 is completely miscible with
supercritical water (Soubiran & Militzer, 2015), so that gaseous solubility is
no longer limited by Henry’s Law (miscibility is expected for helium and other
gases as well.) In consequence, a hydrogen envelope could not remain distinct
from the supercritical water interior, but instead would mix into it and be
diluted into the supercritical ocean. It is an interesting theoretical
question how long it would take for such mixing to occur if the planet started
out with a distinct hydrogen envelope, but the more plausible scenario for a
planet with an initial H2-H2O composition is that the system would start out
with a supercritical H2-H2O mixture. The H2 would never phase-separate into a
distinct layer if the equilibrium radiation balance maintains supercritical
conditions with respect to water. For this reason, we put particular emphasis
on defining planetary parameters for which the ocean surface temperature
approaches the critical point of water. The chief goal of this paper is to
identify conditions in which a subcritical liquid water ocean can coexist with
a hydrogen-rich atmosphere, though we do offer some speculations on the state
a sub-Neptune settles into when those conditions are not met.

###  1.3 Paper Structure

The paper is structured as followed. In Section 2 we will discuss the H2-He
inventory required to produce a strong enough greenhouse effect to drive
surface temperatures to the critical point of water. This will give an upper-
bound on the instellation a planet can receive from its host star before
entering a runaway state. In Section 3 we extend our model to deal with
atmospheres with a water vapour component. We will then discuss our results in
the context of the literature in Section 4.

##  2 Hydrogen-helium atmosphere above a water ocean

###  2.1 Model

To find an upper bound on the hydrogen inventory of a planet that maintains a
liquid water surface, we first model a hydrogen-helium atmosphere above a
water ocean. We assume the temperature structure of the atmosphere is the dry
adiabat for a hydrogen-helium mixture with solar abundances (Asplund et al.,
2009), neglecting the contribution of water vapour to the atmosphere (similar
to the calculation in Pierrehumbert & Gaidos (2011)). This is an unphysical
assumption given warm temperatures will naturally lead to the water
evaporating from the surface ocean and mixing with the hydrogen-helium gas.
However, it serves as an easily calculated upper limit on the hydrogen-helium
content of a planet with temperate instellation since we expect the addition
of water vapour to act as a greenhouse gas, reducing the ability of the planet
to cool. Our choice of a purely adiabatic atmosphere is justified from test
runs using a full radiative-convective iteration (see Section 3). We found
that the tropopause (at $\approx$ 0.05 bars) was typically at lower pressures
than the infra-red photosphere (at $>0.1$ bar).

Radiative fluxes were calculated in the longwave (LW) and shortwave (SW)
regions of the spectrum using the SOCRATES radiative transfer code (Edwards &
Slingo, 1996), the details of which can be found in Appendix A. In the SW
calculation, a surface albedo of 0.12 is specified (Goldblatt et al., 2013).

Since we are interested in calculating the maximum hydrogen inventory possible
before an ocean is driven supercritical, we calculate fluxes for an atmosphere
on a dry adiabat with surface temperature equal to the critical temperature of
water, $647\text{\,}\mathrm{K}$. We vary the surface pressure logarithmically
between 0.5 bar and the critical pressure of water, 220 bar. The SW fluxes
were calculated using a solar instellation value of
$1361\text{\,}\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$, however we note that
since the temperature structure of the atmosphere is fixed, we can simply
multiply our results by a constant factor to retrieve the SW fluxes for any
arbitrary instellation. The surface gravity was kept constant at the Earth
value, $9.81\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-2}$, throughout. We
consider cases where the incoming stellar spectrum has that of a G star or an
M star (details of which are found in Appendix A).

###  2.2 Results

Only moderate pressures of hydrogen and helium are required to force surface
temperatures to supercritical values. Figure 1 shows the OLR and SW absorption
of these atmospheres as a function of surface pressure. The intersection
points between the OLR curve (green) and the SW absorption curves (blue and
orange depending on stellar type) signify atmospheric configurations in global
equilibrium. For a pure H2-He atmosphere with solar instellation, $S_{0}$,
approximately 10 bars of atmosphere will cause a large enough greenhouse
effect due to CIA to drive the surface ocean supercritical. The SW absorption
in the G-star experiments was reduced relative to the M-star thanks to the
enhanced Rayleigh scattering cross sections at low wavelengths. However, the
change in stellar type causes minimal differences in the qualitative behaviour
of our model. Naturally, lower instellation values require larger atmospheric
masses to warm the surface to the critical temperature, and atmospheres
irradiated with less than 1% of solar radiation would require a surface
pressure greater than the critical pressure of water to reach the critical
temperature.

Figure 1: OLR and SW absorption as a function of surface pressure. The green
curve represents the OLR of an H2-He atmosphere on a dry adiabat with surface
temperature equal to the critical temperature of water. The blue and orange
curves show the SW absorption at different instellations relative to solar
($S_{0}$) for an M star and G star respectively.

We note the surface gravity in these experiments was held constant at Earth’s
value. As noted in Pierrehumbert & Gaidos (2011), the OLR depends on the
surface gravity in the combination of $p_{s}^{2}/g$, so changing $g$ can be
accounted for simply by rescaling $p_{s}$. From Figure 1 we see that
increasing the surface gravity from its terrestrial value would be equivalent
to lowering the surface pressure, i.e. increasing the OLR (assuming the change
in Rayleigh scattering, which scales as $p_{s}/g$, does not cause the
planetary albedo to change significantly). Therefore, for the same level of
instellation, a planet with a higher surface gravity can host a higher
pressure atmosphere before going supercritical. Nevertheless, we expect the
surface gravity of temperate sub-Neptunes to be of a similar order of
magnitude to the Earth. The four temperate sub-Neptunes with constrained mass
and radius (Table 2 of Pierrehumbert, 2023) have surface gravities spanning
the range
$4\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-2}24\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-2}$.

##  3 Adding water vapour to the atmosphere

###  3.1 Model Setup

####  3.1.1 Thermodynamics

Water is added to each layer at its saturation pressure value (i.e. the
relative humidity is assumed to be 100%). We use the analytic expressions for
the saturation vapour pressure given by Wagner & Pruß (2002), generating a
lookup table of values for each temperature for faster on-the-fly calculation.

####  3.1.2 Radiative Procedures

As in Section 2, radiative fluxes are calculated with the SOCRATES radiative
transfer code. Unlike in Section 2 where an inverse modelling approach was
taken, we now timestep the model by iterating the temperature according to:

$\partialderivative{T}{t}=A(p)\partialderivative{F_{\text{net}}}{p}$ (1)

where $T$ is the temperature, $F_{\text{net}}$ is the net flux two-stream flux
calculated at the edge of each model layer (defined such that $F_{\text{net}}$
is positive if the net flux is upwards, in the $-p$ direction). The
coefficient $A(p)$ varies in each model layer in order to approach radiative
equilibrium (where $\partial_{p}F_{\text{net}}=0$) as fast as possible. We use
the method described in Malik et al. (2017), where $A(p)$ is decreased if
temperature oscillations are detected in that temperature layer and increased
on each iteration otherwise. For our purposes we only need to compute the
correct equilibrium profile, and do not need to compute the actual time series
as the system adjusts to equilibrium.

Using this approach, the independent variable is the instellation, $S$, as
opposed to the surface temperature, $T_{s}$, which is permitted to vary until
the atmosphere reaches local and global radiative equilibrium.

####  3.1.3 Convective Adjustment

After each iteration, the temperature-pressure profile is checked for
instability to convection. We use the instability criterion modified for use
in atmospheres where water vapour is at saturation with a lower mean molecular
weight background gas (Leconte et al., 2017):

$\displaystyle(\nabla_{\text{am}}-\nabla_{\text{ad}})(1-\beta q\varpi)>0$ (2a)
$\displaystyle\text{with}\quad\nabla_{x}\equiv\derivative{\ln T_{x}}{\ln p},$
(2b) $\displaystyle\beta\equiv\derivative{\ln p_{\text{sat}}}{\ln T},$ (2c)
$\displaystyle\varpi\equiv(1-\mu_{d}/\mu_{v}).$ (2d)

The “am” and “ad” suffixes denote “ambient” and “adiabatic” respectively,
$p_{\text{sat}}$ is the saturation vapour pressure and $\mu_{d}$ and $\mu_{v}$
are the MMWs of the dry and condensing vapour gases respectively. The first
term in Equation 2a is the Schwarzschild criterion (Schwarzschild, 1906) for
convective instability when the ambient lapse rate is greater than the moist
adiabatic lapse rate. The second term represents the effect of mean molecular
weight gradients in an environment where water vapour is held at its
saturation vapour pressure value. When the water vapour concentration, $q$,
exceeds a critical value, $q_{c}\equiv 1/\beta\varpi$, the parcel is no longer
unstable to convection despite any super-adiabatic lapse rates. The value of
the moist adiabatic lapse rate, $\nabla_{\text{ad}}$ is calculated using the
expression in Ding & Pierrehumbert (2016), which accounts for water vapour
being a non-negligible component of the atmosphere but neglects the effects of
retained condensates. We replace occurences of $L/R_{v}T$ in the lapse rate
formula with the more appropriate $\beta$ factor obtained directly from a
lookup table, and any other factors of $L$ are calculated using a lookup table
using fits from Wagner & Pruß (2002). This accounts for vanishing $L$ as $T$
approaches the critical temperature of water ($647\text{\,}\mathrm{K}$).
However, we note that at temperatures approaching the critical point of water,
the non-idealness of water vapour will become important but is neglected in
this current work for simplicity. Convective adjustment is performed pairwise
on layers working upwards from the bottom of the atmosphere. If Equation 2a is
satisfied, we adjust the temperature of both layers to the moist adiabat. This
procedure is repeated until convergence is reached.

Note that when $q>q_{c}$, a saturated layer with lapse rate
$\nabla_{\text{am}}$ less than the adiabatic lapse rate – even an isothermal
layer – is unstable, and in our scheme is convectively adjusted to the
saturated adiabat; radiative cooling may then cause further steepening of the
lapse rate.

The implementation of Equation 2a leads to two radiative zones forming in the
atmosphere when temperatures near the surface are high enough – a traditional
stratosphere in the upper atmosphere with low lapse rates and a moisture-
inhibited radiative zone below the first convective region. Since the
atmosphere is optically thick in this region, the radiative lapse rate is
usually high and requires higher vertical resolution. We ran the model with
200 layers, with the bottom 100 layers dedicated to resolving this radiative
region of the atmosphere. In the traditional stratosphere we implement a cold
trap for moisture, setting the moisture concentration to its minimum value at
pressures lower than the moisture minimum.

####  3.1.4 Surface Physics

At the surface we implement a rudimentary heat transfer scheme which transfers
sensible heat from the surface to the lowest model layer. This prevents the
surface temperature becoming unphysically hotter than the lowest model layer
and is a $0^{\text{th}}$ order attempt to represent some of the heat exchange
processes in the turbulent boundary layer. The surface energy equation is:

$\rho_{s}c_{s}h\derivative{T_{s}}{t}=-F_{\text{net}}(p_{s})+c_{p}\rho_{\text{air}}C_{D}U(T(p_{s})-T_{s})$
(3)

where values for the surface density ($\rho_{s}$) and heat capacity ($c_{s}$)
are taken to be those of liquid water and $h$ (the bucket depth) is taken as
$1\text{\,}\mathrm{m}$. The drag coefficient $C_{D}$ and the characteristic
drag velocity $U$ are taken to be 0.001 and
$10\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ respectively, from
Pierrehumbert (2010). The effect of varying these parameters had a negligible
effect on the overall temperature structure of the atmosphere and the final
conclusions of this work.

In traditional calculations of the runaway greenhouse limit using an inverse
modelling approach (e.g. Kasting et al., 1993), the surface pressure $p_{s}$
is taken to be the sum of a constant dry component, $p_{0}$, and the
saturation vapour pressure at the specified surface temperature $T_{s}$.
Surface pressure can increase with temperature, allowing an increase in total
atmospheric mass as more water vapor is added to the atmosphere. However,
holding $p_{0}$ fixed as temperature increases leads to an implied change in
the dry mass of the atmosphere as temperature changes, because the mean
molecular weight of the atmosphere varies and the concentration of the dry
mass is not uniform over the profile. For inverse climate modeling, it is
straightforward to allow $p_{s}$ to vary with temperature, but with radiative-
convective calculations that compute equilibria using time-stepping or related
iterations, the re-gridding needed to allow $p_{s}$ to vary becomes unwieldy
and can lead to numerical issues. For that reason, we make the additional
simplification of holding $p_{s}$ itself fixed, until conditions make such an
assumption physically inconsistent. The expedient of holding $p_{s}$ fixed was
also used in the calculations of Figure 9 in Piette & Madhusudhan (2020),
which were carried out with a water-saturated atmosphere. In such
calculations, there is an additional reduction in implied dry mass as
temperature increases, as $p_{0}$ needs to be reduced in order to compensate
for the increase with temperature of surface water vapor partial pressure. As
long as this change isn’t drastic, it is of little consequences for our
purpose, as we are not attempting to track the actual time evolution of an
atmosphere. It only means that the dry air mass in the equilibrium state is
somewhat different from what was specified for the initial condition of the
calculation. The magnitude of the adjustment will be quantified in Section
4.3.

When the surface saturation vapor pressure approaches the specified $p_{s}$,
$p_{0}\rightarrow 0$ and it is no longer possible to keep surface pressure
fixed as surface temperature is further increased. To deal with this case, we
introduce a pure-steam layer for $p>p_{s}$, extending to a greater surface
pressure $p_{s}^{\prime}$; we then time-step the radiative-convective model
only for $p<p_{s}$, computing the radiation from the pure steam layer assuming
it to lie on the pure-steam (i.e dewpoint) moist adiabat, as discussed in
Pierrehumbert & Ding (2016). Since water vapour is heavy and non-buoyant in a
hydrogen-dominated background atmosphere, it will remain at the bottom of the
atmosphere.

A radiative layer in the lower levels of the atmosphere tends to dry out the
upper atmosphere (due to its steep lapse rates), leading to a layered
structure with pure water vapor at the bottom, nearly pure hydrogen-helium at
the top, and a sharp transition layer between the two. All of the mass of
hydrogen and helium resides in the upper layer, and the opacity of this layer
(which in turn depends on the dry hydrogen-helium mass) strongly affects the
conditions for the surface temperature to enter a runaway state.

In cases which require a pure steam layer at the bottom, we set the bottom of
the atmosphere to a fixed surface temperature (which fixes the surface
pressure as the saturation vapour pressure at this temperature). The approach
in this case is to determine the instellation which is compatible with the
specified surface temperature, requiring multiple runs of the time-stepped
model. Unlike the usual inverse climate modeling approach, the required
instellation cannot be determined by just computing the OLR corresponding to a
given $T(p)$ profile, since the profile is affected by the instellation
through stellar absorption within the atmosphere. Instead, we need to guess an
instellation, time-step the model (subject to a lower boundary condition
provided by the steam layer) until $T(p)$ reaches equilibrium, and then check
the top of atmosphere balance. The instellation is then adjusted until top of
atmosphere balance is achieved. When the surface temperature is too high, the
OLR becomes decoupled from surface temperature owing to the optically thick
steam layer, and so equilibrium cannot be reached when instellation exceeds a
threshold value, which defines the runaway condition. For instellation above
the runaway threshold, the temperature increases until some process intervenes
to allow OLR to increase again, as discussed in Boukrouche et al. (2021);
Pierrehumbert (2023). In this paper, we do not compute the super-runaway
equilibrated state, but for a water-rich sub-Neptune with a deep water layer
it is certain to be above the critical point of water.

###  3.2 Experimental Procedure

To model the effect of different dry mass paths, we run the model with two
surface pressures – 1 bar and 10 bar. As described above, this surface
pressure is held constant until $q(T_{s},p_{s})=1$ is reached at the bottom of
the atmosphere, after which the surface temperature and pressure are increased
on the pure steam adiabat. Since steep radiative lapse rates keep the
atmosphere relatively dry at all pressures except very close to the surface,
this acts to keep the dry mass path of the atmosphere relatively constant at
$\approx$$1\text{\times}{10}^{4}\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-2}$
for the 1 bar initial condition and
$\approx$$1\text{\times}{10}^{5}\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-2}$
for the 10 bar case (we will quantify this in Section 4.3). Extending the
atmosphere on the pure steam adiabat does not add any dry mass to the
atmosphere (and instead assumes that the extra mass comes from evaporation
from the surface ocean). We will refer to these two cases as the “1 bar” and
“10 bar” cases interchangeably with
“$1\text{\times}{10}^{4}\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-2}$” and
“$1\text{\times}{10}^{5}\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-2}$”
cases.

For both the 1 bar and 10 bar atmospheres, we perform separate runs with the
G-star and M-star spectra described in Appendix A. For each set of runs, we
begin with a low instellation that gives surface temperatures between
$270\text{\,}\mathrm{K}300\text{\,}\mathrm{K}$ and run the model to radiative-
convective equilibrium. After the equilibrium state is found, we reinitialise
the model with an incrementally higher instellation, chosen as a balance
between numerical stability and computational efficiency. Increments of
$5\text{\,}\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$ and
$1\text{\,}\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$ were used for the 1 bar and
10 bar runs respectively. The initial temperature profile of the first run is
a dry adiabat with an isothermal stratosphere – runs at higher instellations
are then initialised on the final temperature profile of the previous run.
When $q(T_{s},p_{s})=1$ is reached at the bottom of the atmosphere (i.e.
$p_{\text{sat}}(T_{s})=p_{s}$), we switch to the procedure described above
where the surface temperature is held fixed and the instellation is iterated
until both local and global radiative equilibrium are attained. We run the
model until the surface temperature reaches 600 K, giving us a range of
surface temperatures between $270\text{\,}\mathrm{K}300\text{\,}\mathrm{K}$
and $600\text{\,}\mathrm{K}$.

###  3.3 Results

####  3.3.1 Temperature and Humidity Profiles

Figure 2: Sample of temperature-pressure profiles (left column) and specific
humidity profiles (right column) for the M star experiments with dry mass
paths of $10^{4}$ $\mathrm{kg}\text{\,}{\mathrm{m}}^{-2}$ (top row) and
$10^{5}$ $\mathrm{kg}\text{\,}{\mathrm{m}}^{-2}$ (bottom row). The
introduction of radiative layers in the lower atmosphere causes a sharp
increase in the surface temperature before the lower atmosphere becomes pure
steam. The radiative layers have lower lapse rates in the $10^{5}$
$\mathrm{kg}\text{\,}{\mathrm{m}}^{-2}$ case because less SW radiation
penetrates to the lower atmosphere. Figure 3: Same as Figure 2 but for the G
star experiments.

Figures 2 and 3 show sample temperature-pressure profiles and specific
humidity profiles for the M star and G star experiments. At pressure levels
where $q>q_{c}$, there is a sharp increase in lapse rate which corresponds to
increased surface temperatures with respect to convecting lower atmospheres.
Once the atmosphere reaches the pure steam limit in the lower atmosphere
($q=1$) the temperature profile follows a pure steam adiabat. In the M star 1
bar experiment and both G star experiments, enough SW radiation penetrates the
lower atmosphere to make the radiative layers have an extremely steep lapse
rate compared to both the moist adiabat and the pure steam adiabat. However,
in the M star 10 bar experiment (bottom row of Figure 2) attenuation of SW
radiation results in a smoother transition between adiabatic and radiative
regions.

####  3.3.2 Surface Temperature vs. Instellation

Figure 4: Surface temperature as a function of incoming instellation for
different hydrogen-helium inventories. The introduction of layers with
convective inhibition when the moisture content at the bottom of the
atmosphere (BOA) is high enough (orange region) causes a rapid increase in
surface temperature (and moisture content). Once the bottom of the atmosphere
is pure steam (pink region), the atmosphere becomes optically thick and
increasing surface temperature is no longer linked with an increase in
cooling, causing a runaway state.

Figure 4 shows the surface temperature, $T_{s}$, as a function of the incoming
instellation. We divide our graphs into three regions. In the region where
$q(T_{s},p_{s})<q_{c}$ (green in Figure 4), the atmospheric $T$-$p$ structure
is a moist adiabat in the lower atmopshere with a radiative stratosphere at
low pressures. Once $q(T_{s},p_{s})>q_{c}$ (orange in Figure 4), there is a
radiative layer in the lower atmosphere. The steep lapse rates in this region
tend to increase the surface temperature sharply with instellation than in the
lower temperature cases. Lastly, when $q(T_{s},p_{s})=1$ (pink in Figure 4),
the bottom of the atmosphere is pure steam and lies on a pure steam adiabat.
At this point the atmosphere becomes optically thick at all wavelengths and
the surface temperature decouples from the OLR. Further increases in the
instellation cannot be compensated by additional cooling and the surface
temperature increases until the ocean reservoir is depleted or the critical
point of water is reached. In the latter case (applicable to sub-Neptunes with
a significant water inventory), the hydrogen atmosphere is miscible with the
supercritical water envelope. The possible equilibrated states of a super-
runaway pure steam atmosphere are discussed in Pierrehumbert (2023). The
addition of hydrogen in this scenario is beyond the scope of this work.

When our initial condition is a 1 bar hydrogen-helium atmosphere, the M-star
runaway limit instellation is
$435\text{\,}\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$, and our G-star runaway
limit is $530\text{\,}\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$. The higher
runaway limit for G-stars is due to the Rayleigh scattering cross-section
being larger at shorter wavelengths, which leads to an increased SW albedo for
the G-star experiment where more of the instellation is at low wavelengths.

Interesting behaviour is caused by the sudden drop in maximum possible
instellation received in both cases when the atmosphere becomes pure steam (at
the boundary between the pink and orange regions in Figure 4). For global
equilibrium, we require:

$S=\frac{4\text{OLR}}{1-\alpha}$ (4)

where $\alpha$ is the albedo. Inspecting the relevant terms, the drop in
absorbed instellation at this boundary is caused by a sudden decrease in the
albedo with the introduction of the steam layer. This initial drop is caused
by a sudden increase in SW absorption from the sharp increase in water vapour
at the bottom of the atmosphere. Less radiation reflectsfrom the surface,
decreasing $\alpha$ and therefore $S$ in Equation 4 (assuming OLR remains
approximately constant). Adding more steam at the bottom of the atmosphere
(increasing surface temperature) eventually causes the albedo to increase
again, since the albedo of a thick pure steam layer is greater than the
surface albedo, 0.12. In the G-star case, more SW radiation penetrates into
the lower atmosphere than the M-star case, leading to a larger increase in
albedo and hence a larger increase in instellation.

The drop in absorbed SW radiation leads to a narrow range of instellations
with multiple equilibrium surface temperatures, some of which are unstable
(depending on the sign of $\derivative*{S}{T_{s}}$, (Koll & Cronin, 2019)). In
the M-star case with 1 bar H2-He, the drop in absorbed radiation leads to the
atmosphere abruptly entering a runaway state at any instellation above 0.32
$S_{0}$.

With 10 bars of H2-He as the initial condition, we see similar overall
behaviour as the 1 bar case. The boundaries between the three regimes
identified above have shifted to higher surface temperatures due to the higher
surface pressure. Due to the much higher CIA optical depth of the dry gas
inventory compared to the 1 bar case, there is more SW absorption in these
atmospheres. This mutes the varying albedo effect at the pure steam boundary
described above, and leads to the surface temperature increasing monotonically
with instellation towards the runaway limit. This limit is
$92\text{\,}\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$ for the G-star and
$116\text{\,}\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$ for the M-star case. We
note that the M-star limit is now higher than the G-star, despite the higher
albedos of the G-star irradiated atmospheres. In the M-star radiative regions
at $p\approx 10$ bar, there is very little SW radiation penetrating the
radiative layer. In the optically thick limit, radiative flux can be
approximated as radiative diffusion (Pierrehumbert, 2010; Heng et al., 2014),
with corresponding radiative equilibrium lapse rate:

$\partialderivative{T}{p}=\frac{3}{16}\frac{\kappa}{g\sigma
T^{3}}S_{\text{net}}$ (5)

where $\kappa$ is the Rosseland mean opacity and $S_{\text{net}}$ is the net
SW flux penetrating the region, which is also the flux that must be carried
upward through the region by radiative transfer. If $S_{\text{net}}$ is small
in the radiative region (as is the case with the M star 10 bar atmospheres),
then the lapse rate is also relatively small, meaning increasing $S$ has less
of an effect on the surface temperature than in cases where more radiation
penetrates into the deeper layers.

The curves in Figure 4 exhibit discrete stepping in regions where there are
superadiabatic radiative layers in the atmosphere. This numerical artefact
arises from the sensitivity of the surface temperature to the structure of the
radiative layer. Discrete stepping is worse if the vertical resolution is low
and there are few model levels in the radiative region. In this case, there is
a large jump in surface temperature when a new model level becomes inhibited
to convection. This behaviour motivated the use of higher vertical resolution
in the lower atmosphere (discussed in Section 3.1.3) which reduces the size of
the jumps in $T_{s}$ but does not completely remove them. The structure of the
stepping also changes as the instellation increment is increased or decreased,
suggesting the surface temperature may be somewhat dependent on the initial
state of the model (which is initialised from the final temperature profile of
the previous run). We do not believe the numerical artefacts affect the
conclusions of our work.

##  4 Discussion

###  4.1 Comparison of Runaway Limits to Canonical Values

Table 1 summarizes the runaway greenhouse instellations found in the previous
section. We can compare these to the classical runaway greenhouse limit
calculated for hydrogen atmospheres assumed to be on a simple moist adiabat
(e.g. Koll & Cronin, 2019). Since at high $T_{s}$ the moist adiabat approaches
the pure steam limit smoothly (Koll & Cronin, 2019), this is the same as
asking the maximum instellation a pure steam atmosphere can receive and remain
in global radiative equilibrium. The maximal OLR (sometimes called the
Simpson-Nakajima limit), is approximately
$280\text{\,}\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$ for Earth’s surface
gravity. In global equilibrium, this must be equal to $(1-\alpha)S/4$ by
Equation 4. Our estimation of the maximum instellation therefore depends on
the calculated albedo, $\alpha$, which varies by spectral type. We calculate
the runaway limit for a pure steam atmosphere to be
$1410\text{\,}\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$ and
$1150\text{\,}\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$ for our G-star and M-star
cases respectively. Comparing these numbers to our model results, for 1 bar of
solar H2-He mixture, the maximum instellation is less than half of the
Simpson-Nakajima limit and is less than 10% of the classical limit with 10
bars of H2-He mixture. This affects the placement of the inner-edge of the
habitable zone, which is given by:

Table 1: Summary of runaway instellations Expermiment | Runaway Limit[$\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$] ([$S_{0}$]) | HZ inner edge [AU]
---|---|---
G-star, 1 bar H2-He | 530 (0.389) | 1.60
M-star, 1 bar H2-He | 435 (0.320) | 0.280
G-star, 10 bar H2-He | 92 (0.0676) | 3.85
M-star, 10 bar H2-He | 116 (0.0852) | 0.543
Simpson-Nakajima, G-star | 1410 (1.04) | 0.982
Simpson-Nakajima, M-star | 1150 (0.847) | 0.172

$d=\quantity(\frac{L/L_{\odot}}{S/S_{0}})^{1/2}\text{AU}$ (6)

where $L/L_{\odot}$ is the luminosity of the star normalized by the solar
value. We take $L/L_{\odot}=1$ for the G-star and
$\log_{10}(L/L_{\odot})=-1.6$ for the M-star, the same as K2-18’s luminosity
(Benneke et al., 2019).

Our new inner-edge estimates are presented in Table 1 and are further from the
host star than previous calculations. For a 10 bar dry mass inventory, our
values lie outside the traditional outer edge of the habitable zone
(instellations between 0.2 and 0.4 $S_{0}$ (Kopparapu et al., 2013)). As a
reality check, we note that the low runaway thresholds are compatible with the
observed H2-dominated outer atmosphere of Uranus, even if the interior is
primarily composed of water, since the instellation of Uranus is only 0.0027
that of Earth – more than an order of magnitude lower than the runaway
threshold. Neptune is even further below the threshold.

###  4.2 Why Convective Inhibition Lowers the Runaway Limit

In this section we explore why, for a given instellation, the surface
temperature is much hotter in our experiments than in traditional calculations
of the inner edge of the habitable zone. Consider a planet with a given
instellation, $S$. The instellation roughly sets the stratospheric temperature
as $T_{\text{strat}}\sim(S(1-\alpha)/4\sigma)^{1/4}$ and the temperature of
the radiating layer where the characteristic LW optical depth, $\tau$, is
unity. With increasing pressure, the $T$-$p$ profile will follow a radiative
layer, followed by a moist adiabatic layer until it reaches the level where
$q=q_{c}$. In our simulations, the atmosphere then follows a radiative lapse
rate, which in general is much steeper than the equivalent moist adiabatic
lapse rate so long as the atmosphere is opaque enough and has enough SW
radiation penetrating to that level. This increased lapse rate leads to much
higher surface temperatures for equivalent levels of instellation. This effect
is illustrated in Figure 5(a). Once the bottom of the atmosphere becomes pure
steam, the surface temperature increases steeply with small increases in
instellation, since the bottom of the atmosphere becomes optically thick and
decoupled from the OLR.

Figure 5: Three $T$-$p$ profiles demonstrating the effect of superadiabatic
layers on the radiative balance of the atmosphere. (a) We compare an
atmosphere with convective inhibition (2) to one without (1), assuming an
identical instellation. Although the radiating temperature (and OLR) is the
same between the two cases, the atmosphere with the superadiabatic layer has a
higher surface temperature. (b) We again consider an atmosphere with (2) and
without (1) a superadiabatic layer, starting from a fixed surface temperature.
In this case, the superadiabatic layer causes the radiation temperature (and
therefore OLR) of the atmosphere to decrease, which in turn reduces the
maximum instellation it can receive. (c) We consider two atmospheres with
inhibited layers with different dry mass paths. The pressure at $\tau=1$ is
approximately constant between the two cases. The atmosphere with the greater
dry mass path, (2), has a greater average lapse rate between the surface and
the radiating level than (1), where the lower atmosphere is pure steam. The
radiating level is therefore colder, reducing the maximum OLR of the
atmosphere.

Equivalently, we can imagine the temperature-pressure profile in both the
classical and inhibited scenarios starting from the same surface temperature,
$T_{s}$. If $q(T_{s},p_{s})>q_{c}$ then our modelled atmospheres will follow a
steep radiative lapse rate, compared to the shallower moist adiabat. Once
$q<q_{c}$, the inhibited atmosphere will again follow a moist adiabat (albeit
one with a much steeper lapse rate than the one departing from $(T_{s},p_{s})$
owing to $q$ now being more dilute). The resulting upper atmosphere
temperature of the inhibited atmosphere will be much lower, leading to a much
lower OLR. In global equilibrium, $\text{OLR}=S(1-\alpha)/4$, and so the
maximum allowed instellation for a given surface temperature will be much
lower. This is illustrated in Figure 5(b).

###  4.3 Why Increased Dry Mass Lowers the Runaway Limit

It is clear from Figure 4 that the increased dry mass path in the 10 bar H2-He
experiments lowers the maximum instellation limit. Figure 5(c) shows two
different dry mass path atmospheres with the same surface temperature. The
pressure level at which $\tau=1$ is relatively constant between the two cases,
because in the upper atmosphere $q<<1$ so $\tau=1$ when $p(\tau=1)\approx
g/\kappa_{d}$ where $\kappa_{d}$ is some characteristic grey opacity of the
dry gas inventory. The temperature at which $q=q_{c}$ is approximately
constant with temperature (see Section 4.4). Increasing the dry mass path of
the atmosphere also increases the pressure level at which $q$ becomes unity
and the atmosphere transitions to a steam layer. As seen in Figure 5(c), since
the pure steam layer (pink) has a much lower lapse rate than the radiative
layer (orange), the average lapse rate between the surface and the $\tau=1$
level increases with increasing mass path. A higher average lapse rate
decreases the radiating temperature with increasing dry mass path, leading to
a drop in the OLR (and hence maximum instellation) in the runaway limit.

###  4.4 Discussion of Analytic OLR

In Appendix B, we calculate the OLR from an inhibited atmosphere to be:

$\text{OLR}=\Gamma(1+4\alpha)\quantity(\frac{\kappa_{d}m_{d}}{\bar{\theta}})^{-4\alpha}\quantity(q_{0}\frac{\kappa_{v}m_{d}}{\varepsilon\bar{\theta}})^{4/(\beta_{0}-1)}\sigma
T_{0}^{4}$ (7)

where $T_{0}\approx 273$ K, $\beta_{0}\equiv L/(R_{v}T_{0})$, $q_{0}\equiv
1/(\beta_{0}\varpi)$, $\alpha\equiv R_{d}/c_{p}$, $\bar{\theta}=3/5$,
$\varepsilon\equiv\mu_{v}/\mu_{d}$, $\Gamma$ is the standard gamma function
and $\kappa_{v}$ and $\kappa_{d}$ are the characteristic grey opacities of the
moist and dry components of the atmosphere respectively. Choosing
$\kappa_{v}=0.01$ ${\mathrm{m}}^{2}\text{\,}{\mathrm{kg}}^{-1}$ to match the
Simpson-Nakajima limit for a pure steam atmosphere, this leaves the OLR as a
function of the dry opacity, $\kappa_{d}$ and the dry mass path, $m_{d}$.
Table 2 shows estimates of the analytic OLR when $\kappa_{d}=$
$1.6\text{\times}{10}^{-4}\text{\,}{\mathrm{m}}^{2}\text{\,}{\mathrm{kg}}^{-1}$,
a sensible value for the dry opacity which corresponds to an H2-He atmosphere
that becomes optically thick at approximately 0.6 bars

The term in $m_{d}^{-4\alpha}$ represents the total optical depth of the dry
component of the atmopshere – as this term increases, the radiating
temperature drops and the OLR decreases. The second term in
$m_{d}^{4/(\beta_{0}-1)}$ represents how increasing the dry mass of the
atmosphere increases the temperature at which the atmosphere becomes
radiative, which increases the OLR. In general, $\alpha\gg(\beta_{0}-1)^{-1}$,
and we estimate $4(\beta_{0}-1)^{-1}-4\alpha\approx-8/7$, so the OLR decreases
with dry mass, explaining the trend in Table 1. Moreover, Table 2 shows that
the 10 bar M-star case is has a reduced dry mass compared to the equivalent
G-star case, explaining its higher instellation limit since from Equation 7
this atmosphere will cool more efficiently. More care could be taken to ensure
that the dry mass path of the atmosphere is conserved across our different
simulations. With the current model setup, this would involve iterating the
surface pressure so that the integral:

$\int_{0}^{p_{s}}(1-q)\frac{\differential{p}}{g}$ (8)

is conserved. This was deemed too computationally expensive for the current
study. Alternatively, one could implement a self-consistent moisture scheme
that keeps track of the mass of the vapour phase, which would naturally
conserve dry mass.

Moreover, in Appendix B we show that the ratio of this OLR to the Simpson-
Nakajima limit can be written as:

$\frac{\text{OLR}}{\text{OLR}_{\text{SN}}}=\frac{\Gamma(1+4\alpha)}{\Gamma(1+4/\beta_{0})}\quantity(q_{0}\frac{\kappa_{v}m_{d}}{\varepsilon\bar{\theta}})^{4/(\beta_{0}-1)}\quantity(\frac{\kappa_{d}m_{d}}{\bar{\theta}})^{-4\alpha}$
(9)

Since the ratio of gamma functions is of order unity, and
$4(\beta_{0}-1)^{-1}\ll 1$, so long as $\kappa_{d}m_{d}/\bar{\theta}>1$, i.e.
the dry inventory of the atmosphere is optically thick, we should expect the
OLR to be lower than the classical Simpson-Nakajima limit.

The analytic expression provides a relatively good estimate of the OLR but
decreases slightly too steeply with increasing $m_{d}$ – more experiments at
different dry paths would need to be run to establish the limitations of the
power-law formulation.

Table 2: Comparison of model OLR with analytical calculations with a dry gas opacity of $1.6\text{\times}{10}^{-4}\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-2}$ Expermiment | OLR [$\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$] | Dry mass [$\mathrm{kg}\text{\,}{\mathrm{m}}^{-2}$] | Analytic OLR [$\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$]
---|---|---|---
G-star, 1 bar H2-He | 100 | $9.92\text{\times}{10}^{4}$ | 108
M-star, 1 bar H2-He | 102 | $9.01\text{\times}{10}^{4}$ | 109
G-star, 10 bar H2-He | 13.8 | $9.74\text{\times}{10}^{3}$ | 11.8
M-star, 10 bar H2-He | 26.4 | $5.99\text{\times}{10}^{4}$ | 19.9

###  4.5 What Does a Super-Runaway State Look Like?

Having discussed how the runaway greenhouse threshold changes for inhibited
atmospheres, it is natural to wonder what a super-runaway atmosphere would
look like. For Earth-like planets with a finite water reservoir, eventually
the water inventory will be entirely in the atmosphere, and the lower
atmosphere will not be in liquid-vapour phase equilibrium, allowing it to lie
on a dry adiabat and increase its OLR (Boukrouche et al., 2021). However, for
a Hycean world with an almost limitless supply of water this is not possible.
As discussed in Pierrehumbert (2023), a super-runaway state is likely to
consist of supercritical water vapour in the lower atmosphere. Since the
supercritical phase is not constrained to lie on the phase equilibrium
boundary between liquid and vapour (in contrast to the condensing layers
above), a deep layer of the interior heats up until the supercritical water
layer penetrates to high enough altitudes that radiation to space can increase
beyond the runaway limit. In the pure water case discussed in Pierrehumbert
(2023), this generally leaves a thin condensing region near the top of the
atmosphere, but in the case with a substantial H2 layer at the top, the
radiating level is in the H2 layer, so the condensing layer is eliminated
entirely. The warming proceeds until the H2 layer becomes hot enough to
increase the OLR, but once supercritical water is in contact with the H2, that
layer would mix into the supercritical water and largely disappear, because H2
(and presumably also He) is completely miscible in supercritical water
(Soubiran & Militzer, 2015). A full thermal evolution model would be needed to
determine how long this process would take. A curious possibility emerges
because the runaway instellation threshold with an H2 layer is considerably
below the pure-steam limit. If the instellation lies between the two
thresholds, once the H2 layer is diluted into the supercritical water
interior, a liquid ocean could form again, forcing some H2 back into the
atmosphere. One possibility is that the mixing between the layers results in
just enough H2 remaining in the outer layer for radiative balance to be
achieved, predicting a self-regulation in the thickness of the H2 layer.

The scenario we have modelled in this paper corresponds to a cold-start, in
which the planet begins in a sub-runaway state and then undergoes a runaway as
the stellar luminosity increases. This is a possible scenario for an F or G
star, but an alternate scenario for the evolution is a hot-start, in which the
planet begins in a super-runaway state, either because of the heat of
formation of the planet or because of intense illumination in the extended pre
main-sequence stage of low mass stars. In a hot-start, the initial H2-H2O
inventory would most likely begin in a supercritical mixed state. While the
instellation remains above the runaway threshold, a significant separated H2
layer would never form (unless H2 is a major proportion of the initial
composition). Once the planet cools down enough to be sub-runaway, though, a
liquid ocean will form, leading an H2 layer to effervesce out of the
subcritical liquid ocean.

###  4.6 Mixing of H2-He into the Pure Steam Layer

The layered structure that occurs in our model at high surface temperatures,
with a pure steam layer at the bottom, is a self-consistent solution of the
radiative-convective equations. It is a peculiarity of the compositional
stability criterion (Equation 2a) that when $q>q_{c}$, $q\rightarrow 1$
appears to be a singular limit. When $q=1$ exactly, the moist stability
criterion is the usual criterion that the lapse rate be steeper than the
adiabat. However, if even an infinitesimal amount of hydrogen or helium mixes
into the pure steam layer, so $q=1-\delta$, with $\delta\ll 1$, the
compositional stability criterion then nominally applies according to which
lapse rate steeper than the moist adiabat are stable. Radiative cooling would
then be expected to generate a steep radiative layer within the nearly pure
steam layer, no matter how small $\delta$ may be. However, since the
stabilizing compositional buoyancy becomes exceedingly week for small
$\delta$, many other mixing processes could intervene, so we find the
generation of a radiative layer under these circumstances to be implausible.
We cannot rule out the possibility, though; it is a matter that will need to
be resolved by future resolved-convection modelling.

###  4.7 Implications for Observations

Our main result is that the runaway limit for sub-Neptune water worlds is
greatly reduced with even 1 bar of hydrogen. Due to observational biases
favouring the detection of low semi-major axis planets, most of the current
observational candidates for Hycean worlds have relatively high equilibrium
temperatures. From Table 1, the inner edge of the habitable zone is around
0.28 AU for our formulation. Table 1 of Madhusudhan et al. (2021) lists
potential Hycean world candidates, all of which lie within our inner edge
estimate and therefore would only be able to sustain a liquid water ocean if
hydrogen were not detected in large abundances in the atmosphere. The well-
studied sub-Neptune K2-18 b, originally thought to lie right on the inner edge
of the classical habitable zone, is well beyond the inner edge by this
measure.

In this case, the planets most likely to host liquid water oceans on close-in
orbits are either terrestrial planets with a higher mean-molecular weight
background gas, or pure “water-worlds” with little to no H2 envelope and an
atmosphere predominantly composed of steam. The latter type of planets and
their evolution have been studied previously in the highly-irradiated regime
(Mousis et al., 2020; Aguichine et al., 2021), finding that water worlds may
fit the observed mass-radius distribution of small radius planets. The super-
Earth sized planet Kepler-138 d (Piaulet et al., 2023) is a candidate
volatile-rich planet. It has a 1.5 $R_{\oplus}$ radius and low density making
it not dense enough to be predominantly rocky, but also too dense to sustain a
significant hydrogen envelope that would not be lost through atmospheric
escape. Although this particular planet is above the runaway greenhouse
instellation threshold, similar cool planets may be able to host liquid water
oceans at near-Earth instellations.

The ability to observationally distinguish liquid water surfaces and super-
runaway mixtures of H2-H2O would allow us to verify some of the predictions of
this paper. Mapping the transition from sub-runaway planets to super-runaway
planets as a function of instellation would give us the runaway instellation
limit, which could be compared to the predictions in Table 1 to provide
evidence for or against robust super-adiabatic layers. The combined non-
detection of ammonia and detection of methanol in a sub-Neptune atmosphere has
been proposed as a method of distinguishing a shallow water surface (Tsai et
al., 2021). However, this method requires around 20 transits with the James
Webb Space Telescope which may prove unfeasible given time allocation
restrictions. Our results also predict a sharp transition between sub-runaway
atmospheres where the upper atmosphere is very dry (see Figures 2 and 3) to
super-runaway atmospheres that are moist due to the mixing of H2 and
supercritical H2O (Pierrehumbert, 2023). In contrast, the moistening of the
upper atmosphere in the classical runaway greenhouse limit without convective
inhibition is more smooth as instellation is increased and would occur at
higher instellations. Further complications arise from the presence of non-
Hycean sub-Neptunes, where the presence of water in the atmosphere is not
necessarily correlated with its interior structure or the presence of a
surface. More work is needed to understand how we can observationally
disentangle the various possible atmospheric structures of habitable zone sub-
Neptunes.

###  4.8 Robustness of Calculations and Caveats

####  4.8.1 Day-Night Averaging

Our model is one-dimensional, and makes the assumption that the stellar
radiation is redistributed evenly over the dayside and nightside of the
planet. The superadiabatic layers in our model are sustained by the need to
remove the stellar flux deposited at the surface of our model. If day-night
heat redistribution is not efficient, then the nightside of these planets may
be able to sustain shallower lapse rates which would aid with radiative
cooling. However, for thick hydrogen atmospheres on temperate sub-Neptunes,
general circulation models (Charnay et al., 2021; Innes & Pierrehumbert, 2022)
have shown that the combination of slow rotation rate and low mean molecular
weight atmospheres produces globally weak temperature gradients thanks to
dynamical redistribution of heat. This suggests heat deposited near the
surface would be transported horizontally to the nightside, maintaining the
steep lapse rates globally. Scaling relations for shallow atmospheres also
suggest low MMW atmospheres should have efficient redistribution of heat
(Koll, 2022).

####  4.8.2 Assumption of Saturation

Our model also assumes 100% relative humidity in the column (except above the
stratospheric cold trap). One major effect of 3D dynamics is to cause
subsiding regions (e.g. due to the descending branch of a Hadley-like
circulation, or night side subsidence on a tidally locked exoplanet).
Descending dry air causes compressional heating and undersaturation, and can
be responsible for regions where the OLR is locally greater the result of a
globally-averaged calculation (Pierrehumbert & Swanson, 1995; Leconte et al.,
2013). If an atmosphere is undersaturated to the point it lies below the
critical water vapour mixing ratio $q_{c}$, then superadiabatic layers
responsible for surface heating may not form. Moreover, on Earth moist
convection is the main mechanism by which water vapour is transported
vertically in the atmosphere. Within the inhibited superadiabatic layer,
mixing by convection is suppressed and our assumption of 100% relative
humidity may break down aloft. However, we note that if the near-surface
layers are saturated (due to being close to the ocean surface) and the layers
aloft are undersaturated, this induces an even greater mean-molecular weight
gradient to stabilise the atmosphere to convection than before. Moreover,
since undersaturated lofted parcels would travel on the dry adiabat, which has
a steeper lapse rate than the moist adiabat, this would again help stabilise
the atmosphere to convection. One could argue that decreasing relative
humidity with height would affect the radiative calculations. However, from
Section 4.3 and Appendix B we can see that the main driver of lower cooling is
the radiative effect of the dry mass of the atmosphere. A decrease in relative
humidity with height would likely decrease the thickness of the superadiabatic
layer, in which case the reduction in OLR may not be as severe as in the fully
saturated scenario.

####  4.8.3 Clouds and Hazes

Cloud and haze opacities were not included in our model. We expect their
introduction to affect our results in three main ways. Firstly, the increase
in LW opacity due to cloud water or hazes would exacerbate the greenhouse
effect and, if taken independently from other cloud and haze radiative
feedbacks, decreases the runaway instellation. Secondly, the SW scattering
properties of clouds and hazes increase the effective planetary albedo. This
opposes the greenhouse effect and increases the value of the runaway
instellation. Thirdly, clouds and hazes may reduce the magnitude of shortwave
radiation penetrating the lower atmosphere. This would reduce the radiative
lapse rates in the inhibited layers by reducing $S_{\text{net}}$ in Equation
5. The reduction of radiative lapse rates will decrease the surface
temperature for any given instellation, raising the runaway instellation. A
similar effect was modelled in Piette & Madhusudhan (2020), who demonstrated
that surface water oceans were possible on K2-18 b if the haze scattering
opacity was high enough. In this case, the lower atmosphere becomes
isothermal, allowing for temperate oceans at high pressures. However, their
model neglected the effect of convective inhibition. Moreover, as discussed in
Section 3.1.3, saturated subadiabatic radiative layers are unstable to
convection when $q>q_{c}$, implying that above the critical moisture
threshold, the moist adiabat is the minimum possible lapse rate.

If the cooling effects of clouds and hazes dominate their potential warming
effect, then our runaway instellations in Table 1 are likely too pessimistic
and the inner edge of the habitable zone could be at lower orbital distances.
The magnitude of the cloud and haze radiative effects is likely to be strongly
dependent on the particles’ microphysical properties and 3D spatial
distribution (e.g., Yang et al., 2013; Turbet et al., 2021). These effects are
beyond the scope of our simplified 1D model, though we encourage future
efforts to quantify the impact of clouds and hazes on our results.

####  4.8.4 Other Heat Transport Mechanisms

We also need to consider other mechanisms which may be able to transport heat
through the stabilised layers. For example, how efficient is thermal
conduction at transporting flux deposited in the lower atmosphere? We can
compare the efficiency of thermal conduction to radiation by comparing the
thermal diffusion coefficient for conduction and radiation, as in Markham et
al. (2022). For an ideal gas, the thermal diffusivity of conduction is:

$k_{\text{cond}}=\rho\lambda c_{v}\sqrt{\frac{2k_{B}T}{\pi
m}}\approx$1\text{\,}\mathrm{W}\text{\,}\mathrm{m}\text{\,}{\mathrm{K}}^{-1}$$
(10)

where $\lambda$ is the mean free path, $m$ is the weight of an average gas
molecule, $k_{B}$ is the Boltzmann constant and $c_{v}$ its specific heat
capacity at constant volume. We choose characteristic values to give an upper
limit on $k_{\text{cond}}$. We use $T=300$ K, $\lambda=k_{B}T/(\sqrt{2}\pi
d^{2}p)$ with $d=290$ pm (the kinetic diameter of a hydrogen molecule, Mehio
et al., 2014) and $m=2$ amu. Density is calculated using the ideal gas law and
we note the final result is independent of pressure.

The radiative diffusivity of an optically thick gas is approximately

$k_{\text{rad}}=\frac{16}{3}\frac{\sigma
T^{3}}{\kappa\rho}\approx$1\text{\times}{10}^{4}\text{\,}\mathrm{W}\text{\,}\mathrm{m}\text{\,}{\mathrm{K}}^{-1}$$
(11)

where we have used $\kappa=0.01$ $\mathrm{kg}\text{\,}{\mathrm{m}}^{-2}$ and
used the same $T$ and $\rho$ in our calculation of $k_{\text{cond}}$. We
conclude that energy transport via radiation is much more efficient than
thermal conduction. Other sources of heat transport that could be considered
are advective heat transports. Although the superadiabatic regions are
statically stable, eddy heat transport could play a role in transporting heat
vertically, especially if the vertical wind shear is high. Moreover, we have
also neglected latent heat fluxes, which could play a significant role if
there is significant condensation or re-evaporation of condensates around the
region of interest. Within the framework of a 1D model with no dynamics, it is
very difficult to get an accurate estimate of the magnitude of these fluxes.
Studying this system with a cloud-resolving model (e.g. Lefèvre et al., 2021;
Tan et al., 2021) is key to understanding the robustness of the superadiabatic
layer against other mechanisms of heat transport. These models would also
indicate the radiative effect of water clouds on the system and derived
runaway greenhouse limits.

####  4.8.5 Possibility of Multiple Equilibria

Lastly, apart from in the region near the runaway limit (see Section 3.3), our
model doesn’t consider the possibility of multiple equilibrium states and
hysteresis. Given that the phase structure of a water-world can change
drastically on either side of the runaway limit from steam above a surface
ocean to a supercritical envelope (Pierrehumbert, 2023), there is a
possibility that much warmer surface temperatures could be achieved with a
similar instellation if we modelled the atmosphere with a supercritical water
layer mixed with hydrogen gas. Our models represent a “cold start”, i.e.
warming a planet up that initially starts with a surface water ocean. However,
realistically a sub-Neptune will form hot and cool down from a state where the
water is supercritical (Misener & Schlichting, 2022; Markham et al., 2022). If
there are multiple equilibrium solutions, it is possible that even when the
instellation lies below the runaway threshold, the water will still be in a
supercritical state. Moreover, in our model we have neglected the effect of
internal heating from the residual heat of formation or tidal heating. Given
that our model only requires
$\approx$1\text{\,}\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$$ of solar flux
penetrating the lower layers to drive steep superadiabatic lapse rates, a
similar level of internal flux could equally sustain very high surface
temperatures in the absence of a significant instellation.

##  5 Conclusions

The aim of this study was to determine the sustainabilty of a liquid water
ocean on a Hycean world with a significant H2-He inventory. Our major findings
are:

  1. 1.

Neglecting water vapour feedbacks, 10-20 bars of solar H2-He mixture will
drive a surface ocean supercritical when forced with solar instellation. A
planet receiveing 10 times solar instellation would have to have less than 1
bar of hydrogen to sustain a liquid water ocean.

  2. 2.

Including water vapour feedbacks, the presence of superadiabatic layers where
convection is inhibited in the lower atmosphere reduces the runaway greenhouse
instellation limit from the Simpson-Nakajima limit significantly. For a solar
H2-He inventory of around
$1\text{\times}{10}^{4}\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-2}$, the
runaway greenhouse limit to an instellation of approximately
$530\text{\,}\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$ for a G-star and
$435\text{\,}\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$ for an M-star. This reduces
further to around $100\text{\,}\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$ for an
H2-He inventory of around
$1\text{\times}{10}^{5}\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-2}$.

  3. 3.

The reduced instellation limits correspond to moving the inner edge of the
habitable zone to around 1.6 AU (3.85 AU) for a planet orbiting a G-star with
1 bar (10 bar) of H2-He and equivalently 0.280 AU (0.543 AU) for a planet
orbiting an M-star (c.f. 0.982 AU and 0.172 AU for a G-star and M-star
respectively from previous models).

  4. 4.

Analytical models of the OLR show the key parameter responsible for the
reduction in the OLR is the total optical depth of the dry inventory, given
that steep superadiabatic lapse rates in the inhibited layers dry the
atmosphere aloft. A higher dry optical depth reduces the radiating temperature
of the atmosphere and caps the maximum cooling from a H2-H2O atmosphere. If we
model the atmosphere as having a constant, gray opacity for the dry gases,
then the limiting OLR scales roughly as the inverse of the dry mass path.

  5. 5.

Our results suggest that most of the current Hycean world targets are within
the inner limit of the habitable zone and unlikely to host liquid water
oceans. The most promising targets for observing a liquid water ocean on a
close-in orbit are therefore traditional terrestrial-like planets with a high-
mean molecular weight background atmosphere or “water worlds” with negligible
H2-He envelopes.

We conclude by encouraging the use of 3D cloud resolving models to study the
robustness of the inhibited, superadiabatic radiative layers to 3D dynamics
and other sources of heat flux.

This paper is supported by funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(Grant agreement No. 740963). S.-M.T. thanks Matej Malik and Daniel Kitzmann
for the discussion on the water continuum absorption. S.-M.T. acknowledges
support from the University of California at Riverside and NASA Exobiology
grant 80NSSC20K1437.

##  Appendix A Calculation of radiative fluxes

We calculate the radiative fluxes using the SOCRATES code, based on Edwards &
Slingo (1996). SOCRATES uses the correlated-$k$ method to efficiently
calculate two-stream fluxes in both the longwave (LW) and shortwave (SW)
regions of the spectrum. Gaseous overlap is treated with the equivalent
extinction method with resorting and rebinning (Lacis & Oinas, 1991). In
Section 2, we calculate LW fluxes in 300 bands equally spaced in wavenumber
between $1\text{\,}{\mathrm{cm}}^{-1}$ and $5000\text{\,}{\mathrm{cm}}^{-1}$.
In the SW region, we perform two separate calculations for a G-type star
(using the Lean & DeLand (2012) solar spectrum) and an M-type star (calculated
using PHOENIX (Husser et al., 2013) for a $3500\text{\,}\mathrm{K}$ star with
$\log(g)=5.0$ and solar metallicities and alpha element abundances). In
Section 2 we use 300 bands equally spaced in wavenumber in each of the LW and
SW regions of the spectrum. For the calculations in Section 3 we use 30 bands
in each of the LW and SW regions to speed the convergence of the numerical
iteration. We found no significant differences in the temperature-pressure
profiles in benchmark models run with 300 bands in each region compared to 30
bands, justifying this reduction in wavenumber resolution. For the M-star we
use the range $250\text{\,}{\mathrm{cm}}^{-1}$ to
$40\,000\text{\,}{\mathrm{cm}}^{-1}$ and for the G-star we use
$250\text{\,}{\mathrm{cm}}^{-1}$ to $50\,000\text{\,}{\mathrm{cm}}^{-1}$. We
consider the collision-induced absorption (CIA) due to H2-H2 and H2-He
interactions as the main sources of absorption. We calculate H2-He opacities
using the HITRAN database (Karman et al., 2019) and use the HITRAN database
with additional data from Borysow (2002) for the calculation of H2-H2
opacities. In the SW calculation we include the effects of Rayleigh scattering
by both hydrogen and helium, calcualting the cross sections using fits of
refractive indices taken from Peck & Huang (1977) and Cuthbertson &
Cuthbertson (1932) respectively. A SW surface albedo of 0.12 was specified
(Goldblatt et al., 2013) and the surface temperature was assumed to be
identical to the temperature of the lowest layer of the atmosphere.

Water vapour absorption $k$-coefficients are calculated using HITRAN data
(Karman et al., 2019) and continuum absorption is calculated with the MT_CKD
model (Mlawer et al., 2012). The refractive index of water used to calculate
the Rayleigh scattering coefficients was taken from Ciddor (1996).

##  Appendix B Analytic OLR calculations

We can understand our results better by constructing a simple grey analytic
framework for calculating the maximum OLR of our model atmospheres.

In the optically thick limit ($\tau\gg 1$), the OLR of a grey atmosphere is
given by:

$\text{OLR}\approx\int_{0}^{\infty}\sigma
T(\tau)^{4}e^{-\tau}\differential{\tau}$ (B1)

To calculate this integral, we require the temperature profile $T(\tau)$ and
the optical depth:

$\tau(p)=\int_{0}^{p}\frac{\differential{p^{\prime}}}{\bar{\theta}g}\quantity(\kappa_{d}(1-q)+q\kappa_{v})$
(B2)

where $\kappa_{d}$ and $\kappa_{v}$ are characteristic grey opacities for the
dry and condensible phases of the atmosphere respectively and $\bar{\theta}$
is the average zenith angle, accounting for the angular distribution of the LW
radiation. To simplify the problem, we split our atmospheres into two regions
separated at $q=q_{c}=1/(\beta\varpi)$, since we expect the lapse rates in
these two regions to be very different. We assume that the moisture at this
level is dilute such that we can make the approximation:

$q\approx\varepsilon\frac{p_{\text{sat}}}{p}$ (B3)

This approximation is approximately valid if $q_{c}<0.1$, i.e.
$T<0.1L\varpi/R_{v}=480$ K assuming $\beta=L/R_{v}T$. We can verify that this
is the case in all of our simulations. We then approximate $p_{\text{sat}}(T)$
in the same way as Koll & Cronin (2019), writing:

$p_{\text{sat}}(T)=p_{0}\quantity(\frac{T}{T_{0}})^{\beta_{0}}$ (B4)

where $\beta_{0}\equiv L/(R_{v}T_{0})$ and $(p_{0},T_{0})$ is some reference
point on the vapour-liquid phase curve of water chosen to be close to our
region of interest. Like in Koll & Cronin (2019), we will choose $p_{0}$ as
the pressure where a pure steam atmosphere has unity optical thickness:

$p_{0}=\frac{g\bar{\theta}}{\kappa_{v}}$ (B5)

This gives $(p_{0},T_{0})=(588.6~{}\text{Pa},272.6~{}\text{K})$ for $g=9.81$
$\mathrm{m}\text{\,}{\mathrm{s}}^{-2}$, $\kappa_{v}=0.01$
${\mathrm{m}}^{2}\text{\,}{\mathrm{kg}}^{-1}$, $\bar{\theta}=3/5$ \-
remarkably close to the triple point of water vapour. If Equation B3 holds,
then the temperature and pressure $(T_{*},p_{*})$ at which $q=q_{c}$
satisfies:

$\displaystyle\varepsilon\frac{p_{0}}{p_{*}}\quantity(\frac{T_{*}}{T_{0}})^{\beta_{0}}$
$\displaystyle=\frac{R_{v}T_{0}}{p}$ (B6) $\displaystyle T_{*}$
$\displaystyle=T_{0}\quantity(\frac{1}{\beta_{0}\varpi}\frac{p_{*}}{\varepsilon
p_{0}})^{1/(\beta_{0}-1)}$ (B7)

By noting that $\beta_{0}\gg 1$ in our range of $T_{0}$ (Koll & Cronin, 2019),
we can readily verify our dilute approxmation:

$\frac{T_{*}}{L\varpi/R_{v}}\approx\frac{1}{\beta_{0}\varpi}=0.06$ (B8)

for $T_{0}=272.6$ K. We will also define $q_{0}=1/(\beta_{0}\varpi)$ as the
moisture inhibition threshhold at $T=T_{0}$

We then want to relate our pressure $p_{*}$ to the dry mass path. We note that
since $q_{c}$ is relatively dilute, the atmosphere will quickly dry on the
moist adiabat extending upwards from this point, leaving most of the upper
atmosphere dry. For a dry mass path $m_{d}$, we can then write $p_{*}\approx
m_{d}g$.

We assume that the temperature structure of the upper atmosphere ($p<p_{*}$)
is approximately a dry adiabat emanating from $(p_{*},T_{*})$. This neglects
the effect of moisture on the lapse rate, which can be significant but quickly
leads to intractable solutions since lapse rate $\derivative*{T}{p}$ depends
on $q(p,T)$. This dry adiabat has the form:

$T=T_{*}\quantity(\frac{p}{p_{*}})^{\alpha},\quad\alpha\equiv\frac{R_{d}}{c_{p}}$
(B9)

This allows us to write an equation for how $\tau$ varies with $\ln T$:

$\displaystyle\derivative{\tau}{\ln T}$
$\displaystyle=\quantity(\derivative{\ln T}{\ln p})^{-1}p\derivative{\tau}{p}$
(B10)
$\displaystyle=\alpha^{-1}\quantity[\frac{\kappa_{d}p}{\bar{\theta}g}+(\kappa_{v}-\kappa_{d})\frac{\varepsilon
p_{\text{sat}}}{\bar{\theta}g}]$ (B11)
$\displaystyle=\alpha^{-1}\quantity[\frac{\kappa_{d}p_{*}}{\bar{\theta}g}\quantity(\frac{T}{T_{*}})^{1/\alpha}+(\kappa_{v}-\kappa_{d})\frac{\varepsilon
p_{0}}{\bar{\theta}g}\quantity(\frac{T}{T_{0}})^{\beta_{0}}]$ (B12)

Integrating this relation and letting $\tau=0$ when $(T,p)=(0,0)$ yields:

$\tau=\alpha^{-1}\quantity[\frac{\alpha\kappa_{d}p}{\bar{\theta}g}+(\kappa_{v}-\kappa_{d})\frac{p}{\beta_{0}\bar{\theta}g}q]$
(B13)

To calculate the OLR, ideally we would invert Equation B13 to find $T(\tau)$
in Equation B1. However, due to the mixed powers of $T$ in Equation B13 (one
term in $T^{1/\alpha}$ and one in $T^{\beta_{0}}$), this cannot be done
analytically. To proceed, we compare the magnitude of the two terms and argue
that the second term can be neglected so long as:

$q\ll\frac{\beta_{0}\alpha}{\kappa_{v}/\kappa_{d}-1}$ (B14)

Let us take characteristic values of $\beta_{0}$ at $T_{0}=273$ K,
$\alpha\approx 2/7$, and $\kappa_{v}=0.01$
${\mathrm{m}}^{2}\text{\,}{\mathrm{kg}}^{-1}$ and $\kappa_{d}=$
$1.6\text{\times}{10}^{-4}\text{\,}{\mathrm{m}}^{2}\text{\,}{\mathrm{kg}}^{-1}$
(this value of $\kappa_{d}$ yields good agreement between our final analytical
OLR and simulations and is consistent with H2 becoming optically thick in the
infra-red between 0.1 and 1 bar). In this case, Equation B14 gives:

$q\ll 0.1$ (B15)

We assert this to be true around the $\tau\sim 1$ region of the atmosphere,
even if it doesn’t strictly hold at $p=p_{*}$.

We then have:

$T=T_{*}\quantity(\frac{\tau}{\tau_{*}})^{\alpha},\quad\tau_{*}\equiv\frac{\kappa_{d}p_{*}}{\bar{\theta}g}=\frac{\kappa_{d}m_{d}}{\bar{\theta}}$
(B16)

which when inserted into B1 yields

$\text{OLR}=\Gamma(1+4\alpha)\tau_{*}^{-4\alpha}\sigma T_{*}^{4}$ (B17)

where $\Gamma$ is the standard gamma function. Expanding $T_{*}$ from Equation
B7 and $\tau_{*}$ from Equation B16 yields:

$\text{OLR}=\Gamma(1+4\alpha)\quantity(\frac{\kappa_{d}m_{d}}{\bar{\theta}})^{-4\alpha}\quantity(q_{0}\frac{\kappa_{v}m_{d}}{\varepsilon\bar{\theta}})^{4/(\beta_{0}-1)}\sigma
T_{0}^{4}$ (B18)

We immediately see that this OLR limit does not depend on the surface
temperature, which is characteristic of a runaway greenhouse atmosphere. The
first bracket corresponding to $\tau_{*}$ represents the effect of increasing
the dry opacity of the atmosphere. This shifts the radiating temperature up
the adiabat with exponent $\alpha$, and therefore reduces the OLR of the
atmosphere. The second bracket traces back to Equation B7 and represents how
increasing the dry mass of the atmosphere increases the temperature at which
the atmosphere becomes inhibited, $T_{*}$. Increasing the base temperature of
the adiabat is associated with an increase in the OLR, albeit with a weak
dependence of $m_{d}^{4/(\beta_{0}-1)}$. We can see that the OLR depends on
the dry mass path $m_{d}$ with an exponent of:

$4(\beta_{0}-1)^{-1}-4\alpha\approx-4\alpha$ (B19)

since $\beta_{0}\gg 1$. For a diatomic ideal gas $-4\alpha\approx-8/7$, so we
would expect the OLR to drop of sharply with dry mass path.

We can compare Equation B18 with Equation 24 in (Koll & Cronin, 2019), for a
pure steam atmosphere, which was an estimate of the classical Simpson-Nakajima
limit:

$\text{OLR}_{\text{SN}}=\Gamma(1+4/\beta_{0})\quantity(\frac{\kappa_{v}p_{0}}{\bar{\theta}g})^{-4/\beta_{0}}\sigma
T_{0}^{4}=\Gamma(1+4/\beta_{0})\sigma T_{0}^{4}$ (B20)

Taking the ratio of this equation and Equation B18 gives:

$\frac{\text{OLR}}{\text{OLR}_{\text{SN}}}=\frac{\Gamma(1+4\alpha)}{\Gamma(1+4/\beta_{0})}\quantity(q_{0}\frac{\kappa_{v}m_{d}}{\varepsilon\bar{\theta}})^{4/(\beta_{0}-1)}\quantity(\frac{\kappa_{d}m_{d}}{\bar{\theta}})^{-4\alpha}$
(B21)

In general, since $4/(\beta_{0}-1)\ll 1$ and the gamma function ratio is of
order unity, so long as $\kappa_{d}m_{d}/\bar{\theta}>1$ (i.e. the dry mass of
the atmosphere is optically thick), the OLR limit will be lower than the
classical Simpson-Nakajima limit.

Lastly, if we note that if one were to take the opposite limit of Equation B14
(i.e. the moist optical depth at $\tau\sim 1$ dominates the OLR, we would
find:

$\text{OLR}=\Gamma(1+4/\beta_{0})\quantity(\frac{c_{p}}{LT_{0}})^{-4/\beta_{0}}\sigma
T_{0}^{4}$ (B22)

which is identical to the “dilute limit” found in Koll & Cronin (2019). This
limit does not depend on the dry mass path and is only moderately lower than
the Simpson-Nakajima limit. Since our results vary greatly with dry mass path
and are much lower than this limit, this should reassure us that Equation B14
is a good assumption.

## References

  * Aguichine et al. (2021) Aguichine, A., Mousis, O., Deleuil, M., & Marcq, E. 2021, The Astrophysical Journal, 914, 84, doi: 10.3847/1538-4357/abfa99
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, Annual Review of Astronomy and Astrophysics, 47, 481, doi: 10.1146/annurev.astro.46.060407.145222
  * Batalha (2014) Batalha, N. M. 2014, Proceedings of the National Academy of Sciences of the United States of America, 111, 12647, doi: 10.1073/pnas.1304196111
  * Benneke et al. (2019) Benneke, B., Wong, I., Piaulet, C., et al. 2019, The Astrophysical Journal, 887, L14, doi: 10.3847/2041-8213/ab59dc
  * Borysow (2002) Borysow, A. 2002, Astronomy and Astrophysics, 390, 779, doi: 10.1051/0004-6361:20020555
  * Boukrouche et al. (2021) Boukrouche, R., Lichtenberg, T., & Pierrehumbert, R. T. 2021, arXiv.  http://arxiv.org/abs/2107.14150
  * Charnay et al. (2021) Charnay, B., Blain, D., Bézard, B., et al. 2021, Astronomy and Astrophysics, 646, doi: 10.1051/0004-6361/202039525
  * Ciddor (1996) Ciddor, P. E. 1996, Applied Optics, 35, 1566, doi: 10.1364/AO.35.001566
  * Cloutier et al. (2017) Cloutier, R., Astudillo-Defru, N., Doyon, R., et al. 2017, Astronomy and Astrophysics, 608, A35, doi: 10.1051/0004-6361/201731558
  * Cloutier et al. (2019) —. 2019, Astronomy and Astrophysics, 621, doi: 10.1051/0004-6361/201833995
  * Cuthbertson & Cuthbertson (1932) Cuthbertson, C., & Cuthbertson, M. 1932, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 135, 40, doi: 10.1098/rspa.1932.0019
  * Ding & Pierrehumbert (2016) Ding, F., & Pierrehumbert, R. T. 2016, The Astrophysical Journal, 822, 24, doi: 10.3847/0004-637X/822/1/24
  * Edwards & Slingo (1996) Edwards, J. M., & Slingo, A. 1996, Quarterly Journal of the Royal Meteorological Society, 122, 689, doi: 10.1002/qj.49712253107
  * Fulton et al. (2017) Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, The Astronomical Journal, 154, 109, doi: 10.3847/1538-3881/aa80eb
  * Goldblatt et al. (2013) Goldblatt, C., Robinson, T. D., Zahnle, K. J., & Crisp, D. 2013, Nature Geoscience, 28, doi: 10.1038/NGEO1892
  * Guillot (1995) Guillot, T. 1995, Science, 269, 1697, doi: 10.1126/science.7569896
  * Gupta & Schlichting (2019) Gupta, A., & Schlichting, H. E. 2019, Monthly Notices of the Royal Astronomical Society, 487, 24, doi: 10.1093/mnras/stz1230
  * Heng et al. (2014) Heng, K., Mendonça, J. M., & Lee, J. M. 2014, Astrophysical Journal, Supplement Series, 215, doi: 10.1088/0067-0049/215/1/4
  * Husser et al. (2013) Husser, T. O., Wende-Von Berg, S., Dreizler, S., et al. 2013, Astronomy and Astrophysics, 553, 6, doi: 10.1051/0004-6361/201219058
  * Ingersoll (1969) Ingersoll, A. 1969, Journal of the Atmospheric Sciences, 26, 1191, doi: 10.1175/1520-0469(1969)026<1191:trgaho>2.0.co;2
  * Ingersoll & Kanamori (1995) Ingersoll, A. P., & Kanamori, H. 1995, Nature, 374, 706, doi: 10.1038/374706a0
  * Innes & Pierrehumbert (2022) Innes, H., & Pierrehumbert, R. T. 2022, The Astrophysical Journal, 927, 38, doi: 10.3847/1538-4357/ac4887
  * Irwin et al. (2022) Irwin, P. G. J., Teanby, N. A., Fletcher, L. N., et al. 2022, Journal of Geophysical Research: Planets, 127, e2022JE007189, doi: 10.1029/2022JE007189
  * Karman et al. (2019) Karman, T., Gordon, I. E., van der Avoird, A., et al. 2019, Icarus, 328, 160, doi: 10.1016/j.icarus.2019.02.034
  * Kasting et al. (1993) Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icarus, 101, 108, doi: 10.1006/icar.1993.1010
  * Koll (2022) Koll, D. D. B. 2022, The Astrophysical Journal, 924, 134, doi: 10.3847/1538-4357/ac3b48
  * Koll & Cronin (2019) Koll, D. D. B., & Cronin, T. W. 2019, The Astrophysical Journal, 881, 120, doi: 10.3847/1538-4357/ab30c4
  * Komabayasi (1967) Komabayasi, M. 1967, Journal of the Meteorological Society of Japan. Ser. II, 45, 137, doi: 10.2151/jmsj1965.45.1_137
  * Kopparapu (2013) Kopparapu, R. K. 2013, Astrophysical Journal Letters, 767, 8, doi: 10.1088/2041-8205/767/1/L8
  * Kopparapu et al. (2013) Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013, Astrophysical Journal, 765, 131, doi: 10.1088/0004-637X/765/2/131
  * Kuchner (2003) Kuchner, M. J. 2003, The Astrophysical Journal, 596, L105, doi: 10.1086/378397
  * Lacis & Oinas (1991) Lacis, A. A., & Oinas, V. 1991, Journal of Geophysical Research, 96, 9027, doi: 10.1029/90JD01945
  * Lean & DeLand (2012) Lean, J. L., & DeLand, M. T. 2012, Journal of Climate, 25, 2555, doi: 10.1175/JCLI-D-11-00571.1
  * Leconte et al. (2013) Leconte, J., Forget, F., Charnay, B., Wordsworth, R., & Pottier, A. 2013, Nature, 504, 268, doi: 10.1038/nature12827
  * Leconte et al. (2017) Leconte, J., Selsis, F., Hersant, F., & Guillot, T. 2017, Astronomy and Astrophysics, 598, 98, doi: 10.1051/0004-6361/201629140
  * Lefèvre et al. (2021) Lefèvre, M., Turbet, M., & Pierrehumbert, R. 2021, The Astrophysical Journal, 913, 101, doi: 10.3847/1538-4357/abf2c1
  * Li & Ingersoll (2015) Li, C., & Ingersoll, A. P. 2015, Nature Geoscience, 8, 398, doi: 10.1038/ngeo2405
  * Lopez et al. (2012) Lopez, E. D., Fortney, J. J., & Miller, N. 2012, Astrophysical Journal, 761, 59, doi: 10.1088/0004-637X/761/1/59
  * Léger et al. (2004) Léger, A., Selsis, F., Sotin, C., et al. 2004, Icarus, 169, 499, doi: 10.1016/j.icarus.2004.01.001
  * Madhusudhan et al. (2020) Madhusudhan, N., Nixon, M. C., Welbanks, L., Piette, A. A. A., & Booth, R. A. 2020, The Astrophysical Journal Letters, 891, L7, doi: 10.3847/2041-8213/ab7229
  * Madhusudhan et al. (2021) Madhusudhan, N., Piette, A. A. A., & Constantinou, S. 2021, The Astrophysical Journal, 918, 1, doi: 10.3847/1538-4357/abfd9c
  * Malik et al. (2017) Malik, M., Grosheintz, L., Mendonça, J. M., et al. 2017, The Astronomical Journal, 153, 56, doi: 10.3847/1538-3881/153/2/56
  * Marcy et al. (2014) Marcy, G. W., Weiss, L. M., Petigura, E. A., et al. 2014, Proceedings of the National Academy of Sciences of the United States of America, 111, 12655, doi: 10.1073/pnas.1304197111
  * Markham et al. (2022) Markham, S., Guillot, T., & Stevenson, D. 2022, Astronomy & Astrophysics, 665, A12, doi: 10.1051/0004-6361/202243359
  * Mehio et al. (2014) Mehio, N., Dai, S., & Jiang, D.-e. 2014, The Journal of Physical Chemistry A, 118, 1150, doi: 10.1021/jp412588f
  * Misener & Schlichting (2022) Misener, W., & Schlichting, H. E. 2022, Monthly Notices of the Royal Astronomical Society, 514, 6025, doi: 10.1093/mnras/stac1732
  * Mlawer et al. (2012) Mlawer, E. J., Payne, V. H., Moncet, J. L., et al. 2012, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 370, 2520, doi: 10.1098/rsta.2011.0295
  * Mousis et al. (2020) Mousis, O., Deleuil, M., Aguichine, A., et al. 2020, The Astrophysical Journal, 896, L22, doi: 10.3847/2041-8213/ab9530
  * Nakajima et al. (1992) Nakajima, S., Hayashi, Y.-Y., & Abe, Y. 1992, Journal of the Atmospheric Sciences, 49, 2256, doi: 10.1175/1520-0469(1992)049<2256:ASOTGE>2.0.CO;2
  * Nixon & Madhusudhan (2021) Nixon, M. C., & Madhusudhan, N. 2021, Monthly Notices of the Royal Astronomical Society, 505, 3414, doi: 10.1093/mnras/stab1500
  * Owen & Wu (2013) Owen, J. E., & Wu, Y. 2013, Astrophysical Journal, 775, 105, doi: 10.1088/0004-637X/775/2/105
  * Peck & Huang (1977) Peck, E. R., & Huang, S. 1977, Journal of the Optical Society of America, 67, 1550, doi: 10.1364/josa.67.001550
  * Petigura et al. (2013) Petigura, E. A., Howard, A. W., & Marcy, G. W. 2013, Proceedings of the National Academy of Sciences of the United States of America, 110, 19273, doi: 10.1073/pnas.1319909110
  * Piaulet et al. (2023) Piaulet, C., Benneke, B., Almenara, J. M., et al. 2023, Nature Astronomy, 7, 206, doi: 10.1038/s41550-022-01835-4
  * Pierrehumbert & Gaidos (2011) Pierrehumbert, R., & Gaidos, E. 2011, Astrophysical Journal Letters, 734, 13, doi: 10.1088/2041-8205/734/1/L13
  * Pierrehumbert (2010) Pierrehumbert, R. T. 2010, Principles of Planetary Climate (Cambridge: Cambridge University Press), doi: 10.1017/cbo9780511780783
  * Pierrehumbert (2023) —. 2023, The Astrophysical Journal, 944, 20, doi: 10.3847/1538-4357/acafdf
  * Pierrehumbert & Ding (2016) Pierrehumbert, R. T., & Ding, F. 2016, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 472, doi: 10.1098/rspa.2016.0107
  * Pierrehumbert & Swanson (1995) Pierrehumbert, R. T., & Swanson, K. L. 1995, Annu. Rev. Fluid Mech, 27, 419, doi: 10.1146/annurev.fluid.27.1.419
  * Piette & Madhusudhan (2020) Piette, A. A., & Madhusudhan, N. 2020, The Astrophysical Journal, 904, 154, doi: 10.3847/1538-4357/abbfb1
  * Rogers (2015) Rogers, L. A. 2015, Astrophysical Journal, 801, doi: 10.1088/0004-637X/801/1/41
  * Rogers et al. (2011) Rogers, L. A., Bodenheimer, P., Lissauer, J. J., & Seager, S. 2011, Astrophysical Journal, 738, 59, doi: 10.1088/0004-637X/738/1/59
  * Rogers & Seager (2010) Rogers, L. A., & Seager, S. 2010, Astrophysical Journal, 716, 1208, doi: 10.1088/0004-637X/716/2/1208
  * Schwarzschild (1906) Schwarzschild, K. 1906, Nachrichten von der Königlichen Gesellschaft der Wissenschaften zu Göttingen. Math.-phys. Klasse, 195, 41.  https://ui.adsabs.harvard.edu/abs/1906WisGo.195...41S
  * Simpson (1929) Simpson, G. 1929, Quarterly Journal of the Royal Meteorological Society, 55, 73, doi: 10.1002/qj.49705522908
  * Soubiran & Militzer (2015) Soubiran, F., & Militzer, B. 2015, The Astrophysical Journal, 806, 228, doi: 10.1088/0004-637X/806/2/228
  * Tan et al. (2021) Tan, X., Lefèvre, M., & Pierrehumbert, R. T. 2021, The Astrophysical Journal Letters, 923, L15, doi: 10.3847/2041-8213/ac3e69
  * Tsai et al. (2021) Tsai, S.-M., Innes, H., Lichtenberg, T., et al. 2021, The Astrophysical Journal Letters, 922, L27, doi: 10.3847/2041-8213/ac399a
  * Tsiaras et al. (2019) Tsiaras, A., Waldmann, I. P., Tinetti, G., Tennyson, J., & Yurchenko, S. N. 2019, Nature Astronomy, 3, 1086, doi: 10.1038/s41550-019-0878-9
  * Turbet et al. (2021) Turbet, M., Bolmont, E., Chaverot, G., et al. 2021, Nature, 598, doi: 10.1038/s41586-021-03873-w
  * Wagner & Pruß (2002) Wagner, W., & Pruß, A. 2002, Journal of Physical and Chemical Reference Data, 31, 387, doi: 10.1063/1.1461829
  * Winn & Fabrycky (2015) Winn, J. N., & Fabrycky, D. C. 2015, Annu. Rev. Astron. Astrophys, 53, 409, doi: 10.1146/annurev-astro-082214-122246
  * Yang et al. (2013) Yang, J., Cowan, N. B., & Abbot, D. S. 2013, Astrophysical Journal Letters, 771, doi: 10.1088/2041-8205/771/2/L45